Ren et al. (2024). Projected soil carbon loss with warming in constrained Earth system models

, Pete Smith

4

& Shilong Piao

1,5

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 dif

ﬁ

cult.

Soil is the largest reservoir of terrestrial organic carbon

, 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

4

5

, either in terms of the magnitude of the

effect or even its sign. One signi

ﬁ

cant, and poorly understood, com-

ponent of the system is soil carbon turnover

, which is de

ﬁ

ned as the

average time it takes for a carbon atom to enter and leave the soil
system

7

. 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 con-
stants, which are modi

ﬁ

ed only by abiotic factors (temperature and

moisture) with few spatially-uniform parameters

9

–

. The models vary

greatly in the assumed values of constants and parameters

, and the

observational constraints on these values are generally lacking

,

4

,

13

. This

de

ﬁ

ciency 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

14

and from radiocarbon-based

observations

,

15

, reducing con

ﬁ

dence 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 (

τ

i

) across the globe

16

,

. Notably, soil

τ

i

is representative of the theoretical carbon turnover time under opti-
mal 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

τ

i

in the real-

world settings

(Methods). Using global soil

τ

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)

18

–

21

, to

Received: 31 March 2023

Accepted: 13 December 2023

Check for updates

1

State Key Laboratory of Tibetan Plateau Earth System, Environment and Resources (TPESER), Institute of Tibetan Plateau Research, Chinese Academy of

Sciences, Beijing, China.

2

University of Chinese Academy of Sciences, Beijing, China.

3

Laboratoire de Géologie, École normale supérieure, CNRS, PSL

University, IPSL, Paris, France.

4

Institute of Biological and Environmental Sciences, School of Biological Sciences, University of Aberdeen, Aberdeen AB24

3UU, UK.

5

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

Nature Communications

| (2024) 15:102

1

12

3

4

5

6

78

9

0

()

:,;

12

3

4

5

6

7

8

9

0

()

:,;

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

2

) ef

ﬂ

ux data from

temperature-controlled soil incubation studies conducted at 178 sites
covering eight biomes, ranging from Arctic permafrost to dry Medi-
terranean forests (Fig.

a). Due to the large heterogeneity of soil car-

bon, we represent the time evolution of soil CO

2

ef

ﬂ

ux (

SR

) using a

three-pool model which partitions the total soil carbon (

C

tot

) into three

reservoirs, with each decaying at its own intrinsic turnover rate (that is,
the inverse of turnover time

τ

)

16

,

.

SR

=

X

3

p = 1

1

τ

p

C

tot

f

p

ð

1

Þ

where

f

is the partitioning coef

ﬁ

cient of

C

tot

for each pool (

p

). Using

deconvolution analysis, the pool-speci

ﬁ

c

τ

values were quanti

ﬁ

ed 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

23

, suggesting that the inverted carbon pool-speci

ﬁ

c

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

τ

i

within each biome and

