Soil models

Overview

SoilModel is a model of soil physics that couples the relevant processes governing energy, water, and carbon in natural soils. It is the primary model type for simulating heat conduction, freeze-thaw processes (e.g. permafrost), and variably saturated water flow in the subsurface.

arch = CPU()
grid = ColumnGrid(arch, Float32, ExponentialSpacing(N = 10)) # 10 soil layers
model = SoilModel(grid)
integrator = initialize(model)

Integrator of SoilModel{Float32, Oceananigans.Grids.RectilinearGrid{Float32, Oceananigans.Grids.Periodic, Oceananigans.Grids.Flat, Oceananigans.Grids.Bounded, Oceananigans.Grids.StaticVerticalDiscretization{Oceananigans.Utils.StaticStridesOffsetArray{Float32, 1, (17,), (-3,), Vector{Float32}}, Oceananigans.Utils.StaticStridesOffsetArray{Float32, 1, (16,), (-3,), Vector{Float32}}, Oceananigans.Utils.StaticStridesOffsetArray{Float32, 1, (17,), (-4,), Vector{Float32}}, Oceananigans.Utils.StaticStridesOffsetArray{Float32, 1, (16,), (-3,), Vector{Float32}}}, Float32, Float32, Oceananigans.Utils.StaticStridesOffsetArray{Float32, 1, (3,), (-1,), StepRangeLen{Float32, Float64, Float64, Int64}}, Oceananigans.Utils.StaticStridesOffsetArray{Float32, 1, (3,), (-1,), StepRangeLen{Float32, Float64, Float64, Int64}}, Float32, Float32, Nothing, Nothing, CPU, Oceananigans.Grids.GridSize{1, 1, 10, 1, 0, 3}}, SoilEnergyWaterCarbon{Float32, SoilStratigraphy{Float32, 1, Tuple{ConstantSoilHorizon{Float32, :soil, ConstantSoilPorosity{Float32}}}}, SoilThermodynamics{Float32, SoilEnergyTemperatureClosure, SoilThermalProperties{Float32, FreeWater, InverseQuadratic, SoilThermalConductivities{Float32}}}, SoilHydrology{Float32, RichardsEq, SoilSaturationPressureClosure, SoilHydraulicsSURFEX{Float32, BrooksCorey{FreezeCurves.SoilWaterVolume{Float32, Float32, Float32}, Float32, Float32}, UnsatKBrooksCorey{Float32}}, Nothing}, ConstantSoilCarbonDensity{Float32}}, DefaultInitializer{Float32}, IMEX{Float32, ForwardEuler{Float32}, BackwardEuler{Float32, Picard}}} with timestepper IMEX{Float32, ForwardEuler{Float32}, BackwardEuler{Float32, Picard}}(ForwardEuler{Float32}(3600.0f0), BackwardEuler{Float32, Picard}(Picard(2)))
├── Current time: 0.0
├── StateVariables{Float32}(clock = Clock{Float32, Float64}(time=0 seconds, iteration=0, last_Δt=Inf days), prognostic = (:internal_energy, :saturation_water_ice), auxiliary = (:temperature, :liquid_water_fraction, :pressure_head, :ground_temperature, :hydraulic_conductivity, :water_table, :water_flux), inputs = (), namespaces = (:soil,), timestepper_cache = IMEXCache)

mutable struct SoilModel{NF, GridType<:Oceananigans.Grids.AbstractGrid, Soil<:Terrarium.AbstractSoil, Initializer<:Terrarium.AbstractInitializer, Timestepper<:Terrarium.AbstractTimeStepper} <: Terrarium.AbstractSoilModel{NF, GridType}

The default timestepper of SoilModel depends on the soil hydrology. With Richards-equation water flow (the default), it is IMEX(ForwardEuler(NF; Δt = 3600), BackwardEuler(NF)), which treats soil heat conduction and water flow implicitly with a 1-hour time step. Without water flow (NoFlow), it is ForwardEuler(NF). Value boundary conditions such as PrescribedSurfaceTemperature require an explicit timestepper, e.g. SoilModel(grid; soil, timestepper = ForwardEuler(NF)). See Time stepping for details.

variables(model)

