Diffusion limitations and Michaelis–Menten kinetics as drivers of combined temperature and moisture effects on carbon fluxes of mineral soils

CO₂ production in soils responds strongly to changes in temperature and moisture, but the magnitude of such responses at different timescales remains difficult to predict. Knowledge of the mechanisms leading to the often observed interactions in the effects of these drivers on soil CO₂ emissions is especially limited. Here we test the ability of different soil carbon models to simulate responses measured in soils incubated at a range of moisture levels and cycled through 5, 20, and 35 ^∘C. We applied parameter optimization methods while modifying two structural components of models: (1) the reaction kinetics of decomposition and uptake and (2) the functions relating fluxes to soil moisture. We found that the observed interactive patterns were best simulated by a model using Michaelis–Menten decomposition kinetics combined with diffusion of dissolved carbon (C) and enzymes. In contrast, conventional empirical functions that scale decomposition rates directly were unable to properly simulate the main observed interactions. Our best model was able to explain 87 % of the variation in the data. Model simulations revealed a central role of Michaelis–Menten kinetics as a driver of temperature sensitivity variations as well as a decoupling of decomposition and respiration C fluxes in the short and mid-term, with general sensitivities to temperature and moisture being more pronounced for respiration. Sensitivity to different model parameters was highest for those affecting diffusion limitations, followed by activation energies, the Michaelis–Menten constant, and carbon use efficiency. Testing against independent data strongly validated the model (R²=0.99) and highlighted the importance of initial soil C pool conditions. Our results demonstrate the importance of model structure and the central role of diffusion and reaction kinetics for simulating and understanding complex dynamics in soil C.

1 Introduction

Soils are a main component of the global carbon (C) cycle, storing ca. 2200 Pg of C in the top 100 cm according to recent estimates (Batjes, 2014). This soil C pool is dynamic and often exists in a non-equilibrium state as the result of an imbalance between input and output C fluxes, in which case it will act either as a C sink or as a C source over time. Changes in the speed at which soil organisms decompose soil organic matter (SOM) and mineralize soil organic carbon (SOC) into CO₂ are one way in which an imbalance can occur, producing a net sink or source of atmospheric CO₂.

It is well known that SOC mineralization and resulting CO₂ fluxes are highly sensitive to variations in soil temperature and moisture (Hamdi et al., 2013; Moyano et al., 2013). As a result, feedback effects, either positive or negative, are expected to occur from the interaction between climate change and global soil C stocks (Crowther et al., 2016; Davidson and Janssens, 2006; Kirschbaum, 2006). However, the direction and magnitude of such feedbacks at the global scale remain uncertain. Increased soil respiration with a resulting net loss of soil C, and thus a positive climate feedback, is expected with the warming of permafrost soils and the drying of wetland soils. But there is still large uncertainty and a lack of consensus regarding the long-term response to climate variability of soils that are non-saturated, non-frozen, and dominated by a mineral matrix (Crowther et al., 2016), i.e. soils found under most forests, grasslands, and agricultural lands.

Future predictions of soil C dynamics require the use of mathematical models. Early soil C models and most still in use are based on first-order decay of multiple C pools, with temperature and moisture having general non-interactive effects on decay rates (Rodrigo et al., 1997). When appropriately calibrated, these models do well at simulating soil respiration fluxes of soils under relatively stable conditions. They are often developed to approximate long-term steady-state conditions under specific land uses. They are also capable of fitting long-term trends of soil C loss, such as data from long-term bare fallow where all litter input has stopped (Barré et al., 2010). However, they lack a theoretical basis justifying their basic assumptions of pool partitioning and decay mechanisms. They also generally need calibration for specific soil types or land cover types and often fail to properly simulate observed short- and mid-term variability in soil respiration. Some of the most relevant observations these models have failed to reproduce include changes (typically a dampening) of temperature sensitivities of decomposition over time (Hamdi et al., 2013), non-linear responses to soil moisture content (Borken and Matzner, 2009), and changes in decomposition rates in response to variations in concentrations of organic matter (Blagodatskaya and Kuzyakov, 2008). Such model shortcomings, which reflect missing or wrongly simulated processes, create a difficult-to-quantify uncertainty in global long-term predictions of soil C and its feedback to climate change. It is therefore unclear if first-order models can predict long-term changes in C stocks under more dynamic (and therefore realistic) environmental conditions.

Second-order models have a more realistic basic structure compared to conventional first-order models, since they simulate organic matter decomposition as a reaction between SOC and decomposers (i.e. a microbial or enzyme pool). This single but fundamental change in decomposition kinetics strongly affects predicted long-term changes in soil C, largely as a result of the dynamics of the decomposer pool, which itself can respond to temperature in a number of ways (Wutzler and Reichstein, 2008). Second-order models also lead to more complex dynamics of short- to mid-term soil respiration, with apparent temperature sensitivities that vary over time, more in line with many observations.

The temporal variability in the response of decomposition to moisture is most evident in the strong respiration pulses after dry soils are re-wetted, known as the Birch effect (Birch, 1958). But studies have shown that a successful simulation of such pulses requires the incorporation of additional mechanisms, namely the explicit representation of a bio-available C pool, such as dissolved organic matter (DOC), and a moisture regulation of decomposer's access to this pool that may differ from the moisture regulation on the decomposition reaction itself (Lawrence et al., 2009; Zhang et al., 2014).

The response of soil respiration to temperature and moisture is highly dynamic, both spatially and temporally (Hamdi et al., 2013; Moyano et al., 2012). Moisture and temperature interactions have been observed in a number of experimental studies (Craine and Gelderman, 2011; Rey et al., 2005; Suseela et al., 2012; Wickland and Neff, 2008), but neither consistent trends nor general explanatory theories have been identified. Improving our understanding of these interactions is a crucial step towards increasing confidence in models and important for interpreting modelling and experimental results (Crowther et al., 2016; Tang and Riley, 2014). Identifying the model structures and parameterizations that can best represent these interactive effects has been attempted by very few studies (Sierra et al., 2015, 2017).

The objectives of this study are to compare the ability of different soil C modelling approaches to reproduce temperature and moisture interactive effects on soil carbon fluxes and thus to gain insight into mechanisms underlying the observed responses. With the hypothesis that a more mechanistic model will be better capable of simulating such interactions, we compare different model structures, testing first-order, second-order, and Michaelis–Menten reaction kinetics in combination with an explicit simulation of diffusion fluxes. We then compare the best diffusion model with versions based on common empirical moisture relationships.

2 Observational data

Measurements of the interaction effects of temperature and moisture on soil respiration fluxes were obtained by incubating a crop field soil at several fixed levels of soil moisture and variable levels of temperature over a period of ca. 6 months, as detailed in the following.

Soils from 0 to 20 cm depth were sampled at Versailles, France, from the “Le Closeaux” experimental field plot, cultivated with wheat until 1992 and with maize since 1993. Mean annual temperature and rainfall are 10 ^∘C and 640 mm, respectively. The soil is classified as Luvisol (FAO) silt loam (26 % sand, 59 % silt, 15 % clay) containing no carbonates. Organic carbon contents at the start of the incubation were 1.2 % in weight. Soil samples were prepared for elemental analysis (C, N) using a planetary ball mill (3 min at 500 rpm). C concentrations were measured using a CHN auto-analyser (NA 1500, Carlo Erba).

Sampled soils were thoroughly mixed, sieved at 2 mm, and stored moist at 4 ^∘C in plastic bags with holes for aeration for 10 days. Soils were then put in small plastic cylinders containing the equivalent of 90 g dry soil. To ensure a high and equal water conductivity, all samples were compacted to a bulk density of 1.4 g cm⁻³. The resulting soil porosity was 0.45.

All samples were brought to a pF of 4.2, corresponding to about 12.5 % volumetric moisture. Three replicate samples were then adjusted to each of the moisture levels (2 %, 5.5 %, 12.5 %, 16 %, 20 %, 23.5 %, 27 %, 30.5 %, 34 %, 38 %, and 45 %) by adding water or air drying. These values range from air-dry to saturated. Immediately thereafter, the plastic cylinders were put in 500 mL jars containing a small amount of water on the bottom (except for the 2 % and 5.5 % moisture) to prevent soil drying and equipped with a lid and a rubber septum for gas sampling. Because of the extremely low respiration rates, samples with 2 % moisture were placed in 125 mL jars containing 170 g of soil.

To minimize post-disturbance effects, samples were pre-incubated at 4 ^∘C for 10 days. The samples were then cycled through incubation temperatures following the sequence 5–20–35–5–20–35 ^∘C, thus applying two temperature cycles to each sample. This was done in order to capture possible hysteresis of temperature effects and to reduce the covariance between a temperature response and substrate depletion (helping constrain model parameters). Soil respiration was calculated at every temperature step by measuring the amount of CO₂ accumulated in flask headspace. For this, samples were flushed with CO₂-free air and left to accumulate CO₂ for 3 to 74 days. The variable accumulation time was chosen so that sufficient CO₂ accumulated for the micro gas chromatographer measurements (at least 100 ppm), thus depending on the soil temperature and moisture content. After the accumulation time, an air sample was taken from each soil sample headspace, and respiration rates were calculated as the accumulated amount over the accumulation time. Samples were incubated for a total period of ca. 6 months (Fig. 1).

As shown in Fig. 1, the timing of temperature treatments was not equal for all samples, with some temperature steps missing at low moisture levels. This was partly due to the time required for CO₂ concentrations in the flask headspace to reach detectable limits, the time necessary for carrying out measurements and human error. However, while important for a statistical comparison between treatments, such differences are of little consequence when looking at model performance and the fit between model and data, which is the focus of this study.

Figure 1 Graphical representation of the incubated soil samples showing the fixed levels of moisture content and the times at different temperatures.

3 Modelling approach

3.1 Structure and state variables

We started with a basic soil C model with the following state variables: a bio-unavailable polymeric C pool (C_P); a bio-available dissolved C pool (C_D); a microbial C pool (C_M); and two extracellular enzyme C pools (C_E), one representing the enzyme fraction at the decomposition site ($C_{{\mathrm{E}}_{\mathrm{D}}}$) and one the fraction at the microbial site ($C_{{\mathrm{E}}_{\mathrm{M}}}$). With this model we assume two conceptual soil spaces that are separated by a diffusion barrier, one being the site of decomposition and the other the site of microbial uptake and enzyme production (Fig. 2). This model thus closely follows Manzoni et al. (2016) and otherwise builds on other published microbial models (Allison et al., 2010; Schimel and Weintraub, 2003). We refer to those studies for general assumptions and application of this type of model. Aspects specific to this study are described below.

The rates of change of the model state variable were defined as

$$\begin{array}{ll}\text{(1)}& {\displaystyle }& {\displaystyle \frac{\mathrm{d}{C}_{\mathrm{P}}}{\mathrm{d}t}}=F_{{\mathrm{L}}_{\mathrm{S}}\mathrm{P}}+F_{\mathrm{MP}}-F_{\mathrm{PD}},{\displaystyle }& {\displaystyle \frac{\mathrm{d}{C}_{\mathrm{D}}}{\mathrm{d}t}}=F_{{\mathrm{L}}_{\mathrm{M}}\mathrm{D}}+F_{\mathrm{PD}}+F_{{\mathrm{E}}_{\mathrm{D}}\mathrm{D}}+F_{{\mathrm{E}}_{\mathrm{M}}\mathrm{D}}-F_{\mathrm{DM}}\\ \text{(2)}& {\displaystyle }& {\displaystyle }-F_{{\mathrm{DR}}_{\mathrm{G}}}-F_{{\mathrm{DE}}_{\mathrm{M}}},\text{(3)}& {\displaystyle }& {\displaystyle \frac{\mathrm{d}{C}_{\mathrm{M}}}{\mathrm{d}t}}=F_{\mathrm{DM}}-F_{\mathrm{MP}}-F_{{\mathrm{MR}}_{\mathrm{M}}},\text{(4)}& {\displaystyle }& {\displaystyle \frac{\mathrm{d}{C}_{{\mathrm{E}}_{\mathrm{D}}}}{\mathrm{d}t}}=F_{{\mathrm{E}}_{\mathrm{M}}{\mathrm{E}}_{\mathrm{D}}}-F_{{\mathrm{E}}_{\mathrm{D}}\mathrm{D}},\text{(5)}& {\displaystyle }& {\displaystyle \frac{\mathrm{d}{C}_{{\mathrm{E}}_{\mathrm{M}}}}{\mathrm{d}t}}=F_{{\mathrm{DE}}_{\mathrm{M}}}-F_{{\mathrm{E}}_{\mathrm{M}}{\mathrm{E}}_{\mathrm{D}}}-F_{{\mathrm{E}}_{\mathrm{M}}\mathrm{D}},\end{array}$$

