Groundwater flow

Single layer groundwater flow is defined in the struct GroundwaterFlow, which contains the following fields:

aquifer
connectivity
constanthead
boundaries

Base.@kwdef struct GroundwaterFlow
    aquifer::A where A <: Aquifer
    connectivity::Connectivity
    constanthead::ConstantHead
    boundaries::Vector{B} where B <: AquiferBoundaryCondition
end

The fields of struct GroundwaterFlow are described in more detail below.

Aquifer types

Groundwater flow contains the Abstract typeAquifer, that can be either a confined (ConfinedAquifer) or unconfined (UnconfinedAquifer) aquifer. Groundwater flow is solved forward in time and central in space.

Wflow.Aquifer — Type

Aquifer

Abstract type representing an aquifer, either confined or unconfined.

The vertically averaged governing equation of an unconfined, inhomogeneous and isotropic aquifer in one dimension can be written as:

S * ∂ϕ / ∂t = ∂ / ∂x * (kH * ∂ϕ / ∂x) + Q

with:

S: storativity (or storage coefficient)
ϕ: hydraulic head (groundwater level)
t: time
k: conductivity
H: H the (saturated) aquifer height: groundwater level - aquifer bottom elevation
η: elevation of aquifer bottom
Q: fluxes from boundary conditions (e.g. recharge or abstraction)

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

Sᵢ * (ϕᵢᵗ⁺¹ - ϕᵢᵗ) / Δt = -Cᵢ₋₁ * (ϕᵢ₋₁ - ϕᵢ) -Cᵢ * (ϕᵢ₊₁ - ϕᵢ) + Qᵢ

with:

ᵢ as cell index
ᵗ as time index
Δt as step size
Cᵢ₋₁ as the intercell conductance between cell i-1 and i
Cᵢ as the intercell conductance between cell i and i+1

Conductance is defined as:

C = kH * w / l

with:

w the width of the cell to cell connection
l 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. See the documentation below.

There is only one unknown, ϕᵢᵗ⁺¹. Reshuffling terms:

ϕᵢᵗ⁺¹ = ϕᵢᵗ + (Cᵢ₋₁ * (ϕᵢ - ϕᵢ₋₁) + Cᵢ * (ϕᵢ₊₁ - ϕᵢ) + Qᵢ) * Δt / Sᵢ

This can be generalized to two dimensions, for both regular and irregular cell connectivity.

See this paper for more details: Chu, W. S., & Willis, R. (1984). An explicit finite difference model for unconfined aquifers. Groundwater, 22(6), 728-734.

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.

source

Wflow.ConfinedAquifer — Type

ConfinedAquifer{T} <: 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 over the aquifer and transmissivity kH (m d⁻¹) is a constant.

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 1E-5 (stiff) and 0.01 (weak).

NOTA BENE: specific storage is per m of aquifer (conf. specific weight). Storativity or (storage coefficient) is for the entire aquifer (conf. transmissivity).

source

Wflow.UnconfinedAquifer — Type

UnconfinedAquifer{T} <: Aquifer

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).

source

Connectivity

The connectivity between cells is defined as follows.

Wflow.Connectivity — Type

Connectivity{T}

Stores connection data between cells. Connections are stored in a compressed sparse column (CSC) adjacency matrix: only non-zero values are stored. Primarily, this consist of two vectors:

the row value vector holds the cell indices of neighbors
the column pointers marks the start and end index of the row value vector

This matrix is square: n = m = ncell. nconnection is equal to nnz (the number of non-zero values).

ncell: the number of (active) cells in the simulation
nconnection: the number of connections between cells
length1: for every connection, the length in the first cell, size nconnection
length2: for every connection, the length in the second cell, size nconnection
width: width for every connection, (approx.) perpendicular to length1 and length2, size nconnection
colptr: CSC column pointer (size ncell + 1)
rowval: CSC row value (size nconnection)

source

Constant head

Dirichlet boundary conditions can be specified through the field constanthead.

Aquifer boundary conditions

River

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

\[ Q_{riv} = \Bigg\lbrace{C_{i} \,\text{min}(h_{riv} - B_{riv}, h_{riv} - \phi), \,h_{riv} > \phi \atop C_{e} (h_{riv} - \phi) , \,h_{riv} \leq \phi}\]

where $Q_{riv}$ is the exchange flux from river to aquifer [L$^3$ T$^{-1}$], $C_i$ [L$^2$ T$^{-1}$] is the river bed infiltration conductance, $C_e$ [L$^2$ T$^{-1}$] is the river bed exfiltration conductance, $B_{riv}$ the bottom of the river bed [L], $h_{riv}$ is the rive stage [L] and $\phi$ is the hydraulic head in the river cell [L].

Drainage

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

\[Q_{drain} = C_{drain} \text{min}(0, h_{drain} - \phi)\]

where $Q_{drain}$ is the exchange flux from drains to aquifer [L$^3$ T$^{-1}$], $C_{drain}$ [L$^2$ T$^{-1}$] is the drain conductance, $h_{drain}$ is the drain elevation [L] and $\phi$ is the hydraulic head in the cell with drainage [L].

Recharge

The recharge flux $Q_{r}$ to the aquifer is calculated as follows:

\[Q_{r} = R \, A\]

with $R$ the recharge rate [L T$^{-1}$] and $A$ the area [L$^2$ ] of the aquifer cell.

Head boundary

This boundary is a fixed head with 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. The flux from the boundary $Q_{hb}$ [L$^3$ T$^{-1}$] is calculated as follows:

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

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

Well boundary

A volumetric well rate [L$^3$ T$^{-1}$] can be specified as a boundary condition.

Note

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