Soil Biogeochemistry: CO₂ Production and O₂ Consumption

This section describes the coupled soil biogeochemistry model implemented in ClimaLand, which simulates the production and diffusion of CO₂, consumption of O₂, and microbial respiration sourced by a fixed pool of soil organic carbon ($Y_{C_{som}}$). The model combines the Dual Arrhenius and Michaelis-Menten (DAMM) kinetics framework for microbial respiration with gas diffusion equations that account for soil structure, temperature, and moisture, and with Henry's-law partitioning of dissolved CO₂ and O₂ between the gas and aqueous phases of the pore space.

Model Overview

The biogeochemistry model carries three prognostic fields that evolve in space (with soil depth) and time. The variables, listed in the units used internally, are:

• the total CO₂ stored per unit soil volume (gas phase plus dissolved), $Y_{CO_2}$ (kg C m⁻³ soil), with $Y_{CO_2} = \theta_{\rm eff}\, c_g$ where $c_g$ is the air-phase carbon mass concentration (kg C m⁻³ air) and $\theta_{\rm eff}$ is the effective porosity, defined as the sum of the gas-filled pore volume and the dissolved-equivalent volume of the liquid pore water (see Effective Porosity below);
• the volumetric O₂ fraction in soil air, $Y_{\theta_{O_2}}$ (dimensionless);
• the soil organic carbon density, $Y_{C_{som}}$ (kg C m⁻³), held fixed at its initial value (no tendency, see below).

The model captures the fundamental processes of aerobic decomposition: microbes consume soluble carbon and oxygen to produce carbon dioxide. CO₂ and O₂ diffuse through the soil pore space, with diffusion rates controlled by soil structure (porosity, tortuosity), moisture, temperature, and pressure.

\[Y_{C_{som}}\]

is initialized from the SoilGrids organic carbon density product and is treated as a stationary spatial parameter: we assume that the input of organic matter (e.g., from litter and root turnover) is equal to the rate of decomposition, so the actual $Y_{C_{som}}$ is constant in time. The model thus captures geographic heterogeneity of carbon stocks but does not integrate carbon depletion on simulation timescales — a deliberate simplification appropriate for the multi-decade runs targeted by ClimaLand. In the future, we may model the rate of decomposition explicitly and track $Y_{C_{som}}$ over time; expressing $Y_{C_{som}}$ as a prognostic variable now allows that change to be made easily.

Governing Equations

CO₂ Transport and Production

