Field-warmed soil carbon changes imply high 21st-century modeling uncertainty

The feedback between planetary warming and soil carbon loss has been the focus of considerable scientific attention in recent decades, due to its potential to accelerate anthropogenic climate change. The soil carbon temperature sensitivity is traditionally estimated from short-term respiration measurements – either from laboratory incubations that are artificially manipulated or from field measurements that cannot distinguish between plant and microbial respiration. To address these limitations of previous approaches, we developed a new method to estimate soil temperature sensitivity (Q₁₀) of soil carbon directly from warming-induced changes in soil carbon stocks measured in 36 field experiments across the world. Variations in warming magnitude and control organic carbon percentage explained much of field-warmed organic carbon percentage (R² = 0.96), revealing Q₁₀ across sites of 2.2 [1.6, 2.7] 95 % confidence interval (CI). When these field-derived Q₁₀ values were extrapolated over the 21st century using a post hoc correction of 20 Coupled Model Intercomparison Project Phase 5 (CMIP5) Earth system model outputs, the multi-model mean soil carbon stock changes shifted from the previous value of 88 ± 153 Pg carbon (weighted mean ± 1 SD) to 19 ± 155 Pg carbon with a Q₁₀-driven 95 % CI of 248 ± 191 to −95 ± 209 Pg carbon. On average, incorporating the field-derived Q₁₀ values into Earth system model simulations led to reductions in the projected amount of carbon sequestered in the soil over the 21st century. However, the considerable parameter uncertainty led to extremely high variability in soil carbon stock projections within each model; intra-model uncertainty driven by the field-derived Q₁₀ was as great as that between model variation. This study demonstrates that data integration should capture the variation of the system, as well as mean trends.

1 Introduction

The flux of carbon dioxide between the soil and atmosphere is a major control on atmospheric carbon dioxide concentrations. Warming temperatures, driven by increases in atmospheric carbon dioxide, have the potential to stimulate carbon decomposition, accelerating its release into the atmosphere (Davidson and Janssens, 2006). If this is not counterbalanced by an equivalent increase in primary productivity (the opposing carbon flux), then it has the potential to drive a land carbon–climate feedback that will accelerate anthropogenic climate change. Recent global compilations of data from ecosystem warming experiments lend support to this idea (Carey et al., 2016), suggesting that warming alone could drive a loss of carbon from the upper soil horizons (Crowther et al., 2016). However, these studies addressed the impact of warming in isolation, and it remains unclear how this process will interact with the variety of other global change drivers to affect the global soil carbon stock over the rest of this century. Reflective of such uncertainty, soil carbon changes projected for 2100 under the Representative Concentration Pathway (RCP) 8.5, business-as-usual, scenario for the Coupled Model Intercomparison Project Phase 5 (CMIP5) range from −70 to 250 Pg carbon across different Earth system models (ESMs; Todd-Brown et al., 2014), making the land–carbon feedback one of the largest sources of uncertainty in future climate projections (Friedlingstein et al., 2014). Improving the soil carbon component of the Earth system models is essential to predicting the future evolution of the Earth system and thus establishing meaningful greenhouse gas emissions targets.

A fundamental parameter describing soil temperature sensitivity in soil carbon models is the Q₁₀ – the factor of the change in decomposition rate associated with 10 ^∘C of warming from a reference temperature (Davidson and Janssens, 2006; Lloyd and Taylor, 1994). Traditional laboratory incubations have found a wide range of Q₁₀ values, varying from 1.4 (Townsend et al., 1997) to > 3 (Davidson et al., 1998, 2006), with 2 being the most commonly accepted value. Complicating this, theoretical analyses based on chemical kinetics suggest Q₁₀ is itself dependent on temperature (Davidson and Janssens, 2006; Lloyd and Taylor, 1994), though these values are typically very close to 2 in most environmental temperature ranges (Lloyd and Taylor, 1994). More recently, large-scale analyses of field respiration have converged on Q₁₀ estimates of 1.4 to 1.5 (Bond-Lamberty and Thomson, 2010; Hashimoto et al., 2015; Mahecha et al., 2010). Unsurprisingly, this temperature response is also critical in Earth system models, where the temperature sensitivity parameter is known to be a major driver of variation (Booth et al., 2012; Jones and Cox, 2001; Jones et al., 2006). However, it is unclear what is driving the lower Q₁₀ estimates in these field-based syntheses compared to the average lab-based estimates from single-site studies, and there appears to be a relatively wide range of “typical” Q₁₀ values in the literature. Nevertheless, most Earth system models use values that range from 1.5 (Oleson et al., 2013; Raddatz et al., 2007) to 2 (Bonan, 1996; Cox, 2001).

