Back to report · Publisher record

Sanderman et al. (2017). Soil carbon debt of 12,000 years of human land use

Page 1

Figures and page graphics, page 1

SUSTAINABILITY

SCIENCE

Soil carbon debt of 12,000 years of human land use

Jonathan Sanderman

a,1,2

, Tomislav Hengl

b,1

, and Gregory J. Fiske

a

a

Woods Hole Research Center, Falmouth MA 02540; and

b

The International Soil Resource and Information Center – World Soil Information, 6708 PB

Wageningen, The Netherlands

Edited by William H. Schlesinger, Cary Institute of Ecosystem Studies, Millbrook, NY, and approved July 14, 2017 (received for review April 12, 2017)

Human appropriation of land for agriculture has greatly altered
the terrestrial carbon balance, creating a large but uncertain car-
bon debt in soils. Estimating the size and spatial distribution of
soil organic carbon (SOC) loss due to land use and land cover
change has been difficult but is a critical step in understand-
ing whether SOC sequestration can be an effective climate mit-
igation strategy. In this study, a machine learning-based model
was fitted using a global compilation of SOC data and the His-
tory Database of the Global Environment (HYDE) land use data
in combination with climatic, landform and lithology covariates.
Model results compared favorably with a global compilation of
paired plot studies. Projection of this model onto a world with-
out agriculture indicated a global carbon debt due to agricul-
ture of 133 Pg C for the top 2 m of soil, with the rate of loss
increasing dramatically in the past 200 years. The HYDE classes
“grazing” and “cropland” contributed nearly equally to the loss
of SOC. There were higher percent SOC losses on cropland but
since more than twice as much land is grazed, slightly higher
total losses were found from grazing land. Important spatial pat-
terns of SOC loss were found: Hotspots of SOC loss coincided
with some major cropping regions as well as semiarid grazing
regions, while other major agricultural zones showed small losses
and even net gains in SOC. This analysis has demonstrated that
there are identifiable regions which can be targeted for SOC
restoration efforts.

agriculture

|

soil organic matter

|

climate change

|

soil degradation

T

he incredible rise of human civilizations and the contin-
uing sustainability of current and future human societies

are inextricably linked to soils and the wide array of services
soils provide (1–3). Human population and economic growth
has led to an exponential rise in use of soil resources. Roughly
50 million km

2

of soils are currently being managed to some

degree by humans for food, fiber, and livestock production (4),
leading to the declaration that we live on a “used planet” (5).
The consequences of human domination of soil resources are
far ranging (6, 7): accelerated erosion, desertification, saliniza-
tion, acidification, compaction, biodiversity loss, nutrient deple-
tion, and loss of soil organic matter (SOM).

Of these soil threats, loss of SOM has received the most atten-

tion, due to the critical role SOM plays in the contemporary car-
bon cycle (8, 9) and as a key component of sustaining food pro-
duction (10, 11). Despite the intense research interest in SOM
and soil organic carbon (SOC) as the dominant component of
SOM, there remain many unknowns (12) that impede progress
in implementing sound land management strategies to rebuild
SOC stocks (13).

Conversion of native soil to agricultural uses typically leads to

a decline in SOC levels (14–16). The rate and extent of decline
in SOC stocks should vary greatly across the globe, due to dif-
ferences in soil properties, climate, type of land-use conversion,
and, importantly, the specific management implementation of a
given form of land use. Loss of SOC under agricultural land use
is not universal; modest gains are seen when soil of naturally
low fertility is improved and the previous constraint (e.g., mois-
ture, fertility, and hardpan) on plant growth is alleviated (17–
19). However, for the vast majority of land, SOC loss is more

common. In fact, in a metaanalysis of the available literature (

SI

Appendix

), we found median SOC loss values of 26% for the

upper 30 cm and 16% for the top 100 cm of soil, but ranges of

−

36 to 78% and

−

25 to 61%, respectively, have been reported

for these two depth increments (

SI Appendix

, Fig. S2

). Scaling

these limited point measurements to calculate a cumulative SOC
loss for the world’s agricultural land has been difficult, with esti-
mates ranging from 40 Pg C to over 500 Pg C (20). Recent esti-
mates from dynamic global vegetation models run with actual
land use versus with potential natural vegetation have put this
figure at 30 Pg C to 62 Pg C for the industrial post-1850 period
(21, 22).

A credible estimate of the global total and spatial distribu-

tion of SOC loss is a critical step in understanding the potential
for soil carbon sequestration to be an effective climate abate-
ment strategy. To quantify the cumulative impact of human land
use on changes in SOC at the global scale, we have developed
a machine learning-based data-driven statistical model, which
is based on a global compilation (

N

= 158,147) of soil pro-

file observations and samples collated at The International Soil
Resource and Information Center (ISRIC) – World Soil Infor-
mation over the last decade (23). We predicted current SOC
stocks as a function of climatic, topographic, geologic, and land-
use characteristics. By using the History Database of the Global
Environment (HYDE) v3.2 spatially explicit database of historic
pattern of human land use (4), we projected the predictive SOC
model back in time, holding climate constant, to estimate his-
toric SOC levels, with the difference between historic and current
SOC levels representing the amount of SOC lost due to human
land use.