across biomes for each of the three carbon reservoirs (Fig.

b

–

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

τ

i

variations.

Of the tested predictors, mean annual temperature (MAT) was the
most important variable in explaining cross-site variability of

τ

i

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

τ

i

variability in the fast soil carbon pool (

C

fast

) and in the

slow (

C

slow

) and passive (

C

passive

) carbon pools. Speci

ﬁ

cally, soil

τ

i

had

positive correlations with MAT (Supplementary Figs. 2

–

3), with lower

MAT being associated with shorter

τ

i

.

Based on the empirical relationships between environmental

predictors and soil

τ

i

, we generated a predictive model that could

explain more than 75% of the cross-site variability in

τ

i

of each soil

carbon pool without bias (Supplementary Fig. 1). This predictive model
enabled us to interpolate

τ

i

of each soil carbon pool to a depth of 1 m

across global soils (Fig.

a, d, g). We derived a global

τ

i

of 0.3 yr for

C

fast

,

with a spatial variation ranging between 0.06 yr (

ﬁ

rst percentile) and

0.64 yr (99th percentile) (Fig.

a and Supplementary Table 1). The

global

C

slow

and

C

passive

τ

i

are 6.68 (2.3

–

12.1) and 398 (89

–

696) yr,

respectively, with the values in parentheses denoting the spatial range
between the

ﬁ

rst and 99

th

percentiles (Fig.

d, g and Supplementary

Table 1). We also calculated carbon-weighted

τ

i

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 speci

ﬁ

c pools

and environmental drivers across sites (see Methods). The global
carbon-weighted

τ

i

is 316 yr with a spatial variation ranging between

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 (

τ

i

) of

C

fast

(

b

),

C

slow

(

c

) and

C

passive

(

d

) that are inverted from all

experiments within each of the eight biomes, respectively. The eight biomes are

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.

Article

https://doi.org/10.1038/s41467-023-44433-2

Nature Communications

| (2024) 15:102

2

81 yr (

ﬁ

rst percentile) and 609 yr (99th percentile) (Supplementary

Table 1).

The longest

τ

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

ﬁ

nd that the

τ

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

τ

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.

a, 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

ﬁ

nd that soil

τ

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

τ

i

speci

ﬁ

ed in the

ﬁ

ve 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

τ

i

in ESMs are

generally described as global constants without any spatial variability
(Supplementary Table 2). We found that ESMs overestimated soil

τ

i

in

high-latitude ecosystems for all carbon pools (Fig.

and Supplemen-

tary 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

τ

i

are overestimated to a

less extent, and even underestimated in some regions (Supplementary
Fig. 5). Overall, the model ensemble overestimated soil

τ

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.

). The model-data

bias of soil

τ

i

could be attributed to the omission of critical microbial

processes from ESMs, such as thermal adaptation

24

. Speci

ﬁ

cally,

microbial turnover rates have been shown to adjust to temperature
changes via biochemical trade-offs in enzyme and cell membrane
structure and function

25

26

. Low temperatures typically select for

Fig. 2 | Global distributions of intrinsic soil carbon turnover times.
a

,

d

,

g

Observation-derived intrinsic soil carbon turnover times (

τ

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

τ

i

estimates derived from different combinations of data sets.

c

,

f

,

i

Latitudinal

pro

ﬁ

les of soil

τ

i

, aggregated at 0.5° latitudinal resolution. The black and red lines

indicate observation-derived and ESMs-averaged soil

τ

i

, respectively. The pink

dashed lines are soil

τ

i

speci

ﬁ

ed in each model (Supplementary Table 2). The red

shading represents the standard deviation. The x-axis data in

c

,

f

,

i

are log-

transformed.

Article

https://doi.org/10.1038/s41467-023-44433-2

Nature Communications

| (2024) 15:102

3

$Figures and page graphics, page 4$

enzymes and/or membranes that are highly

ﬂ

exible to ef

ﬁ

ciently alter

conformation and facilitate interactions

27

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

ing microbial metrics (such as thermal adaptation

24

, species

composition

28

) into ESMs to build con

ﬁ

dence in predicting soil carbon-

climate feedback.

The large deviation of modeled

τ

i

from the observations could

lead to substantial biases in projections of soil carbon stock changes.
To remove simulation biases due to misinterpretations of

τ

i

in ESMs,

we

ﬁ

rstly 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 re

ﬁ

ned reduced-complexity model to con-

strain ESM projections under three different emissions scenarios
(SSP1-2.6, SSP2-4.5 and SSP5-8.5; Supplementary Fig. 7). In the re

ﬁ

ned

reduced-complexity model, we described the soil carbon decay as a
function of

τ

i

that is not only modi

ﬁ

ed by climate factors, but also by

the two competing processes (mineral protection and the rhizosphere
priming effect) (see Methods). We forced the model with climatolo-
gical mean (2000

–

2014) climate variables and satellite-based net pri-

mary productivity (NPP), and calibrated the model against both the
gridded historical soil carbon pool and upscaled

ﬁ

eld determinations

of soil carbon age (as a proxy of apparent soil carbon turnover times)
derived from radiocarbon measurements

15

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 Map-
ping Studies) (Supplementary Table 3 and Supplementary Fig. 8).

The re

ﬁ

ned 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 sce-
narios (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.

and Supplementary Fig. 9). In contrast, constrained

ESM simulations using the re

ﬁ

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

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

ﬁ

nding is gen-

erally consistent with the expectation that the intrinsic turnover time
was overestimated in complex ESMs (Fig.

). 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 con-
straints 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 inter-
mediate 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

29

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

).