Traditionally, these Q₁₀ values have been calculated from warming-induced changes in soil respiration rates. However, this approach has two main limitations: (1) respiration rates measured under idealized laboratory conditions fail to reflect the structure, heterogeneity, and variability of natural systems, and (2) field measurements cannot directly isolate heterotrophic soil respiration from autotrophic root respiration without substantially altering the system. Estimating Q₁₀ directly from warming-induced changes in soil carbon stocks could be a valuable approach to addressing these limitations, but the variability and relative imprecision of soil carbon stock data necessitate a large sample size to adequately describe variation at the global scale (Bradford et al., 2016). Yet results from a recent Earth system model meta-analysis indirectly suggest that, with enough sample coverage, it may be possible to infer Q₁₀ directly from changes in soil carbon stocks (Todd-Brown et al., 2014).

Here we present a new approach to estimating the global Q₁₀ value from net changes in soil carbon stocks under warming, rather than soil respiration measurements, and examine the consequences of these estimates – with associated uncertainty – on CMIP5 Earth system model projections of global carbon storage over the rest of the 21st century. To do this, we use a global database of soil carbon stock data from 36 field-warming experiments around the world, each of which includes control (ambient) plots, and those which have been warmed for extended (years to decades) periods of time (Crowther et al., 2016) (Table S1 in the Supplement), and outputs from 20 Earth system models in the CMIP5 RCP 8.5, business-as-usual, scenario (Taylor et al., 2011; see Table 1 and S3 for model and output details). These field data were used previously to derive Earth system model independent estimates of global soil carbon temperature sensitivity where the effect of warming was isolated from other global change drivers or the interacting climate system (Crowther et al., 2016). In this study we develop a novel approach that enables us to explore these field results in the context of the temperature sensitivity function (Q₁₀) used in an integrated Earth system model. We then examine the consequences of the data-driven Q₁₀ estimates, and the associated uncertainty, for CMIP5 Earth system model projections of global carbon storage over the rest of the 21st century using a novel post hoc modification of the CMIP5 simulation outputs.

Table 1 This is a summary of the soil decomposition sub-models for the ESMs used in this study and includes the number of pools, structure of the carbon exchange between those pools, temperature sensitivity function, and citations. Temperature sensitivity function is denoted as f(T), borrowed from the Century model (Parton et al., 1987, 1988); Lloyd and Taylor, taken from the recommended form from Lloyd and Taylor (1994); Q₁₀(T), a temperature-dependent Q₁₀ as defined by Arora (2003); or a Q₁₀ parameter for the Q₁₀ function as defined in this paper.

2 Materials and methods

2.1 Field sites

The field sites were drawn from a previous analysis (Crowther et al., 2016). From this initial database of 48 paired case–control studies, we selected 36 studies that were run longer than 2 years to match the metastable state assumption articulated below. Eighteen of these sites were temperate grasslands, savannas, and shrublands; 10 were temperate broadleaf and mixed forests; six were tundra; one boreal forests or taiga; and one site was in a Mediterranean forest, woodland, and scrub. A traditional statistical analysis of the sites is provided by Crowther et al. (2016). For this study, we used the increase in soil temperature due to warming, length of the study, and the percent of soil organic carbon in paired warmed and control plots (Table S1).

2.2 Q₁₀ calculations

We calculated traditional Q₁₀ estimates based on these warming-induced soil carbon losses, enabling us to embed this temperature sensitivity information into a soil decomposition model framework. Traditional soil decomposition models follow a first-order linear decay framework:

$$\begin{array}{}\text{(1)}& {\displaystyle \frac{\mathrm{d}\mathit{C}\left(t\right)}{\mathrm{d}t}}=u_{\mathrm{in}}\left(t\right)\mathit{b}-\left({\mathbf{Q}}_{\mathrm{10}}(T,t)\mathbf{KA}\right)\mathit{C}\left(t\right){\mathbf{KQ}}_{\mathrm{10}}(T,t)\mathbf{A},\end{array}$$

where C is a vector of soil carbon pools with unique turnover times, t is time, u_in is a scaler of soil carbon inputs, b is an allocation vector describing how the inputs are divided between the soil carbon pools, K is a diagonal matrix representing the decomposition rates of the pools, Q₁₀ is a diagonal matrix with entries of the form $q_i^{(T\left(t\right)-T_{\mathrm{0}})/\mathrm{10}}$ representing the temperature sensitivity factor, T is a scalar describing the soil temperature and T₀ an arbitrary reference temperature, and A is the transfer matrix representing movement of carbon between soil carbon pools.

