GitHub

Soil energy balance

Overview

Heat transfer in the subsurface can be represented according to the heat equation, with the upper boundary set to surface temperature and the lower boundary set to a constant positive heat flux representing heat produced by the inner earth [14, 15]. If both the upper and lower boundaries are assumed to be constant over time, the steady-state temperature profile takes the form of a continuous piecewise linear function increasing over depth with the slope determined by the thermal properties of the ground material. The instantaneous temperature field can then be generally represented as

\[\begin{equation} T(z,t) = T_0 + \frac{Q_{\text{geo}}}{\kappa_{\text{h}}(z)}z + \Delta T(z,t) \end{equation}\]

where $T(z,t)$ is the temperature field (K) over depth $z$ (m) and time $t$ (s), $T_0$ is the mean annual GST (K), $Q_{\text{geo}}$ is the geothermal heat flux (W/m²), and $\kappa_{\text{h}}(z)$ is the thermal conductivity which may vary with depth depending on the material (W/m K). The last term $\Delta T(z,t)$ represents transient disturbances to the steady state temperature profile due to both seasonal and long-term fluctuations in the upper and lower boundary conditions of the vertical domain. Simulating the impacts of these transient changes is one of the primary objectives of most numerical permafrost and land surface models.

Diffusive heat flow in a solid medium is governed by Fourier's law,

\[\begin{equation} \mathbf{j}_\text{h}^{\mathrm{v}} = \mathbf{j}_\text{h} \cdot \mathbf{n}_z = -\kappa_{\text{h}}\frac{\partial T}{\partial z}\,, \end{equation}\]

where $\mathbf{j}_\text{h}$ is the diffusive heat flux vector (W/m²) and $\mathbf{n}_z$ is the upward facing normal vector along the vertical $z$ axis.

Phase change of pore water/ice

Since ground materials are often porous, i.e., there exists void space between the solid particles, it is necessary to consider the potential presence of water and/or ice in this void space, which is hereafter referred to as pore space, or simply, soil pores. The thermal effects of water and ice can be accounted for by considering not only the temperature of the material but rather the total internal energy of the elementary volume. Combining the diffusive flux with a potential advective heat flux $\mathbf{j}_{\text{h}}^{\text{w}}$ due to water flow yields the energy conservation law,

\[\begin{equation} \label{eq:energyconservation} \frac{\partial U(T,\theta)}{\partial t} = - \boldsymbol{\nabla} \cdot \left(\mathbf{j}_\text{h} + \mathbf{j}_h^{\text{w}}\right) + F_h(z,t), \end{equation}\]

where $U(T,\theta)$ (J/m³) is the volumetric internal energy as a function of temperature and total water/ice content $\theta$ (m³/m³), and $F_h(z,t)$ is an inhomogeneous heat source/sink (forcing) term.

The advective heat flux $\mathbf{j}_{\text{h}}^{\text{w}}$ can be represented as,

\[\begin{equation} \mathbf{j}_{\text{h}}^{\text{w}} = \left( c_{\text{w}} T + L_{\text{sl}} \right) \mathbf{j}_{\text{w}} \rho_{\text{w}} \end{equation}\]

where $L_{\text{sl}}$ (J/kg) and $c_{\text{w}}$ (J/(kg K)) represent the specific latent heat of fusion (solid → liquid) and heat capacity of liquid water respectively. Note that $\mathbf{j}_{\mathrm{w}}$ is the water flux vector as defined in Soil hydrology. This flux term accounts for the energy transferred by the movement of water within the soil matrix. In model configurations that neglect subsurface water flow, this flux term is implicitly assumed to be zero.

Terrarium.SoilThermodynamics — Type
struct SoilThermodynamics{NF, EnergyClosure<:Terrarium.AbstractEnergyClosure, ThermalProps<:SoilThermalProperties{NF}} <: Terrarium.AbstractSoilThermodynamics{NF}

Soil thermal dynamics with freezing and thawing of pore water. The internal energy U evolves by vertical heat conduction and advection of sensible heat by the vertical water flux j_w (m/s, positive upward):

∂U/∂t = ∂z(κ ∂z T) - ∂z(C_w T j_w),

with temperature T (°C), thermal conductivity κ (W/m/K), and volumetric heat capacity of liquid water C_w (J/m³/K). The closure relates temperature and internal energy, $U(T,\phi)$, with additional parameters $\phi$ set by the model configuration.

Properties:

  • closure::Terrarium.AbstractEnergyClosure: Closure relating energy and temperature

  • thermal_properties::SoilThermalProperties: Soil thermal properties

source

As for Soil hydrology, only the vertical fluxes are currently considered. This simplifies the notation of $\eqref{eq:energyconservation}$ to

\[\begin{equation} \frac{\partial U(T,\theta)}{\partial t} = \frac{\partial}{\partial z}\left[\kappa_{\text{h}} \frac{\partial T}{\partial z}\right] + F_h(z,t)\,. \end{equation}\]

Process interface

