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

Abstract

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_{\mathrm{10}}$) 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^{\mathrm{2}}$ $=$ 0.96), revealing $Q_{\mathrm{10}}$ across sites of 2.2 [1.6, 2.7] 95 % confidence interval (CI). When these field-derived $Q_{\mathrm{10}}$ 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 $\pm$ 153 Pg carbon (weighted mean $\pm$ 1 SD) to 19 $\pm$ 155 Pg carbon with a $Q_{\mathrm{10}}$-driven 95 % CI of 248 $\pm$ 191 to $-$95 $\pm$ 209 Pg carbon. On average, incorporating the field-derived $Q_{\mathrm{10}}$ 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_{\mathrm{10}}$ 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.

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_{\mathrm{10}}$ – the factor of the change in decomposition rate associated with 10 ${}^{\circ }$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_{\mathrm{10}}$ 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_{\mathrm{10}}$ 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_{\mathrm{10}}$ 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_{\mathrm{10}}$ 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_{\mathrm{10}}$ 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_{\mathrm{10}}$ 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_{\mathrm{10}}$ 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_{\mathrm{10}}$ directly from changes in soil carbon stocks (Todd-Brown et al., 2014).

Here we present a new approach to estimating the global $Q_{\mathrm{10}}$ 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_{\mathrm{10}}$) used in an integrated Earth system model. We then examine the consequences of the data-driven $Q_{\mathrm{10}}$ 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.

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\left(T\right)$, borrowed from the Century model (Parton et al., 1987, 1988); Lloyd and Taylor, taken from the recommended form from Lloyd and Taylor (1994); $Q_{\mathrm{10}}\left(T\right)$, a temperature-dependent $Q_{\mathrm{10}}$ as defined by Arora (2003); or a $Q_{\mathrm{10}}$ parameter for the $Q_{\mathrm{10}}$ function as defined in this paper.

Model center	Earth system model	Soil/land carbon sub-model	Number of soil carbon pools	Pool structure	Temperature sensitivity	Citations
BCC	BCC-CSM1.1	BCC_AVIM1.0; AVIM2; CEVSA	8	Cascade	$f\left(T\right)$	Cao and Woodward (1998), Ji et al. (2008), Wu et al. (2013)
CCCma	CanESM2	CTEM	2	Cascade	$Q_{\mathrm{10}}\left(T\right)$	Arora (2003), Arora et al. (2011), Arora and Boer (2010)
NCARNSF-DOE-NCAR  NCC	CCSM4CESM1(BGC)CESM1(CAM5)CESM1(WACCM)NorESM1-MNorESM1-ME	CLM 4.0; CN	7	Cascade	Lloyd andTaylor	Gent et al. (2011), Lawrence et al. (2011), Oleson et al. (2010), Thornton et al. (2007), Thornton and Rosenbloom (2005), Tjiputra et al. (2013)
NOAA GFDL	GFDL-ESM2GGFDL-ESM2M	LM3.0; ED	2	Independent	$f\left(T\right)$	Dunne et al. (2013), Moorcroft et al. (2001), Shevliakova et al. (2009)
MOHC${}^{*}$	HadGEM2-CCHadGEM2-ES	ROTHC	4	Feedback	$Q_{\mathrm{10}}=$2	Coleman and Jenkinson (1999), Collins et al. (2011)
INM	INM-CM4IPSL-CM5A-MRIPSL-CM5B-LR	LSM 1.0	1	Independent	$Q_{\mathrm{10}}$ $=$ 2	Bonan (1996), Volodin (2007)
IPSL	IPSL-CM5A-LR	ORCHIDEE	5	Feedback	$f\left(T\right)$	Dufresne et al. (2013), Krinner et al. (2005)
MIROC	MIROC-ESMMIROC-ESM-CHEM	SEIB-DGVM	3	Cascade	Lloyd andTaylor	Sato et al. (2007), Watanabe et al. (2011)
MPI-M	MPI-ESM-MR	JSBACH	5	Cascade	Lloyd andTaylor	Giorgetta et al. (2013), Schneck et al. (2013)
MRI	MRI-ESM1	LPJ-DGBM	3	Cascade	Lloyd andTaylor	Adachi et al. (2013), Obata and Shibata (2012), Sitch et al. (2003)