Such a reversal of the soil carbon balance is most prominent in the
arctic tundra and boreal regions (Supplementary Fig. 10), and occur-
red there even without invoking warming-induced permafrost soil
carbon degradation

30

31

. The arctic tundra and boreal region have the

shortest soil-carbon intrinsic turnover time of all the regions (Fig.

),

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

32

. Fur-

thermore, 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

,

33

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

dence suggests that this CO

2

fertilization effect in ESMs may have been

overestimated

34

, 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

35

, respectively. We further esti-

mated 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 car-
bon budget should be reduced by nearly 15

–

66% over the course of

this century to achieve the warming goals, rendering climate mitiga-
tion much more dif

ﬁ

cult.

In summary, we synthesized long-term incubated-soil CO

2

ﬂ

ux

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

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.

Article

https://doi.org/10.1038/s41467-023-44433-2

Nature Communications

| (2024) 15:102

4

feedback to climate change in a warmer world

, 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 para-
meterize 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 car-
bon dynamics

, 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 re

ﬂ

ects kinetic properties of

various soil organic compounds under optimal conditions

. 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 freez-
ing,

ﬂ

ooding and physical protection, can dampen decomposition

processes

, 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

2

ﬂ

uxes under controlled conditions

,

39

. Results from such soil incuba-

tion experiments are invaluable for informing process-based ecosys-
tem 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

2

ﬂ

uxes 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

Here, we assembled a raw dataset of soil CO

2

ﬂ

uxes from soil

incubation experiments to assess the intrinsic decomposability of
different soil carbon pools. We included only soil incubations which
satis

ﬁ

ed the following criteria: (1) the length of the incubation

experiment was longer than 6 months; (2) the aerobic CO

2

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 incuba-
tion 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

; (3) the initial

soil carbon concentration must be available; (4) the soils were incu-
bated 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

ﬁ

gures using

the GetData (v.2.25) software. In total, we gathered 102 peer-reviewed
publications covering 178 sites (Fig.

), resulting in a total of >5000

time-series of soil CO

2

ﬂ

ux data. The soil CO

2

ﬂ

ux 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 pre-
cipitation), 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 classi

ﬁ

ed

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.

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

2

ﬂ

ux is composed of contributions from

the three different carbon pools (that is

C

fast

,

C

slow

and

C

passive

). For

each sampling site, we

ﬁ

tted the time series of total soil CO

2

ﬂ

uxes

using the three-pool carbon decomposition model

described below.

In this model, the total soil CO

2

ﬂ

ux (

R

, in mg C g

-1

dw d

-1

) was the sum of

the respiration rates (

r

) derived from the three carbon pools (

p

) with

different sizes and turnover rates (Eq. (

)).

R

=

X

3

p = 1

r

p

=

X

3

p = 1

k

p

×

C

tot

×

f

p

ð

2

Þ

f

p

=

C

p

C

tot

,

X

3

p = 1

f

p

= 1

ð

3

Þ

where the pool-speci

ﬁ

c respiration rate (

r

p

) was computed as the

product of intrinsic pool-speci

ﬁ

c decay rate (

k

p

, in d

−

1

; the inverse of

τ

i

),

the total initial carbon pool (

C

tot

, in mg C g

−

1

dw) and a partitioning

coef

ﬁ

cient (

f

p

). The partitioning coef

ﬁ

