Simulating long-term responses of soil organic matter turnover to substrate stoichiometry by abstracting fast and small-scale microbial processes: the Soil Enzyme Steady Allocation Model (SESAM; v3.0)

Understanding the coupling of nitrogen (N) and carbon (C) cycles of land ecosystems requires understanding microbial element use efficiencies of soil organic matter (SOM) decomposition. Whereas important controls of those efficiencies by microbial community adaptations have been shown at the scale of a soil pore, a simplified representation of those controls is needed at the ecosystem scale. However, without abstracting from the many details, models are not identifiable; i.e. they cannot be fitted without ambiguities to observations. There is a need to find, implement, and validate abstract simplified formulations of theses processes.

Therefore, we developed the Soil Enzyme Allocation Model (SEAM). The model explicitly represents community adaptation strategies of resource allocation to extracellular enzymes and enzyme limitations on SOM decomposition. They thus provide an abstraction from several microbial functional groups to a single holistic microbial community. Here we further simplify SEAM using a quasi-steady-state assumption for extracellular enzyme pools to derive the Soil Enzyme Steady Allocation Model (SESAM) and test whether SESAM can provide the same decadal-term predictions as SEAM.

SESAM reproduced the priming effect, the SOM banking mechanism, and the damping of fluctuations in carbon use efficiency with microbial competition as predicted by SEAM and other more detailed models. This development is an important step towards a more parsimonious representation of soil microbial effects in global land surface models.

1 Introduction

Soil organic matter (SOM) dynamics constitutes a strong link of global nutrient cycles because the microbial decomposer community has a rather strict homeostatic regulation of their stoichiometry (Sterner and Elser, 2002; Zechmeister-Boltenstern et al., 2015). Hence, understanding and properly modelling SOM dynamics is required to understand the linkage of the global element cycles of nitrogen (N) and carbon (C) (Thornton et al., 2007; Janssens et al., 2010; Zaehle and Dalmonech, 2011; Todd-Brown et al., 2012; Xu et al., 2014). The discussion about microbial carbon use efficiency (CUE), the ratio between microbial carbon uptake and microbial growth, has received attention in the literature (Wieder et al., 2013; Bradford et al., 2016; Hagerty et al., 2018; Liu et al., 2018; Fatichi et al., 2019) as an important predictor for building up SOM stocks. However, relevant processes happen at the pore scale, and processes strongly vary between microsites and sites. On the other hand we are concerned with dynamics at the pedon or ecosystem scale and strive for general predictions that do not depend on the varying details. Hence, there is a discrepancy in detail and scale between process understanding and what is feasible to implement in global models.

CUE is an emergent value that depends on microbial traits, such as maintenance requirements, and stoichiometric imbalances of the substrates that microbes feed on (Sect. 2.1.1). At the micro-pore scale, Kaiser et al. (2014) have shown that competition between microbes allows for a dampening of CUE fluctuations despite changing stoichiometry of the decomposed SOM. Specifically, by using a detailed individual-based model they showed how the succession of different microbial populations with different substrate preferences sustained high community CUE despite changing substrate stoichiometry.

In a first abstraction step models represent different microbial groups or guilds instead of single microbes (Allison, 2014; Perveen et al., 2014; Huang et al., 2018). This matches the scale of inputs and outputs to the SOM model, making it potentially feasible to implement them as submodels in global models. However, for each microbial group properties and parameters have to be determined. The increase in number of model parameters inflates the risk of equifinality, i.e. simulating the observed data for the wrong reasons, and renders model–data-integration studies difficult. For example Huang et al. (2018) used 11 microbe-related model parameters for each group.

An alternative model abstraction represents a single microbial community with adapting properties by, for example, optimizing microbial growth in the model. This concept was applied in the Soil Enzyme Allocation Model (SEAM) (Wutzler et al., 2017), which was successful in reproducing the simulated patterns of the group-based SYMPHONY model (Perveen et al., 2014).

In a further model abstraction, the effect of changing CUE can be modelled in conventional pool-based models without explicit microbial community by a growth-adapted humification coefficient and stoichiometry-dependent decomposition rates as in the PRIM model (Guenet et al., 2016). However, such a model cannot represent the changes in decomposition rate with apparent CUE depending on biomass properties (Tang and Riley, 2014).

Hence, there is a need for a model abstraction with fewer details that is still able to represent effects of stoichiometry such as priming effects due to N mining from SOM under N limitation (Kuzyakov, 2010) and the banking mechanism (Perveen et al., 2014; Wutzler et al., 2017), where excess available N is stored by SOM build-up and made available again by increased SOM decomposition during periods of N limitation. To meet this need, we previously developed SEAM (Wutzler et al., 2017) and, here, present the Soil Enzyme Steady Allocation Model (SESAM), a further simplification of SEAM.

SESAM is intended to capture the longer-term, i.e. decadal, dynamics of SOM decomposition and abstracts from short-term dynamics, i.e. shorter than seasonal, by applying the quasi-steady-state assumption (Wutzler and Reichstein, 2013) to soil extracellular enzyme pools.

The aim of this paper is to present SESAM without pre-knowledge of SEAM and show that it can reproduce the predictions of SEAM at a decadal timescale and is therefore able to simulate the priming effect due to N mining and the banking mechanism.

2 Methods

2.1 SESAM

The dynamic Soil Enzyme Steady Allocation Model (SESAM) allows exploration of the consequences of soil microbial stoichiometry for SOM cycling at the soil core to the ecosystem scale at a decadal timescale. The modelled system consists of C and N pools in SOM in a volume of soil. Such a system can be a layer of a soil profile or a laboratory incubation. It can be integrated into a larger model that explicitly represents the entire soil column and vertical transport.

Table 1 State variables and model drivers. Values refer to the reference state of the sensitivity analysis and were adopted from Wutzler et al. (2017).

^∗ Arbitrary high value so that plant uptake is constraint by $u_{{\mathrm{I}}_{\mathrm{N}},\text{max}}\left(t\right)$ (Eq. A7c).

SESAM models several SOM pools containing C and N by differential equations for their mass fluxes. Model drivers are inputs of C and N by plant litter (both aboveground litter and rhizodeposition), input of inorganic N from deposition and fertilizers, and prescribed root uptake of inorganic N (Table 1).

SESAM represents several SOM fractions by several pools that differ by their stoichiometry, and it represents changes in microbial community structure by changing preferences in degrading specific SOM pools. The litter pool, L, is rich in C; the residue pool, R, is richer in N (Fig. 1, Table 1). While Appendix A provides all the detailed model equations, the following paragraph summarizes the most important assumptions. Table A1 explains the symbols used.

Figure 1 SESAM structure: microbial biomass B produces enzymes that depolymerize substrate pools (L and R) that differ in their elemental ratios. Microbial community enzyme allocation α determines which part of the microbial community depolymerizes L versus R by producing respective enzymes E_L and E_R. Microbes take up dissolved organic matter (DOM) and use it for synthesizing new biomass or new enzymes or for catabolic respiration. A part of microbial turnover (tvr) adds to the residue pool, another part is mineralized, and another part adds to DOM and is recycled into microbial biomass. Stoichiometric imbalance between DOM and B is resolved by mineralizing the excess element or immobilizing the required element (Φ_B) from inorganic N (I) (further detailed in Fig. A1). Boxes correspond to pools, disks to fluxes. Solid lines represent fluxes of both C and N, while dotted and dashed lines represent separate C or N fluxes, respectively.

This paper presents model version 3, which in addition to the enzyme steady-state assumption differs from published SEAM (Wutzler et al., 2017) by explicitly tracking community composition, α, as a state variable; by using a modified formula for weighting optimal states near co-limitation; and by using a modified calculation of the revenue, which is used in the optimality assumption. The enzymes are assumed to be in a steady state; i.e. production of enzymes equals their turnover. This assumption allows computation of the size of the enzyme pools based on other pools, simplification of equations, and use of fewer parameters in the overall model. The microbial community is assumed to self-organize in ways to optimize growth of the entire community. This involves adjusting the relative proportion of groups that preferentially decompose litter and groups that preferentially decompose microbial residues. The enzymes, which are produced in different shares, in turn affect the decomposition of the respective pools.

Decomposition of the litter and residue pools is modelled by reverse Michaelis–Menten kinetics (Schimel and Weintraub, 2003). $\mathrm{C}/\mathrm{N}$ ratios, β_N, of fluxes are equal to the $\mathrm{C}/\mathrm{N}$ ratios of the source pools. While the $\mathrm{C}/\mathrm{N}$ ratios of the substrate pools may change over time due to changing $\mathrm{C}/\mathrm{N}$ ratio of total influxes to these pools, the $\mathrm{C}/\mathrm{N}$ ratios of biomass and enzymes are assumed to be fixed. Total enzyme allocation is modelled as a fraction, a_E, of the microbial biomass, B, per time. SESAM assumes the DOM pool to be in a quasi-steady state (Wutzler and Reichstein, 2013) and does not explicitly simulate its changes because the dynamics of the DOM pool is usually much faster than the dynamics of the other pools. Therefore, microbial uptake equals the sum of all influxes to the DOM pool, i.e. decomposition plus part of the enzyme turnover. The microbial carbon balance can be negative (i.e. uptake is smaller than maintenance and enzyme synthesis), and then microbial biomass declines.

