Quantifying aggregation bias in marsh carbon flux estimates caused by rhizosphere oxygen heterogeneity

Oxygen (O₂) in marsh rhizospheres varies strongly at millimetre scales due to radial oxygen loss (ROL), re-oxidation of reduced compounds, hydrology, and sediment structure. Yet many models estimate heterotrophic (microbial) respiration at a single mean O₂, which can bias CO₂ fluxes because the O₂ response of respiration is curved (nonlinear). We calibrated an enzyme kinetic respiration function to marsh-soil incubations, then compared (I) respiration evaluated at a single mean O₂ value with (II) the mean of different respiration fluxes over a prescribed rhizosphere O₂ distribution. To isolate oxygen effects, temperature and the relation between substrate and soil moisture were held fixed; incubations were used only to calibrate carbon kinetics. For a gas-phase, air-equivalent O₂ proxy, the single mean approach overestimates respiration by ∼ 12 % on average. This bias increases with increasing variation of O₂ and with the ratio of mean O₂ to the O₂ half-saturation constant $K_{\mathrm{M},{\mathrm{O}}_{\mathrm{2}}}$. In the dissolved phase, mean O₂ ≪ $K_{\mathrm{M},{\mathrm{O}}_{\mathrm{2}}}$ and variation in O₂ is small, so the bias is negligible due to linear behaviour of the relation in these conditions; when physiology is unit-matched to a similar relative saturation, however, the dissolved case shows a negative bias comparable to the gas-phase one. Our findings suggest reporting the mean and variation of O₂ together with a calibrated $K_{\mathrm{M},{\mathrm{O}}_{\mathrm{2}}}$ allows modellers to anticipate or correct the bias associated with replacing heterogeneous rhizosphere conditions with a single mean.

1 Introduction

Marsh soils are major coastal carbon reservoirs where high primary production and efficient sediment trapping coincide with comparatively slow decomposition (Barbier et al., 2011; Duarte et al., 2013; Granse et al., 2024; Kirwan et al., 2016). Whether these systems store or release carbon depends not only on mean environmental conditions but also on the availability of oxygen (O₂) to microbes in the rhizosphere, which regulates aerobic heterotrophic respiration.

Rhizosphere O₂ is inherently patchy at millimetre scales. Radial oxygen loss (ROL) from roots varies across species and genotypes and with plant status (growth stage, stress), while aerenchyma formation and barriers to ROL modulate internal gas pathways. Hydrologic dynamics, such as tidal flooding and drainage, platform elevation, and episodic drying interact with sediment texture to shape air-filled porosity and diffusion pathways. O₂ that enters the root zone is, in addition, consumed by the re-oxidation of reduced compounds such as dissolved sulfide, iron sulfides and ammonium, which competes with aerobic respiration for O₂ and confines oxic zones to the immediate vicinity of root surfaces (Howes and Teal, 1994; Koop-Jakobsen et al., 2017). In situ microelectrode and optical measurements consistently show steep radial O₂ gradients around roots and strong small scale variability within and across root neighbourhoods (Kühl et al., 2021; Mueller et al., 2016; Van der Nat and Middelburg, 1998; Visser et al., 2000). Thus, microbes experience a mosaic of microhabitats ranging from well oxygenated to effectively anoxic over centimetres or less. In such a setting, the curved, nonlinear shape of the O₂ response of respiration may lead to differing estimates depending on the level of aggregation of O₂ conditions.

Despite this heterogeneity, many ecosystem and land-surface models represent belowground redox drivers with a single bulk state, for example a mean O₂ level or a uniform moisture or redox scalar within a grid cell or plant functional type (Alizad et al., 2016; Kirwan et al., 2016; Morris et al., 2002; Rietl et al., 2021; Swanson et al., 2014; Zhang et al., 2002). Under a concave O₂-limitation function (saturating, Michaelis–Menten-type), evaluating respiration at the mean O₂ generally differs from averaging respiration across a heterogeneous O₂ field, a classic aggregation (upscaling) bias arising from response curvature (Jensen's inequality) (Rastetter et al., 1992). For marsh soils, strong O₂ heterogeneity is well documented and concave limitation is physiologically expected, yet the magnitude, direction, and conditions of the resulting bias remain poorly quantified (Bernal et al., 2017; Bhattacharyya and Furtak, 2022; Mueller et al., 2020; Spivak et al., 2019).

Here we quantify the aggregation bias associated with spatially heterogeneous rhizosphere O₂, focusing explicitly on aerobic heterotrophic (microbial) respiration in the root zone. While soil respiration includes both autotrophic (root) and heterotrophic (microbial) components, we focus on the heterotrophic component since the aggregation bias arises from the nonlinear O₂ limitation of microbial respiration. The drivers of autotrophic respiration, including climatic conditions and plant physiology, are outside our scope. Moreover, we do not model explicitly plant internal gas transport, or pore scale advection-diffusion. Instead, we treat ROL and hydrology as drivers that shape local O₂ distributions and ask how that heterogeneity propagates through a standard enzyme kinetic respiration formulation with an explicit O₂-limitation term (cf. Davidson et al., 2012 for the general structure). To isolate oxygen effects, we calibrate carbon kinetic parameters to marsh soil incubations and hold temperature and the substrate–moisture linkage fixed.

We represent rhizosphere O₂ using two complementary proxies: (I) a gas-phase, air-equivalent proxy that captures aeration via air-filled porosity; and (II) a dissolved-phase proxy derived from Henry's law. For like-for-like comparisons, we also analyse a diagnostic in which the microbial O₂ half-saturation constant $K_{\mathrm{M},{\mathrm{O}}_{\mathrm{2}}}$ is unit-matched so that the relative saturation r = mean(O₂) $/$ $K_{\mathrm{M},{\mathrm{O}}_{\mathrm{2}}}$ equals that of the gas case. We define the uniform baseline in O₂ space as f(mean(O₂)) (rather than mapping mean moisture to O₂), then quantify bias as the difference between this uniform evaluation and the mean of respiration over the full O₂ distribution.

We address three questions. (Q1) How large, and in which direction, is the bias when O₂ heterogeneity is replaced by a single mean? (Q2) How does the bias scale with relative saturation and with the variance of O₂ at fixed mean? (Q3) Under what conditions is the bias negligible, providing a simple rule for when single mean treatments are acceptable? Answering these provides a practical diagnostic for modellers: reporting the mean and variance of O₂ together with a prescribed $K_{\mathrm{M},{\mathrm{O}}_{\mathrm{2}}}$ is sufficient to anticipate, and if needed correct, aggregation error in marsh carbon balance models.

2 Materials and methods

2.1 Model overview: DAMM with explicit O₂ limitation

We use the DAMM framework (Davidson et al., 2012) in which volumetric aerobic CO₂ production R_sx (mg C cm⁻³ h⁻¹) depends on temperature, soluble-carbon substrate, and oxygen:

$$\begin{array}{}\text{(1)}& R_{\mathrm{s}x}\left(C_{\mathrm{org}},{\mathrm{O}}_{\mathrm{2}},T\right)=\mathit{\alpha }{e}^{\frac{-E_{\mathrm{a}}}{R_{\mathrm{gas}}T}}{\displaystyle \frac{S_x}{K_{\mathrm{M},\mathrm{S}}+S_x}}{\displaystyle \frac{\left[{\mathrm{O}}_{\mathrm{2}}\right]}{K_{\mathrm{M},{\mathrm{O}}_{\mathrm{2}}}+\left[{\mathrm{O}}_{\mathrm{2}}\right]}}\end{array}$$

Here C_org is bulk organic C (g C cm⁻³); S_x is the soluble substrate concentration (g C cm⁻³); O₂ is an air-equivalent oxygen proxy expressed as cm³ O₂ cm⁻³; T is absolute temperature (K); α (mg C cm⁻³ h⁻¹) and E_a (J mol⁻¹) are the pre-exponential factor and the activation energy; K_M,S (g C cm⁻³) and $K_{\mathrm{M},{\mathrm{O}}_{\mathrm{2}}}$ (cm³ O₂ cm⁻³) are half-saturation constants; and R_gas is the universal gas constant. Symbols, units, and baseline values are listed in Table 1.

To isolate O₂ effects in the aggregation experiment, soluble carbon is linked to bulk carbon via a fixed moisture proxy:

$$\begin{array}{}\text{(2)}& S_x=C_{\mathrm{org}}p{D}_{\mathrm{liq}}{{\mathit{\theta }}_{\mathrm{ref}}}^{\mathrm{3}}\end{array}$$

where p (–) converts bulk to soluble carbon, D_liq (–) scales aqueous availability, and θ_ref is the volumetric water content, used only inside the substrate term (to isolate oxygen effects). The oxygen-limitation term O₂ $/$ $(K_{\mathrm{M},{\mathrm{O}}_{\mathrm{2}}}+{\mathrm{O}}_{\mathrm{2}})$ is concave for O₂ > 0; this curvature is the source of aggregation bias examined below.

2.2 Incubation medians and calibration

We calibrated α and K_M,S against median aerobic respiration rates derived from long-term marsh-soil incubations summarized with a two-pool introductory carbon balance model (ICBM; Andrén and Kätterer, 1997; Beer et al., 2022; Knoblauch et al., 2013). Incubation rates were converted from µmol CO₂ g⁻¹ d⁻¹ to mg C cm⁻³ h⁻¹ using soil bulk density and molar-mass factors. Parameters were estimated by bounded constrained optimization (SQP; MATLAB fmincon) in log ₁₀ space with box constraints, minimizing the RMSE between model predictions and the binned medians, at a fixed reference oxygen O_2,ref and a fixed water content θ_ref. The reference oxygen was computed from the gas-phase proxy defined in Sect. 2.3 evaluated at θ_ref. Uncertainty was quantified by non-parametric bootstrap resampling of the (C_org, R) pairs (B = 500). We report 95 % pointwise bands for the calibration curve, confidence intervals for α and K_M,S, and the Spearman correlation ρ between log ₁₀α and log ₁₀K_M,S.

2.2.1 Incubation experiments

Marsh soils from the Elbe estuary (Germany) were gently sieved to 2 mm. For each incubation, 20 g dry-mass equivalent soil was weighed (wet masses corrected to dry-mass equivalence). The samples were then adjusted to 60 % water-holding capacity and placed in 1000 mL glass flasks. Incubations were conducted in the dark at 20 °C for up to 465 d. At regular intervals, 1 mL headspace samples were withdrawn and analyzed for CO₂ by gas chromatography (JAS 6890N, Germany). To maintain aerobic conditions, we kept a slight overpressure and added synthetic air whenever headspace CO₂ exceeded 2.5 %, preventing O₂ depletion. In the context of our study, these incubations provide aerobic respiration rates across a broad C_org range for Elbe estuary marsh soils; they do not impose or resolve O₂ heterogeneity. We use the rates only to calibrate carbon-kinetic parameters; O₂ effects are controlled separately in the aggregation analyses (Sect. 2.3 to 2.6).

2.2.2 From CO₂ time series to characteristic rates (ICBM fits)

CO₂ time series were fit with the two-pool ICBM (fast and slow pools, first-order kinetics, humification from fast to slow) using gradient-descent optimization to estimate pool decay and humification parameters against the gas chromatography data. For each of seven C_org bins, we computed the median CO₂ rate (µmol CO₂ g⁻¹ d⁻¹), then converted it to mg C cm⁻³ h⁻¹ as above. These bin medians are the only empirical constraints used to fit α and K_M,S at fixed O_2,ref and θ_ref, thereby isolating oxygen effects in subsequent experiments.

2.3 Prescribed rhizosphere O₂ distributions

We prescribe two rhizosphere O₂ drivers and summarize each by its mean and variance.

For the gas-phase proxy (primary analysis), air-filled porosity a is approximated as

$$\begin{array}{}\text{(3)}& a=\mathrm{1}-{\displaystyle \frac{\text{BD}}{\text{PD}}}-\mathit{\theta }\end{array}$$

where θ is volumetric water content and BD and PD are bulk and particle densities (g cm⁻³), respectively. Air-equivalent oxygen is then

$$\begin{array}{}\text{(4)}& {{\mathrm{O}}_{\mathrm{2}}}^{\mathrm{gas}}=D_{\mathrm{gas}}\times \mathrm{0.209}\times a^{\frac{\mathrm{4}}{\mathrm{3}}}\end{array}$$

where D_gas (–) represents effective gas-phase diffusivity and 0.209 the atmospheric O₂ mole fraction, which for an ideal gas is numerically equal to its volume fraction. Because Eq. (4) multiplies this fraction by the air-filled porosity, O${}_{\mathrm{2}}^{\mathrm{gas}}$ is an amount of O₂ per unit volume of soil (cm³ O₂ cm⁻³), the same unit as $K_{\mathrm{M},{\mathrm{O}}_{\mathrm{2}}}$. This proxy captures aeration in the rhizosphere and is used for all main analyses.

For the aqueous-phase proxy, dissolved O₂ is derived from Henry's law:

$$\begin{array}{}\text{(5)}& {{\mathrm{O}}_{\mathrm{2}}}^{\mathrm{aq}}=\mathit{\theta }{D}_{{\mathrm{O}}_{\mathrm{2}}}{C}_{{\mathrm{O}}_{\mathrm{2}}}^{*}{V}_{\mathrm{m}}\end{array}$$

where $D_{{\mathrm{O}}_{\mathrm{2}}}$ is the diffusivity of O₂ in water, written here as a dimensionless scaling factor, consistently with D_gas in Eq. (4), and $C_{{\mathrm{O}}_{\mathrm{2}}}^{*}$ = H_cp $P_{{\mathrm{O}}_{\mathrm{2}},\mathrm{atm}}$ is the Henry's-law equilibrium concentration (mol m⁻³) at atmospheric partial pressure. This concentration is converted using the molar volume V_m, so that O${}_{\mathrm{2}}^{\mathrm{aq}}$ has units of volume fraction (cm³ O₂ cm⁻³), the same as O${}_{\mathrm{2}}^{\mathrm{gas}}$. We do not rescale dissolved O₂; we retain its native magnitude relative to $K_{\mathrm{M},{\mathrm{O}}_{\mathrm{2}}}$.

In a diagnostic used for Sect. 3.4 and Appendix A, we keep the dissolved proxy unaltered and instead unit-match the microbial half-saturation so that relative saturation is identical across phases:

$$\begin{array}{}\text{(6)}& K_{\mathrm{M},{\mathrm{O}}_{\mathrm{2}}}^{\mathrm{aq},\mathrm{UM}}=K_{\mathrm{M},{\mathrm{O}}_{\mathrm{2}}}^{\mathrm{gas}}{\displaystyle \frac{\mathit{\mu }\left({{\mathrm{O}}_{\mathrm{2}}}^{\mathrm{aq}}\right)}{\mathit{\mu }\left({{\mathrm{O}}_{\mathrm{2}}}^{\mathrm{gas}}\right)}}\end{array}$$

This aligns ${\mathit{\mu }}_{{\mathrm{O}}_{\mathrm{2}}}$ $/$ $K_{\mathrm{M},{\mathrm{O}}_{\mathrm{2}}}$ between gas and dissolved cases without artificially rescaling O₂.

We map a uniform grid of θ values to O₂ using Eqs. (3) to (5) to obtain discrete O₂ samples. In the main text we use θ ∈ $[\mathrm{0.02},\mathrm{1}-\text{BD}/\text{PD}-\mathrm{0.02}]$; as a high-θ benchmark, we also consider a narrow saturated θ ∈ [0.65,0.68] (both mapped through the same proxy). Each driver is summarized by its first two moments,

$$\begin{array}{}\text{(7)}& {\mathit{\mu }}_{{\mathrm{O}}_{\mathrm{2}}}=E\left[{\mathrm{O}}_{\mathrm{2}}\right],\text{Var}\left({\mathrm{O}}_{\mathrm{2}}\right)=E\left[{\left({\mathrm{O}}_{\mathrm{2}}-{\mathit{\mu }}_{{\mathrm{O}}_{\mathrm{2}}}\right)}^{\mathrm{2}}\right]\end{array}$$

which are used later in the Jensen diagnostic.

2.4 Uniform versus heterogeneous experiment (aggregation error)

To isolate the oxygen effect, the substrate S_x is held fixed by evaluating it at a single θ_ref in all runs. For each C_org on a prescribed grid we compute two configurations of the same calibrated DAMM model. Uniform configuration: R_uni(C_org) = $f(C_{\mathrm{org}},{\mathit{\mu }}_{{\mathrm{O}}_{\mathrm{2}}})$, where ${\mathit{\mu }}_{{\mathrm{O}}_{\mathrm{2}}}$ = E[O₂] is the arithmetic mean of the imposed O₂ distribution. Heterogeneous configuration: R_het(C_org) = E[f(C_org,O₂)], that is, the arithmetic mean of the respiration response over the same O₂ distribution, with C_org and θ_ref held fixed.

Aggregation bias (percent) is calculated as

$$\begin{array}{}\text{(8)}& {\mathrm{\Delta }}_{\mathrm{s}}\left(C_{\mathrm{org}}\right)={\displaystyle \frac{R_{\mathrm{het}}\left(C_{\mathrm{org}}\right)-R_{\mathrm{uni}}\left(C_{\mathrm{org}}\right)}{max\left(R_{\mathrm{uni}}\left(C_{\mathrm{org}}\right),\mathit{ϵ}\right)}}\times \mathrm{100}\mathit{\%}\end{array}$$

with ϵ = 10⁻¹² mg C cm⁻³ h⁻¹ used only to avoid division by zero (in practice R_uni ≫ ϵ). In the implementation, the expectation E[⋅] is a discrete average over the vector of O₂ values obtained by mapping the θ grid through the chosen proxy.

Crucially, the uniform baseline is defined in O₂ space as $f(C_{\mathrm{org}},{\mathit{\mu }}_{{\mathrm{O}}_{\mathrm{2}}})$, to avoid baseline reversal when the θ to O₂ mapping is nonlinear. Unless stated otherwise, main text analyses use the gas-phase driver; dissolved-phase analyses appear as sensitivity/diagnostics, and where indicated also with unit-matched $K_{\mathrm{M},{\mathrm{O}}_{\mathrm{2}}}$ to equalize relative saturation ${\mathit{\mu }}_{{\mathrm{O}}_{\mathrm{2}}}/K_{\mathrm{M},{\mathrm{O}}_{\mathrm{2}}}$.

2.5 Dependence on O₂ affinity and transport distributions

Since $K_{\mathrm{M},{\mathrm{O}}_{\mathrm{2}}}$ sets the relative O₂-saturation regime, we performed a structured $K_{\mathrm{M},{\mathrm{O}}_{\mathrm{2}}}$ variation (multipliers × 0.25, × 0.5, × 1, × 2, × 4) for the gas proxy and for the aqueous proxy with unit-matched $K_{\mathrm{M},{\mathrm{O}}_{\mathrm{2}}}$. We further summarize bias behaviour on a response surface as a function of mean$\left({\mathrm{O}}_{\mathrm{2}}\right)/K_{\mathrm{M},{\mathrm{O}}_{\mathrm{2}}}$ and the coefficient of variation (CV) using a symmetric two-point distribution at fixed mean (reported in Appendix A).

2.6 Jensen diagnostic (mechanistic check)

To link bias magnitude to heterogeneity directly, we fixed the mean oxygen at μ and constructed symmetric two-point mixtures with values μ − d and μ + d (equal probability; range constrained to positive values) to vary Var[O₂] while holding the mean constant. For each variance level we computed Δ from Eq. (8) and compared it to the quadratic approximation implied by the curvature of the O₂-limitation term:

$$\begin{array}{}\text{(9)}& \mathrm{\Delta }\left(d\right)\approx {\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\displaystyle \frac{f^{\prime \prime }\left(\mathit{\mu }\right)\text{Var}\left[{\mathrm{O}}_{\mathrm{2}}\left(d\right)\right]}{f\left(\mathit{\mu }\right)}}\times \mathrm{100}\mathit{\%}\end{array}$$

where f^′′(μ) is the second derivative of f with respect to O₂ at fixed C_org and T.

2.7 Implementation details

All analyses were performed with consistent grid resolution for C_org (300 points across 0.001–0.30 g C cm⁻³). Calibration used bounded least squares; bootstrap size was 500; a fixed random seed was used for reproducibility.

3 Results

3.1 Model calibration and parameter uncertainty

The respiration formulation calibrated well against incubation medians across the observed range of C_org (Fig. 1). Best-fit parameters are α = 1.48 × 10¹¹ mg C cm⁻³ h⁻¹ (95 % CI: 1.04 × 10¹¹–9.93 × 10¹³) and K_M,S = 5.98 × 10⁻⁶ g C cm⁻³ (RMSE = 6.04 × 10⁻⁴ mg C cm⁻³ h⁻¹), bootstrap band (B = 500) around the calibration curve widens at high C_org where constraints are sparse. A few medians fall outside the ribbon because it reflects uncertainty in the fitted curve, not confidence intervals on binned medians. Incubations were aerobic at 20 °C and span the observed C_org range; these medians calibrate the carbon-kinetic terms only. In all aggregation experiments, temperature and the substrate–moisture linkage are held fixed to isolate oxygen effects; only the treatment of O₂ variability differs.

Figure 1 Calibration of the DAMM respiration model to incubation-derived respiration medians under aerobic conditions. Points show median rates for seven bins of bulk organic carbon (C_org) obtained from a two-pool ICBM fitted to Elbe estuary incubations (units converted to mg C cm⁻³ h⁻¹). The solid line is the best DAMM fit evaluated at the air-equivalent O₂ corresponding to θ_ref = 0.30; the shaded envelope is the 95 % bootstrap band (B = 500). Best-fit parameters: α = 1.48 × 10¹¹ mg C cm⁻³ h⁻¹, K_M,S = 5.98 × 10⁻⁶ g C cm⁻³; RMSE = 6.04 × 10⁻⁴ mg C cm⁻³ h⁻¹. Parameter confidence intervals and provenance are summarized in Table 1.

Table 1 Parameters used in the DAMM analyses (fixed and fitted). All constants are baseline values unless noted; fitted parameters carry 95 % bootstrap confidence intervals.

3.2 Rhizosphere O₂ distributions

Both prescribed O₂ drivers are right-tailed on a log scale (Fig. 2). Because the dissolved proxy is about 450 times smaller in magnitude than the gas proxy, we show them in separate panels to preserve visibility. For the gas-phase proxy [air-equivalent units], the mean and variance are μ_gas = 0.0891 cm³ O₂ cm⁻³ and Var_gas = 3.53 × 10⁻³ (cm³ O₂ cm⁻³)². For the dissolved proxy [volume fraction], μ_aq = 1.97 × 10⁻⁴ cm³ O₂ cm⁻³ and Var_aq = 1.16 × 10⁻⁸ (cm³ O₂ cm⁻³)² (Fig. A1). We retain each proxy in its native magnitude (no cross-rescaling) so its position relative to $K_{\mathrm{M},{\mathrm{O}}_{\mathrm{2}}}$ is preserved; these moments are used in the Jensen diagnostic (Sect. 3.5). A high-θ benchmark narrows the gas-phase distribution by mapping θ ∈ [0.65,0.68] through the same proxy.

Figure 2 Prescribed rhizosphere O₂ distributions used to drive the DAMM model. (a) Gas-phase proxy in air-equivalent units (cm³ O₂ cm⁻³), obtained from air-filled porosity and an effective gas diffusivity. (b) Dissolved-phase proxy based on Henry's law, expressed as a volume fraction (cm³ O₂ cm⁻³). Both panels show probability density functions on a logarithmic y axis; in-panel annotations report the mean (μ) and variance (Var) of each distribution, which are later used in the Jensen diagnostic.

3.3 Aggregation bias: uniform versus heterogeneous O₂

Across the full C_org range, evaluating respiration at a single mean O₂ (uniform) yields higher fluxes than averaging the model over the full O₂ distribution (heterogeneous) for the gas-phase driver (Fig. 3). The mean aggregation bias is −12 %, i.e., using a single mean O₂ overestimates aerobic heterotrophic respiration by about one tenth on average under observed gas-phase heterogeneity. For the unscaled dissolved proxy, μ_aq ≪ $K_{\mathrm{M},{\mathrm{O}}_{\mathrm{2}}}$ and Var_aq is small, so the response is near-linear and the mean bias is negligible (−0.05 %; Fig. A2). Tightening the gas-phase distribution in the high-θ benchmark reduces variance and yields a small bias (−0.5 % at C_org = 0.15 g C cm⁻³; Fig. A5).

The ∼ 12 % overestimate is model-structural (Jensen concavity) rather than statistical: the comparison uses the same calibrated model in both configurations and differs only by how O₂ heterogeneity is treated. Hence, the bias is operationally meaningful for budgeting and benchmarking aerobic CO₂ fluxes.

Figure 3 Aerobic respiration as a function of C_org when evaluated at a single mean O₂ (“Uniform”) versus averaged over the full gas-phase O₂ distribution (“Heterogeneous”). Both curves use the calibrated DAMM parameters and the gas-phase O₂ proxy. The reported mean aggregation bias (about 12 %) quantifies the relative difference between the heterogeneous and uniform cases over the C_org range considered.

3.4 Dependence on O₂ affinity and relative saturation

To separate transport from microbial physiology, we varied $K_{\mathrm{M},{\mathrm{O}}_{\mathrm{2}}}$ while holding the O₂ distributions fixed and plotted the aggregation bias against the relative saturation ratio mean$\left({\mathrm{O}}_{\mathrm{2}}\right)/K_{\mathrm{M},{\mathrm{O}}_{\mathrm{2}}}$, which collapses cases onto an interpretable axis (Fig. 4). For the gas proxy, varying $K_{\mathrm{M},{\mathrm{O}}_{\mathrm{2}}}$ by factors of × 0.25, × 0.5, × 1, × 2, × 4 yields a bias range from −13 % to −6 %. For the dissolved proxy with unit-matched $K_{\mathrm{M},{\mathrm{O}}_{\mathrm{2}}}$ (that is, chosen so that its mean O₂ concentration corresponds to the same relative saturation as the gas case), the range is −9 % to −4 %. The baseline points occur at mean$\left({\mathrm{O}}_{\mathrm{2}}\right)/K_{\mathrm{M},{\mathrm{O}}_{\mathrm{2}}}$ ≈ 0.74 for both proxies; biases differ in magnitude because the underlying O₂ PDFs differ in spread.

In short, the magnitude of bias increases as the system approaches the nonlinear mid-saturation regime and as O₂ variability increases, and it does not decrease again over the range of relative saturation examined here.

Figure 4 Aggregation bias (%) as a function of mean(O₂) $/$ $K_{\mathrm{M},{\mathrm{O}}_{\mathrm{2}}}$ for the gas-phase O₂ proxy (circles) and the dissolved-phase proxy with unit-matched $K_{\mathrm{M},{\mathrm{O}}_{\mathrm{2}}}$ (squares). For each proxy, $K_{\mathrm{M},{\mathrm{O}}_{\mathrm{2}}}$ is varied via multiplicative factors around the calibrated gas-phase value. Bias is computed as the percent difference between heterogeneous (averaged over the O₂ distribution) and uniform (evaluated at the mean O₂) respiration at fixed C_org. The bias is most pronounced at intermediate relative saturation (mean(O₂) $/$ $K_{\mathrm{M},{\mathrm{O}}_{\mathrm{2}}}$ ≈ 1). Over the range shown it becomes smaller towards low relative saturation, whereas at high relative saturation it remains almost constant, because the prescribed O₂ distributions retain appreciable density close to zero.

3.5 Jensen diagnostic

Holding the mean O₂ fixed and increasing Var(O₂) with symmetric two-point mixtures produces a near-linear growth in bias (Fig. A3). A second order (quadratic) approximation based on the local curvature of the Michaelis–Menten O₂ term closely tracks the simulations, confirming that the bias is a structural consequence of concavity: averaging a heterogeneous O₂ field yields E[f(O₂)] < f(E[O₂]). A response-surface view in Appendix A (Fig. A6) shows that bias remains non-positive across the explored domain and peaks when mean(O₂) is comparable to $K_{\mathrm{M},{\mathrm{O}}_{\mathrm{2}}}$ and the coefficient of variation is large.

4 Discussion

4.1 Principal findings

Our analysis shows that spatial heterogeneity in rhizosphere oxygen can systematically bias modelled aerobic heterotrophic respiration when the process is evaluated at a single mean oxygen concentration. Using a Dual Arrhenius–Michaelis–Menten (DAMM) scheme calibrated to marsh-soil incubations, we compared a “uniform” configuration, in which respiration is evaluated at the spatial mean oxygen, with a “heterogeneous” configuration, where respiration is averaged over a prescribed distribution of oxygen in the root zone. For gas-phase (air-equivalent) oxygen distributions representative of tidally influenced marsh soils, the uniform configuration overestimates aerobic CO₂ production by roughly 10 %–15 % across a wide range of organic carbon concentrations. This aggregation bias arises because the oxygen-limitation term in the DAMM model is a saturating Michaelis–Menten function of oxygen, which is concave over the range where respiration is partially oxygen-limited.

The magnitude of the bias depends on two simple descriptors of the oxygen field: the ratio of the mean oxygen concentration to the microbial half-saturation constant for oxygen, and the variance of oxygen in the rhizosphere. When mean oxygen is far below the half-saturation constant, microbial respiration responds nearly linearly to oxygen and the bias is negligible. When mean oxygen is far above the half-saturation constant, the oxygen term is close to saturation and again becomes insensitive to spatial variation. The largest aggregation errors occur at intermediate mean oxygen, where respiration is strongly limited but not saturated, and where the variance in oxygen is substantial. Our two point and lognormal variance experiments, which hold the mean oxygen fixed while increasing its variance, confirm that the magnitude of the negative bias increases monotonically with oxygen variance and that the sign of the bias does not change within the Michaelis–Menten framework (Figs. A3, A4).

The comparison between gas-phase and dissolved-phase oxygen drivers illustrates the same mechanism from a complementary angle. When the dissolved-phase proxy is used without rescaling, its mean and variance are both much smaller than the half-saturation constant, so respiration is effectively in the linear regime of the oxygen response, and the aggregation bias is close to zero (Fig. A2). When we adjust the microbial half-saturation constant so that the dissolved-phase proxy has a similar degree of oxygen limitation to the gas-phase proxy (that is, a similar ratio of mean oxygen to half-saturation), the dissolved-phase simulations exhibit biases comparable in sign and magnitude to the gas-phase case (Fig. 4). The difference between the gas and aqueous results therefore reflects how each proxy is positioned relative to microbial oxygen affinity, rather than a change in model structure.

Overall, the marsh-soil case study provides a concrete example of a general mechanism: when a saturating oxygen-limitation term is applied to a spatially heterogeneous rhizosphere, evaluating respiration at a single bulk oxygen level tends to overestimate aerobic heterotrophic CO₂ production. The size of that overestimate can be understood and, in principle, approximated from simple statistics of the oxygen distribution and plausible ranges of the half-saturation constant.

4.2 Relation to process understanding and prior work

Tidal marsh soils are characterized by strong small-scale redox heterogeneity. Hydrologic forcing through tidal inundation, drainage, and groundwater flow drives rapid switches between oxic and anoxic conditions, while plant traits such as radial oxygen loss from aerenchymatous roots create localized oxic zones around roots embedded in largely reduced sediment. Microelectrode and planar optode studies in flooded and marsh soils consistently show millimetre-centimetre oxygen mosaics, with steep gradients around roots and patches where oxygen fluctuates over diel and tidal cycles (Koop-Jakobsen et al., 2018). Our imposed gas-phase oxygen distributions are designed to be consistent with this process understanding: they arise from variation in water content and air-filled porosity, both of which are strongly affected by tidal inundation and drainage, and they are intended to bracket realistic degrees of rhizosphere oxygen heterogeneity in tidally influenced marshes.

Previous work has emphasized that tidal marshes can be strong carbon sinks because high primary production and efficient sediment trapping are combined with slow anaerobic decomposition. At the same time, a growing literature highlights how root-mediated oxygen supply can locally accelerate aerobic decomposition and alter greenhouse gas production, especially in systems with high root biomass and pronounced aerenchyma development. However, most marsh carbon models represent belowground redox conditions using a single bulk variable (such as average water content, a redox class, or an empirical soil-moisture scalar), and they typically apply this scalar uniformly within a grid cell or plant functional type. In doing so, they effectively assume that respiration responds to a uniform oxygen environment, even though observations show that roots and microbes experience a patchwork of oxic and anoxic microhabitats (Pezeshki and DeLaune, 2012).

The broader ecological modelling literature has long recognized that using mean environmental conditions in nonlinear response functions can lead to systematic upscaling biases. Rastetter et al. (1992), for example, showed that aggregating nonlinear responses to water and nutrients can misrepresent ecosystem responses to climate change. Yet, to our knowledge, the specific consequences of small-scale rhizosphere oxygen heterogeneity for marsh carbon flux estimates have not been quantified using a physiologically based respiration scheme calibrated to marsh soils. Our results therefore connect a well established mathematical principle (Jensen's inequality) to a specific, process rich context in which oxygen heterogeneity is both expected and extensively documented. They complement previous work on micro-scale redox and organic-matter dynamics in soils and sediments by translating that heterogeneity into explicit consequences for ecosystem-scale flux calculations (Wilson and Gerber, 2021), and recent three dimensional reactive transport simulations that resolve plant mediated oxygen inputs and fine scale porewater biogeochemistry in coastal wetlands (Zhou et al., 2024) by translating those insights into an explicit correction for marsh carbon flux calculations.

Finally, by framing the bias in terms of the mean and variance of oxygen and the microbial half-saturation constant, our results fit naturally alongside DAMM-type applications in terrestrial and wetland systems. DAMM has been used to represent soil heterotrophic respiration in land surface and ecosystem models, where it replaces simple Q₁₀ moisture scalars with a mechanistic treatment of temperature, substrate, and diffusion constraints. Our analysis extends this line of work by showing that, once oxygen limitation is represented explicitly, the spatial distribution of oxygen within a model grid cell becomes a first-order control on predicted respiration, rather than a secondary detail (Moyano et al., 2018).

4.3 Implications for modelling and observation

The most direct implication for marsh process models is that replacing heterogeneous rhizosphere oxygen with a single mean value will, under realistic conditions, tend to overestimate aerobic heterotrophic CO₂ production. In our marsh case study, gas-phase oxygen heterogeneity produces a mean aggregation bias of −12 %, that is, a uniform-oxygen configuration overestimates respiration by roughly one tenth. Because net carbon sequestration reflects the balance between inputs and heterotrophic losses, this implies a systematic underestimation of the marsh carbon sink. If a bias of this magnitude were broadly representative of tidal marshes, then, given published global salt-marsh carbon burial estimates of about 10–21 Tg C yr⁻¹ (Alongi, 2020; Ouyang and Lee, 2014), correcting for aggregation bias would increase inferred global marsh sequestration by on the order of 1–2.5 Tg C yr⁻¹, corresponding to roughly 0.1–0.25 Gt C over a century. Given the wide uncertainty in global salt-marsh carbon burial, with estimates based on carbon accumulation rates spanning roughly 0.9–31.4 Tg C yr⁻¹ (Chastain et al., 2022), these numbers should be interpreted as order of magnitude extrapolations rather than precise forecasts. Nevertheless, they indicate that aggregation bias is non negligible for coastal carbon budgets and for greenhouse gas inventories that include blue carbon ecosystems, even though the corrected flux remains small compared with anthropogenic fossil fuel emissions (IPCC, 2023).

At the land surface and Earth system model scale, most soil biogeochemical schemes still employ relatively simple representations of heterotrophic respiration, often combining a temperature response with an empirical soil moisture scalar. More mechanistic schemes, including DAMM and related enzyme kinetic formulations, are now being incorporated into regional and global models to capture interactions among temperature, moisture, and substrate supply more explicitly. Recent work using a reduced complexity carbon–climate model shows that replacing a first order respiration formulation with a Michaelis–Menten formulation can substantially increase the strength of the terrestrial carbon–climate feedback and raise the remaining carbon budget to meet a given temperature target (Beer, 2025). This underlines that large scale projections are structurally sensitive to the mathematical form chosen for respiration. Our results add that, when such schemes include an explicit oxygen limitation term but still operate at coarse spatial resolution, they also inherit an additional structural error that depends on how oxygen is aggregated within each grid cell. Regions with shallow water tables, strong hydrologic variability, and high root biomass, such as tidal marshes, riparian wetlands, peatlands and some rice systems, are likely to be particularly sensitive.

In practical terms, the modelling framework developed here can be used to diagnose where aggregation bias is likely to matter and to design simple mitigation strategies. First, modellers can characterize or bracket the relative saturation of oxygen, that is, the ratio between the mean oxygen in the rhizosphere and a plausible microbial half-saturation constant, together with a measure of subgrid variability. A convenient measure of variability is the coefficient of variation, defined as the standard deviation of oxygen divided by its mean. When relative saturation is very low or very high, or when variability is small, our results suggest that aggregation error will be negligible. When relative saturation is intermediate and variability in oxygen is substantial, a negative aggregation bias should be expected if a single bulk oxygen value is used.

Second, if subgrid oxygen distributions cannot be resolved explicitly, our Jensen-type diagnostics offer a way to correct for at least part of the bias or to represent it as an uncertainty range. For small to moderate variability, a Taylor series approximation using the local curvature of the oxygen-limitation term can be used to derive a multiplicative correction factor for the uniform-oxygen flux, which we demonstrate for the gas-phase proxy (Fig. A7). When variability is large or poorly constrained, it may be more appropriate to propagate a plausible range of oxygen variance into a range of predicted fluxes, thereby making the aggregation error explicit in model uncertainty budgets.

Finally, we note that a conceptually similar issue arises when oxygen limitation is parameterized as a function of another state variable, such as water content or water-table depth. If oxygen is first expressed as a nonlinear function of water content, and then the model uses mean water content within a grid cell, the resulting mean-water baseline can differ from the baseline defined by mean oxygen. In extreme cases, this can even alter the apparent sign of the aggregation effect. Our analysis suggests that it is preferable to define the uniform reference case directly in oxygen space (that is, to evaluate respiration at the mean oxygen concentration) rather than in terms of the mean of a proxy variable such as water content.

4.4 Uncertainties and limits

Several limitations of our approach should be kept in mind when interpreting the results. First, the gas-phase and dissolved-phase oxygen distributions we impose are proxies that reflect hydrologic and structural controls (through water content, porosity, and effective diffusivity) but do not resolve root architecture, advective transport, or biogeochemical feedback between reactions and oxygen supply. They are designed to constrain plausible degrees of rhizosphere heterogeneity rather than to reproduce site specific oxygen probability density functions. Fully coupled transport–reaction models that explicitly resolve root geometry, gas-phase diffusion, and solute flows would be needed to predict detailed oxygen distributions for particular marshes or plant communities.

A related simplification is that aerobic heterotrophic respiration is treated as the only sink for O₂ in the root zone. In marsh soils, the re-oxidation of reduced compounds, including dissolved sulfide, iron sulfides and ammonium, also consumes O₂ and competes with respiration for it (Howes and Teal, 1994; Koop-Jakobsen et al., 2017; Pezeshki and DeLaune, 2012). Because the O₂ distributions are prescribed here rather than obtained from a coupled transport–reaction model, this competition is not resolved, and it affects how representative the imposed distributions are of a given marsh rather than the comparison between the uniform and the heterogeneous configuration.

Second, we treat microbial oxygen affinity through a single half-saturation constant and explore its influence via sensitivity sweeps, rather than by refitting oxygen kinetics for each hydrologic or redox regime. In reality, microbial communities and enzyme systems can acclimate or adapt to chronic hypoxia or frequent redox oscillations, shifting effective half-saturation constants and potentially altering the degree of concavity in the oxygen response. Our sensitivity analysis indicates that the sign of the aggregation bias is robust as long as the oxygen-limitation term retains its saturating form, but the exact magnitude of the bias in any specific marsh will depend on the local microbial community and redox history.

Third, our calibration focuses on aerobic heterotrophic respiration derived from long term soil incubations under controlled temperature and water content. Autotrophic root respiration and anaerobic pathways (for example fermentation, sulfate reduction, and methanogenesis) are outside the scope of the present study. In the field, these processes coexist and interact; root respiration can covary with oxygen and hydrology, while anaerobic pathways contribute substantially to CO₂ and CH₄ production under persistent anoxia. Our results should therefore be interpreted as applying to the aerobic heterotrophic component of belowground CO₂ production, not to total soil respiration or full greenhouse gas budgets.

Lastly, we use a single marsh system and a fixed temperature as our calibration context. While the DAMM formulation is designed to be portable across sites and temperature regimes, additional calibration and validation against independent datasets from other marsh types, climates, and plant communities would be needed to confirm the quantitative magnitude of aggregation bias across the global diversity of coastal wetlands.

4.5 Observational priorities and data needs

Our diagnostics highlight a specific set of observational needs that would allow the aggregation mechanism to be better constrained in marshes and other wetland ecosystems. First, direct measurements of oxygen distributions in marsh rhizospheres are still rare relative to the spatial and temporal variability expected from hydrology and plant traits. Fine scale microelectrode and planar optode measurements (Koop-Jakobsen and Wenzhöfer, 2015) have proven capable of resolving steep oxygen gradients and dynamic oxic–anoxic interfaces in sediments and wetland soils. Extending these techniques across tidal elevations, plant functional types, and seasons would provide the empirical oxygen means, variances, and coefficients of variation needed to apply our bias diagnostics more directly.

Second, joint observations of oxygen, water content, and CO₂ fluxes in the root zone would permit explicit tests of the DAMM based aggregation framework. For example, combining high resolution oxygen imaging with controlled manipulations of water level and plant species could reveal how hydrologic state and radial oxygen loss map into the statistical properties of rhizosphere oxygen distributions, and how those properties in turn control aerobic respiration.

Third, improved constraints on microbial oxygen kinetics in marsh soils through laboratory respiration assays across oxygen gradients, in combination with molecular or enzymatic measurements, would narrow the plausible range of half saturation constants and inform when and where the concavity driven aggregation bias is likely to be large. Together, such observational efforts would reduce structural uncertainty in marsh carbon models, enable more targeted use of correction factors in DAMM type schemes, and provide a foundation for representing rhizosphere oxygen heterogeneity in land surface and Earth system models.

5 Conclusions

We examined how small-scale heterogeneity in rhizosphere oxygen affects modelled aerobic heterotrophic respiration in tidal marsh soils. Using a Dual Arrhenius–Michaelis–Menten formulation calibrated to long term marsh soil incubations, we compared respiration evaluated at a single mean oxygen concentration with respiration averaged over heterogeneous oxygen distributions generated by hydrologic and structural variation in the root zone. For gas-phase oxygen fields representative of tidally influenced marsh soils, evaluating respiration at the mean oxygen led to an overestimate of aerobic CO₂ production on the order of 10 %–15 %.

The magnitude of the bias is controlled by the ratio between mean oxygen and the microbial half-saturation constant for oxygen, together with the variance of oxygen in the rhizosphere. Biases are negligible when oxygen is either extremely scarce or effectively saturating, or when spatial variability is small, and they are largest at intermediate oxygen levels and substantial variability.

These findings have direct implications for marsh carbon budgets and for the representation of wetland soils in biogeochemical and Earth system models. In models that employ mechanistic, oxygen-limited respiration schemes but represent oxygen as a single bulk variable per grid cell, the concavity driven aggregation bias identified here can introduce a structural overestimate of aerobic heterotrophic CO₂ production in systems with strongly heterogeneous rhizospheres. Our diagnostic framework provides a simple way to anticipate when this bias is likely to be important, namely where mean oxygen in the root zone is comparable to microbial half-saturation constants and where variability in oxygen is nontrivial, and to design correction factors or uncertainty ranges accordingly.

Finally, the work points to clear observational and modelling priorities. High resolution measurements of rhizosphere oxygen distributions, coupled with constraints on microbial oxygen kinetics, would allow aggregation effects to be quantified and corrected for specific marsh types. Extending the present framework to include autotrophic root respiration, anaerobic pathways, and fully coupled transport reaction models would further clarify how plant traits, hydrology, and microbial communities jointly shape carbon fluxes in coastal wetlands. By explicitly linking rhizosphere oxygen heterogeneity to aggregation bias in a widely used respiration scheme, this study provides a template for incorporating small scale redox dynamics into ecosystem and Earth system scale assessments of marsh carbon cycling.

Appendix A

Figure A1 Dissolved-phase O₂ proxy distribution with the same units and proxy definition as in Fig. 2b, shown on a logarithmic y axis. The vertical lines indicate the distribution mean (μ) and the unit-matched $K_{\mathrm{M},{\mathrm{O}}_{\mathrm{2}}}$ value used in the dissolved sensitivity sweep in Fig. 4, placing the gas and dissolved cases on a comparable relative O₂ saturation scale.

Figure A2 Sensitivity of aggregation bias to the choice of dissolved-phase O₂ driver. (a) Uniform (mean O₂) versus heterogeneous (full distribution) DAMM respiration curves driven by the unscaled dissolved O₂ proxy. Because dissolved O₂ is orders of magnitude below the calibrated $K_{\mathrm{M},{\mathrm{O}}_{\mathrm{2}}}$, the O₂ limitation is effectively linear and the bias is negligible. (b) Same analysis but with $K_{\mathrm{M},{\mathrm{O}}_{\mathrm{2}}}$ rescaled to match the mean relative saturation of the gas-phase case. Under this unit-matched $K_{\mathrm{M},{\mathrm{O}}_{\mathrm{2}}}$, the dissolved proxy exhibits a negative aggregation bias similar in magnitude to the gas-phase proxy.

Figure A3 Jensen diagnostic for a symmetric two-point O₂ family. Bias in respiration E[f(O₂)] − f(E[O₂]) versus Var[O₂] at fixed mean O₂ and C_org, computed from explicit DAMM evaluations (primary curve) and from the second-order Taylor approximation based on the curvature of O₂ limitation (secondary curve). The close agreement over a broad variance range supports the use of a curvature-based Jensen correction.

Figure A4 Aggregation bias versus Var[O₂] at fixed mean O₂ for two families of O₂ distributions: symmetric two-point mixtures and lognormal distributions. For each family, the mean O₂ is held constant and the variance is varied. In both cases, the bias remains negative and its magnitude increases monotonically with Var[O₂], indicating that the concavity mechanism is robust to the choice of O₂ distribution shape.

Figure A5 Probability density of O₂ obtained by mapping a narrow water-content band, θ ∈ [0.65,0.68], to gas-phase O₂ using the same proxy as in the main analysis; the density is shown on a log scale. This high-θ band has much lower variance than the full rhizosphere distribution. Evaluated at C_org = 0.15 g C cm⁻³, the respiration predicted at the mean O₂ (uniform configuration) is slightly larger than the heterogeneous average E[f(O₂)], giving an aggregation bias of −0.5 % (uniform overpredicts). Numerically, R_uniform = 1.2657 × 10⁻⁴ mg C cm⁻³ h⁻¹ and R_{heterogeneous} = 1.2592 × 10⁻⁴ mg C cm⁻³ h⁻¹ (calibrated DAMM parameters; substrate term evaluated at θ_ref = 0.30).

Figure A6 Response surface of aggregation bias as a function of mean(O₂) $/$ $K_{\mathrm{M},{\mathrm{O}}_{\mathrm{2}}}$ and the coefficient of variation (CV) of O₂, for a symmetric two-point O₂ family at fixed mean O₂. Colour shading indicates bias (%) with negative values only, consistent with concave O₂ limitation. Bias magnitude increases with CV and peaks at intermediate relative saturation, summarising how physiology ($K_{\mathrm{M},{\mathrm{O}}_{\mathrm{2}}}$) and O₂ heterogeneity jointly control aggregation error.

Figure A7 Example of a small variance Jensen correction applied to the gas-phase O₂ proxy. The uniform curve shows respiration evaluated at the mean O₂, the heterogeneous curve shows the explicit average over the prescribed O₂ distribution, and the corrected curve applies a second-order curvature-based factor using the local O₂ variance and curvature. For the modest O₂ variability considered, the correction removes most of the negative aggregation bias, illustrating a lightweight approach for models that cannot explicitly resolve subgrid O₂.