GitHub

Soil hydrology

Overview

Soil hydrology processes characterize the dynamics of ground water in both saturated and unsaturated soil. It defines parameters and methods needed to compute water fluxes between layers and grid cells within the soil domain. Implementations should extend AbstractSoilHydrology and should generally consist of at least four components:

  • A scheme for computing vertical water fluxes between soil layers
  • A closure parameterization linking soil saturation and pressure head
  • A parameterization for soil hydraulic properties
  • A forcing term representing user-defined, internal sources/sinks in the soil domain (not including evapotranspiration)
Terrarium.AbstractSoilHydrology — Type
abstract type AbstractSoilHydrology{NF} <: Terrarium.AbstractProcess{NF}

Base type for soil hydrology implementations. Subtypes should define state variables for saturation_water_ice, hydraulic_conductivity, liquid_water_fraction, and the current water_table level, along with any other implementation-specific state variables.

source
subtypes(Terrarium.AbstractSoilHydrology)
1-element Vector{Any}:
 SoilHydrology

Terrarium currently provides a single general implementation of SoilHydrology following the above interface:

Terrarium.SoilHydrology — Type
struct SoilHydrology{NF, VerticalFlow<:Terrarium.AbstractVerticalFlow, SaturationClosure<:Terrarium.AbstractSoilWaterClosure, SoilHydraulics<:Terrarium.AbstractSoilHydraulics{NF}, VWCForcing<:Union{Nothing, Oceananigans.Forcings.ContinuousForcing{LX, LY, LZ, P} where {P, LX, LY, LZ}, Oceananigans.Forcings.DiscreteForcing}} <: Terrarium.AbstractSoilHydrology{NF}

Properties:

  • vertical_flow::Terrarium.AbstractVerticalFlow: Soil water vertical flow operator

  • closure::Terrarium.AbstractSoilWaterClosure: Closure relation for the soil hydrology state

  • hydraulic_properties::Terrarium.AbstractSoilHydraulics: Soil hydraulic properties parameterization

  • vwc_forcing::Union{Nothing, Oceananigans.Forcings.ContinuousForcing{LX, LY, LZ, P} where {P, LX, LY, LZ}, Oceananigans.Forcings.DiscreteForcing}: Forcing for soil moisture (volumetric water content)

source

Process interface

Dispatches for NoFlow:

Dispatches for RichardsEq:

Terrarium.initialize! — Method
initialize!(
    state,
    grid,
    hydrology::SoilHydrology{NF, RichardsEq, SaturationClosure, SoilHydraulics} where {SaturationClosure<:Terrarium.AbstractSoilWaterClosure, SoilHydraulics<:(Terrarium.AbstractSoilHydraulics{NF})},
    soil::Terrarium.AbstractSoil,
    constants::PhysicalConstants,
    args...
) where NF
source
Terrarium.compute_auxiliary! — Method
compute_auxiliary!(
    state,
    grid,
    hydrology::SoilHydrology{NF, RichardsEq, SaturationClosure, SoilHydraulics} where {SaturationClosure<:Terrarium.AbstractSoilWaterClosure, SoilHydraulics<:(Terrarium.AbstractSoilHydraulics{NF})},
    soil::Terrarium.AbstractSoil,
    args...
) where NF
source
Terrarium.compute_tendencies! — Method
compute_tendencies!(
    state,
    grid,
    hydrology::SoilHydrology{NF, RichardsEq, SaturationClosure, SoilHydraulics} where {SaturationClosure<:Terrarium.AbstractSoilWaterClosure, SoilHydraulics<:(Terrarium.AbstractSoilHydraulics{NF})},
    soil::Terrarium.AbstractSoil,
    constants::PhysicalConstants
) where NF
compute_tendencies!(
    state,
    grid,
    hydrology::SoilHydrology{NF, RichardsEq, SaturationClosure, SoilHydraulics} where {SaturationClosure<:Terrarium.AbstractSoilWaterClosure, SoilHydraulics<:(Terrarium.AbstractSoilHydraulics{NF})},
    soil::Terrarium.AbstractSoil,
    constants::PhysicalConstants,
    evtr::Union{Nothing, Terrarium.AbstractEvapotranspiration}
) where NF
compute_tendencies!(
    state,
    grid,
    hydrology::SoilHydrology{NF, RichardsEq, SaturationClosure, SoilHydraulics} where {SaturationClosure<:Terrarium.AbstractSoilWaterClosure, SoilHydraulics<:(Terrarium.AbstractSoilHydraulics{NF})},
    soil::Terrarium.AbstractSoil,
    constants::PhysicalConstants,
    evtr::Union{Nothing, Terrarium.AbstractEvapotranspiration},
    timestepping::Timestepping
) where NF

Compute the saturation tendency of the Richards-equation hydrology, with the evapotranspiration sink of evtr (if given) and the user-defined vwc_forcing as source terms. timestepping is the resolved Timestepping class of saturation_water_ice (see timestepping(state, ::Val)). Under Explicit the full divergence of the Darcy flux is added. Under Implicit the fluxes through the interior faces are omitted (integrated by the implicit sub-stepper, see BackwardEuler), leaving the fluxes through the two boundary faces (see compute_boundary_darcy_flux_divergence). The flux boundary conditions on saturation_water_ice (e.g. InfiltrationFlux) from compute_boundary_conditions! and the source terms are part of the tendency in both cases.

The Darcy fluxes at tⁿ through all faces are stored in water_flux (compute_water_flux), with which the soil energy advects heat. Under Implicit time stepping, the implicit step overwrites the interior fluxes with those of its conservative update (see implicit_update!).

source

Vertical flow

Static soil hydrology ("No Flow")

The simplest possible soil hydrology scheme is one in which the soil saturation state remains constant over time. This can be appropriate for simulations in regions where the soil is normally waterlogged or where changes in the hydrological state can be otherwise assumed to play a negligible role.

Terrarium.NoFlow — Type
struct NoFlow <: Terrarium.AbstractVerticalFlow

Represents a hydrology scheme where soil water is immobile.

source
variables(SoilHydrology(Float32, NoFlow()))
Variables
├─ Prognostic: 
├─ Auxiliary: 
├── saturation_water_ice [-] on XYZ{Center, Center, Center} of the ground domain
├── water_table [m] on XY{Center, Center, Nothing} of the ground domain
├── hydraulic_conductivity [m s^-1] on XYZ{Center, Center, Face} of the ground domain
├─ Inputs: 
├── liquid_water_fraction [-] on XYZ{Center, Center, Center} of the ground domain
├─ Namespaces:

Richardson-Richards equation for variably saturated flow

The vertical flow of water in porous media, such as soils, can be formulated as following the conservation law

\[ \phi\frac{\partial \xi(\psi)}{\partial t} = - \boldsymbol{\nabla} \cdot \textbf{j}_{\text{w}} + F_{\text{w}}(z,t),\]