2.1.1 Soil microbial stoichiometry

SESAM assumes that stoichiometry is one of the overarching controls of decadal-scale SOM changes. It assumes that microbial community and development of different SOM stocks adapt to changes in drivers without the need to explicitly model all the details of this adaptation.

There are three principle ways (Mooshammer et al., 2014b) in which decomposers can adapt to stoichiometric imbalance, i.e. differences in elemental composition between the requirement of feeders and its food (Sterner and Elser, 2002). First, decomposer communities can change their nutrient use efficiencies (Sinsabaugh et al., 2013) by overflow respiration or mineralizing N. For example, if there is more N in DOM uptake for biomass synthesis than constrained by other ways, such as available C, excess N will be mineralized. Such regulation of nutrient use efficiencies has large consequences for loss of nutrients from the ecosystem (Mooshammer et al., 2014a), plant nutrition and soil–plant feedback (Rastetter, 2011), and carbon sequestration in SOM (Allison, 2014; Wieder et al., 2013). Second, microbes can adapt their stoichiometric requirements. There are possible shifts in community composition between species with lower $\mathrm{C}/\mathrm{N}$ ratio, such as many bacteria, or species with high $\mathrm{C}/\mathrm{N}$ ratio, such as many fungi (Cleveland and Liptzin, 2007; Xu et al., 2013). However, the range of adjustment is quite constrained. Third, microbial community can alter the stoichiometry of uptake by adapting their preferences of degrading SOM fractions that vary in elemental composition (Moorhead et al., 2012; Li et al., 2021).

SEAM and SESAM assume that a combination of options 1 and 3 is used in a way to optimize growth, and option 2 is negligible. Modelled microbial community develops in a way so that different kinds of enzymes are produced in proportion to their revenue, i.e. the decomposition return per unit of limiting element invested into enzyme production. This microbial enzyme allocation strategy performed better in simulation experiments (Wutzler et al., 2017) than a fixed enzyme allocation or stoichiometrically optimized allocation irrespective of the amount of substrate (Moorhead et al., 2012).

While this adaptive single microbial community was a necessary step in the direction to simplify models, SEAM required two model parameters for the dynamics of the enzyme pools. These parameters are hard to measure and added complexity to model–data integration. The presented SESAM uses only one enzyme-production-related parameter, while the enzyme turnover parameter merges with the half-saturation parameter of the SOM decomposition.

2.1.2 Quasi-steady-state approximation of enzyme dynamics

SESAM abstracts from short-term dynamics of enzyme pools in SEAM by assuming that soil extracellular enzyme pools are in a quasi-steady state (Wutzler and Reichstein, 2013); i.e. pools are at levels where current enzyme production equals enzyme turnover.

It assumes that compared to the intended decadal modelling timescale, the amount of soil extracellular enzyme mass and the composition of the microbial community approach a steady state given the annually smoothed inputs and drivers of the modelled system and current SOM stocks. This yields the enzyme states in Eq. (1).

$$\begin{array}{}\text{(1)}& E_Z^{\ast }={\displaystyle \frac{{\mathit{\alpha }}_Z{a}_{\mathrm{E}}B}{k_{\mathrm{N}}}},\end{array}$$

where $Z\in \mathit{\{}R,L\mathit{\}}$ denote enzymes depolymerizing litter and residue substrate. α_Z is the proportion of biosynthesis invested into the respective enzyme (Sect. 2.1.3); a_E is the enzyme production per microbial biomass, B, per time; and k_N is the enzyme turnover rate.

This steady-state expression is used instead of the explicitly modelled enzyme pool in SEAM to simplify other equations. For example, decomposition of the residue pool now depends on biomass rather than enzyme levels (Eq. 2).

$$\begin{array}{}\text{(2a)}& {\displaystyle }{\text{dec}}_{\mathrm{R}}& {\displaystyle }=k_{\mathrm{R}}R{\displaystyle \frac{E_{\mathrm{R}}}{k_{\mathrm{m}}+E_{\mathrm{R}}}}\text{(2b)}& {\displaystyle }{\text{dec}}_{\mathrm{R}}^{\ast }& {\displaystyle }=k_{\mathrm{R}}R{\displaystyle \frac{{\mathit{\alpha }}_{\mathrm{R}}{a}_{\mathrm{E}}B}{k_{\mathrm{m}}{k}_{\mathrm{N}}+{\mathit{\alpha }}_{\mathrm{R}}{a}_{\mathrm{E}}B}},\end{array}$$

where k_m is the half-saturation constant of the original enzyme-limited decomposition equation. In the steady-state equations, half-saturation k_m and the enzyme turnover rate k_N always occur together as a product. Hence, they can be replaced by a single parameter k_mN.

We further explored two additional potential simplifying model assumptions. First the microbial biomass can be assumed to be in a quasi-steady state (SteadyB; Appendix F). Second, the mass fluxes from microbial biomass to DOM due to production and turnover of enzymes can be neglected (NoEnzFlux; Appendix E). Both assumptions yielded models with the same number of parameters. The SteadyB variant had one fewer state variable; however, the representation of biomass B by the solution of a third polynomial yielded no further simplification of model equations compared to SESAM. The NoEnzFlux variant still required the parameters for enzyme production to compute revenues and did not result in model simplification.

2.1.3 Community enzyme allocation α

Enzyme allocation α is defined as the proportion of total enzyme synthesis that is invested into production of a specific enzyme. With SESAM it is computed to be proportional to the revenue, i.e. return per investment. It is normalized to add up to unity (Eq. 3).

$$\begin{array}{}\text{(3)}& {\mathit{\alpha }}_{Z,\text{Opt}}={\displaystyle \frac{{\text{rev}}_Z}{{\sum }_{\mathit{\zeta }}{\text{rev}}_{\mathit{\zeta }}}},\end{array}$$

where rev_Z is the revenue from a given enzyme Z.

While the original SEAM computed both the return and the investment for each element $E\in (\mathrm{C},\mathrm{N})$ separately, SESAM computes the revenue based on limitation-weighted return divided by a limitation-weighted investment. The return is determined by the mineralization flux of element E from enzyme action, the investment is calculated by the amount of E required for enzyme production, and the weight for E is determined by the difference in potential growth constrained by other elements and growth constrained by E (Appendix A4). Compared to the original SEAM formulation this avoids the unreasonable high revenues when the concentration of the limiting element in enzymes is low, e.g. for an extension of SESAM which includes phosphorus limitation.

The return of an element E for enzyme Z in SESAM is the flux toward the microbial biomass from decomposition of targeted substrate, S, computed with current steady-state enzyme levels, which are determined by current enzyme production (Eq. 2). Alternatively to decomposition flux, the microbial uptake flux can define the return. This option is discussed in Appendix B, where it is argued that this only differs from the simpler decomposition approach in rare cases.

When inserting the steady-state revenue into Eq. (3), one notices that the revenue in turn depends on current α. Nevertheless explicit formulas for optimal α can be obtained for either assuming C limitation or assuming N limitation. However, situations near co-limitations then pose a problem. At the same substrate levels, a community optimized for C usage is N-limited, while a community optimized for N usage then is C-limited. In SEAM this was prevented by the inertia of explicitly modelled enzyme levels. In contrast, SESAM does not use these explicit formulas to compute α but rather computes the revenues based on current instead of the optimum enzyme allocation. To this end it models α as an additional state variable that develops towards an optimal value, α_Opt, that can be computed by Eq. (3) given the current value of α and corresponding decomposition fluxes and revenues (Eq. A13). Since α represents a property of the microbial community, it changes at the timescale of growth and microbial turnover (Eq. A12).

2.2 Increased C input simulations

In order to show the ability of SESAM to reproduce the priming effects due to N mining and the banking mechanism, we repeated the CO₂ fertilization experiment in Wutzler et al. (2017). It studied the effect of increased continuous litter C input on SOM cycling. Such increased C input is expected with elevated atmospheric CO₂ concentration. The simulations started from a steady state corresponding to initial litter C input of 400 $\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$, applied 20 % increased C inputs during years 10 to 60, and applied initial litter inputs again during the next 50 years. The litter N inputs were kept constant over time, implying an increase in the litter $\mathrm{C}/\mathrm{N}$ ratio of 20 % in annual average litter inputs. Maximum plant uptake rate of inorganic N was set to litter input rate, plus N deposition of $\mathrm{1}/\mathrm{200}$ of litterfall N was compensated by a small leaching rate.

We compared results of the following model variants:

• SEAM – baseline model with explicit representation of extracellular enzymes,

• SESAM – enzyme levels assumed in a quasi-steady state,

• SESAM-NoEnzFlux – additionally neglecting the mass flux of N from microbial biomass to DOM via enzyme production and turnover (Appendix E),

• SESAM-SteadyB – additionally microbial biomass assumed in a quasi-steady state (Appendix F).

