# Projected soil carbon loss with warming in constrained Earth system models

Shuai Ren, Tao Wang, Bertrand Guenet, Dan Liu, Yingfang Cao, Jinzhi Ding, Pete Smith, Shilong Piao

► To cite this version:

Shuai Ren, Tao Wang, Bertrand Guenet, Dan Liu, Yingfang Cao, et al. Projected soil carbon loss with warming in constrained Earth system models. Nature Communications, 2024, 15 (1), pp.102. 10.1038/s41467-023-44433-2. hal-04790877

HAL Id: hal-04790877

https://ens.hal.science/hal-04790877v1

Submitted on 19 Nov 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Projected soil carbon loss with warming in constrained Earth system models

Received: 31 March 2023

Accepted: 13 December 2023

Published online: 02 January 2024

<figure><img src="image_3.png" /></figure>

Check for updates

Shuai Ren Ⓔ<fnref n="1" />,², Tao Wang Ⓔ<fnref n="1" />✉, Bertrand Guenet³, Dan Liu<fnref n="1" />, Yingfang Cao<fnref n="1" />,², Jinzhi Ding Ⓔ<fnref n="1" />, Pete Smith Ⓔ⁴ & Shilong Piao Ⓔ<fnref n="1" />,⁵

The soil carbon-climate feedback is currently the least constrained component of global warming projections, and the major source of uncertainties stems from a poor understanding of soil carbon turnover processes. Here, we assemble data from long-term temperature-controlled soil incubation studies to show that the arctic and boreal region has the shortest intrinsic soil carbon turnover time while tropical forests have the longest one, and current Earth system models overestimate intrinsic turnover time by 30 percent across active, slow and passive carbon pools. Our constraint suggests that the global soils will switch from carbon sink to source, with a loss of 0.22–0.53 petagrams of carbon per year until the end of this century from strong mitigation to worst emission scenarios, suggesting that global soils will provide a strong positive carbon feedback on warming. Such a reversal of global soil carbon balance would lead to a reduction of 66% and 15% in the current estimated remaining carbon budget for limiting global warming well below 1.5 °C and 2 °C, respectively, rendering climate mitigation much more difficult.

Soil is the largest reservoir of terrestrial organic carbon<fnref n="1" />, and there is compelling experimental evidence of accelerated soil carbon loss with warming²,³, suggesting that soils may act as a positive carbon feedback on climate change. However, it is still not clear how global warming will affect soil carbon dynamics⁴,⁵, either in terms of the magnitude of the effect or even its sign. One significant, and poorly understood, component of the system is soil carbon turnover<fnref n="1" />,⁶, which is defined as the average time it takes for a carbon atom to enter and leave the soil system⁷. In the sub-models of Earth system models (ESMs), the soil organic carbon is generally viewed as a heterogeneous mix of two to several pools with different degrees of decomposability⁸, and the soil carbon decomposition is commonly dictated by intrinsic decay constants, which are modified only by abiotic factors (temperature and moisture) with few spatially-uniform parameters⁹⁻¹². The models vary greatly in the assumed values of constants and parameters¹², and the observational constraints on these values are generally lacking<fnref n="1" />,⁴,¹³. This

deficiency in model structure and/or parameters is manifested in simulations of soil carbon turnover time at local to regional scales which diverge both from each other¹²,¹⁴ and from radiocarbon-based observations⁶,¹⁵, reducing confidence in forecasts of how soil carbon stock will respond in the long term.

Here, we take advantage of the growing number of long-term soil incubation experiments to generate a spatially-explicit understanding of intrinsic soil carbon turnover ($\tau_i$) across the globe¹⁶,¹⁷. Notably, soil $\tau_i$ is representative of the theoretical carbon turnover time under optimal conditions. While various environmental constraints such as freezing and physical protection could inhibit the achievement of this theoretical value, leading to longer apparent value of $\tau_i$ in the real-world settings⁸ (Methods). Using global soil $\tau_i$ observations, we then develop an observationally-calibrated three-pool model, which includes emerging concepts of controls on soil carbon stabilization (e.g., soil physical-chemical protection and priming effect)¹⁸⁻²¹, to

<footnote n="1">¹State Key Laboratory of Tibetan Plateau Earth System, Environment and Resources (TPESER), Institute of Tibetan Plateau Research, Chinese Academy of Sciences, Beijing, China. ²University of Chinese Academy of Sciences, Beijing, China. ³Laboratoire de Géologie, École normale supérieure, CNRS, PSL University, IPSL, Paris, France. ⁴Institute of Biological and Environmental Sciences, School of Biological Sciences, University of Aberdeen, Aberdeen AB24 3UU, UK. ⁵Institute of Carbon Neutrality, Sino-French Institute for Earth System Science, College of Urban and Environmental Sciences, Peking University, Beijing, China. ✉e-mail: twang@itpcas.ac.cn</footnote>
constrain projected soil carbon changes in ESMs by the end of this century under different emissions scenarios.

## Results and discussion

We assessed the soil carbon turnover time by assembling a global database of aerobic carbon dioxide (CO₂) efflux data from temperature-controlled soil incubation studies conducted at 178 sites covering eight biomes, ranging from Arctic permafrost to dry Mediterranean forests (Fig. 1a). Due to the large heterogeneity of soil carbon, we represent the time evolution of soil CO₂ efflux (SR) using a three-pool model which partitions the total soil carbon (Cₜot) into three reservoirs, with each decaying at its own intrinsic turnover rate (that is, the inverse of turnover time τ)¹⁶,¹⁷.

$$SR = \sum_{p=1}^{3} \frac{1}{\tau_p} \cdot C_{tot} \cdot f_p \quad (1)$$

where f is the partitioning coefficient of Cₜot for each pool (p). Using deconvolution analysis, the pool-specific τ values were quantified at their own incubation temperatures, and then scaled to a common temperature of 15 °C (ref. 22; see Methods). Note that the soils were incubated at constant temperatures, while other environmental factors such as soil moisture were generally maintained at the optimum level²³, suggesting that the inverted carbon pool-specific values of τ can be considered as intrinsic ones without environmental constraints.

## Global patterns of intrinsic soil carbon turnover times

There were considerable variations in soil τᵢ within each biome and across biomes for each of the three carbon reservoirs (Fig. 1b–d). We

then used a boosted regression trees (BRT) model to determine the dominant environmental drivers (e.g., vegetation growth, local climate and soil attributes consisting of bulk density, pH value, organic carbon, total nitrogen, C:N ratio and soil texture) of the cross-site τᵢ variations. Of the tested predictors, mean annual temperature (MAT) was the most important variable in explaining cross-site variability of τᵢ for all three soil carbon pools, with importance values of 27–42% across the different pools (Supplementary Fig. 1). The normalized difference vegetation index (NDVI) and soil organic carbon came second in explaining τᵢ variability in the fast soil carbon pool (Cₐₗt) and in the slow (Cₐₗl) and passive (Cₐₗsiv) carbon pools. Specifically, soil τᵢ had positive correlations with MAT (Supplementary Figs. 2–3), with lower MAT being associated with shorter τᵢ.