Terrarium.initialize! — Method
initialize!(
    state,
    grid,
    energy::SoilThermodynamics,
    soil::Terrarium.AbstractSoil,
    constants::PhysicalConstants,
    args...
)
source
Terrarium.compute_tendencies! — Method
compute_tendencies!(
    state,
    grid,
    energy::SoilThermodynamics,
    soil::Terrarium.AbstractSoil
)
compute_tendencies!(
    state,
    grid,
    energy::SoilThermodynamics,
    soil::Terrarium.AbstractSoil,
    constants::Union{Nothing, PhysicalConstants}
)
compute_tendencies!(
    state,
    grid,
    energy::SoilThermodynamics,
    soil::Terrarium.AbstractSoil,
    constants::Union{Nothing, PhysicalConstants},
    timestepping::Timestepping
)

Compute the explicitly integrated part of the internal energy tendency of energy. The Timestepping argument is the class of internal_energy (see timestepping(state, ::Val)), i.e. how the variable is integrated, not which part of the tendency is computed. Explicit: the whole tendency is integrated explicitly, i.e. the full divergence of the diffusive and advective heat fluxes. Implicit: the fluxes through the interior faces are integrated by the implicit sub-stepper (BackwardEuler), so only the diffusive fluxes and the advective inflow through the boundary faces remain. Flux boundary conditions are included in both cases.

source

Advection of heat by soil water flow

When the soil water flows (a RichardsEq hydrology), it carries its heat along, as described by the advective heat flux $\mathbf{j}_{\text{h}}^{\text{w}}$ above. Since $U$ is measured relative to liquid water at 0 °C, moving liquid water carries no latent heat, and ice is immobile, so only the sensible heat of the liquid water is advected: the vertical advective heat flux is $C_{\text{w}}\, T\, j_{\text{w}}$, where $j_{\text{w}} = \mathbf{j}_{\text{w}} \cdot \mathbf{n}_z$ is the vertical water flux (m/s, positive upward; see Soil hydrology) and $C_{\text{w}} = \rho_{\text{w}} c_{\text{w}}$ the volumetric heat capacity of liquid water (the water heat capacity of SoilThermalProperties, the same value with which water enters the bulk heat capacity). The energy equation then reads

\[\frac{\partial U}{\partial t} = \frac{\partial}{\partial z}\left[\kappa_{\text{h}} \frac{\partial T}{\partial z}\right] - \frac{\partial}{\partial z}\left(C_{\text{w}}\, T\, j_{\text{w}}\right) + F_h(z,t)\,.\]

In the discretization, the water flux is evaluated at the cell faces and the temperature of the advective flux at a face is that of the upwind cell (first-order upwind; see advective_heat_flux).

The flux $j_{\text{w}}$ is the water_flux auxiliary variable of the hydrology, i.e. the Darcy flux that moved the water in the step: the flux at $t^n$ under explicit time stepping of the saturation, and the flux of the conservative update of the implicit Richards step otherwise (see Soil hydrology). Heat thus always moves with the water that moves. Under explicit time stepping of internal_energy, the advective term is part of the explicit tendency and subject to an advective CFL condition, reported by cell_advection_timescale (see Adaptive time stepping). Under implicit time stepping, it is part of the implicit heat operator (see Implicit vertical heat conduction) and has no CFL limit.

Water that crosses the column boundaries through flux boundary conditions on saturation_water_ice, such as the infiltration flux, carries the temperature of the cell it enters or leaves, i.e. infiltrating water enters at the temperature of the top soil cell; this heat is added by compute_boundary_advected_heat! from exactly the same boundary flux as the water. Water that the hydrology closure redistributes between cells, or removes from an oversaturated surface, carries its sensible heat as well (see adjust_saturation_profile! and compute_redistributed_heat!); the heat of excess water routed to the surface excess pool leaves the system. Without water flow (NoFlow) the advective flux is zero.

Implicit vertical heat conduction

The Timestepping trait argument of compute_tendencies! is the resolved timestepping class of internal_energy (see Implicit-explicit time stepping). It is Explicit unless the model is run with an IMEX timestepper, in which case it is Implicit.

  • Under Explicit(), the tendency contains the full divergence of the diffusive and advective heat fluxes, compute_energy_tendency, plus the flux boundary conditions.
  • Under Implicit(), the diffusive and advective heat fluxes through the interior cell faces are omitted, since they are integrated by the implicit sub-stepper. Only the diffusive fluxes through the two boundary faces of the column, via compute_boundary_energy_tendency, and the heat of the water entering the column through them with the boundary Darcy flux, via compute_boundary_advective_energy_tendency, remain. For flux boundary conditions the diffusive boundary fluxes vanish because the halo cells are filled with a zero gradient, and the boundary fluxes themselves enter through the flux boundary conditions on internal_energy. For a gradient boundary condition on temperature, they carry the prescribed flux $-\kappa\, \partial T / \partial z$ through the boundary face. The heat of the water leaving the column through the boundary faces (e.g. by gravity drainage) is part of the implicit operator, at the temperature of the boundary cell at $t^{n+1}$, so that an explicit drainage flux that overdraws the water of the boundary cell cannot remove more heat than the cell holds.