The derivatives of the model variants were implemented in the R programming language (R Core Team, 2021) and solved, i.e. simulated over time, using the function lsoda from the package deSolve (Soetaert et al., 2010).

2.3 Substrate pulse simulations

In order to show the ability of SESAM to reproduce the damping of fluctuations in CUE with adapting microbial community, we simulated an incubation experiment. SESAM models CUE as an emergent property instead of a model parameter. With the substrate pulse scenario we simulated an experiment similar to the one in Kaiser et al. (2014), who showed that competition among microbial groups controlled the emergent CUE in a spatially explicit model.

In this experiment microbial community decomposes a carbon-rich chunk of initial litter whose initial $\mathrm{C}/\mathrm{N}$ ratio varied by simulation scenario. Continuous L input rate was set close to 0 to 40 $\mathrm{g}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$, and N leaching rate was set to zero. CUE was computed as a synthesis of microbial biomass carbon, i.e. excluding enzyme production and turnover, per carbon taken up.

Simulations were compared between SESAM, which has dynamic enzyme allocation, and a model version where we fixed community enzyme allocation α_Z=0.5 to not dynamically change over time.

We used the same R-based computational setup as in the increased C input simulations.

2.4 Sensitivity analysis

In order to explore which parameters most influence the steady-state and transient predictions of the increased C input simulation (Sect. 2.2), we performed a global sensitivity analysis using Sobol sensitivity indices (Saltelli et al., 2008) with the soboltouati method in the sensitivity R package using two samples of N = 5000.

For each parameter we prescribed prior distributions of possible parameter values (Table A1). Next, we computed the cumulative probability of initial parameters and sampled cumulative probabilities ± 10 % around this value. Two samples of cumulative probabilities were used to generate design matrices of the sensitivity method. Next we transformed the cumulative probabilities of the design matrices back to parameter values and simulated model output. Hence, we performed a global uncertainty analysis in a subspace of the entire parameter space that covered 20 % of the parameter range of each parameter. This procedure sampled larger intervals for more uncertain parameters and avoided many unreasonable parameter combinations that would occur with a global sensitivity analysis across the entire prior space not accounting for parameter correlations.

We checked robustness of the setting by repeating the analysis by sampling ± 20 % around the reference parameter values to cover 40 % of each parameter range.

While the increased C input simulations (Sect. 2.2) used a high decomposition rate of the residue pool, k_R, for demonstrating stock changes after 50 years, the sensitivity analysis and fluctuation analysis used a lower decomposition rate farther away from the edge of the parameter's prior distribution (Table A1) and simulated increased inputs from year 0 to 100. This helped to avoid severely truncated intervals of sampled parameter space around the decomposition rate.

For each model run, we computed (a) the steady-state SOM stocks and (b) the change in SOM stocks after 100 years of increased C input. The computed sensitivity indices tell about the proportion of variance in these outputs due to variation in input parameters by either the single effect of a given parameter (first-order index) or the combined effect with interaction with other parameters (total sensitivity index) (Saltelli et al., 2008; Sobol, 1990).

SESAM was recoded using the Julia programming language (Bezanson et al., 2017) using the package ModelingToolkit.jl (Ma et al., 2021) and solved using the method Tsit5 (Tsitouras, 2011) implemented in the package DifferentialEquations.jl (Rackauckas and Nie, 2017). This allowed faster simulation and the application of global sensitivity analysis.

2.5 Fluctuation nonlinearity experiment

In order to explore the possible bias in long-term predictions due to the enzyme steady-state assumption combined with the non-linearity of its effect on decomposition we simulated strong seasonal fluctuations in litter inputs causing seasonal fluctuations in elemental limitation and enzyme community allocation.

SESAM incorporates nonlinear functions. Hence, average decomposition computed with fluctuating enzyme levels will give different results compared to decomposition computed with first averaging enzyme levels. Hence, the difference between explicitly modelled enzyme levels and steady-state enzyme levels has the potential to introduce bias also in the average long-term predictions.

In order to investigate the effect of both the time-averaging and enzyme steady-state assumption we performed an experiment where we ran both SESAM and a version of SEAM that explicitly tracks enzyme pools but otherwise uses the same formulations, with the setting of the increased C input simulation (Sect. 2.2) and seasonally varying litter input. We imposed higher litter input in autumn and lower input during the rest of the year (Fig. 2) and inspected possible bias in the simulation results.

Figure 2 Seasonally fluctuating litter input rate was simulated by assuming 50 % of litter input by aboveground litterfall in autumn only. Physical decay of an aboveground litter pool with turnover time of 2 months then contributed to inputs to the SOM model. At year zero an increase in average annual litter input was prescribed.

The fluctuating litter scenario displayed stiff properties; hence we used the Vern7 method (Verner, 2010) to solve, i.e. simulate the system across time.

Figure 3 SOM stocks, here approximated by the sum of litter, L, and residue, R, develop in the same way as the “increased C input” simulation with the different model version. All model variants could simulate the liberation of N from organic R stocks with the stoichiometric imbalance with increased C inputs.

3 Results

3.1 Increased C input simulations

The imbalance in stoichiometry with increased C input was compensated by shifting enzyme production towards decomposition of the N-rich residue pool. This led to an increase in litter stocks and a decrease in residue stocks over time (Fig. 3). After input stoichiometry returned to initial conditions, the stocks slowly recovered towards the initial state.

Hence, the models simulated microbial N mining, i.e. the behaviour where under increased C input and hence N limitation, N is liberated from SOM. They also simulated the banking mechanism, i.e. the build-up of N-rich SOM stocks when N limitation was relieved.

At this timescale there were no apparent differences between the enzyme explicit SEAM and the quasi-steady-state models SESAM and its SteadyB variant. The NoEnzFlux variant lacked a refuelling of the DOM pool by the N-rich enzyme turnover and hence simulated a residue pool with higher $\mathrm{C}/\mathrm{N}$ ratio and decomposition slightly shifted towards the residue pool already in an initial steady state. However, the variant predicted the same pattern across time.

3.2 Substrate pulse simulations

CUE varied dynamically in the substrate pulse simulations (Fig. 4). During initial stages there was not enough N in the substrate, leading to overflow respiration and low CUE. The lowest CUE was found with the highest initial substrate $\mathrm{C}/\mathrm{N}$ ratios. The subsequent accumulation of microbial residues provided a source of N, which helped increase CUE again.

Figure 4 Variation in carbon use efficiency (CUE) over time with the substrate pulse simulation is more dampened, i.e. is not changing as much with changing substrate stoichiometry, if microbial community can adapt enzyme allocation compared to fixed allocation. This was true both across time (x axis) and across initial substrate $\mathrm{C}/\mathrm{N}$ ratio (line colour and line type).

The differences in CUE across time and across initial litter $\mathrm{C}/\mathrm{N}$ ratios was damped with adaptable enzyme allocation compared to fixed enzyme allocation.

3.3 Sensitivity analysis

SOM stocks and their transient changes in the increased C input simulation were only sensitive to 5 out of 14 model parameters (Fig. 5). As expected, the total stocks in a steady state were most sensitive to the potential decomposition rate, k_R, of the residue pool, which is the largest SOM pool. In addition, they were sensitive to parameters of microbial turnover, τ and ϵ_tvr, and turnover of enzymes, k_mN. Transient changes in SOM stocks were also sensitive to microbial turnover, but not to decomposition rate. Rather they were sensitive to stoichiometric parameters, specifically the $\mathrm{C}/\mathrm{N}$ ratio of microbial biomass, β_NB, and the intrinsic CUE, ϵ.

Figure 5 Simulated stock change with increased C litter input (ΔSOM) was most sensitive to $\mathrm{C}/\mathrm{N}$ ratio of microbial biomass (β_NB), while initial steady-state stocks (SOM) were most sensitive to turnover rate of the residue pool (k_R). Both were sensitive to parameters of microbial turnover. The first-order effects (triangles) were very similar to total order effects (circles) for SOM but not for ΔSOM. This means that effects on steady-state stocks were mostly direct, while the effects on stock changes were mostly sensitive to interactions of parameters.

When repeating the sensitivity analysis on a subspace that included 40 % rather than 20 % of each parameter range, the results were influenced by extreme values due to unusual parameter combinations. We observed similar total effects of SOM stocks, but first-order effects were slightly smaller. SOM stock changes now were additionally sensitive to decomposition rate, k_R, and the sensitivity to $\mathrm{C}/\mathrm{N}$ ratio of microbial biomass, β_NB, was less strong (Fig. C1 in the Appendix).

3.4 Fluctuation analysis

Simulation results differed only marginally between steady-state enzymes (SESAM) and explicit representation of enzyme level (SEAM) in the fluctuating litter input simulations (Fig. 6; sesam_annual overplotting seam_annual). This resulted from steady-state enzymes being close to the explicit enzyme levels (Appendix D). Hence, the enzyme steady-state assumptions did not introduce bias in predictions.

Figure 6 Small differences in transient N-leaching between simulation scenarios of seasonally fluctuating litter input (wiggling curves) and annually averaged litter inputs (smooth curves) were caused by different timing of litter input between the scenarios after the prescribed increase in average annual inputs at year zero. Abstraction of enzyme steady states (sesam) did not cause additional differences between the litter scenarios, as seen by overplotting.

