Groundwater flow equation

Last updated

Used in hydrogeology, the groundwater flow equation is the mathematical relationship which is used to describe the flow of groundwater through an aquifer. The transient flow of groundwater is described by a form of the diffusion equation, similar to that used in heat transfer to describe the flow of heat in a solid (heat conduction). The steady-state flow of groundwater is described by a form of the Laplace equation, which is a form of potential flow and has analogs in numerous fields.


The groundwater flow equation is often derived for a small representative elemental volume (REV), where the properties of the medium are assumed to be effectively constant. A mass balance is done on the water flowing in and out of this small volume, the flux terms in the relationship being expressed in terms of head by using the constituitive equation called Darcy's law, which requires that the flow is laminar. Other approaches are based on Agent Based Models to incorporate the effect of complex aquifers such as karstic or fractured rocks (i.e. volcanic) [1]

Mass balance

A mass balance must be performed, and used along with Darcy's law, to arrive at the transient groundwater flow equation. This balance is analogous to the energy balance used in heat transfer to arrive at the heat equation. It is simply a statement of accounting, that for a given control volume, aside from sources or sinks, mass cannot be created or destroyed. The conservation of mass states that, for a given increment of time (Δt), the difference between the mass flowing in across the boundaries, the mass flowing out across the boundaries, and the sources within the volume, is the change in storage.

Diffusion equation (transient flow)