For the energy equation with advection (see Advection of heat by soil water flow),

\[\frac{\partial U}{\partial t} = \frac{\partial}{\partial z}\left[\kappa \frac{\partial T}{\partial z}\right] - \frac{\partial}{\partial z}\left(C_\text{w}\, T\, j_\text{w}\right) + G_\text{exp}\,,\]

where the second term is the advection of sensible heat by the vertical soil water flux $j_\text{w}$ (zero without water flow), and $G_\text{exp}$ contains the flux boundary conditions, the heat fluxes through the boundary faces, and any other explicit sources, the implicit step is a backward Euler step in temperature, solved with the Picard solver (see Implicit Euler with a Picard solver). Writing $C_\text{app} = \partial U / \partial T$ for the apparent heat capacity and $T^m$, $U^m$ for the latest iterate of $T^{n+1}$, $U^{n+1}$, the energy increment is linearized as $U^{n+1} - U^m \approx C_\text{app}\,(T^{n+1} - T^m)$, giving

\[C_\text{app}\, T^{n+1} - \Delta t\, \frac{\partial}{\partial z}\left[\kappa \frac{\partial T^{n+1}}{\partial z}\right] + \Delta t\, \frac{\partial}{\partial z}\left(C_\text{w}\, T^{n+1} j_\text{w}\right) = C_\text{app}\, T^m - U^m + U^+\,, \qquad U^+ = U^n + \Delta t\, G_\text{exp}\,.\]

Discretized in space, with the temperature of the advective flux at each face taken from the upwind cell (first-order upwind; see Tridiagonal operator), this system is tridiagonal in each column and is solved with a batched tridiagonal solver. The coefficients $\kappa$ (at cell faces) and $C_\text{app}$ (at cell centers) are lagged, i.e. evaluated at the latest iterate rather than at $T^{n+1}$. The water flux $j_\text{w}$ is the flux with which the water was moved in the same step (see Coupling of water and energy); for a fixed $j_\text{w}$, the upwind terms keep the system an M-matrix (see Tridiagonal operator), so the solution has no spurious oscillations and the advection needs no CFL limit. The first iterate is the internal energy at the beginning of the step, $U^n$, stored by update_cache!, with the temperature re-diagnosed from it at the current saturation.

The FreeWater freeze curve has an apparent heat capacity that is a Dirac delta at the melting point, which cannot be used directly in a linear system. Terrarium therefore regularizes it: for cells undergoing phase change (liquid water fraction strictly between zero and one), the latent heat term is replaced by a large constant capacity,

\[C_\text{app} = C + C_\text{pc}\,, \qquad C_\text{pc} = 1.5 \times 10^8~\text{J/m}^3\text{/K}\,,\]

with bulk heat capacity $C$. Outside of phase change, $C_\text{app} = C$. For typical water contents, $C_\text{pc}$ corresponds to spreading the latent heat over an interval of roughly 1 K. Since the energy is updated conservatively (see below), $C_\text{pc}$ only affects how quickly the temperature of a freezing or thawing cell may change between iterates, not the energy balance.

Once the temperature has been solved for, the internal energy is updated conservatively from the implicit diffusive and advective fluxes through the interior faces,

\[U^{n+1} = U^+ + \Delta t\, \frac{\partial}{\partial z}\left[\kappa \frac{\partial T^{n+1}}{\partial z}\right] - \Delta t\, \frac{\partial}{\partial z}\left(C_\text{w}\, T^{n+1} j_\text{w}\right)\,,\]

rather than from the linearization. The temperature and liquid water fraction are then diagnosed from the closure, after every iteration including the last, so that they are consistent with $U^{n+1}$ when the implicit step returns. Since the interior fluxes cancel pairwise between adjacent cells, the change in the column-integrated energy over a step equals $\Delta t$ times the net boundary flux contained in $G_\text{exp}$, independent of the number of iterations and the regularization.

The lagged apparent heat capacity and face conductivity are given by compute_implicit_capacity! and compute_implicit_conductivity!, and, with water flow, the advective face velocity $C_\text{w} j_\text{w}$ by compute_implicit_advective_velocity!. The implicit step is supported for the FreeWater freeze curve, flux boundary conditions on internal_energy, and gradient boundary conditions on temperature; value boundary conditions on temperature are not yet supported.

The implicit operator of internal_energy is a TridiagonalOperator in temperature, with a face velocity for the advection if the water flows, returned by implicit_operator (which also validates the configuration), and each iteration of the implicit solver calls the process hooks of SoilThermodynamics below. The linear system is solved for the temperature, the variable of the SoilEnergyTemperatureClosure (see its varname method below). If saturation_water_ice is integrated implicitly as well, the Richards step is solved first, so that the heat step is evaluated at the saturation $\xi^{n+1}$ at the end of the step and advects heat with the water flux of the step (see Coupling of water and energy); linearize! therefore re-evaluates the energy closure at the current saturation before the coefficients.

