7.5 Applying Governing Equations
The general mathematical governing equations that describe groundwater flow conditions in a groundwater system are the foundation of models developed to characterize groundwater flow and its transport of dissolved substances. The governing equations describe how the hydraulic heads, h, within a problem domain are distributed under a set of specified conditions. The steps involved in formulating a groundwater model include:
- defining the purpose of the model;
- building a conceptual model of the groundwater system including a water balance based on field and laboratory data, and hydrogeologic theory;
- translation of the conceptual model into a format compatible with the appropriate governing equation and solution technique;
- application of the technique to obtain a solution;
- calibration (adjusting hydraulic parameters and boundary conditions until the solution matches heads and fluxes measured at the field site);
- using the calibrated model to make predictions about how the groundwater system will respond under a new stress such as pumping or a drought; and,
- reporting model results and their uncertainty.
The Role of a Water Budget in Formulating Models
In a few places in this book, the concept of a water budget or water balance has been noted as a key component used to understand groundwater systems. When scientists develop a conceptual model of a groundwater system, they examine the characteristics of the physical geologic framework, assess hydrogeological parameters, define problem boundaries, and generate a field-based groundwater budget that accounts for the volumes of groundwater entering and leaving the problem domain. The groundwater budget is a simple accounting of the groundwater inflows, outflows and changes in groundwater storage over a specified time interval for a defined problem domain (volume). In its simplest form it is expressed as Equation 81.
Groundwater Inflow = Groundwater Outflow + Change in Groundwater Storage | (81) |
Given the addition of the change in groundwater storage on the right-hand side of Equation 81, a sign convention is required for the change in groundwater storage. As presented here, “Change in Groundwater Storage” is negative for water coming out of groundwater storage (i.e., falling water levels) and positive for water going into groundwater storage (i.e., rising water levels). The sign convention may vary in other presentations. The convention for any presentation can be deduced by looking at the form of the expression and reasoning out which sign is required for the system to balance.
Groundwater inflow can include flow of water into the domain at boundaries including leakage from a deeper aquifer system; recharge by natural precipitation; agricultural irrigation and artificial basin recharge; and seepage from surface water features such as streams, lakes and wetlands. Outflow can include flow of water out of the domain at boundaries including leakage to a deeper aquifer system; groundwater flowing to natural discharge locations such as streams, lakes, wetlands and springs; loss of groundwater to evapotranspiration, and extraction of groundwater by water supply wells and drains. Changes in storage occur in response to natural or induced changes in the potentiometric surface. A conceptual groundwater budget is illustrated with words in Figure 56 and as images in Figure 57.