The prognostic variable $Y_{CO_2}$ is the total carbon mass stored per unit soil volume, summing the gas phase (occupying the air-filled pore volume $\theta_a$ at concentration $c_g$) and the dissolved phase (held in the liquid pore water $\theta_l$ at concentration $c_{aq}$, related to $c_g$ by Henry's law as $c_{aq} = \beta\, c_g$). The total is therefore

\[\begin{equation} Y_{CO_2} = \theta_a\, c_g + \theta_l\, c_{aq} = (\theta_a + \beta\, \theta_l)\, c_g \equiv \theta_{\rm eff}\, c_g, \end{equation}\]

where $\theta_{\rm eff} \equiv \theta_a + \beta\, \theta_l$ is the effective porosity — the equivalent volume that would hold $Y_{CO_2}$ at concentration $c_g$ if all the carbon were in the gas phase. The gas-phase concentration is recovered as $c_g = Y_{CO_2} / \theta_{\rm eff}$. Storing $Y_{CO_2}$ rather than $c_g$ as the prognostic keeps the formulation well-behaved as the soil saturates ($\theta_a \to 0$): dissolved storage keeps $\theta_{\rm eff} \geq \beta\, \theta_l > 0$ as long as there is any liquid water. Diffusion is driven by gas-phase gradients:

\[\begin{equation} \frac{\partial Y_{CO_2}}{\partial t} = \nabla \cdot \left[D \, \nabla c_g\right] + S_m, \qquad c_g = Y_{CO_2} / \theta_{\rm eff} \end{equation}\]

where:

• the effective diffusivity of CO₂ in soil, $D$ (m² s⁻¹);
• the microbial CO₂ production rate, $S_m$ (kg C m⁻³ s⁻¹).

Here $\theta_a = \nu - \theta_l - \theta_i$ is the volumetric air content, i.e. the fraction of the soil volume occupied by the air phase. Note that we do not currently model fluxes of dissolved CO₂ due to fluxes of the soil water; only gas-phase diffusion is represented.

O₂ Transport and Consumption

Oxygen transport and consumption in soil is described by:

\[\begin{equation} \frac{\partial Y_{\theta_{O_2}}}{\partial t} = \frac{R\, T}{\theta_{\rm eff,O_2}\, P\, M_{O_2}}\, \nabla \cdot \left[D_{O_2}\, \nabla c_{O_2}\right] - \frac{R\, T}{\theta_{\rm eff,O_2}\, P\, M_{C}}\, S_m \end{equation}\]

where:

• the effective porosity for O₂, $\theta_{\rm eff,O_2} = \theta_a + \beta_{O_2}\, \theta_l$ (m³ m⁻³ soil) — same construction as the CO₂ effective porosity defined above, but with the O₂-specific Henry's-law factor. Because O₂ is much less soluble than CO₂ ($\beta_{O_2} \approx 0.032$ vs $\beta_{CO_2} \approx 0.84$ at 298 K), the dissolved-equivalent term contributes much less for O₂;
• the O₂ mass concentration in air, $c_{O_2}$ (kg O₂ m⁻³ air), recovered from $Y_{\theta_{O_2}}$ via the ideal gas law;
• the effective diffusivity of O₂ in soil, $D_{O_2}$ (m² s⁻¹);
• the universal gas constant, $R$ (8.314 J mol⁻¹ K⁻¹);
• the soil temperature, $T$ (K);
• the atmospheric surface pressure, $P$ (Pa);
• the molar mass of O₂, $M_{O_2}$ (0.032 kg mol⁻¹);
• the molar mass of carbon, $M_C$ (0.012 kg mol⁻¹).

The second term represents O₂ consumption by microbial respiration. The stoichiometry follows C + O₂ → CO₂: each 12 g of carbon respired consumes 32 g of O₂. The conversion factor $R T / (\theta_{\rm eff,O_2}\, P\, M_C)$ translates the carbon mass source $S_m$ (kg C m⁻³ s⁻¹ soil) into a tendency on the volumetric fraction $Y_{\theta_{O_2}}$.

Soil Organic Carbon

The $Y_{C_{som}}$ pool is held fixed at its initial value:

\[\begin{equation} \frac{\partial Y_{C_{som}}}{\partial t} = 0 \end{equation}\]

The pool $Y_{C_{som}}$ enters the DAMM substrate kinetics (below) but is not depleted by respiration. Its value is set from the SoilGrids organic carbon density at the start of the simulation, providing the spatial structure of carbon availability.

Microbial Respiration: DAMM Model

The microbial source term $S_m$ is computed using the Dual Arrhenius and Michaelis-Menten (DAMM) kinetics model Davidson et al. [32]:

\[\begin{equation} S_m = V_{\max} \cdot MM_{sx} \cdot MM_{O_2} \end{equation}\]

where:

• the maximum potential rate of respiration (temperature-dependent), $V_{\max}$ (kg C m⁻³ s⁻¹);
• the substrate availability factor, $MM_{sx}$ (dimensionless, 0–1);
• the oxygen limitation factor, $MM_{O_2}$ (dimensionless, 0–1).

Maximum Respiration Rate

The maximum respiration rate follows Arrhenius kinetics, written here in centered form:

\[\begin{equation} V_{\max} = V_{{\rm ref},sx} \, \exp\!\left[-\frac{E_{a,sx}}{R}\!\left(\frac{1}{T} - \frac{1}{T_{{\rm ref},sx}}\right)\right] \end{equation}\]

where:

• the maximum respiration rate at the reference temperature, $V_{{\rm ref},sx}$ (kg C m⁻³ s⁻¹);
• the reference temperature for the centered Arrhenius form, $T_{{\rm ref},sx}$ (K);
• the activation energy, $E_{a,sx}$ (J mol⁻¹).

This is algebraically equivalent to the uncentered form $V_{\max} = \alpha_{sx}\, \exp[-E_{a,sx}/(RT)]$ with $\alpha_{sx} = V_{{\rm ref},sx}\, \exp[E_{a,sx}/(R\, T_{{\rm ref},sx})]$. The centered form is preferred for calibration because $V_{{\rm ref},sx}$ (the rate at $T_{{\rm ref},sx}$) and $E_{a,sx}$ (the temperature sensitivity) are nearly orthogonal under parameter estimation, while $\alpha_{sx}$ and $E_{a,sx}$ are strongly correlated.

Substrate Availability

The concentration of soluble carbon substrate available to microbes depends on diffusion through soil water films:

\[\begin{equation} [S_x] = p_{sx} \cdot Y_{C_{som}} \cdot D_{liq} \cdot \theta_l^{3} \end{equation}\]

where:

• the concentration of soluble substrate, $[S_x]$ (kg C m⁻³);
• the soluble fraction of $Y_{C_{som}}$, $p_{sx}$ (dimensionless);
• the soil organic carbon density, $Y_{C_{som}}$ (kg C m⁻³);
• the dimensionless diffusion coefficient of soluble carbon in water, $D_{liq}$;
• the volumetric liquid water content, $\theta_l$ (m³ m⁻³).

The cubic dependence on moisture reflects the strong constraint that water films place on substrate diffusion to microbial sites.

The substrate limitation factor follows Michaelis-Menten kinetics:

\[\begin{equation} MM_{sx} = \frac{[S_x]}{kM_{sx} + [S_x]} \end{equation}\]

where $kM_{sx}$ is the substrate Michaelis constant (kg C m⁻³).

Oxygen Availability

The oxygen availability for microbial respiration accounts for diffusion limitations in the partially saturated pore space using a Millington–Quirk tortuosity model:

\[\begin{equation} O_{2,{\rm avail}} = D_{oa} \cdot Y_{\theta_{O_2}} \cdot \theta_a^{4/3} \end{equation}\]

where:

• the dimensionless O₂ availability metric, $O_{2,{\rm avail}}$;
• the dimensionless O₂ diffusion coefficient used in the availability scaling, $D_{oa}$;
• the volumetric O₂ fraction in air, $Y_{\theta_{O_2}}$ (dimensionless);
• the Millington–Quirk tortuosity factor, $\theta_a^{4/3}$.

The exponent 4/3 captures how tortuous diffusion pathways limit gas transport in partially saturated porous media.

The oxygen limitation factor follows Michaelis-Menten kinetics:

\[\begin{equation} MM_{O_2} = \frac{O_{2,{\rm avail}}}{kM_{O_2} + O_{2,{\rm avail}}} \end{equation}\]

where $kM_{O_2}$ is the dimensionless O₂ Michaelis constant.

Air–Water Partitioning of Dissolved Gases

CO₂ and O₂ partition between the soil's gas phase and the aqueous phase of pore water. Treating only the gas phase becomes singular as the soil saturates ($\theta_a \to 0$); accounting for the dissolved phase keeps the equations well-posed and is physically more accurate at high water content, especially for CO₂.

Henry's Law

The temperature-dependent Henry's-law solubility is computed using a van 't Hoff form Sander [33]:

\[\begin{equation} K_H(T) = K_H(T_{\rm ref}^{H})\, \exp\!\left[\frac{\partial \ln K_H}{\partial T}\!\left(\frac{1}{T} - \frac{1}{T_{\rm ref}^{H}}\right)\right] \end{equation}\]

where:

• the Henry's-law constant at the Sander reference temperature, $K_H(T_{\rm ref}^{H})$ (mol m⁻³ Pa⁻¹);
• the temperature coefficient (a dimensional analogue of $-\Delta_{\rm sol} H/R$), $\partial \ln K_H / \partial T$ (K);
• the Sander reference temperature, $T_{\rm ref}^{H} = 298.15$ K.

The dimensionless Henry's-law factor $\beta$ is defined by:

\[\begin{equation} \beta(T) = K_H(T)\, R\, T \end{equation}\]

and converts a gas-phase concentration to an equivalent dissolved concentration: $c_{aq} = \beta\, c_g$. At 298 K this gives $\beta_{CO_2} \approx 0.84$ (dissolved CO₂ storage is comparable to gas-phase storage) and $\beta_{O_2} \approx 0.032$ (dissolved O₂ storage is small but stabilizing).

Effective Porosity

The effective porosity used to translate between $Y_{CO_2}$ and $c_g$ — and similarly for the O₂ tendency — is the sum of the gas-filled pore volume and the dissolved-equivalent volume of the liquid pore water:

\[\begin{equation} \theta_{\rm eff} = \max(\theta_a + \beta\, \theta_l,\ \theta_{\rm eff,min}) \end{equation}\]

where the gas-filled pore volume $\theta_a$ and the volumetric liquid water content $\theta_l$ are introduced below. The term $\beta\, \theta_l$ is the "air-equivalent" volume of dissolved gas held in pore water; $\theta_{\rm eff,min} = 10^{-4}$ is a numerical floor that prevents division by zero in extreme conditions where both $\theta_a$ and $\theta_l$ are vanishingly small (see Numerical Safeguards). The same construction with $\beta_{O_2}$ defines $\theta_{\rm eff,O_2}$ used in the O₂ tendency.

Gas Diffusivity in Soil

CO₂ and O₂ Diffusivity

Both CO₂ and O₂ diffusivity in soil are computed using the same parameterization Ryan et al. [34], with separate reference diffusivities for the two gases:

\[\begin{equation} D_{\rm gas} = D_{0,{\rm gas}} \cdot (2\theta_{a100}^{3} + 0.04\, \theta_{a100}) \cdot \left(\frac{\theta_a}{\theta_{a100}}\right)^{2 + 3/b} \end{equation}\]

where:

• the effective diffusivity of the gas in soil, $D_{\rm gas}$ (m² s⁻¹), evaluated for either CO₂ ($D$) or O₂ ($D_{O_2}$);
• the volumetric air content, $\theta_a$ (m³ m⁻³);
• the air-filled porosity at a soil water potential of $-100$ cm H₂O, $\theta_{a100}$ (dimensionless);
• the pore-size distribution parameter, $b$ (dimensionless).

The free-air reference diffusivity is corrected for soil temperature and pressure:

\[\begin{equation} D_{0,{\rm gas}} = D_{{\rm ref},{\rm gas}} \left(\frac{T}{T_{\rm ref}}\right)^{1.75} \left(\frac{P_{\rm ref}}{P}\right) \end{equation}\]

where:

• the gas-specific reference diffusivity at standard conditions, $D_{{\rm ref},{\rm gas}}$ (m² s⁻¹), evaluated as $D_{\rm ref}$ for CO₂ and $D_{{\rm ref},O_2}$ for O₂;
• the standard reference temperature, $T_{\rm ref} = 273$ K;
• the standard reference pressure, $P_{\rm ref} = 101325$ Pa.

The temperature exponent 1.75 reflects increased molecular kinetic energy, while the pressure dependence accounts for the change in gas number density with surface pressure.

Volumetric Air Content

The volumetric air content uses the total water content (liquid plus ice), so frozen water displaces air space:

\[\begin{equation} \theta_a = \nu - \theta_w, \qquad \theta_w = \theta_l + \theta_i \end{equation}\]

where:

• the total soil porosity, $\nu$ (m³ m⁻³);
• the volumetric ice content, $\theta_i$ (m³ m⁻³).

By contrast, the dissolved storage in $\theta_{\rm eff}$ uses only $\theta_l$ — gases dissolve in liquid water, not in ice.

O₂ Concentration Conversions

From Volumetric Fraction to Mass Concentration

The O₂ mass concentration in air is computed from the volumetric fraction using the ideal gas law:

\[\begin{equation} c_{O_2} = Y_{\theta_{O_2}}\, \frac{P\, M_{O_2}}{R\, T} \end{equation}\]

where the O₂ mass concentration in air, $c_{O_2}$ (kg O₂ m⁻³ air), is recovered from the dimensionless volumetric fraction $Y_{\theta_{O_2}}$ together with the surface pressure $P$, the molar mass $M_{O_2} = 0.032$ kg mol⁻¹, the gas constant $R = 8.314$ J mol⁻¹ K⁻¹, and the soil temperature $T$.

From Mass Concentration to Volumetric Fraction

The inverse transformation is:

\[\begin{equation} Y_{\theta_{O_2}} = c_{O_2}\, \frac{R\, T}{P\, M_{O_2}} \end{equation}\]

These conversions are needed because diffusion is most naturally expressed in terms of mass concentration gradients, while the prognostic variable $Y_{\theta_{O_2}}$ and the microbial kinetics use the volumetric fraction.

Boundary Conditions

CO₂ Boundary Conditions

At the top boundary, the diffusive flux drives the soil gas-phase concentration toward the atmospheric value:

\[\begin{equation} F_{CO_2,{\rm top}} = -D\, \frac{c_{g,{\rm atm}} - c_{g,{\rm top}}}{\Delta z_{\rm top}}, \qquad c_{g,{\rm atm}} = \chi_{CO_2}\, P\, M_C / (R\, T_{\rm top}) \end{equation}\]

where:

• the atmospheric CO₂ mole fraction (mol mol⁻¹), $\chi_{CO_2}$, supplied as a driver;
• the gas-phase carbon concentration in the top soil layer, $c_{g,{\rm top}}$ (kg C m⁻³ air);
• the distance from the top layer center to the surface, $\Delta z_{\rm top}$ (m);
• the soil temperature at the top of the column, $T_{\rm top}$ (K).

At the bottom boundary, a zero-flux condition is applied: $F_{CO_2,{\rm bottom}} = 0$.

O₂ Boundary Conditions

At the top boundary, the atmospheric O₂ volumetric fraction $\theta_{O_2,{\rm atm}}$ ($\approx 0.21$) is converted to a mass concentration via the ideal gas law and used as the boundary value:

\[\begin{equation} F_{O_2,{\rm top}} = -D_{O_2}\, \frac{c_{O_2,{\rm atm}} - c_{O_2,{\rm top}}}{\Delta z_{\rm top}}, \qquad c_{O_2,{\rm atm}} = \theta_{O_2,{\rm atm}}\, P\, M_{O_2} / (R\, T_{\rm top}). \end{equation}\]

At the bottom boundary, a zero-flux condition is applied: $F_{O_2,{\rm bottom}} = 0$.

Numerical Safeguards

The continuous PDE above is integrated with an explicit scheme. Several pragmatic safeguards keep the discrete update well-behaved in extreme regimes (e.g. very dry hot columns or saturated ice-rich columns), without changing the equations in their physical regime:

• the floor on the effective porosity, $\theta_{\rm eff,min} = 10^{-4}$, prevents a singular $c_g = Y_{CO_2}/\theta_{\rm eff}$ when both $\theta_a$ and $\theta_l$ are vanishingly small (saturated, frozen soil);
• the diffusivity ratio is capped at $\theta_a / \theta_{a100} \le 5$, preventing the Ryan/Moldrup expression from producing unphysically large $D$ in soils much drier than their reference state ($\theta_a \gg \theta_{a100}$);
• the per-step bracket on the $Y_{\theta_{O_2}}$ tendency clamps $\partial Y_{\theta_{O_2}}/\partial t$ so that one explicit step cannot push $Y_{\theta_{O_2}}$ outside $[0,\ \theta_{O_2,{\rm atm}}]$; this kicks in only when the explicit step would otherwise violate the diffusion CFL bound, and is conservative for multi-stage Runge-Kutta integrators;
• an extra Michaelis-Menten attenuation $Y_{\theta_{O_2}} / (Y_{\theta_{O_2}} + K_{O_2,{\rm lim}})$ with $K_{O_2,{\rm lim}} = 10^{-4}$ multiplies the O₂ consumption term, which together with the tendency bracket prevents the explicit step from driving $Y_{\theta_{O_2}}$ negative in the rare cells where $Y_{\theta_{O_2}}$ becomes very small.

Summary Tables

Model Outputs

Drivers

Parameters

Constants