Terrarium.implicit_operator — Method
implicit_operator(
    energy::SoilThermodynamics,
    hydrology::Terrarium.AbstractSoilHydrology,
    state,
    grid
)

Return the implicit operator of the soil internal energy: a TridiagonalOperator on ground_domain(grid) for heat conduction in temperature, with a face velocity for upwind advection of heat by the water flux if the hydrology lets the water flow (SoilHydrology with RichardsEq). Validates on the host that the freeze curve is FreeWater (the only one with an implemented regularized apparent heat capacity) and that temperature has no value boundary conditions (not yet supported). Flux boundary conditions on internal_energy and gradient boundary conditions on temperature are supported; the latter prescribe the diffusive flux -κ ∂T/∂z through the boundary face, which enters the explicit tendency via compute_boundary_energy_tendency.

The operator is the backward Euler step of ∂U/∂t = ∂z(κ ∂z T) - ∂z(C_w T j_w) + G_exp in temperature, with first-order upwind advection by the vertical water flux j_w (T_up the temperature of the upwind cell of a face) and a conservative update of the energy (see the soil energy docs for the derivation). With the Picard solver, for iterations m = 0, …, M-1, starting from U⁰ = Uⁿ (plus any heat moved by the closure of the water step) and T⁰ ← closure(U⁰):

  1. linearize!: Tᵐ ← closure(Uᵐ), lagged κᵐ (faces), regularized apparent heat capacity C_appᵐ (centers), advective face velocity C_w j_w (faces), and rhs = C_appᵐ Tᵐ - Uᵐ + U⁺ with U⁺ = U⁰ + Δt G_exp;
  2. solve!: C_appᵐ Tᵐ⁺¹ - Δt ∂z(κᵐ ∂z Tᵐ⁺¹) + Δt ∂z(C_w j_w T_upᵐ⁺¹) = rhs;
  3. implicit_update!: Uᵐ⁺¹ = U⁺ + Δt [∂z(κᵐ ∂z Tᵐ⁺¹) - ∂z(C_w j_w T_upᵐ⁺¹)] through the interior faces only;
  4. implicit_closure!: Tᵐ⁺¹ and the liquid water fraction ← closure(Uᵐ⁺¹).

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

source
Terrarium.linearize! — Method
linearize!(
    operator::TridiagonalOperator,
    rhs,
    state,
    grid,
    energy::SoilThermodynamics,
    soil::Terrarium.AbstractSoil,
    constants::PhysicalConstants,
    explicit_energy,
    Δt
)

Linearize the implicit heat operator of the soil internal energy at the current iterate (linearly implicit backward Euler in temperature, conservative update of the internal energy). On entry, explicit_energy holds U⁺ = U⁰ + Δt G_exp, where G_exp contains the flux boundary conditions (e.g. SoilHeatFlux and GeothermalHeatFlux), the diffusive fluxes implied by gradient boundary conditions on temperature, the advective heat fluxes through the boundary faces, and the heat carried by water crossing the boundaries (compute_boundary_advected_heat!), but not the fluxes through the interior faces. This hook

  1. re-evaluates the energy closure, so that temperature and liquid water fraction are consistent with the internal energy and the current saturation (changed by the water step since tⁿ),
  2. evaluates the face conductivity κ (compute_implicit_conductivity!), the regularized apparent heat capacity C_app (compute_implicit_capacity!), and, with water flow, the advective face velocity C_w j_w from the water_flux of the step (compute_implicit_advective_velocity!), storing them in the face_conductivity, diagonal_coefficient, and face_velocity of operator,
  3. and assembles the right-hand side C_app Tᵐ - Uᵐ + U⁺ of C_app Tᵐ⁺¹ - Δt ∂z(κ ∂z Tᵐ⁺¹) + Δt ∂z(C_w j_w T_upᵐ⁺¹) = C_app Tᵐ - Uᵐ + U⁺ (compute_implicit_energy_rhs).

implicit_update! then updates the internal energy conservatively: for any number of iterations, the change of column-integrated energy is exactly Δt times the net boundary flux in G_exp.

source
Terrarium.implicit_update! — Method
implicit_update!(
    state,
    grid,
    energy::SoilThermodynamics,
    soil::Terrarium.AbstractSoil,
    operator::TridiagonalOperator,
    explicit_energy,
    Δt
)

Conservative update of the soil internal energy, Uᵐ⁺¹ = U⁺ + Δt [∂z(κ ∂z Tᵐ⁺¹) - ∂z(C_w j_w T_upᵐ⁺¹)], from the implicit diffusive and advective fluxes through the interior faces (compute_implicit_energy_update), with U⁺ = explicit_energy and the lagged face conductivity κ and face velocity C_w j_w of operator.

source

After each iteration, the temperature and liquid water fraction are re-diagnosed with the regular energy closure, which has no once-per-step effects.

Coupling of water and energy

With a Richards-equation hydrology, the water flux advects sensible heat: the soil energy equation carries the advective term $-\partial_z (C_\text{w} T\, j_\text{w})$, discretized with first-order upwinding (see Advection of heat by soil water flow). For the scheme to keep a uniform temperature uniform, heat must move with exactly the water that moves. Terrarium ensures this as follows.

  • Water before energy. The SplitOperator of SoilModel and LandModel solves saturation_water_ice before internal_energy. The heat step therefore sees the saturation at the end of the step and the water flux of the step, while the Richards step sees the temperature and liquid water fraction at the beginning of the step: the effect of frozen ground on the hydraulic conductivity is lagged by one step, as in most sequentially split land models. The Stefan problem test (test/soil/soil_freeze_thaw_tests.jl) quantifies this against a fully explicit reference.
  • The water flux of the step. The water_flux auxiliary variable holds the Darcy flux through the cell faces that moved the water: the flux at $t^n$, computed with the tendencies, through the boundary faces and, under explicit time stepping of the saturation, through the interior faces; under implicit time stepping, the flux of the conservative update of the last Picard iteration through the interior faces. The advective heat fluxes use this flux, explicitly when internal_energy is explicit, and implicitly (upwind, in the heat operator) when it is implicit. In the implicit case, water that enters the column through a boundary face carries the temperature of the boundary cell at $t^n$ (explicitly), while water that leaves it, e.g. by gravity drainage, carries the temperature of the boundary cell at $t^{n+1}$ (implicitly), so that a large explicit drainage flux cannot remove more heat from the bottom cell than it holds. The heat carried by flux boundary conditions on the saturation (e.g. infiltration, see compute_boundary_advected_heat!) is explicit in both cases, since the water flux is.
  • Heat moved by the hydrology closure. adjust_saturation_profile! moves the sensible heat $C_\text{w} T$ of the water it shifts between cells, at the temperature of the cell that the water leaves, and removes the heat of the excess water that it discards or routes to the surface excess pool. The hydrology closure runs before the energy closure, so the temperature is recomputed from the adjusted energy. When both variables are implicit, the internal energy then holds $U^n$ plus this adjustment, and since the Picard solver forms the explicit stage of the heat step from the current internal energy and the explicit increment, the adjustment is kept, as in a sequential split.

The water heat capacity $C_\text{w}$ is the one with which water enters the bulk heat capacity. Since the bulk heat capacity also includes air, water that replaces air in the pores changes the temperature slightly (by a fraction of the order of $C_\text{air} / C_\text{w} \approx 3 \times 10^{-4}$ of the change of the water content), as the heat of the displaced air is not tracked.

Implicit water with explicit energy

Integrating saturation_water_ice implicitly while internal_energy is explicit (by specializing timestepping) is allowed, but not recommended with water flow: the explicit advection of heat is then computed at $t^n$, before the water is moved by the implicit flux of the step, so heat and water move with different fluxes. During wetting fronts, this can warm or cool the soil substantially and become unstable for large time steps. The default routing of SoilModel and LandModel integrates both implicitly.

Under explicit time stepping of internal_energy, the stability limit of the explicit advection is reported by cell_advection_timescale (see Adaptive time stepping); under implicit time stepping, it is Inf.

Closures

The constitutive relationship between energy and temperature plays a critical role in simulating the thermal state of porous media. This relation can be defined in integral form as

\[\begin{equation} U(T,\theta) = \int_{T_{\text{ref}}}^T \tilde{C}(x,\theta) \, \mathrm{d}x \end{equation}\]

where $\tilde{C}$ is referred to as the effective or apparent heat capacity and $T_{\text{ref}}$ is a reference temperature. The internal energy $U$ is defined relative to a reference temperature $T_{\text{ref}}$ and a reference physical state if phase change is considered. For Terrarium, we choose $T_{\text{ref}} = 273.15$ K. We assume water to be the only constituent undergoing phase change; as such, the total internal energy of the soil volume is computed as,

\[\begin{equation} U(T,\theta) = C(\theta_{\text{w}},\theta) (T - T_{\text{ref}}) - \rho_{\text{w}} L_{\text{sl}} (\theta - \theta_{\text{w}}(\theta, T)), \end{equation}\]

where $\theta_{\text{w}}(T,\theta)$ is the volumetric unfrozen water content as a function of temperature $T$ and total water/ice content $\theta$. Note that, under this formulation, zero internal energy corresponds to zero ice content (fully thawed) at 0°C. $C(\theta_{\text{w}},\theta)$ (J/(K m³)) is the bulk volumetric material heat capacity of the total soil volume as a function of the unfrozen and total ($\theta$) water contents, calculated as $\sum_x \rho_x c_x$. The apparent heat capacity is then defined as the derivative of the energy-temperature relation,

\[\begin{equation} \tilde{C}(T,\theta) := \frac{\partial U}{\partial T} = \overbrace{C(\theta_{\text{w}},\theta) + T \frac{\partial C}{\partial \theta_{\text{w}}}\frac{\partial \theta_{\text{w}}}{\partial T}}^{\text{Sensible}} \, + \, \overbrace{\rho_{\text{w}} L_{\text{sl}} \frac{\partial\theta_{\text{w}}}{\partial T}}^{\text{Latent}}\,, \end{equation}\]