The conceptual budget is defined for a problem domain bordered by a low permeability fault (no groundwater flow into or out of the model); a permeable boundary along the Green Mountains where groundwater inflow occurs; a river boundary (Bright River) along which groundwater is recharged by river losses in the upper reaches and groundwater discharges to the gaining river in the lower reaches; and a lake where groundwater discharges to the lake as shown in Figure 57. Within the problem domain areal recharge occurs, as well as infiltration from an irrigated field. Groundwater leaves the system at a municipal well and by groundwater evapotranspiration at a wetland. The water budget is calculated over a one-year period so the storage is the product of the specific yield and the change in water level over the year (on average for the entire domain or the sum of values if the specific yield varies over the domain).
The average rates of inflow and outflow for the year are computed and summed. The values are obtained in a variety of ways. Recharge is estimated as the spatially averaged precipitation minus evapotranspiration. Groundwater flow in from the Green Mountains and out to Fish Lake are estimated by Darcy’s Law, using field measurements of gradients, hydraulic conductivity and flow areas. River seepage is estimated by noting the amount of gain and loss of stream flow between river gages. Municipal well discharge is obtained from town records. Irrigation recharge is estimated by knowing the amount of water imported for irrigation and subtracting the amount used by the type of crop in the climate of the area. Evapotranspiration from the wetland groundwater is estimated based on the local climate and character of the wetland vegetation. Change in groundwater storage is estimating by determining the average change in water levels in wells over the domain and multiplying by the area of the domain and the specific yield of the subsurface materials.
Many of the budget components are not easily determined, so they are estimates thus, the budget will not balance. If the imbalance is large, it is likely that the estimate for some component of the budget is in error, or an important component was overlooked and is not included. If the imbalance is small it may be distributed over the components of input or output, or considered to be part of the most uncertain component. When such a distribution is made, the water budget will be perfectly balanced, indicating the errors were incorporated into the component estimates. If the most uncertain component is determined by solving the budget equation, then the measurement error for all of the components is incorporated into the computed value. It is recommended that each component be measured or estimated independently and errors for each component estimated. A quantitative water budget provides the foundation for assessing whether analytical equations and numerical groundwater models make sense. A good reference for developing water budgets is available from the United States Geological Survey (Healy et al., 2007).
Boundary Value Problems
To apply the governing equation, a boundary value problem is formulated. It requires selecting a governing equation that describes the appropriate hydrogeological conditions, assigning head or flow rates at boundary conditions and assigning hydrogeologic parameters (e.g., hydraulic conductivity, recharge rates). Also, when conditions are transient, initial heads and storativity values are assigned. The boundary value problem is a model of the hydrogeological system.
Methods for Solving Groundwater Problems
Methods used to obtain solutions to groundwater problems include the following five approaches:
1. Inspection of field data
Inspection of the data for spatial and temporal trends and correlating data with the site-specific geological framework and basic hydrogeological principles allow for an initial analysis. This often includes building a conceptual model of how the groundwater system works and developing a water budget to assess interrelationships. A water budget lists the inputs and outputs (including water going into or coming out of storage) and estimates the magnitude of each budget item. The water budget must balance because mass is conserved. Inspection is also used to determine whether field data collected in simple settings contain trends that make it possible to predict future conditions.
2. Graphical techniques
Graphical methods include flow-net construction, a method discussed in Section 8.
3. Analog models
Analog models include using simple laboratory-scale systems made of materials that emulate groundwater flow under prescribed boundaries. These could be sand tanks or electric circuit boards.
4. Analytical mathematical techniques
Analytical solutions are algebraic equations that are continuous solutions to conditions described by a governing equation and assigned boundaries. A continuous solution is one that provides a value of head at every location in the system domain. Analytical solutions can be solved by hand (e.g., pencil-and-paper or hand calculator) or a spreadsheet. Analytical solutions usually require substantial simplifying assumptions and are typically only applicable to field systems with simple boundary conditions and flow patterns (for example, isotropic and homogeneous systems that can be represented by one- or two-dimensional flow).
5. Numerical mathematical techniques
Most field-scale groundwater modeling efforts require numerical methods to approximate solutions to a particular governing equation and generate values of head at discrete points in the groundwater system. The benefit of using such an approximation is that it can address complex boundary conditions and parameter distributions. Numerical modeling includes methods like finite difference and finite element modeling, which are used in many engineering fields. The United States Geological Survey has developed public domain numerical modeling software (e.g., https://water.usgs.gov/software/lists/groundwater) for simulating various aspects of groundwater systems.
Regardless of the solution technique, the governing equations are used to define the distribution of hydraulic head and fluxes within a system where the groundwater professional defines the spatial distribution of aquifer properties (e.g., K, and S) as well as boundary conditions for the volume of earth material. That volume is called the problem domain or model domain.
Boundary Conditions
Boundary conditions are generalized into three types that are defined along all boundaries of the model domain and can change with time. The three general types of boundary conditions are:
Type 1: Specified head boundary (also known as Dirichlet conditions). Steady or time-varying head values are specified along the boundary (for steady-state models the values are constant). If an unchanging value of head is defined on a boundary then the Type 1 boundary is referred to as a constant head boundary. Such a boundary could represent a large lakeshore where the lake stage elevation is assigned as a groundwater discharge boundary. When the model is solved, the model calculates flow rates at the location where the head is specified.
Type 2: Specified flux boundary (also known as Neumann conditions). Steady or time-varying flux is defined along the boundary (for steady-state models the values are constant). A special case is a no-flow or zero flux boundary. The gradient across a no-flow boundary is 0, so no flow crosses the boundary, Q = 0. A low permeability fault zone bordering a model domain could be assigned as a no-flow boundary. When the model is solved, the model calculates head at the location where the flux is specified.
Type 3: Head-dependent flux boundary (also known as Cauchy conditions; also referred to as a mixed boundary). Steady or time-varying head values external to the model domain are specified along the boundary and a resistance layer is defined between the model domain and the external head (for steady-state models the head values are constant). Flow across the boundary is controlled by the difference between the specified head value located outside of the domain and the head calculated at the boundary of the model domain. The calculated head depends on conditions in the rest of the model domain. As the calculated head at the boundary varies, so does the flux in or out of the boundary. A Type 3 boundary could be used to represent flow into a groundwater system from alluvial sediments underlying a large river. When the model is solved, the model calculates the time-varying flow rate at each location where the head dependent flux condition is specified.
Franke et al. (1987) provide a well-written discussion of how groundwater boundaries are formulated. A few conceptual examples of the application of models to groundwater problems are provided here.
Application of Flow Equations (Unconfined Aquifer Flow Between Water Bodies)
This subsection shows how an algebraic analytical solution is developed. Groundwater professionals usually do not need to develop their own analytical solutions because many solutions have been created and published. To use a solution the professional need only understand the hydrogeological conditions being represented and how mathematical simplifications effect its application to a particular site being investigated.
Consider the steady-state, isotropic and homogeneous, one-dimensional, section of an unconfined aquifer shown in Figure 58. The constant head reservoirs on the left and right can represent any water body, for example, a reservoir, lake or river and will be represented as Type 1 specified head boundaries. To be completely accurate, the water bodies need to fully penetrate the aquifer, however if the partial penetration provides good connection, this requirement is commonly overlooked and reasonable results are obtained.

Equation 80 is applicable to this situation. When unconfined flow is formulated in one dimension, the y term is omitted, the head is measured from the aquifer base (because it also represents aquifer thickness), and the partial differential is changed to an ordinary differential becoming Equation 82. This equation represents groundwater conditions for the problem domain, the area bounded by the two water bodies, the water table and aquifer base.
![]() |
(82) |
Those familiar with differential equations will see that the general solution to Equation 82 by integration is shown as Equation 83. Those who are unfamiliar with differential equations can, for now, accept that this is the case because this step is not vital to understanding and use of the final solution.
hx2 = C1x + C2 | (83) |
where:
hx | = | head at location x (L) |
C1, C2 | = | constants of integration (L, L2, respectively) |
x | = | position on the x axis (L) |
By applying the boundary condition h = h1 at x = 0 on the left face (a Type 1 boundary) and the condition h = h2 at x = L on the right face (also a Type 1 boundary), the values of the constants of integration can be determined. By substituting at x = 0, the first term, C1, of Equation 83 is zero and so, C2 = h12. As shown in Equation 84.
h12 = C2 | (84) |
Knowing the value of C2 and substituting h = h2 at x = L into Equation 83 produces Equation 85.
h22= C1L + h12 | (85) |
![]() |
Returning to Equation 83 and substituting Equation 84 and Equation 85 for C1 and C2, respectively, results in Equation 86, which is then rearranged with respect to the order of h1 and h2, and solved for hx (the head at any distance x from x = 0) in any of the three forms shown in Equation 86.
![]() |
(86) |
![]() |
![]() |
The heads at any x location in the system can be derived using Equation 86. For example, referring to Figure 58, if h1 is 15 m, h2 is 10 m, and L is 1000 m, what would be the value of head (hx) at x = 343 m? Using the third representation of Equation 86, the head at x = 343 m is calculated as follows.
![]() |
For a similar situation, the groundwater professional need not develop an analytical solution from scratch, because an existing solution is available. Once an appropriate existing analytical solution is identified, problem-specific parameters are inserted, and a head is calculated.
The next task is evaluating the flow rate per unit width, q’. Expressing the flow area as the product of the height of the water table and the unit width of the system in the y direction, the flow is shown in Equation 87.
![]() |
(87) |
With Δy = 1 it can be removed from the expression, so that q’ is the flow per unit width (Δy = 1) with units of L2/T (Equation 88). The symbol q’ is used here to differentiate it from specific discharge, q (L/T).
![]() |
(88) |
where:
q’ | = | flow per unit width (Δy = 1) (L2/T) |
Rearranging and integrating from x = 0 to x = L and from h at 0 = h1 to h at L = h2, as shown in Equation 89, results in Equation 90.
![]() ![]() |
(89) |
![]() ![]() |
(90) |
Again, using Figure 58, if h1 is 15 m, h2 is 10 m, L is 1000 m, and K = 25 m/d then the flow per unit width, q’, can be computed using Equation 90 as follows.
![]() |
If the cross-sectional area perpendicular to the boundaries is 100 m wide, the discharge of groundwater though the 100 m wide aquifer would be q’ (100 m) = 156 m3/d.
Box 6 extends the boundary value problem by adding constant recharge to the one-dimensional model. Box 6 also includes a simple spreadsheet that allows the reader to input values for this problem, including recharge rates and quickly see the head distribution and flow rates for the groundwater system. Click here to read Box 6 and experiment with the interactive spreadsheet.
It is important to remember that in this unconfined system the heads not only determine the gradient, but also the saturated thickness so the aquifer base must be used as the datum if a reasonable flow rate is to be obtained. At times, a hydrologist forgets this issue and calculates outrageous values for flow rate because sea level is used as a datum for a system that is at a high elevation with an aquifer base only 10 meters below the surface (e.g., perhaps at 2000 m on a continental divide instead of 10 m). Also, because the head, as measured from the aquifer base, determines the saturated thickness and the gradient, error will be introduced if the Dupuit assumptions are applied to an aquifer with a sloping base.
It should be noted that when the water table slope is small approximations of q and v are often derived by directly using Darcy’s law and assuming confined conditions as a reasonable approximation. In this case the confined conditions state the change in the water table is linear and it is assumed the thickness does not change significantly (T is constant).
Example Numerical Application of Flow Equations to a Dewatering Problem
This subsection shows how a numerical groundwater model, that incorporates a governing equation to describe groundwater conditions, can be applied to a field problem. Groundwater professionals often use software developed by mathematicians and computer programmers to solve problems. There are a number of standard, open-source, groundwater modeling software packages available to groundwater professionals (e.g., https://water.usgs.gov/software/lists/groundwater/). The professional does not need have all the skills needed to develop new software, but does need to understand how to apply models to solve problems.
The model presented here is intended to provide the reader with a general idea of a numerical model application in order to illustrate how groundwater systems are simplified for modeling and appreciate what they can achieve after mastering knowledge of groundwater science.
Often, large construction projects or mining endeavors need to dewater a groundwater system to accomplish their work. This example illustrates a situation in which a mining company wants to construct an exploratory mine shaft and needs to keep the shaft dewatered by maintaining a water level in the shaft at an elevation of 900 m. The shaft will be sealed through overlying formations including a thick shale layer overlying the confined aquifer in which the shaft terminates.
The model formulation is illustrated in Figure 59. To allow flexibility in assigning parameter values to the problem domain, the governing equation for confined three-dimensional, anisotropic, heterogeneous hydraulic conductivity distributions (Kx, Ky, Kz) is selected (Equation 67). A long-term pumping rate is sought, so a steady-state formulation of the problem is selected (left-hand side is set to 0). In addition to the flow equation, assignment of internal hydrogeologic parameters at multiple locations (e.g., Kx, Ky, Kz), and the domain boundaries are specified based on information from field-based hydrogeologic investigations. In this hypothetical example, water levels from boreholes in the field are used to estimate values of head that are specified on the left and right boundaries (Type 1). Recharge is estimated from local precipitation and soil conditions and is applied as a constant flux at the top of the model. The unconfined aquifer is also represented using Equation 67, however the computer code continually changes the unconfined hydraulic properties (T = Kb) as the saturated thickness of the unconfined layer changes. The boundaries on the front and back faces and the base are assigned as Type 2, zero flux, no-flow, boundaries, where water cannot pass into or out of the model domain. By defining the head on the bottom of the shaft as a constant 900 m, the model will calculate the outflow required to maintain that water level (large red arrow in Figure 59).

Finally, to solve the governing equation, a public-domain computer code that uses a numerical mathematical technique to generate heads and a water budget is applied. Figure 60 schematically illustrates the model results including: head distribution within each of the three geologic layers; flow rates into the left and right boundaries; leakage rates between model layers; and the answer to the question that prompted the model formulation (i.e., the steady-state pumping rate required to maintain the water level in the shaft at 900 m). The software also provides a water budget that accounts for the volume of water flowing in and out of every model boundary and from applied sources (recharge) and sinks (water extracted at the shaft). The water budget must balance (for steady-state conditions, inflow must equal outflow). A water budget that does not balance indicates a problem with the model formulation and the modeler must not use the results, but rather determine the causes and modify conditions to obtain a reasonable water balance.

If the mining company wanted to know how long they would need to pump in order to lower the water level to the bottom of the shaft, they could construct a model using the transient formulation as presented in Equation 67 and begin pumping the shaft at a fixed rate to explore how pumping at various rates would impact water levels in the vicinity of the shaft.