Significance

Land use and land cover change has resulted in substantial
losses of carbon from soils globally, but credible estimates of
how much soil carbon has been lost have been difficult to
generate. Using a data-driven statistical model and the His-
tory Database of the Global Environment v3.2 historic land-
use dataset, we estimated that agricultural land uses have
resulted in the loss of 133 Pg C from the soil. Importantly,
our maps indicate hotspots of soil carbon loss, often associ-
ated with major cropping regions and degraded grazing lands,
suggesting that there are identifiable regions that should be
targets for soil carbon restoration efforts.

Author contributions: J.S. and T.H. designed research; J.S., T.H., and G.J.F. performed

research; J.S., T.H., and G.J.F. contributed new reagents/analytic tools; J.S., T.H., and G.J.F.
analyzed data; and J.S. and T.H. wrote the paper.

The authors declare no conflict of interest.

This article is a PNAS Direct Submission.

Freely available online through the PNAS open access option.

Data deposition: Meta-analysis data have been deposited on Harvard University’s website,

https://dataverse.harvard.edu/dataset.xhtml?persistentId=doi:10.7910/DVN/QQQM8V

.

Code and model outputs have been archived on GitHub,

https://github.com/whrc/Soil-

Carbon-Debt

.

1

J.S. and T.H. contributed equally to this work.

2

To whom correspondence should be addressed. Email: jsanderman@whrc.org.

This article contains supporting information online at

www.pnas.org/lookup/suppl/doi:10.

1073/pnas.1706103114/-/DCSupplemental

.

www.pnas.org/cgi/doi/10.1073/pnas.1706103114

PNAS

|

September 5, 2017

|

vol. 114

|

no. 36

|

9575–9580

Downloaded from https://www.pnas.org by 98.155.36.23 on October 19, 2026 from IP address 98.155.36.23.

SEE CORRECTION FOR THIS ARTICLE

Page 2

Figures and page graphics, page 2

Results and Discussion

Model Performance and Predictors.

Models explaining SOC den-

sity distribution as a function of depth, climate, relief, lithol-
ogy, land cover, and land use were evaluated using fivefold
cross-validation with repeated fitting. For the 0- to 2-m depth
interval,

R

-square was 54% with a root mean square error

(RMSE) of 12 kg C

·

m

−

3

based on random forest model, and

45% with a RMSE of 14 kg C

·

m

−

3

based on the gradient boost-

ing model. Correlation plots indicated that the ensemble model
for the 0- to 2-m depth interval was unbiased (no overesti-
mation/underestimation) over the whole span of values. Model
errors were lowest where sample density was high, and were
greatest in regions with sparse data and with high SOC values
(

SI Appendix

, Fig. S4

).

While climate and topographic attributes were the most

important variables in both models (

SI Appendix

, Fig. S3

), mul-

tiple land-use variables have been shown to be also important
for explaining the distribution of current SOC. Individual cor-
relation plots for the land-use variables generally decline with
increasing intensity of land use within a pixel (

SI Appendix

,

Fig. S3

). The relatively high importance of the HYDE land use

classes “grazing” and “cropland” in our SOC models validates
our conceptual approach to estimating historic SOC stocks in the
absence of human land use.

Current SOC Stocks.

Projecting our data-driven statistical model

across the globe for the year 2010 suggested that global SOC
stocks were 863 Pg C, 1,824 Pg C, and 3,012 Pg C in the upper
0.3 m, 1 m, and 2 m of soil, respectively (

SI Appendix

, Fig. S5

A

and Table S3

). On a per hectare basis, land classified as crop-

land [International Geosphere-Biosphere Programme (IGBP)
class 12 (24)] contained an average of 62 [31 to 96, 95%
confidence intervals (CI)], 127 (61 to 200, 95% CI), and 198
(96 to 315, 95% CI) Mg C

·

ha

−

1

for the same depth incre-

ments. The global total SOC stocks and depth distribution are
consistent with most recent estimates of this quantity using a
number of approaches (25, 26). This agreement is particularly
encouraging because we excluded remote sensing data prod-
ucts characterizing current vegetation, which has been found to
be particularly important in explaining SOC distribution (27),
because of the requirement to be able to project this model back
in time.

Historic SOC Stocks.

Reprojection of our SOC model to a no

land-use (NoLU) condition with all other variables held con-
stant resulted in global SOC stocks of 899 Pg C, 1,899 Pg C, and
3144 Pg C in the upper 0.3 m, 1 m, and 2 m of soil, respectively
(

SI Appendix

, Fig. S5

B

and Table S3

), suggesting that human-

driven land-use decisions have resulted in substantial reduc-
tions in global SOC levels. In the absence of accurate SOC data
from past millenia, we have attempted to assess the accuracy
of these historic projections by comparing the modeled NoLU
SOC stocks to SOC measurements taken in remnant patches
of native vegetation which were compiled from the literature
(

SI Appendix

, Fig. S6

). Given the limitations of this comparison