where (as defined in Soil stratigraphy) $\phi$ is the natural porosity (or saturated water content) of the soil volume and $\xi \in [0,1]$ the saturation of pore water/ice. $F_{\text{w}}(z,t)$ is an inhomogeneous source/sink (forcing) term (1/s) and $\textbf{j}_{\text{w}}$ the water flux vector (m/s).

Vertical fluxes in the soil column can be represented by combining gravity-driven advection with Darcy's law

\[\begin{equation} \textbf{j}_{\text{w}}^{\mathrm{v}} =\textbf{j}_{\text{w}} \cdot \mathbf{n}_z = -\kappa_{\text{w}}\frac{\partial \left(\psi + z\right)}{\partial z}, \end{equation}\]

where $\psi$ is the matric potential (m), $\mathbf{n}_z$ the normal vector perpendicular to the surface and $\kappa_{\text{w}}$ the hydraulic conductivity (m/s). Given the positive upwards convention (see Numerical core), the $z-$axis is thus defined positive away from the surface. Substituting this equation into the aforementioned conservation law yields the widely known Richardson-Richards equation for variably saturated flow in porous media [7]. When the divergence of the water flux vector ($\boldsymbol{\nabla} \cdot \textbf{j}_{\text{w}}$) is simplified to only consider vertical water fluxes, this gives:

\[\begin{equation} \phi\frac{\partial \xi(\psi)}{\partial t} = \frac{\partial}{\partial z}\left[\kappa_w \frac{\partial \left(\psi + z\right)}{\partial z}\right] + F_{\text{w}}(z,t) \end{equation}\]

which is called the mixed-form equation, as the rate of change on the left-hand side is in $\xi$ but the vertical gradient on the right-hand side is in $\psi$ [2].

Terrarium.RichardsEq — Type
RichardsEq{PS} <: AbstractVerticalFlow

SoilHydrology flow operator implementing the mixed saturation-pressure form of the Richardson-Richards equation.

State variables defined by the Richards' formulation of SoilHydrology:

  • saturation_water_ice: saturation level of water and ice in the pore space.
  • hydraulic_conductivity: hydraulic conductivity at cell centers (m/s).
  • water_table`: elevation of the water table (m).
  • liquid_water_fraction: fraction of unfrozen liquid water in the pore space (dimensionless).
  • water_flux: Darcy flux (m/s, positive upward) through the cell faces that moved the water in the step (see compute_water_flux).

Excess water that reaches the soil surface is routed to the surface_excess_water pool defined by AbstractSurfaceRunoff process (see adjust_saturation_profile!), if passed as a dependency. For standalone soil hydrology with no surface runoff, surface excess water is discarded.

See also SoilSaturationPressureClosure and AbstractSoilHydraulics for details regarding the closure relating saturation and pressure head.

source
variables(SoilHydrology(Float32, RichardsEq()))
Variables
├─ Prognostic: 
├── saturation_water_ice [-] on XYZ{Center, Center, Center} of the ground domain
├─ Auxiliary: 
├── pressure_head [m] on XYZ{Center, Center, Center} of the ground domain
├── hydraulic_conductivity [m s^-1] on XYZ{Center, Center, Face} of the ground domain
├── water_table [m] on XY{Center, Center, Nothing} of the ground domain
├── water_flux [m s^-1] on XYZ{Center, Center, Face} of the ground domain
├─ Inputs: 
├── liquid_water_fraction [-] on XYZ{Center, Center, Center} of the ground domain
├─ Namespaces:

Implicit Richards equation

The Timestepping trait argument of compute_tendencies! for RichardsEq is the resolved timestepping class of saturation_water_ice (see Implicit-explicit time stepping). It is Explicit unless the model is run with an IMEX timestepper, in which case SoilModel and LandModel route the saturation of a Richards-equation hydrology to the implicit sub-stepper:

Terrarium.default_timestepping — Function
default_timestepping(hydrology::AbstractSoilHydrology)

Default Timestepping class of saturation_water_ice under an AbstractIMEX timestepper: Implicit for the Richards equation (SoilHydrology{NF, RichardsEq}) and Explicit otherwise.

julia> Terrarium.default_timestepping(SoilHydrology(Float64, RichardsEq()))
Implicit()

julia> Terrarium.default_timestepping(SoilHydrology(Float64, NoFlow()))
Explicit()
source
  • Under Explicit(), the tendency contains the full divergence of the darcy_flux, the evapotranspiration sink, the user-defined forcing, and the flux boundary conditions on saturation_water_ice (e.g. InfiltrationFlux).
  • Under Implicit(), the Darcy fluxes through the interior cell faces are omitted, since they are integrated by the implicit sub-stepper. Only the Darcy fluxes through the two boundary faces of the column remain, via compute_boundary_darcy_flux_divergence. These vanish for the default (no-flux) conditions on pressure_head and carry the prescribed flux for a gradient condition on pressure_head (e.g. FreeDrainage).

In both cases, the Darcy fluxes at $t^n$ through all cell faces are stored in the water_flux auxiliary variable (m/s, positive upward, on cell faces; see compute_water_flux). Under Implicit(), the implicit step overwrites the fluxes through the interior faces with those of its conservative update (compute_implicit_water_flux). After a step, water_flux therefore holds the flux that moved the water, with which the soil energy advects heat (see Soil energy balance).

The hydrology closure redistributes over- and undersaturation with adjust_saturation_profile!. It moves the sensible heat $C_\text{w} T$ of the water it shifts between cells in the internal_energy, at the temperature of the cell that the water leaves, and removes that of the excess water it discards (standalone SoilModel) or routes to the surface_excess_water pool (LandModel), which stores no heat. The energy closure, which runs after the hydrology closure, then recomputes the temperature from the adjusted energy.

The implicit step solves the mixed form of the Richards equation (see Richardson-Richards equation for variably saturated flow) with the mass-conservative modified Picard iteration of [3], using the lagged face conductivity of the Darcy flux and the lagged specific moisture capacity given by compute_implicit_conductivity! and compute_implicit_capacity!. With porosity $\phi$, saturation $\xi$ (so that the volumetric water content is $\theta = \phi \xi$), and total hydraulic head $\Psi = \psi_m(\theta) + \psi_z + \psi_h$ (matric, elevation, and hydrostatic head), the backward Euler step reads

\[\phi\, \frac{\xi^{n+1} - \xi^n}{\Delta t} = \frac{\partial}{\partial z}\left[K \frac{\partial \Psi^{n+1}}{\partial z}\right] + \phi\, G_\text{exp}\,,\]

where $G_\text{exp}$ (1/s) is the explicit tendency of the saturation $\xi$; it is multiplied by $\phi$ because the equation is written for the volumetric water content $\theta = \phi \xi$, the quantity whose flux divergence is the right-hand side. $G_\text{exp}$ contains the flux boundary conditions on the saturation (e.g. InfiltrationFlux), the Darcy fluxes through boundary faces implied by gradient conditions on the pressure head, evapotranspiration, and the user-defined forcing. The elevation head and the hydrostatic head below the water table are held fixed within the step, so that $\partial\theta/\partial\Psi = C$, the specific moisture capacity. Linearizing the water content about the latest iterate $m$, $\phi \xi^{m+1} \approx \phi \xi^m + C^m (\Psi^{m+1} - \Psi^m)$, with the hydraulic conductivity $K^m$ of the Darcy flux lagged at the same iterate, gives the tridiagonal system

\[C^m\, \Psi^{m+1} - \Delta t\, \frac{\partial}{\partial z}\left[K^m \frac{\partial \Psi^{m+1}}{\partial z}\right] = C^m\, \Psi^m - \phi\,(\xi^m - \xi^+)\,, \qquad \xi^+ = \xi^n + \Delta t\, G_\text{exp}\,,\]

for the total head, starting from $\xi^0 = \xi^n$. As for heat conduction (see Implicit vertical heat conduction), the saturation is then updated conservatively from the implicit fluxes through the interior faces,

\[\xi^{m+1} = \xi^+ + \frac{\Delta t}{\phi}\, \frac{\partial}{\partial z}\left[K^m \frac{\partial \Psi^{m+1}}{\partial z}\right],\]

and the hydrology closure is applied after every iteration: it redistributes over- and undersaturation (see adjust_saturation_profile!), updates the water table, and recomputes $\Psi^{m+1}$ from $\xi^{m+1}$. Excess water that reaches a saturated surface is routed into the surface_excess_water pool of a LandModel's surface runoff only by the closure of the last iteration (and is discarded by a standalone SoilModel), since every iteration restarts from $\xi^+$. For the same reason, only the closure of the last iteration moves the heat of the redistributed water (see Coupling of water and energy), so the internal energy is adjusted exactly once per step. The column water content therefore changes by exactly $\Delta t$ times the boundary fluxes and sources in $G_\text{exp}$, minus the excess routed to the surface pool, for any number of iterations.

In saturated cells $C = 0$, and the diagonal of the system is carried by the conductivity terms $\Delta t K / \Delta z^2$. Where the conductivity vanishes as well (frozen saturated cells), the system would be singular, so $C$ is bounded from below by the small minimum_moisture_capacity (1/m, default $10^{-5}$) of the hydraulic properties (see ConstantSoilHydraulics), which acts like a small specific storage within the iteration but does not affect mass conservation. [3] recommend at least two iterations, which is the default of Picard.

Flux boundary conditions on saturation_water_ice and gradient boundary conditions on pressure_head are supported; value boundary conditions on pressure_head are not yet supported.

The terms that stay explicit can remove more water from a cell in one step than it holds, e.g. gravity drainage through the bottom face of a draining sandy column. water_outflow_timescale reports the time in which they would drain a cell to its residual saturation, which bounds the step of a TimeStepWizard through cell_advection_timescale (see Adaptive time stepping with IMEX). It evaluates the explicit terms selected by the resolved class of saturation_water_ice, so under Explicit() it includes the interior Darcy fluxes as well.

The implicit operator of saturation_water_ice is a TridiagonalOperator in the total hydraulic head pressure_head, the variable of the SoilSaturationPressureClosure (see its varname method below), returned by implicit_operator (which also validates the configuration), and each iteration of the implicit solver calls the process hooks of SoilHydrology with RichardsEq below:

Terrarium.implicit_operator — Method
implicit_operator(
    hydrology::SoilHydrology{NF, RichardsEq, SaturationClosure, SoilHydraulics} where {SaturationClosure<:Terrarium.AbstractSoilWaterClosure, SoilHydraulics<:(Terrarium.AbstractSoilHydraulics{NF})},
    state,
    grid
) where NF

Allocate the TridiagonalOperator of saturation_water_ice on ground_domain(grid). Throws an ArgumentError unless the closure is the SoilSaturationPressureClosure and pressure_head has no value boundary conditions.

The operator is the backward Euler step of the Richards equation in the total pressure head Ψ, after the mass-conservative modified Picard iteration of [3]: with porosity φ, saturation ξ, and ξ⁺ = ξⁿ + Δt G_exp from the explicit stage, the water content is linearized about the iterate m, φ ξᵐ⁺¹ ≈ φ ξᵐ + Cᵐ (Ψᵐ⁺¹ - Ψᵐ) with the specific moisture capacity Cᵐ and the hydraulic conductivity Kᵐ lagged, giving the tridiagonal system Cᵐ Ψᵐ⁺¹ - Δt ∂z(Kᵐ ∂z Ψᵐ⁺¹) = Cᵐ Ψᵐ - φ (ξᵐ - ξ⁺) (see the soil hydrology docs for the derivation). Each Picard iteration m = 0, …, M-1, starting from ξ⁰ = ξⁿ, calls the hooks:

  1. linearize!: hydraulic conductivity, lagged Kᵐ (faces) and Cᵐ (centers, bounded below by minimum_moisture_capacity so that frozen saturated cells stay regular), and the right-hand side;
  2. solve!: the system above for Ψᵐ⁺¹;
  3. implicit_update!: ξᵐ⁺¹ = ξ⁺ + Δt ∂z(Kᵐ ∂z Ψᵐ⁺¹) / φ through the interior faces, which conserves water exactly, and water_flux ← -Kᵐ ∂z Ψᵐ⁺¹, the flux with which the heat step advects heat;
  4. implicit_closure!: redistribution of over- and undersaturation, water table, Ψᵐ⁺¹, and hydraulic conductivity. Since every iteration restarts from ξ⁺, only the final one routes surface excess water to the runoff pool and moves the heat of redistributed water.

The solver calls the hooks through the SoilModel/LandModel forwarding methods and the SoilEnergyWaterCarbon methods in soil_coupled.jl.

References

  • [3] Celia, Bouloutas, and Zarba, Water Resources Research 26, 1483–1496 (1990)
source
Terrarium.linearize! — Method
linearize!(
    operator::TridiagonalOperator,
    rhs,
    state,
    grid,
    hydrology::SoilHydrology{NF, RichardsEq, SaturationClosure, SoilHydraulics} where {SaturationClosure<:Terrarium.AbstractSoilWaterClosure, SoilHydraulics<:(Terrarium.AbstractSoilHydraulics{NF})},
    soil::Terrarium.AbstractSoil,
    constants::PhysicalConstants,
    explicit_saturation,
    Δt
) where NF

Update the hydraulic conductivity (compute_hydraulics!), fill the lagged capacity and face conductivity of operator (compute_implicit_water_flow_coefficients!), and assemble the right-hand side C Ψᵐ - φ (ξᵐ - ξ⁺) in rhs (compute_implicit_water_flow_rhs), with ξ⁺ = explicit_saturation. [3] recommend at least two iterations.

References

  • [3] Celia et al., Water Resources Research 26, 1483–1496 (1990)
source
Terrarium.implicit_update! — Method
implicit_update!(
    state,
    grid,
    hydrology::SoilHydrology{NF, RichardsEq, SaturationClosure, SoilHydraulics} where {SaturationClosure<:Terrarium.AbstractSoilWaterClosure, SoilHydraulics<:(Terrarium.AbstractSoilHydraulics{NF})},
    soil::Terrarium.AbstractSoil,
    operator::TridiagonalOperator,
    explicit_saturation,
    Δt
) where NF

Update saturation_water_ice conservatively from the solved pressure head (compute_implicit_saturation_update) and store the interior Darcy fluxes in water_flux (compute_implicit_water_flux).

source
Terrarium.implicit_closure! — Method
implicit_closure!(
    state,
    grid,
    hydrology::SoilHydrology{NF, RichardsEq, SaturationClosure, SoilHydraulics} where {SaturationClosure<:Terrarium.AbstractSoilWaterClosure, SoilHydraulics<:(Terrarium.AbstractSoilHydraulics{NF})},
    soil::Terrarium.AbstractSoil,
    constants::PhysicalConstants,
    runoff::Union{Nothing, Terrarium.AbstractSurfaceRunoff},
    ::Terrarium.Iteration{name, target, true}
) where {NF, name, target}

Apply the regular saturation–pressure closure! after the final iteration of the implicit step, which routes the surface excess to runoff and moves the heat of the redistributed water, then update the hydraulic conductivity.

source
Terrarium.implicit_closure! — Method
implicit_closure!(
    state,
    grid,
    hydrology::SoilHydrology{NF, RichardsEq, SaturationClosure, SoilHydraulics} where {SaturationClosure<:Terrarium.AbstractSoilWaterClosure, SoilHydraulics<:(Terrarium.AbstractSoilHydraulics{NF})},
    soil::Terrarium.AbstractSoil,
    constants::PhysicalConstants,
    runoff::Union{Nothing, Terrarium.AbstractSurfaceRunoff},
    ::Terrarium.Iteration{name, target, false}
) where {NF, name, target}

Apply the saturation–pressure closure without its once-per-step effects after an intermediate iteration of the implicit step (compute_saturation_closure! without runoff routing or heat transfer), since every iteration restarts from the explicit stage, then update the hydraulic conductivity.

source
Terrarium.compute_saturation_closure! — Function
compute_saturation_closure!(
    state,
    grid,
    closure::SoilSaturationPressureClosure,
    hydrology::SoilHydrology{NF, RichardsEq, SaturationClosure, SoilHydraulics} where {SaturationClosure<:Terrarium.AbstractSoilWaterClosure, SoilHydraulics<:(Terrarium.AbstractSoilHydraulics{NF})},
    soil::Terrarium.AbstractSoil,
    runoff::Union{Nothing, Terrarium.AbstractSurfaceRunoff},
    energy::Union{Nothing, Terrarium.AbstractSoilThermodynamics}
) where NF

Saturation–pressure closure of a Richards-equation hydrology: redistribute over- and undersaturation (adjust_saturation_profile!), routing surface excess water to runoff and moving the heat of the redistributed water in the internal_energy of energy (either may be nothing to skip the routing or the heat transfer), then update the water table and compute the pressure_head. Shared by closure! and the implicit Richards step, whose intermediate iterations skip the routing and the heat transfer.

source

Hydraulic properties

Match the conductivity to the retention curve

The unsaturated hydraulic conductivity and the soil water retention curve should come from the same family. By default, ConstantSoilHydraulics and SoilHydraulicsSURFEX select the conductivity that matches the retention curve (see Terrarium.default_unsat_hydraulic_cond). In particular, UnsatKLinear is not suitable for dry soil with a Brooks–Corey retention curve: the matric potential grows as a power law of the inverse saturation while the conductivity falls only linearly, so that the suction fluxes in dry soil become unphysically large.

Besides the static soil characteristics mentioned in Soil stratigraphy such as porosity ($\phi$), a key time-varying hydraulic property is the (non-saturated) hydraulic conductivity $\kappa_{\text{w}}$. Different ways of modelling $\kappa_{\text{w}}$ are defined as subyptes of AbstractUnsatK. The simplest formulation, called UnsatKLinear, models $\kappa_{\text{w}}$ as:

\[\begin{equation} \kappa_{\text{w}} = \frac{\theta_{\text{liq}}}{\theta + \theta_{\text{air}}} \kappa_{\text{w,sat}}\, \end{equation}\]

with $\kappa_{\text{w,sat}}$ the saturated hydraulic conductivity (m/s). A more complex formulation taking ice into account, called UnsatKVanGenuchten,follows [8], here presented in the form of [9, Eq. (26)]:

\[\begin{equation} \kappa_{\text{w}} = I_\text{ice}\left(\frac{\theta_{\text{liq}}}{\phi}\right)^{0.5} \left[1 - \left(1 - \left(\frac{\theta_{\text{liq}}}{\phi}\right)^{n/(n+1)}\right)^{(n-1)/n}\right]^2 \kappa_{\text{w,sat}}\,, \end{equation}\]

which adds an ice impedance factor $I_\text{ice} = 10^{-\Omega\theta_{\text{ice}}/\theta}$ to the classic van Genuchten formulation to account for the blocking of water-filled pores by ice.

Terrarium.UnsatKLinear — Type
struct UnsatKLinear{NF} <: Terrarium.AbstractUnsatK{NF}

Simple formulation of hydraulic conductivity as a linear function of the liquid water saturated fraction, i.e. soil.water / (soil.water + soil.ice + soil.air).

source
Terrarium.UnsatKVanGenuchten — Type
struct UnsatKVanGenuchten{NF} <: Terrarium.AbstractUnsatK{NF}

Formulation of hydraulic conductivity as a function of saturated hydraulic conductivity K_sat and volumetric fractions, assumed to include those of water, ice, and air, following the van Genuchten formulation [8] extended with an ice impedance factor [9].

References

  • [8] Van Genuchten, Soil Science Society of America Journal (1980)
  • [9] Westermann et al., Geoscientific Model Development (2023)
source

The conductivity of the Brooks–Corey model, called UnsatKBrooksCorey, is a power law of the liquid water saturation, with the same ice impedance factor:

\[\begin{equation} \kappa_{\text{w}} = I_\text{ice}\left(\frac{\theta_{\text{liq}}}{\phi}\right)^{3 + 2/\lambda} \kappa_{\text{w,sat}}\,, \end{equation}\]

where $\lambda$ is the pore size distribution index of the retention curve. [10] derived the exponent from their retention curve using the pore size distribution model of Burdine; [11] obtained the same law, written with $b = 1/\lambda$ as $(\theta_{\text{liq}}/\phi)^{2b + 3}$, and [12] tabulated $b$ by soil texture. In this form, it is the conductivity used in SURFEX/ISBA [13].

Terrarium.UnsatKBrooksCorey — Type
struct UnsatKBrooksCorey{NF} <: Terrarium.AbstractUnsatK{NF}

Unsaturated hydraulic conductivity of the Brooks–Corey model, extended with an ice impedance factor:

K = I_ice K_sat S^(3 + 2/λ)

where S = θw / θsat is the liquid water saturation, λ is the pore size distribution index of the BrooksCorey retention curve, and I_ice is the Terrarium.ice_impedance_factor of [9].

Brooks and Corey [10] define the power law on the effective saturation (θw − θres) / (θsat − θres). Terrarium evaluates the retention curve with the residual water content of the BrooksCorey curve itself, which is zero by default, and in that case the effective saturation reduces to S. A nonzero residual water content of the retention curve is not yet accounted for here; the residual_saturation of the hydraulic properties is only a numerical lower bound on the saturation.

The power law and its exponent 3 + 2/λ were derived by [10] from their retention curve, using the pore size distribution model of Burdine. [11] obtained the same law, written as S^(2b + 3) with b = 1/λ, and [12] tabulated b by soil texture; in this form, it is the conductivity of SURFEX/ISBA [13]. The exponent λ is taken from the retention curve so that the two cannot be configured inconsistently; this formulation can therefore only be combined with a BrooksCorey retention curve.

References

  • [10] Brooks & Corey, Hydrology Papers, Colorado State University (1964)
  • [11] Campbell, Soil Science (1974)
  • [12] Clapp & Hornberger, Water Resources Research (1978)
  • [13] Boone et al., Journal of Applied Meteorology (2000)
  • [9] Westermann et al., Geoscientific Model Development (2023)
source
Terrarium.default_unsat_hydraulic_cond — Function
default_unsat_hydraulic_cond(
    ::Type{NF},
    ::SWRC
) where NF -> UnsatKLinear

Return the unsaturated hydraulic conductivity formulation that matches the soil water retention curve swrc: UnsatKBrooksCorey for BrooksCorey, UnsatKVanGenuchten for VanGenuchten, and UnsatKLinear otherwise.

julia> Terrarium.default_unsat_hydraulic_cond(Float32, BrooksCorey())
UnsatKBrooksCorey{Float32}(7.0f0)

julia> Terrarium.default_unsat_hydraulic_cond(Float32, VanGenuchten())
UnsatKVanGenuchten{Float32}(7.0f0)
source

Values for $\kappa_{\text{w,sat}}$, wilting point ($\theta_{wp}$), and field capacity ($\theta_{fc}$) can be set as constants (see ConstantSoilHydraulics) or can be calculated from soil texture using pedotransfer functions (PTFs). Currently, PTFs based on [6] for $\theta_{wp}$ and $\theta_{fc}$ are available in SoilHydraulicsSURFEX.

Terrarium.ConstantSoilHydraulics — Type
struct ConstantSoilHydraulics{NF, RC, UnsatK<:Terrarium.AbstractUnsatK{NF}} <: Terrarium.AbstractSoilHydraulics{NF, RC, UnsatK}

Represents a simple case where soil hydraulic properties are given as constant values. This is mostly provided just for testing, although it may be useful in certain cases where direct measurements of hydraulic properites are available.

Properties:

  • swrc::Any: Soil water retention curve

  • unsat_hydraulic_cond::Terrarium.AbstractUnsatK: Unsaturated hydraulic conductivity formulation; defaults to the formulation matching swrc (see default_unsat_hydraulic_cond)

  • saturated_conductivity::Any: Hydraulic conductivity at saturation

  • field_capacity::Any: Constant field capacity

  • wilting_point::Any: Constant wilting point

  • residual_saturation::Any: Residual (minimum) saturation level

  • minimum_moisture_capacity::Any = 1.0e-5: Lower bound (1/m) of the specific moisture capacity ∂θ/∂ψ on the diagonal of the implicit Richards solve (implicit time stepping only)

source
Terrarium.SoilHydraulicsSURFEX — Type
struct SoilHydraulicsSURFEX{NF, RC, UnsatK<:Terrarium.AbstractUnsatK{NF}} <: Terrarium.AbstractSoilHydraulics{NF, RC, UnsatK}

Soil hydraulics parameterization that includes the SURFEX [6, Eq. (28-29)] formulation of field capacity and wilting point as a function of soil texture.

Properties:

  • swrc::Any: Soil water retention curve

  • unsat_hydraulic_cond::Terrarium.AbstractUnsatK: Unsaturated hydraulic conductivity formulation; defaults to the formulation matching swrc (see default_unsat_hydraulic_cond)

  • saturated_conductivity::Any: Hydraulic conductivity at saturation

  • wilting_point_effect::Any: Linear coefficient of wilting point adjustment due to clay content

  • field_capacity_effect::Any: Linear coefficient of field capacity adjustment due to clay content

  • field_capacity_exp::Any: Exponent of field capacity adjustment due to clay content

  • field_capacity_min::Any: Minimum field capacity at 0% clay content

  • residual_saturation::Any: Residual (minimum) saturation level

  • minimum_moisture_capacity::Any = 1.0e-5: Lower bound (1/m) of the specific moisture capacity ∂θ/∂ψ on the diagonal of the implicit Richards solve (implicit time stepping only)

References

  • [6] Noilhan & Mahfouf, Global and Planetary Change (1996)
source

Closures

Critical to the modelling of water transport in the soil, is the relationship between soil water/ice saturation ($\xi$) and the matric potential ($\Psi_m$), which is called the soil-water retention curve (SWRC). Recall from the documentation on the core interfaces that the total hydraulic head $\Psi$ is defined as:

\[\begin{equation} \Psi = \Psi_m(\xi) + \Psi_z + \Psi_h \end{equation}\]

with $\Psi_z$ the elevation head and $\Psi_h$ the hydrostatic head contributed by free water above the water table. Currently, two implementations of the SWRC are available in Terrarium via FreezeCurves.jl: the Brooks-Corey (BrooksCorey) and van Genuchten (VanGenuchten) formulations.

Terrarium.SoilSaturationPressureClosure — Type
struct SoilSaturationPressureClosure <: Terrarium.AbstractSoilWaterClosure

Represents a closure relating saturation of water/ice in soil pores to a corresponding pressure (or hydraulic) head. Note that here "pressure head" is defined to be synonymous with hydraulic head, i.e. including all both elevation and hydrostatic pressure contributions. This relation is typically described by soil property-dependent soil-water retention curve (SWRC) which is here defined in implementations of AbstractSoilHydraulics.

source

Methods

Terrarium.get_swrc — Function
get_swrc(::AbstractUnsatK)

Return the soil water retention curve associated with the given unsaturated hydraulic conductivity scheme.

source
Terrarium.get_closure — Function
get_closure(
    hydrology::SoilHydrology
) -> Terrarium.AbstractSoilWaterClosure

Return the saturation-pressure closure defined by the given hydrology process, or nothing if not defined for the given configuration.

source
get_closure(snow::SingleLayerSnow)
source
Terrarium.compute_water_table! — Function
compute_water_table!(
    water_table,
    i,
    j,
    grid,
    sat,
    ::SoilHydrology{NF, VerticalFlow, SaturationClosure, SoilHydraulics} where {VerticalFlow<:Terrarium.AbstractVerticalFlow, SaturationClosure<:Terrarium.AbstractSoilWaterClosure, SoilHydraulics<:(Terrarium.AbstractSoilHydraulics{NF})}
) where NF

Kernel function that diagnoses the water table at grid cell i, j given the current soil saturation profile.

source
Terrarium.adjust_saturation_profile! — Function
adjust_saturation_profile!(state, grid, hydrology::SoilHydrology, runoff = nothing, strat = nothing, bgc = nothing, energy = nothing)

Redistribute over- and undersaturation in the soil columns of state (see redistribute_saturation_profile!). With runoff = nothing, excess water reaching the surface is discarded; otherwise it is routed into the surface_excess_water pool of runoff as a water depth (m), using the porosity of the topmost layer from strat and bgc.

With an optional SoilThermodynamics energy process, the sensible heat of the water moves with it (compute_redistributed_heat!), so the heat of the excess water leaves the column; this requires strat and bgc.

source
adjust_saturation_profile!(
    out,
    i,
    j,
    grid,
    fields,
    hydrology::SoilHydrology{NF, VerticalFlow, SaturationClosure, SoilHydraulics} where {VerticalFlow<:Terrarium.AbstractVerticalFlow, SaturationClosure<:Terrarium.AbstractSoilWaterClosure, SoilHydraulics<:(Terrarium.AbstractSoilHydraulics{NF})},
    ::Nothing,
    energy,
    strat,
    bgc
) where NF

Kernel function that adjusts the saturation profile at i, j and discards any excess water reaching the surface. Used for standalone soil hydrology, where no surface runoff process owns a water pool.

source
adjust_saturation_profile!(
    out,
    i,
    j,
    grid,
    fields,
    hydrology::SoilHydrology{NF, VerticalFlow, SaturationClosure, SoilHydraulics} where {VerticalFlow<:Terrarium.AbstractVerticalFlow, SaturationClosure<:Terrarium.AbstractSoilWaterClosure, SoilHydraulics<:(Terrarium.AbstractSoilHydraulics{NF})},
    runoff::Terrarium.AbstractSurfaceRunoff,
    energy,
    strat::Terrarium.AbstractStratigraphy,
    bgc::Terrarium.AbstractSoilBiogeochemistry
) where NF

Kernel function that adjusts the saturation profile at i, j and routes any excess water reaching the surface into the surface_excess_water pool owned by the given surface runoff process, as a water depth (m): the excess saturation times the thickness and porosity (porosity_top) of the topmost layer.

source
Terrarium.compute_hydraulics! — Function
compute_hydraulics!(state, grid, hydrology::SoilHydrology, soil::AbstractSoil, args...)

Compute all state-dependent hydraulic auxiliaries such as hydraulic conductivity and field capacity, and wilting point.

source
Terrarium.minimum_moisture_capacity — Function
minimum_moisture_capacity(hydraulics::AbstractSoilHydraulics)

Numerical lower bound (1/m) of the specific moisture capacity ∂θ/∂ψ used on the diagonal of the implicit Richards solve (see compute_implicit_capacity!). It plays the role of a small specific storage: in saturated cells ∂θ/∂ψ = 0, and a saturated cell whose faces have zero hydraulic conductivity (e.g. a frozen saturated cell) would otherwise give a singular row of the tridiagonal system.

source
Terrarium.water_outflow_timescale — Function
water_outflow_timescale(
    state,
    grid,
    hydrology::SoilHydrology{NF, RichardsEq, SaturationClosure, SoilHydraulics} where {SaturationClosure<:Terrarium.AbstractSoilWaterClosure, SoilHydraulics<:(Terrarium.AbstractSoilHydraulics{NF})},
    soil::Terrarium.AbstractSoil,
    constants::PhysicalConstants,
    evtr::Union{Nothing, Terrarium.AbstractEvapotranspiration},
    timestepping::Timestepping
) where NF

Return the minimum over all cells of the explicit water outflow timescale (compute_water_outflow_timescale) of a Richards-equation hydrology: the time in which the explicit terms of the saturation tendency, evaluated at the current state, would drain a cell to its residual saturation. timestepping is the resolved Timestepping class of saturation_water_ice, which selects the explicit terms as in compute_tendencies!.

source

Kernel functions

Terrarium.saturation_water_ice — Function
saturation_water_ice(i, j, k, grid, fields, ::AbstractSoilHydrology)

Compute or retrieve the current saturation level of water + ice in the pore space.

source
Terrarium.hydraulic_conductivity — Function

Compute (variably saturated) hydraulic conductivity based on the given hydraulic properties, soil water retention curve (SWRC), and volumetric fractions.

source
Terrarium.ice_impedance_factor — Function
ice_impedance_factor(Ω, f)

Ice impedance factor I_ice = 10^(-Ω (1 - f)) [9], which reduces the hydraulic conductivity of partially frozen soil to account for the blocking of water-filled pores by ice, given the scaling parameter Ω and the unfrozen fraction f of the pore water.

References

  • [9] Westermann et al., Geoscientific Model Development (2023)
source
Terrarium.liquid_water_fraction — Function
liquid_water_fraction(i, j, k, grid, fields, ::AbstractSoilHydrology)

Compute or retrieve the current fraction of unfrozen water in the pore space.

source
Terrarium.water_table — Function
water_table(i, j, k, grid, fields, ::AbstractSoilHydrology)

Compute or retrieve the current water table level relative to the surface.

source
Terrarium.surface_excess_water — Function
surface_excess_water(i, j, k, grid, fields, ::AbstractSoilHydrology)

Retrieve the current saturation level of water + ice in the pore space.

source
Terrarium.compute_volumetric_water_content_tendency — Function
compute_volumetric_water_content_tendency(
    i,
    j,
    k,
    grid,
    clock,
    fields,
    hydrology::SoilHydrology{NF, VerticalFlow, SaturationClosure, SoilHydraulics} where {VerticalFlow<:Terrarium.AbstractVerticalFlow, SaturationClosure<:Terrarium.AbstractSoilWaterClosure, SoilHydraulics<:(Terrarium.AbstractSoilHydraulics{NF})},
    constants::PhysicalConstants,
    evtr::Union{Nothing, Terrarium.AbstractEvapotranspiration}
) where NF
compute_volumetric_water_content_tendency(
    i,
    j,
    k,
    grid,
    clock,
    fields,
    hydrology::SoilHydrology{NF, VerticalFlow, SaturationClosure, SoilHydraulics} where {VerticalFlow<:Terrarium.AbstractVerticalFlow, SaturationClosure<:Terrarium.AbstractSoilWaterClosure, SoilHydraulics<:(Terrarium.AbstractSoilHydraulics{NF})},
    constants::PhysicalConstants,
    evtr::Union{Nothing, Terrarium.AbstractEvapotranspiration},
    ::Timestepping
) where NF

Compute the volumetric water content (VWC) tendency at grid cell i, j k f. Note that the VWC tendency is not scaled by the porosity and is thus not the same as the saturation tendency.

source
compute_volumetric_water_content_tendency(
    i,
    j,
    k,
    grid,
    clock,
    fields,
    hydrology,
    constants,
    evapotranspiration
)
compute_volumetric_water_content_tendency(
    i,
    j,
    k,
    grid,
    clock,
    fields,
    hydrology,
    constants,
    evapotranspiration,
    timestepping
)

Compute the volumetric water content (VWC) tendency at grid cell i, j k according to the Richardson-Richards equation. Note that the VWC tendency is not scaled by the porosity and is thus not the same as the saturation tendency. Like the energy tendency, it is the explicitly integrated part of the tendency, and timestepping is the class of saturation_water_ice: under Explicit (the default) the full divergence of the Darcy flux is included; under Implicit only the boundary-face fluxes are (compute_boundary_darcy_flux_divergence), the interior ones being integrated by the implicit sub-stepper.

source
Terrarium.compute_boundary_darcy_flux_divergence — Function
compute_boundary_darcy_flux_divergence(
    i,
    j,
    k,
    grid,
    fields
)

Part of the Darcy flux divergence at i, j, k due to the two boundary faces of the column (bottom face of cell k = 1, top face of cell k = Nz), omitting all interior faces:

(δ_{k,Nz} j_{w,k+½} - δ_{k,1} j_{w,k-½}) / Δz_k,

with j_w the darcy_flux. This is the explicit part when the interior fluxes are integrated implicitly (see BackwardEuler). The boundary fluxes are computed from the pressure_head halo: they vanish for the default (no-flux) conditions and carry the prescribed flux for a gradient condition on pressure_head. Flux boundary conditions on saturation_water_ice (e.g. InfiltrationFlux) are added separately by compute_boundary_conditions!.

source
Terrarium.darcy_flux — Function
darcy_flux(i, j, k, grid, ψ, K)

Kernel function for computing the Darcy flux q = -K ∂ψ/∂z (m/s, positive upward) over layer faces from the total pressure head ψ and the face hydraulic conductivity K (see darcy_conductivity).

source
Terrarium.darcy_conductivity — Function
darcy_conductivity(i, j, k, grid, ψ, K)

Hydraulic conductivity (m/s) used by the darcy_flux at face k: the minimum of K (defined on faces) at face k and at the adjacent face in the direction of the gradient of the total pressure head ψ (face k - 1 for ∂ψ/∂z < 0, else k + 1).

source
Terrarium.redistribute_saturation_profile! — Function
redistribute_saturation_profile!(
    out,
    i,
    j,
    grid,
    fields,
    hydrology::SoilHydrology{NF, VerticalFlow, SaturationClosure, SoilHydraulics} where {VerticalFlow<:Terrarium.AbstractVerticalFlow, SaturationClosure<:Terrarium.AbstractSoilWaterClosure, SoilHydraulics<:(Terrarium.AbstractSoilHydraulics{NF})}
) where NF
redistribute_saturation_profile!(
    out,
    i,
    j,
    grid,
    fields,
    hydrology::SoilHydrology{NF, VerticalFlow, SaturationClosure, SoilHydraulics} where {VerticalFlow<:Terrarium.AbstractVerticalFlow, SaturationClosure<:Terrarium.AbstractSoilWaterClosure, SoilHydraulics<:(Terrarium.AbstractSoilHydraulics{NF})},
    energy
) where NF
redistribute_saturation_profile!(
    out,
    i,
    j,
    grid,
    fields,
    hydrology::SoilHydrology{NF, VerticalFlow, SaturationClosure, SoilHydraulics} where {VerticalFlow<:Terrarium.AbstractVerticalFlow, SaturationClosure<:Terrarium.AbstractSoilWaterClosure, SoilHydraulics<:(Terrarium.AbstractSoilHydraulics{NF})},
    energy,
    strat
) where NF
redistribute_saturation_profile!(
    out,
    i,
    j,
    grid,
    fields,
    hydrology::SoilHydrology{NF, VerticalFlow, SaturationClosure, SoilHydraulics} where {VerticalFlow<:Terrarium.AbstractVerticalFlow, SaturationClosure<:Terrarium.AbstractSoilWaterClosure, SoilHydraulics<:(Terrarium.AbstractSoilHydraulics{NF})},
    energy,
    strat,
    bgc
) where NF

Kernel function that adjusts saturation profiles to account for oversaturation and undersaturation arising due to numerical error. This implementation scans over the saturation profiles at each lateral grid cell and redistributes excess water upward layer-by-layer until reaching the topmost layer, and saturation deficits below the residual saturation downward and then back upward, drawing on the water of the neighboring layers. Water is conserved unless the column as a whole holds less water than its residual saturation, in which case the topmost layer is filled up to the residual saturation. The remaining surface excess is returned as the excess saturation of the topmost layer times its thickness (m); the caller multiplies by the topmost porosity to get a water depth (m³/m²) to route into a surface water pool, or discards it. With an energy process, every saturation change also changes the internal energy in out by the sensible heat of the water (compute_redistributed_heat!).

source
Terrarium.compute_redistributed_heat! — Function
compute_redistributed_heat!(out, i, j, k, grid, fields, energy, strat, bgc, k_source, Δs)

Change out.internal_energy of cell i, j, k by the sensible heat C_w T φ Δs of a saturation change Δs by redistributed water, with C_w the water heat capacity of the energy process, φ the porosity of cell k, and T the temperature of the cell k_source that the water leaves. A no-op without a SoilThermodynamics energy process.

source
Terrarium.compute_matric_potential — Function
compute_matric_potential(i, j, k, grid, fields)

Matric potential ψₘ = Ψ - ψ_h - ψ_z (m) at cell i, j, k, reconstructed from the total pressure_head Ψ by removing the hydrostatic head ψ_h = max(0, z₀ - z) below the water table z₀ and the elevation head ψ_z = z - z_ref relative to the soil surface z_ref, as in SoilSaturationPressureClosure.

source
Terrarium.specific_moisture_capacity — Function
specific_moisture_capacity(swrc::SWRC, ψ, θsat)

Specific moisture capacity ∂θ/∂ψ (1/m) of the soil water retention curve swrc at the matric potential ψ (m) for the saturated water content θsat, from the derivative provided by FreezeCurves.jl, in the number type of ψ; zero for saturated conditions.

julia> swrc = VanGenuchten(α = 2.0, n = 2.0);

julia> Terrarium.specific_moisture_capacity(swrc, 0.0, 0.4)
0.0

julia> Terrarium.specific_moisture_capacity(swrc, -0.5, 0.4) ≈ swrc(Terrarium.FreezeCurves.derivative, -0.5; θsat = 0.4)
true
source
Terrarium.compute_water_outflow_timescale — Function
compute_water_outflow_timescale(
    i,
    j,
    k,
    grid,
    clock,
    fields,
    hydrology::SoilHydrology{NF, RichardsEq, SaturationClosure, SoilHydraulics} where {SaturationClosure<:Terrarium.AbstractSoilWaterClosure, SoilHydraulics<:(Terrarium.AbstractSoilHydraulics{NF})},
    strat::Terrarium.AbstractStratigraphy,
    bgc::Terrarium.AbstractSoilBiogeochemistry,
    constants::PhysicalConstants,
    evtr::Union{Nothing, Terrarium.AbstractEvapotranspiration},
    timestepping::Timestepping,
    saturation_bcs::NamedTuple
) where NF

Kernel function returning the explicit water outflow timescale $τ = φ (s - s_res) Δz / q_out$ at cell i, j, k, with porosity φ, saturation s, the residual_saturation s_res, and the outflow $q_out = -Δz G_θ$ (m/s) where the explicit volumetric water content tendency $G_θ$ (compute_volumetric_water_content_tendency plus the flux boundary conditions saturation_bcs on saturation_water_ice, see compute_saturation_boundary_flux) is negative. Cells without net explicit outflow return Inf.

source
Terrarium.compute_saturation_boundary_flux — Function
compute_saturation_boundary_flux(
    bc::Oceananigans.BoundaryConditions.BoundaryCondition{<:Flux},
    i,
    j,
    grid,
    clock,
    fields
)

Kernel function returning the flux (positive upward) of the flux boundary condition bc on saturation_water_ice at column i, j, and zero for other boundary conditions.

source

The following kernel functions implement the implicit Richards step (see Implicit Richards equation):

Terrarium.compute_implicit_conductivity! — Method
compute_implicit_conductivity!(
    face_conductivity,
    i,
    j,
    k,
    grid,
    fields,
    hydrology::SoilHydrology{NF, RichardsEq, SaturationClosure, SoilHydraulics} where {SaturationClosure<:Terrarium.AbstractSoilWaterClosure, SoilHydraulics<:(Terrarium.AbstractSoilHydraulics{NF})}
) where NF

Store the hydraulic conductivity (m/s) of the current iterate at face i, j, k in face_conductivity, as in the explicit darcy_flux (darcy_conductivity).

source
Terrarium.compute_implicit_capacity! — Method
compute_implicit_capacity!(
    diagonal_coefficient,
    i,
    j,
    k,
    grid,
    fields,
    hydrology::SoilHydrology{NF, RichardsEq, SaturationClosure, SoilHydraulics} where {SaturationClosure<:Terrarium.AbstractSoilWaterClosure, SoilHydraulics<:(Terrarium.AbstractSoilHydraulics{NF})},
    strat::Terrarium.AbstractStratigraphy,
    bgc::Terrarium.AbstractSoilBiogeochemistry
) where NF

Store the specific moisture capacity C = ∂θ/∂ψ (1/m) of the current iterate at cell i, j, k in diagonal_coefficient (compute_specific_moisture_capacity), bounded below by minimum_moisture_capacity.

source
Terrarium.compute_implicit_water_flow_coefficients! — Function
compute_implicit_water_flow_coefficients!(
    out,
    i,
    j,
    k,
    grid,
    fields,
    hydrology::SoilHydrology{NF, RichardsEq, SaturationClosure, SoilHydraulics} where {SaturationClosure<:Terrarium.AbstractSoilWaterClosure, SoilHydraulics<:(Terrarium.AbstractSoilHydraulics{NF})},
    strat::Terrarium.AbstractStratigraphy,
    bgc::Terrarium.AbstractSoilBiogeochemistry
) where NF

Store the capacity (compute_implicit_capacity!) and face conductivity (compute_implicit_conductivity!) at i, j, k in out.

source
Terrarium.compute_implicit_water_flow_rhs — Function
compute_implicit_water_flow_rhs(
    i,
    j,
    k,
    grid,
    fields,
    diagonal_coefficient,
    explicit_saturation,
    strat::Terrarium.AbstractStratigraphy,
    bgc::Terrarium.AbstractSoilBiogeochemistry
)

Right-hand side C Ψᵐ - φ (ξᵐ - ξ⁺) at i, j, k, with C = diagonal_coefficient and ξ⁺ = explicit_saturation.

source
Terrarium.compute_implicit_saturation_update — Function
compute_implicit_saturation_update(
    i,
    j,
    k,
    grid,
    fields,
    face_conductivity,
    explicit_saturation,
    strat::Terrarium.AbstractStratigraphy,
    bgc::Terrarium.AbstractSoilBiogeochemistry,
    Δt
)

Saturation ξ = ξ⁺ + Δt Q_diff / φ at i, j, k, with Q_diff = ∂z(K ∂z Ψ) (1/s) through the interior faces (compute_implicit_diffusion_tendency), ξ⁺ = explicit_saturation, and K = face_conductivity.

source
Terrarium.compute_implicit_water_flux — Function
compute_implicit_water_flux(
    i,
    j,
    k,
    grid,
    water_flux,
    face_conductivity,
    pressure_head
)

Darcy flux q = -K ∂z Ψ (m/s, positive upward) at interior face i, j, k, with K = face_conductivity; at the bottom face, the current water_flux (from tⁿ) is kept.

source

References

[2]
G. Bonan. Soil Moisture. In: Climate Change and Terrestrial Ecosystem Modeling (Cambridge University Press, 2019); pp. 115–133.
[3]
M. A. Celia, E. T. Bouloutas and R. L. Zarba. A General Mass-Conservative Numerical Solution for the Unsaturated Flow Equation. Water Resources Research 26, 1483–1496 (1990).
[6]
J. Noilhan and J.-F. Mahfouf. The ISBA land surface parameterisation scheme. Global and Planetary Change 13, 145–159 (1996). Soil Moisture Simulation.
[7]
L. A. Richards. Capillary Conduction of Liquids through Porous Mediums. Physics 1, 318–333 (1931).
[8]
M. T. van Genuchten. A Closed-form Equation for Predicting the Hydraulic Conductivity of Unsaturated Soils. Soil Science Society of America Journal 44, 892–898 (1980).
[9]
S. Westermann, T. Ingeman-Nielsen, J. Scheer, K. Aalstad, J. Aga, N. Chaudhary, B. Etzelmüller, S. Filhol, A. Kääb, C. Renette, L. S. Schmidt, T. V. Schuler, R. B. Zweigel, L. Martin, S. Morard, M. Ben-Asher, M. Angelopoulos, J. Boike, B. Groenke, F. Miesner, J. Nitzbon, P. Overduin, S. M. Stuenzi and M. Langer. The CryoGrid Community Model (Version 1.0) – a Multi-Physics Toolbox for Climate-Driven Simulations in the Terrestrial Cryosphere. Geoscientific Model Development 16, 2607–2647 (2023).
[10]
R. H. Brooks and A. T. Corey. Hydraulic Properties of Porous Media. Hydrology Papers 3 (Colorado State University, Fort Collins, CO, 1964).
[11]
G. S. Campbell. A Simple Method for Determining Unsaturated Conductivity from Moisture Retention Data. Soil Science 117, 311–314 (1974).
[12]
R. B. Clapp and G. M. Hornberger. Empirical Equations for Some Soil Hydraulic Properties. Water Resources Research 14, 601–604 (1978).
[13]
A. Boone, V. Masson, T. Meyers and J. Noilhan. The Influence of the Inclusion of Soil Freezing on Simulations by a Soil–Vegetation–Atmosphere Transfer Scheme. Journal of Applied Meteorology 39, 1544–1569 (2000).