We present analytical and Lattice Boltzmann (LB) solutions for steady-state saturated flows in 2D and 3D anisotropic heterogeneous aquifers. The analytical solution is easy to use and extends the known ones for ground-water whirls to more general combinations of the anisotropic properties of two-layered systems. The Bakker and Hemker's multi-layered' semi-analytical solution and the LB results are compared to the analytical solution for a broad range of anisotropic heterogeneous diffusion tensors. The main components of the LB scheme, the eigenvalues of the linear collision operator and/or the equilibrium functions, become discontinuous when the anisotropy changes between the layers. It is shown that the evolution equation of the LB method needs to be modified at the interfaces in order to satisfy the continuity conditions for the diffusion function and/or its tangential derivatives. The existing LB schemes for anisotropic advectiondispersion equations are formulated in a more general framework in which the leading-order interface corrections are constructed and analyzed for linear and highly nonlinear exact solutions. We also present some stability aspects of these schemes, introduce specified normal gradient boundary conditions and discuss the computation of total and local fluxes. The interface analysis developed here applies to generic LB schemes with discontinuous collision operators.