Varney et al. (2020). A spatial emergent constraint on the sensitivity of soil carbon turnover to global warming

✉

, Sarah E. Chadburn

1

, Pierre Friedlingstein

1,2

, Eleanor J. Burke

3

,

Charles D. Koven

4

, Gustaf Hugelius

5

& Peter M. Cox

1

Carbon cycle feedbacks represent large uncertainties in climate change projections, and the

response of soil carbon to climate change contributes the greatest uncertainty to this. Future

changes in soil carbon depend on changes in litter and root inputs from plants and especially

on reductions in the turnover time of soil carbon (

τ

s

) with warming. An approximation to the

latter term for the top one metre of soil (

Δ

C

s,

τ

) can be diagnosed from projections made with

the CMIP6 and CMIP5 Earth System Models (ESMs), and is found to span a large range even

at 2 °C of global warming (

−

196 ± 117 PgC). Here, we present a constraint on

Δ

C

s,

τ

, which

makes use of current heterotrophic respiration and the spatial variability of

τ

s

inferred from

observations. This spatial emergent constraint allows us to halve the uncertainty in

Δ

C

s,

τ

at

2 °C to

−

232 ± 52 PgC.

https://doi.org/10.1038/s41467-020-19208-8

OPEN

1

College of Engineering, Mathematics and Physical Sciences, University of Exeter, Laver Building, North Park Road, Exeter EX4 4QF, UK.

2

LMD/IPSL, ENS,

PSL Université, École Polytechnique, Institut Polytechnique de Paris, Sorbonne Université, CNRS, 75006 Paris, France.

3

Met Of

ﬁ

ce Hadley Centre, FitzRoy

Road, Exeter EX1 3PB, UK.

4

Earth and Environmental Sciences Division, Lawrence Berkeley National Laboratory, Berkeley, CA 94720, USA.

5

Department of

Physical Geography and Bolin Centre of Climate Research, Stockholm University, Stockholm 10691, Sweden.

✉

email:

r.varney@exeter.ac.uk

NATURE COMMUNICATIONS

| (2020) 11:5544 | https://doi.org/10.1038/s41467-020-19208-8 | www.nature.com/naturecommunications

123456789

0():,;

C

limate

–

carbon cycle feedbacks

must be understood and

quanti

ﬁ

ed if the Paris Agreement Targets are to be met

Changes in soil carbon represent a particularly large

, with the potential to signi

ﬁ

cantly reduce the

carbon budget for climate stabilisation at 2 °C global warming

8

Previous studies have investigated the response of soil carbon to
climate change based on both observational studies

9

and Earth

System Models (ESMs)

10

. ESMs are coupled models which

simulate both climate and carbon cycle processes. Projects such as
the Coupled Model Inter-comparison Project (CMIP)

11

,

12

, have

allowed for consistent comparison of the response of soil carbon
under climate change from existing state-of-the-art ESMs.
However, the uncertainty due to the soil carbon feedback did not
reduce signi

ﬁ

cantly between the CMIP3 and CMIP5 model

generations

6

, or with the latest CMIP6 models (see Fig.

and Supplementary Fig. 1), such that the projected change in
global soil carbon still varies signi

ﬁ

cantly amongst models

13

This study uses an alternative method to obtain a constraint on

the ESM projections of soil carbon change. In previous studies,
emergent constraints based on temporal trends and variations
have been used successfully to reduce uncertainty in climate
change projections

14

. Our approach follows the method used in

Chadburn et al.

15

, where a spatial temperature sensitivity is used

to constrain the future response to climate change

—

which we

term as a

spatial emergent constraint

. Our study combines the

Chadburn et al.

15

method with the soil carbon turnover analysis

of Koven et al.

to get a constraint on the sensitivity of soil

carbon turnover to global warming.

Soil carbon (

C

s

) is increased by the

ﬂ

ux of organic carbon into

the soil from plant litter and roots, and decreased by the break-
down of that organic matter by soil microbes which releases CO

2

to the atmosphere as the heterotrophic respiration

ﬂ

ux (

R

h

). If the

vegetation carbon is at steady-state, litter-fall will equal the Net
Primary Production of plants (NPP). If the soil carbon is also
near to a steady-state

—

and in the absence of signi

ﬁ

cant

ﬁ

re

ﬂ

uxes

and other non-respiratory carbon losses

—

the litter-fall, NPP, and

R

h

will be approximately equal to one another. Even over the

historical period, when atmospheric CO

2

has been increasing and

there has been a net land carbon sink, this approximation holds
well (see Supplementary Fig. 4).

In order to separate the effects of changes in NPP from the

effects of climate change on

R

h

, we de

ﬁ

ne an effective turnover

time

17

for soil carbon as

τ

s

=

C

s

/

R

h

. The turnover time of soil

carbon is known to be especially dependent on temperature

. A

common assumption is that

τ

s

decreases by about 7% per °C of

warming (equivalent to assuming that q

10

=

2)

18

. However, this

sensitivity differs between models, and also between models and
observations.

We can write a long-term change in soil carbon (

Δ

C

s

), as the

sum of a term arising from changes in litter-fall (

Δ

C

s,L

), and a

term arising from changes in the turnover time of soil carbon
(

Δ

C

s,

τ

):

Δ

C

s

¼

Δ

ð

R

h

τ

s

Þ

Δ

C

s

;

L

ð

t

Þ þ

Δ

C

s

;

τ

ð

t

Þ ¼

τ

s

;

0

Δ

R

h

ð

t

Þ þ

R

h

;

0

Δ

τ

s

ð

t

Þ

ð

1

Þ

Model projections of the

ﬁ

rst term (

Δ

C

s,L

) differ primarily

because of differences in the extent of CO

2

-fertilisation of NPP,

and associated nutrient limitations. The second term (

Δ

C

s,

τ

)

differs across models because of differences in the predicted
future warming, and because of differences in the sensitivity of
soil carbon decomposition to temperature (which includes an
in

ﬂ

uence from faster equilibration of fast-turnover compared to

slow-turnover carbon pools under changing inputs

13

). This study

provides an observational constraint on the latter uncertainty. As
the vast majority of the CMIP6 and CMIP5 models do not yet

