Soil hydrology

Overview

The soil hydrology module consists of four prognostic states, comprised of three soil layers and a snow layer. The snow layer is described in reference snow section. The soil layers (\(W_{ly1}\), \(W_{ly2}\), and \(W_{ly3}\)) are characterized functionally distinct entities:

  • \(W_{ly1}\) represents near-surface moisture for bare soil evaporation and plant water uptake. The energy state of this layer (described in reference energy balance section) is used to determine the temperature of heterotrophic respiration (see section reference heterotrophic section).

  • \(W_{ly2}\) represents an additional reservoir accessible by plant roots for water uptake, but sits below \(W_{ly1}\) and therefore does not experience direct evaporation.

  • \(W_{ly3}\) represents water that is below \(W_{ly2}\) and is innaccessible to plant roots.

Functionally, water enters \(W_{ly1}\) through infiltration, water is exchanged between layers as a function of hydraulic gradients, and is removed via (i) bare soil evaporation from \(W_{ly1}\), (ii) plant transpiration from \(W_{ly1}\) and \(W_{ly2}\), and (iii) subsurface runoff each soil layer.

Prognostic and diagnostic water states

Each of the \(n\) soil water layers (\(W_{lyn}\)) consists of a prognistic water mass in units of kgH2O m-2, which is assumed to have a density of 1000 kg m-3 and therefore simplified as being in units of mmH2O. At each time step \(t\), a diagnostic soil moisture state (\(\theta_{lyn}\), in units of m3 m-3) can be derived as

(1)\[\theta_{lyn}(t) = W_{lyn}(t)(10^3p_{lyn}z_{lyn})^{-1}\]

where \(z_{lyn}\) is the depth parameter of layer \(n\) in m, and \(p_{lyn}\) is the porosity parameter of layer \(n\) in m3 m-3.

Each layer has a diagnostic matric water potential (\(\Psi_{lyn}\), in MPa) that is derived as

(2)\[\psi_{lyn}(t) = \psi_{p}\theta_{lyn}(t)^{-b}\]

where \(\psi_{p}\) is set to -1.17 x 10-3 MPa (see Massoud et al. 2022 although we don't use the exact same value) and \(b\) is a unitless retention curve parameter.

Each layer also has a diagnostic temperature (\(T_{lyn}\), in K) and liquid fraction (\(T_{lyn}\), in kg kg-1), the derivation of which is described in section reference energy balance section.

Infiltration into surface

Infiltration of liquid precipitation and snowmelt into the topmost soil layer (\(Qx_{inf}\), in mm d-1) is calculated for each timestep \(t\) as

(3)\[Qx_{inf}(t) = I_{max}\left(1 - \exp\left(-\frac{P_{rain}(t) + Melt(t)}{I_{max}}\right)\right)\]

where \(I_{max}\) is the maximum infiltration parameter, \(P_{rain}\) is the liquid precipitation (\(P_{tot} - P_{snow}\)), and \(Melt\) is snowmelt (defined in equation reference snowmelt equation), all in mm d-1. This formulation ensures \(\approx 100\%\) infiltration when \(P_{rain} + Melt << I_{max}\), and inflitration is \(\approx I_{max}\) \(P_{rain} + Melt >> I_{max}\). Non-infilrated surface water (i.e., \(P_{rain} + Melt - Qx_{inf}\)) is removed as surface runoff (\(Q_{surf}\)).

Inter-layer water transfers

Interlayer water transfers between layers \(n\) and : \(n+1\) (\(Qx_{lyn,n+1}\)) are derived based on the difference in soil water potential (\(\Psi_{lyn}\) and \(\Psi_{lyn+1}\)). The equation used to derive the exchange in mm d:sup:-1 at each time step \(t\) is:

(4)\[Qx_{lyn,n+1}(t) = \rho_l \sqrt{\kappa_{lyn}(t)\kappa_{lyn+1}(t)} \left(\frac{10^{-6}(\Psi_{lyn}(t) - \Psi_{lyn+1}(t))}{\rho_l g 0.5(z_{lyn}+z_{lyn+1})} + 1 \right)\]

where \(\rho_l\) is the density of liquid water (1000 kg m-3), \(g\) is the gravitational acceleration on the surface of Earth (9.8 m s-2), and \(\kappa_{lyn}\) is the hydraulic conductivity of layer \(n\) in m s-1 and is derived as

(5)\[\kappa_{lyn}(t) = \kappa_0\theta_{lyn}(t)^{2b+3}\]

where \(\kappa_0\) is the saturated hydraulic conductivity parameter (in m s-1).

Implementation notes?? We actually calculate potential transfer and check that against the available water in the source layer and available space in the destination layer

Subsurface runoff

Subsurface runoff of layer \(n\) (\(Q_{lyn}\), in mm d-1) is derived based on the difference between the actual soil moisture (\(\theta_{lyn}\)) and the soil moisture at field capacity (\(\theta_{fc}\)), which is derived using the inverse of (2) as

(6)\[\theta_{fc} = \left(\frac{\psi_{p}}{-\psi_{fc}}\right)^{\frac{1}{b}}\]

where \(\psi_{fc}\) is the parameter describing the water potential at field capacity (in -MPa, with the negative sign because CARDAMOM parameters must be positive numbers). The soil moisture difference (\(\Delta\theta_{lyn}\)) is calculated for each time step as

(7)\[\Delta\theta_{lyn}(t) = \max(0,\theta_{lyn}(t) - \theta_{fc})\]

implementation notes: we first make sure :math:`Q_{lyn}` does not exceed 1, and if it does, add any excess to the drainage \(Q_{lyn}\) is then calculated for each time step as

(8)\[Q_{lyn}(t) = \Delta\theta_{lyn}(t) Q_{ex}\left(1 - \frac{\psi_p - \min(\max(\psi_{lyn}(t),-\psi_{fc}),\psi_p))}{\psi_p + \psi_{fc}} \right)\]

References

Massoud, E. C., Bloom, A. A., Longo, M., Reager, J. T., Levine, P. A., and Worden, J. R., 2022. Information content of soil hydrology in a west Amazon watershed as informed by GRACE. Hydrology and Earth System Sciences, 26(5), pp.1407–1423. https://doi.org/10.5194/hess-26-1407-2022