Materials and methods

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).

$Q_{\mathrm{10}}$ calculations

We calculated traditional $Q_{\mathrm{10}}$ 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: $${\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},$$ where $\mathit{C}$ is a vector of soil carbon pools with unique turnover times, $t$ is time, $u_{\mathrm{in}}$ is a scaler of soil carbon inputs, $\mathit{b}$ is an allocation vector describing how the inputs are divided between the soil carbon pools, $\mathbf{K}$ is a diagonal matrix representing the decomposition rates of the pools, ${\mathbf{Q}}_{\mathrm{10}}$ 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_{\mathrm{0}}$ an arbitrary reference temperature, and $\mathbf{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 ${\mathbf{Q}}_{\mathrm{10}}$ 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 $\mathbf{A}$ the identity matrix; cascade, where pools with faster turnover times passed carbon to pools with slower turnover times, making $\mathbf{A}$ a lower triangular matrix; and full-feedback models, where carbon was exchanged between faster and slower pools and vice versa, making $\mathbf{A}$ a fully dense matrix. In all cases $\mathbf{KA}$ is an M matrix, implying there exists an inverse with all positive entries. For the independent and cascade pools $\mathbf{KA}$ is diagonalizable, implying it can be broken down into a diagonal matrix $\mathbf{D}$ and an invertible matrix $\mathbf{P}$ such that $\mathbf{KA}$ $=$ ${\mathbf{P}}^{-\mathrm{1}}\mathbf{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 $\mathbf{KA}$ is an M matrix and this metastable state approximation, we can describe the total soil organic carbon as follows: $$C={\displaystyle \frac{u_{\mathrm{in}}}{k{Q}_{\mathrm{10}}^{(T-T_{\mathrm{0}})/\mathrm{10}}}},$$ 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 $\mathbf{KA}$ and allocation vector of the soil inputs $\mathit{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: $$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).$$ 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 $$C_w=C_c{Q}_{\mathrm{10}}^{-\mathrm{\Delta }T/\mathrm{10}}.$$ 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).

The model–data fits across different $Q_{\mathrm{10}}$ values for random subsets of 34 sites including the root mean squared error, and linear regression metrics $R^{\mathrm{2}}$, 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.

[Figure omitted. See PDF]

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_{\mathrm{10}}$ values from 0.1 to 5 in 0.1 increments. We set the lower bound of the $Q_{\mathrm{10}}$ 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_{\mathrm{10}}$ $>$ 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^{\mathrm{2}}$, slope and intercept).

A $Q_{\mathrm{10}}$ 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^{\mathrm{2}}$ and RMSE were relatively insensitive to the $Q_{\mathrm{10}}$ parameter (see Fig. 1). By selecting the parameter based on model–data fit instead of deriving a direct $Q_{\mathrm{10}}$ 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_{\mathrm{10}}$ distributions were compared using a two-sample Kolmogorov–Smirnov test to test that the distributions were statistical distinct.

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_{\mathrm{10}}$ (INM-CM2 and HadGEM2, $Q_{\mathrm{10}}$ $=$ 2), one model documented a soil-temperature-dependent $Q_{\mathrm{10}}$ (CanESM2), four models documented a soil temperature sensitivity from Lloyd and Taylor which behaves very similarly to $Q_{\mathrm{10}}$ $=$ 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_{\mathrm{10}}$ $=$ 2 at moderate temperatures but declines at high temperatures (Lloyd and Taylor, 1994). The ESMs considered had a single global $Q_{\mathrm{10}}$, or $Q_{\mathrm{10}}$ 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.

The $Q_{\mathrm{10}}$ 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).

