This paper presents a lattice Boltzmann method for the advection and dispersion of solute in three-dimensional variably saturated porous media. The proposed method is based on the BGK model and discretizes the particle velocity space with a cuboid lattice in which the particles move in 19 directions and 7 speeds. In the proposed model a directionally dependent relaxation time is introduced to produce the second-order dispersion tensor, and a modification of the equilibrium distribution functions is given to model the solute in porous media where the volumetric water content varies over space and time. The concentration is calculated from a weighted summation of the particle distribution functions to ensure that the BGK collision does not result in mass change. The accuracy and consistency of the proposed model are verified against benchmark problems and the finite difference solution of solute transport in an unsaturated heterogeneous soil. The results show that the proposed model conserves mass perfectly and gives efficient and accurate solutions for both advection-dominated and dispersion-dominated problems.