# 4. Model Equations¶

In this chapter, we discuss the model equations used by ParFlow for its
fully and variably saturated flow, overland flow, and multiphase flow
and transport models. First, section Steady-State, Saturated Groundwater Flow describes
steady-state, groundwater flow (specified by solver **IMPES**). Next,
section Richards’ Equation describes the Richards’ equation
model (specified by solver **RICHARDS**) for variably saturated flow as
implemented in ParFlow. Section Terrain Following Grid describes the terrain
following grid formulation. Next, the overland flow equations are
presented in section Overland Flow. In section
Multi-Phase Flow Equations we describe the multi-phase flow
equations (specified by solver **IMPES**), and in section
Transport Equations we describe the transport equations.
Finally, section Notation and Units presents some notation
and units and section Water Balance presents some basic water
balance equations.

## 4.1. Steady-State, Saturated Groundwater Flow¶

Many groundwater problems are solved assuming steady-state, fully-saturated groundwater flow. This follows the form often written as:

where \(Q\) is the spatially-variable source-sink term (to represent wells, etc) and \(\textbf{q}\) is the Darcy flux \([L^{2}T^{-1}]\) which is commonly written as:

where \(\textbf{K}\) is the saturated, hydraulic conductivity tensor \([LT^{-1}]\) and \(H\) \([L]\) is the head-potential. Inspection of (4.18) and (4.19) show that these equations agree with the above formulation for a single-phase (\(i=1\)), fully-saturated (\(S_i=S=1\)), problem where the mobility, \({\lambda}_i\), is set to the saturated hydraulic conductivity, \(\textbf{K}\), below. This is accomplished by setting the relative permeability and viscosity terms to unity in (4.20) as well as the gravity and density terms in (4.19). This is shown in the example in Annotated Input Scripts, but please note that the resulting solution is in pressure-head, \(h\), not head potential, \(H\), and will still contain a hydrostatic pressure gradient in the \(z\) direction.

## 4.2. Richards’ Equation¶

The form of Richards’ equation implemented in ParFlow is given as,

where \(\Omega\) is the flow domain, \(p\) is the pressure-head of water \([L]\), \(S\) is the water saturation, \(S_s\) is the specific storage coefficient \([L^{-1}]\), \(\phi\) is the porosity of the medium, \(\textbf{K}(p)\) is the hydraulic conductivity tensor \([LT^{-1}]\), and \(Q\) is the water source/sink term \([L^{3}T^{-1}]\) (includes wells and surface fluxes). The hydraulic conductivity can be written as,

Boundary conditions can be stated as,

where \(\Gamma^D \cup \Gamma^N = \partial \Omega\), \(\Gamma^D \neq \emptyset\), and \({\bf n}\) is an outward pointing, unit, normal vector to \(\Omega\). This is the mixed form of Richards’ equation. Note here that due to the constant (or passive) air phase pressure assumption, Richards’ equation ignores the air phase except through its effects on the hydraulic conductivity, \(K\). An initial condition,

completes the specification of the problem.

## 4.3. Terrain Following Grid¶

The terrain following grid formulation transforms the ParFlow grid to conform to topography [Max13]. This alters the form of Darcy’s law to include a topographic slope component:

where \(\theta_x = \arctan(S_0,x)\) and
\(\theta_y = \arctan(S_0,y)\) which are assumed to be the same as
the **TopoSlope** keys assigned for overland flow, described below. The
terrain following grid formulation can be very useful for coupled
surface-subsurface flow problems where groundwater flow follows the
topography. As cells are distributed near the ground surface and can be
combined with the variable \(\delta Z\) capability, the number of
cells in the problem can be reduced dramatically over the orthogonal
formulation. For complete details on this formulation, the stencil used
and the function evaluation developed, please see Maxwell [Max13]. NOTE: in the original formulation,
\(\theta_x\) and \(\theta_y\) for a cell face is calculated as
the average of the two adjacent cell slopes (i.e. assuming a cell
centered slope calculation). The
**TerrainFollowingGrid.SlopeUpwindFormulation** key provide options to
use the slope of a grid cell directly (i.e. assuming face centered slope
calculations) and removing the sine term from
(4.7). The **Upwind** and **UpwindSine**
options for this key will provide consistent results with
**OverlandKinematic** and **OverlandDiffusive** boundary conditions
while the **Original** option is consistent with the standard
**OverlandFlow** boundary condition.