between point measurements and 10-km model output, model
results for NoLU SOC compared favorably to the measured liter-
ature values (

R

-square = 0.33 to 0.34, with RMSE values of 17 Mg

C

·

ha

−

1

and 37 Mg C

·

ha

−

1

for 0.3- and 1-m depths, respectively).

SOC Loss Due to Land Use.

Subtracting current (2010) SOC

stocks from historic (NoLU) SOC stocks, we found that 37 Pg C,
75 Pg C, and 133 Pg C have been lost due to land-use change in
the upper 0.3 m, 1 m, and 2 m of soil (Fig. 1

B

and

SI Appendix

,

Table S3

). The mean absolute loss due to land use to 2 m across

all pixels with some degree of land use was 17.7 Mg C

·

ha

−

1

(

−

22.7 to 54.0, 95% CI). As a percentage of initial SOC stocks

(

SI Appendix

, Fig. S7

), this represents a mean loss of 8.1% (

−

3.9

to 21.8%, 95% CI). While SOC was lost throughout the soil
profile, consistent with findings from experimental studies (14),
there was an exponential decline in loss with increasing depth (

SI

Appendix

, Fig. S8

).

Comparison of model results with the native remnant database

suggests that modeled SOC loss due to land use was likely a
conservative estimate. First, comparing measured SOC in native
patches to predicted SOC stocks with NoLU, we found a nega-
tive bias of

−

10 Mg C

·

ha

−

1

(

SI Appendix

, Fig. S6

). Second, it was

found that there was 13.8 to 28.0% (95% CI,

n

= 62

) less SOC

to 1 m under agricultural land compared with paired remnant
native vegetation plots (see

SI Appendix

, Fig. S2

and extended

discussion); however, modeled loss (NoLU – 2010) for the same
sites was only 7.7 to 12.8% (95% CI,

n

= 62

) of the NoLU

SOC to 1 m.

This analysis clearly demonstrates that, while, on average, agri-

cultural land use leads to SOC loss, there are important spatial
patterns and contrasts (Fig. 1) suggesting that simple account-
ing methods for SOC change [e.g., Intergovernmental Panel on
Climate Change default emission factors (28)] will misrepre-
sent SOC change in many regions of the world. The major-
ity of data used to generate published emission factors come
from North America and Europe (29). In our own metaanaly-
sis (

SI Appendix

), 82 of the 140 paired comparisons were from

these two regions. Model results from the agricultural heartland
of the United States and much of Europe showed large losses
(Fig. 1

B

), primarily from cropping, that are consistent with the

average documented changes in SOC stocks from field investi-
gations (14, 15).

The largest per pixel losses were found to coincide with crop-

ping regions (Fig. 1

B

); however, grazing, especially in arid and

semiarid regions, with its larger spatial extents (Fig. 1

A

), was

responsible for at least half of the total SOC loss (

SI Appendix

,

Table S3

). The grassland and savanna IGBP land classification

categories collectively lost more SOC than the cropland and
crop/natural vegetation mosaic categories (48 Pg C vs. 35 Pg C).
In particular, the rangelands of Argentina, southern Africa, and
Australia stand out as hotspots of SOC loss when viewed as a
percent of historic SOC (

SI Appendix

, Fig. S7

).

While land use is the underlying anthropogenic driver of SOC

loss, the degree to which land use results in SOC loss is at least
partially dependent upon the degree to which the soil resource
has been exploited (2). Using the Global Land Degradation
Information System biophysical status of land index, a quanti-
tative expression of a given land area’s ability to provide four
categories of ecosystem services (biomass, soil, water, and biodi-
versity) (30), we found a strong correlation between land degra-
dation and SOC loss (

SI Appendix

, Fig. S9

). Results from

SI

Appendix

, Fig. S9

, coupled with the finding that grazing was

the single most important land-use variable in the model (

SI

Appendix

, Fig. S3

), suggest that grazing of relatively unmanaged

rangelands may be a stronger driver of SOC loss than previously
acknowledged.

Agricultural land uses do not always result in large losses

of SOC. For example, the 2-million-km

2

seasonally dry moist

savanna region in Brazil known as “the Cerrado” was thought
to have soils too poor to support intensive agriculture but, over
the past few decades, has been transformed into one of the most
productive agricultural regions of the world, with

∼

750,000 km

2

of crops and 800,000 km

2

of pasture, through liming, fertiliza-

tion, and weed control (31). Because of the naturally infertile
soil state in the Cerrado, agricultural expansion has resulted in
little loss in SOC (Fig. 1

B

), and there are large areas that have

actually accumulated modest amounts of SOC after the advent
of agriculture, which is consistent with results from field studies
(17, 32, 33).

9576

|

www.pnas.org/cgi/doi/10.1073/pnas.1706103114

Sanderman et al.

Downloaded from https://www.pnas.org by 98.155.36.23 on October 19, 2026 from IP address 98.155.36.23.

Page 3

Figures and page graphics, page 3

SUSTAINABILITY

SCIENCE

Fig. 1.