The temperature sensitivity was assumed to be constant across pools. This allows the diagonal Q₁₀ matrix to be collapsed into a single scalar value of the form $Q_{\mathrm{10}}^{(T\left(t\right)-T_{\mathrm{0}})/\mathrm{10}}$. This constant temperature sensitivity assumption is discussed below and follows the structure of the CMIP5 Earth system models.

In general, there are three classes of pool structure for traditional models: independent, where there was no exchange between soil carbon pools, making A the identity matrix; cascade, where pools with faster turnover times passed carbon to pools with slower turnover times, making A a lower triangular matrix; and full-feedback models, where carbon was exchanged between faster and slower pools and vice versa, making A a fully dense matrix. In all cases KA is an M matrix, implying there exists an inverse with all positive entries. For the independent and cascade pools KA is diagonalizable, implying it can be broken down into a diagonal matrix D and an invertible matrix P such that KA = P⁻¹DP.

For most well-developed soils, soil carbon stocks are at a metastable state where soil inputs approximately equal outputs (see Results for discussion of Earth system model outputs). Given that KA is an M matrix and this metastable state approximation, we can describe the total soil organic carbon as follows:

$$\begin{array}{}\text{(2)}& C={\displaystyle \frac{u_{\mathrm{in}}}{k{Q}_{\mathrm{10}}^{(T-T_{\mathrm{0}})/\mathrm{10}}}},\end{array}$$

where C is the total organic carbon stock, u the sum of the soil inputs, and k a bulk decomposition rate that can be constructed from the decay matrix KA and allocation vector of the soil inputs b. For details of this derivation see the Supplement, Sect. “Mathematical Analysis”.

We can now describe the soil carbon stock difference between two soils with the same decay rate but different temperatures and inputs. This could either be two time points from a simulation where the soil output is close (within 10 %) of the soil inputs or a warmed treatment and a control:

$$\begin{array}{}\text{(3)}& C_{\mathrm{2}}=C_{\mathrm{1}}\left({\displaystyle \frac{u_{\mathrm{2}}}{u_{\mathrm{1}}}}{Q}_{\mathrm{10}}^{\left(T_{\mathrm{1}}-T_{\mathrm{2}}\right)/\mathrm{10}}\right).\end{array}$$

For the field sites, we assume that the relative change in inputs due to warming is negligible compared to the effect on the decomposition rate across sites and that the main driver of differences in decomposition rates between control and treatment is the warming treatment, leading us to

$$\begin{array}{}\text{(4)}& C_w=C_c{Q}_{\mathrm{10}}^{-\mathrm{\Delta }T/\mathrm{10}}.\end{array}$$

Finally, we assume that the bulk density of the soil at a given site was unaffected by the warming treatment. This allows us to use the mass percent soil organic carbon instead of the soil organic carbon density for Eq. (4).

Figure 1 The model–data fits across different Q₁₀ values for random subsets of 34 sites including the root mean squared error, and linear regression metrics R², slope, and intercept. The model is take from Eq. (4) (${\mathit{C}}_w={\mathit{C}}_c{Q}_{\mathrm{10}}^{-\mathrm{\Delta }T/\mathrm{10}}$). Slope and intercept values are shown with 2 standard deviation error bars.

2.3 Model–data integration: parameter fitting

Given the relatively small parameter space, we chose a brute-force model–data integration approach where we iteratively calculated the predicted change in soil carbon stock given the control soil carbon (Eq. 4) across a range of Q₁₀ values from 0.1 to 5 in 0.1 increments. We set the lower bound of the Q₁₀ range to 0.1 instead of 1 for two reasons. First, while it is generally accepted that warmer soil temperatures will increase soil respiration (constraining Q₁₀ > 1), it is possible that a warmer soil would result in drier soils and suppress soil respiration. In addition, numerically we wanted to bracket the expected parameter range with our prior. Data–model fits were scored using root mean squared error (RMSE) and linear regression (R², slope and intercept).

A Q₁₀ value was considered a good fit if the resulting model–data linear regression showed a low bias (slopes and intercepts within 2 standard deviations of 1 and 0, respectively); standard fit metrics like the R² and RMSE were relatively insensitive to the Q₁₀ parameter (see Fig. 1). By selecting the parameter based on model–data fit instead of deriving a direct Q₁₀ value for each site and using the distribution, we demonstrate the robustness of the model and have a clear metric to select the parameter range. To test for statistical power, we randomly sampled the data 1000 times with sample sizes from 5 to 34 sites and compared this to samples with randomly assigned control vs. warming (for each study the percent carbon of control and treatment has a 50 % chance of being switched). These random and sample-generated Q₁₀ distributions were compared using a two-sample Kolmogorov–Smirnov test to test that the distributions were statistical distinct.