## 4.4. Flow Barriers¶

The the flow barrier multipliers allow for the reduction in flow across a cell face. This slightly alters Darcy’s law to include a flow reduction in each direction, show here in x:

where \(FB_x\), \(FB_y\) and \(FB_z\) are a dimensionless
multipliers specified by the **FBx**, **FBy** and **FBz** keys. This
creates behavior equivalent to the Hydraulic Flow Barrier (HFB) or
*ITFC* (flow and transport parameters at interfaces) conditions in other
models.

## 4.5. Overland Flow¶

As detailed in Kollet and Maxwell [KM06], ParFlow may simulate fully-coupled surface and subsurface flow via an overland flow boundary condition. While complete details of this approach are given in that paper, a brief summary of the equations solved are presented here. Shallow overland flow is now represented in ParFlow by the kinematic wave equation. In two spatial dimensions, the continuity equation can be written as:

where \({\vec v}\) is the depth averaged velocity vector \([LT^{-1}]\); \(\psi_s\) is the surface ponding depth \([L]\) and \(q_r(x)\) is the a general source/sink (e.g. rainfall) rate \([LT^{-1}]\). If diffusion terms are neglected the momentum equation can be written as:

which is commonly referred to as the kinematic wave approximation. In Equation (4.10) \(S_{o,i}\) is the bed slope (gravity forcing term) \([-]\), which is equal to the friction slope \(S_{f,i}\) \([L]\); \(i\) stands for the \(x\)- and \(y\)-direction. Manning’s equation is used to establish a flow depth-discharge relationship:

and

where \(n\) \([TL^{-1/3}]\) is the Manning’s coefficient. Though complete details of the coupled approach are given in Kollet and Maxwell [KM06], brief details of the approach are presented here. The coupled approach takes Equation eq:kinematic and adds a flux for subsurface exchanges, \(q_e(x)\).

We then assign a continuity of pressure at the top cell of the boundary between the surface and subsurface systems by setting pressure-head, \(p\) in (4.3) equal to the vertically-averaged surface pressure, \(\psi_s\) as follows:

If we substitute this relationship back into Equation (4.13) as follows:

Where the \(\parallel\psi,0\parallel\) operator chooses the greater of the two quantities, \(\psi\) and \(0\). We may now solve this term for the flux \(q_e(x)\) which we may set equal to flux boundary condition shown in Equation (4.5). This yields the following equation, which is referred to as the overland flow boundary condition [KM06]:

This results a version of the kinematic wave equation that is only active when the pressure at the top cell of the subsurface domain has a ponded depth and is thus greater than zero. This method solves both systems, where active in the domain, over common grids in a fully-integrated, fully-mass conservative manner.

The depth-discharge relationship can also be written as

where \(\overline{S_{f}}\) is the magnitude of the friction slope.
This formulation for overland flow is used in the **OverlandKinematic**
and **OverlandDiffusive** boundary conditions. In **OverlandKinematic**
case the friction slope equals the bed slope following Equation
(4.10). For the **OverlandDiffusive** case the
friction slope also includes the pressure gradient. The solution for
both of these options is formulated to do the upwinding internally and
assumes that the user provides face centered bedslopes
(\(S_{o,i}\)). This is different from the original formulation which
assumes the user provides grid cenered bedslopes.

## 4.6. Multi-Phase Flow Equations¶

The flow equations are a set of *mass balance* and *momentum balance*
(Darcy’s Law) equations, given respectively by,

for \(i = 0, \ldots , \nu- 1\) \((\nu\in \{1,2,3\})\), where

Table 5.1 defines the symbols in the above equations, and outlines the symbol dependencies and units.

symbol |
quantity |
units |
---|---|---|

\(\phi({\vec x},t)\) |
porosity |
[] |

\(S_i({\vec x},t)\) |
saturation |
[] |

\({ \vec V}_i({\vec x},t)\) |
Darcy velocity vector |
[\(L T^{-1}\)] |

\(Q_i({\vec x},t)\) |
source/sink |
[\(T^{-1}\)] |

\({\lambda}_i\) |
mobility |
[\(L^{3} T M^{-1}\)] |

\(p_i({\vec x},t)\) |
pressure |
[\(M L^{-1} T^{-2}\)] |

\(\rho_i\) |
mass density |
[\(M L^{-3}\)] |

\({\vec g}\) |
gravity vector |
[\(L T^{-2}\)] |

