Groundwater flow

Single layer groundwater flow requires the four following components:

The aquifer types, constanthead boundary and other boundaries (e.g. river and drainage) are described in more detail below. Which functionality of these components is supported by model type sbm_gwf is also clarified below.

Aquifer types

Groundwater flow can occur either in a confined or unconfined aquifer. Confined aquifers are overlain by a poorly permeable confining layer (e.g. clay). No air can get in to fill the pore space so that the aquifer always remains fully saturated. For a confined aquifer, water will always flow along the complete height \(\SIb{H}{m}\) over the aquifer and transmissivity \(\SIb{kH}{m^2\ d^{-1}}\) is a constant (\(\SIb{k}{m\ d^{-1}}\) is the horizontal hydraulic conductivity). Specific storage is the amount of water an aquifer releases per unit change in hydraulic head, per unit volume of aquifer, as the aquifer and the groundwater itself is compressed. Its value is much smaller than specific yield, between \(10^{-5}\) (stiff) and \(10^{-2}\) (weak).

The upper boundary of an unconfined aquifer is the water table (the phreatic surface). Specific yield (or drainable porosity) represents the volumetric fraction the aquifer will yield when all water drains and the pore volume is filled by air instead. Specific yield will vary roughly between \(0.05\) (clay) and \(0.45\) (peat) (Johnson, 1967).

Groundwater flow is solved forward in time and central in space. The vertically averaged governing equation for an inhomogeneous and isotropic aquifer in one dimension can be written as:

\[ S \frac{\partial \phi}{\partial t} = \frac{\partial}{\partial x} \left(kH \frac{\phi}{\delta x}\right) + Q \]

where \(\SIb{S}{m\ m^{-1}}\) is storativity (or specific yield), \(\SIb{\phi}{m}\) is hydraulic head, \(t\) is time, \(\SIb{k}{m\ t^{-1}}\) is horizontal hydraulic conductivity, \(\SIb{H}{m}\) is the (saturated) aquifer height: groundwater level - aquifer bottom elevation and \(\SIb{Q}{m\ t^{-1}}\) represents fluxes from boundary conditions (e.g. recharge or abstraction), see also Aquifer boundary conditions.

The simplest finite difference formulation is forward in time, central in space, and can be written as:

\[ S_i \frac{\phi_{i}^{t+1} - \phi_i^{t}}{\Delta t} = -C_{i-1} (\phi_{i-1} - \phi_i) - C_i (\phi_{i+1} - \phi_i) + Q_i \]

where \(i\) is the cell index, \(t\) is time, \(\Delta t\) is the step size, \(C_{i-1}\) is the the intercell conductance between cell \(i-1\) and \(i\) and \(C_i\) is the intercell conductance between cell \(i\) and \(i+1\).

Conductance \(C\) is defined as:

\[ C = \frac{kH w}{l} \]

where \(\SIb{w}{m}\) is the width of the cell to cell connection, and \(\SIb{l}{m}\) is the length of the cell to cell connection. \(k\) and \(H\) may both vary in space; intercell conductance is therefore an average using the properties of two cells. For the calculation of the intercell conductance \(C\) the harmonic mean is used (see also Goode and Appel, 1992), here between cell index \(i\) and cell index \(i+1\), in the \(x\) direction:

\[ C_i = w \frac{k_iH_i\cdot k_{i+1}H_{i+1}}{k_iH_i \cdot l_{i+1} + k_{i+1}H_{i+1} \cdot l_i} \]

where \(\SIb{H}{m}\) is the aquifer top - aquifer bottom, and \(k\), \(l_i\) is the length in cell \(i\) (\(0.5 \Delta x_i\)), \(l_{i+1}\) is the length in cell \(i+1\) (\(0.5 \Delta x_{i+1}\)) and \(w\) as previously defined. For an unconfined aquifer the intercell conductance is scaled by using the “upstream saturated fraction” as the MODFLOW documentation calls it. In this approach, the saturated thickness of a cell-to-cell is approximated using the cell with the highest head. This results in a consistent overestimation of the saturated thickness, but it avoids complexities related with cell drying and rewetting, such as having to define a “wetting threshold” or a “wetting factor”. See also the documentation for MODFLOW-NWT (Niswonger et al., 2011) or MODFLOW6 (Langevin et al., 2017) for more background information. For more background on drying and rewetting, see for example McDonald et al. (1991).

For the finite difference formulation, there is only one unknown, \(\phi_i^{t+1}\). Reshuffling terms:

\[ \phi_i^{t+1} = \phi_i^t + (C_{i-1} (\phi_i - \phi_{i-1}) + C_i (\phi_{i+1} - \phi_i) + Q_i) \frac{Δt}{S_i} \]