2.4 Earth system model analysis

Earth system model simulations were drawn from CMIP5, the Coupled Model Intercomparison Project to support the IPCC Fifth Assessment Report (Taylor et al., 2011). We downloaded simulation outputs from the RCP 8.5 scenario, representing the “business-as-usual” scenario, including heterotrophic respiration (rh), soil temperature (tsl), and heterotrophic carbon stock (cSoil and cLitter) from the CMIP5 repository on the Earth System Federation Grid. Ten-year means were taken at the beginning and end of the 21st century for each variable (corresponding to 2006–2015 and 2090–2099). Soil temperature was averaged for the first 10 cm to correspond with experimental soil temperature readings. Soil carbon stock was calculated by adding all heterotrophic-respiration pools (including soil cSoil and litter cLitter) where multiple pools were reported. Soil carbon inputs were calculated from the monthly change in soil carbon stock plus the reported heterotrophic respiration. Model variable summaries can be found in Table S3, and processing code is documented in the Supplement.

These 20 Earth system models are built from previous models which contain 10 distinct soil sub-models (Table 1). The number of soil carbon pools in these ESMs varied from one (INM-CM4) to eight (BCC-CSM1.1), with most models having two to five pools. None of the models reported soil carbon with depth, although GFDL-ESM2G and GFDL-ESM2M document a depth-dependent model. There were three classes of pool structure for these models: independent, where there was no exchange between soil carbon pools; cascade, where pools with faster turnover times passed carbon to pools with slower turnover times; and full-feedback models, where carbon was exchanged between faster and slower pools and vis-versa. In this set of models, two of these soil models were full-feedback models (HadGEM, ISPL-CM), six were cascade pool structure (MRI-ESM1, MIROC-ESM, MPI-ESM, CLM4.0 (CESM1, CCSM4, NorESM1), CanESM2, BCC-CSM1.1), and two were independent pools (GFDL-ESM2, INM-CM4). Only two models documented an explicitly constant Q₁₀ (INM-CM2 and HadGEM2, Q₁₀ = 2), one model documented a soil-temperature-dependent Q₁₀ (CanESM2), four models documented a soil temperature sensitivity from Lloyd and Taylor which behaves very similarly to Q₁₀ = 2 at moderate temperatures (Lloyd and Taylor, 1994), and the remaining three (ISPL-CM5, GFDL-ESM2, BCC-CM1.1) all used a variation of the soil temperature sensitivity proposed in CENTURY (Parton et al., 1987, 1988) which also behaves very similarly to Q₁₀ = 2 at moderate temperatures but declines at high temperatures (Lloyd and Taylor, 1994). The ESMs considered had a single global Q₁₀, or Q₁₀ formula, dependent on soil temperature, uniformly applied to the decay pools. This documented structure should be approached with caution due to frequent lags between model development and documentation; actual values and functions may differ. For details with citations see Table 1.

Figure 2 The Q₁₀ with good one-to-one model–data fits defined in Fig. 1, at 90 % confidence interval (band) with minimum and maximum values (dotted line) and median value (solid line), across 10 different sample sizes ranging from 5 to 34, for the original data set (true: blue) and randomized case–control (random C-C: red).

Soil carbon stocks at the beginning of Earth system model simulations are typically documented to be spun up to close to steady state, and there is numerical support that this holds throughout the simulation (see Results and Fig. S3 in the Supplement). Thus Eq. (3) can be extended to the change in soil carbon stock over the 21st century. This leads to the following explicit calculation for a Q₁₀ value at each grid cell:

$$\begin{array}{}\text{(5)}& \ln \left(Q_{\mathrm{10}}\right)=\left({\displaystyle \frac{\mathrm{10}}{T_{\mathrm{m}}-T_{\mathrm{f}}}}\right)\ln \left({\displaystyle \frac{C_{\mathrm{f}}}{C_{\mathrm{m}}}}{\displaystyle \frac{u_{\mathrm{m}}}{u_{\mathrm{f}}}}\right),\end{array}$$

where the Q₁₀ value is related to the modern soil temperature T_m, future soil temperature at the end of the 21st century T_f, modern soil inputs u_m, future soil inputs u_f, modern soil carbon stock C_m, and future soil carbon stock C_f.

