the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
When and why microbialexplicit soil organic carbon models can be unstable
Samia Ghersheen
Salim Belyazid
Stefano Manzoni
Microbialexplicit soil organic carbon (SOC) cycling models are increasingly being recognized for their advantages over linear models in describing SOC dynamics. These models are known to exhibit oscillations, but it is not clear when they yield stable vs. unstable equilibrium points (EPs) – i.e., EPs that exist analytically but are not stable in relation to small perturbations and cannot be reached by transient simulations. The occurrence of such unstable EPs can lead to unexpected model behavior in transient simulations or unrealistic predictions of steadystate soil organic carbon (SOC) stocks. Here, we ask when and why unstable EPs can occur in an archetypal microbialexplicit model (representing SOC, dissolved OC (DOC), microbial biomass, and extracellular enzymes) and some simplified versions of it. Further, if a model formulation allows for physically meaningful but unstable EPs, can we find constraints in the model parameters (i.e., environmental conditions and microbial traits) that ensure stability of the EPs? We use analytical, numerical, and descriptive tools to answer these questions. We found that instability can occur when the resupply of a growth substrate (DOC) is (via a positive feedback loop) dependent on its abundance. We identified a conservative, sufficient condition in terms of model parameters to ensure the stability of EPs. Principally, three distinct strategies can avoid instability: (1) neglecting explicit DOC dynamics, (2) biomassindependent uptake rate, or (3) correlation between parameter values to obey the stability criterion. While the first two approaches simplify some mechanistic processes, the third approach points to the interactive effects of environmental conditions and parameters describing microbial physiology, highlighting the relevance of basic ecological principles for the avoidance of unrealistic (i.e., unstable) simulation outcomes. These insights can help to improve the applicability of microbialexplicit models, aid our understanding of the dynamics of these models, and highlight the relation between mathematical requirements and (in silico) microbial ecology.
 Article
(2820 KB)  Fulltext XML

Supplement
(1558 KB)  BibTeX
 EndNote
Current Earth system models (ESMs) have very simplified representations of soil organic carbon (SOC) dynamics (Bradford et al., 2016; ToddBrown et al., 2013; Varney et al., 2022). Accuracy in matching observed SOC stocks and turnover times has not significantly improved in the latest ensemble of ESMs used in the Coupled Model Intercomparison Project (CMIP, CMIP6) (Varney et al., 2022), with uncertainty about SOC responses to climate change remaining high (ToddBrown et al., 2013; Varney et al., 2022). Consequently, a need to improve and diversify the description of SOC dynamics in ESMs has been identified (Bradford et al., 2016; ToddBrown et al., 2013; Varney et al., 2022; Wieder et al., 2015, 2018). Current ESMs employ linear degradation kinetics to simulate SOC degradation (ToddBrown et al., 2013). Thereby, they miss to integrate important aspects of our current understanding of major controls on SOC fate and to acknowledge the uncertainties in describing these processes (Abramoff et al., 2018; Abs et al., 2023; Bradford et al., 2016; Wieder et al., 2015, 2018). Nonlinear, microbialexplicit SOC models can improve model–data agreement (Hararuk et al., 2015; Wieder et al., 2013). These models vary in terms of the number and identity of C pools and the degree of nonlinearity (e.g., Allison et al., 2010; Manzoni and Porporato, 2009; Schimel and Weintraub, 2003; Wang et al., 2013, 2015; Wieder et al., 2014, 2015). Among these models, the AWB (Allison–Wallenstein–Bradford) model (Allison et al., 2010) has emerged as an archetypal model structure to study the influence of soil microbial processes on carbon stocks (e.g., Abs et al., 2022; Calabrese et al., 2022; Georgiou et al., 2017; Hararuk et al., 2015; Tao et al., 2023; Wieder et al., 2015). The AWB model explicitly represents pools of microbial biomass, extracellular enzymes produced by microbes, polymeric SOC that is not available for microbial uptake, and a pool of available dissolved organic carbon (DOC) produced from enzymatic depolymerization of SOC (Allison et al., 2010, Fig. 1a). With only four C pools and, commonly, two nonlinear terms, the AWB model retains a comparably simple structure and remains somewhat analytically tractable.
While better at predicting modernday SOC stocks, microbialexplicit SOC models are known to exhibit oscillatory behavior (e.g., Georgiou et al., 2017; Manzoni and Porporato, 2007; Sierra and Müller, 2015; Wang et al., 2014, 2016). Such oscillations can represent carbon–microbe dynamics observed at small spatial scales (see, e.g., discussion in Manzoni and Porporato, 2007) but are unfavorable for application at larger spatial and temporal scales, where such oscillations are generally not observed (Georgiou et al., 2017; Wang et al., 2014). A few studies have analyzed the oscillatory properties of some microbialexplicit SOC models (Georgiou et al., 2017; Manzoni and Porporato, 2007; Wang et al., 2014, 2016). These studies characterized the dynamics exhibited after a perturbation around a model's equilibrium (i.e., when all state variables are at steady state): does a model directly converge back to its previous equilibrium or does it approach the equilibrium with dampened oscillations? Different degrees of oscillatory behavior have been described, but, generally, these models were found to be stable (that is, they do converge back to their previous equilibrium) for given parameterizations or if they follow basic principles such as mass conservation and dependence of fluxes on source pools (Sierra and Müller, 2015; Wang et al., 2014, 2016). Stable oscillatory behavior, however, is only one of the possible dynamics such nonlinear models can exhibit. In fact, these models can also be unstable (that is, after perturbation, a model does not converge back to its previous equilibrium) (Abs et al., 2022; Raupach, 2007; Schimel and Weintraub, 2003; Sierra and Müller, 2015), but the occurrence of unstable equilibria in microbialexplicit SOC models remains largely unexplored. While unstable equilibrium points exist analytically, they can never be reached by transient simulations. Thus, model parameterizations that yield unstable equilibria can lead to unpredictable simulation outcomes as amplifying oscillations can occur; expected equilibrium states may not be reached (because they are unstable), hindering convergence in model spinup; or (some) state variables might collapse (e.g., Fig. 1b, yellow line). Further, if C stocks are predicted based on analytical steadystate solutions, unstable equilibria might lead to unrealistic predictions, causing mismatched outcomes from dynamic simulations. To increase the reliability of model predictions and model applicability, it is important to understand when and for what reasons microbialexplicit SOC models become unstable.
Here, we study an archetypal microbialexplicit SOC model (based on the AWB model of Allison et al., 2010, and some simplified versions of it) to answer the following questions: (1) what mechanisms in microbialexplicit SOC models (model structures, used kinetic formulations, and parameter values) cause unstable equilibrium points to emerge? (2) How can we select model structures and/or constrain model parameters to ensure the stability of equilibrium points?
2.1 Archetypal microbialexplicit SOC model
We start by defining the C mass balance equations for a system encompassing SOC (S), dissolved organic C (DOC, D), microbial biomass C (MBC, B), and extracellular enzyme C (ENZ, E) (Eqs. 1–4). The C compartments and flows are illustrated in Fig. 1a, and symbols for the variables and fluxes are defined in Tables 1 and 3. The C mass balance equations are written as a system of ordinary differential equations (ODEs).
Organic matter enters the system with flux I, which is partitioned between SOC and DOC depending on the fraction f_{I}. SOC is depolymerized by extracellular enzymes at rate P and then is transferred to the DOC pool. DOC is directly available for microbial uptake at rate U. Both P and U are nonlinear functions that can take on various forms (Table 2).
Microbes assimilate the substrate with a maximal efficiency y_{m}≤1 that is limited by physiological and/or thermodynamic constraints (Chakrawal et al., 2022) and use the substrate either for growth (i.e., biomass production at rate y_{B}U) or to produce extracellular enzymes (at rate (y_{m}−y_{B})U). We refer to this uptakedependent pathway of extracellular enzyme production as “inducible” ENZ production. An alternative or complementary mode of ENZ production is the biomassdependent “constitutive” ENZ production at rate R_{E}, given by
where m_{E} is the rate constant of constitutive ENZ production. In both formulations, we assumed that respiratory costs associated with enzyme production are already included in the growth respiration (proportional to 1−y_{m}). Two limiting cases can be derived from this general description of extracellular enzyme production:
Both MBC and ENZ are assumed to decay with a linear decay rate D_{i}:
but we also consider densitydependent microbial decay to be an alternative to the linear kinetic (${D}_{B}^{\prime}={d}_{B}^{\prime}{B}^{b}$, with $\mathrm{1}<b\le \mathrm{2}$; Georgiou et al., 2017). All decayed ENZs are assumed to return to the DOC pool, while only a fraction r_{B} of the decayed microbial biomass is recycled in the system and partitioned between SOC and DOC according to the factor f_{D}. In turn, (1−r_{B})D_{B} represents linear microbial maintenance respiration. SOC, DOC, and ENZs can have abiotic losses L_{i} (e.g., erosion, leaching).
The system of Eqs. (1)–(4) constitutes a model of SOC cycling of varying complexity depending on the chosen kinetics. We refer to this fourpool model version as the SDBE model according to the represented state variables. We use this system as a starting point for our analysis but reduce it to simpler variants to derive specific analytical results (Sect. 2.2).
Commonly, depolymerization of SOC by extracellular enzymes is described by either multiplicative (m), forward Michaelis–Menten (f), reverse Michaelis–Menten (r) (Schimel and Weintraub, 2003), or equilibrium chemistry approximation (ECA, e) (Tang and Riley, 2013) kinetics. Table 2 lists the respective formulations of P_{i} ($i=m,f,r,e$), where ${v}_{i}^{p}$ is the maximal depolymerization rate coefficient, and ${K}_{i}^{p}$ is the respective halfsaturation constant (if applicable). The uptake of DOC by microbes can be described with similar formulations, such as U_{j} ($j=m,f$), simply by replacing S with D and E with B (Table 2), where ${v}_{j}^{u}$ is the maximal uptake rate coefficient, and ${K}_{f}^{u}$ is the respective halfsaturation constant.
Many combinations of depolymerization and uptake kinetics are possible. For model versions with both nonlinear terms, we limit our analysis to only a few combinations of depolymerization and uptake kinetics (indicated by the subscript i×j for the ith depolymerization kinetic and jth uptake kinetic), namely m×m, f×f, and r×f (see the summary of analyzed scenarios in Table 4). The first combination employing only multiplicative kinetics facilitates analytical tractability, the second combination is commonly used in other models (e.g., Allison et al., 2010; Georgiou et al., 2017; Tao et al., 2023), and the third combination is based on the conclusions of Tang and Riley (2019) that r×f might be an appropriate (and analytically tractable) approximation of ECA kinetics.
To improve the analytical tractability of the fourpool model, we neglect abiotic losses of SOC and ENZ by setting ${L}_{S}={L}_{E}=\mathrm{0}$ and by limiting the analysis only to the case of constitutive ENZ production (y_{B}=y_{m}, Eq. 6) (Table 4). The Jacobian matrix of partial derivatives for the fourpool model ${J}_{i\times j}^{\mathrm{SDBE}}$ is given in Eq. (A17). The matrix ${J}_{i\times j}^{\mathrm{SDBE}}$ is expressed in terms of general depolymerization and uptake kinetics P_{i} and U_{j}, which allows for stability analysis of the SDBE model irrespective of the specific kinetics.
2.2 Reduced models for mathematical analysis
To identify which and how structural elements of the fourpool model with two nonlinear kinetics affect model stability, we introduce two reduced model versions:

the SBE (SOC–MBC–ENZ) model, neglecting DOC dynamics, and

the SDB (SOC–DOC–MBC) model, assuming ENZ to be at a quasisteady state.
Both model versions have only three pools but are different as, in the SBE model, only one nonlinear term remains, while the SDB model still has both nonlinear depolymerization and uptake kinetics. We analyze the former, less nonlinear model, for all depolymerization kinetics listed in Table 2 with both constitutive and inducible ENZ production pathways and including abiotic losses of S and E (Table 4). In contrast, we analyze the latter, more nonlinear model, after applying the same simplifying assumptions as for the fourpool model (Table 4) – that is, we set ${L}_{S}={L}_{E}=\mathrm{0}$, y_{B}=y_{m}, and we only consider three combinations of depolymerization and uptake kinetics (m×m, f×f, and r×f).
2.2.1 SBE model
DOC dynamics are neglected in the SBE (SOC–MBC–ENZ) model. Instead, it is assumed that any organic carbon that is made available by depolymerization of SOC is directly taken up by microbes – that is, U=P. The flux of decayed extracellular enzymes D_{E} enters the SOC pool, and the partitioning factors f_{I} and f_{D} are set to 1. The resulting system of equations is given for the ith kinetic formulation for P (Table 2) by
Note that, in this formulation, unless ENZ production is purely constitutive (Eq. 6), ENZ production is (partly) independent of microbial biomass (as P_{i} only depends on E and not on B). Consequently, as soon as there are extracellular enzymes that catalyze depolymerization, further enzyme production follows. The Jacobian matrix of the SBE model ${J}_{i}^{\mathrm{SBE}}$ is a 3×3 matrix given by Eq. (A1).
2.2.2 SDB model
In the SDB (SOC–DOC–MBC) model, the extracellular enzyme pool is assumed to be at a quasisteady state; that is, $\frac{\text{d}E}{\text{d}t}=\mathrm{0}$. With L_{E}=0 and y_{B}=y_{m}, we obtain the quasisteadystate concentration of E from Eqs. (4), (5), and (8) as
The SDB model is obtained by substituting E^{qss} for E in Eqs. (1)–(3), which yields the following after assuming L_{S}=0:
In this model version, two nonlinearities remain. Substituting E^{qss} for E in D_{E} and P_{i}, we obtain ${D}_{E}^{\mathrm{qss}}={R}_{E}$ and ${P}_{i}^{\mathrm{qss}}=f(S,B)$, respectively. The Jacobian matrix of the SDB model, ${J}_{i\times j}^{\mathrm{SDB}}$, is a 3×3 matrix given by Eq. (A10).
2.3 Stability analysis
Here, stability behavior refers to how a model responds to a small perturbation around an equilibrium point (illustrated in Fig. 1b). The equilibrium points are obtained by assuming that the system behavior does not change over time – i.e., by setting the ODEs of all state variables C to be equal to zero ($\frac{\text{d}C}{\text{d}t}=\mathrm{0}$). This yields their steady states C^{*}. For nonlinear systems, stability is determined by the signs of the eigenvalues (λ) of the Jacobian matrix J evaluated at an equilibrium point (denoted as $J{}_{*}$; $\text{eig}\left(J{}_{*}\right)=\mathit{\lambda}$) (e.g., Argyris et al., 2015). Briefly, if the real parts of all eigenvalues are negative (Re(λ)<0), the equilibrium point is stable: the system will converge back to this equilibrium point after a perturbation. Instead, if one or more eigenvalues have positive real parts (Re(λ)>0), the equilibrium point is unstable, and the system will not return to the same state. If the eigenvalues additionally have nonzero imaginary parts (Im(λ)≠0), oscillations around the equilibrium point occur (Fig. 1b). Stability analysis is described in more detail in, e.g., Argyris et al. (2015).
2.3.1 Analytical approach
Because the eigenvalues of the Jacobian matrix can be analytically cumbersome, even in the comparably compact threepool models, we use the Routh–Hurwitz criterion (e.g., Argyris et al., 2015; Horn and Johnson, 1994) for $J{}_{*}$. The Routh–Hurwitz criterion states that all Re(λ) values have negative signs if, and only if,

all coefficients a_{i} of the characteristic polynomial det$(J{}_{*}\mathbf{1}\mathit{\lambda})=\mathrm{0}$ (where 1 is the identity matrix; Eqs. 17 and 18) are positive (i.e., a_{i}>0), and

${a}_{\mathrm{1}}{a}_{\mathrm{2}}{a}_{\mathrm{3}}>\mathrm{0}$ (if $J{}_{*}$ is a 3×3 matrix) or ${a}_{\mathrm{1}}{a}_{\mathrm{2}}{a}_{\mathrm{3}}{a}_{\mathrm{3}}^{\mathrm{2}}{a}_{\mathrm{1}}^{\mathrm{2}}{a}_{\mathrm{4}}>\mathrm{0}$ (if $J{}_{*}$ is a 4×4 matrix).
Thus, by applying the Routh–Hurwitz criterion we can analytically evaluate the stability around the equilibrium points of the nonlinear systems given by the three and fourpool models without directly evaluating λ analytically. The characteristic polynomial for a 3×3 matrix is given by
and for a 4×4 matrix, it is given by
In both cases, a_{1} is the negative trace of $J{}_{*}$ (${a}_{\mathrm{1}}=\text{tr}\left(J{}_{*}\right)$).
2.3.2 Numerical simulations
We also compute λ and the steadystate values of the state variables numerically. If not otherwise specified, 100 000 Monte Carlo simulations were produced by randomly drawing parameter values from (log)uniform distributions using a Latin hypercube sampling algorithm (MATLAB R2022b's lhsdesign function; The MathWorks Inc., 2022). All parameters and their respective ranges are listed in Table 3. Partitioning coefficients were sampled from uniform distributions, while rate constants were log_{10}transformed before sampling.
Following Georgiou et al. (2017) and Sierra and Müller (2015), the stability of equilibria was evaluated using the damping coefficient given by
which ranges between −1 and 1. Note that ζ has positive values only if all Re(λ)<0, indicating a stable equilibrium point, and it has negative values if any Re(λ)>0, indicating an unstable equilibrium point. For Im(λ)=0, ζ is either 1 or −1, indicating no oscillations, while $\mathrm{1}<\mathit{\zeta}<\mathrm{1}$ for Im(λ)≠0 indicates that oscillations occur.
Numerical simulations were carried out in MATLAB 2022b (The MathWorks Inc., 2022), and color maps created by (Crameri, 2018a, b) were used for visualization.
2.3.3 Classification of equilibrium points
Our analyses only consider physically meaningful equilibrium points – that is, only equilibrium points for which all state variables are simultaneously positive and real. Within the physically meaningful equilibrium points, we distinguish three categories:

stable, where all physically meaningful equilibrium points are stable (i.e., stable node or focus points; Argyris et al., 2015);

stable and plausible, where all physically meaningful equilibrium points are stable and also give plausible numerical results;

unstable, where all physically meaningful equilibrium points are not stable (i.e., unstable node or focus points; Argyris et al., 2015).
Based on data synthesized by Wang et al. (2013) and on educated guesses, we applied the following conditions for considering results to be plausible: tOC = SOC + DOC + MBC + ENZ ≤500 mgC g^{−1} (=50 %), DOC/tOC <0.01, MBC/tOC <0.05, and ENZ/MBC <0.1, where tOC indicates the total organic carbon content (the sum of all four carbon pools).
2.3.4 Causal loop analysis
In addition to the mathematical analysis of equilibrium points and their stability, we present causal loop diagrams that qualitatively summarize causal links in a system and the feedbacks they create (Haraldsson, 2004). This analysis can help to understand the behavior a system exhibits after a perturbation around an equilibrium point. In a causal loop diagram, causal connections are depicted by arrows, tying a cause (at the tail of the arrow) to its direct effect (at the head of the arrow) (Haraldsson, 2004). The sign of the causal relation (+ or −) depends on whether an isolated change in one element causes another to change in the same (+) or opposite (−) direction in relation to the initial change (relative to the unchanged state) (Haraldsson, 2004; Richardson, 1986). For example, a decrease in the microbial uptake rate would lead to relatively less microbial biomass (compared to the case where the uptake rate had not changed), describing a positive causal relation. Closed loops with zero or an even number of negative interactions are denoted as positive or reinforcing feedback loops R, and closed loops with a odd number of negative interactions are denoted as negative or balancing feedback loops B (Haraldsson, 2004; Richardson, 1986).
We analyzed three model versions with different model structures (number of state variables and/or nonlinearities; Table 4, Sect. 2.2). We first present analytical results for the simpler threepool models with one (SBE) and two nonlinear terms (SDB), followed by analytical and numerical results for the fourpool SDBE model (Table 4). We use causal loop diagrams to qualitatively interpret these results.
3.1 SBE model: neglecting DOC dynamics
3.1.1 Steadystate solutions
For all kinetic descriptions of the depolymerization rate (Table 2), the threepool SBE model has two equilibrium points (EPs). Of these, one is an “abiotic” equilibrium Q_{0}, where only SOC exists and where microbial biomass and extracellular enzymes are zero; i.e., ${Q}_{\mathrm{0}}=({S}_{\mathrm{0}}^{*},\mathrm{0},\mathrm{0})$. Here, the asterisk indicates a state variable at steady state, and the subscript 0 signifies the abiotic solution. In turn, for each kinetic, there exists an alternative “biotic” equilibrium point with nonzero microbial biomass and extracellular enzymes; i.e., ${Q}_{\mathrm{1},i}=({S}_{\mathrm{1},i}^{*},{B}_{\mathrm{1},i}^{*},{E}_{\mathrm{1},i}^{*})$. The steadystate solutions for these equilibria depend on the ith formulation used to describe P_{i}. All solutions are reported in Table 5, where, for convenience, parameters have been combined into parameter groups: extracellular enzyme turnover α>0, microbial biomass turnover β>0, and parameter groups η≥0 and ω≥0 (Table 3).
While the abiotic equilibrium point is always positive, the biotic one can only be positive (and thus physically meaningful) if
i.e., if the linear SOC loss rate is smaller than SOC inputs. Note that, for f and e kinetics, additional conditions apply for positivity.
If the abiotic loss of SOC is neglected (i.e., l_{S}=0), the abiotic equilibrium point does not exist (SOC would accumulate at the constant rate I), and, in accordance with Eq. (20), the biotic equilibrium is always physically meaningful.
3.1.2 Stability analysis
To analyze whether a physically meaningful equilibrium point is also stable, we apply the Routh–Hurwitz criterion to the Jacobian matrix ${J}_{i}^{\mathrm{SBE}}$ (Eq. A1 in the Appendix) evaluated at the kth equilibrium point (${J}_{i}^{\mathrm{SBE}}{}_{*,k}$) – either the abiotic (k=0) or the biotic equilibrium (k=1) (detailed in Sect. A1 in the Appendix and in Sect. S2.1 in the Supplement). From this, we find that the stability of a physically meaningful abiotic equilibrium is conditional on
whereas all physically meaningful bioticequilibrium points of the threepool SBE model are also stable.
Figure 2a shows a simplified causal loop diagram of the SBE model (sparing all loss and decay terms) that can help to understand the dynamic behavior of the model after a perturbation around an equilibrium point. The reinforcing loop R_{1} describes the increase in microbial biomass with an increasing depolymerization rate (∝ uptake rate), leading to an increased ENZ production rate; an increased ENZ concentration; and, consequently, a further increasing depolymerization rate. This reinforcing effect is dampened by the balancing loops B_{1} (the depletion of SOC by depolymerization) and B_{2} (the carbon cost of ENZ production). The reinforcing loop R_{2} exists only if inducible ENZ production is considered (y_{B}<y_{m}; higher depolymerization stipulates the production of more extracellular enzymes, which promote depolymerization). R_{1} and R_{2} are not independent of each other and have to obey mass balance – i.e., per unit of uptake, an increase in inducible ENZ production ($\propto {y}_{m}{y}_{B}$) can only be achieved by reducing the buildup of microbial biomass (lowering y_{B}). An extreme case of this is y_{B}=0; in this case, no microbial biomass is produced, and only B_{1} and R_{2} remain. From the stability analysis (Sect. A1 and Sect. S2.1 in the Supplement), we obtained no conditionality on the stability of physically meaningful equilibrium points in the SBE model. Thus, for any proportion of constitutive vs. inducible ENZ production, all physically meaningful biotic equilibria are also stable. That is, the dynamic behavior of the model after a perturbation around its equilibrium point is dominated by the balancing feedbacks, ensuring a convergence back to the equilibrium point.
3.1.3 Exclusive stability of either abiotic or biotic equilibrium
We recall that, for the biotic equilibrium to be physically meaningful, it is required that ${S}_{\mathrm{1},i}^{*}<{S}_{\mathrm{0}}^{*}$ (Eq. 20); whereas, for the abiotic equilibrium point to be stable, it is required by Eq. (21) that $\frac{\partial {P}_{i}}{\partial E}{}_{*,\mathrm{0}}<\frac{\mathit{\alpha}\mathit{\beta}}{\mathit{\eta}}$. This condition translates for, e.g., multiplicative kinetics to
This means that, when the biotic equilibrium is physically meaningful, the abiotic equilibrium is unstable and vice versa. Therefore, no region in the parameter space yields a physically meaningful bistability in which biotic and abiotic equilibria are simultaneously physically meaningful and stable. This holds for all evaluated kinetics (see Sect. S2.1 in the Supplement for the remaining analytical derivations).
3.2 SDB model: neglecting ENZ dynamics
3.2.1 Steadystate solutions
Table 6 reports the steadystate solutions of the threepool SDB model, where, for convenience, parameters were further grouped into
and
Because an abiotic equilibrium point exists only for l_{S}>0 (Table 6), we only evaluate the stability of the biotic equilibrium and drop “1” from the subscript for conciseness. Biotic steady states can only be physically meaningful for $I>{l}_{D}{D}_{i\times j}^{*}$, meaning that the DOC leaching flux cannot be larger than the total OC input flux. With f×f and r×f kinetics, it is additionally required that $\mathit{\beta}<{y}_{m}\phantom{\rule{0.125em}{0ex}}{v}_{f}^{u}$, implying that the maximal perbiomass assimilation rate ${y}_{m}\phantom{\rule{0.125em}{0ex}}{v}_{f}^{u}$ must be larger than the microbial biomass turnover ($\mathit{\beta}={m}_{E}+{d}_{B}$; Table 3). For f×f kinetics, it is additionally required that $\frac{{d}_{E}}{{y}_{m}{v}_{f}^{p}{m}_{E}}{\mathit{\gamma}}_{f\times f}<\mathrm{1}$ for steady states to be positive.
In the absence of DOC leaching (for l_{D}=0), with m×m and f×f kinetics, ${S}_{i\times j}^{*}$ becomes independent of I (by Eq. 23), while ${B}_{i\times j}^{*}$ and ${E}_{i\times j}^{qss}$ are linear functions of I. In contrast, for l_{D}=0, ${S}_{r\times f}^{*}$ is linearly dependent on OC input I. For l_{D}>0, ${S}_{m\times m}^{*}$ is a function of $\frac{I}{I{l}_{D}{D}_{m\times m}^{*}}$ (Eq. 23), causing ${S}_{m\times m}^{*}$ to decline with increasing inputs as $\frac{I}{I{l}_{D}{D}_{m\times m}^{*}}\to \mathrm{1}$ for $I\mathit{>>}{l}_{D}{D}_{m\times m}^{*}$. Only in ${S}_{r\times f}^{*}$ does I still appear in a linear term for l_{D}>0 as well.
3.2.2 Stability analysis
The Jacobian matrix of the SDB model around the biotic equilibrium, ${J}_{i\times j}^{\mathrm{SDB}}{}_{*}$, is given by Eq. (A10). Evaluating the respective coefficients of the characteristic polynomial and the requirement that ${a}_{\mathrm{1}}{a}_{\mathrm{2}}{a}_{\mathrm{3}}>\mathrm{0}$ yields a cumbersome sufficient and necessary condition for the stability of the biotic EPs of the SDB model (Eqs. A11–A13; see details in Sect. A2 and Sect. S2.2 in the Supplement). The appearance of such a conditional statement means that, in contrast to the SBE model, the SDB model can have physically meaningful but unstable EPs (i.e., if the conditions described by Eqs. A11–A13 do not hold). A perturbation around such an unstable EP will cause the system to diverge from the EP. In this case, the biotic pools (MBC and quasisteadystate ENZ) will collapse, while DOC will reach a steady state as ${D}_{\mathrm{0},i\times j}^{*}$ (for l_{D}>0), and SOC will accumulate indefinitely.
The obtained sufficient and necessary condition given by Eqs. (A11)–(A13) does not allow for easy interpretation or application. We thus propose a sufficient (i.e., a more conservative or strict) condition for stability that is easier to trace analytically as follows:
This sufficient condition for the stability of the bioticequilibrium points of the SDB model holds for all cases relevant to our analysis (Sect. S2.2 in the Supplement). Described in words, this condition requires the depolymerization rate to be less sensitive to a change in microbial biomass than to a proportional change in SOC. Figure 2b illustrates this relation in a simplified causal loop diagram. The reinforcing loop R_{1} causes the depolymerization rate (${P}_{i}^{\mathrm{qss}}$) to increase as the (quasisteadystate) ENZ concentration increases (quantified by $\frac{\partial {P}_{i}^{\mathrm{qss}}}{\partial B}{}_{*}$). This then increases the DOC concentration and uptake and ultimately causes a further increase in microbial biomass and (quasisteadystate) ENZs. This positive feedback is accelerated by an additional reinforcing feedback loop (R_{2}) that causes uptake to further increase as microbial biomass increases. The balancing loops B_{2} and B_{3}, respectively, describe the depletion of DOC with increasing uptake and the reduction of biomass as more extracellular enzymes are produced. Lastly, the balancing loop B_{1} causes SOC to change in the opposite direction compared to the depolymerization rate (i.e., SOC is depleted as depolymerization increases, and relatively more SOC remains as depolymerization decreases), counteracting the initial change in ${P}_{i}^{\mathrm{qss}}$ (quantified by $\frac{\partial {P}_{i}^{\mathrm{qss}}}{\partial S}{}_{*}$). Therefore, the sufficient stability condition in Eq. (25) can be interpreted in the sense that the negative feedback B_{1} must be quantitatively stronger than the positive feedback R_{1} (by some factor y_{m} and buffered by a constant term; Eq. 25).
In essence, the positive feedback R_{1} can drive the system to overshoot or collapse: e.g., if the microbial biomass or (quasisteadystate) ENZ concentration happens to decrease due to a perturbation, this will reduce the depolymerization rate and, following the positive feedback, will result in further reduced MBC and (quasisteadystate) ENZ concentrations. The biotic pools would collapse, and the system would not be able to recover to its initial equilibrium (i.e., be unstable). Only if the resulting accumulation of SOC increases the depolymerization rate more than it is reduced by the depletion of (quasisteadystate) ENZs will the system be able to recover and retain the biotic components (i.e., be stable).
We note that, for linear uptake kinetics (i.e., ${U}_{\mathrm{l}}={v}_{l}^{u}D$), additional terms in the necessary and sufficient conditions can help to fulfill the Routh–Hurwitz criterion for stability (see details in Sect. A2). In the causal loop diagram, linear uptake kinetics remove the positive feedback between uptake rate and biomass (R_{2} in Fig. 2b vanishes). Although we could not show this analytically, numerical evaluation showed that, for the chosen parameter spaces (Table 3) with linear uptake kinetics U_{l}, the physically meaningful EPs of the SDB model were always stable (Fig. S3c–d in the Supplement).
3.3 SDBE model: full archetypal model
The fourpool SDBE model with ${L}_{S}={L}_{E}=\mathrm{0}$ and y_{B}=y_{m} has the same steadystate solutions as the threepool SDB model, but now the solution for the ENZ pool is denoted as ${E}_{i\times j}^{*}$ because ENZ is not considered to be at a quasisteady state as in the SDB model (Table 6).
3.3.1 Analytical stability analysis
In the fourpool SDBE model, the coefficients of the characteristic polynomial of ${J}_{i\times j}^{\mathrm{SDBE}}{}_{*}$ (Eq. A17) remain analytically tractable (Sect. S2.3 in the Supplement). The trace of ${J}_{i\times j}^{\mathrm{SDBE}}{}_{*}$ is always negative (and, thus, $\text{tr}\left({J}_{i\times j}^{\mathrm{SDBE}}{}_{*}\right)={a}_{\mathrm{1}}>\mathrm{0}$), and its determinant is always positive (det$\left({J}_{i\times j}^{\mathrm{SDBE}}{}_{*}\right)={a}_{\mathrm{4}}>\mathrm{0}$). However, the additional Routh–Hurwitz criterion for the 4×4 matrix ${J}_{i\times j}^{\mathrm{SDBE}}{}_{*}$ (given by ${a}_{\mathrm{1}}{a}_{\mathrm{2}}{a}_{\mathrm{3}}{a}_{\mathrm{3}}^{\mathrm{2}}{a}_{\mathrm{1}}^{\mathrm{2}}{a}_{\mathrm{4}}>\mathrm{0}$) becomes analytically intractable. Despite this additional complexity, we can still draw some conclusions based on similarities between the SDB and SDBE models. Considering that, in the SDBE model, ENZ dynamics are explicitly represented (and, thus, e.g., $\frac{\partial {P}_{i}^{\mathrm{qss}}}{\partial B}\to \frac{{m}_{E}}{{d}_{E}}\frac{\partial {P}_{i}}{\partial E}$), similar conditions emerge for the positivity of the coefficients of the characteristic polynomial, as in the SDB model (Sect. S2.3 in the Supplement). Based on these similarities, we propose that the sufficient condition for stability found for the threepool SDB model might also hold in the SDBE model. This proposed sufficient condition is given by
The simplified causal loop diagram of the SDBE model in Fig. 2b gives rise to the same interpretation of this condition as in the SDB model. In the following, we confirm that this condition holds in the SDBE model via numerical analysis.
3.3.2 Numerical stability analysis
Testing the sufficient condition for stability
We produced 100 000 Monte Carlo simulations and computed the damping coefficient ζ (Eq. 19) to numerically evaluate the stability of equilibrium points in the SDBE model. Within the sampled parameter space (Table 3), physically meaningful equilibrium points for which ${Z}_{i\times j}\ge \mathrm{0}$ (Eq. 26) also always had ζ>0 and were thus stable (Fig. 3a–b, Fig. S1 in the Supplement). Damping coefficients with ζ≤0 were only observed when ${Z}_{i\times j}<\mathrm{0}$ (that is for points above the black line in Fig. 3a or red points in Fig. 3c). While a total of 46 932 evaluated equilibrium points were physically meaningful and stable, fewer than half of these (22 386) also fulfilled the condition ${Z}_{i\times j}\ge \mathrm{0}$. In turn, the majority of these equilibrium points (24 546) were stable despite contradicting this condition – i.e., the condition given by ${Z}_{i\times j}\ge \mathrm{0}$ is very conservative. In case the condition is fulfilled, the value of Z_{i×j} correlates well with the value of ζ, meaning that, for larger Z_{i×j}, oscillations are generally more dampened (Fig. 3b).
Changing environmental conditions, microbial physiology, and stability
Even for the simplified sufficient condition ${Z}_{i\times j}\ge \mathrm{0}$, analysis is cumbersome for kinetics other than m×m (Sect. S2.3 in the Supplement). We thus varied specific parameters individually and evaluated their effect on the numerically computed damping coefficient ζ.
Keeping all microbial and enzymatic parameters constant (set to their baseline value – that is for constitutive ENZ production; Table 3), stability depends on the environmental control parameters l_{D} and I (Fig. 4a). At baseline parameter values, most environmental conditions yield stable EPs, but strong oscillations occur around these EPs (damping coefficient ζ<1). As conditions become less favorable and as either l_{D} increases and/ or I decreases, equilibrium points can become unstable (ζ<0).
Next, we analyzed the influence of individual microbial parameter values on the stability of EPs for a number of scenarios defined by combinations of I and l_{D} (Fig. 4b–d; Table 7). Generally, if DOC leaching is neglected (solid lines, Fig. 4b–d), the variation in just one parameter rarely leads to unstable equilibria (only at very high microbial decay rates d_{B}). In contrast, if DOC leaching occurs, variation in key physiological parameters can lead to a transition from stable to unstable EPs. This happens, e.g., as y_{m} becomes too low or d_{B} becomes too high – i.e., for specific environmental conditions, there are lower threshold values for y_{m} and upper thresholds values for d_{B} beyond which equilibria become unstable (Fig. 4b–c, Table 7). With the exception of f_{D}, the partitioning of decayed microbial biomass between SOC and DOC, all parameters show such a threshold within the explored ranges (Fig. S2 and Sect. S2.3 in the Supplement).
Using the alternative description of inducible production of extracellular enzymes, the stability behavior with respect to changes in y_{B} is more varied (Fig. 4d). As for y_{m}, there are lower y_{B} threshold values below which steady states become unstable. However, there can also be upper thresholds for y_{B} above which too few enzymes are being produced to ensure sufficient C acquisition.
In summary, by varying only individual parameters, instabilities can arise when assimilation, depolymerization, or ENZ production are too low or when abiotic C losses are too high (Fig. 4b–d, Fig. S2 in the Supplement, Table 7). These results are in line with the analytic analysis of the sufficient condition (Eq. 26) for m×m kinetics (Sect. S2.3 in the Supplement).
Densitydependent mortality
Alternatively to the conventional linear decay term, Georgiou et al. (2017) proposed a densitydependent formulation of the microbial decay rate (${D}_{B}^{\prime}={d}_{B}^{\prime}{B}^{b}$, with $\mathrm{1}<b\le \mathrm{2}$) that yields mostly stable nonoscillatory behavior (note that this formulation causes both microbial mortality and maintenance respiration to be densitydependent). For our SDBE model with DOC leaching, we could only find an analytical steadystate solution for m×m kinetics and b=2. In this case, the densitydependent formulation could vastly alleviate the previously observed instability and could result in damping coefficients for plausible equilibrium points close to 1 for most of the explored parameter spaces (Fig. S4a in the Supplement). However, some physically meaningful but unstable EPs were still observed. Only with negligible DOC leaching did physically meaningful but unstable equilibrium points vanish completely. This was numerically tested for m×m, f×f, and r×f kinetics (Fig. S4b–d in the Supplement).
Instability and predicted organic carbon pools
Figure 5a–b illustrate the joint distributions of physically meaningful SOC and MBC pools in the f×f model for scenarios where DOC leaching is either considered (Fig. 5a) or neglected (Fig. 5b) (for 10 000 Monte Carlo simulations within the parameter ranges given in Table 3). By accounting for DOC leaching, just about half of the simulations yield physically meaningful EPs – of which most (4698) were also stable, but more than 10 % (687) were unstable. Neglecting DOC leaching increases the total number of physically meaningful results (6896) and simultaneously reduces the relative share of physically meaningful but unstable results to <5 % (312). In both scenarios, most simulations yield implausible results for steadystate C stocks – e.g., MBC being larger than SOC. However, in both scenarios, unstable EPs largely overlap with stable and plausible outcomes in the SOC–MBC solution space. This is also evident in the empirical probability density functions of all four state variables (Fig. 5c–d). Especially if DOC leaching is considered (l_{D}>0, Fig. 5c), values of plausible and unstable steadystate SOC pools largely overlap. In contrast, the distributions of plausible and unstable steadystate pool sizes of MBC, DOC, and ENZ do not overlap as closely as those for SOC. These distinctions are amplified in the cases where DOC leaching is neglected (l_{D}=0, Fig. 5d), in particular for DOC.
4.1 Model structure matters: standard microbialexplicit SOC models can have unstable equilibria
Microbialimplicit models with linear decomposition kinetics are stable as long as there is no inert pool, mass conservation applies, and rates are proportional to the amount of carbon in the donor pool (Sierra and Müller, 2015). It is more difficult to define general stability criteria for nonlinear models, whose structure and type of nonlinearity affects the model behavior around the equilibrium state. Manzoni and Porporato (2007) and Raupach (2007) showed analytically that the nontrivial steady states of twopool models, consisting of a substrate pool and a microbial pool (denoted as a “harvester” system in Raupach, 2007), are always stable for multiplicative and forward Michaelis–Menten kinetics (but only under the assumption that the input to the substrate (I) is a constant; Raupach, 2007). We show here that the same is also true for all physically meaningful, nontrivial (biotic) equilibrium points if a third pool representing extracellular enzymes is added (SBE model). This result holds irrespective of the kinetic laws used to describe SOC depolymerization and whether ENZ production is considered to be constitutive, inducible, or a combination of both.
Interestingly, by introducing a second nonlinear term, Raupach (2007) found that unstable equilibria could emerge in their twopool model. In contrast, Wang et al. (2014, 2016) demonstrated for several versions of a threepool (litter–SOC–microbes) model with two nonlinearities (microbial degradation and subsequent uptake of litter and SOC) that the equilibrium points of these models were always stable. An underlying assumption in these models was that the available substrate pool (similarly to what we described as DOC) was at a quasisteady state. Our derivation of the SBE model follows a similar simplification – and also does not yield unstable behavior. By contrast, unstable equilibrium points are possible in our three and fourpool model versions with two nonlinearities that explicitly consider DOC (SDB and SDBE models; yellow boxes in Fig. 2). Whether equilibrium points in microbialexplicit SOC models can become unstable is thus not dependent on the number of pools or the number of nonlinearities per se but is rather dependent on the combination of nonlinearities, the coupling of different pools and rates, and what feedbacks they create.
Comparing our three and fourpool models to the simpler twopool model analyzed by Raupach (2007) can help to understand why instability can occur in these models. Briefly, their model describes human consumption of a food resource but is similar in structure to our models (analogous terms in our models are given in brackets): the resource (SOC in SBE and DOC in SDB and SDBE models) is taken up by the human consumer (microbes) and thereby depleted. The uptake process is always described as a nonlinear term, equivalent to our description of U_{j}. Raupach (2007) analyzed two different cases with respect to the resupply of the resource (I in SBE and $(\mathrm{1}{f}_{I})I+{P}_{i}^{\u2020}$ in SDB and SDBE): (1) resupply is independent of the available resource, and (2) resupply is dependent on the resource itself. The first case is similar to our SBE model, where the resource SOC is replenished by the external input I and is thus independent of the SOC availability itself. In these cases, the biotic (resource–human coexistence in Raupach, 2007) equilibrium is always stable if it is physically meaningful. In turn, the second case can be compared to our SDB and SDBE models: unless the external input to DOC is very high (for low f_{I}), the replenishment of the resource DOC is dominated by the depolymerization rate ${P}_{i}^{\u2020}$ – which, via the positive feedback loop R_{1}, is dependent on the available DOC (compare Fig. 2b). These models can have unstable biotic equilibria (Raupach, 2007). Therefore, the (more or less direct) dependency of the resource resupply on the abundance of the resource itself can be identified as the root cause that allows for instability in these models.
We further tested this hypothesis by setting f_{I}=0 in the SDBE model (i.e., all external input goes into DOC directly), bypassing the dependency of DOC replenishment on depolymerization. In line with our expectation (and the analytical solution of Z_{m×m} and Z_{f×f} for f_{I}=0 in Sect. S2.3 in the Supplement), this effectively prevented the occurrence of unstable EPs (Fig. S5 in the Supplement).
4.2 Avoiding instability
Our analysis of different model structures and their stability behavior points to three approaches to avoid unstable EPs in microbialexplicit SOC models:

Model structure. Avoid positive feedback coupling between the microbial growth substrate (here DOC) available for uptake and its resupply.

Kinetic formulations. Avoid accelerated depletion of DOC by reducing the dependency of uptake on microbial biomass.

Parameter values. Choose parameter values so that the sufficient and/or necessary conditions for stability are met.
The first approach is commonly taken in models that assume DOC to be in a quasisteady state (e.g., Wang et al., 2014, 2016) but might have shortcomings in cases where DOC dynamics become important – e.g., if drying–rewetting dynamics or leaching are relevant (also see Wutzler and Reichstein, 2013 for a discussion on timescale dependence of appropriate abstraction level). If DOC leaching is not considered to be a relevant process, neglecting this process but keeping a dynamic description of DOC can already considerably reduce the likelihood of unstable EPs.
The second approach is used, e.g., in models that assume DOC uptake to be independent of microbial biomass but dependent on the availability of DOC (e.g., Manzoni et al., 2014). Alternatively, using, for instance, reverse Michaelis–Menten kinetics to describe microbial uptake can dampen oscillations (Wang et al., 2016). Since uptake kinetics using reverse Michaelis–Menten or the ECA formulation become similar to linear uptake kinetics at relatively high concentrations of microbial biomass, they could also help to alleviate instability issues under these conditions. However, when using forward Michaelis–Menten kinetics for decomposition and DOC assimilation, the halfsaturation constants can attain large values after data assimilation, suggesting that decomposition kinetics might be approximated by multiplicative rather than linear kinetics, at least in largescale model applications (Tao et al., 2024b).
Lastly, the third approach might seem straightforward as we could expect parameter values calibrated with measurement data to yield both stable and plausible EPs. However, our numerical simulations indicated that, especially if DOC leaching is considered, calibrating parameter values with SOC and microbial biomass data alone could still lead to plausible yet unstable EPs (Fig. 5). While data on carbon contents in the extracellular enzyme pool are still not available, combining microbial biomass data with quantitative data on DOC pools (as in, e.g., Wang et al., 2013) could help to avoid calibration according to parameter values that lead to unstable EPs. Moreover, stability criteria can be obeyed in various other ways, e.g., by considering correlations between parameter values or introducing additional constraints on microbial physiology.
4.2.1 Correlations between parameter values
Correlations between parameter values could effectively alleviate the occurrence of unstable EPs by simultaneously changing parameter values that appear on both sides of the inequality given by Eqs. (25) or (26), thereby ensuring that these conditions are always fulfilled even as parameter values change. Some evidence proving this to be effective is provided by, e.g., Hararuk et al. (2015), who used the fourpool AWB model (Allison et al., 2010) (similar to the SDBE model) for predictions of global carbon stocks. They prescribed, e.g., the uptake and depolymerization rate coefficients and the respective halfsaturation constants to positively correlate with temperature. Thus, with the same directional change in temperature, these parameter values change in opposite directions with respect to their threshold values for stability (Table 7) – i.e., with a decrease in temperature, ${v}_{f}^{p}$ decreases, moving closer to its threshold; however, simultaneously, ${K}_{f}^{p}$ also decreases, moving further away from its threshold. Consequently, for a wide range of parameter values, the conditions for stability could be fulfilled. Indeed, Hararuk et al. (2015) reported that they did not observe any unstable equilibria with maximumlikelihood parameters in their global study.
Beyond the qualitative assessment of parameter thresholds (Table 7), explicit analytical expressions of the necessary or sufficient conditions for stability (as for Z_{m×m} in Sect. S2.3 in the Supplement) could be used to quantitatively assess what parameter correlations are required to ensure stability of equilibria across reasonable parameter ranges. However, for other kinetic formulations besides m×m, these terms might become difficult to trace analytically.
4.2.2 Ecoevolutionary constraints on microbial traits
The observation that stability thresholds for parameters shift as environmental conditions change can further be interpreted in the light of expected variations in microbial functional traits. While some kinetic rate parameters might be correlated due to thermodynamics (e.g., temperature response of rate parameters), correlations among other parameters, like the investment into growth or extracellular enzyme production, might rather emerge as outcomes of ecoevolutionary processes that select specific combinations of traits in a given environment (Abs et al., 2023, 2024). These combinations of traits would manifest themselves as microbial life history strategies under different environmental conditions (e.g., Malik et al., 2020). Following this logic, changing environmental conditions could constrain the space for microbial physiological adaptation because microbial traits would need to ensure stability. For example, very inefficient microbes (having a low y_{m}) could not establish a stable equilibrium under very unfavorable conditions (low OC input and/or high DOC leaching, Fig. 4b) unless other traits change simultaneously. This reflects a basic principle of ecology: organisms have to be adapted to the environment they inhabit, and such adaptation is expressed through sets of coordinated traits. Currently, this basic principle has not been integrated into microbialexplicit SOC models (but see Abs et al., 2022, for a recent attempt to address this challenge), which can lead to matching specific environmental conditions with a (modeled) microbial population that is not able to sustain itself under those conditions.
While our results are restricted to equilibrium conditions in the absence of external perturbations (e.g., fluctuating environmental conditions), we can speculate that trait coordination should also emerge in a fluctuating environment (for an example of plant trait coordination under stochastic soil moisture, see Bassiouni et al., 2023). In fact, changes in the microbial community with relevance for SOC cycling can happen on different timescales, from hours (stress response) and weeks (changes in community composition) to years and decades (mutations) (Abs et al., 2024). Environmental fluctuations occur at all these scales, and microbial communities are well equipped to adapt to these stochastically changing environments.
Integrating soil microbial ecological understanding into microbialexplicit SOC models could thus yield alternative mathematical descriptions or parameter relations that could prevent mismatching between parameter values and environmental conditions and ultimately improve model applicability (Georgiou et al., 2017). Evidence regarding the importance of microbial ecology and evolution for SOC cycling is accumulating (Abs et al., 2022, 2023). For instance, microbes have been found to invest more into the production of extracellular enzymes in soils with lower SOC contents (Calabrese et al., 2022; Malik et al., 2019), and densitydependent microbial mortality, a concept derived from ecological considerations, can effectively alleviate oscillatory behavior (Georgiou et al., 2017).
Yet, revisiting the proposed stability criterion from an ecoevolutionary perspective prompts a paradox. In a model that obeys the proposed stability criterion (Eqs. 25 and 26), a reduction in enzyme production will always lead to an increase in the depolymerization rate – a strong disincentive for microbial investment into enzyme production. We can resolve this paradox by realizing that, for the SDB and SDBE models analyzed (and applied by, e.g., Hararuk et al., 2015, but also including the original model formulation by Allison et al., 2010), there is no process (such as competition with other (micro)organisms or abiotic removal of (accessible) SOC through, e.g., leaching or occlusion) that competes with microbes for SOC. Consequently, the only alternative to SOC degradation is its accumulation. In turn, this allows for high depolymerization rates with minimal microbial investment into enzyme production. This outcome also emerges from a mathematical analysis of optimal substrate utilization – without losses of substrate due to abiotic processes or competition, decomposers do not have any reason to invest in resource acquisition (Manzoni et al., 2023). Thus, a rigorous ecoevolutionary optimization approach (that, for instance, aims to maximize microbial growth rate) cannot be readily implemented with the current model structure but would require an extension of the model to account for competition for SOC. Our analysis instead demonstrates how (given the model structure) microbes might adapt to environmental conditions not to maximize their fitness but to attain a stable population. Whether stability and evolutionary fitness maximization are convergent (i.e., an organism at an evolutionary fitness maximum would also establish a stable population) poses an interesting question for future research.
4.2.3 Comparing approaches to avoid instability
It is important to note that, while our analysis shows all three of the approaches listed in Sect. 4.2 to be feasible means to avoid instability, it gives no indication of which of the approaches should be preferred. As different objectives and research questions have differing requirements (e.g., whether the occurrence of instabilities should globally be avoided by removing the positive DOC feedback or whether it suffices to constrain the available parameter space), we cannot give a general recommendation. However, an important distinction between the different approaches is that, while the first two (removal of an explicit DOC representation and the dependency of uptake on microbial biomass) are simplifications of the system that might require further justification for specific use cases (Wutzler and Reichstein, 2013), the third approach can add realism to the model by explicitly considering the interaction between environmental conditions and the microbial community. In line with Abs et al. (2022, 2023, 2024) and Georgiou et al. (2017), we want to highlight the potential of using such ecologically consistent mathematical descriptions to improve current model formulations. In other words, we cannot simply add a biotic component to models without acknowledging that this component has to be “adapted” (as species and communities are in the real world) to the environmental conditions it is exposed to.
4.3 Implications
4.3.1 Mathematical insights into microbially mediated SOC cycling
Currently, the debate on the implications of different model formulations for microbially mediated SOC cycling is ongoing (He et al., 2024; Lennon et al., 2024; Tao et al., 2023, 2024a). Our contribution provides an additional perspective to this discussion by leveraging mathematical requirements for stability. Based on our analysis, we cannot draw any conclusions regarding which (if any) model version is most suitable to realistically represent microbially mediated SOC cycling at the large scale. Microbialexplicit models including a DOC pool (SDB and SDBE models) create a positive feedback loop that allows for instabilities to occur, but this finding by no means indicates a shortcoming of such models. In fact, the positive feedback loop is a direct consequence of our conceptual understanding of the microbially mediated, enzymatically driven degradation of organic carbon substrates in soils (Kuzyakov et al., 2000; Kuzyakov, 2010). At small spatial scales, oscillatory behavior can be realistic (Manzoni and Porporato, 2007); additionally, the collapse of a local microbial population is not difficult to imagine (the alternative – and, perhaps, more realistic – outcome that microbes turn dormant cannot be described by these models). However, at large spatial scales, these behaviors are not observed, indicating that SOC cycling is more stable at these scales (Wang et al., 2014; Georgiou et al., 2017). This can be seen as being indicative of a scaling problem – the same process cannot be described by the same mathematical formulations across scales (Chakrawal et al., 2020; Wilson and Gerber, 2021). The proposed stability criterion, as well as the outlined approaches to avoid instability, could, in turn, guide the development of upscaled model formulations – suitable upscaled kinetics might employ one of the proposed approaches to avoid instability or could be designed to obey the proposed (note: very conservative; see further discussion below) stability criteria.
4.3.2 Applicability of the proposed stability criterion
We could identify a sufficient and necessary condition for the stability of the SDB model (Eqs. A11–A13). However, the condition we found is difficult to interpret and apply. We thus proposed a stricter but simpler sufficient condition for stability (Eq. 25). By comparing the SDB and SDBE model, we proposed that a similar constraint (${Z}_{i\times j}\ge \mathrm{0}$, Eq. 26) would also hold as a sufficient condition for the SDBE model despite the Routh–Hurwitz stability criterion being not fully tractable analytically for this model version. Numerical analyses confirmed that the proposed sufficient condition ensures stability of the SDBE model within the vast parameter space we explored. However, these sufficient conditions are very conservative and can exclude a substantial fraction of the physically meaningful and stable equilibrium points. Further, despite a clear correlation between Z_{i×j} and the damping coefficient ζ, the stability condition does not give direct insights into the oscillation behavior. How useful the stricter sufficient and necessary conditions would be in constraining model parameters – as compared to the simpler sufficient conditions – might depend on the specific model applications. Despite the potential challenges in evaluating these conditions, they can still be useful to understand the processes or parameter interactions that cause unstable EPs to occur and can guide ecologyinformed model developments.
4.3.3 Conditions leading to instability
Our numerical analysis of the SDBE model indicates that instability of equilibrium points becomes more likely with decreasing carbon inputs, increasing DOC leaching, and low process rates (Fig. 4, Table 7). All these conditions are most likely to be met in highaltitude and/or highlatitude environments. This is in line with Hararuk et al. (2015), who observed the strongest oscillations (longest time to dampen oscillations, indicative of diminishing real parts of the eigenvalues) in their calibrated fourpool model in these regions. Therefore, analytical steady states of microbial explicit SOC models applied in highaltitude and/or highlatitude environments could be unstable, and analytical steadystate solutions could thus not be used reliably for the initialization of simulation runs or prediction of SOC stocks.
Microbialexplicit SOC models aim to improve the representation of SOC dynamics by accounting for its biotic controls. At very small spatial and temporal scales, their oscillatory behavior and potential for instability can reflect relevant (micro)ecosystem processes (Manzoni and Porporato, 2007). However, if applied at larger scales, such as in Earth system models, these properties can result in unrealistic simulation outcomes (Georgiou et al., 2017; Wang et al., 2014). Here, we analyzed what processes can lead to instability in these models. By comparing the stability behavior of an archetypal microbialexplicit SOC model (the AWB model; Allison et al., 2010) with some reduced model versions and with the stability analysis of similar models in the literature, we found that instability can occur in models that assume a positive feedback between the resupply of a microbial growth substrate (i.e., DOC) and its abundance. We found that stability is (sufficiently) conditional on the balance between the sensitivity of the depolymerization rate to changes in extracellular enzymes vs. SOC concentration. Based on these analyses, we suggest that instability can be avoided by selecting specific (1) model structures, (2) kinetic formulations, and/or (3) parameter relations or values. While these approaches can differ vastly, an emerging common theme is that acknowledging ecological principles and processes can be leveraged to improve model applicability. These findings have implications for further development of microbialexplicit models and potential upscaling approaches, calling for ecologically consistent model formulations and rigorous mathematical analysis of newly introduced models.
A1 SBE model
The Jacobian matrix of the SBE model is given by
where $\frac{\partial x}{\partial y}$ is the partial derivative of x with respect to y, and $\dot{x}$ is used to denote $\frac{\text{d}x}{\text{d}t}$ so that, e.g., $\frac{\partial \dot{S}}{\partial S}=\frac{\partial}{\partial S}\left(\frac{\text{d}S}{\text{d}t}\right)$.
A1.1 Stability analysis of the abiotic equilibrium
First, we evaluate the parameter space in which the abiotic equilibrium is stable. Substituting the steadystate solutions for Q_{0} given in Table 5 into ${J}_{i}^{\mathrm{SBE}}$ (Eq. A1) yields
For Q_{0} to be stable according to the Routh–Hurwitz criterion, it is required that all the coefficients a_{i} of the characteristic polynomial of ${J}_{i}^{\mathrm{SBE}}{}_{*,\mathrm{0}}$ and, additionally, a_{1}a_{2}−a_{3} be positive. By these means, we find that the stability of Q_{0} is conditional on the following sufficient and necessary condition (Sect. S2.1 in the Supplement):
A1.2 Stability analysis of the biotic equilibrium
Next, we analyze the stability of the biotic equilibrium by evaluating the Jacobian matrix ${J}_{i}^{\mathrm{SBE}}$ (Eq. A1) according to its biotic steady states Q_{1,i} (Table 5):
To evaluate the Routh–Hurwitz criterion, it is convenient to reexpress P_{i} in terms of $\frac{\partial {P}_{i}}{\partial E}$ as follows:
where the factor x_{i} is introduced to maintain the generality of this substitution for any P_{i}, as defined in Table 2:
x_{i} has the convenient property $\mathrm{0}<{x}_{i}\le \mathrm{1}$. The ensuing derivations hold for all values of x_{i} as long as $\mathrm{0}<{x}_{i}\le \mathrm{1}$ and also for other kinetics P_{j} not explored here, for which a $\mathrm{0}<{x}_{j}\le \mathrm{1}$ can be found in order to satisfy Eq. (A5).
Substituting Eq. (A5) and ${B}_{\mathrm{1},i}^{*}=\frac{\mathit{\alpha}{y}_{B}}{\mathit{\eta}}{E}_{\mathrm{1},i}^{*}$ (Table 5) into Eq. (12), evaluated at a steady state, yields
from which, for ${E}_{\mathrm{1},i}^{*}\ne \mathrm{0}$, we find
With this definition, we obtain
From this, it can be seen that the trace of ${J}_{i}^{\mathrm{SBE}}{}_{*,\mathrm{1}}$ (the sum of the diagonal entries) is always negative since $\frac{\partial {P}_{i}}{\partial S}{}_{*,\mathrm{1}}>\mathrm{0}$ and ${x}_{i}{}_{*,\mathrm{1}}\le \mathrm{1}$; thus, a_{1}>0. Likewise, it can be shown that all remaining coefficients of the characteristic polynomial are always positive and that ${a}_{\mathrm{1}}{a}_{\mathrm{2}}{a}_{\mathrm{3}}>\mathrm{0}$ (see detailed analytical derivations in Sect. S2.1 in the Supplement). Thus, all physically meaningful biotic equilibrium points of the threepool SBE model are stable.
A2 SDB model
The Jacobian matrix ${J}_{i\times j}^{\mathrm{SDB}}$ of the SDB model evaluated around the biotic equilibrium is given by
Since only the biotic EP of the SDB model is analyzed, we dropped the subscript 1 here.
Evaluating the Routh–Hurwitz criterion for the SDB model with the chosen kinetic formulations (Tables 2 and 4) yields the sufficient and necessary condition for stability as follows:
with
and
both shown for l_{D}=0 for conciseness (see full expressions with l_{D}≥0 and detailed analysis of all coefficients of the characteristic polynomial in Sect. S2.2 in the Supplement). The sufficient condition for stability given by Eq. (25) in the main text holds for any l_{D}≥0. Beyond this sufficient condition for stability, the additional positive term $\mathit{\pi}\frac{\partial {U}_{j}}{\partial D}{}_{*}$ in Eq. (A12) might indicate the stabilizing influence of the balancing loop B_{2} (Fig. 2b) and the recycling of ENZ and MBC to SOC (compare Eq. 24).
$\frac{\partial \dot{B}}{\partial B}{}_{*,j}$ in Eq. (A13) is the lowerright entry of ${J}_{i\times j}^{\mathrm{SDB}}{}_{*}$ (Eq. A10) and is given by
For any choice of U_{j} that is linear in B (as is the case for U_{m} and U_{f} – compare Table 2), we find from solving $\frac{\text{d}B}{\text{d}t}$ (Eq. 16) at a steady state that ${y}_{m}\frac{\partial {U}_{(m,f)}}{\partial B}{}_{*}={m}_{E}+{d}_{B}$; thus $\frac{\partial \dot{B}}{\partial B}{}_{*,(m,f)}=\mathrm{0}$ so that ${Y}_{i\times (m,f)}=\mathrm{0}$ (Eq. A13). The only necessary and sufficient condition for the stability of the SDB model in these cases is thus ${X}_{i\times j}>\mathrm{0}$ (Eq. A12) (for l_{D}=0; see Sect. S2.2 in the Supplement for the corresponding necessary condition for l_{D}>0).
For linear uptake kinetics U_{l}, that is
Y_{i×l} (Eq. A13) does not vanish from Eq. (A11) (since $\frac{\partial {U}_{\mathrm{l}}}{\partial B}{}_{*}=\mathrm{0}$) or, consequently, from Eq. (A14).
Because of its additional positive components, Y_{i×j} can be positive even if the sufficient condition of Eq. (25) is not fulfilled. Thus, using U_{l} can help to ensure the positivity of all the coefficients of the characteristic polynomial and ${a}_{\mathrm{1}}{a}_{\mathrm{2}}{a}_{\mathrm{3}}>\mathrm{0}$.
A3 SDBE model
The Jacobian matrix of the fourpool SDBE model with ${L}_{S}={L}_{E}=\mathrm{0}$ and y_{B}=y_{m} (Eq. 6) (Table 4) and evaluated around the biotic equilibrium is given by
Further details on the stability analysis of the SDBE model are given in Sect. S2.3 in the Supplement.
The code used for the numerical stability analysis and comparison with the proposed sufficient stability condition is publicly available at Zenodo https://doi.org/10.5281/zenodo.12749207 (Schwarz, 2024).
All data used in this paper are from literature sources (see Table 3, Sect. 2.3 and Sect. S1 in the Supplement).
The supplement related to this article is available online at: https://doi.org/10.5194/bg2134412024supplement.
ES led the study. ES and SM conceptualized the study. ES and SG conducted the formal analysis and investigation. SM and SB assisted in the formal analysis and investigation. All the authors discussed the results together. ES wrote the original draft of the paper and produced the figures, with feedback from SM. All authors reviewed and commented on the original draft of the paper and its revisions.
The contact author has declared that none of the authors has any competing interests.
Publisher’s note: Copernicus Publications remains neutral with regard to jurisdictional claims made in the text, published maps, institutional affiliations, or any other geographical representation in this paper. While Copernicus Publications makes every effort to include appropriate place names, the final responsibility lies with the authors.
We thank Björn Lindahl for the inspiring discussions and valuable comments on the paper and two anonymous reviewers for their insightful comments and reflections.
This research has been supported by the European Research Council (ERC) under the European Union’s Horizon 2020 Research and Innovation Programme (grant no. 101001608), the Vetenskapsrådet (grant no. 202003910), and the Svenska Forskningsrådet FORMAS (grant no. 202102121).
The publication of this article was funded by the Swedish Research Council, Forte, FORMAS, and Vinnova.
This paper was edited by David Medvigy and reviewed by two anonymous referees.
Abramoff, R., Xu, X., Hartman, M., O'Brien, S., Feng, W., Davidson, E., Finzi, A., Moorhead, D., Schimel, J., Torn, M., and Mayes, M. A.: The Millennial model: in search of measurable pools and transformations for modeling soil carbon in the new century, Biogeochemistry, 137, 51–71, https://doi.org/10.1007/s1053301704097, 2018. a
Abs, E., Saleska, S., and Ferriere, R.: Microbial ecoevolutionary responses amplify global soil carbon loss with climate warming, Research Sqaure [preprint], https://doi.org/10.21203/rs.3.rs1984500/v1, 2022. a, b, c, d, e
Abs, E., Chase, A. B., and Allison, S. D.: How do soil microbes shape ecosystem biogeochemistry in the context of global change?, Environ. Microbiol., 25, 780–785, https://doi.org/10.1111/14622920.16331, 2023. a, b, c, d
Abs, E., Chase, A. B., Manzoni, S., Ciais, P., and Allison, S. D.: Microbial evolution–An underappreciated driver of soil carbon cycling, Glob. Change Biol., 30, e17268, https://doi.org/10.1111/gcb.17268, 2024. a, b, c
Allison, S. D., Wallenstein, M. D., and Bradford, M. A.: Soilcarbon response to warming dependent on microbial physiology, Nat. Geosci., 3, 336–340, https://doi.org/10.1038/ngeo846, 2010. a, b, c, d, e, f, g, h
Argyris, J. H., Faust, G., Haase, M., and Friedrich, R.: An Exploration of Dynamical Systems and Chaos: Completely Revised and Enlarged Second Edition, Springer Berlin, Heidelberg, ISBN 9783662460429, https://doi.org/10.1007/9783662460429, 2015. a, b, c, d, e
Bassiouni, M., Manzoni, S., and Vico, G.: Optimal plant water use strategies explain soil moisture variability, Adv. Water Resour., 173, 104405, https://doi.org/10.1016/j.advwatres.2023.104405, 2023. a
Bradford, M. A., Wieder, W. R., Bonan, G. B., Fierer, N., Raymond, P. A., and Crowther, T. W.: Managing uncertainty in soil carbon feedbacks to climate change, Nat. Clim. Change, 6, 751–758, https://doi.org/10.1038/nclimate3071, 2016. a, b, c
Calabrese, S., Mohanty, B. P., and Malik, A. A.: Soil microorganisms regulate extracellular enzyme production to maximize their growth rate, Biogeochemistry, 158, 303–312, https://doi.org/10.1007/s10533022008998, 2022. a, b
Crameri, F.: Geodynamic diagnostics, scientific visualisation and StagLab 3.0, Geosci. Model Dev., 11, 2541–2562, https://doi.org/10.5194/gmd1125412018, 2018a. a
Crameri, F.: Scientific colour maps (4.0), Zenodo [code], https://doi.org/10.5281/zenodo.1243862, 2018b. a
Chakrawal, A., Herrmann, A. M., Koestel, J., Jarsjö, J., Nunan, N., Kätterer, T., and Manzoni, S.: Dynamic upscaling of decomposition kinetics for carbon cycling models, Geosci. Model Dev., 13, 1399–1429, https://doi.org/10.5194/gmd1313992020, 2020. a
Chakrawal, A., Calabrese, S., Herrmann, A. M., and Manzoni, S.: Interacting Bioenergetic and Stoichiometric Controls on Microbial Growth, Front. Microbiol., 13, 859063, https://doi.org/10.3389/fmicb.2022.859063, 2022. a
Cotrufo, M. F. and Lavallee, J. M.: Soil organic matter formation, persistence, and functioning: A synthesis of current understanding to inform its conservation and regeneration, in: Advances in Agronomy, 172, 1–66, Academic Press, ISBN 9780323989534, https://doi.org/10.1016/bs.agron.2021.11.002, 2022. a
Georgiou, K., Abramoff, R. Z., Harte, J., Riley, W. J., and Torn, M. S.: Microbial communitylevel regulation explains soil carbon responses to longterm litter manipulations, Nat. Commun., 8, 1223, https://doi.org/10.1038/s4146701701116z, 2017. a, b, c, d, e, f, g, h, i, j, k, l, m
Haraldsson, H. V.: Introduction to System Thinking and Causal Loop Diagrams, Lund University, Department of Chemical Engineering, ISSN 11042877, 2004. a, b, c, d
Hararuk, O., Smith, M. J., and Luo, Y.: Microbial models with data‐driven parameters predict stronger soil carbon responses to climate change, Glob. Change Biol., 21, 2439–2453, https://doi.org/10.1111/gcb.12827, 2015. a, b, c, d, e, f, g
He, X., Abramoff, R. Z., Abs, E., Georgiou, K., Zhang, H., and Goll, D. S.: Model uncertainty obscures major driver of soil carbon, Nature, 627, E1–E3, https://doi.org/10.1038/s41586023069991, 2024. a
Horn, R. A. and Johnson, C. R.: Topics in Matrix Analysis, Cambrigde University Press, ISBN: 0521467136, 1994. a
Kuzyakov, Y.: Priming effects: Interactions between living and dead organic matter, Soil Biol. Biochem., 42, 1363–1371, https://doi.org/10.1016/j.soilbio.2010.04.003, 2010. a
Kuzyakov, Y., Friedel, J., and Stahr, K.: Review of mechanisms and quantification of priming effects, Soil Biol. Biochem., 32, 1485–1498, https://doi.org/10.1016/S00380717(00)000845, 2000. a
Lennon, J. T., Abramoff, R. Z., Allison, S. D., Burckhardt, R. M., DeAngelis, K. M., Dunne, J. P., Frey, S. D., Friedlingstein, P., Hawkes, C. V., Hungate, B. A., Khurana, S., Kivlin, S. N., Levine, N. M., Manzoni, S., Martiny, A. C., Martiny, J. B. H., Nguyen, N. K., Rawat, M., Talmy, D., ToddBrown, K., Vogt, M., Wieder, W. R., and Zakem, E. J.: Priorities, opportunities, and challenges for integrating microorganisms into Earth system models for climate change prediction, mBio, 15, e0045524, https://doi.org/10.1128/mbio.0045524, 2024. a
Malik, A. A., Puissant, J., Goodall, T., Allison, S. D., and Griffiths, R. I.: Soil microbial communities with greater investment in resource acquisition have lower growth yield, Soil Biol. Biochem., 132, 36–39, https://doi.org/10.1016/j.soilbio.2019.01.025, 2019. a
Malik, A. A., Martiny, J. B. H., Brodie, E. L., Martiny, A. C., Treseder, K. K., and Allison, S. D.: Defining traitbased microbial strategies with consequences for soil carbon cycling under climate change, ISME J., 14, 1–9, https://doi.org/10.1038/s4139601905100, 2020. a
Manzoni, S. and Porporato, A.: A theoretical analysis of nonlinearities and feedbacks in soil carbon and nitrogen cycles, Soil Biol. Biochem., 39, 1542–1556, https://doi.org/10.1016/j.soilbio.2007.01.006, 2007. a, b, c, d, e, f
Manzoni, S. and Porporato, A.: Soil carbon and nitrogen mineralization: Theory and models across scales, Soil Biol. Biochem., 41, 1355–1379, https://doi.org/10.1016/j.soilbio.2009.02.031, 2009. a
Manzoni, S., Schaeffer, S., Katul, G., Porporato, A., and Schimel, J.: A theoretical analysis of microbial ecophysiological and diffusion limitations to carbon cycling in drying soils, Soil Biol. Biochem., 73, 69–83, https://doi.org/10.1016/j.soilbio.2014.02.008, 2014. a
Manzoni, S., Chakrawal, A., and Ledder, G.: Decomposition rate as an emergent property of optimal microbial foraging, Frontiers in Ecology and Evolution, 11, 1094269, https://doi.org/10.3389/fevo.2023.1094269, 2023. a
Raupach, M. R.: Dynamics of resource production and utilisation in twocomponent biospherehuman and terrestrial carbon systems, Hydrol. Earth Syst. Sci., 11, 875–889, https://doi.org/10.5194/hess118752007, 2007. a, b, c, d, e, f, g, h, i
Richardson, G. P.: Problems with causalloop diagrams, Syst. Dynam. Rev., 2, 158–170, https://doi.org/10.1002/sdr.4260020207, 1986. a, b
Schimel, J. P. and Weintraub, M. N.: The implications of exoenzyme activity on microbial carbon and nitrogen limitation in soil: a theoretical model, Soil Biol. Biochem., 35, 549–563, https://doi.org/10.1016/S00380717(03)000154, 2003. a, b, c, d
Schwarz, E.: Numerical stability analysis of an archetypal microbialexplicit soil organic carbon model, Zenodo [code], https://doi.org/10.5281/zenodo.12749207, 2024. a
Sierra, C. A. and Müller, M.: A general mathematical framework for representing soil organic matter dynamics, Ecol. Monogr., 85, 505–524, https://doi.org/10.1890/150361.1, 2015. a, b, c, d, e
Tang, J. and Riley, W. J.: Competitor and substrate sizes and diffusion together define enzymatic depolymerization and microbial substrate uptake rates, Soil Biol. Biochem., 139, 107624, https://doi.org/10.1016/j.soilbio.2019.107624, 2019. a
Tang, J. Y. and Riley, W. J.: A total quasisteadystate formulation of substrate uptake kinetics in complex networks and an example application to microbial litter decomposition, Biogeosciences, 10, 8329–8351, https://doi.org/10.5194/bg1083292013, 2013. a, b
Tao, F., Huang, Y., Hungate, B. A., Manzoni, S., Frey, S. D., Schmidt, M. W. I., Reichstein, M., Carvalhais, N., Ciais, P., Jiang, L., Lehmann, J., Wang, Y.P., Houlton, B. Z., Ahrens, B., Mishra, U., Hugelius, G., Hocking, T. D., Lu, X., Shi, Z., Viatkin, K., Vargas, R., Yigini, Y., Omuto, C., Malik, A. A., Peralta, G., CuevasCorona, R., Di Paolo, L. E., Luotto, I., Liao, C., Liang, Y.S., Saynes, V. S., Huang, X., and Luo, Y.: Microbial carbon use efficiency promotes global soil carbon storage, Nature, 618, 981–985, https://doi.org/10.1038/s41586023060423, 2023. a, b, c, d, e
Tao, F., Houlton, B. Z., Frey, S. D., Lehmann, J., Manzoni, S., Huang, Y., Jiang, L., Mishra, U., Hungate, B. A., Schmidt, M. W. I., Reichstein, M., Carvalhais, N., Ciais, P., Wang, Y.P., Ahrens, B., Hugelius, G., Hocking, T. D., Lu, X., Shi, Z., Viatkin, K., Vargas, R., Yigini, Y., Omuto, C., Malik, A. A., Peralta, G., CuevasCorona, R., Di Paolo, L. E., Luotto, I., Liao, C., Liang, Y.S., Saynes, V. S., Huang, X., and Luo, Y.: Reply to: Model uncertainty obscures major driver of soil carbon, Nature, 627, E4–E6, https://doi.org/10.1038/s41586023070009, 2024a. a
Tao, F., Houlton, B. Z., Huang, Y., Wang, Y.P., Manzoni, S., Ahrens, B., Mishra, U., Jiang, L., Huang, X., and Luo, Y.: Convergence in simulating global soil organic carbon by structurally different models after data assimilation, Global Change Biol., 30, e17297, https://doi.org/10.1111/gcb.17297, 2024b. a
The MathWorks Inc.: MATLAB version: 9.13.0.2105380 (R2022b) Update 2, Natick, Massachusetts, United States, https://www.mathworks.com (last access: 12 July 2024), 2022. a, b
ToddBrown, K. E. O., Randerson, J. T., Post, W. M., Hoffman, F. M., Tarnocai, C., Schuur, E. A. G., and Allison, S. D.: Causes of variation in soil carbon simulations from CMIP5 Earth system models and comparison with observations, Biogeosciences, 10, 1717–1736, https://doi.org/10.5194/bg1017172013, 2013. a, b, c, d
Varney, R. M., Chadburn, S. E., Burke, E. J., and Cox, P. M.: Evaluation of soil carbon simulation in CMIP6 Earth system models, Biogeosciences, 19, 4671–4704, https://doi.org/10.5194/bg1946712022, 2022. a, b, c, d
Wang, G., Post, W. M., and Mayes, M. A.: Development of microbialenzymemediated decomposition model parameters through steadystate and dynamic analyses, Ecol. Appl., 23, 255–272, https://doi.org/10.1890/120681.1, 2013. a, b, c
Wang, G., Jagadamma, S., Mayes, M. A., Schadt, C. W., Megan Steinweg, J., Gu, L., and Post, W. M.: Microbial dormancy improves development and experimental validation of ecosystem model, ISME J., 9, 226–237, https://doi.org/10.1038/ismej.2014.120, 2015. a
Wang, Y. P., Chen, B. C., Wieder, W. R., Leite, M., Medlyn, B. E., Rasmussen, M., Smith, M. J., Agusto, F. B., Hoffman, F., and Luo, Y. Q.: Oscillatory behavior of two nonlinear microbial models of soil carbon decomposition, Biogeosciences, 11, 1817–1831, https://doi.org/10.5194/bg1118172014, 2014. a, b, c, d, e, f, g, h
Wang, Y. P., Jiang, J., ChenCharpentier, B., Agusto, F. B., Hastings, A., Hoffman, F., Rasmussen, M., Smith, M. J., ToddBrown, K., Wang, Y., Xu, X., and Luo, Y. Q.: Responses of two nonlinear microbial models to warming and increased carbon input, Biogeosciences, 13, 887–902, https://doi.org/10.5194/bg138872016, 2016. a, b, c, d, e, f
Wieder, W. R., Bonan, G. B., and Allison, S. D.: Global soil carbon projections are improved by modelling microbial processes, Nat. Clim. Change, 3, 909–912, https://doi.org/10.1038/nclimate1951, 2013. a
Wieder, W. R., Grandy, A. S., Kallenbach, C. M., and Bonan, G. B.: Integrating microbial physiology and physiochemical principles in soils with the MIcrobialMIneral Carbon Stabilization (MIMICS) model, Biogeosciences, 11, 3899–3917, https://doi.org/10.5194/bg1138992014, 2014. a
Wieder, W. R., Allison, S. D., Davidson, E. A., Georgiou, K., Hararuk, O., He, Y., Hopkins, F., Luo, Y., Smith, M. J., Sulman, B., Todd‐Brown, K., Wang, Y., Xia, J., and Xu, X.: Explicitly representing soil microbial processes in Earth system models, Global Biogeochem. Cy., 29, 1782–1800, https://doi.org/10.1002/2015GB005188, 2015. a, b, c, d
Wieder, W. R., Hartman, M. D., Sulman, B. N., Wang, Y., Koven, C. D., and Bonan, G. B.: Carbon cycle confidence and uncertainty: Exploring variation among soil biogeochemical models, Glob. Change Biol., 24, 1563–1579, https://doi.org/10.1111/gcb.13979, 2018. a, b
Wilson, C. H. and Gerber, S.: Theoretical insights from upscaling Michaelis–Menten microbial dynamics in biogeochemical models: a dimensionless approach, Biogeosciences, 18, 5669–5679, https://doi.org/10.5194/bg1856692021, 2021. a
Wutzler, T. and Reichstein, M.: Priming and substrate quality interactions in soil organic matter models, Biogeosciences, 10, 2089–2103, https://doi.org/10.5194/bg1020892013, 2013. a, b