where the chain-rule is applied on $C(\theta_{\text{w}},\theta) (T - T_{\text{ref}})$ as both $C$ and $\theta_{\text{w}}$ are functions of $T$. The grouping of terms on the right-hand side show the partitioning of energy change into sensible and latent heat. The sensible component represents the energy necessary to heat a volume of the material to a particular temperature, whereas the latent component corresponds to the energy required for the phase change of water in the volume from solid (frozen) to liquid (thawed).

In the simplest case where we neglect the effect of capillary action in the soil, the energy-temperature relation can be derived according to that of "free" water (i.e. unbound by the soil matrix),

\[\begin{equation} \theta_{\text{w}}(U) = \begin{cases} 0 & U < -\rho_{\text{w}}L_{\text{sl}}\theta \\ \theta \left( 1 + \frac{U}{\rho_{\text{w}}L_{\text{sl}}\theta} \right) & -\rho_{\text{w}}L_{\text{sl}}\theta \leq U < 0 \\ \theta & U \geq 0\,, \end{cases} \end{equation}\]

with temperature then determined by

\[\begin{equation} T - T_{\text{ref}} = U^{-1}(U(T,\theta)) = \begin{cases} \frac{U + \rho_{\text{w}}L_{\text{sl}}\theta}{C} & U < -\rho_{\text{w}}L_{\text{sl}}\theta \\ 0 & -\rho_{\text{w}}L_{\text{sl}}\theta \leq U < 0 \\ \frac{U}{C} & U \geq 0\,. \end{cases} \end{equation}\]

By expressing $T$ in °C, $T_{ref}$ can be dropped as it equals 0°C. For more information on more complex formulations of these so-called soil freezing characteristic curves, see FreezeCurves.jl.

Terrarium.SoilEnergyTemperatureClosure — Type
struct SoilEnergyTemperatureClosure <: Terrarium.AbstractEnergyClosure

Defines the constitutive relationship between the the internal energy and temperature of a soil volume, i.e.

\[U(T) = T\times C(T) - \rho_w L_{sl} \theta (1 - F(T))\]

where T is temperature (°C), C(T) is the temperature-dependent heat capacity (J/m³/K), ρw L{sl} θ is the volumetric latent heat of fusion (J/m³), and F(T) = θ_w/θ is the constitutive relation between T and the unfrozen fraction of pore water with θ the sum of the volumetric fractions of water and ice. Note that, under this formulation, zero energy corresponds to 0°C with no ice, i.e. all pore water fully thawed.

The closure relation is defined as being a mapping from the conserved quantity (energy) to the continuous quantity (temperature), i.e. the inverse of U(T).

source
Terrarium.varname — Method
varname(::SoilEnergyTemperatureClosure) -> Symbol

Name of the variable returned by the SoilEnergyTemperatureClosure (varname interface for closures): temperature. The closure also diagnoses liquid_water_fraction, but the temperature is its primary variable, in which the implicit heat operator of internal_energy is solved (see BackwardEuler).

source

Methods

Terrarium.compute_boundary_advected_heat! — Function
compute_boundary_advected_heat!(
    state,
    grid,
    energy::SoilThermodynamics,
    hydrology::SoilHydrology{NF, RichardsEq, SaturationClosure, SoilHydraulics} where {SaturationClosure<:Terrarium.AbstractSoilWaterClosure, SoilHydraulics<:(Terrarium.AbstractSoilHydraulics{NF})},
    strat::Terrarium.AbstractStratigraphy,
    bgc::Terrarium.AbstractSoilBiogeochemistry
) where NF

Add the sensible heat of water crossing the column boundaries through flux boundary conditions on saturation_water_ice (e.g. InfiltrationFlux) to the internal energy tendency (compute_boundary_advected_energy_tendency). Must be called right after the hydrology boundary conditions, while the saturation tendency holds only their contributions. No-op without water flow.

source

Kernel functions

Terrarium.compute_boundary_energy_tendency — Function
compute_boundary_energy_tendency(
    i,
    j,
    k,
    grid,
    fields,
    energy::Terrarium.AbstractThermodynamics,
    args...
)

Internal energy tendency ∂U∂t at i, j, k from the diffusive_heat_flux q through the boundary faces of the column only (bottom of k = 1, top of k = Nz):

∂U∂t = -(δ_{k,Nz} q_{k+½} - δ_{k,1} q_{k-½}) / Δz_k.

Zero for flux and no-flux conditions, whose halos have zero gradient (flux conditions are added separately); nonzero, e.g., for a gradient condition on temperature.

source
Terrarium.advective_heat_flux — Function
advective_heat_flux(
    i,
    j,
    k,
    grid,
    fields,
    energy::Terrarium.AbstractThermodynamics,
    args...
)

