A Lattice Boltzmann model with two relaxation times for the 2D/3D advection and anisotropic diffusion equation (AADE) is introduced. The method is applied to Richards' equation for variably saturated flow in isotropic homogeneous media by extending retention curves into the saturated zone in a linear manner. The Darcy velocity is computed locally from the population solution. The method possesses intrinsic mass conservation, it is explicit and especially suitable for parallel computations. Designed for regular grids, the LB approach meets the boundary conditions accurately with an unified “multi-reflexion” technique, introduced to fit pressure head and/or specified flux conditions on static and seepage boundaries. The physical space can assume an uniform rectangular discretization grid which is transformed into the cubic computational grid after proper rescaling  of the AADE. The diffusion term is considered in two forms: the conventional one and the transformed one. The integral transformation may avoid problems encountered with the unbounded diffusion coefficients at the residual and saturated limits. An analytical expression for the transformed diffusion function is obtained for the original and modified VGM retention curves . Analytical instationary solutions for constant flux infiltration with non-linear models [2,15] are revised. An exact unstationary solution  valid for unsaturated, saturated or variably saturated flow is constructed using the BCM hydraulic conductivity function [3,11], moisture tension being fixed at the surface. Stationary infiltration profiles are generated for the BCM and the VGM conductivity functions. The LB method is validated against these and other reference solutions.