Global distribution of cropping and grazing in 2010 from (

A

) HYDE v3.2 and (

B

) modeled SOC change in the top 2 m. In

A

, color gradients indicate

proportion of grid cell occupied by given land use. In

B

, legend is presented as histogram of SOC loss (Mg C

·

ha

−

1

), with positive values indicating loss and

negative values depicting net gains in SOC.

Historic SOC Loss Trend.

By using the HYDE v3.2 land-use

dataset, we predicted the temporal evolution of SOC stocks due
to changes in land use alone (Fig. 2). Globally, SOC loss follows
the exponential rise in used land, but not in a linear fashion (Fig. 2,

Inset

). There were low annual rates of SOC loss (

<

0.05 Pg C

·

y

−

1

)

until AD 1800, followed by a century of losses

>

0.3 Pg C

·

y

−

1

,

with a slight moderation of this rate during the last century
(0.13 Pg C

·

y

−

1

). Comparison of this historic trend among the

10 countries with the greatest losses reveals some interesting
contrasts (

SI Appendix

, Fig. S10

), with old-world countries (e.g.,

India, Kazakhstan, and Germany) showing large losses before AD
1000 but new-world countries (e.g., Brazil and Argentina) only
losing significant amounts of SOC in the past century.

Fig. 2.

Historic reconstruction of loss in SOC relative to 10,000 BC (assumed NoLU). Temporal evolution of cropland and grazing land is given in stacked

area plots. (

Inset

) Biplot of SOC loss (Pg C) v. total used land area (10

6

km

2

) for each predicted time interval.

Limitations of This Study.

While great care has been taken to

ensure that the input data were of the highest quality possi-
ble (23, 27), there remain several limitations in the underlying
datasets and therefore predicted SOC change. First, the train-
ing dataset used to build spatial predictions models was not ideal
for testing the hypotheses. ISRIC’s soil profile dataset is a com-
pilation of national inventories from a large number of nations
with different reasons for undertaking soil surveys and different
methods of laboratory analysis. Data were collected over a 50-y
period, which is likely smoothing out some of the SOC loss in
the model. In addition, the mismatch in scale between a soil
pedon (0.5 m on a side) and the pixel size of the HYDE v3.2
land-use data (10 km) can create situations where the dominant

Sanderman et al.

PNAS

|

September 5, 2017

|

vol. 114

|

no. 36

|

9577

Downloaded from https://www.pnas.org by 98.155.36.23 on October 19, 2026 from IP address 98.155.36.23.

Page 4

Figures and page graphics, page 4

soil properties of a pixel are not represented by the particular
soil pedons that were sampled within that pixel (34). Relatedly,
regions with low sampling density may be overly influenced by a
few data points that may not be representative of that region as
a whole. These last two issues are the major drivers of the spatial
distribution of model error (

SI Appendix

, Fig. S4

).

The HYDE dataset itself also presents a few limitations. First,

these land-use data are a combination of national and subna-
tional level statistics and remotely sensed land-use change (4),
which, in some cases, can create artificial changes in SOC along
political boundaries. This became particularly apparent when we
used the reconstructed land-use histories in prior centuries to
estimate SOC stocks for those time periods (Fig. 2). Second, only
very coarse land-use categories are represented, so management-
specific practices which can influence SOC levels (13), such as
tillage practices, rotations, and cover crops, are not represented.
Relatedly, HYDE does not contain direct information on forest
or wetland loss, both known drivers of SOC loss. Third, in our
model formation, there is no indication of the duration of a given
land use. SOC stocks, while often declining most rapidly in the
first decade after land-use change (35), often take many decades
to over a century to reach a new steady state (36, 37). Finally,
the HYDE dataset describes the extent of land use but not the
intensity. This limitation may be particularly important for the
grazing category, as SOC levels have been shown to decline with
increasing grazing pressure (38), although this effect appears to
be dependent upon grass species composition, with C3 grasses
showing large declines and C4 grasses showing small gains in
SOC with increased grazing pressure (39). It is very likely that,
taking these limitations together, our estimate of soil carbon debt
covers only a smaller fraction of the actual debt due to human
influence.

Implications.

This analysis indicates that the majority of the used

portions of planet Earth have lost SOC, resulting in a cumu-
lative loss of

∼

133 Pg C due to agricultural land use. These

SOC losses are on par with estimates of carbon lost from liv-
ing vegetation primarily due to deforestation (40) and are nearly
100 Pg C higher than earlier estimates of land use and land use
change-driven losses of SOC (41). Importantly, as Fig. 1 demon-
strates, there are hotspots of SOC loss, associated with extensive
cropping regions but also with highly degraded grazing land (

SI

Appendix

, Fig. S9

), suggesting that there are identifiable regions

which should be targets for SOC restoration efforts.

The potential to recover lost SOC may be more limited than

is often assumed. The amount of SOC that has been lost his-
torically can be thought of as the carbon sink potential of the
soil (42). Our analysis has found that this sink potential is

∼

133 Pg C (

SI Appendix

, Table S3

). A widely repeated figure