For soils that are very close to zero soil carbon stocks, have minimal shifts in soil temperature, or have very low soil inputs, the estimated Q₁₀ is not finite. Similarly soils which are not well described by their shift in soil temperature (for example, if there is a significant shift in the moisture regime) may have non-typical Q₁₀ values that are either less than 0.5 or greater than 5. We examined the amount of shift in soil carbon stocks associated with the four categories of Q₁₀ values (nonfinite, less than 0.5, greater than 5, or typical), as well as the spatial patterns associated with these categories.

To support the assertion that the Q₁₀ value can be calculated from relatively short timescales found in the field experiments, we examined the distribution-typical Q₁₀ values associated with similar soil temperature steps experienced by the field experiments on 1-, 5-, 10-, 15-, 20-, 50-, 75-, and 84-year timescales using 10-year mean gridded values of soil carbon stocks, soil inputs, and soil temperature. It should be noted, however, that changes in the moisture conditions over the 21st century may complicate this analysis of the Earth system model simulations; thus it is not an exact proxy for the field experiments where the control and treatment plots experienced similar baseline climate conditions and a more or less constant offset throughout the experiment.

Finally, the Q₁₀ distribution was scaled to reflect the best estimate and uncertainty from the field data. This distribution shift was done by normalizing the Q₁₀ map to the mean of the distribution and multiplying it by the experimentally derived values. The Q₁₀ correction was only applied to grids with typical Q₁₀'s (non-typical Q₁₀'s were considered to have predominately non-temperature driving variables, and their soil carbon stocks were not altered). This normalization shifted the global Q₁₀ distribution within the models to match the most common (geographically likely) Q₁₀ with the data-driven Q₁₀ value, yet by preserving the distribution we preserved other factors affecting changes in decomposition rate (i.e., moisture shifts) in the model. We then recalculated the change in soil carbon for each grid cell with this modified Q₁₀ according to Eq. (3) and calculated the global area-weighted totals.

The full analysis script and those used to generate the figures are available in the Supplement.

3 Results

From the changes in soil carbon stocks across field studies, we find a global Q₁₀ of 2.2 (90 % CI: 1.6, 2.7; R² > 0.95; root mean squared error < 2; Figs. 1, S2). The model–data fit was evaluated using a linear regression and root mean squared error (Fig. 1). While the R² of the model–data comparison was relatively insensitive to the Q₁₀ value, there was a notable improvement in the bias with Q₁₀ (as defined as the slope within 2 standard deviations of 1 and intercept within 2 standard deviations of 0). These bias-driven selection criteria were used to select Q₁₀ values from a prior range of (0.1, 5); see Methods for details.

The Q₁₀ distribution was compared with a random null distribution and was significantly distinct (Kolmogorov-Smirnov D = 0.441, p < 2 × 10⁻¹⁶; see Table S2, Figs. 2 and S1). Five to 34 randomly selected sites from the full data set were compared to a null distribution where control vs. warmed labels were randomized. The quartiles of the data subsets notably converged at a sample size of 25, where the null distribution was relatively invariant across sample size (Fig. 2). The distribution of the Q₁₀ values under null appeared lognormal, centered around 1, demonstrating no temperature effect (Fig. S1). The distribution of the Q₁₀ range for the data subsets converged to around 2.2 (Fig. S1).

The balance between gridded soil inputs and heterotrophic respiration at both the initial and final 10-year mean for the 21st century was within 10 % for over 93 % of the grid cells across all models with half of the grid cells within 0.1 %. Most models had 95 % and two models consistently had 100 % of their grid cells within 10 % – the absolute value of the net flux was within 10 % of the highest primary flux (Fig. S3). Thus, the soil inputs are on the same order of magnitude as the soil outputs. This was reflected in very similar Q₁₀ distributions regardless of whether soil inputs or heterotrophic respiration was used to derive the Q₁₀ value (Fig. S4). A notable exception to this was the MIROC-ESM model, which did see differences in inputs and heterotrophic respiration drive different Q₁₀ distributions (Fig. S4).

Figure 3 Inferred Q₁₀ values from the Earth system models (CMIP5, RCP 8.5). The color scheme is centered around the field-driven Q₁₀ median value of 2.2. Grey indicates non-typical Q₁₀ values that were nonfinite, less than 0.5, or greater then 5.

The inferred Q₁₀ values in the Earth system models derived from 10-year mean changes across different time steps (1, 5, 10, 15, 20, 50, 75, and 84 years) had similar distributions in most of the models (Fig. S4). There were minor shifts in the mode of most models, which could be attributable to changes in the moisture conditions or other (non-temperature or input) environmental variables in the simulation. Models aggregated across common land models showed marked similarity in their Q₁₀ distributions (Fig. S5). There was also an extremely high correlation between Q₁₀ values derived from soil inputs compared to those derived for heterotrophic respiration across all models (Fig. S6).

Table 2 Global model summary with multi-center mean and standard deviation for modern soil organic carbon (SOC) stocks (Pg-C), relative shift in soil inputs (u_f∕u_m), absolute change in soil temperature (dT) (^∘C), inferred mean of Q₁₀ as calculated by grid cell (see Eq. 5), change in soil organic carbon (dSOC) over the 21st century (Pg-C), and change in soil organic carbon with rescaled Q₁₀ values (1.6, 2.2, and 2.7).

The inferred Q₁₀ values in the Earth system models from the decadal average across the 21st century fell into four categories (nonfinite, less than 0.5, greater than 5, or typical; Fig. S7); however most of the change in soil carbon stocks over the 21st century occurred in grid cells with typical Q₁₀ values between 0.5 and 5 (Fig. S6). A notable exception to this trend was the MRI-ESM1 model, where roughly half of the change in carbon stocks occurred in grid cells with Q₁₀ values greater than 5 (Fig. S6). Spatially the Q₁₀ categories showed strong geographical patterns (Fig. S7). The GFDL-ESM2 models were dominated by nonfinite values at high northern latitudes (Fig. S7). MIROC-ESM, CCSM4, CESM1, and NorESM1 models were dominated by Q₁₀ values above 5 at the high northern latitudes (Fig. S7). Unless otherwise noted, only typical Q₁₀ values are addressed for the remainder of this study.

The inferred Q₁₀ values for the decadal average across the 21st century also showed strong geographic patterns (Fig. 3) and were typically unimodal (Fig. S6). MIROC-ESM and MIROC-ESM-CHEM showed the weakest spatial patterns with high grid-to-grid variation (Fig. 3). Mean Q₁₀ values fell within the 90 % CI of the field data Q₁₀, ranging between 1.8 (CESM1(CAM5), HadGEM2-ES, ISPL-CM5A, and MPI-ESM-MR) and 2.6 (MIROC-ESM-CHEM), with the multi-center Q₁₀ values at 2.0 ± 0.2 (Table 2).

When the inferred Q₁₀ values were modified to reflect the data-driven Q₁₀ range, resulting variation in the multi-center mean was almost as large as the variation across model projections (Fig. 4, Table 2). Re-centering the global Q₁₀ distribution to reflect the range of field-driven Q₁₀ values (Fig. S8) resulted in changes in soil carbon stocks over the 21st century of between −452 Pg carbon (MPI-ESM-MR) and 525 Pg carbon (HadGEM2-CC), with a best-estimate Q₁₀ (Q₁₀ = 2.2) resulting in 19 ± 155 Pg carbon (multi-center mean ± 1 SD) and a field-drive bound (Q₁₀ = 1.6, 2.7) of [248 ± 191, −95 ± 209] Pg carbon (Fig. 4, Table 2).

Figure 4 Changes in soil carbon stock (10-year means) over the 21st century from Earth system models (RCP 8.5). Grey dots are the original estimates; the open box is the soil carbon loss after the Q₁₀ is rescaled using the 2.5, 50, and 97.5 % quartiles from the field data.

4 Discussion

By capturing information about warming-induced changes to relatively undisturbed field soil carbon stocks directly rather than inferring this from soil respiration rates, this is the first study to generate field Q₁₀ estimates of soil carbon losses without needing to correct for belowground autotrophic respiration. Using a simplified version of a traditional decomposition model with a soil temperature sensitivity function, we estimate that the global Q₁₀ value is 2.2 ([1.6, 2.7] 95 % CI, Figs. 1, S2). This Q₁₀ is notably higher than previous global estimates based on field soil respiration data (Q₁₀ = 1.4 to 1.5; Bond-Lamberty and Thomson, 2010; Mahecha et al., 2010) yet well within the range of estimates from laboratory-based studies (Davidson and Janssens, 2006) as well as close to documented soil temperature sensitivity parameters (∼ 2) of Earth system models (Table 1).