\({ \bar k}({\vec x},t)\) |
intrinsic permeability tensor |
[\(L^{2}\)] |

\(k_{ri}({\vec x},t)\) |
relative permeability |
[] |

\(\mu_i\) |
viscosity |
[\(M L^{-1} T^{-1}\)] |

\(g\) |
gravitational acceleration |
[\(L T^{-2}\)] |

Here, \(\phi\) describes the fluid capacity of the porous medium,
and \(S_i\) describes the content of phase \(i\) in the porous
medium, where we have that \(0 \le \phi\le 1\) and
\(0 \le S_i\le 1\). The coefficient \({\bar k}\) is considered a
scalar here. We also assume that \(\rho_i\) and \(\mu_i\) are
constant. Also note that in ParFlow, we assume that the relative
permeability is given as \(k_{ri}(S_i)\). The Darcy velocity vector
is related to the *velocity vector*, \({\vec v}_i\), by the
following:

To complete the formulation, we have the following \(\nu\)
*consititutive relations*

where, \(p_{ij} = p_i - p_j\) is the *capillary pressure* between
phase \(i\) and phase \(j\). We now have the \(3 \nu\)
equations, (4.18), (4.19), (4.22), and
(4.23), in the
\(3 \nu\) unknowns, \(S_i, {\vec V}_i\), and \(p_i\).

For technical reasons, we want to rewrite the above equations. First, we
define the *total mobility*, \({\lambda}_T\), and the *total
velocity*, \({\vec V}_T\), by the relations

After doing a bunch of algebra, we get the following equation for \(p_0\):

After doing some more algebra, we get the following \(\nu- 1\) equations for \(S_i\):

The capillary pressures \(p_{ji}\) in (4.27) are rewritten in terms of the constitutive relations in (4.23) so that we have

where by definition, \(p_{ii} = 0\). Note that equations (4.27) are analytically the same equations as in (4.18). The reason we rewrite them in this latter form is because of the numerical scheme we are using. We now have the \(3 \nu\) equations, (4.26), (4.27), (4.25), (4.19), and (4.23), in the \(3 \nu\) unknowns, \(S_i, {\vec V}_i\), and \(p_i\).

## 4.7. Transport Equations¶

The transport equations in ParFlow are currently defined as follows:

where \(i = 0, \ldots , \nu- 1\) \((\nu\in \{1,2,3\})\) is the number of phases, \(j = 0, \ldots , n_c- 1\) is the number of contaminants, and where \(c_{i,j}\) is the concentration of contaminant \(j\) in phase \(i\). Recall also, that \(\chi_A\) is the characteristic function of set \(A\), i.e. \(\chi_A(x) = 1\) if \(x \in A\) and \(\chi_A(x) = 0\) if \(x \not\in A\). Table 5.2 defines the symbols in the above equation, and outlines the symbol dependencies and units. The equation is basically a statement of mass conservation in a convective flow (no diffusion) with adsorption and degradation effects incorporated along with the addition of injection and extraction wells.

symbol |
quantity |
units |
---|---|---|

\(\phi({\vec x})\) |
porosity |
[] |

\(c_{i,j}({\vec x},t)\) |
concentration fraction |
[] |

\({\vec V}_i({\vec x},t)\) |
Darcy velocity vector |
[\(L T^{-1}\)] |

\(\lambda_j\) |
degradation rate |
[\(T^{-1}\)] |

\(\rho_{s}({\vec x})\) |
density of the solid mass |
[\(M L^{-3}\)]] |

\(F_{i,j}({\vec x}, t)\) |
mass concentration |
[\(L^{3} M^{-1}\)] |

\(n_{I}\) |
number of injection wells |
[] |

\(\gamma^{I;i}_{k}(t)\) |
injection rate |
[\(T^{-1}\)] |

\(\Omega^{I}_{k}({\vec x})\) |
injection well region |
[] |

\({\bar c}^{k}_{ij}()\) |
injected concentration fraction |
[] |

\(n_{E}\) |
number of extraction wells |
[] |

\(\gamma^{E;i}_{k}(t)\) |
extraction rate |
[\(T^{-1}\)] |

\(\Omega^{E}_{k}({\vec x})\) |
extraction well region |
[] |

These equations will soon have to be generalized to include a diffusion term. At the present time, as an adsorption model, we take the mass concentration term (\(F_{i,j}\)) to be instantaneous in time and a linear function of contaminant concentration :