is that, with adoption of best management practices, two thirds
of lost SOC can be recovered (42). If the two-thirds figure is
accurate, then SOC sequestration has the potential to offset
88 Pg C (322 Pg CO

2

) of emissions. However, bottom-up esti-

mates of the maximum biophysical potential for carbon seques-
tration on cropping and grazing land range from 0.4 Pg C

·

y

−

1

to

1.4 Pg C

·

y

−

1

(20, 43). Assuming SOC reaches a new steady state

in 20 y (35, 44), this calculation suggests that 8 Pg C to 28 Pg C
can be recaptured. Even the range of 8 Pg C to 28 Pg C is likely
overly ambitious given the various social, economic, and techni-
cal constraints on universal adoption of best management prac-
tices (45), suggesting that the amount of the carbon sink that can
be filled is on the order of, at best, 10 to 30% globally and may
well be

<

10%.

Conclusions

Our data-driven statistical analysis confirms that agricultural
land use is a significant driver of SOC levels. Importantly, we
have generated estimates for the global cumulative loss of SOC

which potentially represent a maximum estimate of the SOC sink
capacity, and have demonstrated that there are hotspots of SOC
loss which are closely associated with land that has been identi-
fied as highly degraded. This analysis also demonstrated that not
all land use is associated with large losses in SOC, particularly in
regions with naturally infertile soils. These results provide a basis
for national and international policies to target SOC restoration
efforts but also suggest that more effort needs to be put into col-
lecting, integrating, and using legacy soil profile data, especially
historic data 50+ y old, so that even more reliable models of SOC
dynamics can be produced.

Materials and Methods

SOC varies in complex but mostly predictable ways across the landscape. It
can be best modeled as a function of climate, potential vegetation, topo-
graphic relief, soil parent material, and time (46, 47). Numerous spatially
explicit data layers now exist that cover most of these state factors of soil
formation. This is the foundation of the current state of the art in predictive
soil mapping (27, 48),

SOC

=

f

(

CL

,

R

,

L

,

LC

,

LU

)

,

[1]

where CL represents climate, R = relief, L = lithology, LC = land cover, and
LU = land use. Our interest here is in applying such soil

×

environment +

land-use models to historic land use while holding all other factors constant.
Once we have estimated a model for SOC (

m

) that relates SOC to environ-

mental + land-use data (using current data), we can use it to predict global
distribution of both SOC also for historical land cover/land-use data,

ˆ

SOC(

y

)

=

E

[

m

|

CL

,

R

,

L

,

LC

,

LU(

y

)

]

,

[2]

where

m

is the statistical model used that relates SOC to environmental

conditions and LU(

y

) are the historical land-use images for year

y

. By using

the HYDE v3.2 dataset (4), we have a spatially explicit representation of land
use [LU(

y

)] over the past 12,000 y, so that we can determine the difference

in SOC between the current and historical land use,

∆ ˆ

SOC

=

ˆ

SOC(

y

2

)

−

ˆ

SOC(

y

1

)

.

[3]

In this work, we have used the year 2010 to represent current conditions,
and, to estimate the cumulative impact of human land use on SOC, we have
replaced current land use with NoLU, which corresponds also to the preagri-
cultural times, for example, 10,000 BC.

Soil Profile Data.

ISRIC) (

www.isric.org

) curates the largest repository of spa-

tially explicit soil profile observations and samples (

N

= 158,147) covering all

climatic zones and biomes. These data come from both public and private
collections, and ISRIC has gone to great lengths to try to harmonize differ-
ences among different collections (23). Given that multiple depths are sam-
pled in each profile, there are over 850,000 measurements of SOC content
(in grams per kilogram), bulk density (in kilograms per cubic meter), and
coarse fragments (27). Based on these, we derive SOC density (in kilograms
C per cubic meter) for each horizon using

SOCD

=

OC

1000

·

BD

·

100

−

CF

100

,

[4]

where OC is SOC mass fraction in permilles, BD is soil bulk density in kilo-
grams per cubic meter, and CF is volumetric fraction of coarse fragments
(

>

2 mm) in percent. Unfortunately, BD is not available for many of the soil

horizons and needs to be determined either using global SoilGrids predic-
tions (27) or, for organic soils (C

>

8%), using a simple pedotransfer rule

relating BD to C, for example, from K ¨ochy et al. (49),

ˆ

BD

= [

1

.

38

−

0

.

31

·

log(C

/

10)

]

·

1000

.

[5]

After all missing values for BD have been filled, we can determine SOC den-
sity (SOCD) for every soil horizon.

Spatial Prediction Model for Organic Carbon Density.

We overlay the train-

ing points and environmental covariates and fit spatial prediction models
following the formula

SOCD(

xyd

)

=

d

+

X

1

(

xy

)

+

X

2

(

xy

)

+

. . .

+

X

p

(

xy

)

,

[6]

where

d

is the depth of observation and

X

p

(

xy

) are covariates. We model

spatial distribution of SOCD in three dimensions (soil depth used as a pre-
dictor) using all soil horizons layers at different depths. The derived spatial
prediction model is then used to predict SOCD at standard depths 0 cm,
30 cm, 100 cm, and 200 cm, so that the SOC stock for 0 m to 2 m can be

9578

|

www.pnas.org/cgi/doi/10.1073/pnas.1706103114

Sanderman et al.

Downloaded from https://www.pnas.org by 98.155.36.23 on October 19, 2026 from IP address 98.155.36.23.

Page 5

Figures and page graphics, page 5

SUSTAINABILITY

SCIENCE

derived as a cumulative sum of the four layers. Aggregated SOCD can be
converted to stock by using

SOCS

[

kg m

−

2

] =

SOCD

[

kg m

−

3

]

·

h

[

m

]

,

[7]

where

h

is the layer thickness in meters.

In addition to spatial prediction of OCD values, we also modeled associ-

ated spatial prediction uncertainty using the following procedure: (

i

) Run

fivefold cross-validation using the machine learning framework and derive
prediction residuals; (

ii

) overlay prediction residuals (absolute values) and

covariates and fit a new spatial prediction model using the same machine
learning framework; and (

iii

) apply this error model over the whole area of

interest to produce maps of absolute errors.

Environmental Covariates.

For modeling purposes, we use a large stack of

spatially explicit covariate raster data layers. These covariate layers have
been compiled from numerous sources representing the major state fac-
tors mentioned above, including the following: (

i

) The HYDE 3.2 Historic

land-use dataset (

ftp://ftp.pbl.nl/hyde/hyde3.2/

) (4) contains the distribution

of main agricultural systems from 10,000 BC (prehistoric NoLU condition)
to present time. Each raster layer represents the area (square kilometers)
of each pixel occupied by a given land-use category, with the 10 cate-
gories being as follows: total cropping, total grazing, pasture (improved
grazing land), rangeland (unimproved grazing land), total rainfed cropping,
and total irrigated cropping with further subdivisions for rice and nonrice
cropping systems for both rainfed and irrigated cropping. (

ii

) The Climate

Research Unit TS2.1 climatic surfaces for period 1960 to 1990 (

www.ipcc-

data.org/observ/clim/

) (50) include precipitation, ground frost frequency,

near-surface air temperature (daily minimum, daily maximum, and mean),
water vapor pressure, wet day frequency, cloud area fraction, near-surface
air temperature, and diurnal range. (

iii

) Current and historic forest extent

was estimated from the United Nations Environment Programme World
Conservation Monitoring Centre Generalised Original and Current Forest
cover map (

www.unep-wcmc.org/resources-and-data/generalised-original-

and-current-forest

). (

iv

) Topographic parameters were derived from Global

EarthEnv-Digital Elevation Model 90 (51) including slope, curvature, topo-
graphic index, topographic openness, valley depth, and multiresolution

valley bottom index. Topographic properties were derived using the Sys-
tem for Automated Geoscientific Analyses Geographic Information System
(52) at finer resolution (250 m) and then resampled to 10-km resolution. (

v

)

Landform and geologic substrate were determined from the United States
Geologic Survey Global Ecophysiography landform classification and litho-
logical map (

rmgsc.cr.usgs.gov/outgoing/ecosystems/Global/

) (53).

If not already available at 10-km resolution, spatial data layers were

resampled to 10-km resolution using the Geospatial Data Abstraction
Library software (54).

Statistical Modeling and Prediction.

The statistical modeling was accom-

plished using machine learning techniques implemented in R environment
for statistical computing (55). We used an ensemble prediction of two algo-
rithms: (

i

) random forest as implemented in the package ranger (56) and

(

ii

) gradient boosting as implemented in the package xgboost (57). For

model fitting, we used all soil profiles, then used this model to predict
SOCD at all grid nodes and all depths for current and historic land use. Dur-
ing the prediction, we derived ensemble estimates as a weighted average
between the two models, with weights estimated based on overall model
performance/accuracy. In most cases, random forest and gradient boost-
ing resulted in almost equal prediction accuracy; hence the weights are al-
most equal.

All computing was run on ISRIC High Performance Computing servers

with 48 cores and 256 GB RAM. Total computing time required to pro-
duce all outputs from scratch is about 18 h of optimized computing (or
about 1,000 central processing unit hours). All derived maps (as geo-
tiff files) and R code used in this analysis can be found at Woods Hole
Research Center’s github repository (

https://github.com/whrc/Soil-Carbon-

Debt

). Metaanalysis data can be found at

https://dataverse.harvard.edu/

dataset.xhtml?persistentId=doi:10.7910/DVN/QQQM8V

.

ACKNOWLEDGMENTS.

We thank Emily Cheney for compiling the remnant

native vegetation soil carbon database and Rebecca McCulley for contribut-
ing data to the remnant native vegetation soil database. Funding was pro-
vided by The Nature Conservancy and the Doris Duke Charitable Founda-
tion. ISRIC – World Soil Information is a nonprofit foundation primarily
funded by The Netherlands government.

1. Koch A, et al. (2013) Soil security: Solving the global soil crisis.

Glob Pol

4:434–441.