Fluxes based on averaging litter inputs also roughly matched the average of the fluxes based on fluctuating litter inputs at a steady state (Fig. 6 before increase in C input at time t<0). Hence, also averaging litter inputs did not introduce bias at a steady state.

The largest simulated differences due to averaging litter input were observed in transient changes in the fast pools, e.g. the inorganic N pool and associated leaching (Fig. 6 after increase in C input at time t>0). This is expected because immediate transient changes depend on the timing of litter inputs within the year.

4 Discussion

4.1 Abstract from microbial details at the seasonal and pedon scale

Soil organic matter (SOM) science has experienced a paradigm shift from understanding persistence of SOM formerly on chemical SOM properties towards understanding persistence as an interactive effect of environmental conditions (Schmidt et al., 2011). SOM formation and turnover are expected to be controlled by microbial processes, such as depolymerization of SOM, respiration, mineralization/immobilization, and microbial turnover. Hence, a more mechanistic representation of microbial processes should improve dynamic SOM models and their applicability across different sites. Indeed, microbial models have shown to be able to better represent global patterns of SOM stocks (Wieder et al., 2013). However, differences in representing these microbial processes lead to qualitatively different predictions. Furthermore, it is hard to get observations as required to constrain their parameters (Marschmann et al., 2019).

Microbial processes work on pore spatial scale and hourly to daily temporal scales. In our work we pursue the hypothesis that at the pedon scale and longer decadal-term scale, stoichiometry provides one of the most important constraints (Buchkowski et al., 2019) and that at this scale we can abstract from many detailed processes while keeping the effects that are most important when averaged across pedons and years.

The presented SESAM employs the simplifying assumption of enzyme levels being close to steady-state (Sect. 2.1.2). Further, both SEAM and SESAM abstract from detailed microbial strategies to cope with varying litter stoichiometry such as internal reserves (Manzoni et al., 2021) or energetic trade-offs (Dufour et al., 2021) and neglect smoothing dynamics that occurs when explicitly modelling DOM and enzyme pools. Hence, SESAM was expected to not precisely reproduce seasonal dynamics with litter inputs due to nonlinear effects of enzymes on decomposition. It may predict overflow respiration or excess N mineralization due to stoichiometric imbalance in uptake of DOM due to short-term fluctuations in model drivers.

4.2 Possible bias due to abstraction

Because of the just-explained problems of omitting buffering capabilities of soil microbes to sudden environmental changes, we recommend driving SESAM with annually averaged model drivers. However, averaging inputs together with nonlinear functions can cause bias (Chakrawal et al., 2020; Graham et al., 2019), and also the timing of litter inputs can be important (Luu et al., 2022; Zhou et al., 2021). Specifically, the mostly concave functions of decomposition according to Michaelis–Menten kinetics yield a higher decomposition flux of the average stock (input argument to the decomposition function) compared to the average of the fluxes computed using the fluctuating stocks. Hence, we expected slightly higher decomposition rates and lower stocks with the average litter input scenario. The fluctuation analysis revealed, fortunately, that due to inherent dynamics of community pools, there was no apparent bias due to the enzyme steady-state assumption. Rather, there were only small transient deviations in averaged predictions by driving SESAM with annually averaged litter inputs compared to supplying seasonally fluctuating litter inputs (Fig. 6). This shows that SESAM may be applicable also at shorter timescales under some conditions, although the model has been primarily designed to predict the effect of decadal-term changes considering annually averaged model inputs.

4.3 Simulate patterns emerging from microbial processes

Competition between microbial groups and adaptation of the microbial community is one of the detailed processes that have been shown to exert strong control on decadal-term SOM dynamics (Kaiser et al., 2014). SESAM represents such competition by modelling adaptation of community traits, specifically the allocation of enzyme production. We performed simulation experiments to see whether the model can predict similar patterns to more detailed models. The “substrate pulse simulations” (Fig. 4) reproduced the pattern shown by Kaiser et al. (2014) of more dampened CUE with microbial adaptation, i.e. smaller range of fluctuating CUE on fluctuating substrate stoichiometry. Avoiding the low CUE with litter pulses by shifting the community to grow more on the N-rich substrates helps sequester more C to the soil pools. The “increased C input simulations” (Fig. 3) reproduced the pattern shown by Perveen et al. (2014) of liberating N from SOM stocks during conditions of imbalanced higher C inputs. Making this N available for plants helps them to avoid or delay progressive N limitation (Averill et al., 2015).

Whether increased C inputs lead to decrease in SOM due to priming or to increase in SOM due to larger input of microbial turnover to SOM is still an open question in soil science. As it probably depends on the combination of other constraints such as nutrient availability (Hicks et al., 2021; Vain et al., 2021; Feng and Zhu, 2021), we might gain new insights using the SESAM concept in model–data-integration studies.

Model–data-integration studies require observations at the modelled timescale. SESAM predicts a change in proportion of different SOM pools in response to shifting nutrient limitations. While the relative changes in SOM pools are so small that they are very hard to directly measure, changes can potentially be detected by observing changing $\mathrm{C}/\mathrm{N}$ ratio of the total SOM (vandenEnden et al., 2021; Melillo, 2019).

4.4 Model complexity and identifiability

Optimal detail or complexity of models depends on the purpose of the model and on the available data to constrain the models (Jakeman et al., 2006). Microbial processes are strongly determined by the pore-scale environment (Pot et al., 2021; Kaiser et al., 2014), and available data are mostly from the laboratory and from hourly to seasonal timescales. Larger-scale compilations are available for stoichiometry of litter, microbial biomass, and enzyme activities (Xu et al., 2013; Zechmeister-Boltenstern et al., 2015) but scarce for microbial growth and turnover rates (Spohn et al., 2016) and element use efficiencies (Manzoni et al., 2018) and very rare for different groups of microbes. Hence, microbial parameters need to be constrained by calibrating models to larger-scale observations.

There are many attempts to directly implement microbial processes into global models with introducing many free parameters (Campbell and Paustian, 2015). For example the ORCHIMIC (Huang et al., 2018) explicitly models several microbial functional types and active and dormant fractions. It introduces n_Group ⋅ 11 microbial parameters. More parameters raise the problem of model identifiability in model inversions. The model may fit the observations for the wrong reasons, and predictions then likely fail for conditions different from that of the calibration.

SESAM aims at reducing model complexity. There are in total 14 model parameters, and long-term SOM stock changes were sensitive to only a few of them. This is a more tractable number of parameters for model inversions, although there will be more parameters for temperature and moisture dependencies and transport when integrated into larger models. Because SESAM targets the decadal-term scale, decadal-term drivers and observations should also be used in SESAM model–data-integration studies. Currently, the free air enrichment experiments are running about 20 years. Thus, obtained observations are getting long enough to calibrate and test models at a decadal timescale.

4.5 Sensitivity of predictions to microbial properties

$\mathrm{C}/\mathrm{N}$ ratio of microbial biomass, β_BN, was the most sensitive parameter for SOM changes (Fig. 5). Hence, the assumption that it can be fixed because its range is rather constrained has to be revisited. We think of future developments to let it change together with other microbial properties of enzyme allocation.

The turnover rate of microbial biomass, τ, was also among the most sensitive parameters. One reason is that it contributes an important source for inorganic N when microbes are N-limited (Joly et al., 2020; Clarholm, 1985). Without predation of microbes and associated mineralization, N would be largely locked in the organic loop. In SESAM microbial turnover it is currently modelled first-order to microbial biomass and predation as a fixed factor of this turnover. A part of the turnover feeds to the DOM pool, and a part adds to the residue pool. Developments of SESAM focus on refining those processes in several ways. Predation rate is smaller at low microbial biomass when predators have to move farther between encountering microbes. Further, the stoichiometry of the parts feeding to DOM and the residue pool will be different because there is different cell material in lysed cells and cell walls. Differences in stoichiometry of microbial decay becoming DOM and decay becoming residue will have a large impact on modelled SOM stocks (Yu et al., 2020).

5 Conclusions

The application of the quasi-steady-state assumption for extracellular enzyme pools simplified a model of microbial adaptation to substrate stoichiometry. The simplified SESAM could reproduce important effects of microbial stoichiometry on SOM dynamics at a decadal timescale, specifically the priming effect, microbial N mining, the SOM banking mechanism, and dampened CUE. This development is an important step towards more parsimonious representation of soil microbial effects in global land surface models.

Appendix A: SESAM equations

For an overview of symbol definitions see Tables 1, A1, and A2.

Table A1 Model parameters and distributions used in the sensitivity analysis.

Table A2 Further symbols of quantities derived within the system.

A1 Carbon dynamics