where \(K_{d;j}\) is the distribution coefficient of the component ([\(L^{3} M^{-1}\)]). If (4.30) is substituted into (4.29) the following equation results (which is the current model used in ParFlow) :

## 4.8. Notation and Units¶

In this section, we discuss other common formulations of the flow and transport equations, and how they relate to the equations solved by ParFlow.

We can rewrite equation (4.19) as

where

Table 5.3 defines the symbols and their units.

symbol |
quantity |
units |
---|---|---|

\({\vec V}_i\) |
Darcy velocity vector |
[\(L T^{-1}\)] |

\({\bar K}_i\) |
hydraulic conductivity tensor |
[\(L T^{-1}\)] |

\(h_i\) |
pressure head |
[\(L\)] |

\(\gamma\) |
constant scale factor |
[\(M L^{-2} T^{-2}\)] |

\({\vec g}\) |
gravity vector |
[\(L T^{-2}\)] |

We can then rewrite equations (4.26) and (4.27) as

Note that \({\bar K}_i\) is supposed to be a tensor, but we treat it as a scalar here. Also, note that by carefully defining the input to ParFlow, we can use the units of equations (4.34) and (4.35). To be more precise, let us denote ParFlow input symbols by appending the symbols in table 5.1 with \((I)\), and let \(\gamma= \rho_0 g\) (this is a typical definition). Then, we want:

By doing this, \({\bar k}(I)\) represents hydraulic conductivity of the base phase \({\bar K}_0\) (e.g. water) under saturated conditions (i.e. \(k_{r0} = 1\)).

## 4.9. Water Balance¶

ParFlow can calculate a water balance for the Richards’ equation,
overland flow and `clm`

capabilities. For a schematic of the water
balance in ParFlow please see Maxwell [Max10]. This water balance is computes
using `pftools`

commands as described in Manipulating Data: PFTools.
There are two water balance storage components, subsurface and surface,
and two flux calculations, overland flow and evapotranspiration.
The storage components have units [\(L^3\)] while the fluxes may be
instantaneous and have units [\(L^3T^{-1}\)] or cumulative over an
output interval with units [\(L^3\)]. Examples of water balance
calculations and errors are given in the scripts `water_balance_x.tcl`

and `water_balance_y.tcl`

. The size of water balance errors
depend on solver settings and tolerances but are typically very
small, \(<10^{-10}\) [-]. The water balance takes the form:

where \(Vol_{subsurface}\) is the subsurface storage [\(L^3\)]; \(Vol_{surface}\) is the
surface storage [\(L^3\)]; \(Q_{overland}\) is the overland flux [\(L^3 T^{-1}\)];
\(Q_{evapotranspiration}\) is the evapotranspiration flux passed
from `clm`

or other LSM, etc, [\(L^3 T^{-1}\)]; and
\(Q_{source sink}\) are any other source/sink fluxes specified in
the simulation [\(L^3 T^{-1}\)]. The surface and subsurface
storage routines are calculated using the ParFlow toolset commands `pfsurfacestorage`

and `pfsubsurfacestorage`

respectively. Overland flow out of the domain is calculated
by `pfsurfacerunoff`

. Details for the use of these commands are given in PFTCL Commands
and Common examples using ParFlow TCL commands (PFTCL). \(Q_{evapotranspiration}\) must be written out by ParFlow as a
variable (as shown in Code Parameters) and only contains the external fluxes passed
from a module such as `clm`

or WRF. Note that these volume and flux quantities are calculated
spatially over the domain and are returned as array values, just like any other quantity in ParFlow.
The tools command `pfsum`

will sum these arrays into a single value for the enrite domain.
All other fluxes must be determined by the user.

The subsurface storage is calculated over all active cells in the domain, \(\Omega\), and contains both compressible and incompressible parts based on Equation (4.3). This is computed on a cell-by-cell basis (with the result being an array of balances over the domain) as follows:

The surface storage is calculated over the upper surface boundary cells in the domain, \(\Gamma\), as computed by the mask and contains based on Equation [eq:kinematic]. This is again computed on a cell-by-cell basis (with the result being an array of balances over the domain) as follows:

For the overland flow outflow from the domain, any cell at the top boundary that has a slope that points out of the domain and is ponded will remove water from the domain. This is calculated, for example in the y-direction, as the multiple of Equation [eq:manningsy] and the area: