Linear two-pool models are insufficient to infer soil organic matter decomposition temperature sensitivity from incubations

Terrestrial carbon (C)-climate feedbacks depend strongly on how soil organic matter (SOM) decomposition responds to temperature. This dependency is often represented in land models by the parameter Q10, which quantifies the relative increase of microbial soil respiration per 10 °C temperature increase. Many studies have conducted paired laboratory soil incubations and inferred “active” and “slow” pool Q10 values by fitting linear two-pool models to measured respiration time series. Using a recently published incubation study (Qin et al. in Sci Adv 5(7):eaau1218, 2019) as an example, here we first show that the very high parametric equifinality of the linear two-pool models may render such incubation-based Q10 estimates unreliable. In particular, we show that, accompanied by the uncertain initial active pool size, the slow pool Q10 can span a very wide range, including values as high as 100, although all parameter combinations are producing almost equally good model fit with respect to the observations. This result is robust whether or not interactions between the active and slow pools are considered (typically these interactions are not considered when interpreting incubation data, but are part of the predictive soil carbon models). This very large parametric equifinality in the context of interpreting incubation data is consistent with the poor temporal extrapolation capability of linear multi-pool models identified in recent studies. Next, using a microbe-explicit SOM model (RESOM), we show that the inferred two pools and their associated parameters (e.g., Q10) could be artificial constructs and are therefore unreliable concepts for integration into predictive models. We finally discuss uncertainties in applying linear two-pool (or more generally multiple-pool) models to estimate SOM decomposition parameters such as temperature sensitivities from laboratory incubations. We also propose new observations and model structures that could enable better process understanding and more robust predictive capabilities of soil carbon dynamics.