where F represents the flux of C from one pool to another as indicated by the subscripts, so that F_PD is the flux from the polymeric pool to the dissolved pool. The subscripts L_S and L_M denote input of structural and metabolic litter (as defined by Parton et al., 1987), which for simulating the incubated soils were set to zero, and R_M and R_G are microbial growth and maintenance respiration.

Figure 2 Diagram showing C pools and fluxes, as well as the points of diffusion limitations. Second-order decay may refer also to Michaelis–Menten reaction kinetics. Variations of this scheme were tested in this study.

3.2 Decomposition and microbial uptake

The flux of C_P to C_D, F_PD, represents decomposition of organic matter, a process that in soils is largely driven by the activity of microorganisms, which produce exo-enzymes that catalyse the decomposition reaction. U_D represents the total uptake flux by microbes of the water-soluble decomposed pool C_D (microbes being the reaction “catalysers”). Conventional soil C models simulate decomposition as a first-order decay reaction. However, more realistic models can be built by using either simple second-order or Michaelis–Menten reaction kinetics. Thus, optional ways of modelling both F_PD and U_D include

$$\begin{array}{}\text{(6)}& {\displaystyle }& {\displaystyle }F=V\left[R\right],\text{(7)}& {\displaystyle }& {\displaystyle }F=V\left[R\right]\left[C\right],\text{(8)}& {\displaystyle }& {\displaystyle }F={\displaystyle \frac{V\left[R\right]\left[C\right]}{K+\left[R\right]}},\text{(9)}& {\displaystyle }& {\displaystyle }F={\displaystyle \frac{V\left[R\right]\left[C\right]}{K+\left[C\right]}},\end{array}$$

where F is the flux, V a base reaction rate, K the half-saturation constant, R the reactant, and C the catalyst. The “reverse” Michaelis–Menten (Eq. 9) was applied by Schimel and Weintraub (2003) as an alternative for improving model stability and is included here for completeness.

The value for V is not equivalent among these equations, differing by several orders of magnitude. As a result, different parameters were used for V in each case, namely V_Dm, V_Dmr, V_D1, and V_D2. Similarly, parameters K_D and K_De were used for K in Eqs. (8) and (9), respectively. The terms [R] and [C] are concentrations of C_P and $C_{{\mathrm{E}}_{\mathrm{D}}}$. In the case of uptake, the parameters are respectively V_U, K_U, C_D, and C_M. The four approaches for reaction kinetics were tested in order to find the best fit between model and data, as described in Sect. 4.

3.3 Diffusive fluxes

Diffusion fluxes depend on a concentration difference, a diffusivity term, and the distance over which diffusion occurs (Manzoni et al., 2016). For the purpose of modelling diffusion in soils, values of diffusivity and diffusion distances are required that best average or represent the actual underlying soil complexity. For practical purposes, we combined these two values into a single calibrated parameter, a conductance (g₀), representing the compound effects of diffusivity and distance. This was done because the values of the latter are unconstrained (from lack of information), and their effects are inversely correlated, so simultaneous calibration would lead to a problem of parameter identifiability. The moisture-scaled conductance (g), which in our model is assumed equal for the C_D and C_E pools, is then given by

$$\begin{array}{}\text{(10)}& {\displaystyle }g=g_{\mathrm{0}}{d}_{\mathit{\theta }},\end{array}$$

where d_θ is a function of soil volumetric water content (VWC or θ):

$$\begin{array}{}\text{(11)}& {\displaystyle }{d}_{\mathit{\theta }}=(\mathit{\varphi }-{\mathit{\theta }}_{\mathrm{th}}{)}^m{\left({\displaystyle \frac{\mathit{\theta }-{\mathit{\theta }}_{\mathrm{th}}}{\mathit{\varphi }-{\mathit{\theta }}_{\mathrm{th}}}}\right)}^n,\end{array}$$

where ϕ is pore space, and n and m are calibrated parameters (Hamamoto et al., 2010; Manzoni et al., 2016), which are variable and were also calibrated in this study. θ_th is the percolation threshold for solute diffusion, for which Manzoni and Katul (2014) reported a value of −15 MPa. This value was not optimal in our case, so θ_th was also calibrated. The diffusive flux of enzyme C between the microbial and the decomposition spaces is then calculated as

$$\begin{array}{}\text{(12)}& {\displaystyle }{F}_{{\mathrm{E}}_{\mathrm{M}}{\mathrm{E}}_{\mathrm{D}}}=g(C_{{\mathrm{E}}_{\mathrm{M}}}-C_{{\mathrm{E}}_{\mathrm{D}}}).\end{array}$$

Diffusion limitations also affect the amount of the dissolved pool (C_D) available for microbial uptake. Instead of dividing C_D into a pool for each space, the conductance, g, was used as a multiplier of the base uptake rate, V_U (Eqs. 6–9). This served to reduce the number of model pools and parameters while still retaining a diffusivity limitation on this flux.

3.4 Microbial and enzyme dynamics

U_D is split into F_DM, $F_{{\mathrm{DR}}_{\mathrm{G}}}$, and $F_{{\mathrm{DE}}_{\mathrm{M}}}$, representing the fluxes of C_D going to C_M, R_G, and $C_{{\mathrm{E}}_{\mathrm{M}}}$, respectively. These fluxes are defined as

$$\begin{array}{}\text{(13)}& {\displaystyle }& {\displaystyle }{F}_{\mathrm{DM}}=U_{\mathrm{D}}{f}_{\mathrm{ug}}\left(\mathrm{1}-f_{\mathrm{ge}}\right),\text{(14)}& {\displaystyle }& {\displaystyle }{F}_{{\mathrm{DR}}_{\mathrm{G}}}=U_{\mathrm{D}}\left(\mathrm{1}-f_{\mathrm{ug}}\right),\text{(15)}& {\displaystyle }& {\displaystyle }{F}_{{\mathrm{DE}}_{\mathrm{M}}}=U_{\mathrm{D}}{f}_{\mathrm{ug}}{f}_{\mathrm{ge}},\end{array}$$

where f_ug represents the fraction of uptake going to growth, otherwise known as microbial growth efficiency or carbon use efficiency, and f_ge is the fraction of growth going to enzyme production. Enzyme production thus depends here on uptake rather than on microbial biomass. This approach follows the assumption that microbes produce enzymes only when new carbon is available and save resources otherwise. C_M goes to either maintenance respiration or the C_P pool according to