Advective heat flux (W/m², positive upward) at the cell face i, j, k. The default is zero, i.e. no heat is transported by a moving medium; see the method for soil water flow following the Richards equation.

source
advective_heat_flux(
    i,
    j,
    k,
    grid,
    fields,
    energy::SoilThermodynamics,
    hydrology::SoilHydrology{NF, RichardsEq, SaturationClosure, SoilHydraulics} where {SaturationClosure<:Terrarium.AbstractSoilWaterClosure, SoilHydraulics<:(Terrarium.AbstractSoilHydraulics{NF})},
    args...
) where NF

Upwind advective heat flux C_w T_up j_w (W/m², positive upward) at face i, j, k, the first-order upwind discretization of C_w T j_w, with j_w (m/s) the water_flux of the step (compute_advective_water_flux) and C_w the water heat capacity of the SoilThermalProperties, consistent with the energy closure. T_up is the temperature of cell k for j_w < 0 and of cell k - 1 otherwise, clamped to the column at the boundary faces. Only the sensible heat of liquid water is advected; ice is immobile.

Under Implicit time stepping, only the boundary inflow is explicit (compute_boundary_advective_energy_tendency); the rest is part of the implicit heat operator (compute_implicit_advective_velocity!). The heat of water entering through flux boundary conditions on saturation_water_ice (e.g. InfiltrationFlux) is added by compute_boundary_advected_heat!.

source
Terrarium.compute_boundary_advective_energy_tendency — Function
compute_boundary_advective_energy_tendency(
    i,
    j,
    k,
    grid,
    fields,
    energy::Terrarium.AbstractThermodynamics,
    args...
)

Compute the part of the internal energy tendency ∂U∂t at i, j, k due to the advective heat flux through the boundary faces of the column that remains explicit when the interior fluxes are integrated implicitly (see BackwardEuler). The default is zero, i.e. no heat is transported by a moving medium; see the method for soil water flow following the Richards equation.

source
compute_boundary_advective_energy_tendency(
    i,
    j,
    k,
    grid,
    fields,
    energy::SoilThermodynamics,
    hydrology::SoilHydrology{NF, RichardsEq, SaturationClosure, SoilHydraulics} where {SaturationClosure<:Terrarium.AbstractSoilWaterClosure, SoilHydraulics<:(Terrarium.AbstractSoilHydraulics{NF})},
    args...
) where NF

Explicit advective energy tendency at i, j, k under Implicit time stepping: the heat of water entering the column through the boundary faces with the Darcy flux at tⁿ, at the temperature of the boundary cell. Outflow is implicit (compute_implicit_advection_tendency), so that it cannot remove more heat than a cell holds.

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

Internal energy tendency C_w T φ G_s (W/m³) at i, j, k of water entering (G_s > 0) or leaving (G_s < 0) the cell through flux boundary conditions on saturation_water_ice, with G_s = saturation_tendency[i, j, k], porosity φ, and cell temperature T. Infiltration thus enters at the temperature of the top cell.

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

Kernel function returning the advective timescale $τ = C Δz / (C_w |j_w|)$ at cell i, j, k, with C the bulk volumetric heat capacity of the cell, C_w that of liquid water, and |j_w| the larger magnitude of the darcy_flux through the two faces of the cell. This is the stability limit of the explicit upwind advection, which changes the cell temperature at the rate C_w j_w / C (m/s). Cells without flow return Inf.

source
Terrarium.compute_energy_tendencies! — Function
compute_energy_tendencies!(
    tendencies,
    i,
    j,
    k,
    grid,
    fields,
    energy::SoilThermodynamics,
    ::Explicit,
    args...
)

Add the explicitly integrated part of the internal energy tendency at i, j, k when internal_energy is integrated explicitly, i.e. the full tendency (compute_energy_tendency).

source
compute_energy_tendencies!(
    tendencies,
    i,
    j,
    k,
    grid,
    fields,
    energy::SoilThermodynamics,
    ::Implicit,
    args...
)

Add the explicitly integrated part of the internal energy tendency at i, j, k when internal_energy is integrated implicitly: the interior fluxes belong to the implicit operator, so only the diffusive fluxes through the boundary faces (compute_boundary_energy_tendency) and the heat of water entering through them (compute_boundary_advective_energy_tendency) remain.

source

The following kernel functions implement the implicit heat conduction and advection step (see Implicit vertical heat conduction):

Terrarium.compute_implicit_conductivity! — Function
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
compute_implicit_conductivity!(
    face_conductivity,
    i,
    j,
    k,
    grid,
    fields,
    energy::SoilThermodynamics,
    hydrology::Terrarium.AbstractSoilHydrology,
    strat::Terrarium.AbstractStratigraphy,
    bgc::Terrarium.AbstractSoilBiogeochemistry
)

Store the lagged thermal conductivity at face i, j, k in face_conductivity, interpolated as in the explicit diffusive_heat_flux.