2. Amundson R, et al. (2015) Soil and human security in the 21st century.

Science

348:1261071.

3. Montanarella L, et al. (2016) World’s soils are under threat.

SOIL

2:79–82.

4. Klein Goldewijk K, Beusen A, Van Drecht G, De Vos M (2011) The HYDE 3.1 spatially

explicit database of human-induced global land-use change over the past 12,000
years.

Glob Ecol Biogeogr

20:73–86.

5. Ellis EC, Fuller DQ, Kaplan JO, Lutters WG (2013) Dating the anthropocene: Towards

an empirical global history of human transformation of the terrestrial biosphere.

Elem Sci Anthropocene

1:000018.

6. Food and Agriculture Organization of the United Nations and Intergovernmental

Technical Panel on Soils (2015)

Status of the World’s Soil Resources (SWSR) — Main

Report

(Food Agric Org, Rome).

7. Keesstra S, et al. (2016) Forum paper: The significance of soils and soil science towards

realization of the UN sustainable development goals (SDGS).

Soil Discuss

2016:

1–28.

8. Schlesinger W (2005)

Biogeochemistry,

Treatise on Geochemistry (Elsevier, New York),

Vol 8.

9. Crowther T, et al. (2016) Quantifying global soil carbon losses in response to warming.

Nature

540:104–108.

10. Tilman D, Cassman KG, Matson PA, Naylor R, Polasky S (2002) Agricultural sustainabil-

ity and intensive production practices.

Nature

418:671–677.

11. Oldfield EE, Wood SA, Palm CA, Bradford MA (2015) How much SOM is needed for

sustainable agriculture?

Front Ecol Environ

13:527–527.

12. Stockmann U, et al. (2013) The knowns, known unknowns and unknowns of seques-

tration of soil organic carbon.

Agric Ecosyst Environ

164:80–99.

13. Paustian K, et al. (2016) Climate–Smart soils.

Nature

532:49–57.

14. Guo LB, Gifford R (2002) Soil carbon stocks and land use change: A meta analysis.

Glob Change Biol

8:345–360.

15. Don A, Schumacher J, Freibauer A (2011) Impact of tropical land-use change on soil

organic carbon stocks—A meta-analysis.

Glob Change Biol

17:1658–1670.

16. Wei X, Shao M, Gale W, Li L (2014) Global pattern of soil carbon losses due to the

conversion of forests to agricultural land.

Sci Rep

4:4062.

17. Neill C, et al. (1997) Soil carbon and nitrogen stocks following forest clearing for

pasture in the southwestern Brazilian Amazon.

Ecol Appl

7:1216–1225.

18. Hoyle F, D’Antuono M, Overheu T, Murphy D (2014) Capacity for increasing soil

organic carbon stocks in dryland agricultural systems.

Soil Res

51:657–667.

19. Churchman GJ, Noble A, Bailey G, Chittleborough D, Harper R (2014)

Clay addition

and redistribution to enhance carbon sequestration in soils, Soil Carbon

(Springer,

New York), pp 327–335.

20. Lal R (2001) World cropland soils as a source or sink for atmospheric carbon.

Adv

Agron

71:145–191.

21. Pugh T, et al. (2015) Simulated carbon emissions from land-use change are sub-

stantially enhanced by accounting for agricultural management.

Environ Res Lett

10:124008.

22. Smith P, et al. (2016) Global change pressures on soils from land use and management.

Glob Change Biol

22:1008–1028.

23. Batjes NH, et al. (2017) Wosis: Providing standardised soil profile data for the world.

Earth Syst Sci Data

9:1–14.

24. Friedl MA, et al. (2010) Modis collection 5 global land cover: Algorithm refinements

and characterization of new datasets.

Remote sens Environ

114:168–182.

25. Scharlemann JP, Tanner EV, Hiederer R, Kapos V (2014) Global soil carbon: Under-

standing and managing the largest terrestrial carbon pool.

Carbon Manage

5:81–91.

26. Batjes N (2016) Harmonized soil property values for broad-scale modelling (wise30sec)

with estimates of global soil carbon stocks.

Geoderma

269:61–68.

27. Hengl T, et al. (2017) Soilgrids250m: Global gridded soil information based on

machine learning.

PLoS one

12:e0169748.

28. Intergovernmental Panel on Climate Change (2006)

2006 IPCC Guidelines for Nation-

al Greenhouse Gas Inventories, Prepared by the National Greenhouse Gas Inventories
Programme

, eds Eggleston HS, Buendia L, Miwa K, Ngara T, Tanabe K (Institute for

Global Environmental Strategies, Hayama, Japan).

29. Ogle SM, Breidt FJ, Paustian K (2005) Agricultural management impacts on soil

organic carbon storage under moist and dry climatic conditions of temperate and
tropical regions.

Biogeochemistry

72:87–121.

30. Nachtergaele F, et al. (2010)

Global Land Degradation Information System (GLADIS).

An

I

nformation

D

atabase for

L

and

D

egradation

A

ssessment at

G

lobal

L

evel Land

D

egradation

A

ssessment in

D