This Q₁₀ range is statistically significant. Resampling the 36-study data set demonstrates the need for over 25 sites to distinguish the Q₁₀ range from random (Figs. 2 and S1). While the Q₁₀ distribution for the 34-study subset is distinct from the null distribution (Kolmogorov-Smirnov D = 0.441; p < 2 × 10⁻¹⁶), there appears to be some minor drift in the range, suggesting that more study sites could be informative, and we hope future studies will include data recently identified (van Gestel et al., 2018).

Inferring decadal-scale environmental sensitivity from an annual-scale experiment is generally controversial. However, in this case, traditional model structures assume a temperature sensitivity function that is invariant across space and time, and numerical trends in the Earth system model reflect this. In the traditional model structure the soil temperature sensitivity function is applied as a single scaler to multi-pool models, causing the relative decomposition response in both fast and slow pools to be the same (for example, Parton et al., 1987). Examining the inferred gridded Q₁₀ values from annual means across timescales from 1 to 84 years in Earth system models shows a strong similarity in the distribution of most models (Fig. S4). Similarly using soil inputs as opposed to heterotrophic respiration did not affect the distribution of the gridded Q₁₀ values, with the notable exception of MIROC-ESM, which is explained by unusual differences in soil inputs and outputs (Figs. S3, S4). Differences in the Q₁₀ distribution across timescales are likely driven then by interaction with other sensitivity functions like moisture or in shifts in the allocation of dead vegetation to different pools as the plant type distribution changes over time.

If soils are more sensitive to warming than previously expected, then how would this affect future soil carbon stocks over the 21st century? To address this question, we turned to the CMIP5 Earth system models run under RCP 8.5 (Taylor et al., 2011). In order to modify the Earth system model output to reflect the data-driven Q₁₀, we applied similar assumptions used in the field data analysis. We first examine the soil temperature sensitivity of soil carbon stocks simulated by CMIP5 Earth system models. In contrast to the field data, we take into account the effect of the change in soil inputs on soil carbon stocks in the Earth system models because these coupled simulations include CO₂ fertilization and other climate effects known to influence primary production (see Methods, Eq. 5). Though these inferred Q₁₀ values (Q₁₀ = 1.8, 2.6) fall within the uncertainty of the field-derived Q₁₀ values (Q₁₀ = 1.6, 2.7), most ESM Q₁₀ means fell under the median data Q₁₀ value of 2.2 (Table 2), implying ESMs were, on average, less sensitive to soil temperature shifts than the field-warmed data would imply. It should be noted that this inferred Q₁₀ value is not exactly the parameterized Q₁₀ value but instead a combination of the soil temperature sensitivity and other environmental sensitivities. If there were, on average, an additional constraint on respiration (such as moisture), we might expect the inferred Q₁₀ parameter to be lower than the model-parameterized Q₁₀.

There were notable regional patterns across all but two of the Earth-system-model-inferred Q₁₀'s (Figs. 3, S7). High northern latitudes tended to have either large or nonfinite Q₁₀ values, suggesting that something other than soil temperature and input shifts were driving changes in soil carbon stock. This alternative driver could be a shift in moisture regimes or dynamics driven by thaw thresholds, which could similarly affect the analysis of the field data. Additional drivers of soil decomposition dynamics, beyond temperature and inputs considered here, have the potential to explain some of the variation in the Q₁₀ range, and new model structures are being explored to take some of these mechanisms into account (Luo et al., 2015; Wieder et al., 2015a). This remains an active area of research.

Propagating this field Q₁₀ range into the ESM projections resulted in greater carbon losses from the soil by the end of the 21st century (multi-center means of the soil carbon change shifted from 88 to 19 Pg carbon) with large uncertainties; ESM multi-center standard deviation was initially 152 Pg carbon, which is half of the range of the multi-model mean attributed to Q₁₀ 95 % CI [248, −95] Pg carbon (Fig. 4, Table 2). To calculate these modified projections, means of the model-specific Q₁₀ distributions were re-centered to reflect the best-estimate Q₁₀ and associated 95 % CI from the field data analysis. By preserving the distribution within the model, we attempted to propagate soil moisture sensitivities and other model-specific effects into the modified projections. We also did not modify grid cells with non-typical Q₁₀ values (nonfinite, below 0.5, or above 5) since those grids likely governed other non-temperature drivers. The large range of carbon shifts in each ESM driven by this Q₁₀ CI confirms the importance of considering parameter uncertainty in the land carbon component of Earth system model projections. The post hoc correction that we present provides an innovative way to account for this parameter variation without the computational burden of additional ensemble runs.