$$\begin{array}{}\text{(16)}& {\displaystyle }& {\displaystyle }{F}_{\mathrm{MP}}=C_{\mathrm{M}}{r}_{\mathrm{md}}\left(\mathrm{1}-f_{\mathrm{mr}}\right),\text{(17)}& {\displaystyle }& {\displaystyle }{F}_{{\mathrm{MR}}_{\mathrm{M}}}=C_{\mathrm{M}}{r}_{\mathrm{mr}},\end{array}$$

where r_md is the rate of microbial death and r_mr is the rate of microbial maintenance respiration. The breakdown of enzymes going to the C_D pool is determined by the rate of enzyme decay, r_ed, as

$$\begin{array}{}\text{(18)}& {\displaystyle }& {\displaystyle }{F}_{{\mathrm{E}}_{\mathrm{D}}\mathrm{D}}=C_{{\mathrm{E}}_{\mathrm{D}}}{r}_{\mathrm{ed}},\text{(19)}& {\displaystyle }& {\displaystyle }{F}_{{\mathrm{E}}_{\mathrm{M}}\mathrm{D}}=C_{{\mathrm{E}}_{\mathrm{M}}}{r}_{\mathrm{ed}}.\end{array}$$

3.5 Temperature effects

Reaction rates ($V_{\mathrm{U}},V_{\mathrm{D}},K_{\mathrm{U}},and{K}_{\mathrm{D}}$ in Eqs. 6–9), decay and respiration rates (r_ed, r_md, r_mr) are temperature sensitive and calculated from their reference values following an Arrhenius type temperature response:

$$\begin{array}{}\text{(20)}& {\displaystyle }r=r_{\mathrm{ref}}\exp \left(-{\displaystyle \frac{E_{\mathrm{a}}}{R}}\left({\displaystyle \frac{\mathrm{1}}{T}}-{\displaystyle \frac{\mathrm{1}}{T_{\mathrm{ref}}}}\right)\right),\end{array}$$

where r is the temperature-modified rate, r_ref the reference rate at temperature T_ref, T temperature in kelvin, E_a the activation energy, and R the universal gas constant. Three parameters were used for E_a: E_{a_m} and E_{a_e} for microbial and enzyme decay rates, respectively, and E_{a_V} for other reaction rates. Temperature thus affects the rates of decomposition and uptake, the half-saturation constant in the Michaelis–Menten equation, and the rates of microbial and enzyme decay. Apparent activation energies – describing the observed temperature relationship, both in measurements and model data – were obtained by fitting an Arrhenius equation to the temperature–flux relationship at each level of moisture. E_a was calculated for measured respiration, modelled respiration (R_G+R_M), and modelled decomposition (F_PD).

4 Model calibration and comparisons

Calibrated and non-calibrated parameters for all models are given in the Supplement (Tables S1, S2, and S3). Whenever possible, fixed parameters as well as lower and upper bounds for calibrated parameters (Table S1) were set according to values reported in the literature (e.g. Hagerty et al., 2014; Li et al., 2014; Price and Sowers, 2004). Equilibrium conditions were not assumed at the start of the experimental procedure, as such a state is unlikely for samples that have been processed and disturbed. Therefore, initial conditions were obtained by also optimizing the fractions of initial carbon pool sizes (f_P, f_D, f_M). Total organic C was set equal to the measured value. Models were calibrated by optimizing parameters to best fit the measured soil respiration data described in Sect. 2. Each model was calibrated by fitting a single set of parameters simultaneously to all the incubation data (Table S3). For this, the model was run to reproduce each sample treatment, i.e. the applied incubation times and temperatures for each level of moisture (Fig. 1). Accumulated soil respiration amounts were then calculated to match those from the observed data. Measured and simulated data from all samples were then combined, and an overall model cost was calculated using the root mean square error (RMSE) and a weighting term, as described below.

Parameter optimization was carried out in two steps. We first explored parameter spaces using a Latin hypercube of parameter values. For this we randomly selected unique parameter sets from a uniform distribution over each parameter range (R function randomLHS, package lhs; Stein, 1987) to obtain 30 000 parameter sets. Model costs were then obtained by running models with each set. In the second step we used the Nelder–Mead algorithm (as implemented in the function modFit in package FME of the R programming language; R Development Core Team, 2016; Soetaert and Petzoldt, 2010) with initial parameter values being the set from the previous step with the lowest model cost. For the cost calculations we used an error term (“err” argument to FME function modCost) to weight the residuals. The error was calculated as the normalized (0–1) standard deviation of measured values at each combination of temperature and moisture, with 0.1 added to avoid an unreasonable weighting of measurements with near-zero errors.

For a visual inspection of the model–data fits, we plotted both the measured and the model relationship between soil respiration and moisture, soil respiration and temperature, and apparent activation energy (E_a) and moisture content.

Table 1 Different model versions with their weighted and unweighted root mean squared errors (RMSEs, in units mg C kg soil⁻¹ h⁻¹) and R² after parameter calibration. F_PD: decomposition flux; U_D: dissolved C uptake flux; 1: first-order kinetics; 2: second-order kinetics; M: Michaelis–Menten kinetics; M_r: reverse Michaelis–Menten kinetics.

4.1 Comparison of reaction kinetics

Models were named according to their decomposition kinetics followed by the uptake kinetics and the moisture function, using the abbreviations: 1 for first order, 2 for second order, M for Michaelis–Menten, M_r for reverse Michaelis–Menten, dif for diffusion, psi for water potential function, and sat for water saturation function. Alternative reaction kinetics leading to fluxes F_PD and U_D were compared in diffusion-based models using different combinations of Eqs. (6)–(9). Specifically, we compared first order for decomposition and uptake (11-dif), second order for decomposition and uptake (22-dif), and Michaelis–Menten decomposition with all combinations of uptake (M1-dif, M2-dif, and MM-dif). In addition, we tested reverse Michaelis–Menten decomposition with second-order uptake (M_r2-dif). We then evaluated the model–data fit based on RMSE values as well as on a visual inspection of the plotted relationships. A “best” model was then selected for further analysis.

4.2 Comparison of moisture regulations: diffusion vs. empirical

