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...
)
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...
)
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
)
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}
)
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},
    runoff::Union{Nothing, Terrarium.AbstractSurfaceRunoff},
    args...
)
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}
├── water_table [m] on XY{Center, Center}
├── hydraulic_conductivity [m s^-1] on XYZ{Center, Center, Face}
├─ Inputs: 
├── liquid_water_fraction [-] on XYZ{Center, Center, Center}
├─ 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 [4]. 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$ [5].

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.
  • surface_excess_water: excess water at the soil surface (m^3/m^2).
  • 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).

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}
├── surface_excess_water [m] on XY{Center, Center}
├─ Auxiliary: 
├── pressure_head [m] on XYZ{Center, Center, Center}
├── hydraulic_conductivity [m s^-1] on XYZ{Center, Center, Face}
├── water_table [m] on XY{Center, Center}
├─ Inputs: 
├── liquid_water_fraction [-] on XYZ{Center, Center, Center}
├─ Namespaces:

Hydraulic properties

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 [6], here presented in the form of [7, 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 [6] extended with an ice impedance factor [7].

References

  • [6] Van Genuchten, Soil Science Society of America Journal (1980)
  • [7] Westermann et al., Geoscientific Model Development (2023)
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 [3] 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<:Terrarium.AbstractUnsatK{NF}}

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 sat_hydraulic_cond

  • sat_hydraulic_cond::Any: Hydraulic conductivity at saturation

  • field_capacity::Any: Constant field capacity

  • wilting_point::Any: Constant wilting point

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

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

Soil hydraulics parameterization that includes the SURFEX [3, 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 sat_hydraulic_cond

  • sat_hydraulic_cond::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

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

References

  • [3] 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
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})}
)

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!(
    out,
    i,
    j,
    grid,
    hydrology::SoilHydrology{NF, VerticalFlow, SaturationClosure, SoilHydraulics} where {VerticalFlow<:Terrarium.AbstractVerticalFlow, SaturationClosure<:Terrarium.AbstractSoilWaterClosure, SoilHydraulics<:(Terrarium.AbstractSoilHydraulics{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, where any remaining excess water is added to the surface_excess_water pool.

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

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

References

[3]
J. Noilhan and J.-F. Mahfouf. The ISBA land surface parameterisation scheme. Global and Planetary Change 13, 145–159 (1996). Soil Moisture Simulation.
[4]
L. A. Richards. Capillary Conduction of Liquids through Porous Mediums. Physics 1, 318–333 (1931).
[5]
G. Bonan. Soil Moisture. In: Climate Change and Terrestrial Ecosystem Modeling (Cambridge University Press, 2019); pp. 115–133.
[6]
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).
[7]
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).