This can be generalized to two dimensions, for both regular and irregular cell connectivity. Finally, a stable time step size can be computed given the forward-in-time, central in space scheme, based on the following criterion from Chu and Willis (1984):

\[ \frac{\Delta t k H}{\Delta x \Delta y S} \le \alpha \]

where \(\SIb{\Delta t}{d}\) is the stable time step size, \(\SIb{\Delta x}{m}\) is the cell length in the \(x\) direction and \(\SIb{\Delta y}{m}\) is the cell length in the \(y\) direction, \(\SIb{k}{m^2\ d^{-1}}\) is the horizontal hydraulic conductivity, \(\SIb{H}{m}\) is the saturated thickness of the aquifer and \(\alpha\) is a coefficient to enhance the stability of the simulation. For each cell \(\frac{\Delta x \Delta y S}{k H}\) is calculated, the minimum of these values is determined, and multiplied by \(\alpha\), to get the stable time step size. The stability coefficient \(\alpha\) (default = 0.25 as recommended by Chu and Willis (1984)) can be set as follows in the TOML file:

[model]
subsurface_water_flow__alpha_coefficient = 0.2 # coefficient for model stability (default = 0.25)

For more details about the finite difference formulation and the stable time step size criterion we refer to the paper of Chu and Willis (1984).

Note

Using the default value for the stability criterion from Chu and Willis (1984) can cause in some cases numerical instabilities, for example source and sink flux terms are not included. We recommend to check for possible numerical instabilities, especially when fluxes from aquifer boundaries (e.g. river) can be large.

Boundary conditions can be classified into three categories:

  • specified head (Dirichlet)
  • specified flux (Neumann)
  • head-dependent flux (Robin)

Neumann and Robin conditions are implemented by adding to or subtracting from a net (lumped) cell flux. Dirichlet conditions are special cased, since they cannot (easily) be implemented via the flux, but the head is set directly instead.

The groundwater flow component of model type sbm_gwf consists of a single layer unconfined aquifer. The list of input parameters for an unconfined aquifer can be found here. Hydraulic head \(\SIb{\phi}{m}\) is a state variable and output variable.

Below an example of setting the conductivity_profile and input parameters in the TOML configuration file for an unconfined aquifer:

[model]
conductivity_profile = "exponential"  # saturated hydraulic conductivity depth profile, default is "uniform"

[input.static]
land_surface__elevation = "wflow_dem"
subsurface_surface_water__horizontal_saturated_hydraulic_conductivity = "kh_surface"
subsurface_water__specific_yield = "specific_yield"
subsurface__horizontal_saturated_hydraulic_conductivity_scale_parameter = "gwf_f"

Constant head

For model type sbm_gwf Dirichlet boundary conditions can be specified (optional) in the TOML file as follows:

[model]
constanthead__flag = true # optional, default is "false"

[input.static]
"model_boundary_condition~constant_hydraulic_head" = "constant_head"

Aquifer boundary conditions

River

The flux between river and aquifer is calculated using Darcy’s law following the approach in MODFLOW:

\[ \subtext{Q}{riv} = \begin{align*} \begin{cases} C_i \min \left\{\subtext{h}{riv} - \subtext{B}{riv}, \subtext{h}{riv} - \phi\right\} &\text{ if }\quad \subtext{h}{riv} > \phi \\ C_e (\subtext{h}{riv} - \phi) &\text{ if }\quad \subtext{h}{riv} \le \phi \end{cases} \end{align*} \]

where \(\SIb{\subtext{Q}{riv}}{L^3 T^{-1}}\) is the exchange flux from river to aquifer, \(\SIb{C_i}{L^2 T^{-1}}\) is the river bed infiltration conductance, \(\SIb{C_e}{L^2 T^{-1}}\) is the river bed exfiltration conductance, \(\SIb{\subtext{B}{riv}}{L}\) the bottom of the river bed, \(\SIb{\subtext{h}{riv}}{L}\) is the river stage and \(\SIb{\phi}{L}\) is the hydraulic head in the river cell.

The list of input parameters for the river boundary of groundwater flow that can be provided through the TOML file can be found here.

The exchange flux (river to aquifer) \(\subtext{Q}{riv}\) is an output variable, and is used to update the total flux in a river cell. For model type sbm_gwf, the water depth \(\SIb{h}{m}\) of the river routing in combination with the river bottom is used to update the river stage each time step.

Drainage

The flux from drains to the aquifer is calculated as follows:

\[ \subtext{Q}{drain} = \subtext{C}{drain} \min(0, \subtext{h}{drain} - \phi) \]

where \(\SIb{\subtext{Q}{drain}}{L^3\ T^{-1}}\) is the exchange flux from drains to aquifer, \(\SIb{\subtext{C}{drain}}{L^2\ T^{-1}}\) is the drain conductance, \(\SIb{\subtext{h}{drain}}{L}\) is the drain elevation and \(\SIb{\phi}{L}\) is the hydraulic head in the cell with drainage.

The list of input parameters for the drainage boundary of groundwater flow that can be provided through the TOML file can be found here.

The exchange flux (drains to aquifer) \(\subtext{Q}{drain}\) is an output variable, and is used to update the total flux in a cell with drains. For model type sbm_gwf this boundary condition is optional, and can be specified in the TOML file as follows:

[model]
drain__flag = true # optional, default is "false"

[input.static]
land_drain_location__mask = "drain"
land_drain__conductance = "cond_drain"
land_drain__elevation = "elev_drain"

Recharge

The net recharge flux \(\subtext{Q}{net recharge}\) to the aquifer is calculated as follows:

\[ \subtext{Q}{net recharge} = \subtext{R}{net} \, A \]

with \(\SIb{\subtext{R}{net}}{L\ T^{-1}}\) the net recharge rate and \(\SIb{A}{L^2}\) the area of the aquifer cell.

The recharge flux \(\subtext{Q}{net recharge}\) is an output variable, and is used to update the total flux in a cell where net recharge occurs. For model type sbm_gwf, the recharge rate \(\SIb{}{mm\ t^{-1}}\) from the SBM soil model is used to update the recharge rate \(\subtext{R}{net}\) of the recharge boundary each time step.

Head boundary

This boundary is a fixed head with time (not affected by the model stresses over time) outside of the model domain, and is generally used to avoid an unnecessary extension of the model domain to the location of the fixed boundary (for example a large lake). The flux from the boundary \(\SIb{Q_{hb}}{L^3\ T^{-1}}\) is calculated as follows:

\[ Q_{hb} = C_{hb} (\phi_{hb} - \phi) \]

with \(\SIb{C_{hb}}{L^2\ T^{-1}}\) the conductance of the head boundary, \(\SIb{\phi_{hb}}{L}\) the head of the head boundary and \(\phi\) the head of the aquifer cell.

The head boundary flux \(Q_{hb}\) is an output variable, and is used to update the total flux in a cell where this type of boundary occurs. The parameter Head \(\phi_{hb}\) can be specified as a fixed or time dependent value.

Note

This boundary is not (yet) part of model type sbm_gwf.

Well boundary

A volumetric well rate \(\SIb{}{L^3\ T^{-1}}\) can be specified as a boundary condition.

The volumetric well rate \(\subtext{Q}{well}\) can be can be specified as a fixed or time dependent value. If a cell is dry, the actual well flux is set to zero (see also the last note on this page).

Note

This boundary is not (yet) part model type sbm_gwf.

Note

For an unconfined aquifer the boundary fluxes are checked, in case of a dry aquifer cell a negative flux is not allowed.

References

  • Chu, W. S., & Willis, R. (1984). An explicit finite difference model for unconfined aquifers. Groundwater, 22(6), 728-734.
  • Goode, D. J., & Appel, C. A. (1992). Finite-Difference Interblock Transmissivity for Unconfined Aquifers and for Aquifers having Smoothly Varying Transmissivity Water-resources investigations report, 92, 4124.
  • Johnson, A. I. (1967), Specific yield: compilation of specific yields for various materials, Water Supply Paper 1662-D, Washington, D.C.: U.S. Government Printing Office, p. 74, doi:10.3133/wsp1662D.
  • Langevin, C.D., Hughes, J.D., Banta, E.R., Niswonger, R.G., Panday, Sorab, and Provost, A.M., 2017, Documentation for the MODFLOW 6 Groundwater Flow Model: U.S. Geological Survey Techniques and Methods, book 6, chap. A55, 197 p., https://doi.org/10.3133/tm6A55.
  • McDonald, M.G., Harbaugh, A.W., Orr, B.R., and Ackerman, D.J., 1991, A method of converting no-flow cells to variable-head cells for the U.S. Geological Survey modular finite-difference groundwater flow model: U.S. Geological Survey Open-File Report 91-536, 99 p.
  • Niswonger, R.G., Panday, Sorab, and Ibaraki, Motomu, 2011, MODFLOW-NWT, A Newton formulation for MODFLOW-2005: U.S. Geological Survey Techniques and Methods 6-A37, 44 p.
Back to top