cient describes the ratio of

carbon pool

p

to the total carbon pool, with the sum of the three

coef

ﬁ

cients equal to unity (Eq. (

)).

The model parameters (

m

), consisting of

k

and

f

for each pool,

were optimized against measured soil CO

2

ﬂ

uxes (

O

) using a Bayesian

probabilistic inversion approach

45

. This approach states that the pos-

terior 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

j

O

Þ /

p

ð

O

j

m

Þ

p

ð

m

Þ

ð

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

ﬁ

rst

speci

ﬁ

ed 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 mea-

sured CO

2

ﬂ

ux data. Through assuming that the errors of measured

CO

2

ﬂ

ux data follow a normal distribution with zero mean, the like-

lihood function

p

(

O

|

m

) is given by,

p

ð

O

j

m

Þ /

exp

1

2

σ

2

X

t

2

obs

ð

S

i

Þ

S

i

ð

t

Þ

O

i

ð

t

Þ

2

8

<
:

9

=
;

ð

5

Þ

where O

i

(t) and S

i

(t) are the measured and simulated CO

2

ﬂ

uxes, and

σ

is the standard deviation of the measured CO

2

ﬂ

uxes.

To derive the posterior distribution of model parameters analy-

tically, 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

Article

https://doi.org/10.1038/s41467-023-44433-2

Nature Communications

| (2024) 15:102

5

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 con

ﬁ

dence intervals for all parameters

from their posterior distributions. Overall, the performance of the
three-pool model in

ﬁ

tting the measured soil CO

2

ﬂ

uxes was good

(Supplementary Figs. 12 and 13).

In the compiled data, soil incubation temperatures ranged from 4

to 35 °C, making it dif

ﬁ

cult to compare

τ

i

across studies in a quanti-

tative manner. Therefore, we adjusted

τ

i

at different incubation tem-

peratures to a reference temperature of 15 °C by using the following
equation

22

49

.

τ

i, 15

=

τ

i,T

× Q

T

15

10

ð

6

Þ

where

τ

i,15

and

τ

i,T

are soil intrinsic turnover times at the reference

temperature of 15 °C and the incubation temperature T (°C), respec-
tively. Q

10

is a temperature sensitivity parameter de

ﬁ

ned as the factor

by which soil respiration increases with a 10 °C increase in tempera-
ture. Here, we assumed that Q

10

varies with temperature, and their

empirical function was derived from a previous synthesis analysis of
Q

10

and temperature from laboratory studies across various ecosys-

tems at the global scale

23

(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 impor-
tance 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

τ

i

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

50

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

51

. Our approach,

which is superior to most traditional methods, can handle different
types of predictor variables and interaction effects between pre-
dictors, and does not require data transformation or outlier elimina-
tion. 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

τ

i

data were log10 transformed before

starting the analysis.

The constructed BRT model was able to explain more than 75% of

the variances in

τ

i

across different sites (Supplementary Fig. 1), and was

subsequently used to produce global maps of

τ

i

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 com-
binations of different sources of climate and edaphic data (Fig.

). For

climate predictors, we used mean annual temperature and mean
annual precipitation from both the Climate Research Unit (CRU) ver-
sion 4.01

52

and WorldClim version 2.0

43

; for edaphic factors (soil phy-

sicochemical properties), we considered both the Global Soil Dataset
for Earth System Modeling (GSDE)

53

and the SoilGrids

44

data (Supple-

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

). Thus, more long-term soil

incubations are urgently needed in these speci

ﬁ

c regions.

To calculate soil carbon-weighted

τ

i

(

τ

i,w

in Eq. (

7

)), we used the

same method as used to obtain global maps of fractions in different
soil carbon pools (

f

p

in Eq. (

7

); Supplementary Fig. 16).

τ

i,w

=

X

3

p = 1

τ

i,p

×

f

p

ð

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

ﬁ

ve ESMs.

This three-pool model partitions soil organic carbon into components
with different intrinsic turnover rates (Eqs. (

–

10

)).

d

C

f

dt

=

NPP

ð

t

Þ

k

f

×

C

f

ð

t

Þ

d

C

s

dt

=

k

f

×

C

f

ð

t

Þ

×

r

f

k

s

×

C

s

ð

t

Þ

d

C

p

dt

=

k

s

×

C

s

ð

t

Þ

×

r

s

k

p

×

C

p

ð

t

Þ

8

>

<
>

>

:

ð

8

Þ

where

r

f

and

r

s

are transfer coef

ﬁ

cients for carbon

ﬂ

owing 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

-1

) cal-

