Fluid flow in the Earth's crust can be accurately described by a diffusion equation, with highly anisotropic, heterogeneous and time dependent diffusivities. The so‐called Bhatnagar‐Gross‐Krook (BGK) models provide a simple and efficient way to solve the scalar diffusion equation with complex boundary conditions. We show that it can be adapted to include an anisotropic, space and time dependent diffusivity, and a gravity body force, hence providing a suitable numerical model for realistic simulations of incompressible pore fluid diffusion in heterogeneous, anisotropic rocks in the presence of a gravity field. A stability condition limits in some cases the magnitude of anisotropy that can be simulated. Numerical results are in very good agreement with a known analytical sol