[Figure omitted. See PDF]

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_{\mathrm{10}}$ value at each grid cell: $$\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),$$ where the $Q_{\mathrm{10}}$ value is related to the modern soil temperature $T_{\mathrm{m}}$, future soil temperature at the end of the 21st century $T_{\mathrm{f}}$, modern soil inputs $u_{\mathrm{m}}$, future soil inputs $u_{\mathrm{f}}$, modern soil carbon stock $C_{\mathrm{m}}$, and future soil carbon stock $C_{\mathrm{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_{\mathrm{10}}$ 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_{\mathrm{10}}$ 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_{\mathrm{10}}$ 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_{\mathrm{10}}$ value can be calculated from relatively short timescales found in the field experiments, we examined the distribution-typical $Q_{\mathrm{10}}$ 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_{\mathrm{10}}$ distribution was scaled to reflect the best estimate and uncertainty from the field data. This distribution shift was done by normalizing the $Q_{\mathrm{10}}$ map to the mean of the distribution and multiplying it by the experimentally derived values. The $Q_{\mathrm{10}}$ correction was only applied to grids with typical $Q_{\mathrm{10}}$'s (non-typical $Q_{\mathrm{10}}$'s were considered to have predominately non-temperature driving variables, and their soil carbon stocks were not altered). This normalization shifted the global $Q_{\mathrm{10}}$ distribution within the models to match the most common (geographically likely) $Q_{\mathrm{10}}$ with the data-driven $Q_{\mathrm{10}}$ 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_{\mathrm{10}}$ 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.

Results

From the changes in soil carbon stocks across field studies, we find a global $Q_{\mathrm{10}}$ of 2.2 (90 % CI: 1.6, 2.7; $R^{\mathrm{2}}$ $>$ 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^{\mathrm{2}}$ of the model–data comparison was relatively insensitive to the $Q_{\mathrm{10}}$ value, there was a notable improvement in the bias with $Q_{\mathrm{10}}$ (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_{\mathrm{10}}$ values from a prior range of (0.1, 5); see Methods for details.

The $Q_{\mathrm{10}}$ distribution was compared with a random null distribution and was significantly distinct (Kolmogorov-Smirnov $D$ $=$ 0.441, $p$ $<$ 2 $\times$ 10${}^{-\mathrm{16}}$; 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_{\mathrm{10}}$ values under null appeared lognormal, centered around 1, demonstrating no temperature effect (Fig. S1). The distribution of the $Q_{\mathrm{10}}$ 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_{\mathrm{10}}$ distributions regardless of whether soil inputs or heterotrophic respiration was used to derive the $Q_{\mathrm{10}}$ 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_{\mathrm{10}}$ distributions (Fig. S4).

Inferred $Q_{\mathrm{10}}$ values from the Earth system models (CMIP5, RCP 8.5). The color scheme is centered around the field-driven $Q_{\mathrm{10}}$ median value of 2.2. Grey indicates non-typical $Q_{\mathrm{10}}$ values that were nonfinite, less than 0.5, or greater then 5.

[Figure omitted. See PDF]

The inferred $Q_{\mathrm{10}}$ 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_{\mathrm{10}}$ distributions (Fig. S5). There was also an extremely high correlation between $Q_{\mathrm{10}}$ values derived from soil inputs compared to those derived for heterotrophic respiration across all models (Fig. S6).

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_{\mathrm{f}}/u_{\mathrm{m}}$), absolute change in soil temperature ($\mathrm{d}T$) (${}^{\circ }$C), inferred mean of $Q_{\mathrm{10}}$ 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_{\mathrm{10}}$ values (1.6, 2.2, and 2.7).

	SOC	Rel.	d$T$	$Q_{\mathrm{10}}$	dSOC	dSOC	dSOC	dSOC
	[Pg-C]	Inputs	[${}^{\circ }$C]		[Pg-C]	$Q_{\mathrm{10}}$ $=$ 1.6	$Q_{\mathrm{10}}$ $=$ 2.2	$Q_{\mathrm{10}}$ $=$ 2.7
BCC-CSM1.1	1050	1.40	3.7	2.2	198	312	198	134
CanESM2	1541	1.29	7.1	2.0	$-$53	239	$-$158	$-$354
CCSM4	515	1.32	4.2	1.9	6	34	$-$16	$-$45
CESM1(BGC)	515	1.29	3.8	1.9	8	29	$-$9	$-$31
CESM1(CAM5)	553	1.30	4.6	1.8	$-$1	17	$-$30	$-$56
CESM1(WACCM)	502	1.32	3.9	1.9	5	25	$-$12	$-$33
GFDL-ESM2G	1422	1.41	5.1	1.9	$-$2	25	$-$23	$-$49
GFDL-ESM2M	1278	1.38	4.5	2.0	$-$8	36	$-$24	$-$56
HadGEM2-CC	1122	1.55	8.4	1.9	285	525	118	$-$71
HadGEM2-ES	1129	1.56	8.3	1.8	259	417	41	$-$133
INM-CM4	1688	1.27	3.3	2.3	69	238	88	2
IPSL-CM5A-LR	1361	1.48	8.2	1.8	28	192	$-$205	$-$394
IPSL-CM5A-MR	1403	1.43	7.6	1.8	7	158	$-$209	$-$387
IPSL-CM5B-LR	1274	1.41	7.6	1.9	85	289	$-$63	$-$236
MIROC-ESM	2586	1.35	7.2	2.5	$-$105	363	11	$-$170
MIROC-ESM-CHEM	2588	1.30	7.3	2.6	$-$89	467	75	$-$123
MPI-ESM-MR	3110	1.31	6.3	1.8	212	461	$-$150	$-$452
MRI-ESM1	1452	1.52	4.4	2.0	415	521	374	294
NorESM1-M	547	1.31	3.7	1.9	$-$21	$-$4	$-$34	$-$51
NorESM1-ME	553	1.32	3.6	2.0	5	31	$-$6	$-$27
Multi-center mean	1403	1.37	5.4	2.0	88	248	19	$-$95
Multi-center SD	793	0.09	1.8	0.2	153	191	155	209

The inferred $Q_{\mathrm{10}}$ 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_{\mathrm{10}}$ 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_{\mathrm{10}}$ values greater than 5 (Fig. S6). Spatially the $Q_{\mathrm{10}}$ 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_{\mathrm{10}}$ values above 5 at the high northern latitudes (Fig. S7). Unless otherwise noted, only typical $Q_{\mathrm{10}}$ values are addressed for the remainder of this study.

The inferred $Q_{\mathrm{10}}$ 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_{\mathrm{10}}$ values fell within the 90 % CI of the field data $Q_{\mathrm{10}}$, 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_{\mathrm{10}}$ values at 2.0 $\pm$ 0.2 (Table 2).

When the inferred $Q_{\mathrm{10}}$ values were modified to reflect the data-driven $Q_{\mathrm{10}}$ 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_{\mathrm{10}}$ distribution to reflect the range of field-driven $Q_{\mathrm{10}}$ 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_{\mathrm{10}}$ ($Q_{\mathrm{10}}$ $=$ 2.2) resulting in 19 $\pm$ 155 Pg carbon (multi-center mean $\pm$ 1 SD) and a field-drive bound ($Q_{\mathrm{10}}$ $=$ 1.6, 2.7) of [248 $\pm$ 191, $-$95 $\pm$ 209] Pg carbon (Fig. 4, Table 2).

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_{\mathrm{10}}$ is rescaled using the 2.5, 50, and 97.5 % quartiles from the field data.

[Figure omitted. See PDF]

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_{\mathrm{10}}$ 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_{\mathrm{10}}$ value is 2.2 ([1.6, 2.7] 95 % CI, Figs. 1, S2). This $Q_{\mathrm{10}}$ is notably higher than previous global estimates based on field soil respiration data ($Q_{\mathrm{10}}$ $=$ 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 ($\sim$ 2) of Earth system models (Table 1).

This $Q_{\mathrm{10}}$ range is statistically significant. Resampling the 36-study data set demonstrates the need for over 25 sites to distinguish the $Q_{\mathrm{10}}$ range from random (Figs. 2 and S1). While the $Q_{\mathrm{10}}$ distribution for the 34-study subset is distinct from the null distribution (Kolmogorov-Smirnov $D$ $=$ 0.441; $p$ $<$ 2 $\times$ 10${}^{-\mathrm{16}}$), 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_{\mathrm{10}}$ 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_{\mathrm{10}}$ 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_{\mathrm{10}}$ 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_{\mathrm{10}}$, 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${}_{\mathrm{2}}$ fertilization and other climate effects known to influence primary production (see Methods, Eq. 5). Though these inferred $Q_{\mathrm{10}}$ values ($Q_{\mathrm{10}}$ $=$ 1.8, 2.6) fall within the uncertainty of the field-derived $Q_{\mathrm{10}}$ values ($Q_{\mathrm{10}}$ $=$ 1.6, 2.7), most ESM $Q_{\mathrm{10}}$ means fell under the median data $Q_{\mathrm{10}}$ 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_{\mathrm{10}}$ value is not exactly the parameterized $Q_{\mathrm{10}}$ 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_{\mathrm{10}}$ parameter to be lower than the model-parameterized $Q_{\mathrm{10}}$.

There were notable regional patterns across all but two of the Earth-system-model-inferred $Q_{\mathrm{10}}$'s (Figs. 3, S7). High northern latitudes tended to have either large or nonfinite $Q_{\mathrm{10}}$ 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_{\mathrm{10}}$ 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_{\mathrm{10}}$ 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_{\mathrm{10}}$ 95 % CI [248, $-$95] Pg carbon (Fig. 4, Table 2). To calculate these modified projections, means of the model-specific $Q_{\mathrm{10}}$ distributions were re-centered to reflect the best-estimate $Q_{\mathrm{10}}$ 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_{\mathrm{10}}$ 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_{\mathrm{10}}$ 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_{\mathrm{10}}$ value, while a decrease in soil inputs would cause an overestimation of the $Q_{\mathrm{10}}$ 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_{\mathrm{10}}$ 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_{\mathrm{10}}$ 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.

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_{\mathrm{10}}$ 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.

The code for this study is included in the Supplement.

The supplement related to this article is available online at: https://doi.org/10.5194/bg-15-3659-2018-supplement.

KTB and TWC both contributed to the writing of the manuscript. KTB was responsible for the analysis. TWC coordinated the field data set. BZ was responsible for the mathematical analysis and made contributions to the revisions of the manuscript.

The authors declare that they have no conflict of interest.

Acknowledgements

We thank Ben Bond-Lamberty and Vanessa Bailey for their reviews of the manuscript before submission, Will Wieder and Mark Bradford for discussions in formulating this analysis, and Donald Jorgensen for his work on the graphical abstract. We would also like to thank the reviewers for their feedback.

Part of the research described in this paper was conducted under the Laboratory Directed Research and Development Program at the Pacific Northwest National Laboratory, a multiprogram national laboratory operated by Battelle for the U.S. Department of Energy; Katherine Todd-Brown is grateful for the support of the Linus Pauling Distinguished Postdoctoral Fellowship program. TWC was supported by a Marie Skłodowska-Curie Action fellowship. Bin Zheng was supported in part by the U.S. Department of Energy, Office of Science, Office of Advanced Scientific Computing Research, as part of the Uncertainty Quantification for Complex Systems project.

We thank the experimentalists who generously made their data available for this analysis.

We acknowledge the World Climate Research Programme's Working Group on Coupled Modelling, which is responsible for CMIP, and we thank the climate modeling groups (listed in Table S3) for producing and making available their model output. For CMIP the U.S. Department of Energy's Program for Climate Model Diagnosis and Intercomparison provides coordinating support and led development of software infrastructure in partnership with the Global Organization for Earth System Science Portals. Edited by: Michael Weintraub Reviewed by: Chris Jones, Melanie Mayes, and three anonymous referees