culated using Eq. (

9

).

k

=

1

τ

a

=

1

τ

i

×

F

ð

T

Þ

×

F

ð

P

Þ

ð

9

Þ

F

ð

T

Þ

= Q

T

Tref

10

;

F

ð

P

Þ

=

P

b

ð

10

Þ

Where

τ

a

denotes the actual-, and

τ

i

the intrinsic, carbon turnover time

(yr) at a reference temperature (

T

ref

).

T

and

P

are the ambient tem-

perature (°C) and precipitation (mm), respectively.

F

(

T

) is a Q

10

-based

standard exponential function to represent the temperature modi

ﬁ

er

of

τ

i

.

F

(

P

) represents a moisture modi

ﬁ

er 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,

ﬁ

ve of eight parameters, including

intrinsic soil carbon turnover time of each pool, Q

10

and associated

T

ref

were from original ESMs (Supplementary Table 2). In particular, for
ESMs with more than three carbon pools, we aggregated

τ

i

of

cLitte

and/or

cCwd

into that of

Cfast

. The transfer coef

ﬁ

cients (

r

f

and

r

s

) and

the environmental dependency parameter (

b

) were diagnosed using a

Bayesian global optimization algorithm in each grid cell

54

, and their

uniform priors are shown in Supplementary Table 8. Speci

ﬁ

cally, we

initialized the reduced complexity model with pool-speci

ﬁ

c 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 diag-
nosed parameters was found to reproduce ESM soil carbon

Article

https://doi.org/10.1038/s41467-023-44433-2

Nature Communications

| (2024) 15:102

6

simulations very well (Supplementary Fig. 6), suggesting that the
three-pool model was an almost perfect approximation of the ESMs.

Notably, the

ﬁ

ve models used in this study may not be repre-

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

55

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

).

Here, to remove soil carbon simulation biases due to misinterpreta-
tions 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 moist-
ure availability

9

–

. In fact, emerging processes, such as mineral pro-

tection and the rhizosphere priming effect, have the potential to affect
soil carbon turnover time

, but which were mostly absent in current

ESMs. We, therefore, re

ﬁ

ned the reduced-complexity model by

including the impacts of climatic factors and these emerging processes
(Eqs. (

–

)) on soil carbon turnover time.

k

=

1

τ

i

×

F

ð

T

Þ

×

F

ð

P

Þ

×

F

ð

M

Þ

×

F

ð

RP

Þ

ð

11

Þ

Where

k

is the actual soil carbon decay rate (yr

-1

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

ﬁ

ers described by

Eq. (

10

). Speci

ﬁ

cally, Q

10

in

F

(

T

) was spatially-heterogeneous and pool-

speci

ﬁ

c, rather than a constant as assumed in ESMs

56

. Furthermore, the

Q

10

values of different carbon pools follow the carbon quality-

temperature hypothesis

38

, 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. (

);

Supplementary Fig. 17a)

57

.

F

ð

M

Þ

= 24

:

2 ×

clay

< 0

:

033

F

ð

M

Þ

=

2

:

1 ×

clay

2

+ 6

:

2 ×

clay

+ 0

:

6

clay

≥

0

:

033

ð

12

Þ

F

(

RP

), a scalar representing the plant rhizosphere priming effect,

is given as a function of root respiration (

R

root

)

32

.

F

ð

RP