represent vertically resolved deep soil carbon in permafrost or
peatlands, we focus our constraint on carbon change in the top 1
metre of soil. To ensure a fair like-for-like comparison we also
exclude the two CMIP6 models that do represent vertically-
resolved soil carbon (CESM2 and NorESM2), although this has a
negligible effect on our overall result. Our study therefore applies
to soil carbon loss in the top 1 metre of soil only. Below we show
that it is possible to signi

ﬁ

cantly reduce the uncertainty in this key

feedback to climate change using current-day spatial data to
constrain the sensitivity to future warming.

Results and discussion

Proof of concept

. For each ESM, we begin by calculating the

effective

τ

s

using time-averaged (1995

–

2005) values of

C

s

and

R

h

at each grid-point, and applying our de

ﬁ

nition of

τ

s

=

C

s

/

R

h

. We

do likewise for observational datasets of soil carbon in the top 1
metre

19

,

20

and

time-averaged

(2001

–

2010)

heterotrophic

respiration

, as shown in Fig.

. Figure

c shows the map of

inferred values of

τ

s

from these observations, with a notable

increase from approximately 7 years in the warm tropics to over
100 years in the cooler high northern latitudes.

Similar maps can be diagnosed for each of year of data, for each

ESM, and for each future scenario, giving time and space varying
values of

τ

s

for each model run. This allows us to estimate

Δ

C

s,

τ

,

via the last term on the right of Eq. (

). For each ESM, the

R

h,0

value is taken as the mean over the decade 1995

–

2005, to overlap

with the time period of the observations and to maintain
consistency across CMIP generations. Individual grid-point

τ

s

values are calculated for each year before calculating area-
weighted global totals of

Δ

C

s,

τ

. The uncertainty of

Δ

C

s,

τ

stems

from the uncertainty in soil carbon turnover (

τ

s

), and the

uncertainty due to differing climate sensitivities of the models. In
this study, we aim to quantify and constrain the uncertainty in

τ

s

.

To isolate the latter uncertainty, we consider

Δ

C

s,

τ

for differing

levels of global mean warming in each model. The resulting
dependence of global total

Δ

C

s,

τ

on global warming is shown in

Fig.

a, for each of the ESMs considered in both CMIP6 and

CMIP5 (seven CMIP6 ESMs and nine CMIP5 ESMs), and for
three Shared Socioeconomic Pathways (SSP): SSP126, SSP245 and
SSP585 (CMIP6)

22

, or the equivalent Representative Concentra-

tion Pathways (RCP): RCP2.6, RCP4.5 and RCP8.5 (CMIP5)

23

. In

all cases

Δ

C

s,

τ

is negative, which is consistent with the soil carbon

turnover time decreasing with warming. The more surprising
thing to note is the huge range in the projections, with a spread at
2 °C global mean warming of approximately 400 PgC, regardless
of future SSP/RCP scenario. Figure

b plots the fractional change

in soil carbon

Δ

C

s,

τ

/

C

s,0

, showing that there is a large range of

effective q

10

sensitivities between the model projections.

Unfortunately, we do not have time-varying observational

datasets of

C

s

and

R

h

that might allow us to directly constrain this

projection uncertainty. Instead we explore whether the observed
spatial variability in

τ

s

(as shown in Fig.

c) provides some

observational constraint on the sensitivity of

τ

s

to temperature. In

doing so, we are motivated by Chadburn et al.

15

who used the

correlation between the observed geographical distributions of
permafrost and air temperature to constrain projections of future
permafrost area under global warming. Similarly, we use ESMs to
test whether the spatial variation in

τ

s

reveals the sensitivity of soil

carbon turnover to temperature. The spatial patterns of

τ

s

in

CMIP5 simulations and observations were previously shown in
Koven et al.

, and here we test whether such relationships can be

used to estimate the response of soil carbon to future climate
change, using a combination of CMIP6 and CMIP5 models.

Figure

a is a scatter plot of log

τ

s

against temperature, using

the

τ

s

values shown in Fig.

c and mean temperatures from the

ARTICLE

NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-19208-8

2

NATURE COMMUNICATIONS

| (2020) 11:5544 | https://doi.org/10.1038/s41467-020-19208-8 | www.nature.com/naturecommunications

ACCESS-ESM1-5

Model colours

BCC-CSM2-MR

CanESM5

CNRM-ESM2-1

IPSL-CM6A-LR

MIROC-ES2L

UKESM1-0-LL

BNU-ESM

CanESM2

CESM1-CAM5

GFDL-ESM2G

GISS-E2-R

HadGEM2-ES

IPSL-CM5A-LR

MIROC-ESM

NorESM1-M

200

0

–200

–400

–600

–800

Model

Δ

C

s,

(PgC)

0

1

2

3

4

5

Δ

T

(°C)

0

–0.5

q

10

= 1.5

q

10

= 2

q

10

= 3

q

10

= 6

–0.4

–0.3

–0.2

–0.1

0.0

SSP126/RCP2.6

SSP245/RCP4.5

SSP585/RCP8.5

Model fr

actional change (

Δ

C

s,

/

C

s,

0

)

1

2

3

4

5

Δ

T

(°C)

a

b

Fig. 1 Uncertainty in future changes in soil carbon due to reduction in turnover time.

Δ

C

s,

τ

vs.

Δ

T

plot diagnosed from sixteen Earth System Models

(seven CMIP6 ESMs and nine CMIP5 ESMs), for three different future scenarios: SSP126, SSP245, SSP585, or RCP2.6, RCP4.5, RCP8.5, respectively.

a

The

change in soil carbon due to the change in soil carbon turnover time against change in global mean temperatures;

b

The fractional change in soil carbon due

to the change in soil carbon turnover time against change in global mean temperatures, and compared to different effective q

10

sensitivities.

80°N

60°N

40°N

20°N

0°

20°S

40°S

180°W

120°W

60°W

60°E

120°E

180°E

0°

80°N

60°N

40°N

20°N

0°

20°S

40°S

180°W

120°W

60°W

60°E

120°E

180°E

0°

80°N

60°N

40°N

20°N

0°

20°S

40°S

180°W

120°W

60°W

60°E

120°E

180°E

3

7

20

55

Inf

erred

s

(yr)

R

h

(kg C m

–2

yr

–1

)

C

S

(kg C m

–2

)

150

400

0.0

0.2

0.4

0.6

0.8

1.0

1.2

1.4

0

10

20

30

40

50

60

70

0°

a

b

c

Fig. 2 Spatial variability of soil carbon turnover time inferred from observations.

Maps of

a

observed soil carbon (

C

s

) to a depth of 1 m (kg C m

−

2

,

b

observed heterotrophic respiration (

R

h,0

) (kg C m

−

2

yr

−

1

)

, and

c

inferred soil carbon turnover time (log

τ

s

) (yr).

NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-19208-8

ARTICLE

NATURE COMMUNICATIONS

| (2020) 11:5544 | https://doi.org/10.1038/s41467-020-19208-8 | www.nature.com/naturecommunications

WFDEI dataset over the period 2001

–

2010

24

. The thick black-

dotted line is a quadratic

ﬁ

t through these points. Also shown for

comparison are equivalent quadratic

ﬁ

ts for each model (coloured

lines), using the model log

τ

s

and mean near-surface air

temperature (

T

) values for each grid-point, over an overlapping

period with the observations (1995

–

2005). There is a spread in

the individual data points due to variation in soil moisture, soil
type, and other soil parameters

25

. The model speci

ﬁ

c spread in

the data can be seen for the CMIP6 and CMIP5 models in
Supplementary Figs. 2 and 3, respectively. Although models do
not account for every possible factor contributing to this spread,
the spread of points in the models is generally similar to the
observations. However, differences between the best-

ﬁ

t functions

relating

τ

s

to

T

are evident between the models, and between the

models and the observations

.

This suggests that we may be able to constrain

Δ

C

s,

τ

using the

observed

τ

s

vs.

T

ﬁ

t from the observations, but only if we can

show that such functions can be used to predict

Δ

C

s,

τ

under

climate change. In order to test that premise, we attempt to
reconstruct the time-varying

Δ

C

s,

τ

projection for each model

using the time-invariant

τ

s

vs.

T

ﬁ

t across spatial points (Fig.

a),

and the time-invariant

R

h,0

ﬁ

eld. The change in soil carbon

turnover time (

Δ

τ

s

(t)) for a given model run is estimated at each

point based-on the

τ

s

vs.

T

curve, and the time-varying

projection of

T

at that point. A local estimate of the subsequent

change in soil carbon can then be made based-on the farthest
right-hand term of Eq. (

) (

R

h,0

Δ

τ

s

), which can be integrated up

to provide an estimated change in global soil carbon in the top 1
metre (

Δ

C

s,

τ

).

Figure

b shows the result of this test for all models and all

respective SSP/RCP scenarios. The axes of this plot show
equivalent variables which represent the global

Δ

C

s,

τ

between

the mean value for 2090

–

2100 and the mean value for 1995

–

2005.

The

y

-axis represents the actual values for each model as shown

in Fig.

, and the

x

-axis represents our estimate derived from

spatial variability (as in Fig.

a). As hoped, actual vs. estimated

values cluster tightly around a one-to-one line with an

r

2

correlation coef

ﬁ

cient value of 0.90. Although some hot-climate

regions will inevitably experience temperatures beyond those
covered by current-day spatial variability, these tend to be regions
with low soil carbon, so this does not have a major impact on the
success of our method.

Spatial emergent constraint

. This gives us con

ﬁ

dence to use the

τ

s

vs.

T

ﬁ

t and

R

h,0

from observations to constrain future pro-

jections of

Δ

C

s,

τ

. To remove the uncertainty in future

Δ

C

s,

τ

due to

the climate sensitivity of the models, we investigate a common
amount of global mean warming in each model. Figure

a is

similar to Fig.

b but instead for the more policy-relevant case of

2 °C of global warming. As before, the

y

-axis represents the

modelled

Δ

C

s,

τ

, and the

x

-axis is our estimate derived from spatial

variability. Once again, the actual and estimated values of

Δ

C

s,

τ

cluster around the one-to-one line (with

r

2

=

0.87). The model

range arises partly from differences in the initial

ﬁ

eld of hetero-

trophic respiration (

R

h,0

), and partly from differences in

Δ

τ

s

(compare

ﬁ

rst row to penultimate row of Table

).

The vertical green line in Fig.

a represents the mean estimate

when the

τ

s

vs.

T

relationship and the

R

h,0

ﬁ

eld from the model are

replaced with the equivalents from the observations. The spread
shown by the shaded area represents the relatively small impact on

Δ

C

s,

τ

of differences in modelled spatial climate change patterns at 2

°C of global warming. In order to estimate the remaining
uncertainty in

Δ

C

s,

τ

, we treat this spread as equivalent to an

observational uncertainty in an emergent constraint approach

26

.

We apply a standard statistical approach

27

,

28

to estimate the

probability density function of the

y

-axis variable (model

Δ

C

s,

τ

),

ACCESS-ESM1-5

Observations

Model colours

9

0

SSP126/RCP2.6
SSP245/RCP4.5
SSP585/RCP8.5

–100

–200

–300

–400

Model

Δ

C

s,

(PgC)

Relationship-derived

Δ

C

s,

(PgC)

–500

–600

–700

–600

–500

–400

–300

–200

–100

0

8

7

6

5

log(

s

) (log(yr))

4

3

2

1

T

(°

C

)

–20

–10

0

10

20

30

BCC-CSM2-MR

CanESM5

CNRM-ESM2-1

IPSL-CM6A-LR

MIROC-ES2L

UKESM1-0-LL

BNU-ESM

CanESM2

CESM1-CAM5

GFDL-ESM2G

GISS-E2-R

HadGEM2-ES

IPSL-CM5A-LR

MIROC-ESM

NorESM1-M

b

a

Fig. 3 Using spatial variability of soil carbon turnover time to estimate

Δ

C

s,

τ

. a

Scatter-plot of the relationship between log

τ

s

and mean air temperature

from observations

–

(black points), and a quadratic

ﬁ

t (black-dotted line) representing the observational temperature sensitivity of log

τ

s

. The

equivalent quadratic

ﬁ

ts for the ESMs are shown by the coloured lines;

b

The proof of principle for our method, showing an actual vs. estimated

comparison, representing the modelled versus the relationship-derived values of the

Δ

C

s,

τ

, where the change is considered between the start (1995

–

2005)

and the end (2090

–

2100) of 21st century and is assumed to relate to the top 1 metre of soil.

Table 1 Sensitivity study of spatial emergent constraint.

Constrained

Δ

C

s,

τ

at 2 °C global mean warming

Combined

CMIP6

CMIP5

CARDAMOM

R

h

−

232 ± 52

−

238 ± 62

−

227 ± 48

MODIS NPP

−

201 ± 53

−

206 ± 63

−

196 ± 49

Raich 2002

R

s

−

243 ± 50

−

249 ± 59

−

238 ± 46

CARDAMOM

R

h

(Observational

τ

s

v

T

ﬁ

t, model

R

h,0

)

−

227 ± 95

−

220 ± 75

−

230 ± 109

Unconstrained

Δ

C

s,

τ

−

196 ± 117

−

216 ± 109

−

180 ± 121

The table presents the sensitivity of the emergent constraint on

Δ

C

s,

τ

to model ensemble:

CMIP5, CMIP6 or CMIP5 and CMIP6 combined (columns), and to the observational dataset for
heterotrophic respiration (rows). The penultimate row presents the constraint using the
observational

τ

s

v T

ﬁ

t and model

R

h,0

, opposed to observational

R

h,0

, to isolate the uncertainty

reduction from these different components. For comparison, the last row shows the mean and
standard deviation of the unconstrained model ensemble.

ARTICLE

NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-19208-8

4

NATURE COMMUNICATIONS

| (2020) 11:5544 | https://doi.org/10.1038/s41467-020-19208-8 | www.nature.com/naturecommunications

accounting for both this observational spread and the quality of the
emergent relationship. To test the robustness to the choice of
observations we have repeated the analysis with different datasets
that represent heterotrophic respiration, which produces strongly-
overlapping emergent constraints, and completing the analysis with
both CMIP6 and CMIP5 models shows that the result is also robust
to the choice of model ensemble (see Table

).

Figure

b shows the resulting emergent constraint (blue line),

and compares to the unweighted histogram of model values (grey
blocks), and a Gaussian

ﬁ

t to that prior distribution (black line).

The spatial emergent constraint reduces the uncertainty in

Δ

C

s,

τ

at

2 °C of global warming from

−

196 ± 117 PgC to

−

232 ± 52 PgC

(where these are mean values plus and minus one standard
deviation for the top 1 metre). This same method can be applied to

ﬁ

nd constrained values of

Δ

C

s,

τ

for other values of global warming.

Figure

c shows the constrained range of

Δ

C

s,

τ

as a function of

global warming. This rules out the most extreme projections
but nonetheless suggests substantial soil carbon losses due to
climate change even in the absence of losses of deeper permafrost
carbon.

Methods

Obtaining spatial relationships

. In this section we explain how the quadratic

relationships representing the spatial log

τ

s

-temperature sensitivity shown in

Fig.

a (and Supplementary Figs. 2, 3 and 6) were derived, for both the Earth

System Models (ESMs) in CMIP6 and CMIP5, and using the observational data.
This is similar to the method used in Koven et al.

Obtaining spatial relationships for CMIP models

. The CMIP6 models used in

this study are shown in the Table

, and the CMIP5 models used in this study are

shown in Table

.

To obtain model speci

ﬁ

c spatial log

τ

s

-temperature relationships, the following

method was used. A reference time period was considered (1995

–

2005), this was

taken as the end of the CMIP5 historical simulation to be consistent across CMIP
generations and to best match the observational data time frame considered. Then,
monthly model output data was time averaged over this period, for the output
variables

‘

soil carbon content

’

(

C

s

) in kg m

−

2

,

‘

heterotrophic respiration carbon

ﬂ

ux

’

(

R

h

) in kg m

−

2

s

−

1

, and

‘

air temperature

’

in K. The variables

C

s

and

R

h

were

used to obtain values for soil carbon turnover time (

τ

s

) in years, using the equation

τ

s

=

C

s

/(

R

h

× 86400 × 365). The model temperature variable units were converted

from K to °C.

For each model, these values of log

τ

s

were plotted against the corresponding

spatial temperature data to obtain the spatial log

τ

s

-temperature plot. Then,

quadratic

ﬁ

ts (using the python package numpy poly

ﬁ

t) are calculated for each

model, which represent the spatial log

τ

s

relationship and sensitivity to

temperature. These model speci

ﬁ

c relationships are shown by the coloured lines in

Fig.

a in the main manuscript, and in Supplementary Fig. 2 for CMIP6 and in

Supplementary Fig. 3 for CMIP5.

Obtaining spatial relationships for observations

. Following Koven et al.

, we

estimated observational soil carbon data (to a depth of 1 m) by combining the
Harmonized World Soils Database (HWSD)

19

and Northern Circumpolar Soil

Carbon Database (NCSCD)

20

soil carbon datasets, where NCSCD was used where

overlap occurs. To calculate soil carbon turnover time,

τ

s

, using the following

equation:

τ

s

=

C

s

/

R

h

, we require a global observational dataset for heterotrophic

respiration. In the main manuscript, CARDAMOM (2001

–

2010) heterotrophic

respiration (

R

h

) is used

. We completed a sensitivity study on the choice of

observational heterotrophic respiration dataset, see below. The WFDEI dataset is

0

–100

–200

–300

Model

Δ

C

s,

(PgC)

Model

Δ

C

s,

(PgC)

Relationship-derived

Δ

C

s,

(PgC)

–400

–500

–600

0.000

0.001

0.002

0.003

Probab

lity density P

er PgC

0.004

0.005

0.006

0.007

0

–50

–100

–150

–200

–250

Magnitude of global

C

s

loss due to

Δ

s

(PgC)

–300

–350

–400

0.0

0.5

1.0

Global mean temperature warming (°C)

1.5

2.0

2.5

3.0

Emergent constraint
Gaussian model fit

q

10

= 2

q

10

= 2.5

q

10

= 3

Constrained Values
Confidence Limits

Model range

SSP126/RCP2.6
SSP245/RCP4.5
SSP585/RCP8.5
Observational Constraint
Emergent Constraint

–500 –400 –300 –200 –100

0

100

200

–500

–400

–300

–200

–100

0

c

b

a

Fig. 4 Emergent constraint on

Δ

C

s,

τ

as a function of global warming. a

Actual vs. estimated scatter plot for

Δ

C

s,

τ

for 2 °C of global warming. The vertical

green line de

ﬁ

nes the observational constraint which is derived using observational data and the future spatial temperature

ﬁ

eld of each model (decadal

average), and the shaded region represents the corresponding uncertainty (±1 standard deviation). The horizontal blue line represents our emergent
constraint, with the shaded region showing the corresponding uncertainty (±1 standard deviation) which results from the differing future spatial warming
patterns seen in the future spatial temperature

ﬁ

elds across the ESMs, and the emergent relationship between the model data points (black line).

b

Probability density function showing the Gaussian distribution of

Δ

C

s,

τ

values from the unweighted prior model ensemble (black line) and the emergent

constraint (blue line).

c

Constrained

Δ

C

s,

τ

values at different levels of global warming (blue), including the likely (±1 standard deviation) uncertainty bounds

(shaded blue). Different effective global q

10

values shown for comparison; our emergent constraint is consistent with an effective q

10

≈

2.5 ± 0.6.

Table 2 CMIP6 models.

Model

Institute

ACCESS-ESM1-5

Australian Community Climate and Earth
Systems Simulator, Australia

BCC-CSM2-MR

The Beijing Climate Center, China

CanESM5

Canadian Centre for Climate Modelling and
Analysis, Canada

CNRM-ESM2-1

CNRM/CERFACS, French Centre National
de la Recherche Scienti

ﬁ

que, France

IPSL-CM6A-LR

Institut Pierre-Simon Laplace, France

MIROC-ES2L

Atmosphere and Ocean Research Institute
and Japan Agency for Marine-Earth Science
and Technology, Japan

UKESM1-0-LL

NERC and Met Of

ﬁ

ce Hadley Centre, UK

Table 3 CMIP5 models.

Model

Institute

BNU-ESM

College of Global Change and Earth System
Science, China

CanESM2

Canadian Centre for Climate Modelling and
Analysis, Canada

CESM1-CAM5

National Science Foundation, Department of Energy,
NCAR, USA

GFDL-ESM2G

NOAA Geophysical Fluid Dynamics Laboratory, USA

GISS-ES-R

NASA Goddard Institute for Space Studies, USA

HadGEM2-ES

Met Of

ﬁ

ce Hadley Centre, UK

IPSL-CM5A-LR

Institut Pierre-Simon Laplace, France

MIROC-ESM

Atmosphere and Ocean Research Institute and
Japan Agency for Marine-Earth Science and
Technology, Japan

NorESM-M

Norwegian Climate Centre, Norway

NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-19208-8

ARTICLE

NATURE COMMUNICATIONS

| (2020) 11:5544 | https://doi.org/10.1038/s41467-020-19208-8 | www.nature.com/naturecommunications

5

used for our observational air temperatures (2001

–

2010)

24

. Then, these datasets

can be used to obtain the observational log

τ

s

-temperature relationship, using the

same quadratic

ﬁ

tting as with the models. This represents the

‘

real world

’

spatial

temperature sensitivity of log

τ

s

, and is shown by the thick-dotted-black line in

Fig.

a of the main manuscript. A comparison of the derived observational rela-

tionships can be seen in Supplementary Fig. 6.

Observational sensitivity study

. We completed a sensitivity study to investigate

our constraint dependence on the choice of observational heterotrophic respira-
tion dataset (CARDAMOM (2001

–

2010)

). The other observational datasets

considered are as follows: NDP-08

‘

Interannual Variability in Global Soil

Respiration on a 0.5 Degree Grid Cell Basis

’

dataset (1980

–

1994)

29

,

‘

Global spa-

tiotemporal distribution of soil respiration modelled using a global database

’

30

, and

MODIS net primary productivity (NPP) (2000

–

2014)

31

. Supplementary Fig. 4

shows scatter plots showing one-to-one comparisons of these observational data-
sets against one another, and Supplementary Fig. 5 shows the corresponding
comparisons of the equivalent log

τ

s

values calculated from each dataset.

The CARDAMOM

R

h

dataset is used in the main manuscript for the following

two main reasons:

ﬁ

rstly, we calculate

τ

s

using heterotrophic respiration which

allows for consistency between models and observations, and secondly, the dataset
does not use a prescribed q

10

sensitivity

. Instead, the CARDAMOM

R

h

dataset

was derived by explicitly assimilating observations into a process-based diagnostic
land-surface model. To test the robustness of our results, we also repeated our
analysis with MODIS NPP and Raich 2002, for both CMIP6 and CMIP5 together,
and as separate model ensembles. Supplementary Fig. 6 shows the observational
log

τ

s

-temperature relationships, derived using each of these observational datasets.

The results are presented in Table

which shows the constrained values of

Δ

C

s,

τ

at

2 °C global mean warming.

We decided not to complete the paper analysis using the Hashimoto dataset

since not only is it inconsistent with the three other datasets considered, it also
shows an arbitrary maximum respiration level (Supplementary Fig. 4), which likely
results from the assumed temperature-dependence of soil respiration in this dataset
which takes a quadratic form

30

. The quadratic form is justi

ﬁ

ed based on a site-level

study in which it is used to

ﬁ

t temporal dynamics. However, the parameters for the

quadratic function that are

ﬁ

tted in the Hashimoto study are very different from

those in the site-level study, which therefore suggests that the same relationship
does not apply to the global distribution of mean annual soil respiration.

Equation for the soil carbon turnover time component of soil carbon change

.

The equation used in this study for the component of the change in soil carbon (

Δ

C

s

)

due to the change in soil carbon turnover time (

Δ

τ

s

) was derived in the following way.

Starting with the equation for soil carbon (based on the de

ﬁ

nition of

τ

s

):

C

s

¼

R

h

τ

s

ð

2

Þ

As discussed in the main manuscript, we can write this change in soil carbon

(

Δ

C

s

), as the sum of a term arising from changes in litter-fall (

Δ

C

s,L

), and a term

arising from changes in the turnover time of soil carbon (

Δ

C

s,

τ

):

Δ

C

s

¼

Δ

ð

R

h

τ

s

Þ

Δ

C

s

;

L

ð

t

Þ þ

Δ

C

s

;

τ

ð

t

Þ ¼

τ

s

;

0

Δ

R

h

ð

t

Þ þ

R

h

;

0

Δ

τ

s

ð

t

Þ

ð

3

Þ

Hence, the equation for the component of soil carbon change due to the change

in

τ

s

is:

Δ

C

s

;

τ

¼

R

h

;

0

Δ

τ

s

ð

4

Þ

In this study we use

R

h

from the reference period (

‘

present day

’

), which we call

R

h,0

,

to allow us to investigate the response of

Δ

C

s,

τ

as a result of the response of

τ

s

to

climate change.

Modelled future temperature

. The proof of principle

ﬁ

gure (Fig.

b) con-

siders

Δ

C

s,

τ

between the end of the 21st century (2090

–

2100), for each future SSP

scenario (SSP126, SSP245, SSP585)

22

or equivalent future RCP scenario (RCP2.6,

RCP4.5 and RCP8.5)

23

, and our reference period from the historical simulation

(1995

–

2005), for each CMIP6 ESM and CMIP5 ESM, respectively.

To consider speci

ﬁ

c °C of global warming (Fig.

), the future spatial

temperature pro

ﬁ

les at these speci

ﬁ

c global mean warming levels, for example:

1 °C, 2 °C and 3 °C global mean warming, were calculated as follows. The
temperature change is calculated from our reference period (1995

–

2005), and then

a 5-year rolling mean of global mean temperature is taken to remove some of the
interannual variability. Once the year that the given temperature increase has been
reached is obtained, a time average including

−

5 and

+

5 years is taken, and the

spatial temperature distribution of that model averaged over the deduced time
period is used for the calculations of future

τ

s

.

Anomaly correction for future temperature projections

. To remove uncertainty

due to errors in the models

’

historical simulation, a spatial future temperature

anomaly was projected using each model and each respective future SSP/RCP
scenario separately. To calculate this, the temperature at the reference time frame
(1995

–

2005), which overlaps the WFDEI observational temperature data time

frame (2001

–

2010), is subtracted from the future temperature pro

ﬁ

le for each

model (as calculated above), to calculate the temperature change. Then, this
temperature anomaly is added onto the observational temperature dataset to give a
model-derived future

‘

observational

’

temperature for each model.

Proof of concept for our method

. Our method relies on the idea that the spatial

temperature sensitivity can be used to project and constrain the temporal sensi-
tivity of

τ

s

to temperature, and subsequently global warming. To test the robustness

of this method,

Δ

C

s,

τ

calculated using model

Δ

τ

s

, and temperature sensitivity

relationship-derived

Δ

τ

s

, are compared.

The change in soil carbon turnover time (

Δ

τ

s

) was either calculated using model

output data to obtain model-derived

Δ

τ

s

as follows:

Δ

τ

s

¼

τ

f
s

τ

h
s

ð

5

Þ

where,

τ

s

¼

C

s

=

R

h

ð

6

Þ

Or calculated using the derived quadratic log

τ

s

-temperature relationships to

obtain relationship-derived

Δ

τ

s

, which is based on the following equation:

Δ

τ

s

¼

exp

ð

p

ð

T

f

ÞÞ

exp

ð

p

ð

T

h

ÞÞ

ð

7

Þ

where,

T

is near surface air temperature, and

T

f

represents a future temperature,

and

T

h

represents historical (present day) temperature from our reference period

(1995

–

2005). The exponentials (exp) are taken to turn log

τ

s

values to

τ

s

values.

p(

T

) represents the quadratic log

τ

s

-temperature relationship as a function of

temperature to obtain our estimated log

τ

s

.

These

Δ

τ

s

values are then put back into the Eq. (

) (with model-speci

ﬁ

c

R

h,0

) to

obtain the corresponding

Δ

C

s,

τ

values. The proof of principle

ﬁ

gure (Fig.

b)

investigates the robustness of our method, where projections of model and
relationship-derived values of

Δ

C

s,

τ

are compared, and an

r

2

value of 0.90 is

obtained. The correlation of the data was also tested when investigating different
levels of global mean warming to obtain the constrained values (Fig.

). The

r

2

values for were as follows: 1 °C is 0.84, 2 °C is 0.87 and 3 °C is 0.87.

Calculating constrained values

. To obtain the constrained values of

Δ

C

s,

τ

, the

model-derived future

‘

observational

’

temperature for each model is used together

with the observational derived log

τ

s

-temperature relationship, to project values for

future

τ

s

. Then this together with relationship-derived historical

τ

s

deduced using

the observational temperature dataset, can be used to calculate

Δ

τ

s

. Finally global

Δ

C

s,

τ

can be obtained by multiplying

Δ

τ

s

by the observational dataset for

R

h,0

(using Eq. (

)), and then calculating a weighted-global total. As each model-derived

future

‘

observational

’

temperature is considered separately, we obtain a range of

projected observational-constrained

Δ

C

s,

τ

values.

We have now obtained a set of

x

and

y

values, corresponding to the

relationship-derived and modelled values of

Δ

C

s,

τ

, respectively, for each ESM.

Where we have an

x

and

y

value for each model, representing the modelled

Δ

C

s,

τ

(

y

values), and the model speci

ﬁ

c relationship-derived

Δ

C

s,

τ

(

x

values). We also

have an

x

obs

value representing the mean observational-constrained

Δ

C

s,

τ

value,

and a corresponding standard deviation due to the uncertainty in the modelled
spatial pro

ﬁ

les of future temperatures. We follow the method used in Cox et al.

2018, which can be seen in the

‘

Least-squares linear regression

’

section and the

‘

Calculation of the PDF for ECS

’

section of the methods from this study

32

. Using

this method, we obtain an emergent relationship between our

x

and

y

data points,

which we can use together with our x

obs

and corresponding standard deviation to

produce a constraint on our

y

-axis. This is shown in Fig.

a. From this we obtain a

constrained probability density function on

Δ

C

s,

τ

, with a corresponding

uncertainty bounds which we consider at the 68% con

ﬁ

dence limits (±1 standard

deviation). Figure

b show the probability density functions representing the

distribution of the range of projections, before and after the constraint.

This method allows us to calculate a constrained probability density function on

Δ

C

s,

τ

at each °C of global mean warming, using the data seen in Fig.

a for 2 °C

warming, and our corresponding constrained values for 1 °C and 3 °C warming.
Figure

c shows the resultant constrained mean value of

Δ

C

s,

τ

obtained for each °C

of global mean warming, and the corresponding uncertainty bounds at the 68%
con

ﬁ

dence limits (±1 standard deviation).

Calculating effective q

10

for change in soil carbon

. Simple models of soil carbon

turnover are often based on just a q

10

function, which means that

τ

s

depends on

temperature as follows:

τ

s

¼

τ

s

;

0

exp

ðð

0

:

1 log q

10

Þ

Δ

T

Þ

ð

8

Þ

We compared the results for

Δ

C

s,

τ

that would be derived from a simple q

10

function with our emergent constraint results for

Δ

C

s,

τ

, to estimate an effective q

10

sensitivity of heterotrophic respiration.

To do this, we can obtain an equation for

Δ

τ

s

derived from Eq. (

8

). This is done by

considering the following, where

τ

s,0

is an initial

τ

s

, we can substitute in

τ

s

in temperature

sensitivity form to obtain an equation for

Δ

τ

s

in temperature sensitivity form:

Δ

τ

s

¼

τ

s

τ

s

;

0

ð

9

Þ

ARTICLE

NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-19208-8

6

NATURE COMMUNICATIONS

| (2020) 11:5544 | https://doi.org/10.1038/s41467-020-19208-8 | www.nature.com/naturecommunications

Δ

τ

s

¼

τ

s

;

0

exp

ðð

0

:

1 log q

10

Þ

Δ

T

Þ

τ

s

;

0

ð

10

Þ

Then, we can substitute this

Δ

τ

s

into Eq. (

) and simplify to obtain an equation

relating

Δ

C

s,

τ

and

Δ

T

:

Δ

C

s

;

τ

¼

R

h

;

0

τ

s

;

0

½

exp

ðð

0

:

1 log q

10

Þ

Δ

T

Þ

1

ð

11

Þ

Δ

C

s

;

τ

¼

C

s

;

0

½

exp

ðð

0

:

1 log q

10

Þ

Δ

T

Þ

1

ð

12

Þ

This equation was used to calculate different

Δ

C

s,

τ

–

Δ

T

sensitivity curves based on

different values on q

10

, for example q

10

=

2, with different amounts of global mean

warming to represent

Δ

T

, and initial observational soil carbon stocks

C

s,0

. These curves

can be seen on Figs.

b and

c. Note that there is no direct relationship between the

effective q

10

for soil carbon change shown in Figs.

b and

c, and the spatial

τ

s

–

T

relationships in Fig.

https://esgf-node.llnl.gov/search/CMIP5/

a. Our q

10

value is an effective q

10

value that indicates the

sensitivity of global soil carbon (in the top 1 metre) to global mean temperature.

Data availability

The datasets analysed during this study are available online: CMIP5 model output
[

], CMIP6 model output [

https://esgf-node.llnl.

gov/search/cmip6/

], The WFDEI Meteorological Forcing Data [

https://rda.ucar.edu/

datasets/ds314.2/

], CARDAMOM Heterotrophic Respiration [

https://datashare.is.ed.ac.

uk/handle/10283/875

], MODIS Net Primary Production [

https://lpdaac.usgs.gov/

products/mod17a3v055/

], Raich et al. 2002 Soil Respiration [

https://cdiac.ess-dive.lbl.

gov/epubs/ndp/ndp081/ndp081.html

], Hashimoto et al. 2015 Heterotrophic Respiration

[

http://cse.ffpri.affrc.go.jp/shojih/data/index.html

], and the datasets for observational Soil

Carbon [

https://github.com/rebeccamayvarney/soiltau_ec

].

Code availability

The Python code used to complete the analysis and produce the

ﬁ

gures in this study is

available in the following online repository [

https://github.com/rebeccamayvarney/

soiltau_ec

].

Received: 13 January 2020; Accepted: 17 September 2020;

References

1.

Ciais, P. et al. Carbon and other biogeochemical cycles. In

Climate Change

2013: The Physical Science Basis. Contribution of Working Group I to the Fifth
Assessment Report of the Intergovernmental Panel on Climate Change

, Book

section 6, (eds. Stocker, T. et al.) 465

–

570 (Cambridge University Press,

Cambridge, New York, 2013).

2.

UNFCCC.

Adoption of the Paris Agreement

. fccc/cp/2015/l.9/rev.1.

http://

unfccc.int/resource/docs/2015/cop21/eng/l09r01.pdf

(2015).

3.

Jenkinson, D. S., Adams, D. & Wild, A. Model estimates of CO

2

emissions

from soil in response to global warming.

Nature

351

, 304 (1991).

4.

Cox, P. M., Betts, R. A., Jones, C. D., Spall, S. A. & Totterdell, I. J. Acceleration
of global warming due to carbon-cycle feedbacks in a coupled climate model.

Nature

408

, 184 (2000).

5.

Bradford, M. A. et al. Managing uncertainty in soil carbon feedbacks to
climate change.

Nat. Clim. Change

6

, 751 (2016).

6.

Todd-Brown, K. et al. Changes in soil organic carbon storage predicted by earth
system models during the 21st century.

Biogeosciences

11

, 2341

–

2356 (2014).

7.

Arora, V. K. et al. Carbon

–

concentration and carbon

–

climate feedbacks in

CMIP5 earth system models.

J. Clim.

26

, 5289

–

5314 (2013).

8.

Jones, C. et al. Twenty-

ﬁ

rst-century compatible CO

2

emissions and airborne

fraction simulated by CMIP5 earth system models under four representative
concentration pathways.

J. Clim.

26

, 4398

–

4413 (2013).

9.

Knorr, W., Prentice, I. C., House, J. & Holland, E. Long-term sensitivity of soil
carbon turnover to warming.

Nature

433

, 298 (2005).

10. Todd-Brown, K. et al. Causes of variation in soil carbon simulations from

CMIP5 earth system models and comparison with observations.

Biogeosciences

10

, 1717

–

1736 (2013).

11. Eyring, V. et al. Overview of the coupled model intercomparison project phase

6 (CMIP6) experimental design and organization.

Geosci. Model Dev.

9

,

1937

–

1958 (2016).

12. Taylor, K. E., Stouffer, R. J. & Meehl, G. A. An overview of CMIP5 and the

experiment design.

Bull. Am. Meteorol. Soc.

93

, 485

–

498 (2012).

13. Koven, C. D. et al. Controls on terrestrial carbon feedbacks by productivity

versus turnover in the CMIP5 earth system models.

Biogeosciences

12

,

5211

–

5228 (2015).

14. Hall, A., Cox, P., Huntingford, C. & Klein, S. Progressing emergent constraints

on future climate change.

Nat. Clim. Change

9

, 269

–

278 (2019).

15. Chadburn, S. et al. An observation-based constraint on permafrost loss as a

function of global warming.

Nat. Clim. Change

7

, 340 (2017).

16. Koven, C. D., Hugelius, G., Lawrence, D. M. & Wieder, W. R. Higher

climatological temperature sensitivity of soil carbon in cold than warm
climates.

Nat. Clim. Change

7

, 817 (2017).

17. Carvalhais, N. et al. Global covariation of carbon turnover times with climate

in terrestrial ecosystems.

Nature

514

, 213 (2014).

18. Davidson, E. A. & Janssens, I. A. Temperature sensitivity of soil carbon

decomposition and feedbacks to climate change.

Nature

440

, 165 (2006).

19. FAO, I. & ISRIC, I. Jrc: Harmonized world soil database (version 1.2). FAO,

Rome, Italy and IIASA, (Laxenburg, Austria, 2012).

20. Hugelius, G. et al. The northern circumpolar soil carbon database: spatially

distributed datasets of soil coverage and soil carbon storage in the northern
permafrost regions.

Earth System Science Data

5

, 3

–

13 (2013).

21. Bloom, A. et al. CARDAMOM 2001-2010 global carbon model-data fusion

(mdf) analysis. (2015).

22. O

’

Neill, B. C. et al. The scenario model intercomparison project (scenariomip)

for CMIP6 (2016).

23. Meinshausen, M. et al. The rcp greenhouse gas concentrations and their

extensions from 1765 to 2300.

Clim. Change

109

, 213 (2011).

24. Weedon, G. P. et al. The wfdei meteorological forcing data set: watch forcing

data methodology applied to era-interim reanalysis data.

Water Res. Res.

50

,

7505

–

7514 (2014).

25. Exbrayat, J.-F., Pitman, A. J., Zhang, Q., Abramowitz, G. & Wang, Y.-P.

Examining soil carbon uncertainty in a global model: response of microbial
decomposition to temperature, moisture and nutrient limitation.

Biogeosciences

10

, 7095

–

7108 (2013).

26. Eyring, V. et al. Taking climate model evaluation to the next level.

Nat. Clim.

Change

9

, 102

–

110 (2019).

27. Cox, P. M. et al. Sensitivity of tropical carbon to climate change constrained by

carbon dioxide variability.

Nature

494

, 341 (2013).

28. Wenzel, S., Cox, P. M., Eyring, V. & Friedlingstein, P. Projected land

photosynthesis constrained by changes in the seasonal cycle of atmospheric co
2.

Nature

538

, 499 (2016).

29. Raich, J. W., Potter, C. S. & Bhagawati, D. Interannual variability in global soil

respiration, 1980

–

94.

Global Change Biol.

8

, 800

–

812 (2002).

30. Hashimoto, S. et al. Global spatiotemporal distribution of soil respiration

modeled using a global database.

Biogeosci.

12

, 4121

–

4132 (2015).

31. Zhao, M., Heinsch, F. A., Nemani, R. R. & Running, S. W. Improvements of

the modis terrestrial gross and net primary production global data set.

Remote

Sens. Environ.

95

, 164

–

176 (2005).

32. Cox, P. M., Huntingford, C. & Williamson, M. S. Emergent constraint on

equilibrium climate sensitivity from global temperature variability.

Nature

553

, 319 (2018).

Acknowledgements

This work was supported by the European Research Council

‘

Emergent Constraints on

Climate-Land feedbacks in the Earth System (ECCLES)

’

project, grant agreement

number 742472 (R.M.V. and P.M.C.). S.E.C. was supported by a Natural Environment
Research Council independent research fellowship, grant no. NE/R015791/1. P.F. was
supported by the European Union

’

s Framework Programme Horizon 2020 for Research

and Innovation under Grant Agreement No. 821003, Climate-Carbon Interactions in the
Current Century (4C) project. E.J.B. was funded by the European Commission

’

s Horizon

2020 Framework Programme, under Grant Agreement number 641816, the

‘

Coordinated

Research in Earth Systems and Climate: Experiments, Knowledge, Dissemination and
Outreach (CRESCENDO)

’

project (11/2015

–

10/2020) and the Met Of

ﬁ

ce Hadley Centre

Climate Programme funded by BEIS and Defra. C.D.K. acknowledges support from the
U.S. Department of Energy, Of

ﬁ

ce of Science, Of

ﬁ

ce of Biological and Environmental

Research, through the RUBISCO SFA and the Early Career Research Program. G.H. was
supported by the Swedish Research Council (2018-04516) and the EU H2020 project
Nunataryuk (773421). We also acknowledge the World Climate Research Programme

’

s

Working Group on Coupled Modelling, which is responsible for CMIP, and we thank the
climate modelling groups (listed in Tables

and