Introduction
Predicting the fate of the large amount of soil carbon stored globally [more than twice as in the current atmosphere (Ciais et al. 2013)] is critical for Abstract Terrestrial carbon (C)-climate feedbacks depend strongly on how soil organic matter (SOM) decomposition responds to temperature. This dependency is often represented in land models by the parameter Q 10 , which quantifies the relative increase of microbial soil respiration per 10 °C temperature increase. Many studies have conducted paired laboratory soil incubations and inferred "active" and "slow" pool Q 10 values by fitting linear two-pool models to measured respiration time series. Using a recently published incubation study (Qin et al. in Sci Adv 5(7):eaau1218, 2019) as an example, here we first show that the very high parametric equifinality of the linear two-pool models may render such incubationbased Q 10 estimates unreliable. In particular, we show that, accompanied by the uncertain initial active pool size, the slow pool Q 10 can span a very wide range, including values as high as 100, although all parameter combinations are producing almost equally good model fit with respect to the observations. This result quantifying feedbacks between terrestrial biogeochemistry and the Earth's climate. In particular, land models need accurate quantifications of soil organic matter (SOM) decomposition sensitivity to soil temperature under a warming climate. To address this need, thousands of empirical studies have inferred SOM decomposition temperature sensitivities, often characterized by the parameter Q 10 that measures the relative increase in microbial SOM respiration per 10 °C temperature increase (e.g., Fang and Moncrieff 2001;Haddix et al. 2011;Hamdi et al. 2013;Lloyd and Taylor 1994;Schadel et al. 2013). However, debate continues regarding the laboratory-based Q 10 values for incorporation into SOM models, or whether alternative model structures are more appropriate Hemingway et al. 2019;Tang and Riley 2015).
The use of Q 10 (or alternatively activation energy in the Arrhenius representation of SOM turnover) is inherent in the substrate quality (SQ) conceptual model of SOM decomposition (e.g., Hartley and Ineson 2008;Parton et al. 1988;Wetterstedt et al. 2010). In these SQ-based models, SOM is assumed to consist of several abstract (i.e., not measurable) pools, each of which has its own decomposition turnover time that varies as a function of temperature, moisture, and other edaphic conditions. Despite their popularity, SQ-based models cannot explain the observation that old carbon (which by assumption is recalcitrant and therefore of long turnover time) can be decomposed rapidly when made accessible to active microbes (e.g., Kleber et al. 2011;Nowinski et al. 2010;Schuur et al. 2009).
Another problem with SQ-based models is their usually very high parametric equifinality, i.e., for a given set of observations, many parameter combinations are able to produce model outputs that fit the observations equally well (e.g., Beven and Freer 2001;Tang and Zhuang 2008), while their predictions often diverge wildly. For instance, using the APSIM model, which has six SQ-based SOM pools, Luo et al. (2015) showed that "convergent modeling of past soil organic carbon stocks" is followed by "divergent projections" of twenty-first century SOM dynamics. In a subsequent study, Luo et al. (2017) deduced the same conclusion for a SQ-based linear two-transfer-pool model (that considers interactions between active and slow pools; see supplemental material for a description and solutions for these two types of linear two-pool models). These findings therefore motivate an analysis of whether the high parametric equifinality of SQ-based models may lead to large uncertainties in temperature sensitivities inferred from incubation experiments.
Notwithstanding these conceptual defects, the simplicity of SQ-based models has led them to be widely applied to interpret laboratory soil incubations. For instance, Liang et al. (2015) compared results from several linear multiple-pool models to infer temperature sensitivity parameters from respiration time series measured during SOM incubations. They found that the two-discrete-pool model, which does not consider interactions between active and slow pools, was able to fit the soil respiration data very well under many conditions, and Schadel et al. (2013) drew the same conclusion from a similar analysis. In a recent study, which we analyze here, Qin et al. (2019; hereafter Q2019) used a linear two-discrete-pool model to fit their 330day incubation of topsoil (0-10 cm) and subsoil (30-50 cm) samples collected from three sites on the northeastern Tibetan Plateau. They inferred that topsoil is more sensitive than subsoil to temperature warming and concluded that SOM models can be improved by incorporating such a contrast in temperature sensitivity.
Meanwhile, to address the conceptual difficulties with SQ-based models (e.g., "old carbon" can often behave just like "active" pool carbon), new model structures have been developed to attempt explicit representations of the myriad interactions between microbes, enzymes, substrates, and various abiotic soil processes (Abramoff et al. 2019;Dwivedi et al. 2017Dwivedi et al. , 2019Riley et al. 2014;Sulman et al. 2018;Tang and Riley 2015;Wang et al. 2013;Wieder et al. 2013). These new models then predict that SOM decomposition emerges from these interactions. In general, including more processes requires more model parameters and will potentially increase a model's parametric equifinality. Nevertheless, these new models enable more direct correspondence between their simulations and empirical measurements, and provide mechanistic insights that are otherwise not available from SQ-based models. Moreover, with improved mathematical formulations (Tang and Riley 2013, 2019b, these new models may enable first-principle based parameterization that is not feasible with SQ-based models (e.g., Tang and Riley 2019a, b).
Here, by reanalyzing the Q2019 incubation data using the linear two-discrete-pool model (e.g., Liang et al. 2015;Qin et al. 2019;Schadel et al. 2013) and the linear two-transfer-pool model (Luo et al. 2017), we first show that the high parametric equifinality makes it almost impossible to obtain a robust inference of Q 10 values and other related parameters for SQ-based SOM pools. Next, using a microbe-explicit model, i.e., the REaction-network-based model of Soil Organic Matter and microbes (RESOM) (Tang and Riley 2015), we show that apparent active and slow pools can emerge from the decomposition of a single substrate with a single temperature sensitivity (i.e., activation energy). We then discuss why mechanistic SOM decomposition models like RESOM should be more extensively explored and how laboratory incubations can be better used to inform the development of these new models.
Large parametric equifinality in linear two-pool models leads to uncertain Q 10 and other parameters Q2019 conducted a long-term (330-day) laboratory incubation for topsoil (0-10 cm) and subsoil (30-50 cm) samples collected from three sites on the northeastern Tibetan Plateau. They then applied a linear two-discrete-pool model to derive the temperature sensitivities of active and slow SOM pools, and asserted that topsoil has higher temperature sensitivity than subsoil for bulk soil, active, and slow pools. We digitized their respiration time series data using the MATLAB script grabit (https:// www. mathw orks. com/ matla bcent ral/ filee xchan ge/ 7173grabit), and conducted Markov chain Monte Carlo (MCMC) based inversions (Vrugt 2016) for a linear two-discrete-pool model using these respiration data. We note that Q2019 set their Q 10 parameter prior ranges to [1, 2] and [2, 4] for the active and slow pool, respectively (reported in Table S4 of their Supplemental Material), thereby strongly limiting the range of their possible posterior values. Analogously Schadel et al. (2013) imposed upper and lower limits to the decay rates of active and slow pools. However, some observationally-based studies have reported bulk soil Q 10 values as high as a few hundred (e.g., Hamdi et al. 2013). Therefore, if we interpret their approach correctly, the Q2019 parameter inversion results could have been biased methodologically, because these unjustified prior constraints force their inferred Q 10 values to fall within an unrealistically small range. We further note that emergent Q 10 values have been observed to vary across a wide range and be hysteretic with respect to temperature over time (e.g., Pingintha et al. 2010), and microbe-explicit models are better at resolving such hysteresis by accounting for hysteresis in microbial biomass and kinetics (e.g., Tang and Riley 2015).
We performed seven MCMC inversions to quantify how sensitive the inferred Q 10 values are to the imposed parameter ranges for both active and slow pools. For both pools, the seven inversions used maximum Q 10 values of 4, 10, 20, 40, 50, 70, and 100; the experiments are named accordingly as IQ4, IQ10, IQ20, IQ40, IQ50, IQ70, and IQ100, where IQx refers to an inversion with Q 10 bounded at a maximum value of x. Additionally, since Q 10 has been reported to be as low as 0.5 (Hamdi et al. 2013), we allowed Q 10 prior values to extend to 0, but none of the MCMC inversions resulted in Q 10 values smaller than 1.
When the maximum Q 10 is bounded at 4 (for both pools), our MCMC inversions ( Fig. 1) corroborate Q2019s inference that Q 10 values of both the active and slow pools in topsoil ( 2.06 ± 0.012 and 3.89 ± 0.092 , respectively) are higher than those in subsoil ( 1.51 ± 0.031 and 2.86 ± 0.23 , respectively), although these Q 10 values differ from those in Q2019, who reported mean values of about 1.8 and 2.4 for topsoil active and slow pools, respectively; and 1.6 and 2.2 for subsoil active and slow pools, respectively. These differences may exist because (1) we did not exactly recover their data and the inversion is sensitive to data uncertainty characterization, and (2) we used different, but still restrictive, parameter ranges. Further, when we imposed the same upper prior ranges as Q2019 did, the inversion inferred topsoil Q 10 values for the active and slow pool of 2.00 ± 0.003 and 3.95 ± 0.045 , respectively, and for the subsoil, of 1.53 ± 0.030 and 3.69 ± 0.26 , respectively. However, for inversions with larger maximum Q 10 values, we find very wide ranges of slow pool Q 10 values (i.e., a manifestation of high parametric equifinality; Fig. 1c, d), while all predicted respiration time series match the observations equally well (Fig. 1a, b). Moreover, this large parametric equifinality is higher for topsoil than for subsoil. For example, the most uncertain inference occurs for topsoil when Q 10 is bounded at 100, where the posterior maximum and minimum slow pool Q 10 values are 97.5 and 4.6, respectively; whereas for subsoil the most uncertain inference occurs when Q 10 is bounded at 70, where the posterior maximum and minimum slow pool Q 10 values are 40.9 and 2.80, respectively. In summary, our results are consistent with those of Luo et al. (2017), even though they instead used time series of soil carbon stocks as constraints in their parameter inversions for the linear two-transfer-pool model, which has one more parameter than does the linear two-discrete-pool model (see supplemental material).
When the posterior distributions of slow pool decay rates at different temperatures and initial active pool fractions are compared for inversions with Q 10 values bounded at 4 and 100 (i.e., IQ4 and IQ100 in Fig. 2), we find that the wider range of topsoil slow pool Q 10 for IQ100 (than for IQ4; Fig. 2a, c, e) occurs because (1) the inferred slow pool decay rates at 10 °C and 20 °C are left edge-hitting and more center-peaked, respectively, and (2) the initial active pool is larger in IQ100 than in IQ4 (Fig. 2g, h). For subsoil, the posterior differences between IQ100 and IQ4 are smaller, but high parametric equifinality remains (Fig. 2b, d, f).
Our analyses above show that more than 12,500 combinations (out of the 500,000 MCMC samples for each inversion) of posterior parameters fit the incubation soil respiration observations equally well. This high parametric equifinality and the resultant wide range of Q 10 values indicate that (1) it is difficult to use a linear two-discrete-pool model to robustly derive the temperature sensitivity of SOM decomposition from laboratory incubation soil respiration data alone, and (2) the initial active and slow pool SOM stocks are difficult to robustly quantify. For assertion (2), although the differences in the fraction of initial SOM as active pool between inversions (Fig. 2g, h) are relatively small (about 0.09 vs. 0.11 on average), they force significant uncertainty into the slow pool decay rates, making those small differences consequential (which is reflected in the wide range of slow pool Q 10 ).
Parameters inferred from linear two-discrete-pool models are often applied to inform more complex transfer-pool models (e.g., CENTURY; Parton et al.  1988)). However, the difference in these two model structures makes it difficult to use the discrete-pool model to demonstrate the impacts of the high parametric equifinality on model projections. To facilitate such a demonstration, we conducted an inversion using a linear two-transfer-pool model (see supplemental material for model formulation) for the topsoil with maximum Q 10 values of 4 and 100, and performed forward simulations using the inferred posterior parameters. We found that the full set of inversions using the linear two-transfer-pool model produced very similar results (not shown) to those shown in Figs. 1 and 2. We find the MCMC inversions with the linear twotransfer-pool model are equally uncertain (Fig. 3) as those for the linear two-discrete-pool model (Figs. 1,  2). Moreover, the associated 25-year projection with no warming demonstrates a very strong divergence through amplifying effects from the carbon input flux on the posterior parameters after year 2 (Fig. 4). If warming is considered, we expect that the divergence shown in Fig. 4 will become even larger because it will include additional uncertainty from Q 10 (Feller 1968). Given these results and previous studies Tang and Zhuang 2008), we conclude that linear two-pool models (and probably multiple-pool models) are of such high parametric equifinality that they cannot be well-constrained by incubation time series observations of either respiration or soil carbon stocks.
We also investigated whether additional data could help better constrain the linear two-pool models. We tested this idea by dramatically increasing the sampling frequency and reducing the observational error of respired CO 2 . This synthetic inversion experiment with the linear two-transfer-pool model indicates that, although sampling respiration daily at 2% relative error (which is rarely achievable in reality) does constrain Q 10 values quite well, it is still insufficient to robustly infer other model parameters (i.e., initial active pool fraction and fraction of decomposed active pool into slow pool) that also must be used for predictions (Fig. 5). Increasing the relative error of daily observations to a more realistic yet still challenging-to-achieve 5% would have resulted in an even more biased posterior of Fig. 2 For topsoil, when the MCMC inversion applies a maximum Q 10 value of 100 (i.e., IQ100), the slow pool Q 10 distribution (panel a) has a longer tail (than IQ4) because the distribution of IQ100 slow pool decay rate ( k s ) at 10 °C (panel c) does not align well with the distribution at 20 °C (panel e) and the initial active pool is larger (panel g). Differences between IQ100 and IQ4 for subsoil also exist but are comparatively smaller (panels b, d, f, h) the initial fraction of SOM as active pool (0.027 for IQ4, and 0.21 for IQ100), and transfer coefficient (0.02 for IQ4, and 0.83 for IQ100) between active and slow pools during decomposition ( Figure S1). Analogously, we suspect that adding multi-site bulk SOM stocks or radiocarbon respiration data as observational constraints is unlikely to improve the inversion, because of the weak signal of the slow pool as can be triggered in an incubation experiment (which will lead to biased estimation of the fraction of decomposed active pool that enters slow pool, as shown in Fig. 5g).  We next conducted incubation experiments with the microbe-explicit RESOM model (Abramoff et al. 2019;Tang and Riley 2015) to show that apparent active and slow pools may be inferred from incubation data, even though a single substrate with a single activation energy of 55 kJ (mol C) −1 was actually undergoing microbial decomposition. For this demonstration of commonly-observed SOM decomposition patterns, we configured RESOM with one microbial population (with a substrate assimilation efficiency of 0.27 for topsoil, and 0.50 for subsoil, respectively), one extracellular enzyme, one mineral soil surface (with surface area of 575 g C m −2 and 1800 g C m −2 for topsoil and subsoil, respectively), and one substrate (which can exist in both aqueous and solid phases). The substrate assimilation efficiency and mineral soil surface were obtained through trial-anderror by comparing the simulated respiration data to that of Q2019, and other parameters are from previous publications (Abramoff et al. 2019;Tang and Riley 2015). Following Tang and Riley (2015), we first generated initial conditions from a spinup to equilibrium and then ran RESOM for 330 days at the two temperatures to mimic the incubation. After normalizing the respiration to initial SOM, these RESOM simulations reproduced the Q2019 respiration time series for both topsoil and subsoil within the observational uncertainty (Fig. 6a, b), suggesting that the predictive capability of RESOM is scientifically reasonable to demonstrate the potential to mis-identify multiple substrate pools from incubation data. We also note that when the temperature forcing has higher temporal variability, like that observed in field conditions, the emergent temperature sensitivity is dynamic and the corresponding instantaneous Q 10 can vary across a wide range. This dynamics occurs because temperature sensitivity of SOM decomposition depends on all the interactions between microbes, substrates, enzymes, and mineral surfaces represented in the model, as discussed extensively in Tang and Riley (2015).
When the natural logarithm of the RESOM-predicted respiration is plotted against time (Fig. 6c, d), c, d, e, f and g are probability distributions of active pool Q 10 , slow pool Q 10 , initial active pool fraction, active pool k f at 10 °C, slow pool k s at 10 °C, and the fraction of decomposed active pool into slow pool nonlinear curves are found for both topsoil and subsoil at 10 °C and 20 °C. This commonly-observed pattern, which can also be produced with other model parameters or by other microbial models (e.g., Grant et al. 1993;Smith 1979;Wieder et al. 2014), is a primary reason that researchers assume SOM has multiple pools with different turnover times. In this one-substrate-pool RESOM version, these nonlinearcurves occur because, as the incubation progresses, (1) substrate is gradually depleted, and (2) microbial biomass is reduced accordingly. Indeed, by conducting linear two-discrete-pool model parameter MCMC inversions with the RESOM predicted respiration time series data, the inferred model is found to fit almost perfectly to the RESOM subsoil simulation (Fig. 6d), and relatively accurately for the topsoil (more than two pools are required to fit the RESOM simulation for topsoil) (Fig. 6c). This contrast of model fitting to RESOM simulations for topsoil and subsoil suggests that definitions of SQ pools are not unique even for a single substrate. In real soils, substrate diversity is likely very high, and the forcing temperature is much more variable. Thus, depending on the relative strength of interactions (and likely the imposed temperature forcing in incubation experiments), SQ-based linear multiple-pool models may succeed in one situation (still at the expense of uncertain parameters) while fail in another, even for similar overall system properties. Improving new SOM model structures from better use of observations and cross-disciplinary knowledge Our analyses demonstrate that, because of their high parametric equifinality, SQ-based linear two-pool models provide inaccurate inferences of SOM decomposition temperature sensitivity (and other related parameters). We also showed that the oft-inferred active and slow pools could be artificial constructs that emerge even when only a single substrate is being decomposed. Therefore, improving SOM decomposition models and temperature sensitivity inferences from observations requires more explicit representation of the interactions between the many biotic and abiotic processes that contribute to SOM dynamics. Intuitively, one could argue that models with more mechanistic representations may require more parameters and may not resolve the curse of high parametric a b c d Fig. 6 The one-substrate RESOM model accurately fits the laboratry incubation data (panels a and b). The model predictions and observations show why more than one apparent SOM pool could be mistakenly inferred (panels c and d). In particu-lar, calibrating (from MCMC inversion) the two-discrete-pool model fits the RESOM simulations for subsoil almost perfectly (panel d) and reasonably well for topsoil (panel c). Panels a and c are for topsoil, b and d are for subsoil equifinality. However, more mechanistic SOM models have at least several advantages, including that they: (1) can incorporate first-principle based representations of kinetic parameters, e.g., substrate affinity parameter, maximum substrate uptake rate, or even substrate assimilation efficiency, whose variability is unlikely to be resolved through multiplier functions as often used in SQ-based models (Tang and Riley 2019a, b); (2) can provide mechanistic insights into important systems dynamics, including why some old carbon (i.e., with relatively negative Δ 14 C values) can relatively easily be decomposed when exposed to microbes and how mineral protection affects SOM persistence (Hemingway et al. 2019;Tang and Riley 2015); (3) can produce more realistic dynamic predictions of how temperature and moisture affect SOM dynamics, which have been statically represented using multiplier functions in SQ-based models (Tang and Riley 2019b), allowing for observed bulk soil Q 10 values to vary over a wide range or be hysteretic (Pingintha et al. 2010); (4) enable consistent representation of the myriad biogeochemical processes that control the cycling of various important chemicals (e.g., CO 2 , CH 4 , N 2 O, etc.), a need evidenced by the SQ-based models' ever-increasing introduction of new response functions for new processes (Sierra et al. 2015;Zhuang et al. 2004); (5) enable a direct link between models and the vast amount of knowledge accumulated in related fields like microbiology, microbial ecology, aqueous chemistry, physical chemistry, etc., so that more diverse data (e.g., microbial cell size) can be used for model constraints [e.g., Jin and Kirk (2018a, b), la Cecilia et al. (2019), LaRowe and Van Cappellen (2011), and Tang and Riley (2019a, b)]; (6) may allow the use of parametric equifinality to reduce predictive divergence. This counter-intuitive idea is analogous to the approach used to link kinetic gas theory and Navier-Stokes equations in fluid mechanics (Batchelor 1967), except that in microbial models parametric equifinality could be represented with bio-diversity. When biodiversity is further integrated with known ecolog-ical principles (e.g., competition and symbiosis), the many parameters can compensate each other to reduce predictive uncertainty and meanwhile enhance the models' resilience to perturbations (see Sect. 6 of Tang and Riley (2017) for a discussion of this concept).
Finally, with the advantages described above, we emphasize that mechanistic models will benefit from laboratory respiration measurements when biogeochemical and biogeophysical processes affecting respiration are also measured. Although we demonstrated that the inferred temperature sensitivities in Q2019 are likely biased because of their assumed SOM model structure and inversion approach, their other measurements characterizing soil aggregates, mineral soil protection, and relative microbial abundances and composition are valuable resources for future analyses. As experiments begin collecting such critical information, a suite of process-rich predictive models can be developed which in turn enable better interpretation of laboratory incubation studies and improve predictions of terrestrial biogeochemical interactions at larger scales.

Conclusions
SQ-based SOM models have been argued to be selfconsistent (Blankinship et al. 2018) and have provided useful understanding of SOM dynamics, yet we show here that their high parametric equifinality prevents them from reliably being used to infer temperature sensitivity and other related parameters of assumed SOM pools from laboratory incubation data. This problem is serious, since many global land models used in climate assessments have SOM decomposition models parameterized from these types of observations. New and more mechanistic models of SOM decomposition should be explored to represent the myriad biotic and abiotic entities and interactions that affect temperature [and moisture; e.g., Tang and Riley (2019b)] sensitivity of SOM decomposition. By combining observations and model structures that more explicitly represent the underlying processes affecting SOM decomposition, we can expect improved understanding and predictions of how SOM cycling may respond to climate change.