Þ

= 1

.

1 +

2

:

47 ×

R

root

13

:

01 +

R

root

ð

13

Þ

This Michaelis

–

Menten function is derived from a global meta-

analysis, which showed a positive relationship between the rhizo-
sphere 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

32

. Here, the rhizosphere priming effect was

assumed to only occur in slow and passive pools with

“

older

”

soil

organic carbon

18

58

. The rhizosphere priming effect could also increase

the release of soil nutrients, which would, in turn, stimulate plant
growth

59

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 re

ﬁ

ned reduced-

complexity models due to a general absence of carbon-nitrogen cou-
pling 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 re

ﬁ

ned 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-speci

ﬁ

c temperature scalar Q

10

, as well as the transfer

coef

ﬁ

cients

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 modi

ﬁ

ed 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

15

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 con-
sidered the four combinations of the two gridded soil carbon data sets
(SoilGrids and the Harmonized World Soil Database

60

) and the two net

primary productivity data sets (MODIS

61

and Global Inventory Model-

ing and Mapping Studies

37

) in the optimization algorithm (Supple-

mentary Table 3).

To initialize the soil carbon pool for the start of the future simu-

lations, the re

ﬁ

ned 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 de

ﬁ

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

62

63

.

Following the steady-state run, we then ran the re

ﬁ

ned 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

64

65

to correct biases in the future

climate and NPP at the monthly timescale (Eq. (

14

)).

T

fut,cor

=

T

bas,obs

+

ð

T

fut,mod

T

bas,mod

Þ

P

fut,cor

=

P

bas,obs

×

P

fut,mod

P

bas,mod

NPP

fut,cor

=

NPP

bas,obs

+

ð

NPP

fut,mod

NPP

bas,mod

Þ

8

>

<
>

>

:

ð

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

Article

https://doi.org/10.1038/s41467-023-44433-2

Nature Communications

| (2024) 15:102

7

relationship in our default simulation to scale site-level

τ

i

at their own

incubation temperatures to that at 15 °C and then further obtained
global estimates of soil

τ

i

in a machine learning algorithm linking these

τ

i

to environmental variables across sites. We then prescribed

τ

i

using

these data-driven estimates based on the Arrhenius function in re

ﬁ

ned

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 uncer-
tainties in the quanti

ﬁ

cation of

τ

i

of slow-cycling carbon pool especially

C

passive

. In SE2, we prescribed

τ

i

of

C

passive

as the ensemble mean of

ESM

’

s own values rather than the data-driven estimates in our default

simulation (Fig.

).

Third, in SE3-SE6, we used the reduced-complexity model (Eq.

),

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 re

ﬁ

ned model (Eq.

) to constrain soil

carbon projections in ESMs. In order to evaluate the relative impor-
tance of soil

τ

i

and Q

10

in soil carbon projections, we performed the

following SE3-SE5 tests. For SE3, we replaced soil

τ

i

with our data-

driven estimates, and replaced Q

10

for each pool with that derived

from the re

ﬁ

ned reduced-complexity model (Eq.

https://esgf-node.llnl.gov/projects/cmip6/

), in which Q

10

for

each pool was obtained through calibration against observed soil
carbon stocks (Supplementary Table 9). For SE4, we replaced soil

τ

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

τ

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

τ

i

is

more important than Q

10

in determining projections of soil carbon

dynamics. In addition, since the use of laboratory incubation experi-
ments would have uncertainties in the quanti

ﬁ

cation of

τ

i

of

C

passive

, for

SE6, we used the ensemble mean of ESM

’

s own soil

τ

i

of

C

passive

, but

assigned

τ

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

35

. 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 tem-
perature and the observationally-constrained changes in global soil
carbon stock between the end of this century and the historical period
across the

ﬁ

ve 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 calcu-

lation 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

(

).

The

WorldClim and CRU climate data are available at

http://www.

worldclim.com/version2

and

https://crudata.uea.ac.uk/cru/data/hrg/