Groundwater Flow

HydroGeoSphere assumes that the subsurface flow equation in a porous medium is always solved during a simulation, either for fully saturated or variably saturated flow conditions. The subsurface flow equation can be expanded to incorporate discrete fractures, a second interacting porous continuum, wells, tile drains and surface flow. The following assumptions are made for subsurface flow:

• The fluid is essentially incompressible.
• The porous medium and fractures (or macropores), if present, are non-deformable.
• The system is under isothermal conditions.
• The air phase is infinitely mobile.

The Richards’ equation is used to describe three-dimensional transient subsurface flow in a variably saturated porous medium:

${\displaystyle -\nabla \cdot (w_{m}{\textbf {q}})+\sum \Gamma _{eq}\pm Q=w_{m}{\frac {\partial }{\partial t}}(\theta _{s}S_{w})}$

The fluid flux, ${\displaystyle {\textbf {q}}}$, is represented by the Darcy's law shown as:

${\displaystyle {\textbf {q}}=-{\textbf {K}}\cdot k_{r}\nabla (\psi +z)}$

where ${\displaystyle w_{m}}$ is the volumetric fraction of the total porosity occupied by the porous medium, ${\displaystyle \Gamma _{ex}}$ is the internal fluid exchange rate (e.g. surface water, wells, and tile drains), ${\displaystyle Q}$ is the external fluid outside of the model domain, ${\displaystyle \theta _{s}}$ is the saturated water content, ${\displaystyle S_{w}}$ is the degree of saturation, ${\displaystyle {\textbf {K}}}$ is the hydraulic conductivity tensor, ${\displaystyle k_{r}}$ is the relative permeability of the medium calculated as a function of saturation, ${\displaystyle \psi }$ is the pressure head, and ${\displaystyle z}$ is the elevation head.

Surface water flow

Areal surface water flow is represented in HydroGeoSphere by a two-dimensional depth-averaged flow equation, which is the diffusion-wave approximation of the Saint Venant equation for surface water flow. HydroGeoSphere's surface water flow component is implemented with the following assumptions:

• Depth-averaged flow velocities
• Hydrostatic pressure distribution vertically
• Mild slope
• Dominant bottom shear stresses.

The surface flow components are solved by the following three equations, which are given by the following mass balance equation:

${\displaystyle {\frac {\partial \phi _{o}h_{o}}{\partial t}}+{\frac {\partial {\bar {v}}_{xo}d_{o}}{\partial x}}+{\frac {\partial {\bar {v}}_{yo}d_{o}}{\partial y}}+d_{o}\Gamma _{o}\pm Q_{o}=0}$

coupled with the momentum equations, neglecting inertia terms, for the x-direction:

${\displaystyle {\bar {v}}_{ox}=-K_{ox}{\frac {\partial h_{o}}{\partial x}}}$

and for the y-direction:

${\displaystyle {\bar {v}}_{oy}=-K_{oy}{\frac {\partial h_{o}}{\partial y}}}$

where ${\displaystyle \phi _{o}}$ is the surface flow domain porosity, ${\displaystyle h_{o}}$ is the water surface elevation, ${\displaystyle {\bar {v}}_{xo}}$ and ${\displaystyle {\bar {v}}_{yo}}$ are the vertically averaged flow velocities in the x and y directions, ${\displaystyle d_{o}}$ is the depth of surface water flow, ${\displaystyle \Gamma _{o}}$ is the internal fluid exchange, and ${\displaystyle Q_{o}}$ is the external fluid exchange. The surface conductances, ${\displaystyle K_{ox}}$ and ${\displaystyle K_{oy}}$ are approximated by either the Manning or Chezy equation.

Solute transport

Three-dimensional transport of solutes is described by the modified reactive transport advective-dispersion equation:

${\displaystyle -\nabla \cdot w_{m}({\textbf {q}}C-\theta _{s}S_{w}{\textbf {D}}\nabla C)+[w_{m}\theta _{s}S_{w}R\lambda C]_{par}+\sum {\Omega _{ex}+Q_{c}}=w_{m}\left[{\frac {\partial (\theta _{s}S_{w}RC)}{\partial t}}+\theta _{s}S_{w}R\lambda C\right]}$

where ${\displaystyle C}$ is the solute concentration, ${\displaystyle \lambda }$ is the first-order decay constant, ${\displaystyle Q_{c}}$ is the external source or sink term, ${\displaystyle \Omega }$ is the internal solute transfer between domains, ${\displaystyle R}$ is the retardation factor, ${\displaystyle {\textbf {D}}}$ is the diffusion coefficient, and ${\displaystyle par}$ designates parent species for the case of a decay chain.

Heat transport

Graf [2005] incorporated heat transport within the saturated-zone flow regime into HydroGeoSphere together with temperature-dependent fluid properties, such as viscosity and density. The model’s capability was successfully demonstrated for the case of thermohaline flow and transport in porous and fractured porous media [Graf and Therrien, 2007]. This work extends the model’s capability to include thermal energy transport in the unsaturated zone and in the surface water, which is considered a key step in the linkage between the atmospheric and hydrologic systems. Surface heat fluxes from atmospheric inputs are an important source/sink of thermal energy, especially to the surface water system. As such, surface heat fluxes across the land surface were also incorporated into HydroGeoSphere. A complete description of the physical processes and governing flow and solute transport equations that form the basis of HydroGeoSphere can be found in Therrien et al. [2007] and therefore will not be presented here.

The general equation for variably saturated subsurface thermal energy transport following Molson et al. [1992] is given by:

${\displaystyle -\nabla \cdot {\Big (}{\textbf {q}}\rho _{w}c_{w}T-(k_{b}+c_{b}\rho _{b}{\textbf {D}})\nabla T{\Big )}+\Omega _{o}\pm Q_{T}={\frac {\partial \rho _{b}c_{b}T}{\partial t}}}$

where ${\displaystyle \rho }$ is the density, ${\displaystyle c}$ is the heat capacity, ${\displaystyle T}$ is the temperature of the bulk subsurface, ${\displaystyle k}$ is the thermal conductivity, ${\displaystyle D}$ is the thermal dispersion term, ${\displaystyle Q_{T}}$ is the thermal source/sink, ${\displaystyle \Omega _{o}}$ is the thermal interactions between the surface and subsurface, and ${\displaystyle Q_{T}}$ is the external thermal interactions.