$$\begin{array}{}\text{(A1a)}& {\displaystyle \frac{\mathrm{d}B}{\mathrm{d}t}}& {\displaystyle }={\text{syn}}_{\mathrm{B}}-{\text{tvr}}_{\mathrm{B}}\text{(A1b)}& {\displaystyle \frac{\mathrm{d}L}{\mathrm{d}t}}& {\displaystyle }=-{\text{dec}}_{\mathrm{L}}+i_{\mathrm{L}}\left(t\right)\text{(A1c)}& {\displaystyle \frac{\mathrm{d}R}{\mathrm{d}t}}& {\displaystyle }=-{\text{dec}}_{\mathrm{R}}+{\mathit{ϵ}}_{\text{tvr}}{\text{tvr}}_{\mathrm{B}}+(\mathrm{1}-{\mathit{\kappa }}_{\mathrm{E}}){\text{syn}}_{\text{Enz}},\end{array}$$

where ϵ_tvr is the fraction of microbial turnover C that is respired by predators, i_L is the litter C input to the system, and κ_E is the fraction of enzyme turnover that is transferred to the DOM instead of the R pool. Details of the specific fluxes are given below.

Total enzyme production syn_Enz, maintenance respiration r_M, and microbial turnover tvr_B are modelled as a first-order kinetics of biomass:

$$\begin{array}{}\text{(A2a)}& {\displaystyle }{\text{syn}}_{\text{Enz}}& {\displaystyle }=a_{\mathrm{E}}B\text{(A2b)}& {\displaystyle }{r}_{\mathrm{M}}& {\displaystyle }=mB\text{(A2c)}& {\displaystyle }{\text{tvr}}_{\mathrm{B}}& {\displaystyle }=\mathit{\tau }B.\end{array}$$

With assuming enzyme production and turnover to be in a quasi-steady state and reverse Michaelis–Menten kinetics for substrate decomposition (Schimel and Weintraub, 2003; Tang and Riley, 2019), substrate depolymerization is

$$\begin{array}{}\text{(A3a)}& {\displaystyle }{\text{dec}}_{S,\text{Pot}}& {\displaystyle }=k_{S}S\text{(A3b)}& {\displaystyle }{\text{dec}}_S& {\displaystyle }={\text{dec}}_{S,\text{Pot}}{\displaystyle \frac{{\mathit{\alpha }}_S{\text{syn}}_{\text{Enz}}}{k_{\text{mN}}+{\mathit{\alpha }}_S{\text{syn}}_{\text{Enz}}}},\end{array}$$

where dec_S,Pot is the decomposition flux when enzymes are saturated, $S\in \mathit{\{}L,R\mathit{\}}$ is one of the substrate pools, k_S is the maximum decomposition rate of S, α_S is the proportion of total investment into enzymes that is allocated to S (Sect. A4), and k_mN is a lumped parameter of enzyme half-saturation constant and enzyme turnover.

We assume a quasi-steady state of the DOM pool, and hence, microbial uptake equals the sum of all influxes to the DOM pool (decomposition + part of the enzyme turnover).

$$\begin{array}{}\text{(A4)}& u_{\mathrm{C}}={\text{dec}}_{\mathrm{L}}+{\text{dec}}_{\mathrm{R}}+{\mathit{\kappa }}_{\mathrm{E}}{\text{syn}}_{\text{Enz}}\end{array}$$

With C limitation, C_synBC, the carbon available for synthesis of new biomass, equals the C uptake minus the expenses for enzyme synthesis (Eq. A2a) and maintenance respiration (Eq. A2b). If this balance is positive, then a part is used for growth respiration.

