Mathematical model of groundwater head variation processes in heterogeneous porous media

This study develops a mathematical formulation and numerical solution for determining the spatial–temporal distribution of the groundwater head in heterogeneous porous media, and illustrates them with a transient multi-well simulation. The relevance of the study stems from the need for efficient management of groundwater resources, for ecological sustainability, and for forecasting hydrological processes under a changing climate. Starting from Darcy’s law and the law of mass conservation, a second-order parabolic partial differential equation with variable coefficients is derived that describes filtration in anisotropic heterogeneous media; anisotropy is represented through a diagonal hydraulic-conductivity tensor and heterogeneity through spatially variable parameters. Boundary conditions are formulated for all six faces of the domain, including river–aquifer interaction, infiltration, evapotranspiration, and irrigation return flow. The differential problem is discretized by an explicit finite-difference scheme in which the interface conductivities in the heterogeneous medium are evaluated by the harmonic mean, which correctly reproduces series-connected hydraulic resistances. The stability of the scheme is analysed by the von Neumann method and a Courant–Friedrichs–Lewy stability condition is obtained; convergence follows from the Lax equivalence theorem. A transient numerical experiment is carried out for an aquifer containing two extraction wells and one injection well: the computations reproduce the progressive formation and interaction of the depression cones and the recharge mound, the transient head response at the wells, and the three-dimensional structure of the head surface. The computational algorithm parallelizes naturally and can be applied to groundwater-resource management, to forecasting the hydrodynamic state of aquifers, and to the design of wells and drainage systems.