rylands (Food Agric Org, Rome), Tech Rep 17 Beta Vers.

31. Spera SA, Galford GL, Coe MT, Macedo MN, Mustard JF (2016) Land-use change affects

water recycling in Brazil’s last agricultural frontier.

Glob Change Biol

22:3405–3413.

32. Conant RT, Paustian K, Elliott ET (2001) Grassland management and conversion into

grassland: Effects on soil carbon.

Ecol Appl

11:343–355.

33. Fujisaki K, et al. (2015) From forest to cropland and pasture systems: A critical review

of soil organic carbon stocks changes in Amazonia.

Glob Change Biol

21:2773–2786.

34. Pachepsky Y, Hill RL (2017) Scale and scaling in soils.

Geoderma

287:4–30.

35. West TO, et al. (2004) Carbon management response curves: Estimates of temporal

soil carbon dynamics.

Environ Manage

33:507–518.

36. Sanderman J, Baldock JA (2010) Accounting for soil carbon sequestration in national

inventories: A soil scientist’s perspective.

Environ Res Lett

5:034003.

37. Senthilkumar S, Basso B, Kravchenko A, Robertson G (2009) Contemporary evidence

of soil carbon loss in the US corn belt.

Soil Sci Soc Am J

73:2078–2086.

38. Han G, et al. (2008) Effect of grazing intensity on carbon and nitrogen in soil and

vegetation in a meadow steppe in inner Mongolia.

Agric Ecosyst Environ

125:21–32.

39. McSherry ME, Ritchie ME (2013) Effects of grazing on grassland soil carbon: A global

review.

Glob Change Biol

19:1347–1357.

Sanderman et al.

PNAS

|

September 5, 2017

|

vol. 114

|

no. 36

|

9579

Downloaded from https://www.pnas.org by 98.155.36.23 on October 19, 2026 from IP address 98.155.36.23.

Page 6

Figures and page graphics, page 6

40. Houghton RA, et al. (2012) Carbon emissions from land use and land-cover change.

Biogeosciences

9:5125–5142.

41. Houghton R (1999) The annual net flux of carbon to the atmosphere from changes in

land use 1850–1990.

Tellus B

51:298–313.

42. Lal R (2004) Soil carbon sequestration impacts on global climate change and food

security.

Science

304:1623–1627.

43. Smith P, et al. (2008) Greenhouse gas mitigation in agriculture.

Philos Trans R Soc

Lond B Biol Sci

363:789–813.

44. Eggleston H, Buendia L, Miwa K, Ngara T, Tanabe K (2006)

IPCC Guidelines for

National Greenhouse Gas Inventories, Prepared by the National Greenhouse Gas
Inventories Programme

(Intergov Panel Clim Change, Hayama, Japan).

45. Smith P, et al. (2005) Carbon sequestration potential in European croplands has been

overestimated.

Glob Change Biol

11:2153–2163.

46. Jenny H (1994)

Factors of Soil Formation: A System of Quantitative Pedology

(Courier

Corp, North Chelmsford, MA).

47. McBratney AB, Santos MM, Minasny B (2003) On digital soil mapping.

Geoderma

117:3–52.

48. Scull P, Franklin J, Chadwick O, McArthur D (2003) Predictive soil mapping: A review.

Prog Phys Geogr

27:171–197.

49. K ¨ochy M, Don A, van der Molen MK, Freibauer A (2015) Global distribution of soil

organic carbon-part 2: Certainty of changes related to land use and climate.

Soil

1:367–380.

50. Harris I, Jones P, Osborn T, Lister D (2014) Updated high-resolution grids of monthly

climatic observations–the CRU TS3.10 dataset.

Int J Climatol

34:623–642.

51. Robinson N, Regetz J, Guralnick RP (2014) Earthenv-dem90: A nearly-global, void-

free, multi-scale smoothed, 90m digital elevation model from fused ASTER and SRTM
data.

ISPRS J Photogramm Remote Sens

87:57–67.

52. Conrad O, et al. (2015) System for automated geoscientific analyses (SAGA) v. 2.1.4.

Geosci Model Dev

8:1991–2007.

53. Sayre R, et al. (2014)

A New Map of Global Ecological Land Units—An Ecophysio-

graphic Stratification Approach

(US Geol Surv, Washington, DC).

54. Mitchell T, GDAL Developers (2014)

Geospatial Power Tools: GDAL Raster & Vector

Commands

(Locate Press, Chugiak, AK).

55. R Development Core Team (2009)

R: A Language and Environment for Statistical Com-

puting

(R Found Stat Comput, Vienna).

56. Wright MN, Ziegler A (2016) Ranger: A fast implementation of random forests for

high dimensional data in C++ and R.

J Stat Softw

77:18.

57. Chen T, Guestrin C (2016) XGBoost: A scalable tree boosting system. arXiv:1603.02754.

9580

|

www.pnas.org/cgi/doi/10.1073/pnas.1706103114

Sanderman et al.

Downloaded from https://www.pnas.org by 98.155.36.23 on October 19, 2026 from IP address 98.155.36.23.