Based on the empirical relationships between environmental predictors and soil τᵢ, we generated a predictive model that could explain more than 75% of the cross-site variability in τᵢ of each soil carbon pool without bias (Supplementary Fig. 1). This predictive model enabled us to interpolate τᵢ of each soil carbon pool to a depth of 1 m across global soils (Fig. 2a, d, g). We derived a global τᵢ of 0.3 yr for Cₐₗt, with a spatial variation ranging between 0.06 yr (first percentile) and 0.64 yr (99th percentile) (Fig. 2a and Supplementary Table 1). The global Cₐₗl and Cₐₗsiv τᵢ are 6.68 (2.3–12.1) and 398 (89–696) yr, respectively, with the values in parentheses denoting the spatial range between the first and 99th percentiles (Fig. 2d, g and Supplementary Table 1). We also calculated carbon-weighted τᵢ to a depth of 1 m using the fraction of each soil carbon pool. The fraction of each soil carbon pool at the global scale was extrapolated from that derived at site level, using the empirical relationships between fraction of specific pools and environmental drivers across sites (see Methods). The global carbon-weighted τᵢ is 316 yr with a spatial variation ranging between

<figure><img src="image_4.png" /><figcaption>**Fig. 1 | Distribution of sample locations and intrinsic soil carbon turnover data.**
a The spatial distribution of long-term temperature-controlled soil incubation experiments. b–d The boxplots showing the distributions of intrinsic soil carbon turnover times (τᵢ) of Cₐₗt (b), Cₐₗl (c) and Cₐₗsiv (d) that are inverted from all experiments within each of the eight biomes, respectively. The eight biomes are</figcaption></figure>

tropical forest (n = 53), temperate forest (n = 73), boreal forest (n = 33), cropland (n = 59), grassland (n = 69), shrubland (n = 8), tundra (n = 64) and wetland (n = 15), respectively. On each box, the central black line marks the median, the edges of the box correspond to the 25th and 75th percentiles, the whiskers extend to the range of the data, and the outliers are shown as dots.
<figure><img src="image_5.png" /><figcaption>**Fig. 2 | Global distributions of intrinsic soil carbon turnover times.**
a, d, g Observation-derived intrinsic soil carbon turnover times ($\tau_i$) at 1 m depth that are upscaled from 374 data points at 15 °C using the machine learning model for $C_{fast}$ (a), $C_{slow}$ (d) and $C_{passive}$ (g), respectively. b, e, h Standard deviations of global soil $\tau_i$ estimates derived from different combinations of data sets. c, f, i Latitudinal</figcaption></figure>

profiles of soil $\tau_i$, aggregated at 0.5° latitudinal resolution. The black and red lines indicate observation-derived and ESMs-averaged soil $\tau_i$, respectively. The pink dashed lines are soil $\tau_i$ specified in each model (Supplementary Table 2). The red shading represents the standard deviation. The x-axis data in c, f, i are log-transformed.

81 yr (first percentile) and 609 yr (99th percentile) (Supplementary Table 1).

The longest $\tau_i$ values were found in tropical forests (0.43, 9.41 and 601 yr for $C_{fast}$, $C_{slow}$ and $C_{passive}$, respectively), shrublands (0.31, 7.8 and 467 yr) and grasslands (0.28, 6.93, and 418 yr), whereas tundra (0.25, 4.48 and 288 yr) and boreal forests (0.32, 5.58, and 326 yr) have the shortest values (Supplementary Table 1). We also find that the $\tau_i$ values of the three carbon pools vary considerably with latitude, with higher values in the tropical zone between 20° N and 20° S (mean $\tau_i$ is 0.35, 9.1 and 543 yr) but lower values in the northern high latitudes (0.27, 4.99 and 269 yr above 50° N) (Fig. 2a, d, g). To explore the extent to which the intrinsic turnover times translate into the actual ones due to environmental constraints®, we compared our estimates with the radiocarbon-derived carbon age as a surrogate of realized or apparent turnover times (Supplementary Table 1). We find that soil $\tau_i$ is more than 16 times shorter than the actual one across the globe (316 yr versus 5,238 yr), and this value is largest in tundra (58.6) and boreal forests (18) but smallest in tropical forests (4.4) (Supplementary Table 1; Supplementary Fig. 4).

## Observational constraints on projected soil carbon dynamics in ESMs

Next, we evaluated soil $\tau_i$ specified in the five ESMs which have three soil carbon pools ($C_{fast}$, $C_{slow}$ and $C_{passive}$) archived in the Coupled Model Intercomparison Project Phase 6 (CMIP6). Soil $\tau_i$ in ESMs are generally described as global constants without any spatial variability (Supplementary Table 2). We found that ESMs overestimated soil $\tau_i$ in high-latitude ecosystems for all carbon pools (Fig. 2 and Supplementary Fig. 5), with a factor of 1.2, 2 and 1.7 for $C_{fast}$, $C_{slow}$ and $C_{passive}$, respectively. By contrast, in the tropics, soil $\tau_i$ are overestimated to a less extent, and even underestimated in some regions (Supplementary Fig. 5). Overall, the model ensemble overestimated soil $\tau_i$ by about 30% across the three different pools globally, with larger biases in $C_{slow}$ (58%) than those in $C_{fast}$ (12%) and $C_{passive}$ (17%) (Fig. 2). The model-data bias of soil $\tau_i$ could be attributed to the omission of critical microbial processes from ESMs, such as thermal adaptation²⁴. Specifically, microbial turnover rates have been shown to adjust to temperature changes via biochemical trade-offs in enzyme and cell membrane structure and function²⁵,²⁶. Low temperatures typically select for
<figure><img src="image_6.png" /><figcaption>**Fig. 3 | Projected changes in global soil carbon stock.** Changes in global soil carbon stock between the current period (2005–2014) and the end of the century (2080–2099) from the original and constrained CMIP6 models under SSP1-2.6 (blue), SSP2-4.5 (orange) and SSP5-8.5 (red) emissions scenarios, respectively.</figcaption></figure>

enzymes and/or membranes that are highly flexible to efficiently alter conformation and facilitate interactions²⁷. As a result, cold-adapted microbial communities have faster growth and respiration rates than the warm-adapted when compared at common temperatures⁵,²⁶. Future modeling efforts should seek a spatial representation of soil intrinsic turnover parameters especially for $C_{slow}$, e.g., by incorporating microbial metrics (such as thermal adaptation²⁴, species composition²⁵) into ESMs to build confidence in predicting soil carbon-climate feedback.

The large deviation of modeled $\tau_i$ from the observations could lead to substantial biases in projections of soil carbon stock changes. To remove simulation biases due to misinterpretations of $\tau_i$ in ESMs, we firstly developed a reduced-complexity three-pool model, which could partition soil organic carbon into components with different intrinsic turnover rates, to emulate soil carbon simulations in much more complex ESMs (see Methods). Although our constructed reduced-complexity model is capable of mimicking ESMs soil carbon dynamics and making projections (Supplementary Fig. 6), it is limited to rudimentary processes without integrating emerging knowledge of controls on soil carbon turnover time. To improve the realism of the model, we constructed a refined reduced-complexity model to constrain ESM projections under three different emissions scenarios (SSPI-2.6, SSP2-4.5 and SSP5-8.5; Supplementary Fig. 7). In the refined reduced-complexity model, we described the soil carbon decay as a function of $\tau_i$ that is not only modified by climate factors, but also by the two competing processes (mineral protection and the rhizosphere priming effect) (see Methods). We forced the model with climatological mean (2000–2014) climate variables and satellite-based net primary productivity (NPP), and calibrated the model against both the gridded historical soil carbon pool and upscaled field determinations of soil carbon age (as a proxy of apparent soil carbon turnover times) derived from radiocarbon measurements¹⁵ to minimize the model-data mismatch for the historical period (see Methods). The parameter uncertainties were obtained from different combinations of soil organic carbon (SoilGrids and the Harmonized World Soil Database) and NPP data sets (MODIS and Global Inventory Modeling and Mapping Studies) (Supplementary Table 3 and Supplementary Fig. 8).

The refined reduced-complexity model, constructed in this way, was run forward to 2100 using bias-corrected NPP and climate forcing (temperature and precipitation) under three different emissions scenarios (see Methods). An ensemble of original ESMs projected an increase of 18.6–29.8 petagrams of carbon (PgC) in global soils by the end of this century relative to the period 2005–2014 across the three scenarios (Fig. 3 and Supplementary Fig. 9). In contrast, constrained

ESM simulations using the refined reduced-complexity model showed that global soils will lose 19.1 Pg of carbon under SSP1-2.6, and this loss will be further escalated to 24 and 45.2 PgC under SSP2-4.5 and SSP5-8.5, respectively (Fig. 3). These projections translate into an estimated annual soil carbon loss rate of 0.22–0.53 PgC from now until the end of this century under different emissions scenarios. This finding is generally consistent with the expectation that the intrinsic turnover time was overestimated in complex ESMs (Fig. 2). This is because faster intrinsic soil carbon turnover would lead to a more rapid decline in soil carbon stock in response to warming, although environmental constraints such as climatic limitation and physical protection would inhibit the achievement of this intrinsic value³. Such soil carbon losses could largely offset increases in plant biomass due to the fertilization effect of rising $CO_2$, thereby reducing the potential capacity of land biosphere carbon sequestration in future. For example, in an intermediate emission scenario, the increase in $CO_2$ expected by 2100 would enhance the global plant biomass by 59 PgC using an empirical upscaling of $CO_2$ fertilization effect²⁹. However, this enhancement in biomass carbon stock would be substantially counterbalanced by soil carbon losses (24 PgC).

In most regions, the soil carbon stock switched from a sink in the original simulations to a source in the constrained simulations (Fig. 3). Such a reversal of the soil carbon balance is most prominent in the arctic tundra and boreal regions (Supplementary Fig. 10), and occurred there even without invoking warming-induced permafrost soil carbon degradation³⁰,³¹. The arctic tundra and boreal region have the shortest soil-carbon intrinsic turnover time of all the regions (Fig. 2), and the release of environmental constraints *via* warming and the plant rhizosphere priming effect will largely decrease the apparent soil carbon turnover time, thereby inducing large soil carbon losses. For example, a recent synthesis showed that rhizosphere priming could induce an additional 40 PgC loss from northern soils by 2100³². Furthermore, the entire tropical forest soils, including that in Amazonia, will also become a net carbon source of 8.2–21.4 PgC (Supplementary Fig. 10). Our results are in broad agreement with previous in situ experiments which showed that warming could accelerate the rate of soil carbon losses in tropical forests³,³³. However, the magnitude of soil carbon losses from tropical forests under SSP5-8.5 (21.4 PgC) is lower than a recent estimate (65 PgC) that is simply extrapolated from a two-year warming experiment in a tropical forest under 4 °C warming³. The relatively low soil carbon loss is also due to the large increase in soil carbon input due to $CO_2$ fertilization in ESMs. While emerging evidence suggests that this $CO_2$ fertilization effect in ESMs may have been overestimated³⁴, and our estimates of soil carbon losses are then likely underestimated.

## Implications for the global remaining carbon budget

The vanished capacity of soils to sequester $CO_2$ suggests that a more aggressive strategy toward emissions reduction is required to realize the pledges of the Paris Agreement. A likely (50%) chance of keeping warming well below the 1.5 °C and 2 °C temperature targets requires that the maximum permitted carbon emissions remain below 68 PgC and 327 PgC from the start of 2022³⁵, respectively. We further estimated the amount of carbon absorbed by global soils to be less than the original ESM projections at the end of this century (Supplementary Fig. 11). The results imply that the currently estimated remaining carbon budget should be reduced by nearly 15–66% over the course of this century to achieve the warming goals, rendering climate mitigation much more difficult.

In summary, we synthesized long-term incubated-soil $CO_2$ flux measurements to provide observationally-constrained estimates of future soil carbon stock changes, using an optimized three-pool soil carbon model which includes the emerging concept of controls on soil carbon turnover time. Our observational constraints supported the current idea that global soils will create a strong positive carbon
feedback to climate change in a warmer world$^{4,5,36}$, and indicated that such soil carbon-climate feedback is substantially underestimated in current ESMs due to poor parameterizations of soil carbon turnover processes. The observationally-calibrated soil carbon model, with its spatially-varying estimates of parameters, is effective at capturing soil carbon turnover across diverse terrestrial biomes. It indicates the areas which require improvement in conventional models which parameterize carbon turnover processes in a spatially-uniform manner. The study provides important insights into the potential use of an observationally-constrained soil carbon model on projected soil carbon-climate feedback. However, our constrained projection is still subject to uncertainties, due to the omission of deep permafrost carbon dynamics$^{11}$, a divergent projection of soil carbon inputs$^{37}$, and a lack of feedbacks between soil carbon dynamics and nutrients. By including more data, particularly from under-sampled regions, such as Africa, central and southern Asia and some high latitudes, similar constrained projection studies are likely to provide further value to this area of research.

## Methods

### Collection of long-term soil incubation experiments

The intrinsic soil carbon turnover time reflects kinetic properties of various soil organic compounds under optimal conditions$^8$. Generally, molecular structures with greater complexity exhibit an increased resistance to decomposition$^{38}$, resulting in a longer $τ_i$. However, in real-world conditions, various environmental constraints, including freezing, flooding and physical protection, can dampen decomposition processes$^8$, frequently leading to a longer apparent value of $τ_i$. Laboratory soil incubation experiments, with incubation duration varying from days to years, are a widely-used approach to estimate the intrinsic decomposability of soil carbon pools by measuring soil CO₂ fluxes under controlled conditions$^{17,39}$. Results from such soil incubation experiments are invaluable for informing process-based ecosystem models about carbon pool sizes and their intrinsic turnover rates$^{40}$. In contrast to short-term incubations (a few days to weeks), where measured soil CO₂ fluxes mostly originate from respiration in the fast-cycling carbon pool, long-term incubations can provide more information on the decomposability of slow-cycling carbon pools$^{41,42}$.

Here, we assembled a raw dataset of soil CO₂ fluxes from soil incubation experiments to assess the intrinsic decomposability of different soil carbon pools. We included only soil incubations which satisfied the following criteria: (1) the length of the incubation experiment was longer than 6 months; (2) the aerobic CO₂ production must have been measured over the parallel-warming incubation experiments. If more than one incubated temperature was available, we selected the high-end temperature treatment, as a higher incubation temperature generally leads to a quicker depletion of the fast-cycling carbon pool and therefore enables more information on the slow-cycling carbon decomposability to be retrieved$^{17}$; (3) the initial soil carbon concentration must be available; (4) the soils were incubated without any substrate addition; (5) soil respiration rates in the initial phase must be higher than those at the end of incubation; and (6) data can be taken directly from tables or extracted from figures using the GetData (v.2.25) software. In total, we gathered 102 peer-reviewed publications covering 178 sites (Fig. 1), resulting in a total of >5000 time-series of soil CO₂ flux data. The soil CO₂ flux data were taken from a range of depths from 5 cm to >1 m, with incubation temperatures ranging from 4 to 35 °C and incubation durations ranging from 180 d to >10 yr. The metadata for each sampling site included location (latitude and longitude), climate (mean annual temperature and precipitation), soil physicochemical attributes (bulk density, pH, soil organic carbon, total nitrogen, C:N ratio, and sand and silt content), vegetation types, sampling depth, incubation temperatures and duration. The samples from which data were obtained can be classified into tropical forest (23° S–23° N; n = 53), temperate forest (23–50° N

and S; n = 73), boreal forest (> 50° N; n = 33), grassland (n = 69), cropland (n = 59), shrubland (n = 8), tundra (n = 64) and wetland (n = 15) (Fig. 1). For sampling sites where climate and soil attributes were not available, we extracted data from WorldClim version 2.0$^{43}$ and SoilGrids$^{44}$, respectively, based on the published location and depth information. The normalized difference vegetation index from the MODIS was also obtained for each sampling site as a measure of site productivity.

### Inversion of intrinsic carbon turnover times from soil incubation experiments

We assumed that the total CO₂ flux is composed of contributions from the three different carbon pools (that is $C_{fast}$, $C_{slow}$ and $C_{passive}$). For each sampling site, we fitted the time series of total soil CO₂ fluxes using the three-pool carbon decomposition model$^{17}$ described below. In this model, the total soil CO₂ flux (R, in mg C g⁻¹dw d⁻¹) was the sum of the respiration rates (r) derived from the three carbon pools (p) with different sizes and turnover rates (Eq. (2)).

$$R = \sum_{p=1}^{3} r_p = \sum_{p=1}^{3} k_p \times C_{tot} \times f_p \quad (2)$$

$$f_p = \frac{C_p}{C_{tot}}, \sum_{p=1}^{3} f_p = 1 \quad (3)$$

where the pool-specific respiration rate ($r_p$) was computed as the product of intrinsic pool-specific decay rate ($k_p$, in d⁻¹; the inverse of $τ_i$), the total initial carbon pool ($C_{tot}$, in mg C g⁻¹dw) and a partitioning coefficient ($f_p$). The partitioning coefficient describes the ratio of carbon pool p to the total carbon pool, with the sum of the three coefficients equal to unity (Eq. (3)).

The model parameters (m), consisting of k and f for each pool, were optimized against measured soil CO₂ fluxes (O) using a Bayesian probabilistic inversion approach$^{45}$. This approach states that the posterior probability density function (PDF) of the model parameters $p(m|O)$, conditional on the data (O), can be obtained by applying Bayes theorem$^{46}$.

$$p(m|O) \propto p(O|m)p(m) \quad (4)$$

where $p(m)$ is the prior PDF and $p(O|m)$ is the PDF of O conditional on m, also called the likelihood function. For constructing $p(m)$, we first specified the range of each parameter according to values obtained from the literature (Supplementary Table 4) and then assumed that they are uniformly distributed over this range. The likelihood function $p(O|m)$ summarizes the difference between the simulated and measured CO₂ flux data. Through assuming that the errors of measured CO₂ flux data follow a normal distribution with zero mean, the likelihood function $p(O|m)$ is given by,

$$p(\mathbf{O}|\mathbf{m}) \propto \exp \left\{ -\frac{1}{2\sigma^2} \sum_{t \in \text{obs}(S_i)} [S_i(t) - O_i(t)]^2 \right\} \quad (5)$$

where $O_i(t)$ and $S_i(t)$ are the measured and simulated CO₂ fluxes, and σ is the standard deviation of the measured CO₂ fluxes.

To derive the posterior distribution of model parameters analytically, we resorted to the Metropolis-Hastings (M-H) algorithm, a powerful Markov Chain Monte Carlo (MCMC) technique for simulating complex and nonstandard multivariant distributions, which would iteratively search for the optimum feasible solution$^{47,48}$. During the inversion process, the sum of the fractions (f) of all three carbon pools must be equal to unity, and the intrinsic decay rate (k) should be largest in the fast carbon pool and smallest in the passive one. We then
computed the maximum likelihood estimates for well-constrained parameters (e.g., the decay rates from the fast and slow carbon pools), the means of the poorly-constrained parameters (e.g., the decay rate from the passive pool), and the confidence intervals for all parameters from their posterior distributions. Overall, the performance of the three-pool model in fitting the measured soil CO₂ fluxes was good (Supplementary Figs. 12 and 13).

In the compiled data, soil incubation temperatures ranged from 4 to 35 °C, making it difficult to compare τᵢ across studies in a quantitative manner. Therefore, we adjusted τᵢ at different incubation temperatures to a reference temperature of 15 °C by using the following equation²²,⁴⁹.

$$ \tau_{i,15} = \tau_{i,T} \times Q_{10}^{T_{15}} \quad (6) $$

where τᵢ,₁₅ and τᵢ,T are soil intrinsic turnover times at the reference temperature of 15 °C and the incubation temperature T (°C), respectively. Q₁₀ is a temperature sensitivity parameter defined as the factor by which soil respiration increases with a 10 °C increase in temperature. Here, we assumed that Q₁₀ varies with temperature, and their empirical function was derived from a previous synthesis analysis of Q₁₀ and temperature from laboratory studies across various ecosystems at the global scale²³ (Supplementary Fig. 14).

## Upscaling of site-level intrinsic carbon turnover times to the global level

We used boosted regression trees (BRT) to assess the relative importance of the independent variables, including local climate (mean annual temperature and mean annual precipitation (MAP)), vegetation productivity (NDVI) and edaphic properties (bulk density, pH value, organic carbon content, total nitrogen content, C:N ratio, silt and sand content), on spatial variability in τᵢ for each carbon pool (at 15 °C; Supplementary Table 5). Note that we selected local climate variables because they can affect the structure of soil microbial communities and then the soil carbon decomposition⁵⁰. BRT is an ensemble method that combines the strengths of regression trees (tree-based models that relate a predicant to predictors through recursive partitioning) and boosting algorithms (using large numbers of relatively simple tree models to give improved predictive performance)⁵¹. Our approach, which is superior to most traditional methods, can handle different types of predictor variables and interaction effects between predictors, and does not require data transformation or outlier elimination. We used a grid-search procedure, using ten-fold cross-validation, to select the best hypermeter combination of BRT modeling with the lowest cross-validation root mean square error (Supplementary Table 6). The BRT analysis was performed using the gbm and caret packages in R 4.0.5. The soil τᵢ data were log10 transformed before starting the analysis.

The constructed BRT model was able to explain more than 75% of the variances in τᵢ across different sites (Supplementary Fig. 1), and was subsequently used to produce global maps of τᵢ for each pool at a spatial resolution of 0.1°. The uncertainties of these parameters were further generated by forcing the constructed BRT models with combinations of different sources of climate and edaphic data (Fig. 2). For climate predictors, we used mean annual temperature and mean annual precipitation from both the Climate Research Unit (CRU) version 4.01³² and WorldClim version 2.0⁴³; for edaphic factors (soil physicochemical properties), we considered both the Global Soil Dataset for Earth System Modeling (GSDE)³³ and the SoilGrids⁴⁴ data (Supplementary Table 5). We also tested the extent of our extrapolations and found that our soil samples spanned most environmental conditions around the globe (Supplementary Fig. 15). Despite this, certain regions (such as Africa, central and southern Asia and some high latitudes) are underrepresented by our samples (Fig. 1). Thus, more long-term soil incubations are urgently needed in these specific regions.

To calculate soil carbon-weighted τᵢ (τᵢ,w in Eq. (7)), we used the same method as used to obtain global maps of fractions in different soil carbon pools (fₚ in Eq. (7); Supplementary Fig. 16).

$$ \tau_{i,w} = \sum_{p=1}^{3} \tau_{i,p} \times f_{p} \quad (7) $$

## Development of a reduced-complexity three-pool model

From the CMIP6 historical and future simulations for three different emissions scenarios (SSP1-2.6, SSP2-4.5 and SSP5-8.5), we selected only those ESMs which reported the total soil carbon stock partitioned into the three discrete soil carbon components: cFast, cMedium and cSlow. For ESMs with litter or woody debris carbon pools (cLitter and cCwd), we combined these pools with cFast to form the fast-cycling carbon pool⁶. Since ESMs do not provide depth information for soil carbon, we assumed that the carbon was stored within the top one meter of the soil⁶,¹². The models used here are ACCESS-ESM1-5, CESM2-WACCM, IPSL-CM6A-LR, NorESM2-LM and TaiESM1 (Supplementary Table 7).

We developed a reduced-complexity three-pool model to approximate soil carbon simulations in each grid cell of the five ESMs. This three-pool model partitions soil organic carbon into components with different intrinsic turnover rates (Eqs. (8–10)).

$$ \left\{ \begin{aligned} \frac{dC_f}{dt} &= NPP(t) - k_f \times C_f(t) \\ \frac{dC_s}{dt} &= k_f \times C_f(t) \times r_f - k_s \times C_s(t) \\ \frac{dC_p}{dt} &= k_s \times C_s(t) \times r_s - k_p \times C_p(t) \end{aligned} \right. \quad (8) $$

where r_f and r_s are transfer coefficients for carbon flowing from fast to slow pools and from slow to passive pools, respectively. C_f, C_s, C_p represent the soil carbon stock of the fast, slow and passive pools, respectively. k_f, k_s and k_p are the actual carbon decay rates (yr⁻¹) calculated using Eq. (9).

$$ k = \frac{1}{\tau_a} = \frac{1}{\tau_i} \times F(T) \times F(P) \quad (9) $$

$$ F(T) = Q_{10}^{\left(\frac{T-r_{ref}}{10}\right)}; \quad F(P) = P^b \quad (10) $$

Where τₐ denotes the actual-, and τᵢ the intrinsic, carbon turnover time (yr) at a reference temperature (T ref). T and P are the ambient temperature (°C) and precipitation (mm), respectively. F(T) is a Q₁₀-based standard exponential function to represent the temperature modifier of τᵢ, F(P) represents a moisture modifier that increases with annual precipitation, normalized to maximal annual precipitation for each ESM, using an exponential function where b is greater than zero¹².

In the three-pool model, five of eight parameters, including intrinsic soil carbon turnover time of each pool, Q₁₀ and associated T ref were from original ESMs (Supplementary Table 2). In particular, for ESMs with more than three carbon pools, we aggregated τᵢ of cLitter and/or cCwd into that of cFast. The transfer coefficients (r_f and r_s) and the environmental dependency parameter (b) were diagnosed using a Bayesian global optimization algorithm in each grid cell⁵⁴, and their uniform priors are shown in Supplementary Table 8. Specifically, we initialized the reduced complexity model with pool-specific carbon stock from historical ESM simulations (2000–2014), and then used the variables npp, tas and pr from different emissions scenarios as inputs to force the reduced model from 2015 to 2100 for each ESM. The objective function was constructed based on the root mean square error between the actual and modeled soil carbon stock for each of the three carbon pools. The reduced-complexity model with these diagnosed parameters was found to reproduce ESM soil carbon
simulations very well (Supplementary Fig. 6), suggesting that the three-pool model was an almost perfect approximation of the ESMs.

Notably, the five models used in this study may not be representative of the broader CMIP6 ensemble because they draw heavily from the CENTURY model and then have similar structures (e.g., three different soil carbon pools)⁵⁵. Other CMIP models derived from models other than CENTURY may behave differently and deserve further exploration.

### Constrained projections of soil carbon stock changes

The model-data comparison suggested that $τ_i$ in ESMs has a large deviation from interpolated observations at the global scale (Fig. 2). Here, to remove soil carbon simulation biases due to misinterpretations of $τ_i$ in ESMs, we constrained the ESM projections of soil carbon sink by 2100 (Supplementary Fig. 7). Our constructed reduced-complexity model was demonstrated to well mimic soil carbon dynamics in ESMs (Supplementary Fig. 6), but it generally assumes that decomposition rates are only constrained by temperature and moisture availability⁹⁻¹¹. In fact, emerging processes, such as mineral protection and the rhizosphere priming effect, have the potential to affect soil carbon turnover time⁸, but which were mostly absent in current ESMs. We, therefore, refined the reduced-complexity model by including the impacts of climatic factors and these emerging processes (Eqs. (11–12)) on soil carbon turnover time.

$$k = \frac{1}{\tau_i} \times F(T) \times F(P) \times F(M) \times F(RP) \quad (11)$$

Where k is the actual soil carbon decay rate (yr⁻¹) and $τ_i$ denotes the intrinsic carbon turnover time (yr) at a reference temperature of 15 °C. $F(T)$ and $F(P)$ are temperature and precipitation modifiers described by Eq. (10). Specifically, $Q_{10}$ in $F(T)$ was spatially-heterogeneous and pool-specific, rather than a constant as assumed in ESMs⁵⁶. Furthermore, the $Q_{10}$ values of different carbon pools follow the carbon quality-temperature hypothesis³⁸, with higher values for more recalcitrant carbon pools. $F(M)$, a scalar representing the mineral protection of soil carbon, is described as a function of soil clay content (clay in Eq. (12); Supplementary Fig. 17a)⁵⁷.

$$ \begin{cases} F(M) = 24.2 \times \text{clay} \\ F(M) = -2.1 \times \text{clay}^2 + 6.2 \times \text{clay} + 0.6 \end{cases} \quad \begin{array}{l} \text{clay} < 0.033 \\ \text{clay} \ge 0.033 \end{array} \quad (12) $$

$F(RP)$, a scalar representing the plant rhizosphere priming effect, is given as a function of root respiration ($R_{root}$)³².

$$F(RP) = 1 / \left( 1 + \frac{2.47 \times R_{\text{root}}}{13.01 + R_{\text{root}}} \right) \quad (13)$$

This Michaelis–Menten function is derived from a global meta-analysis, which showed a positive relationship between the rhizosphere priming effect and root respiration ($R_{root}$ in Eq. (13)) across all studies (Supplementary Fig. 17b). The root respiration was estimated to be around 7% of NPP³². Here, the rhizosphere priming effect was assumed to only occur in slow and passive pools with “older” soil organic carbon¹⁸,⁵⁸. The rhizosphere priming effect could also increase the release of soil nutrients, which would, in turn, stimulate plant growth³⁹ and thereby create a positive feedback loop that further decreases soil carbon turnover time. On the other hand, enhanced plant growth due to the rhizosphere priming effect might partly offset soil carbon losses due to enhanced soil carbon turnover rates. The net effect due to this positive feedback on soil carbon stock changes might not be so large. This effect was not included in our refined reduced-complexity models due to a general absence of carbon-nitrogen coupling processes in most of ESMs. To fully resolve this question, we require next generation of Earth system models that explicitly

incorporate the rhizosphere priming effect within the coupled carbon-nitrogen cycle framework.

In the refined reduced-complexity model, the parameter $τ_i$ of each pool was taken from our interpolated observations at the grid cell level (Supplementary Fig. 7). Other parameters, including precipitation scalar b, pool-specific temperature scalar $Q_{10}$, as well as the transfer coefficients $r_f$ and $r_s$, were optimized through minimizing errors between the observed and modeled carbon content of the three pools in each grid cell (Supplementary Table 9). These parameters were empirically calibrated so as to ensure that the input rate for each observed soil carbon pool was coupled to $τ_i$ which is modified by climatic factors, mineral protection and rhizosphere priming. The observed soil carbon content for each pool at the grid cell scale was derived from gridded soil carbon data multiplied by our estimated carbon content fraction for each pool (Supplementary Fig. 16). We used data-driven estimates of soil carbon turnover times from global radiocarbon measurements¹⁵ as a surrogate for observed $τ_a$ at the grid scale. In the optimization process, the 15-year mean (2000–2014) of climate and satellite-based NPP were used as inputs. In addition, to account for uncertainties in these empirical parameters, we considered the four combinations of the two gridded soil carbon data sets (SoilGrids and the Harmonized World Soil Database⁶⁰) and the two net primary productivity data sets (MODIS⁶¹ and Global Inventory Modeling and Mapping Studies³⁷) in the optimization algorithm (Supplementary Table 3).

To initialize the soil carbon pool for the start of the future simulations, the refined three-pool model, with parameters calibrated using mean climate and NPP from 2000–2014, was further run repeatedly, using yearly climate and satellite-based NPP from the period 2000–2014 as inputs, for more than 30 cycles until it reached a steady-state condition. The steady-state condition is defined as the 15-year mean difference between NPP and simulated carbon decomposition rate from the three pools being approximately equal to zero. The inclusion of such a steady-state run is to avoid introducing an artefact into the projection of future soil carbon change. After reaching the steady-state condition, the size of the global soil carbon pool reached 2433 ± 325 PgC over the four ensemble predictions: a value that is well within the range of observation-based global soil C stock estimates (–1500–3000 PgC)⁶²,⁶³.

Following the steady-state run, we then ran the refined three-pool model forward to 2100 using bias-corrected climate (temperature and precipitation) and NPP data simulated by CMIP6 ESMs under three different emissions scenarios (Supplementary Fig. 18). We applied the following delta or change method⁶⁴,⁶⁵ to correct biases in the future climate and NPP at the monthly timescale (Eq. (14)).

$$ \begin{cases} T_{\text{fut,cor}} = T_{\text{bas,obs}} + (T_{\text{fut,mod}} - T_{\text{bas,mod}}) \\ P_{\text{fut,cor}} = P_{\text{bas,obs}} \times \left( \frac{P_{\text{fut,mod}}}{P_{\text{bas,mod}}} \right) \\ \text{NPP}_{\text{fut,cor}} = \text{NPP}_{\text{bas,obs}} + (\text{NPP}_{\text{fut,mod}} - \text{NPP}_{\text{bas,mod}}) \end{cases} \quad (14) $$

where $T_{fut,cor}$, $P_{fut,cor}$ and $NPP_{fut,cor}$ are bias-corrected temperature, precipitation and NPP projections for 2015–2100, respectively. $T_{bas,obs}$, $P_{bas,obs}$ and $NPP_{bas,obs}$ are CRU-derived temperature and precipitation, and satellite-derived NPP during the baseline period 2000–2014, respectively. $T_{bas,mod}$, $P_{bas,mod}$ and $NPP_{bas,mod}$ are ESM-simulated temperature, precipitation and NPP for the baseline period 2000–2014, respectively, and $T_{fut,mod}$, $P_{fut,mod}$ and $NPP_{fut,mod}$ are ESM-simulated temperature, precipitation and NPP for 2015–2100, respectively.

**Sensitivity analysis.** To assess the robustness of our constrained results, we conducted the following sensitivity experiments (SE) (Supplementary Fig. 19: Supplementary Table 10). First, we used the Arrhenius function instead of the empirical $Q_{10}$-temperature
relationship in our default simulation to scale site-level $\tau_i$ at their own incubation temperatures to that at 15 °C and then further obtained global estimates of soil $\tau_i$ in a machine learning algorithm linking these $\tau_i$ to environmental variables across sites. We then prescribed $\tau_i$ using these data-driven estimates based on the Arrhenius function in refined reduced-complexity model to constrain soil carbon projections (SE1). Second, the use of laboratory incubation experiments, albeit with the length of the period longer than six months, would still have uncertainties in the quantification of $\tau_i$ of slow-cycling carbon pool especially $C_{passive}$. In SE2, we prescribed $\tau_i$ of $C_{passive}$ as the ensemble mean of ESM's own values rather than the data-driven estimates in our default simulation (Fig. 2).

Third, in SE3-SE6, we used the reduced-complexity model (Eq. 8), which only considered climate controls on soil carbon turnover time and used model inputs such as NPP, MAT and MAP directly from the original ESMs, instead of the refined model (Eq. 11) to constrain soil carbon projections in ESMs. In order to evaluate the relative importance of soil $\tau_i$ and $Q_{10}$ in soil carbon projections, we performed the following SE3-SE5 tests. For SE3, we replaced soil $\tau_i$ with our data-driven estimates, and replaced $Q_{10}$ for each pool with that derived from the refined reduced-complexity model (Eq. 11), in which $Q_{10}$ for each pool was obtained through calibration against observed soil carbon stocks (Supplementary Table 9). For SE4, we replaced soil $\tau_i$ with our data-driven estimates and used ESM's own $Q_{10}$ and associated reference temperatures. For SE5, we used ESM's own soil $\tau_i$, but assigned $Q_{10}$ for each pool to the calibrated one. The magnitude of soil carbon loss in SE4 (19.5 PgC averaged across scenarios) is much closer to that in SE3 (26.3 PgC) than in SE5 (13 PgC), suggesting that soil $\tau_i$ is more important than $Q_{10}$ in determining projections of soil carbon dynamics. In addition, since the use of laboratory incubation experiments would have uncertainties in the quantification of $\tau_i$ of $C_{passive}$, for SE6, we used the ensemble mean of ESM's own soil $\tau_i$ of $C_{passive}$, but assigned $\tau_i$ of $C_{fast}$ and $C_{slow}$, and $Q_{10}$ to be the same with those in SE4. The results of SE6 showed that global soils would be a source of carbon to the atmosphere (10 PgC), albeit at a lower magnitude than SE4.

## The impact of changes in global soil carbon stock on the remaining carbon budget

Changed soil carbon sequestration potential due to observational constraints could affect the remaining carbon budget for limiting global warming below 1.5 °C and 2 °C³⁵. To estimate the global soil carbon sequestration potential under the warming targets, we relied on a strong linear relationship between the change in global air temperature and the observationally-constrained changes in global soil carbon stock between the end of this century and the historical period across the five CMIP6 models under the three emissions scenarios (Supplementary Fig. 11). The change in global air temperature was calculated using the 1850–1900 mean as a baseline, while the calculation of the change in soil carbon storage used the 2020 value as a baseline.

## Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.

## Data availability

The outputs of the Earth system models can be downloaded from the CMIP6 website (https://esgf-node.llnl.gov/projects/cmip6/). The WorldClim and CRU climate data are available at http://www.worldclim.com/version2 and https://crudata.uea.ac.uk/cru/data/hrg/, respectively. Soil physicochemical attributes of the GSDE and SoilGrids data sets can be obtained from http://globalchange.bnu.edu.cn/research/soilw and https://soilgrids.org/, respectively. The soil carbon content of HWSD can be obtained from http://www.fao.org/soils-portal/data-hub/soil-maps-and-databases/harmonized-world-soil-

database-v12/en/. The global NPP databases of MODIS and GIMMS3g are available at http://files.ntsg.umt.edu/data/NTSG Products/MOD17/ and https://wkolby.org/data-code/, respectively. The collected metadata and gridded maps of soil $\tau_i$ and fractions of different carbon pools have been deposited in the Figshare data repository (https://doi.org/10.6084/m9 figshare.19641759.v1)⁶⁶. Source data are provided with this paper.

## Code availability

Data analysis was carried out using R v.4.0.5 and MATLAB R2016a. The code used in this study is available at the Figshare data repository (https://doi.org/10.6084/m9 figshare.19641759.v1)⁶⁶.

## References

1. Bradford, M. A. et al. Managing uncertainty in soil carbon feedbacks to climate change. *Nat. Clim. Change* **6**, 751–758 (2016).

2. Melillo, J. M. et al. Long-term pattern and magnitude of soil carbon feedback to the climate system in a warming world. *Science* **358**, 101–105 (2017).

3. Nottingham, A. T., Meir, P., Velasquez, E. & Turner, B. L. Soil carbon loss by experimental warming in a tropical forest. *Nature* **584**, 234–237 (2020).

4. Crowther, T. W. et al. Quantifying global soil carbon losses in response to warming. *Nature* **540**, 104–108 (2016).

5. Garcia-Palacios, P. et al. Evidence for large microbial-mediated losses of soil carbon under anthropogenic warming. *Nat. Rev. Earth Environ.* **2**, 507–517 (2021).

6. He, Y. et al. Radiocarbon constraints imply reduced carbon uptake by soils during the 21st century. *Science* **353**, 1419–1424 (2016).

7. Luo, Z., Wang, G. & Wang, E. Global subsoil organic carbon turnover times dominantly controlled by soil properties rather than climate. *Nat. Commun.* **10**, 3688 (2019).

8. Davidson, E. A. & Janssens, I. A. Temperature sensitivity of soil carbon decomposition and feedbacks to climate change. *Nature* **440**, 165–173 (2006).

9. Doetterl, S. et al. Soil carbon storage controlled by interactions between geochemistry and climate. *Nat. Geosci.* **8**, 780–783 (2015).

10. Haaf, D., Six, J. & Doetterl, S. Global patterns of geo-ecological controls on the response of soil respiration to warming. *Nat. Clim. Change* **11**, 623–627 (2021).

11. Koven, C. D., Lawrence, D. M. & Riley, W. J. Permafrost carbon-climate feedback is sensitive to deep soil carbon decomposability but not deep soil nitrogen dynamics. *Proc. Natl Acad. Sci. USA* **112**, 3752–3757 (2015).

12. Todd-Brown, K. E. O. et al. Causes of variation in soil carbon simulations from CMIP5 Earth system models and comparison with observations. *Biogeosciences* **10**, 1717–1736 (2013).

13. Jian, S. et al. Multi-year incubation experiments boost confidence in model projections of long-term soil carbon dynamics. *Nat. Commun.* **11**, 1–9 (2020).

14. Varney, R. M. et al. A spatial emergent constraint on the sensitivity of soil carbon turnover to global warming. *Nat. Commun.* **11**, 5544–5544 (2020).

15. Shi, Z. et al. The age distribution of global soil carbon inferred from radiocarbon measurements. *Nat. Geosci.* **13**, 555–559 (2020).

16. Schaedel, C., Luo, Y., Evans, R. D., Fei, S. & Schaeffer, S. M. Separating soil CO₂ efflux into C-pool-specific decay rates via inverse analysis of soil incubation data. *Oecologia* **171**, 721–732 (2013).

17. Schaedel, C. et al. Circumpolar assessment of permafrost C quality and its vulnerability over time using long-term incubation data. *Glob. Change Biol.* **20**, 641–652 (2014).

18. Fontaine, S. et al. Stability of organic carbon in deep soil layers controlled by fresh carbon supply. *Nature* **450**, 277–280 (2007).
19. Kuzyakov, Y., Friedel, J. & Stahr, K. Review of mechanisms and quantification of priming effects. *Soil Biol. Biochem.* **32**, 1485-1498 (2000).

20. Schmidt, M. W. I. et al. Persistence of soil organic matter as an ecosystem property. *Nature* **478**, 49–56 (2011).

21. Hartley, I. P., Hill, T. C., Chadburn, S. E. & Hugelius, G. Temperature effects on carbon storage are controlled by soil stabilisation capacities. *Nat. Commun.* **12**, 1–7 (2021).

22. Ye, J. S., Bradford, M. A., Maestre, F. T., Li, F. M. & Garcia-Palacios, P. Compensatory thermal adaptation of soil microbial respiration rates in global croplands. *Glob. Biogeochem. Cycles* **34**, e2019GB006507 (2020).

23. Hamdi, S., Moyano, F., Sall, S., Bernoux, M. & Chevallier, T. Synthesis analysis of the temperature sensitivity of soil respiration from laboratory studies in relation to incubation methods and soil conditions. *Soil Biol. Biochem.* **58**, 115–126 (2013).

24. Bradford, M. A. Thermal adaptation of decomposer communities in warming soils. *Front. Microbiol.* **4**, 333 (2013).

25. Bradford, M. A. et al. Cross-biome patterns in soil microbial respiration predictable from evolutionary theory on thermal adaptation. *Nat. Ecol. Evol.* **3**, 223–231 (2019).

26. Dacal, M., Bradford, M. A., Plaza, C., Maestre, F. T. & Garcia-Palacios, P. Soil microbial respiration adapts to ambient temperature in global drylands. *Nat. Ecol. Evol.* **3**, 232–238 (2019).

27. Crowther, T. W. & Bradford, M. A. Thermal acclimation in widespread heterotrophic soil microbes. *Ecol. Lett.* **16**, 469–477 (2013).

28. Sáez-Sandino, T. et al. The soil microbiome governs the response of microbial respiration to warming across the globe. *Nat. Clim. Change* **13**, 1382–1387 (2023).

29. Terrer, C. et al. Nitrogen and phosphorus constrain the CO₂ fertilization of global plant biomass. *Nat. Clim. Change* **9**, 684–689 (2019).

30. Schuur, E. A. G. et al. Climate change and the permafrost carbon feedback. *Nature* **520**, 171–179 (2015).

31. Turetsky, M. R. et al. Permafrost collapse is accelerating carbon release. *Nature* **569**, 32–34 (2019).

32. Keuper, F. et al. Carbon loss from northern circumpolar permafrost soils amplified by rhizosphere priming. *Nat. Geosci.* **13**, 560–565 (2020).

33. Nottingham, A. T. et al. Microbial responses to warming enhance soil carbon loss following translocation across a tropical forest elevation gradient. *Ecol. Lett.* **22**, 1889–1899 (2019).

34. Winkler, A. J. et al. Slowdown of the greening trend in natural vegetation with further rise in atmospheric CO₂. *Biogeosciences* **18**, 4985–5010 (2021).

35. Lamboll, R. D. et al. Assessing the size and uncertainty of remaining carbon budgets. *Nat. Clim. Change* **13**, 1–8 (2023).

36. Wang, M. et al. Global soil profiles indicate depth-dependent soil carbon losses under a warmer climate. *Nat. Commun.* **13**, 5514 (2022).

37. Smith, W. K. et al. Large divergence of satellite and Earth system model estimates of global terrestrial CO₂ fertilization. *Nat. Clim. Change* **6**, 306–310 (2016).

38. Craine, J. M., Fierer, N. & McLauchlan, K. K. Widespread coupling between the rate and temperature sensitivity of organic matter decay. *Nat. Geosci.* **3**, 854–857 (2010).

39. Knorr, W., Prentice, I. C., House, J. I. & Holland, E. A. Long-term sensitivity of soil carbon turnover to warming. *Nature* **433**, 298–301 (2005).

40. Schaedel, C. et al. Decomposability of soil organic matter over time: the Soil Incubation Database (SIDb, version 1.0) and guidance for incubation procedures. *Earth Syst. Sci. Data* **12**, 1511–1524 (2020).

41. Elberling, B. et al. Long-term CO₂ production following permafrost thaw. *Nat. Clim. Change* **3**, 890–894 (2013).

42. von Luetzow, M. & Koegel-Knabner, I. Temperature sensitivity of soil organic matter decomposition-what do we know? *Biol. Fertil. Soils* **46**, 1–15 (2009).

43. Fick, S. E. & Hijmans, R. J. WorldClim 2: new 1-km spatial resolution climate surfaces for global land areas. *Int. J. Climatol.* **37**, 4302–4315 (2017).

44. Hengl, T. et al. SoilGrids250m: Global gridded soil information based on machine learning. *PLoS One* **12**, e0169748 (2017).

45. Xu, T., White, L., Hui, D. F. & Luo, Y. Q. Probabilistic inversion of a terrestrial ecosystem model: Analysis of uncertainty in parameter estimation and model prediction. *Glob. Biogeochem. Cycles* **20**, 1–15 (2006).

46. Tarantola, A. & Valette, B. Generalized nonlinear inverse problems solved using the least squares criterion. *Rev. Geophys.* **20**, 219–232 (1982).

47. Metropolis, N., Rosenbluth, A. W., Rosenbluth, M. N., Teller, A. H. & Teller, E. Equation of state calculations by fast computing machines. *J. Chem. Phys.* **21**, 1087–1092 (1953).

48. Hastings, W. K. Monte Carlo sampling methods using Markov chains and their applications. *Biometrika* **57**, 97–109 (1970).

49. Xu, X. et al. Global pattern and controls of soil microbial metabolic quotient. *Ecol. Monogr.* **87**, 429–441 (2017).

50. Moyano, F. E. et al. The moisture response of soil heterotrophic respiration: interaction with soil properties. *Biogeosciences* **9**, 1173–1182 (2012).

51. Elith, J., Leathwick, J. R. & Hastie, T. A working guide to boosted regression trees. *J. Anim. Ecol.* **77**, 802–813 (2008).

52. Harris, I., Osborn, T. J., Jones, P. & Lister, D. Version 4 of the CRU TS monthly high-resolution gridded multivariate climate dataset. *Sci. Data* **7**, 1–18 (2020).

53. Shangguan, W., Dai, Y., Duan, Q., Liu, B. & Yuan, H. A global soil data set for earth system modeling. *J. Adv. Model. Earth Syst.* **6**, 249–263 (2014).

54. Duan, Q. Y., Sorooshian, S. & Gupta, V. K. Optimal use of the SCE-UA global optimization method for calibrating watershed models. *J. Hydrol.* **158**, 265–284 (1994).

55. Berardi, D. et al. 21st-century biogeochemical modeling: challenges for Century-based models and where do we go from here? *GCB Bioenergy* **12**, 774–788 (2020).

56. Zhou, T., Shi, P., Hui, D. & Luo, Y. Global pattern of temperature sensitivity of soil heterotrophic respiration (Q₁₀) and its implications for carbon-climate feedback. *J. Geophys. Res. Biogeoosci.* **114**, G02016 (2009).

57. Sulman, B. N., Phillips, R. P., Oishi, A. C., Shevliakova, E. & Pacala, S. W. Microbe-driven turnover offsets mineral-mediated storage of soil carbon under elevated CO₂. *Nat. Clim. Change* **4**, 1099–1102 (2014).

58. Chen, L. et al. Regulation of priming effect by soil organic matter stability over a broad geographic scale. *Nat. Commun.* **10**, 1–10 (2019).

59. Drake, J. E. et al. Stoichiometry constrains microbial response to root exudation-insights from a model and a field experiment in a temperate forest. *Biogeosciences* **10**, 821–838 (2013).

60. Wieder, W. R., Boehnert, J., Bonan, G. B. & Langseth, M. Reggridded Harmonized World Soil Database v1.2. ORNL DAAC. https://doi.org/10.3334/ORNLDAAC/1247 (2014).

61. Zhao, M. & Running, S. W. Drought-induced reduction in global terrestrial net primary production from 2000 through 2009. *Science* **329**, 940–943 (2010).

62. Scharlemann, J. P., Tanner, E. V., Hiederer, R. & Kapos, V. Global soil carbon: understanding and managing the largest terrestrial carbon pool. *Carbon Manag.* **5**, 81–91 (2014).

63. Tifafi, M., Guenet, B. & Hatté, C. Large differences in global and regional total soil carbon stock estimates based on SoilGrids, HWSD, and NCSCD: Intercomparison and evaluation based on field data from USA, England, Wales, and France. *Glob. Biogeochem. Cycles* **32**, 42–56 (2018).
64. Wang, T. et al. Data-driven mapping of the spatial distribution and potential changes of frozen ground over the Tibetan Plateau. *Sci. Total Environ.* **649**, 515–525 (2019).

65. Wang, T. et al. Permafrost thawing puts the frozen carbon at risk over the Tibetan Plateau. *Sci. Adv.* **6**, eaaz3513 (2020).

66. Ren, S. et al. Data and code for 'Projected soil carbon loss with warming in constrained Earth system models'. Figshare. https://doi.org/10.6084/m9 figshare.19641759.v1 (2023).

## Acknowledgements

This study was supported by the National Natural Science Foundation of China (42230411), the NSFC project Basic Science Centre for Tibetan Plateau Earth System (41988101), the Second Tibetan Plateau Scientific Expedition and Research Program (2019QZKK0606), the Science and Technology Plan Project of Tibet Autonomous Region (XZ202201ZY0015G), and Innovation Program for Young Scholars of TPESER (TPESER-QNCX2022ZD-02). We also acknowledge the support of Kathmandu Center for Research and Education, Chinese Academy of Sciences—Tribhuvan University.

## Author contributions

T.W. designed the research; T.W. and S.R. wrote the paper; S.R. performed the data analysis; all authors contributed to the interpretation of the results and to the text.

## Competing interests

The authors declare no competing interests.

## Additional information

**Supplementary information** The online version contains supplementary material available at https://doi.org/10.1038/s41467-023-44433-2.

**Correspondence** and requests for materials should be addressed to Tao Wang.

**Peer review information** *Nature Communications* thanks David Armstrong McKay and Daniel Ricciuto for their contribution to the peer review of this work. A peer review file is available.

**Reprints and permissions information** is available at http://www.nature.com/reprints

**Publisher's note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.

© The Author(s) 2024