A second model comparison was carried out to test the impact of different approaches for modelling moisture effects. For this we modified the model M2-dif (Table 1), removing diffusion fluxes and adding empirical moisture functions. This consisted in removing all diffusion effects (so that $C_{{\mathrm{E}}_{\mathrm{M}}}$ and $C_{{\mathrm{E}}_{\mathrm{D}}}$ were replaced by a single C_E pool and the uptake rate, V_U, was no longer modified by g) and adding a function to scale (i.e. multiply) the decomposition flux, F_PD. This approach is equivalent to the conventional way used to model moisture effects on soil C fluxes. Two alternative moisture scaling functions were tested (Moyano et al., 2013), one based on relative water saturation (M2-sat) and the other on water potential (M2-wp):

$$\begin{array}{}\text{(21)}& {\displaystyle }& {\displaystyle }f\left({\mathit{\theta }}_{\mathrm{S}}\right)=a{\mathit{\theta }}_{\mathrm{S}}-b{\mathit{\theta }}_{\mathrm{S}}^{\mathrm{2}},\text{(22)}& {\displaystyle }& {\displaystyle }f\left(\mathrm{\Psi }\right)=max\left\{\begin{array}{l}min\left\{\begin{array}{l}\mathrm{1}-{\displaystyle \frac{\left[{\log }_{\mathrm{10}}\left(\mathrm{\Psi }\right)-{\log }_{\mathrm{10}}\left({\mathrm{\Psi }}_{\mathrm{opt}}\right)\right]}{\left[{\log }_{\mathrm{10}}\left({\mathrm{\Psi }}_{\mathrm{th}}\right)-{\log }_{\mathrm{10}}\left({\mathrm{\Psi }}_{\mathrm{opt}}\right)\right]}}\\ \mathrm{1}\end{array}\right.\\ \mathrm{0}\end{array}\right.,\end{array}$$

where θ_S is relative water saturation; Ψ is soil water potential; and a, b, Ψ_opt, and Ψ_th are fitted parameters. The latter two represent the optimal water potential for decomposition and a threshold water potential with values of −0.03 and −15 MPa, respectively. Water potential was calculated based on Campbell (1974) and Cosby et al. (1984). a and b are empirical parameters and were calibrated.

4.3 Model steady state, sensitivity analysis, and validation

Equations for steady state were derived by setting the rate of change in the state variables to zero in Eqs. (1)–(5) (where the flux terms are replaced by their respective equations) and then solving for the state variables. This was performed in Python using the “SymPy” package (Meurer et al., 2017).

A sensitivity analysis was carried out on model parameters using the “sensFun” function from the R package FME, which perturbs each parameter individually by a small amount. We ran the model as above, i.e. simulating the incubation and using daily output. Daily sensitivities were then averaged to obtain an overall value. Sensitivity values were calculated for the C_P pool alone, as this pool represents the largest fraction of soil C.

For model validation, we used soil respiration data from the study by Rey et al. (2005), where a Mediterranean oak forest soil was incubated for 1 month in a full factorial design at 100 %, 80 %, 60 %, 40 %, and 20 % of water holding capacity and at 30, 20, 10, and 4 ^∘C. This soil differed from the one used for model calibration in at least three aspects: the amount of organic C (7 %), soil pore space (65 %), and texture (classified as silty clay loam). The optimized set of parameters from the model M2-dif was used with the exception of the initial fraction of C pools (f_P, f_D, f_M) and the percolation threshold (θ_th), which we calibrated against the new data (Nelder–Mead calibration). The former was required since we had no information with which to estimate the microbial, dissolved, and enzyme C for this study and information regarding an initial soil steady state was also lacking. In the case of θ_th, we assumed that this parameter is especially sensitive to variations in soil texture and structure. Although in previous studies it has been determined to be equal to a water potential of −15 MPa (Manzoni and Katul, 2014), this value did not provide a good fit when applied to the validation data.

5 Results

5.1 Reaction kinetics

The calibrated values for all models are shown in Table S3. Using different reaction kinetics resulted in variations in model performance as measured by RMSE (Table 1). Changes in RMSE were more sensitive to the kinetics of decomposition (F_PD), with models using M and M_r decomposition kinetics resulting in lower RMSE values compared to first- and second-order kinetics. In terms of uptake kinetics, both first- and second-order kinetics performed better than Michaelis–Menten kinetics.

Models M1-dif, M2-dif, and M_r2-dif all showed a good fit to the data, with the first two having a slightly higher R². Thus, selecting a “best” model necessarily remains partially subjective. A visual comparison shows some weaknesses and strengths in each case. M1-dif and M_r2-dif better captured the variability in the data along the respiration axis at 35 ^∘C (Fig. S1 in the Supplement), while M2-dif more closely captured the relationship at 20 ^∘C and thus the temperature sensitivities (Fig. S2). We selected the model M2-dif (R²=0.87, Fig. 3) as the “best” model, since it better represents the actual mediation of uptake by microbial biomass when compared to the model M1-dif. We also had no theoretical reason to prefer M_r to M decomposition. The decomposition and uptake equations of the model M2-dif are then

$$\begin{array}{}\text{(23)}& {\displaystyle }& {\displaystyle }{F}_{\mathrm{PD}}=V_{\mathrm{D}}{C}_{{\mathrm{E}}_{\mathrm{D}}}{C}_{\mathrm{P}}/(K_{\mathrm{D}}+C_{\mathrm{P}})\text{(24)}& {\displaystyle }& {\displaystyle }{U}_{\mathrm{D}}=C_{\mathrm{D}}{C}_{\mathrm{M}}{V}_{\mathrm{U}}g.\end{array}$$

Figure 3 Model vs. measured accumulated CO₂ of incubated soil samples. Colour depicts the range of volumetric water content (VWC). The model R² is 0.87.

5.2 Moisture regulation

Replacing diffusion effects with empirical moisture scalars followed by recalibration decreased model performance compared to a diffusion-based model, both when using relative water saturation (M2-sat) and when using water potential (M2-wp) functions (Table 1). Although empirical functions were able to approximate the shape of the respiration–moisture relationship at 20 ^∘C, they were unable to capture the variation of this response at higher and colder temperatures, as seen in the measurements and simulated by diffusion base models (Fig. 4). Diffusion-based models more accurately simulated a linear relationship between respiration and moisture at lower temperatures and a steep increase followed by a plateau at high temperatures, with an intermediate response seen at 20 ^∘C.

Figure 4 Relationship of soil respiration with volumetric soil moisture. Results shown over three temperature levels (5, 20, 35 ^∘C) for the observed data (obs) and three model versions (M2-dif, M2-wp, and M2-sat). Lines are smooth loess fits depicting the mean relationship.

5.3 Temperature sensitivities

Figures 5 and 6 show the apparent temperature sensitivities fitted to observations and modelled fluxes at different moisture levels and for two temperature ranges, 5–20 and 20–35 ^∘C. Figure 5 compares different reaction kinetics, and Fig. 6 different moisture functions. Michaelis–Menten decomposition outperformed first- and second-order kinetics when simulating the variability in E_a observed along the moisture axis as well as the differences observed between colder (5–20 ^∘C) and warmer (20–35 ^∘C) temperature ranges. The model M2-dif closely followed the observed E_a values, which were near 100 kJ at colder temperatures and in the 30–70 kJ range at warmer temperatures. Models M2-sat and M2-wp captured the large differences between temperature ranges but did not simulate the variability along the moisture axis as well as diffusion-based models.

Figure 5 Temperature sensitivities of respiration and decomposition fluxes, showing activation energy (E_a) fitted using two temperature ranges (5–20 and 20–35 ^∘C) and the equivalent Q₁₀ derived for a 10 ^∘C range. Plotted are observed respiration data (R-obs) and three models with different reaction kinetics (R-M2-dif, R-11-dif, and R-22-dif). The sensitivities of the decomposition flux from the model M2-dif is included for comparison (D-M2-dif).

6 Model steady state, sensitivity analysis, and validation

Model steady-state equations are provided in the Supplement. For 20 ^∘C, 30 % VWC, 1.2 g d⁻¹ C input, and 30 cm soil depth (z), the equilibrium sizes of the model pools are 2560, 37, 120, and 4 g C for the C_P, C_D, C_M, and $C_{{\mathrm{E}}_{\mathrm{D}}}$ pools, respectively. These values are stable over most of the moisture range and increase exponentially only at very low soil moisture (data not shown). A similar pattern was observed for temperature, with the C_P pool increasing towards high values only at temperatures near 0 ^∘C. The same pool showed little sensitivity to changes in C input.

Figure 6 Equivalent to Fig. 5 but showing observational data (R-obs) next to models with different moisture functions (R-M2-dif, R-M2-wp, and R-M2-sat).

Table 2 shows the averaged values from the sensitivity analysis done on the model C_P pool. High sensitivities were found for g₀ and n, indicating the importance of diffusion fluxes. Large effects were also seen for the activation energy parameters, denoting a strong general effect of temperature. Also high were the sensitivities to K_D and f_ug, reflecting the importance of Michaelis–Menten kinetics for decomposition and carbon use efficiently, respectively. Low sensitivities were found for rates of microbial and enzyme decay.

Table 2 Parameters of the model M2-diff, calibrated and non-calibrated, with results of a sensitivity analysis (sens). Sens shows a relative measure of the sensitivity of the model C_P pool to small perturbations in the parameter values. Values are rounded to two significant digits.

Simulation of the incubated soil from the study of Rey et al. (2005) resulted in a very high fit to the validation data after calibration of initial SOC fractions and θ_th, with an RMSE of 0.09 in fluxes that were almost an order of magnitude higher than those used for calibration and a model R² of 0.99 (Fig. 7). This was reflected in a generally good agreement between the relationships of model and observations with moisture (Fig. 8) and temperature (Fig. 9).

7 Discussion

The interaction often observed in the effects of temperature and moisture on the cycling of soil C is an indicator of the complex nature of soil systems. Such responses are often ignored, particularly by modellers trying to minimize model complexity and derive functions that are easy to parameterize, but also by experimentalists focusing on finding an invariable response to a single factor. But a careful consideration of the nature of soils suggests that interactions should be expected, something that becomes evident in multi-factorial experiments as well as in field measurements. Here we found clear interactive effects in our experimental observations, adding to the evidence that fixed empirical temperature and moisture scalars, as used in conventional soil C models, are inappropriate for simulating the variability often found in natural conditions.

Figure 7 Model vs. measured accumulated CO₂ after simulating the experiment from Rey et al. (2005). Colour depicts the range of volumetric water content (VWC). The model R² is 0.99.

Figure 8 Relationship of soil respiration with volumetric soil moisture shown for the model M2-dif and observations from the validation data (obs). Results are shown over four temperature levels.

Since the total amount of soil C was equal among samples and its relative change in the 6 months of incubation was small, we expected that second-order kinetics would do as well as Michaelis–Menten kinetics. But using Michaelis–Menten increased the R² by ca. 5 % compared to second- and first-order kinetics. This combined with the fact that the model was highly sensitive to a change in K_D, more than to V_D, would indicate that Michaelis–Menten kinetics are in fact important for explaining soil C flows. Indeed, even in this case where the C_P pool is relatively invariant, the outcome of a strong temperature effect modifying K_D (E_a of 94 kJ) cannot be reproduced by second-order kinetics.

Figure 9 Apparent activation energy (E_a) and equivalent Q₁₀ (for the temperature range 15–25 ^∘C) during the validation step. Values fitted to observed respiration (R-obs) as well as modelled respiration (R-M2-dif) and decomposition (D-M2-dif).

The relative importance of different processes was also shown by the model parameter sensitivity values. It is perhaps not surprising that some of the highest values were related to diffusion and temperature, since these were the two factors that varied in our experiment. However, these factors also vary considerably in natural ecosystems and largely drive changes in decomposition rates. No strong correlations between the effects of different parameters were found, with most being below 0.6 (Fig. S4), thus giving a degree of confidence in the estimated values. While we did not obtain statistical confidence intervals, kernel density estimations (Figs. S5–S12) suggest differing degrees of likelihood for different parameters. Activation energies in particular showed narrow ranges of optimal values with a strong dependence on model structure.

Since optimizing all parameters against our data resulted in an R² of 0.87, it was surprising to obtain an R² of 0.99 during model validation. We note that few studies were found with data on moisture and temperature interactions under controlled conditions, and this was the only validation attempt carried out. This very high R² is partially thanks to the recalibration of initial pool sizes and may have to do with the reduced amount of data coming from a simpler experimental design compared to our study. There were only 20 data points in the validation data, one for each temperature and moisture combination. In contrast, we had 3 replicates, 11 moisture levels, and 2 temperature cycles, and therefore more data and associated variability. Despite these points and this being an initial validation step, such a close agreement using independent data and a soil that differed considerably in C content provides strong support for the model structure we used.

Model steady state or equilibrium is attained when the rate of change of all state variables equals zero, reflecting the state towards which the system will tend under invariant input and forcing conditions. Even though external drivers are in constant change in natural systems, steady-state information can indicate the approximate model behaviour under specific average conditions. Results here showed that the model M2-dif gives realistic values in the range of temperature for which it was calibrated but leads to unrealistic values under colder conditions. In addition, the C_P pool shows little sensitivity to changes in C input. While the model fitted well the validation data, it may not be suitable when applied outside the conditions used for development and may need further changes for field applications. The limitations encountered are characteristic of non-linear microbial models and mark their current limitations as predictive tools. However, such limitations are most likely the result of missing processes that still need to be adequately represented. For example, recent work has shown that a density-dependent mortality rate of the microbial pool can lead to much more realistic long-term simulations (Georgiou et al., 2017).

It is important to point out that our approach was to use a simple model with few processes and C pools and modify only those components we tested. This allowed us to distinguish the effects of each modification and minimize parameter identifiability problems arising from having too many parameters with effects that may correlate. While this allowed us to focus on specific processes, it also meant that important mechanisms were left out. Some of these mechanisms are oxygen limitations in saturated conditions, leaching of C_D, the coupling of the C and N cycles (introducing SOC quality and microbial stoichiometry limitations), and organo-mineral interactions. Our model thus needs further development to extend its application and general predictive capacity. In its current form, it cannot be extended to litter decomposition (Cotrufo et al., 2015) or organic soils, which will be much more dependent on substrate quality and less affected by carbon diffusion (Manzoni et al., 2012b). Also, peatlands and other saturated soils (Clymo, 1984; Frolking et al., 2001) will show different dynamics, reflecting the critical role of oxygen as a limiting factor. We did not include mineral adsorption of carbon as an active mechanism in this study. This is contrary to recent studies that used adsorption–desorption fluxes to explain the variability in temperature responses (Tang and Riley, 2014). However, some values of mineral desorption rates found in the literature (Ahrens et al., 2015) suggest that these rates, although important in the long term, are too slow to have a noticeable impact on the timescale of this or similar experiments, and thus on most estimates of soil respiration temperature sensitivities. Finally, nitrogen requirements will impose limits on the growth of microbial communities, which in models with microbial driven uptake and/or decomposition will also regulate C fluxes (Grant et al., 1993; Manzoni et al., 2012a). Despite such limitations, we demonstrated the effects and relevance of combining Michaelis–Menten kinetics with diffusion in mineral soils, with model results being well supported by the data.

7.1 Temperature effects

Unlike other calibrated parameters, the activation energy values for microbial (E_{a_m}) and enzyme (E_{a_e}) decay were fixed at 10 kJ, representing a positive but low temperature sensitivity. This value was used in order to be consistent with two main observations:

• a.

  The effect of E_{a_m} on the amount of microbial carbon. A high E_{a_m} results in large changes of microbial biomass C with temperature. However, observations often show a negative but moderate effect of temperature on microbial biomass (Grisi et al., 1998; Salazar-Villegas et al., 2016).

• b.

  The effect of E_{a_e} on carbon decomposition rates. High E_{a_e} values result in increasing accumulations of soil C with warming (Allison et al., 2010; Tang and Riley, 2014) as a consequence of a decrease in the enzyme pool caused by accelerated turnover. This is a critical aspect of enzyme-driven soil carbon models and largely determines simulated responses to long-term warming. Experimental evidence for E_{a_e} is lacking, but the latest observations of mid-term responses to warming are compatible with low values (Crowther et al., 2016).

The optimized E_{a_V} value of models with first- and second-order decomposition kinetics were in the range 40–50 kJ, translating to a Q₁₀ of ca. 2. In contrast, for all but one model using M decomposition, values were above 90 kJ, translating to a Q₁₀ of nearly 4. This high value was apparent in the modelled respiration fluxes only at lower temperatures, while at temperatures higher than 20 the apparent Q₁₀ approximated the more commonly observed value of 2. Such results followed closely our observations and agree well with general trends in Q₁₀ along the temperature axis reported by Hamdi et al. (2013). These values were mostly stable at high levels of soil moisture but increased sharply under drier conditions. This moisture relationship, however, is not necessarily the norm and seems to depend on initial conditions and/or pool dynamics, as demonstrated by the validation step (Fig. 9), where the apparent E_a remained close to 90 kJ and thus near the parameterized value. Also the change in E_a with moisture content followed a different trend in the validation data, although again values increased with lower moisture.

The difference between prescribed and observed temperature sensitivities may be related to two factors. First, the apparent sensitivities do not represent the instantaneous sensitivities dictated by the prescribed values but reflect also the effects of other limiting factors that change with time. Pool sizes, including C_M and C_E, may differ from the initial conditions as time progresses, making measurements at different temperatures not strictly comparable. The observation that Q₁₀ values from studies using short incubation times (hours to days) are higher compared to those using longer incubation times (Hamdi et al., 2013) is consistent with this idea. The second factor is related to the temperature sensitivity of the K constant of Michaelis–Menten kinetics. Our results are well in line with the theory discussed by Davidson and Janssens (2006), who stated that “because the K_m of most enzymes increases with temperature, the temperature sensitivities of K_m and V_max can neutralize each other, creating very low apparent Q₁₀ values”. Indeed, this seems to be the most important effect of introducing Michaelis–Menten kinetics in our simulations – not, as first assumed, the effects of concentrations of either the C_P or $C_{{\mathrm{E}}_{\mathrm{D}}}$ pools, since the choice of using M vs. M_r kinetics had only a small impact on the results.

The model results described above are thus emergent effects leading to apparent temperature sensitivities that vary in time but are based on constant model parameter (E_a) values. These results demonstrate how apparent sensitivities are the result of the offsetting effects of different processes (e.g. sensitivities of Michaelis–Menten parameters V vs. K_m) and how different values can be measured when soil pool dynamics change (e.g. through changes in diffusion limitations) even when the underlying temperature sensitivities are the same. Much of the variability in reported temperature sensitivities of soil respiration, and in particular its relationship with soil moisture (Craine and Gelderman, 2011), may be the result of the changing dynamics in microbial, enzyme, and dissolved C pools during measurement times. Clearly, misleading conclusions regarding an intrinsic temperature sensitivity of soil C decomposition are often reached by the usual procedure of fitting a simple function to respiration vs. temperature data.

Decomposition, which was only modelled, consistently showed a lower apparent temperature sensitivity than respiration, with a Q₁₀ between 1 and 2 for our experiment and just below 3 for the validation study. These values may be the most relevant for predicting long-term changes, since uptake and respiration ultimately depend on C made available by decomposition. These rather low sensitivities are consistent with some integrative studies at the ecosystem level (Mahecha et al., 2010) and again likely respond to the temperature sensitivities of K_m and V_max neutralizing each other. Such results raise the question of what E_a or Q₁₀ values – i.e. those apparent for respiration, those apparent for decomposition, or those parameterized – are best suited for conventional first-order empirical soil models. Since these models will tend to have similar apparent and intrinsic behaviour, the answer is not clear and will require further research. Ultimately, the best option may be to abandon such models and develop better validated mechanistic alternatives for prediction purposes.

7.2 Moisture effects and diffusion limitations

Diffusion fluxes are a function of water content, diffusivity coefficients, and pool concentrations. Different equations have been used to calculate diffusion as a function of water content in soils (Hamamoto et al., 2010; Hu and Wang, 2003). All these equations generally predict a strong positive near-exponential effect of water content on diffusion. Following previous studies (Manzoni et al., 2016), we chose the function from Hamamoto et al. (2010). This equation allows for an adjustment of the percolation threshold (θ_th) in different soils. We note that when using the θ_th obtained during calibration (0.063) we also obtained a high fit to the validation data (R²=0.97, data not shown), so the recalibration of θ_th led to a noticeable but small improvement. While the value 0.063 for our soil came close to the water potential of −15 MPa found in previous studies (Manzoni and Katul, 2014), this relationship did not hold for the soil used for validation, where we assumed a higher clay and silt content from its classification. Thus, a prerequisite for applying our model to other soils is finding a relationship between θ_th and soil type that holds in all cases.

Diffusion regulations can be implemented either by simulating two separate pools between which diffusion takes place or by determining the available amount of a pool as a function of diffusivity (or conductance in our case) at each time step. In our model we used a combination, simulating a diffusion flux between enzyme pools and calculating how much C_D is available for uptake at each time step. We did not assume a diffusion regulation of available polymeric C, an approach that is closer to empirical functions scaling the decomposition flux directly and that has been implemented in other microbial models (Davidson et al., 2012).

In our study especially, but also in the validation data, the moisture response tended to become less linear and have a larger plateau at higher temperatures. The mechanisms leading to such interactions are still unclear, but our model comparison indicates that solute diffusion limitations play a central role. The plateau behaviour, a decrease near saturation, and even near-linear responses all contrast with the near-exponential relationship between moisture and conductance given by Eq. (11) and with the fact that no oxygen limitations at high moisture levels were modelled. They may, however, result from a faster depletion of available carbon at high moisture levels and at high temperatures, driving down the accumulated fluxes over time.

While a low supply of O₂ usually limits respiration rates in saturated soils under field conditions, O₂ seemed to have a negligible effect in our study. At 35 ^∘C, where fluxes were highest, no clear drop in respiration was observed near saturation, as is expected when O₂ becomes limiting. Rather, the general behaviour was well simulated by our models using solute diffusion limitations only. Schurgers et al. (2006) found that the anaerobic fraction in soils with air O₂ concentrations over 10 % is low until very close to saturation. The minimum flask air O₂ concentrations (corresponding to 56 000 ppm of CO₂, the maximum accumulated by a sample before changing headspace air) was over 15 % O₂, which next to the small sample sizes would not indicate an O₂ limitation.

In models where decomposition and respiration are separated processes, these fluxes can show different responses. This decoupling is especially evident when diffusion limitations come into play. Plots of modelled fluxes against temperature and moisture (Fig. S3) showed a different relationship when comparing respiration and decomposition. Figure 10 shows modelled decomposition against respiration (using M2-dif) as accumulated values, each line being a sample at a different water content. Without any diffusion limitation, the relationship follows a slope of ca. 0.3, determined by 1−f_ug, where f_ug is the fraction of uptake to growth (the C use efficiency). This slope, however, changes as diffusion becomes limiting, with temperature also playing a role as evidenced by the shifts in the slope occurring at various intervals. With time these fluxes will tend to equilibrate as the C_D and $C_{{\mathrm{E}}_{\mathrm{D}}}$ pools adjust. But the proportionality between these fluxes is not constant and will depend on moisture, temperature, and time, even after months of incubation. These results show that, without a proper modelling framework and when assuming a constant proportionality, interpretations based only on respiration activity may lead to wrong conclusions about the dynamics of organic matter decomposition, especially at low moisture content levels and in short- and mid-term experiments.

Figure 10 Modelled decomposed vs. respired C shown as accumulated values over the entire simulated incubation, including temperature steps. Each line is a sample at a different moisture content level.

8 Conclusions

As the main mechanism linking water content with the movement of substrates, microbes, and enzymes, diffusion plays a central role in soil organic matter decomposition. We here showed that integrating it into models can significantly improve our understanding of soil C dynamics. Diffusion-based models were better at simulating the effects of moisture and improved the simulated temperature responses, thus allowing for a better interpretation of the observed temperature sensitivities. This and similar studies indicate that measured temperature sensitivities cannot be generalized or correctly interpreted without having a full understanding of the relevant mechanisms, their interactions, and the state of soil carbon and microbial pools.

We also found evidence that Michaelis–Menten kinetics plays an important role in soil C dynamics, explaining the strong differences in temperature sensitivities across temperature ranges. Our results are consistent with relatively high activation energies for both the V and K Michaelis–Menten parameters and generally lower apparent values.

Creating models that capture the variability in the response of C dynamics across different soils and at different levels of driving factors remains challenging. However, process-based models are of central importance for establishing confidence in C cycle predictions and soil–climate feedbacks. As seen here, the structure and process representation of models can be critical for simulating the complex response of soil C fluxes to combined changes in temperature and moisture. Diffusion as a moisture regulation of soil C fluxes has not been used in large scale predictions, which still rely on empirical scaling functions. Evidence of interactions seen in experiments and presented from a mechanistic model perspective indicates that these simpler approaches do not always hold. Further research should focus on more extensive validation and finding the relationships necessary for extending the application of models to diverse soil types.