This analysis includes several basic assumptions and caveats. Specifically, we assume that the difference between treatment and control is driven entirely by the soil warming effect, and those warming effects are uniform across soil carbon quality. Though warming-induced changes in soil inputs are, on average, relatively small, they are have been shown to be highly variable at similar sites (Lu et al., 2013). The analysis of field data could be extended to account for these changes in inputs in follow-up studies (Eq. 3). A large increase in soil inputs would cause an underestimation of the Q₁₀ value, while a decrease in soil inputs would cause an overestimation of the Q₁₀ value (see Eq. 3). While there is some evidence to support soil temperature sensitivity dependency on soil carbon quality (Knorr et al., 2005), there is also evidence for a uniform soil temperature sensitivity (Hicks Pries et al., 2017), as is represented in the Earth system models considered in this study (Table 1). A quality-dependent Q₁₀ would not be separable from the bulk decay term, and thus a one-pool model would be inappropriate in this case (see Supplement). In addition, the data set has acknowledged biases (see Crowther et al., 2016), which are typical of field studies.

One-pool simplification

We find that multi-year soil carbon dynamics can be well described by a one-pool model at a specific timescale in both the Earth system models and field experiment. If we restrict the decomposition models to those with either independent or cascade pool structures (that is, no carbon passed from the slow to the fast pools), then the temporal dynamics of the total soil carbon of the system at a specific timescale can be approximated by a single pool due to the fact that the lower triangular decomposition–transfer matrix is diagonalizable (see Supplement for details). While this diagonalizable property does not hold for full-feedback models where carbon is transferred from the slower to faster carbon pools, all decomposition–transfer matrices are M matrices. If we combine the positive-inverse properties resulting from this M matrix structure and assume that the soils are close to metastable state (that is, soil inputs are roughly equal to the heterotrophic respiration outputs, as we show in Fig. S3 for the Earth system models considered and would expect for soils from intact systems), then the total soil carbon can be described by a bulk decay rate that is a linear combination of the transfer coefficients, decay rates, and input allocations of the component pools (see Supplement for analytical details). This provides analytical support for the one-pool simplification seen numerically in the Earth system models in the CMIP5 project (Todd-Brown et al., 2013, 2014).

The one-pool simplifications described above are controversial assertions. The one-pool model has proven inadequate to describe laboratory incubations where heterotrophic respiration over time is compared to the soil carbon stock (Thornton, 1998; Weng and Luo, 2011). This is due to the multiple timescales considered (daily, monthly, and annual) and, more importantly, the fact that these laboratory incubations are by their nature not at steady state since any inputs to the system are generally removed. Thus this analysis would not be expected to hold for laboratory incubation, and we would further expect the bulk decay rate to change with timescales for sites undergoing rapid changes in inputs (in other words, the bulk decay rate inferred at a 1-year time step would not match the 100-year time step at a site undergoing transition from grassland to forest). Another key assumption is that soil organic carbon of different quality responds the same to warming. However, the scalar multiplier representing environmental sensitivities is independent of pools in most models (e.g., Parton et al., 1988). These scalar multiplies (like the Q₁₀ temperature sensitivity examined in this study) would be invariant to the timescale if this modeling assumption were applied to the field analysis. Finally shifts in the allocation of dead vegetation to the different soil pools would shift the bulk decay rate of the one-pool approximation (see Supplement: Analytical proofs). With these caveats in mind, we feel that the one-pool approximation is extremely valuable in analyzing soil carbon models and data.

5 Conclusion

It is still unclear how the terrestrial carbon cycle in general, and soils in particular, will respond to climate change over the 21st century. The CMIP5 models, representing our best coupled climate models to date, have a wide range of soil carbon responses over the 21st century (Todd-Brown et al., 2014). While it would be nice to have all the models agree on a tightly bound answer, the question we should be asking scientifically is, does the variation in the models reflect our best scientific understanding? Models must capture not only mean trends but also system variance and must accurately represent scientific uncertainty.

Post hoc correction of simulation results can provide some insight into known gaps in Earth system models without the computational hurdle of re-running simulation results. Previous studies have applied post hoc corrections to address nutrient limitations on net primary production (Wieder et al., 2015b), and this study demonstrates the high level of uncertainty that can be driven by the soil temperature response parameter. This study suggests that soil carbon response to warming is highly variable in the field and ESMs should increase their variability to reflect this field variation. Future studies increasing the number of field-warmed studies (van Gestel et al., 2018), as well as extending the field data to include changes in plant productivity in response to warming, would inform the field-derived Q₁₀ analysis explored here. In addition, explaining field moisture and applying that understanding to a post hoc Earth system model analysis is a logical next step.