6.1 Governing Equations

Vincent E.A. Post; Craig T. Simmons

6.1 Governing Equations

For variable-density flow, the mass conservation equations for fluid (groundwater), mass and heat are the basis for numerical models. The mass balance equation for a compressible fluid in a saturated porous medium is shown in Equation 32.

$\displaystyle \frac{\partial (n\rho )}{\partial t}=-\nabla \cdot (\rho \vec{q})+\rho _{ss}q_{ss}$

(32)

where:

t	=	time (T), e.g., s
$\vec{q}$	=	specific discharge as given by Equation 2 (L/T), e.g., m s⁻¹
ρ_ss	=	density of water source, or sink (M/L³), e.g., kg m⁻³
q_ss	=	discharge rate per unit of volume of water source or sink (1/T), e.g., s⁻¹

The term on the left of the equal sign represents the change in fluid mass per unit of time. The first term to the right of the equal sign represents the change in fluid mass due to differences between inflow and outflow, while the second term is the change in fluid mass due to a source or sink. In groundwater, this can be recharge, inflow from an adjacent aquifer or an extraction/injection well, for example.

The mass balance equation for a solute that instantaneously and linearly partitions between a dissolved phase at concentration C in the pores and an adsorbed phase on the rock surface is Equation 33.

$\displaystyle \left ( 1+\frac{\rho _{b}K_{d}}{n} \right )\frac{\partial (nC)}{\partial t}=-\nabla\cdot (\vec{q}C)+\nabla \cdot (n\overline{D_{h,C}}\cdot \nabla C)+q_{s}C_{s}$

(33)

where:

ρ_b	=	dry bulk density of the rock (M/L³), e.g., kg m⁻³
K_d	=	distribution coefficient (L³/M), e.g., m³ kg⁻¹, which is the ratio of adsorbed mass and solute mass
$\overline{D_{h,C}}$	=	solute dispersion tensor (L²/T), e.g., m²s⁻¹, which includes the effect of both molecular diffusion and mechanical dispersion
C_s	=	concentration associated with a source (M/L³), e.g., kg m⁻³

The heat balance equivalent for Equation 33 for temperature being at local instantaneous equilibrium between the pore water and rock matrix is Equation 34.

$\displaystyle \left ( 1+\frac{1-n}{n}\frac{\rho _{r}c_{s}}{\rho c} \right )\frac{\partial (nT)}{\partial t}=-\nabla \cdot (\vec{q}T)+\nabla \cdot (\overline{D_{h,T}}\cdot \nabla T)+q_{s}T_{s}$

(34)

where:

ρ_r	=	density of the rock (M/L³), e.g., kg m⁻³
c and c_s	=	specific heat capacities of the groundwater and rock, respectively (ML²/(T²Θ)), e.g., J K⁻¹ kg⁻¹
$\overline{D_{h,T}}$	=	tensor which includes the combined effects of thermal conduction and mechanical dispersion (L²/T), e.g., m²s⁻¹
T_s	=	temperature associated with a source (Θ), e.g., K

The terms on the left of the equal sign in Equations 33 and 34 represent the change in solute mass or energy, respectively, per unit of time. For both equations, the terms on the right-hand side of the equal sign represent the change due to:

advective solute/heat transfer;
mechanical dispersion coupled with diffusion of a solute or conduction of heat; and,
fluid entering or exiting the system.

Most numerical codes assume that density is a linear function of C and T, as well as pressure p for the range of conditions within the modeled region such that all of the derivatives have constant values (Equation 35).

$\displaystyle \rho =\rho _{0}+\frac{\partial \rho }{\partial C}(C-C_{0})+\frac{\partial \rho }{\partial T}(T-T_{0})+\frac{\partial \rho }{\partial p}(p-p_{0})$

(35)

where:

ρ₀	=	density of the fluid at initial concentration, temperature and pressure (M/L³), e.g., kg m⁻³
C₀	=	initial concentration of the groundwater (M/L³), e.g., kg m⁻³
T₀	=	initial temperature of the groundwater (Θ), e.g., K
p₀	=	initial fluid pressure (M/(L¹T²)), e.g., kg m⁻¹ s⁻²

It should be noted that the slope $\frac{\partial \rho }{\partial C}$ is constant over a large concentration and temperature range, but $\frac{\partial \rho }{\partial T}$ is temperature-dependent over the range of temperatures of hydrogeological interest. The dependency of ρ on p is very small due to the low compressibility of water such that $\frac{\partial \rho }{\partial p}\ll \frac{\partial \rho }{\partial C}$ and $\frac{\partial \rho }{\partial p}\ll \frac{\partial \rho }{\partial T}$ . In fact, it is much smaller than the compressibility of most of the natural materials that make up aquifers, so it is usually not entered explicitly but combined with the elastic storage properties of the rock in the specific storage coefficient S_p. With the use of Equation 35, Equation 32 can be transformed into Equation 36 (Guo and Langevin, 2002; Jiao and Post, 2019).

$\displaystyle \rho S_{p}\frac{\partial p}{\partial t}+\frac{\partial \rho }{\partial C}\frac{\partial C}{\partial t}+\frac{\partial \rho }{\partial T}\frac{\partial T}{\partial t}=-\nabla \cdot (\rho \vec{q})+\rho _{ss}q_{ss}$

(36)

Equations 33, 34 and 36, or their equivalents, form the governing equations of numerical codes for variable-density groundwater flow. In addition to these, numerical models rely on the pressure-based form of Darcy’s law (Equation 2) to describe groundwater flow, Fick’s law (Equation 21) for dispersive and diffusive solute mass transfer and/or Fourier’s law for conductive heat transfer (Equation 22), as well as an equation of state like Equation 35. An inspection of the equations shows that they are coupled. That is, the solute concentration and temperature affect the density, the density affects the flow field, and the flow field alters the solute concentrations and temperatures. The governing equations for groundwater flow and solute transport, therefore, need to be solved simultaneously for variable-density flow problems. Codes for variable-density flow are therefore inherently more complex than those for constant-density systems.

License