Variables
├─ Prognostic: 
├── internal_energy [J m^-3] on XYZ{Center, Center, Center} of the ground domain
├── saturation_water_ice [-] on XYZ{Center, Center, Center} of the ground domain
├─ Auxiliary: 
├── temperature [°C] on XYZ{Center, Center, Center} of the ground domain
├── liquid_water_fraction [-] on XYZ{Center, Center, Center} of the ground domain
├── pressure_head [m] on XYZ{Center, Center, Center} of the ground domain
├── ground_temperature [°C] on XY{Center, Center, Terrarium.Coordinate{typeof(lastindex), 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: 
├─ Namespaces:
├── Terrarium.Namespace{:soil, Terrarium.Variables}

Components

Each component is summarized briefly below. Follow the linked process pages for full theoretical background, available concrete types, state variables, and method signatures.

Stratigraphy

The strat component parameterizes the vertical distribution of soil material properties (texture, porosity, organic content). It provides kernel functions used by the energy and hydrology sub-processes to look up spatially varying material properties at each grid cell. By default HomogeneousSoilStratigraphy is used, which specifies a single uniform material throughout all vertical profiles. See Soil stratigraphy for details.

Energy balance

The energy component represents heat conduction in the soil column, including the latent heat of freeze-thaw phase change. The default implementation is SoilThermodynamics, which evolves the volumetric internal energy $U$ (J m⁻³) as the prognostic variable and derives temperature via the SoilEnergyTemperatureClosure. See Soil energy balance for details.

Hydrology

The hydrology component governs the vertical movement of water in the soil column. The default implementation is SoilHydrology, which evolves total saturation (liquid + ice) as the prognostic variable. By default, the vertical flow follows the Richards equation (RichardsEq); selecting NoFlow as the vertical_flow operator keeps the soil water immobile. See Soil hydrology for details.

Biogeochemistry

The biogeochem component simulates the spatial distribution of soil organic carbon and associated biogeochemical fluxes. The default implementation is ConstantSoilCarbonDensity, which prescribes a spatially homogeneous organic carbon density and does not define any prognostic variables.

Initializers

Terrarium provides a hierarchy of initializers for SoilModel. The top-level initializer is SoilInitializer, which composes separate sub-initializers for energy, hydrology, and biogeochemistry. It is the recommended initializer for most use cases:

struct SoilInitializer{NF, EnergyInit<:Terrarium.AbstractInitializer{NF}, HydrologyInit<:Terrarium.AbstractInitializer{NF}, BGCInit<:Terrarium.AbstractInitializer{NF}} <: Terrarium.AbstractInitializer{NF}

Passing a SoilInitializer to SoilModel overrides the DefaultInitializer:

initializer = SoilInitializer(Float32;
    energy    = QuasiThermalSteadyState(Float32; T₀ = 2.0f0),
    hydrology = SaturationWaterTable(Float32; water_table_depth = 3.0f0),
)
model = SoilModel(grid; initializer)

SoilModel{Float32} on CPU()
├── grid:  1×1×10 RectilinearGrid{Float32, Periodic, Flat, Bounded} on CPU with 1×0×3 halo
├── soil:  SoilEnergyWaterCarbon{Float32, SoilStratigraphy{Float32, 1, Tuple{ConstantSoilHorizon{Float32, :soil, ConstantSoilPorosity{Float32}}}}, SoilThermodynamics{Float32, SoilEnergyTemperatureClosure, SoilThermalProperties{Float32, FreeWater, InverseQuadratic, SoilThermalConductivities{Float32}}}, SoilHydrology{Float32, RichardsEq, SoilSaturationPressureClosure, SoilHydraulicsSURFEX{Float32, BrooksCorey{FreezeCurves.SoilWaterVolume{Float32, Float32, Float32}, Float32, Float32}, UnsatKBrooksCorey{Float32}}, Nothing}, ConstantSoilCarbonDensity{Float32}}
├── constants:  PhysicalConstants{Float32}
├── initializer:  SoilInitializer{Float32, QuasiThermalSteadyState{Float32, Float32, Float32}, SaturationWaterTable{Float32, Float32, Float32}, DefaultInitializer{Float32}}
├── timestepper:  IMEX{Float32, ForwardEuler{Float32}, BackwardEuler{Float32, Picard}}

Spatially varying initial values

The parameters of the energy and hydrology initializers are declared as input variables of the model, with the values given to the initializer as their defaults:

Each value may therefore be a number, a function of the horizontal node coordinates, an array, a Field, or an AbstractFieldInitializer. These defaults are re-applied at every initialization, so changes to their parameters take effect, and an InputSource with the same name and units always takes precedence over them. Spatial data, such as a regridded climatology, should preferably be supplied this way rather than as a Field stored in the initializer: the initializer is part of the model, which should stay independent of any particular grid and free of state. Field initializers such as LatitudinalClimatology expose their own parameters as model parameters. Note that geothermal_heat_flux is the same variable read by the GeothermalHeatFlux bottom boundary condition, so the initial profile and the boundary condition stay consistent.

column_grid = ColumnGrid(arch, Float32, ExponentialSpacing(N = 10), 3) # three columns
energy = QuasiThermalSteadyState(Float32; T₀ = x -> 2.0f0 * x, Qgeo = 0.05f0)
model = SoilModel(column_grid; initializer = SoilInitializer(Float32; energy))
integrator = initialize(model)
interior(integrator.state.initial_surface_temperature)

3×1×1 view(::Array{Float32, 3}, 4:6, 1:1, 1:1) with eltype Float32:
[:, :, 1] =
 0.33333334
 1.0
 1.6666666

Energy initializers

struct ConstantSoilTemperature{NF, T0} <: Terrarium.AbstractInitializer{NF}

struct QuasiThermalSteadyState{NF, T0, QG} <: Terrarium.AbstractInitializer{NF}

T(z) = T₀ - Qgeo / k_eff * z

struct PiecewiseLinearInitialSoilTemperature{NF, N}

initializer = PiecewiseLinearInitialSoilTemperature(
    0.0u"m" => 5.0, # always in °C!
    0.5u"m" => 2.0,
    1.0u"m" => 1.0,
    10.0u"m" => 1.5,
    ...
)

Hydrology initializers

struct SaturationWaterTable{NF, S, D} <: Terrarium.AbstractInitializer{NF}

julia> init = SaturationWaterTable(Float32; vadose_zone_saturation = VolumetricWaterContent(0.2f0));

julia> init.vadose_zone_saturation
VolumetricWaterContent{Float32}(0.2f0)

The vadose zone water content of SaturationWaterTable can be given either as saturation or as volumetric water content by wrapping the value in one of the following marker types; plain values are interpreted as saturation. Volumetric water content is converted to saturation by dividing by the soil porosity at each grid cell, so that it can be initialized directly from, e.g., soil moisture observations.

hydrology = SaturationWaterTable(Float32; vadose_zone_saturation = VolumetricWaterContent(0.2f0), water_table_depth = 3.0f0)
model = SoilModel(grid; initializer = SoilInitializer(Float32; hydrology))
integrator = initialize(model)
extrema(interior(integrator.state.saturation_water_ice))

(0.40816325f0, 1.0f0)

struct Saturation{T} <: Terrarium.AbstractWaterContent

struct VolumetricWaterContent{T} <: Terrarium.AbstractWaterContent

struct ConstantSaturation{NF, S} <: Terrarium.AbstractInitializer{NF}

Kernel functions

compute_quasi_steady_state_temperature(
    i,
    j,
    k,
    grid,
    fields,
    init::QuasiThermalSteadyState
)

compute_constant_temperature(
    i,
    j,
    k,
    grid,
    fields,
    ::ConstantSoilTemperature
)

compute_water_table_saturation(
    i,
    j,
    k,
    grid,
    fields,
    init::SaturationWaterTable,
    args...
)

compute_vadose_zone_saturation(
    i,
    j,
    k,
    grid,
    fields,
    saturation
)

compute_constant_saturation(
    i,
    j,
    k,
    grid,
    fields,
    ::ConstantSaturation
)

Fallback

Boundary conditions

Terrarium provides a set of named boundary condition aliases for the most common SoilModel configurations. These return NamedTuples in the format expected by the boundary_conditions keyword argument of initialize.

Energy boundary conditions

Hydrology boundary conditions