$$\begin{array}{}\text{(A5a)}& {\displaystyle }{C}_{\text{synBCt}}& {\displaystyle }=u_{\mathrm{C}}-{\text{syn}}_{\text{Enz}}/\mathit{ϵ}-r_{\mathrm{M}}\text{(A5b)}& {\displaystyle }{C}_{\text{synBC}}& {\displaystyle }=\left\{\begin{array}{ll}\mathit{ϵ}{C}_{\text{synBCt}},& \text{if }{C}_{\text{synBCt}}>\mathrm{0}\\ C_{\text{synBCt}},& \text{otherwise}.\end{array}\right.\end{array}$$

A part of syn_B (Eq. A10), the C balance for biomass synthesis is used for catabolic growth respiration r_G to support biomass synthesis, and the remaining fraction ϵ, the anabolic carbon use efficiency, is used for synthesis of biomass and enzymes. ϵ is assumed to be equal for all substrates for simplicity. SESAM assumes that requirements for maintenance and enzyme synthesis must be met first. Therefore, the microbial C balance can become negative, resulting in starvation of and decline in microbial biomass.

A6a$$\begin{array}{}\text{(A6a)}& {\displaystyle }{r}_{\mathrm{G}}& {\displaystyle }=\left\{\begin{array}{ll}\frac{\mathrm{1}-\mathit{ϵ}}{\mathit{ϵ}}{\text{syn}}_{\mathrm{B}},& \text{if }{\text{syn}}_{\mathrm{B}}>\mathrm{0}\\ \mathrm{0},& \text{otherwise},\end{array}\right.\end{array}$$

where syn_B is the C balance for biomass synthesis and is given below by Eq. (A10).

A2 Nitrogen dynamics

Nitrogen fluxes are computed by dividing the respective C fluxes (Eq. A1) by the $\mathrm{C}/\mathrm{N}$ ratio, β_N, of their source.

We assumed fixed $\mathrm{C}/\mathrm{N}$ ratios ${\mathit{\beta }}_{{\mathrm{N}}_{\mathrm{B}}}$ and ${\mathit{\beta }}_{{\mathrm{N}}_{\text{Enz}}}$ of the microbial biomass and enzymes. However, substrate N pools are modelled explicitly because their $\mathrm{C}/\mathrm{N}$ ratio of the substrate pools may change over time.

$$\begin{array}{ll}\text{(A7a)}& {\displaystyle \frac{\mathrm{d}{L}_{\mathrm{N}}}{\mathrm{d}t}}& {\displaystyle }=-{\text{dec}}_{\mathrm{L}}/{\mathit{\beta }}_{\mathrm{N}\mathrm{L}}+i_{\mathrm{L}}\left(t\right)/{{\mathit{\beta }}_{\mathrm{N}}}_i{\displaystyle \frac{\mathrm{d}{R}_{\mathrm{N}}}{\mathrm{d}t}}& {\displaystyle }=-{\text{dec}}_{\mathrm{R}}/{\mathit{\beta }}_{\mathrm{N}\mathrm{R}}+{\mathit{ϵ}}_{\text{tvr}}{\text{tvr}}_{\mathrm{B}}/{\mathit{\beta }}_{{\mathrm{N}}_{\mathrm{B}}}+\\ \text{(A7b)}& {\displaystyle }& {\displaystyle }(\mathrm{1}-{\mathit{\kappa }}_{\mathrm{E}}){\text{syn}}_{\text{Enz}}/{{\mathit{\beta }}_{\mathrm{N}}}_E\text{(A7c)}& {\displaystyle }{u}_{\text{Plant,N}}& {\displaystyle }=\text{min}\left(u_{\text{Plant,N},\text{max}}\right(t),k_{\text{PlantN}}(t\left)I\right)\text{(A7d)}& {\displaystyle \frac{\mathrm{d}{I}_{\mathrm{N}}}{\mathrm{d}t}}& {\displaystyle }=+i_{{\mathrm{I}}_{\mathrm{N}}}\left(t\right)-u_{\text{Plant,N}}-l_{\mathrm{N}}{I}_{\mathrm{N}}+{\mathrm{\Phi }}_{\mathrm{N}}\text{(A7e)}& {\displaystyle }{\mathrm{\Phi }}_{\mathrm{N}}& {\displaystyle }={\mathrm{\Phi }}_{\text{Nu}}+{\mathrm{\Phi }}_{\text{NB}}+r_{\text{tvr}}/{\mathit{\beta }}_{{\mathrm{N}}_{\mathrm{B}}}\text{(A7f)}& {\displaystyle }{\mathrm{\Phi }}_{\text{Nu}}& {\displaystyle }=(\mathrm{1}-{\mathit{\nu }}_{\mathrm{N}})u_{\mathrm{N},\text{OM}}.\end{array}$$

The inorganic N pool I_N balances external inputs ($i_{{\mathrm{I}}_{\mathrm{N}}}$), leaching (l_NI_N), plant uptake (u_Plant,N), and the exchange flux with soil microbial biomass (Φ_N) (Fig. A1).

Figure A1 Several component fluxes sum to total mineralization flux in SESAM: $\mathrm{\Phi }={\mathrm{\Phi }}_{\mathrm{u}}+{\mathrm{\Phi }}_{\mathrm{B}}+{\mathrm{\Phi }}_{\text{tvr}}$. A negative mineralization flux, Φ_B, corresponds to N immobilization of inorganic pool, I, by microbial biomass, B. In addition there is mineralization during microbial turnover, Φ_tvr, and mineralization during uptake, Φ_u. Fluxes across the system boundary are inputs of inorganic N (i_I), inputs of organic N with litter, plant uptake of inorganic N, and leaching.

In addition to the mineralization–immobilization imbalance flux, Φ_NB (Eq. A11c), microbes exchange N with the inorganic pools by apparent mineralization due to soil heterogeneity, Φ_Nu (Manzoni et al., 2008), and mineralization of a part of microbial turnover, $r_{\text{tvr}}/{\mathit{\beta }}_{{\mathrm{N}}_{\mathrm{B}}}$ (Eq. A16, Sect. A5).

Potential N uptake by microbes, u_N,Pot (Eq. A8), is the sum of effective organic N uptake and potential inorganic uptake. The parallel scheme (PAR) (Manzoni et al., 2008) is used to model effective organic uptake. It accounts for an apparent mineralization at the soil core scale caused by sub-scale soil spots with high N concentration in DOM. With PAR, a part of the total organic N uptake, $(\mathrm{1}-{\mathit{\nu }}_{\mathrm{N}})u_{\mathrm{N},\text{OM}}$, is mineralized. Uptake from DOM, u_N,OM, equals the influxes to DOM multiplied by the apparent N use efficiency ν_N.

$$\begin{array}{}\text{(A8a)}& {\displaystyle }{u}_{\mathrm{N},\text{Pot}}={\mathit{\nu }}_{\mathrm{N}}{u}_{\mathrm{N},\text{OM}}+u_{\text{immN,Pot}}\text{(A8b)}& {\displaystyle }{u}_{\mathrm{N},\text{OM}}={\text{dec}}_{\mathrm{L}}/{\mathit{\beta }}_{\mathrm{N}\mathrm{L}}+{\text{dec}}_{\mathrm{R}}/{\mathit{\beta }}_{\mathrm{N}\mathrm{R}}+{\mathit{\kappa }}_{\mathrm{E}}{\text{syn}}_{\text{Enz}}/{\mathit{\beta }}_{{\mathrm{N}}_{\text{Enz}}}\text{(A8c)}& {\displaystyle }{u}_{\text{immN,Pot}}=i_{\mathrm{B}\mathrm{N}}{I}_{\mathrm{N}},\end{array}$$

where $\mathrm{C}/\mathrm{N}$ ratios β_NL and β_NR are calculated using current C and N substrate pools: ${\mathit{\beta }}_{\mathrm{N}\mathrm{L}}=L/L_{\mathrm{N}}$. Note that u_N is the potential microbial N uptake using the potential immobilization flux. The actual net inorganic flux, Φ_NB, is computed taking other limitations into account (Eq. A11c).

The N constraint on biomass synthesis is potential microbial N uptake minus enzyme synthesis. Converted to C units this reads $C_{\text{synB}}\le {\mathit{\beta }}_{{\mathrm{N}}_{\mathrm{B}}}{N}_{\text{synBN}}$.

A9a$$\begin{array}{}\text{(A9a)}& {\displaystyle }{N}_{\text{synBN}}& {\displaystyle }=u_{\mathrm{N},\text{Pot}}-{\text{syn}}_{\text{Enz}}/{\mathit{\beta }}_{{\mathrm{N}}_{\text{Enz}}}\end{array}$$

Equation (A9a) assumes that the N taken up is only used for enzyme production and biomass synthesis. A potential contribution of N to maintenance processes is neglected.

A3 Imbalance fluxes of microbes limited by C or N

There is a constraint on the synthesis of new biomass by each chemical element. In SESAM synthesis follows the minimum of these constraints (Eq. A10).

$$\begin{array}{}\text{(A10)}& {\text{syn}}_{\mathrm{B}}=\text{min}(C_{\text{synBC}},{\mathit{\beta }}_{{\mathrm{N}}_{\mathrm{B}}}{N}_{\text{synBN}})\end{array}$$

The elements in excess then are lost by imbalance fluxes (Eq. A11) so that the mass balance is closed. Excess N is mineralized (M_Imb), and excess C is respired by overflow respiration (r_O).

$$\begin{array}{}\text{(A11a)}& {\displaystyle }{r}_{\mathrm{O}}& {\displaystyle }=u_{\mathrm{C}}-({\text{syn}}_{\mathrm{B}}+r_{\mathrm{G}}+{\text{syn}}_{\text{Enz}}/\mathit{ϵ}+r_{\mathrm{M}})\text{(A11b)}& {\displaystyle }{M}_{\text{ImbN}}& {\displaystyle }=u_{\mathrm{N}}-({\text{syn}}_{\mathrm{B}}/{\mathit{\beta }}_{{\mathrm{N}}_{\mathrm{B}}}+{\text{syn}}_{\text{Enz}}/{\mathit{\beta }}_{{\mathrm{N}}_{\text{Enz}}})\text{(A11c)}& {\displaystyle }{\mathrm{\Phi }}_{E\mathrm{B}}& {\displaystyle }=M_{\text{ImbE}}-u_{\text{immE,Pot}}\end{array}$$

The actual mineralization–immobilization flux Φ_EB for element E is the difference between excess mineralization of E and the potential immobilization flux. With substrate N limitation, Φ_NB will be negative (N immobilization), whereas if microbes are limited by C availability, Φ_NB will be positive (N mineralization). If the required immobilization is larger than potential immobilization ($-{\mathrm{\Phi }}_{\text{NB}}>u_{\text{immN,Pot}}$) then stoichiometry must be balanced by overflow respiration.

A4 Community composition

Microbes in SESAM allocate a proportion α_Z of their total enzyme investments to the synthesis of enzyme Z. This enzyme allocation coefficient, α_Z, reflects the community composition, i.e. the mix of species growing on different portions of the organic matter.

SESAM models composition, α_Z, as a vector of state variables that develop over time towards the optimal composition, α_ZOpt. The timescale of this change is assumed to be equal to the timescale of microbial turnover, τ, and biomass synthesis, syn_B (Eq. A12).

$$\begin{array}{}\text{(A12)}& {\displaystyle \frac{\mathrm{d}{\mathit{\alpha }}_Z}{\mathrm{d}t}}=({\mathit{\alpha }}_{Z\text{Opt}}-{\mathit{\alpha }}_Z)(\mathit{\tau }+|{\text{syn}}_{\mathrm{B}}|/B)\end{array}$$

Community can change fast either if it is growing fast or if it is decaying fast. Hence, both terms are considered in Eq. (A12).

SESAM3 adopts the revenue strategy where investment in enzyme synthesis is proportional to its revenue (Wutzler et al., 2017). It differs from SEAM by considering the weights also in the computation of the investment.

$$\begin{array}{ll}\text{(A13a)}& {\displaystyle }{\mathit{\alpha }}_{Z\text{Opt}}& {\displaystyle }={\displaystyle \frac{{\text{rev}}_Z}{{\sum }_{\mathit{\zeta }}{\text{rev}}_{\mathit{\zeta }}}}\approx {\displaystyle \frac{{\text{rev}}_Z^{\ast }}{{\sum }_{\mathit{\zeta }}{\text{rev}}_{\mathit{\zeta }}^{\ast }}}\text{(A13b)}& {\displaystyle }{\text{rev}}_Z^{\ast }& {\displaystyle }={\displaystyle \frac{\text{limitation-weighted\hspace{0.17em} return}}{\text{limitation-weighted\hspace{0.17em} investment}}}\text{(A13c)}& {\displaystyle }{\text{investment}}_Z^{\ast }& {\displaystyle }={\mathit{\alpha }}_Z{a}_{\mathrm{E}}B(w_{\mathrm{C}}+w_{\mathrm{N}}/{\mathit{\beta }}_{{\mathrm{N}}_{\text{Enz}}}){\displaystyle }{\text{return}}_Z^{\ast }& {\displaystyle }=\sum _E{\text{return}}_{Z\mathrm{E}}^{\ast }{w}_{\mathrm{E}}\\ \text{(A13d)}& {\displaystyle }& {\displaystyle }=(w_{\mathrm{C}}+w_{\mathrm{N}}/{{\mathit{\beta }}_{\mathrm{N}}}_Z){\text{dec}}_{Z_{\mathrm{B}}},\end{array}$$

where α_Z is the current community enzyme allocation, rev_Z is the revenue from given enzyme depolymerizing substrates $Z\in \mathit{\{}L,R\mathit{\}}$ of organic matter and ${\text{rev}}_Z^{\ast }$ its version computed with enzymes assumed in a steady state, and ${\text{dec}}_{Z_{\mathrm{B}}}$ is the decomposition due to biomass-produced enzymes (Eq. A3). The elemental limitation weights, w_E, occur in both the return and the investment; therefore they do not need to be normalized to one.

$$\begin{array}{}\text{(A14)}& {\text{rev}}_Z^{\ast }={\displaystyle \frac{k_{Z}Z}{k_{\text{mN}}+{\mathit{\alpha }}_Z{a}_{\mathrm{E}}B}}{\displaystyle \frac{w_{\mathrm{C}}+w_{\mathrm{N}}/{{\mathit{\beta }}_{\mathrm{N}}}_Z}{w_{\mathrm{C}}+w_{\mathrm{N}}/{\mathit{\beta }}_{{\mathrm{N}}_{\text{Enz}}}}}\end{array}$$

The unnormalized weight of an element limitation, w_E, decreases exponentially with the excess of biosynthesis flux constrained by the given element only over the actual biosynthesis flux constrained by all elements. This excess in potential biosynthesis flux is normalized by microbial turnover in order to derive a unitless quantity. Compared to the SEAM version 2 formulation of weights (Wutzler et al., 2017), this formulation works better for starving microbial community with negative biomass synthesis and can be extended to more than two limiting elements.

$$\begin{array}{}\text{(A15)}& w_{\mathrm{E}}=\text{exp}\left(-\mathit{\delta }{\displaystyle \frac{C_{\text{synBE}}-{\text{syn}}_{\mathrm{B}}}{{\text{tvr}}_{\mathrm{B}}}}\right),\end{array}$$

where δ, arbitrarily set to 40, controls the steepness of the transition between states limited by different elements. C_synBE denotes the available biosynthesis flux in microbial biomass carbon equivalents given the limitation of element E (Eqs. A5 and A9a), e.g. for nitrogen, $C_{\text{synBN}}={\mathit{\beta }}_{{\mathrm{N}}_{\mathrm{B}}}{N}_{\text{synBN}}$.

Compared with SEAM, already a small C limitation causes an increased preference for the C-rich labile pool, i.e. lower α_R. This is because the elemental N limitation is divided by $\mathrm{C}/\mathrm{N}$ ratio in Eq. (A13).

A5 Fate of microbial turnover

During microbial turnover, a part (1−ϵ_tvr) of microbial biomass is mineralized, e.g. by grazing.

$$\begin{array}{}\text{(A16)}& r_{\text{tvr}}=(\mathrm{1}-{\mathit{ϵ}}_{\text{tvr}}){\text{tvr}}_{\mathrm{B}}\end{array}$$

A respective proportion of N ($r_{\text{tvr}}/{\mathit{\beta }}_{{\mathrm{N}}_{\mathrm{B}}}$) is also mineralized and transferred to the inorganic N pools.

The remainder of the microbial turnover goes to the residue pool. The current SESAM version ignores the part that enters the DOM pool and is taken up again by living microbial biomass. This corresponds to an effective uptake rate, assuming that the effects of this DOM flux on pools cancel in their parameterizations. This shortcut leads to a joint small underestimation of microbial turnover, uptake, and CUE. Investigating the effect of this simplifying assumption on isotopic tracers is an outlook.

Appendix B: Alternative return of microbial uptake

Instead of taking the entire decomposition flux as return, one could account for the mineralization–immobilization pathway and the fact that during this path, part of the decomposition flux is routed away from microbial biomass.

$$\begin{array}{}\text{(B1a)}& {\displaystyle }{\text{return}}_{Z\mathrm{E}}& {\displaystyle }={\text{dec}}_{SZ}{\mathit{\nu }}_{\text{TE}}\text{(B1b)}& {\displaystyle }{\mathit{\nu }}_{\text{TE}}& {\displaystyle }={\mathit{\nu }}_E+(\mathrm{1}-{\mathit{\nu }}_E)p_{\text{immo},E},\end{array}$$

where ν_TE is the total nutrient efficiency that includes the uptake by the mineralization–immobilization pathway, and p_immo,E is the ratio of microbial uptake from the pool of the inorganic form of the element E to the sum of all losses from this pool.

This leads to updated equations of return, revenue, and community composition (Eq. B2).

$$\begin{array}{}\text{(B2)}& {\text{return}}_S^{\ast }={\text{dec}}_{S_{\mathrm{B}}}(w_{\mathrm{C}}+w_{\mathrm{N}}{\mathit{\nu }}_{\text{TN}}/{{\mathit{\beta }}_{\mathrm{N}}}_S)\end{array}$$

Notice that ν_TN depends via its dependence on p_immo,E on many aspects of the current state. This makes reasoning about the system more difficult.

When one element is clearly limiting, then the returns in both the numerator and the denominator in the computation of α_Z (Eq. A13a) are dominated by only a single term. In this case the ν_TN can be factored out and cancels. Hence, the computed α_Z equals the version computed by adopting the decomposition flux as the return.

We argue that the case of clear co-limitation is quite rare. Depending on fluctuations in litter input and soil heterogeneity, the microbial community at a given time and a given spot is usually limited by one of the elements. Therefore, SESAM currently adopts the simpler version of the return (Eq. A13).

Appendix C: Robustness of sensitivity parameters

Results of repeated sensitivity analysis on a larger parameter subspace (Fig. C1) were very similar to the original sensitivity analysis (Fig. 5).

Figure C1 Modification of Fig. 5 for sensitivity analysis repeated on a larger subspace encompassing 40 % of parameter ranges.

Appendix D: Supplementary figures of the fluctuation analysis

The following figures help to understand the result of Sect. 3.4 that enzyme steady-state assumptions did not introduce bias in predictions despite the non-linearity of decomposition with enzyme levels.

The aboveground litter inputs in autumn caused time-lagged responses and smoothed responses in the modelled soil properties (Fig. D1). This smoothing behaviour is similar to a daily temperature signal travelling down from top soil to deeper soil layers where the signal is delayed and smoothed.

Figure D1 Model responses to a spike of aboveground litterfall in autumn, i_Lagr, are damped. There is a cascade of delays in litter input from litter layer and roots to soil, i_L; change in litter substrate pool, L; increased share of microbial community depolymerizing this pool, α_L; and an increase in enzyme concentration, E_L. All fluxes were normalized by dividing by their maximum, except α_L, whose range additionally was scaled to 0.5 because it varied only marginally across seasons.

Due to this smoothing and lagging behaviour, the simulated steady-state enzyme levels closely tracked the explicit enzyme levels (Fig. D2).

Figure D2 Enzyme levels, here shown for residue degrading enzyme, E_R, differ only marginally across simulations of explicit enzymes (seam) and steady-state enzymes (sesam).

Appendix E: No mass flux by enzymes

SESAM does not explicitly represent enzyme pools. However, the mass fluxes across the enzyme pool from biomass to DOM and to the residue pool are represented.

A model variant “NoEnzFlux” has been implemented, where the enzyme pools are still part of the revenue computation, but mass fluxes across the enzyme pools are neglected. This has been accomplished by using syn_Enz=0 instead of Eq. (A2a) and using a_EB directly in computation of decomposition (Eq. A3) and subsequently in computation of revenue (Eq. A13).

Appendix F: Steady state of microbial biomass

Here we derive equations for microbial biomass in a quasi-steady state.

By setting dB=0, syn_B=tvr_B.

Carbon available for biomass synthesis, syn_B, is the minimum for each element at limitation, by either C or N. Hence, we compute the steady-state biomass for all limitations and then take the minimum $B_{\mathrm{s}}=\text{min}(B_{\mathrm{s}\mathrm{C}},B_{\mathrm{s}\mathrm{N}})$.

F1 Carbon limitation

With enzymes in a quasi-steady state, the uptake from enzyme turnover equals enzyme production, a_EB. While enzyme investments are subtracted from uptake for computing biomass synthesis, a part κ_E occurs in uptake.

$$\begin{array}{}\text{(F1a)}& {\displaystyle }\mathit{ϵ}{C}_{\text{syn}\mathrm{B}\mathrm{C}}=\mathit{\tau }B\text{(F1b)}& {\displaystyle }{u}_{\mathrm{C}}-{\text{syn}}_E/\mathit{ϵ}-r_{\mathrm{M}}=\mathit{\tau }/\mathit{ϵ}B\text{(F1c)}& {\displaystyle }\begin{array}{rl}& {\displaystyle \frac{{\mathrm{d}}_{\mathrm{L}}(\mathrm{1}-\mathit{\alpha })a_{\mathrm{E}}B}{k_{\text{mN}}+(\mathrm{1}-\mathit{\alpha })a_{\mathrm{E}}B}}+{\displaystyle \frac{{\mathrm{d}}_{\mathrm{R}}\mathit{\alpha }{a}_{\mathrm{E}}B}{k_{\text{mN}}+\mathit{\alpha }{a}_{\mathrm{E}}B}}+{\mathit{\kappa }}_{\mathrm{E}}{a}_{\mathrm{E}}B\\ & -{\displaystyle \frac{a_{\mathrm{E}}B}{\mathit{ϵ}}}-mB=\mathit{\tau }/\mathit{ϵ}B\end{array}\text{(F1d)}& {\displaystyle }\begin{array}{rl}& {\mathrm{d}}_{\mathrm{L}}(\mathrm{1}-\mathit{\alpha })a_{\mathrm{E}}(k_{\text{mN}}+\mathit{\alpha }{a}_{\mathrm{E}}B)+{\mathrm{d}}_{\mathrm{R}}\mathit{\alpha }{a}_{\mathrm{E}}(k_{\text{mN}}\\ & +(\mathrm{1}-\mathit{\alpha })a_{\mathrm{E}}B)=\left[\mathit{\tau }/\mathit{ϵ}+m+(\mathrm{1}/\mathit{ϵ}-{\mathit{\kappa }}_{\mathrm{E}})a_{\mathrm{E}}\right]c_{\mathrm{1}}\end{array},\end{array}$$

with abbreviations

$$\begin{array}{}\text{(F2a)}& {\displaystyle }{\mathrm{d}}_{\mathrm{L}}& {\displaystyle }={\text{dec}}_{\text{L,Pot}}=k_{\mathrm{L}}L\text{(F2b)}& {\displaystyle }{\mathrm{d}}_{\mathrm{R}}& {\displaystyle }={\text{dec}}_{\text{R,Pot}}=k_{\mathrm{R}}R\text{(F2c)}& {\displaystyle }{c}_{\mathrm{1}}& {\displaystyle }=(k_{\text{mN}}+(\mathrm{1}-\mathit{\alpha }\left)a_{\mathrm{E}}B\right)(k_{\text{mN}}+\mathit{\alpha }{a}_{\mathrm{E}}B)\text{(F2d)}& {\displaystyle }{\mathit{\tau }}_{\mathit{ϵ}m}& {\displaystyle }=\mathit{\tau }/\mathit{ϵ}+m+(\mathrm{1}/\mathit{ϵ}-{\mathit{\kappa }}_{\mathrm{E}})a_{\mathrm{E}}.\end{array}$$

This results in a square equation. If there is no real positive solution, biomass cannot be sustained, otherwise the maximum of the two roots gives the required steady-state biomass, B_sC.

$$\begin{array}{}\text{(F3a)}& {\displaystyle }\mathrm{0}& {\displaystyle }=a{B}^{\mathrm{2}}+bB+c\text{(F3b)}& {\displaystyle }a& {\displaystyle }=-{\mathit{\tau }}_{\mathit{ϵ}m}\mathit{\alpha }(\mathrm{1}-\mathit{\alpha })a_{\mathrm{E}}^{\mathrm{2}}\text{(F3c)}& {\displaystyle }b& {\displaystyle }=a_{\mathrm{E}}^{\mathrm{2}}\mathit{\alpha }(\mathrm{1}-\mathit{\alpha })({\mathrm{d}}_{\mathrm{L}}+{\mathrm{d}}_{\mathrm{R}})-{\mathit{\tau }}_{\mathit{ϵ}m}{k}_{\text{mN}}{a}_{\mathrm{E}}\text{(F3d)}& {\displaystyle }c& {\displaystyle }=k_{\text{mN}}{a}_{\mathrm{E}}\left[(\mathrm{1}-\mathit{\alpha }){\mathrm{d}}_{\mathrm{L}}+\mathit{\alpha }{\mathrm{d}}_{\mathrm{R}}\right]-{\mathit{\tau }}_{\mathit{ϵ}m}{k}_{\text{mN}}^{\mathrm{2}}\end{array}$$

F2 Nitrogen limitation

For N limitation, the potential immobilization flux does not cancel, and a complex cubed polynomial equation results.

$$\begin{array}{}\text{(F4a)}& {\displaystyle }\mathit{ϵ}{C}_{\text{synBN}}=\mathit{\tau }B\text{(F4b)}& {\displaystyle }{N}_{\text{synBN}}=\mathit{\tau }B/{\mathit{\beta }}_{\mathrm{B}}\text{(F4c)}& {\displaystyle }\begin{array}{rl}& {\mathit{\nu }}_{\mathrm{N}}\left({\text{dec}}_{\mathrm{L}}/{\mathit{\beta }}_{\mathrm{L}}+{\text{dec}}_R/{\mathit{\beta }}_{\mathrm{R}}+{\mathit{\kappa }}_{\mathrm{E}}{\text{tvr}}_E/{\mathit{\beta }}_E\right)\\ & +u_{\text{imm,Pot}}-{\text{syn}}_E/{\mathit{\beta }}_E=\mathit{\tau }B/{\mathit{\beta }}_{\mathrm{B}}\end{array}\text{(F4d)}& {\displaystyle }\begin{array}{rl}& {\displaystyle \frac{{\mathrm{d}}_{\text{LN}}(\mathrm{1}-\mathit{\alpha })a_{\mathrm{E}}B}{k_{\text{mN}}+(\mathrm{1}-\mathit{\alpha })a_{\mathrm{E}}B}}+{\displaystyle \frac{{\mathrm{d}}_{\text{RN}}\mathit{\alpha }{a}_{\mathrm{E}}B}{k_{\text{mN}}+\mathit{\alpha }{a}_{\mathrm{E}}B}}+{\displaystyle \frac{{\mathit{\kappa }}_{\mathrm{E}}{a}_{\mathrm{E}}B}{{\mathit{\beta }}_E}}\\ & +{\displaystyle \frac{u_{\text{imm,Pot}}}{{\mathit{\nu }}_{\mathrm{N}}}}-{\displaystyle \frac{a_{\mathrm{E}}B}{{\mathit{\nu }}_{\mathrm{N}}{\mathit{\beta }}_E}}={\displaystyle \frac{\mathit{\tau }}{{\mathit{\nu }}_{\mathrm{N}}{\mathit{\beta }}_{\mathrm{B}}}}B\end{array}\text{(F4e)}& {\displaystyle }\begin{array}{rl}& {\mathrm{d}}_{\text{LN}}(\mathrm{1}-\mathit{\alpha })a_{\mathrm{E}}B(k_{\text{mN}}+\mathit{\alpha }{a}_{\mathrm{E}}B)+{\mathrm{d}}_{\text{RN}}\mathit{\alpha }{a}_{\mathrm{E}}B\\ & \times (k_{\text{mN}}+(\mathrm{1}-\mathit{\alpha }\left)a_{\mathrm{E}}B\right)+u_{\mathit{\nu }\mathrm{N}}{c}_{\mathrm{1}}={\mathit{\tau }}_{\mathrm{N}}B{c}_{\mathrm{1}}\end{array},\end{array}$$

with abbreviations

$$\begin{array}{}\text{(F5a)}& {\displaystyle }{\mathit{\tau }}_{\mathrm{N}}& {\displaystyle }={\displaystyle \frac{\mathit{\tau }}{{\mathit{\nu }}_{\mathrm{N}}{\mathit{\beta }}_{\mathrm{B}}}}+\left({\displaystyle \frac{\mathrm{1}}{{\mathit{\nu }}_{\mathrm{N}}}}-{\mathit{\kappa }}_{\mathrm{E}}\right){\displaystyle \frac{a_{\mathrm{E}}}{{\mathit{\beta }}_E}}\text{(F5b)}& {\displaystyle }{u}_{\mathit{\nu }\mathrm{N}}& {\displaystyle }=u_{\text{imm,Pot}}/{\mathit{\nu }}_{\mathrm{N}}\text{(F5c)}& {\displaystyle }{\mathrm{d}}_{\text{LN}}& {\displaystyle }={\text{dec}}_{\text{L,Pot}}/{\mathit{\beta }}_{\mathrm{L}}\text{(F5d)}& {\displaystyle }{\mathrm{d}}_{\text{RN}}& {\displaystyle }={\text{dec}}_{\text{R,Pot}}/{\mathit{\beta }}_{\mathrm{R}}\text{(F5e)}& {\displaystyle }{c}_{\mathrm{1}}& {\displaystyle }=(k_{\text{mN}}+(\mathrm{1}-\mathit{\alpha }\left)a_{\mathrm{E}}B\right)(k_{\text{mN}}+\mathit{\alpha }{a}_{\mathrm{E}}B).\end{array}$$

This results in a cubic equation. Its second root is real and gives the steady-state biomass, B_sN.

$$\begin{array}{}\text{(F6a)}& {\displaystyle }\mathrm{0}& {\displaystyle }=a{B}^{\mathrm{3}}+b{B}^{\mathrm{2}}+cB+d\text{(F6b)}& {\displaystyle }a& {\displaystyle }=-{\mathit{\tau }}_{\mathrm{N}}\mathit{\alpha }(\mathrm{1}-\mathit{\alpha })a_{\mathrm{E}}^{\mathrm{2}}\text{(F6c)}& {\displaystyle }b& {\displaystyle }=a_{\mathrm{E}}^{\mathrm{2}}\mathit{\alpha }(\mathrm{1}-\mathit{\alpha })\left({\mathrm{d}}_{\text{LN}}+{\mathrm{d}}_{\text{RN}}+u_{\mathit{\nu }\mathrm{N}}\right)-{\mathit{\tau }}_{\mathrm{N}}{k}_{\text{mN}}{a}_{\mathrm{E}}\text{(F6d)}& {\displaystyle }c& {\displaystyle }=a_{\mathrm{E}}{k}_{\text{mN}}\left((\mathit{\alpha }-\mathrm{1}){\mathrm{d}}_{\text{LN}}+\mathit{\alpha }{\mathrm{d}}_{\text{RN}}+u_{\mathit{\nu }\mathrm{N}}\right)-{\mathit{\tau }}_{\mathrm{N}}{k}_{\text{mN}}^{\mathrm{2}}\text{(F6e)}& {\displaystyle }d& {\displaystyle }=k_{\text{mN}}^{\mathrm{2}}{u}_{\mathit{\nu }\mathrm{N}}\end{array}$$

While steady-state biomass can be computed and passed to other equations that involve biomass, these other equations are not simplified.