source
Terrarium.compute_implicit_capacity! — Function
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
compute_implicit_capacity!(
    diagonal_coefficient,
    i,
    j,
    k,
    grid,
    fields,
    energy::SoilThermodynamics,
    ::FreeWater,
    hydrology::Terrarium.AbstractSoilHydrology,
    strat::Terrarium.AbstractStratigraphy,
    bgc::Terrarium.AbstractSoilBiogeochemistry,
    constants::PhysicalConstants
)

Store the regularized apparent heat capacity C_app (J/m³/K) at cell i, j, k in diagonal_coefficient, for the FreeWater freeze curve,

C_app = C + C_pc    if 0 < F < 1 (phase change),
C_app = C           otherwise,

where C is the bulk heat capacity, F the liquid water fraction of the current iterate, and C_pc = 1.5 × 10⁸ J/m³/K a fixed regularization of the latent heat term (a Dirac delta at 0 °C for FreeWater). C_pc only damps the temperature change of phase-change cells between iterates and does not enter the conservative energy update. It is not applied to fully frozen or thawed cells, where it would slow the lagged temperature of cells that keep cooling (warming) and overestimate their conductive fluxes.

source
Terrarium.compute_implicit_advective_velocity! — Function
compute_implicit_advective_velocity!(face_velocity, i, j, k, grid, fields, energy, hydrology)

Store the advective face velocity C_w j_w (W/m²/K) of the implicit heat operator at face i, j, k in face_velocity, and at the top boundary face for k = Nz, with j_w the water flux (m/s) of the step through the face (compute_advective_water_flux) and C_w the volumetric heat capacity of liquid water, as in the explicit advective_heat_flux. Through the boundary faces, the operator uses only the outflow (see TridiagonalOperator). A no-op without advection (face_velocity = nothing).

source
Terrarium.compute_implicit_heat_coefficients! — Function
compute_implicit_heat_coefficients!(
    out,
    i,
    j,
    k,
    grid,
    fields,
    energy::SoilThermodynamics,
    fc::FreeWater,
    hydrology::Terrarium.AbstractSoilHydrology,
    strat::Terrarium.AbstractStratigraphy,
    bgc::Terrarium.AbstractSoilBiogeochemistry,
    constants::PhysicalConstants
)

Store the lagged apparent heat capacity (compute_implicit_capacity!), face conductivity (compute_implicit_conductivity!), and, with water flow, advective face velocity (compute_implicit_advective_velocity!) at i, j, k in out.diagonal_coefficient, out.face_conductivity, and out.face_velocity.

source
Terrarium.compute_implicit_energy_rhs — Function
compute_implicit_energy_rhs(
    i,
    j,
    k,
    grid,
    fields,
    diagonal_coefficient,
    explicit_energy
)

Right-hand side C_app Tᵐ - Uᵐ + U⁺ of the linearized system at iterate m of the implicit solver, with C_app = diagonal_coefficient and U⁺ = explicit_energy.

source
Terrarium.compute_implicit_energy_update — Function
compute_implicit_energy_update(
    i,
    j,
    k,
    grid,
    face_conductivity,
    face_velocity,
    temperature,
    explicit_energy,
    Δt
)

Conservative energy update U = U⁺ + Δt (Q_diff + Q_adv) at i, j, k from the implicit heating rates (W/m³) of conduction, Q_diff = ∂z(κ ∂z T) through the interior faces (compute_implicit_diffusion_tendency), and of advection, Q_adv = -∂z(C_w j_w T_up) through the interior faces and the boundary outflow (compute_implicit_advection_tendency), with U⁺ = explicit_energy, κ = face_conductivity, and C_w j_w = face_velocity (nothing without water flow).

source
Terrarium.compute_implicit_diffusion_tendency — Function
compute_implicit_diffusion_tendency(
    i,
    j,
    k,
    grid,
    face_conductivity,
    ϕ
)

Diffusive tendency [∂z(κ ∂z ϕ)]_k, the rate of change of the conserved quantity from the diffusive fluxes through the interior faces of cell k only (consistent with the masked coefficients), used for the conservative update after each solve. Mirrors the explicit diffusive_heat_flux with the lagged face_conductivity.

source
Terrarium.compute_implicit_advection_tendency — Function
compute_implicit_advection_tendency(
    i,
    j,
    k,
    grid,
    face_velocity,
    ϕ
)

Advective tendency -[∂z(u ϕ_up)]_k, the rate of change of the conserved quantity from the upwind fluxes (compute_implicit_advective_flux) through the faces of cell k, used for the conservative update after each solve. Through the boundary faces of the column only the outflow u ϕ_k is included; the inflow enters through the right-hand side. Zero without advection (face_velocity = nothing).

source

References

[14]
J. C. Jaeger. Application of the Theory Of Heat Conduction to Geothermal Measurements. In: Terrestrial Heat Flow (American Geophysical Union (AGU), 1965); Chapter 2, pp. 7–23.
[15]
A. H. Lachenbruch and B. V. Marshall. Changing Climate: Geothermal Evidence from Permafrost in the Alaskan Arctic. Science 234, 689–696 (1986).