Mass can be represented as density times volume, and under most conditions, water can be considered incompressible (density does not depend on pressure). The mass fluxes across the boundaries then become volume fluxes (as are found in Darcy's law). Using Taylor series to represent the in and out flux terms across the boundaries of the control volume, and using the divergence theorem to turn the flux across the boundary into a flux over the entire volume, the final form of the groundwater flow equation (in differential form) is:

This is known in other fields as the diffusion equation or heat equation, it is a parabolic partial differential equation (PDE). This mathematical statement indicates that the change in hydraulic head with time (left hand side) equals the negative divergence of the flux (q) and the source terms (G). This equation has both head and flux as unknowns, but Darcy's law relates flux to hydraulic heads, so substituting it in for the flux (q) leads to

Now if hydraulic conductivity (K) is spatially uniform and isotropic (rather than a tensor), it can be taken out of the spatial derivative, simplifying them to the Laplacian, this makes the equation

Dividing through by the specific storage (Ss), puts hydraulic diffusivity (α = K/Ss or equivalently, α = T/S) on the right hand side. The hydraulic diffusivity is proportional to the speed at which a finite pressure pulse will propagate through the system (large values of α lead to fast propagation of signals). The groundwater flow equation then becomes

Where the sink/source term, G, now has the same units but is divided by the appropriate storage term (as defined by the hydraulic diffusivity substitution).

Rectangular cartesian coordinates

Three-dimensional finite difference grid used in MODFLOW MODFLOW 3D grid.png
Three-dimensional finite difference grid used in MODFLOW

Especially when using rectangular grid finite-difference models (e.g. MODFLOW, made by the USGS), we deal with Cartesian coordinates. In these coordinates the general Laplacian operator becomes (for three-dimensional flow) specifically

MODFLOW code discretizes and simulates an orthogonal 3-D form of the governing groundwater flow equation. However, it has an option to run in a "quasi-3D" mode if the user wishes to do so; in this case the model deals with the vertically averaged T and S, rather than k and Ss. In the quasi-3D mode, flow is calculated between 2D horizontal layers using the concept of leakage.

Circular cylindrical coordinates

Another useful coordinate system is 3D cylindrical coordinates (typically where a pumping well is a line source located at the origin parallel to the z axis causing converging radial flow). Under these conditions the above equation becomes (r being radial distance and θ being angle),


This equation represents flow to a pumping well (a sink of strength G), located at the origin. Both this equation and the Cartesian version above are the fundamental equation in groundwater flow, but to arrive at this point requires considerable simplification. Some of the main assumptions which went into both these equations are:

Despite these large assumptions, the groundwater flow equation does a good job of representing the distribution of heads in aquifers due to a transient distribution of sources and sinks.

Laplace equation (steady-state flow)

If the aquifer has recharging boundary conditions a steady-state may be reached (or it may be used as an approximation in many cases), and the diffusion equation (above) simplifies to the Laplace equation.

This equation states that hydraulic head is a harmonic function, and has many analogs in other fields. The Laplace equation can be solved using techniques, using similar assumptions stated above, but with the additional requirements of a steady-state flow field.

A common method for solution of this equations in civil engineering and soil mechanics is to use the graphical technique of drawing flownets; where contour lines of hydraulic head and the stream function make a curvilinear grid, allowing complex geometries to be solved approximately.

Steady-state flow to a pumping well (which never truly occurs, but is sometimes a useful approximation) is commonly called the Thiem solution.

Two-dimensional groundwater flow

The above groundwater flow equations are valid for three dimensional flow. In unconfined aquifers, the solution to the 3D form of the equation is complicated by the presence of a free surface water table boundary condition: in addition to solving for the spatial distribution of heads, the location of this surface is also an unknown. This is a non-linear problem, even though the governing equation is linear.

An alternative formulation of the groundwater flow equation may be obtained by invoking the Dupuit–Forchheimer assumption, where it is assumed that heads do not vary in the vertical direction (i.e., ). A horizontal water balance is applied to a long vertical column with area extending from the aquifer base to the unsaturated surface. This distance is referred to as the saturated thickness, b. In a confined aquifer, the saturated thickness is determined by the height of the aquifer, H, and the pressure head is non-zero everywhere. In an unconfined aquifer, the saturated thickness is defined as the vertical distance between the water table surface and the aquifer base. If , and the aquifer base is at the zero datum, then the unconfined saturated thickness is equal to the head, i.e., b=h.

Assuming both the hydraulic conductivity and the horizontal components of flow are uniform along the entire saturated thickness of the aquifer (i.e., and ), we can express Darcy's law in terms of integrated groundwater discharges, Qx and Qy:

Inserting these into our mass balance expression, we obtain the general 2D governing equation for incompressible saturated groundwater flow:

Where n is the aquifer porosity. The source term, N (length per time), represents the addition of water in the vertical direction (e.g., recharge). By incorporating the correct definitions for saturated thickness, specific storage, and specific yield, we can transform this into two unique governing equations for confined and unconfined conditions:

(confined), where S=Ssb is the aquifer storativity and

(unconfined), where Sy is the specific yield of the aquifer.

Note that the partial differential equation in the unconfined case is non-linear, whereas it is linear in the confined case. For unconfined steady-state flow, this non-linearity may be removed by expressing the PDE in terms of the head squared:

Or, for homogeneous aquifers,

This formulation allows us to apply standard methods for solving linear PDEs in the case of unconfined flow. For heterogeneous aquifers with no recharge, Potential flow methods may be applied for mixed confined/unconfined cases.

See also

Related Research Articles

Laplaces equation Second order partial differential equation

In mathematics and physics, Laplace's equation is a second-order partial differential equation named after Pierre-Simon Laplace who first studied its properties. This is often written as

Navier–Stokes equations Equations describing the motion of viscous fluid substances

In physics, the Navier–Stokes equations are certain partial differential equations which describe the motion of viscous fluid substances, named after French engineer and physicist Claude-Louis Navier and Anglo-Irish physicist and mathematician George Gabriel Stokes. They were developed over several decades of progressively building the theories, from 1822 (Navier) to 1842–1850 (Stokes).

Potential flow

In fluid dynamics, potential flow describes the velocity field as the gradient of a scalar function: the velocity potential. As a result, a potential flow is characterized by an irrotational velocity field, which is a valid approximation for several applications. The irrotationality of a potential flow is due to the curl of the gradient of a scalar always being equal to zero.

In mathematics, the Laplace operator or Laplacian is a differential operator given by the divergence of the gradient of a scalar function on Euclidean space. It is usually denoted by the symbols , , or . In a Cartesian coordinate system, the Laplacian is given by the sum of second partial derivatives of the function with respect to each independent variable. In other coordinate systems, such as cylindrical and spherical coordinates, the Laplacian also has a useful form. Informally, the Laplacian Δf(p) of a function f at a point p measures by how much the average value of f over small spheres or balls centered at p deviates from f(p).

Stream function

The stream function is defined for incompressible (divergence-free) flows in two dimensions – as well as in three dimensions with axisymmetry. The flow velocity components can be expressed as the derivatives of the scalar stream function. The stream function can be used to plot streamlines, which represent the trajectories of particles in a steady flow. The two-dimensional Lagrange stream function was introduced by Joseph Louis Lagrange in 1781. The Stokes stream function is for axisymmetrical three-dimensional flow, and is named after George Gabriel Stokes.

Euler equations (fluid dynamics)

In fluid dynamics, the Euler equations are a set of quasilinear hyperbolic equations governing adiabatic and inviscid flow. They are named after Leonhard Euler. The equations represent Cauchy equations of conservation of mass (continuity), and balance of momentum and energy, and can be seen as particular Navier–Stokes equations with zero viscosity and zero thermal conductivity. In fact, Euler equations can be obtained by linearization of some more precise continuity equations like Navier–Stokes equations in a local equilibrium state given by a Maxwellian. The Euler equations can be applied to incompressible and to compressible flow – assuming the flow velocity is a solenoidal field, or using another appropriate energy equation respectively. Historically, only the incompressible equations have been derived by Euler. However, fluid dynamics literature often refers to the full set – including the energy equation – of the more general compressible equations together as "the Euler equations".

Darcy's law is an equation that describes the flow of a fluid through a porous medium. The law was formulated by Henry Darcy based on results of experiments on the flow of water through beds of sand, forming the basis of hydrogeology, a branch of earth sciences.

In differential geometry, the four-gradient is the four-vector analogue of the gradient from vector calculus.

Hydraulic conductivity, symbolically represented as , is a property of vascular plants, soils and rocks, that describes the ease with which a fluid can move through pore spaces or fractures. It depends on the intrinsic permeability of the material, the degree of saturation, and on the density and viscosity of the fluid. Saturated hydraulic conductivity, Ksat, describes water movement through saturated media. By definition, hydraulic conductivity is the ratio of velocity to hydraulic gradient indicating permeability of porous media.

Hydraulic head Specific measurement of liquid pressure above a vertical datum

Hydraulic head or piezometric head is a specific measurement of liquid pressure above a vertical datum.

Covariant formulation of classical electromagnetism

The covariant formulation of classical electromagnetism refers to ways of writing the laws of classical electromagnetism in a form that is manifestly invariant under Lorentz transformations, in the formalism of special relativity using rectilinear inertial coordinate systems. These expressions both make it simple to prove that the laws of classical electromagnetism take the same form in any inertial coordinate system, and also provide a way to translate the fields and forces from one frame to another. However, this is not as general as Maxwell's equations in curved spacetime or non-rectilinear coordinate systems.

In fluid mechanics, potential vorticity (PV) is a quantity which is proportional to the dot product of vorticity and stratification. This quantity, following a parcel of air or water, can only be changed by diabatic or frictional processes. It is a useful concept for understanding the generation of vorticity in cyclogenesis, especially along the polar front, and in analyzing flow in the ocean.

The Richards equation represents the movement of water in unsaturated soils, and is attributed to Lorenzo A. Richards who published the equation in 1931. It is a nonlinear partial differential equation, which is often difficult to approximate since it does not have a closed-form analytical solution. The equation is based on Darcy's law for groundwater flow, which was developed for saturated rather than unsaturated flow in porous media. To do this, an additional continuity requirement is added. The transient-state form of the Richards equation in just the vertical direction is

Groundwater discharge is the volumetric flow rate of groundwater through an aquifer.

Mathematical descriptions of the electromagnetic field Formulations of electromagnetism

There are various mathematical descriptions of the electromagnetic field that are used in the study of electromagnetism, one of the four fundamental interactions of nature. In this article, several approaches are discussed, although the equations are in terms of electric and magnetic fields, potentials, and charges with currents, generally speaking.

The intent of this article is to highlight the important points of the derivation of the Navier–Stokes equations as well as its application and formulation for different families of fluids.


HydroGeoSphere (HGS) is a 3D control-volume finite element groundwater model, and is based on a rigorous conceptualization of the hydrologic system consisting of surface and subsurface flow regimes. The model is designed to take into account all key components of the hydrologic cycle. For each time step, the model solves surface and subsurface flow, solute and energy transport equations simultaneously, and provides a complete water and solute balance.

ZOOMQ3D is a numerical finite-difference model, which simulates groundwater flow in aquifers. The program is used by hydrogeologists to investigate groundwater resources and to make predictions about possible future changes in their quantity and quality. The code is written in C++, an object-oriented programming language and can compile and run on Windows and Unix operating systems.

Hydrological model

A hydrologic model is a simplification of a real-world system that aids in understanding, predicting, and managing water resources. Both the flow and quality of water are commonly studied using hydrologic models.

Relativistic heat conduction refers to the modelling of heat conduction in a way compatible with special relativity. This article discusses models using a wave equation with a dissipative term.


  1. Corona, Oliver López; Padilla, Pablo; Escolero, Oscar; González, Tomas; Morales-Casique, Eric; Osorio-Olvera, Luis (2014-10-16). "Complex groundwater flow systems as traveling agent models". PeerJ. 2: e557. doi: 10.7717/peerj.557 . ISSN   2167-8359.

Further reading