Morais et al. (2019). Detailed global modelling of soil organic carbon in cropland, grassland and forest soils

3

]. Land use and land use

change (LU/LUC) are the second most relevant anthropogenic source of C into the atmo-
sphere, after emissions from fossil fuel combustion [

4

]. Long-term transformations from forest

and grassland to cropland are known to deplete soil organic carbon (SOC) stocks [

5

7

] releas-

ing C into the atmosphere. SOC can be replenished through C sequestration [

8

9

]. SOC is

strongly linked with soil management practices (e.g. mulching), soil proprieties (e.g. texture),
climate (e.g. temperature and rainfall) [

]. These factors display high spatial variability

] and make terrestrial C fluxes the most uncertain in the global C cycle [

13

]. Besides

PLOS ONE |

a1111111111
a1111111111
a1111111111
a1111111111
a1111111111

September 19, 2019

1 / 27

OPEN ACCESS

Citation: Morais TG, Teixeira RFM, Domingos T

(2019) Detailed global modelling of soil organic

carbon in cropland, grassland and forest soils.

PLoS ONE 14(9): e0222604.

10.1371/journal.pone.0222604

Editor: Debjani Sihi, Oak Ridge National

Laboratory, UNITED STATES

Received: March 6, 2019

Accepted: September 3, 2019

Published: September 19, 2019

Copyright:

©

2019 Morais et al. This is an open

access article distributed under the terms of the

Creative Commons Attribution License

, which

permits unrestricted use, distribution, and

reproduction in any medium, provided the original

author and source are credited.

Data Availability Statement: Data are available on

Zenodo (

http://doi.org/10.5281/zenodo.3387665

).

Funding: This work was supported by FCT/MCTES

(PIDDAC) through project UID/EEA/50009/2019

and by project “Animal Future - Steering Animal

Production Systems towards Sustainable Future",

funded by the Fundac¸

ã

o para a Ci

ê

ncia e

Tecnologia (SusAn/0001/2016). T.G. Morais was

supported by grant SFRH/BD/115407/2016 and R.

F.M. Teixeira by grant SFRH/BPD/111730/2015

from Fundac¸

ã

o para a Ci

ê

ncia e Tecnologia.

temporal and spatial variability, SOC measurement is also significantly uncertain. Even in
unperturbed systems, this means that changes in the soil C pool can take 5 or more years to be
detectable [

14

15

]. Linking SOC and LU is additionally difficult as global LU maps and data-

bases typically have a low number of classes or lack relevant information [

16

]. The European

Land Use/Land Cover Area Frame Survey (LUCAS) [

17

] database jointly reports SOC and LU

class for c. 19,000 field measurements that will continue to be repeated into the future. How-
ever, LUCAS does not include any information regarding how long after the last LUC event
the measurement was taken. The dynamics of SOC after LUC and the legacy effects of prior
LU are relevant for the determination of potential SOC stocks and how fast SOC is mineralized
[

18

].

In 2006, the Intergovernmental Panel on Climate Change (IPCC) proposed a method [

]

with three tiers of detail to account for changes in soil C stocks due to LU/LUC. This method
has been applied at local [

], regional [

] and global scales [

23

]. Tier 1 has default factors

of emissions and sequestration of LU/LUC. Tier 2 can potentially incorporate country-specific
management systems and biophysical data to calculate SOC stocks for mineral and organic
soils. Emission/sequestration C flows between atmosphere and soil are then calculated as the
yearly change in SOC stocks over 20 years. Nevertheless, Tier 2 approaches still rely heavily on
default factors for unspecific LU classes within major LU types (cropland, forest, pasture) [

],

as shown in the inventories to the United Nations Framework Convention on Climate Change
and submitted to the Kyoto Protocol by countries such as Canada [

24

], France [

25

] and Portu-

gal [

26

]. Consequently, comparisons of the IPCC method with regionalized models have

revealed large discrepancies regarding SOC changes [

The alternative for global and regional large-scale assessments of SOC change is to use pro-

cess-based models. Process-based soil models consider biogeochemical processes formulated
according to mathematical-ecological theory. They are capable of simulating SOC turnover
according to specific site conditions and relating it to management practices [

27

29

]. They

address user-defined temporal and spatial scales based on scenarios that characterize intra and
inter-annual dynamics. The application of these models transfers the need for highly regional-
ized data from the output variable (i.e. SOC) to the input variables. The Rothamsted Carbon
Model (RothC) [

] model is one of the most commonly used soil models today [

–

35

]. The

set of inputs required by RothC is a crucial advantage when compared with other process-
based soil models that require much larger datasets [

]. The range of dynamic SOC processes

captured by RothC is lower, but the fact that it is a parsimonious model makes RothC a good
candidate for global modelling exercises due to its manageability. However, few examples exist
of continental or global-scale applications of RothC [

,

36

], which had either low spatial reso-

lution or a low number of classes. For example, Gottschalk et al. [

] applied it globally to gain

insight on future SOC stocks under scenarios of LU distribution, but including only three LU
classes.

As the past few years have seen a surge in computational power and in global data availabil-

ity for modelling [

29

], here we propose new Tier 2-equivalent defaults obtained from global

SOC modelling with a modified version of RothC. The main innovative features of the mod-
elled data obtained here are (a) the level of spatial differentiation of the results (and respective
uncertainty), and (b) the number of land classes, particularly crop types. Regarding (a), the
basic unit of analysis were around 17,000 unique homogeneous territorial units (UHTU),
which result from the spatial combination of thermal zones, land cover, soil type, soil texture
and countries. Regarding (b), we calculated global SOC dynamics for 80 specific LUs within
broad LU classes (croplands, forests and grassland). The main results targeted with this study
were in particular a new set of highly spatially differentiated and land class-specific attainable
SOC stocks and C mineralization rates. Attainable SOC is the maximum attainable SOC with

Global modelling of soil organic carbon

PLOS ONE |

September 19, 2019

2 / 27

Competing interests: The authors have declared

that no competing interests exist.

occupation by a given vegetation type and climate conditions [

37

]. We calculate it as a present

potential due to the fact that we assume constant near-past climate and soil conditions, as well
as current crop yields, in its determination.

Materials and methods

The RothC model

The Rothamsted Carbon Model version 26.3 (RothC) is a model of C turnover in non-water-
logged soils [

]. It was initially developed to model C turnover in arable soils and later

expanded to grasslands and forests [

]. It takes into account the effects of temperature,

moisture content and soil type. The model uses a monthly step. SOC is divided in five com-
partments or pools, depending on decomposability: inert organic matter (IOM), easily decom-
posable plant material (DPM), resistant plant material (RPM), microbial biomass (BIO) and
humified organic matter (HUM). The IOM pool is resistant to decomposition and does not
receive C inputs [

38

]. Each compartment except IOM decomposes according to a first-order

decomposition process. RothC has been applied using data from long-term experiments across
several ecosystems, climate conditions, and land use (LU) classes [

].

Calculation procedure

The calculation procedure is depicted in

Fig 1

and summarized as follows. The details for each

of these steps are presented in the ensuing sub-sections. First, we collected the source data for
all input variables presented in

Fig 1

(i.e. climate and management data, initial SOC stock and

clay content) and respective uncertainty, assuming normal distributions. These data involve all
climate variables (temperature, precipitation and evaporation), management variables (water
input, residues production per LU and country as function of the yield, farmyard manure and
organic fertilizer application per LU) and initial SOC stock and clay content. Then we defined
c. 17,000 Unique Homogeneous Territorial Units (UHTUs) that are assumed to be similar in
terms of soil type, climate regime and current LU. We assigned all input variables to each
UHTU. We then selected priority LU classes and a set of management practices for cropland
production. In each UHTU only a subset of the list of LU classes is biophysically feasible. The
next step was to insert all variables into a MATLAB version (produced in this paper) of RothC.
We ran model initialization to divide the initial SOC stock under current land use into the five
C pools using the method proposed by Weihermu

¨ller et al. [

39

]. Then, we used RothC one

hundred times in each UHTU for each feasible LU class and for 86 years, in order to simulate
results that take into consideration intra-UHTU variability and obtain uncertainty measures
for average UHTU-level results. Results were then used to fit a two-parameter saturating expo-
nential model for each UHTU. We calculated attainable SOC stocks and mineralization rates
per LU class and UHTU using those two parameters. Finally, we compared results to long-
term SOC measurements and data collected in international projects for a wide range of geo-
graphical regions.

Input data.

Fig 1

indicates the land, soil and climate data inputs required by RothC [

].

Crop parameters involved in calculations are soil cover period, monthly input of plant residues
(t C/ha), farmyard manure (stored and sprayed manure/slurry and direct deposition of excre-
tions during grazing) C inputs (t C/ha) and irrigation (mm). The clay content of soils (%) is
used to determine the initial distribution of SOC between pools. The starting SOC content (t
C/ha) is used to initialize the model and also the CO

2

flux from the soil in the first simulation

period. As for climate data, the model requires average monthly mean air temperature (˚C),
precipitation (mm) and open pan evaporation (mm). We considered a “static world” scenario,
i.e. there was no change in climate in the future. This provides potential SOC dynamics in

Global modelling of soil organic carbon

PLOS ONE |

September 19, 2019

3 / 27

near-past conditions but cannot be extrapolated for forecasting the role of climate change in
SOC. We assigned uncertainty to all input data using mean values and standard deviations.
We assumed a normal probability distribution for all parameters.

Soil cover period is required due to higher decomposability of C in uncovered soils [

].

This is a binary monthly variable, where 1 means that the soil was covered with vegetation dur-
ing that month and 0 means that the soil was uncovered. Soil cover period was set to 1 for per-
manent crops (orchards, olive groves and vineyards), forests, grasslands and shrublands. For
croplands, this parameter was obtained from Chapagain et al. [

40

] and depends on the thermal

zones used by the Global Agro-ecological Zones (GAEZ) Project [

41

], whose definition is

based on temperature and precipitation. They provide sowing/planting dates and the duration
of the vegetation cycles. In the months between sowing/plantation and harvesting we consid-
ered the soil cover parameter to be 1 and in the other months, in the case where residues are
not left on the field and so the soil is bare, we assigned it the value 0. For cereals, when residues
are left on the field we used the value 1 for all months.

To estimate annual plant residue C input, we used the IPCC methods [

,

]. Perma-

nent and annual crop residues are equal to a fraction of above-ground crop yield plus the
below-ground residues [

https://doi.org/10.1371/journal.pone.0222604.g001

]. Below-ground residues are also a fixed fraction of the above-

ground yield. Yield data at country-level was obtained from the Food and Agriculture Organi-
zation of the United Nations (FAO) [

44

]. We used the average and standard deviation of yields

between 2004 and 2014. For forests and grasslands, residue inputs are the biomass left to decay
on the ground (e.g. litter), which depends on the vegetation type and climate region. The full

Fig 1. Workflow procedure for data preparation to be implemented in the RothC model. LU–Land use; SOC–Soil organic carbon; UHTU—Unique Homogeneous
Territorial Unit.

Global modelling of soil organic carbon

PLOS ONE |

September 19, 2019

4 / 27

method description and data used for crop and forest and grassland is in

. The method

used for croplands determines C residues for an entire year. We then distributed the residues
monthly. This distribution depends on the distribution of net primary production per month
and the life stages of plants. We used the method proposed by Jebary et al. and Morais et al.
[

,

] to distribute residue C input to soil between months. According to this method, for

cereal crops, 50% of residues were allocated to the harvesting month and the remaining 50%
were equally distributed among the three months before. For permanent crops, 70% of resi-
dues were allocated to the pruning months and the remaining 30% to the prior four months.
This distribution depends on the monthly distribution of net primary production and life
stages of plants. The months for harvesting and pruning were obtained from Chapagain et al.
[

40

]. For forests and grasslands, residues were equally distributed during the year, thus assum-

ing that there are no management activities.

Farmyard manure is an extra C input in croplands and grasslands. In grasslands, the source

of this input is both direct dung deposition during grazing and (possibly) application of
manure, while in croplands only the latter is applicable. For dung excretions during grazing
we used livestock units on grasslands from Gridded Livestock of the World v2.0 [

47

] combined

with manure excretion from IPCC [

] for cattle, sheep and goat livestock. For manure appli-

cation in grass and croplands there is no known data available. Organic fertilization is highly
variable, as it depends on local manure availability and price, its source and type (solid, liquid)
and its chemical composition (humidity, C:N ratio, etc.). We used data from Mueller et al. [

48

]

for total aggregate consumption of nitrogen (N) fertilizers from synthetic and animal origin,
and considered three scenarios. In the first scenario, the total consumption of N from Mueller
et al. [

48

] was assumed to be synthetic fertilizer only, and thus there is no manure input. In the

second, we considered that only 50% of the total amount of fertilizer was applied as manure.
In the third scenario, we considered that all N used as fertilizer was manure (organic fertilizer).
To convert the N content of manure to C, we used the C:N ratio of farmyard manure provided
by FAO [

49

]. We excluded any vegetable material applied as part of manure (such as straw and

other materials from animal beddings) to avoid double-counting, as plant-based C input is
already included in cases where residues are left on the field. The detailed method is in

.

RothC does not distinguish between water input types (i.e. water from irrigation and water

from precipitation). Water provision during dry seasons can make a difference in SOC miner-
alization [

38

]. We used the method of Pfister et al. [

50

] to calculate water input as the sum of

precipitation and irrigation. The method assumes that, for each crop and attainable yield, the
water requirements of that crop must be attained through precipitation alone or with addi-
tional irrigation. Irrigation is either zero (if precipitation is higher than the water needed for a
given month) or the difference between the water needed by the crop and the precipitation for
each month. Crop water requirements were obtained by multiplying specific crop coefficients
depending on a monthly growth stage factor and potential evapotranspiration. The growth
stage factors were obtained from Chapagain et al. [

40

] for all crops in each thermal zone.

Potential evapotranspiration was obtained using the Thornthwaite equation [

51

], which uses

monthly average air temperature, average day length, in hours, and number of days per month
obtained from MODIS [

52

]. The detailed method description and data used for croplands is in

Average monthly precipitation and mean air temperature were both obtained from

National Aeronautics and Space Administration’s (NASA). Precipitation was obtained from
the database of the “Global Precipitation Climatology Project (GPCP)” [

53

] We used data for

thirteen years (2000–2013). To establish normal climatic conditions, 30 years are usually nec-
essary. However, the data from NASA does not go so far back. Additionally, given recent
changes in climate that produce many anomalous meteorological years, we believe a shorter

Global modelling of soil organic carbon

PLOS ONE |

September 19, 2019

5 / 27

period can nevertheless more accurately represent near-past conditions [

54

56

]. Further, we

also used this period to ensure consistency with other data, as yield data is also valid for the
same period only. Monthly open pan evaporation was calculated as 75% of the potential evapo-
ration, an assumption used by previous studies using RothC [

46

]. Potential evaporation was

calculated using the Thornthwaite equation [

51

].

SOC initialization values were obtained from the European Soil Data Centre (ESDAC) [

57

].

Clay content was obtained from the Harmonized World Soil Database [

58

]. The topsoil layer

considered was 30 cm deep.

Definition of unique homogeneous territorial units. In this study, we used the concept

of UHTU to establish the basic units of analysis. UHTUs are geographical regions where soil
type and texture, climate type and LU are uniform at the scale of analysis. RothC was used for
each unique region individually. UHTUs were defined as a geographical combination of five
layers, namely thermal zones, land cover, soil type, soil texture and country. Thermal zones
were obtained from GAEZ [

41

], which divides the world into 12 zones (

). We excluded

arctic and desert regions. We used data from the Land Processes Distributed Active Archive
Center (LP DAAC) [

59

] for attribution of land classes. This source considers 16 LU classes,

out of which we excluded “water” and “unclassified” regions from the analysis. Soil type and
texture data were obtained from the World Reference Base (WRB), as depicted in Fischer et al.
[

60

]. We considered the 8 types of soil and 13 types of texture from WRB. Finally, it was neces-

sary to use a country layer in the definition of UHTUs because crop yield data from FAO is
only available at country-level. We obtained the country borders from the World Borders
Dataset [

61

]. We used FAO yield data for 202 countries and for the period 2004–2014. The full

list of countries is in

. The combination of these data sets resulted in 17,203 UHTUs (

S3

File

). An individual set of RothC simulations was then carried out for each UHTU. The aver-

age area of each UHTU is 6,400 km2.

Choice of land use classes. We performed simulations of SOC dynamics for 43 LU clas-

ses, including 28 agricultural classes (under different management practices), 16 forest classes
and 1 grassland class. The agricultural classes selected were the 28 most produced and traded
(which are the result of intersecting the top twenty most produced and traded products) dur-
ing the period 2004–2014 (most recent year available) [

44

] in the World (using FAOSTAT

nomenclature): apples; bananas; barley; cabbages and other brassicas; carrots and turnips;
cocoa, beans; coconuts; coffee, green; grapes; groundnuts, with shell; maize; oil, palm fruit;
olives; onions, shallots, green; oranges; potatoes; rapeseed; rice, paddy; seed cotton; sorghum;
soybeans; sugar beet; sugar cane; sunflower seed; sweet potatoes; tobacco, unmanufactured;
tomatoes; and wheat. Regarding management practices, besides the three organic fertilization
scenarios explained previously, we considered rainfed and irrigated production (except cab-
bage, carrot and onion, in which cases we only used irrigated production because rainfed areas
are insignificant compared with irrigated areas). We considered two management options for
residue management of cereals: (1) residues are left on the field; (2) residues are removed from
the field. Forest and grassland classes were obtained from the IPCC classification [

broadleaf deciduous forest and needleleaf evergreen forest in the dry boreal zone, moist boreal
zone, moist cold temperate zone, dry cold temperate zone, dry warm temperate zone, moist
warm temperate zone, sub-tropical zone and tropical zone. The IPCC method [

43

] for forest

residues provides annual C residues and transition time for all forests. Thus, we used the linear
annual litter increase rate during forest growth and after maturity a constant annual C residue
input. IPCC considers only one type of grassland. In total there are 80 possible combinations,
combining cropland with management practices (irrigation/rainfed and residues left/removed
on/from the field) plus the forest classes and the grassland class. To determine the UHTUs
where each agricultural LU class is feasible, country-level FAOSTAT yield data is insufficient.

Global modelling of soil organic carbon

PLOS ONE |

September 19, 2019

6 / 27

Instead, we used production capacity maps from GAEZ [

41

]. These maps integrate climatic

and edaphic variables to determine feasible areas for each agricultural LU class. For each crop,
if there is at least one pixel inside of the UHTU with non-null production capacity in the pres-
ent, we assumed that the crop is feasible in that UHTU. Forest classes are divided by thermal
zone by definition. For grasslands, we used directly the areas classified as “grassland” from LP
DAAC maps [

59

].

MATLAB implementation. RothC can be obtained under a free for research license from

the Rothamsted Research website [

62

]. However, that version of the model only runs for one

region (UHTU) and one LU class at a time. To expedite calculations for all UHTUs simulta-
neously we implemented a Matlab version of RothC. We implemented also a Monte Carlo
approach [

63

] for determining the uncertainty of results. We ran the model 100 times for each

LU class in each UHTU. Each iteration thus used a unique set of input parameters drawn ran-
domly from their respective normal probability distributions.

Calculation of attainable SOC stocks and mineralization rates. We obtained the annual

SOC stock from the model for 86 years (simulation between 2014 and 2100). As 2000–2013 cli-
mate data was used, this estimation considers climate stability and does not take the role of cli-
mate change into account.

RothC is a multiple pool model but the overall shape of the dynamic SOC curves for each

pool is approximately exponential. For grassland and agricultural classes, this approximation
is a good fit for all periods simulated because the C input from residues is constant, but for for-
est classes it is only true after maturity (i.e. when annual C residues are constant). In the transi-
tion period, dynamic SOC curves are approximately a fourth-degree polynomial. For
grasslands, croplands and mature forests, we used a simple two-parameter mass balance
model for SOC dynamics to fit the annual SOC stocks obtained by the application of RothC
for each land class in each UHTU. The SOC balance is the difference between C input and
mineralization, described by

dSOC

dt

¼

K

a

SOC

ð

1

Þ

where SOC is the SOC stock (t C/ha) at time t, K is the C input to soil at time t, and

α

is the C

mineralization rate. C input (parameter K) is a function of the LU. Therefore, K and

α

are

time-invariant and vary with location (i.e. UHTU) and LU class. Integrating between 0 and t,
we obtain

SOC

t

¼

K

a

1

e

a

t

ð

Þ þ

e

a

t

SOC

0

ð

2

Þ

SOC is limited by an upper bound, i.e. SOC reaches a maximum attainable equilibrium,

which is given by K/

α

. The curve fitting procedure thus provides attainable SOC stocks and

mineralization rates that are LU and site-specific.

Comparison of results and assumptions with data and previous uses of the
model

Attainable SOC stocks calculated in this paper are difficult to validate for two main reasons.
First, because they are regional averages valid at UHTU scale, which are relative large and
approximately homogeneous areas. Our results should therefore be an approximation of mea-
surements made at point/plot scale within UHTUs, but not an exact match. Second, they are
difficult to compare with benchmarks because they are potentials rather than actual estimated

Global modelling of soil organic carbon

PLOS ONE |

September 19, 2019

7 / 27

SOC stocks. Comparison, for example, with available global SOC maps would be incorrect as
those depict current SOC levels rather than potentials.

To assess the agreement between our results and field-level measurements of SOC, we

therefore used a mixed strategy of indirect verification. Namely, we triangulated multiple lines
of evidence to inquire whether literature or database sources lend support to our main results.
This analysis included three separate parts.

First, we compared SOC stocks and mineralization rates obtained here with multiple field

measurements taken from multiple croplands and forests in different locations of the world.
In the absence of crop-specific global meta-analyses of SOC stock variation with LUC, we
compare our results with prior local/regional applications of RothC. We compare results from
those in-situ experiments to our results for the same regions, to assess if our approach, which
is global and not tailored for farm-level assessments, can nevertheless adjust well to observa-
tions. However, to avoid comparing potentials (our work) with current SOC levels, we used
only results from long-term experiments where, in principle, SOC stocks were closer to maxi-
mum attainable.

In this first step we used the search engine Google Scholar and the keywords “long-term

experiment” + “soil organic carbon” and “long-term experiment” + “mineralization”. The
search was conducted on July 30, 2018. We obtained 4,480 and 610 matches, respectively. We
then selected only peer-reviewed studies published in books and journals with impact factor
according to the 2017 InCites “Journal Citation Reports” from Clarivate Analytics [

64

]. We

excluded project reports, Master and PhD dissertations and other grey literature sources. We
also excluded experiments involving crop rotations or other LU changes within the study
period, as these are not included in our analysis. As the goal was not to perform a comprehen-
sive meta-analysis, we considered only the first 10 Google Scholar pages, i.e. 200 papers. We
only used studies that explicitly report SOC stock (t C/ha) or SOC concentration (g C/kg soil)
and bulk density (g/cm

3

) (to enable conversion of SOC concentration into stock). Using these

criteria, we obtained a final list of 12 papers (out of 200).

The second comparison we made in this paper involved the model parameters chosen in

local/regional and global applications of RothC regarding the DPM/RPM ratio and the decom-
position rate of each pool. RothC was calibrated using SOC data from long-term experiments
by the original developers of the model [

]. We used the default parameterization of the

model for the global application. As parameters could be UHTU-specific, we checked whether
other applications of the model have used non-default sets of parameters and, if so, if they
were significantly different from the default. We again used the search engine Google Scholar
and the keywords "Rothamsted Carbon Model" + "DPM" OR "RPM" OR "Hum" OR "Bio" OR
"IOM". The search was conducted on December 2, 2018. We obtained 521 matches. We
selected only peer-reviewed studies. Papers that used older versions of RothC than version
26.3 [

] and papers that only mentioned RothC without applying it were also excluded. Here,

studies involving crop rotations were included. We included local/regional applications
(mainly with measured soil C inputs) of RothC and continental/global applications (with esti-
mated and aggregated data for soil C inputs). As the goal was not to perform a comprehensive
meta-analysis, we considered only the first 50 papers (out of 521). Those 50 papers include a
total of 233 applications of the RothC to particular points/regions. The full list of references, as
well as parameters used in each of them, is included in

S2 File

The final step of this section of the work was a comparison of attainable SOC stocks

obtained in this study with SOC concentration measurements from the LUCAS Project [

17

].

In this project SOC concentrations were measured at the 0–20 cm topsoil layer. To compare
results, we assumed a uniform distribution on SOC in the topsoil profile, which means that 2/
3 of SOC is in the 0–20 cm layer. Finally, to convert SOC concentrations (kg C/kg soil) to SOC

Global modelling of soil organic carbon

PLOS ONE |

September 19, 2019

8 / 27

stock (t C/ha), we used bulk density maps produced in the same project [

65

]. These site mea-

surements covered 1,171 European UHTUs. As LUCAS assessed current SOC stocks, we
could not compare values directly. Instead, we assume that higher present SOC should be, on
average, correlated with higher potential SOC within the same land classes. We performed
geospatial correlation analysis using Pearson’s r using 21 LU classes whose classification
matches the system used in LUCAS. We excluded 20 LU classes because they had less than 100
observations each. To calculate these correlations, we used software SPSS version 23. We cre-
ated a correspondence between the classification used in this study and the LUCAS Project
classification, which is described in

Results and discussion

Dynamics of SOC loss and accumulation

We ran RothC for 86 years (from 2014 to 2100) and simulated SOC dynamics in all LU classes
that are biophysically feasible in each of the c. 17,000 UHTUs, starting from present-day mea-
sured SOC stocks.

Fig 2

displays selected results for the four largest UHTUs. SOC stocks (0–30

cm depth) are highly LU class-dependent even for classes that are traditionally aggregated
such as different types of arable land. For example, in the region along the northern border of
the United States of America (USA) with Canada the SOC stock after 86 years of irrigated
maize with residues left on the field is 89

±

3 t C/ha while for irrigated soybean it is 56

±

2 t C/

ha (

), despite the fact that both are examples of typical agricultural classes. For each LU

class, the SOC stock is highly dependent on the region under consideration: for example the
average SOC stock for maize in Northwest India (

Fig 2D

) after 86 years is about 10

±

1 t C/ha,

which is considerably lower than the previously cited stock in North America (

).

Forest LU classes such as “broadleaf deciduous forests” in North America (

) and Asia

(

Fig 2C

) need more time to reach the maximum attainable level. Immediately after transition

to forest, SOC stocks decrease due to the low input of C residues as the forest is starting to
grow. When forests reach a certain level of biomass stock, SOC begins to increase (

and

Fig 2D

) as observed in other studies [

5

]. Climatic features such as temperature and precipita-

tion and different soil types drive differences in soil C mineralization and also feasible yields,
which in turn drives plant C inputs to soil.

In general, agricultural LU classes reach lower SOC levels than forests at maturity. Some

boreal forests in western Scandinavia are exceptions where croplands can have higher SOC
stocks than forests. Potential C input from residues kept on the field is high in these regions
(e.g. tomato and maize), while boreal forests have low tree biomass and productivity and thus
low soil C input [

66

]. Additionally, conversions to agricultural classes typically reach new SOC

equilibria quicker than transitions to forest and grassland LU classes. Consequently, SOC loss
is typically faster than SOC recovery. Grasslands generally reach higher SOC levels than crop-
lands (

Fig 2B

and

Fig 2D

) due to plant shoot and roots and animal C inputs to soils [

with notable exceptions in some regions for agricultural classes with high plant residues, such
as irrigated maize (

and

shows this effect clearly. When residues are

removed from the field, croplands almost always reach lower SOC stocks, but if residues
remain on the field there are large areas of North America, Western Europe, Central Asia and
Australia where croplands accumulate more SOC than grasslands. Two effects justify this
result–residues increase C input to soil, and soil cover with those residues reduces C minerali-
zation. Out of the two, plant residue C input is the key explanatory variable of the results, but
it is also the main source of uncertainty.

In order to test the significance of results in

, we performed ANOVA for all UHTU and

LU classes where attainable SOC for croplands were higher than attainable SOC in grasslands.

Global modelling of soil organic carbon

PLOS ONE |

September 19, 2019

9 / 27

On average, and for each LU class, the difference of the means is statistically significant at 5% in
30% of the UHTUs (15–67% depending on the agricultural LU class). Irrigated wheat with resi-
dues left on the field was the LU class with the highest number of UHTUs with p-value lower
than 5% (165 out of 673 UHTUs), while rainfed sorghum with residues left on the field was the
LUC class with the lowest number (2 out of 3). Irrigated maize with residues left on the field was
the LU with the highest number of UHTU where attainable SOC was higher (625 UHTU), but
the difference is only statistically significant in 20% of those UHTUs. These tests reassure statisti-
cally that for a fraction of the cases identified in

https://doi.org/10.1371/journal.pone.0222604.g002

some croplands indeed attain higher SOC.

Global distribution of long-term SOC equilibria and mineralization rates

To aid in interpretation, we fitted an exponential model to each 86-year SOC time series,
depicting inter-year changes in SOC stocks as a simple balance between fixed C inputs to soil

Fig 2. Soil organic carbon (SOC) stock (t C/ha) in four land use classes, in the scenario without manure application, in four unique territorial units studied (dark grey),
namely: a, 911 (North America). b, 11867 (South America). c, 3499 (Europe). d, 7916 (Asia). Marks in the graphs represent mean annual SOC stock obtained in RothC,
lines represent fitted exponential curves and the shaded area represents the 95% confidence interval of each fitted year.

Global modelling of soil organic carbon

PLOS ONE |

September 19, 2019

10 / 27

and fixed rate C mineralization and emission as CO

2

through soil respiration; this is effectively

the most simplified version of RothC possible, considering a single C pool in the soil and a
mean rate of mineralization for the average pool. The upper limit of the exponential model is
the maximum SOC stock obtainable, or attainable SOC. The full list of Curve fitting results
(i.e. long-term SOC stock, mineralization rates and the parameter K) are presented in

S4 File

.

depicts the global distribution of attainable SOC stocks for six particular LU classes,

while

shows C mineralization for the same classes. Note that none of these classes exist in

all UHTUs, because of biophysical feasibility constraints. Irrigated tomato production, for
example, is feasible in 11,985 UHTUs (the highest number for all LU classes), and rainfed oil
palm in 73 UHTUs (the lowest number).

Regardless of LU class, high latitude regions and (for some land classes) Equatorial

regions display higher attainable SOC stocks (

), but for different reasons. Equatorial

regions are marked by high temperature and high rainfall, causing higher soil respiration
(

) through the effect on soil microbial activity [

69

https://doi.org/10.1371/journal.pone.0222604.g003

70

], but this effect is countered by

Fig 3. Visual representation of the regions where soil organic carbon (SOC) stock (t C/ha) is higher for grasslands
than agricultural land uses (in green) or lower (in red). Agricultural land uses and practices represented: a, b, c, d–
maize, e, f, g, h–wheat; a, c, e, g–irrigated, b, d, f, h–rainfed; a, b, e, f—residues left on the field, c, d, g, h—residues
removed from the field.

Global modelling of soil organic carbon

PLOS ONE |

September 19, 2019

11 / 27

higher yields and consequently higher plant C input to soil. High latitude regions have high
water availability and low temperature, and therefore low mineralization rates (

https://doi.org/10.1371/journal.pone.0222604.g004

). The

Mediterranean and Australian semi-arid regions show high mineralization rates, but attain-
able SOC stocks are higher along the Mediterranean due to higher plant C inputs to soils
(e.g.

Fig 4D

and

Fig 5D

).

Coniferous forests in northern boreal regions (

Fig 4E

) have the highest attainable SOC out

of all types of forest anywhere in the World. Partially this is the result of climatic conditions in
the north boreal regions, but also of high plant C inputs. Southern boreal forests have lower
attainable SOC than northern regions as mineralization rates increase towards the Equator.
Despite lower mineralization in high latitude regions (boreal thermal zones), these forests in
the Southern border tend to have lower attainable SOC stock compared with temperate forests
due to higher annual litter production of the latter. In UHTUs where boreal coniferous forest
and temperate continental forest can coexist, the latter have on average 75 t C/ha higher attain-
able SOC.

Fig 4. Attainable soil organic carbon stock (t C/ha) obtained by curve fitting per unique homogeneous territorial unit, in the scenario without manure
application. a, irrigated maize, residues left on the field. b, irrigated soybean. c, irrigated sunflower. d, irrigated tomato. e, broadleaf deciduous forest in dry warm
temperate zone, and. f, grassland.

Global modelling of soil organic carbon

PLOS ONE |

September 19, 2019

12 / 27

Regarding grasslands, the input of C from dung deposition depends on grazing intensity,

which introduces a high degree of spatial variability in attainable SOC. Lower SOC stocks are
found in arid and dry areas of Oceania (particularly Australia). In general all LU classes in this
region have low attainable SOC stocks (

) and high mineralization rates (

https://doi.org/10.1371/journal.pone.0222604.g005

). Grasslands

have particularly high SOC stocks in equatorial regions (Sub-Saharan regions, Central Amer-
ica and Southern Asia regions), which is 7 t C/ha higher than in croplands on average. How-
ever, croplands with high residue production crops (e.g. maize, wheat) in mid and high
latitude regions, can have stocks up to 20 t C/ha higher than grasslands.

For maize (

Fig 4A

) there is a downward North-South gradient. In North Canada and the

USA SOC stocks are higher than in southern USA and Mexico, although moving towards the
South of the continent SOC increases slightly for some crops and grasslands. Typically, attain-
able SOC stock in Africa is lower than in the other regions due to high turnover rates of
organic material combined with low yields [

71

]. Irrigated sugarcane is the technically feasible

crop most found in African UHTUs, but it is also the agricultural LU class with the highest

Fig 5. Average organic carbon mineralization rate (% per year) obtained by curve fitting per unique homogeneous territorial unit, in the scenario without
manure application. a, irrigated maize, residues left on the field. b, irrigated soybean. c, irrigated coconut. d, irrigated tomato. e broadleaf deciduous forest in dry warm
temperate zone, and. f, grassland.

Global modelling of soil organic carbon

PLOS ONE |

https://doi.org/10.1371/journal.pone.0222604.g006

September 19, 2019

13 / 27

mean mineralization rate (about 5% per year). Tropical and subtropical forest classes have
lower mean SOC mineralization rates (about 1% per year) than arable land classes.

Fig 6

depicts the impact on attainable SOC stock of management practices in grain maize

production. We tested: a) the effects of irrigation against rainfed maize, b) leaving plant resi-
dues on the field against exporting crop leftovers, and c) organic fertilization. The distribution
is regionally asymmetrical, but in general leaving residues on the field (

Fig 6B

) produces higher

increases of attainable SOC stock than irrigation (

Fig 6A

) and organic fertilization (

Fig 6C

).

The average effect of irrigation on all agricultural classes is an increase in attainable SOC of
4

±

2 t C/ha. Irrigation has less effect on permanent crops (0.9

±

0.1 t C/ha). Wheat is the agri-

cultural class where irrigation has a higher positive effect (10

±

5 t C/ha). SOC increased the

most with irrigation in temperate thermal zones. For example, the SOC increase due to irriga-
tion at temperate oceanic regions is 3–8 t C/ha for UHTUs. On average organic fertilization
increases the attainable SOC stock of agricultural classes by 3

±

1 t C/ha. Potato (irrigated and

rainfed) is the class where organic fertilizer application increases attainable SOC stock the
most (about 8

±

3 t C/ha). Temperate thermal zones are where organic fertilizer application

increases SOC stocks the most, between 4 and 6 t C/ha. Maintenance of residues on the field
increases attainable SOC stock for cereals classes by 12 t C/ha. The increase is highest (18

±

4 t

C/ha) for maize and wheat.

Sugar crops (sugarcane and sugar beet) and tomato have the highest attainable SOC stocks

among crops, due to higher yields and production of residues. These crops are estimated to
produce residues that are twice the yield. Others crops, like tobacco, also produce more-than-
proportional residues but have low yields and as such do not accumulate as much SOC. Irri-
gated sugar beet leads to, on a global average, 20 t C/ha higher stabilisation SOC compared
with irrigated maize with residues left on the field. In this case, although maize is more produc-
tive, the fraction that is plant residue is much lower.

Fig 6. Soil organic carbon (SOC) stock difference (t C/ha) for the effect of management practices on attainable soil organic carbon stock for maize. a, SOC for
irrigated maize minus SOC for rainfed maize. b, difference in SOC for maize with residues left on the field and with full extraction of residues. c, difference in SOC for
100% organic fertilization of maize and 0% organic fertilization.

Global modelling of soil organic carbon

PLOS ONE |

September 19, 2019

14 / 27

Additional management practices such as tillage can potentially affect SOC stocks, but are

not explicitly modelled by RothC. In the literature, the effect on SOC stocks of replacing con-
ventional tillage with no-till is on average 2–10 t C/ha [

72

,

73

]. According to the IPCC [

],

tilled cropland soils have 2–15% less SOC than no-till soils. This is less than half (31%) of the
difference between attainable SOC stocks of grasslands and cereals found here. Our conclu-
sions about regions where croplands reach higher SOC stocks than other land uses are there-
fore not affected by omitting the effects of tillage.

Comparison of results with data

Global modelling exercises provide highly uncertain results, in particular for SOC, which is
highly sensitive to geography and dependent of LU class. Our application of RothC provides
results depicting trends at UHTU level. At higher scales, results are well aligned with prior litera-
ture, despite high variation due to the crop-specificity of our results. For example, in

Fig 2C

the

SOC gain due to transformation to forests depicted for temperate zones is about 200%. Prior
research estimated SOC stock gains between 75% and 200% for transitions from cropland to for-
est (after 100 years), also for temperate zones [

5

]. Nevertheless, despite the high level of regionali-

zation, these estimates are not necessarily farm-specific and are not applicable locally. This raises
the question of whether the outputs of the model realistically depict SOC as measured.

In general, attainable SOC stocks and mineralization rates reported here are similar to values

reported in long-term field measurements for specific sites reported in other studies. In the
Huang-Huai-Hai Plain of China [

74

] the long term SOC stock (without organic fertilizer appli-

cation) at 20 cm depth for maize production is 18 t C/ha (27 t C/ha, assuming homogenous dis-
tribution in the first 30 cm depth). Our results for the same region were 20

±

2 t C/ha (30 cm).

In another Chinese region (Jiangxi) [

75

] long-term maize production without fertilization

results in a SOC stock of 20 t C/ha, which changes to between 23 and 31 t C/ha with different
fertilization rates. For the same region we obtained 25

±

5 t C/ha and 38

±

8 t C/ha, respectively.

In a Canadian region [

76

] SOC in wheat fields stabilizes at 35 t C/ha (at 15 cm depth– 70 t C/ha

at 30 cm) when straw is left on the field, while we report 51

±

10 t C/ha for the same case and

region. In the Harpenden region (United Kingdom) [

] SOC stocks of about 25 (at 15 cm

depth– 50 t C/ha at 30 cm) and 60 t C/ha (at 23 cm depth– 78 t C/ha at 30 cm) were reported
for arable land and mixed deciduous forest following a calibration procedure of RothC using
field-level data. We obtained for the same region 35

±

9 and 66

±

15 t C/ha. In the Ruzyneˇ

region (Czech Republic) [

] arable land has 32 t C/ha in the first 20 cm (48 t C/ha at 30 cm),

according also to a calibration study with the same model, while we obtained 29

±

11 t C/ha. In

the Halle region (Germany) [

77

] maize production without fertilizer applications leads to the

accumulation of 59.7 t C/ha and 43.2 t C/ha (at 35 cm depth– 51 and 37 t C/ha at 30 cm), while
in the same region we obtained 49.6 t C/ha and 48.6 t C/ha at 30 cm, respectively.

This assessment shows that, discounting differences in sampling depth, our results are gen-

erally in agreement with measurements. In particular, there is good agreement for European
temperate regions, as the default parameterization of RothC used here was obtained for those
conditions (originally for measurements of SOC in Harpenden, United Kingdom [

]). How-

ever, this is not the case for all regions of the World. Despite the notable exceptions mentioned
earlier, the estimation agreement is lower in regions with significant different climate condi-
tions to those where RothC were calibrated. Estimated attainable SOC for wheat in three
Indian regions (Barrackpore, Ranchi and Akola) [

78

] were 40% lower on average. SOC stocks

for wheat in those regions (at 30 cm depth) of 39, 24 and 19 t C/ha were reported. Respectively
for the same three regions we obtained 17

±

4, 15

±

5 and 14

±

2 t C/ha. In the Sidney region

(Australia), long-term grasslands reached between 58 t C/ha and 82 t C/ha (depending on the

Global modelling of soil organic carbon

PLOS ONE |

September 19, 2019

15 / 27

tillage practice). In our work, we obtained 38

±

3 t C/ha. Long-term rice cultivation in South

Korea leads to a SOC stock of 60t C/ha (without fertilization) and between 60 and 69 t C/ha,
depending on the fertilization rate, while we obtained 30

±

6 t C/ha and 32 t C/ha, respectively.

Regarding the mineralization rate, an average range reported for sorghum production in Bur-
kina Faso [

79

] is 1.5 to 2.6%, which is similar to 3

±

1% obtained in our results. The differences

can be justified due the differences in climatic variables between studies (e.g. the average rain-
fall used in this study is 200 mm higher). The average rate for fertilized grasslands in England
[

80

] are 3.0%, compared with 3.7% in our study. For oak forests in Italy [

81

], mineralization

rates reportedly range between 1 and 6%, while for the same region we obtained an estimated
5

±

1% for temperate oceanic forests. As for the role of management practices on SOC, for

maize production in a Loess Plateau (China) [

82

] leaving straw on the field increases annual

SOC stock in 0.150 t C/ha.year. For the same region we obtained 0.350

±

0.035 t C/ha.year,

considering the difference in SOC between the first year and at the end of the 86-year simula-
tion. A study in central Ohio (USA) [

83

] calculated that wheat residues left on the field

increase SOC by about 3 t C/ha.year. For the same region, we obtained 2

±

2 t C/ha.year on

average for all LU classes.

A systematic comparison of our results with individual studies would require a large meta-

analysis of SOC measurement studies, which is outside the scope of this study. Our results are
valid over relatively aggregated regions that are approximately homogeneous (i.e. UHTUs), but
the studies in the literature are valid for specific points within those regions. It would be neces-
sary to aggregate the results of those studies into broader areas in order to compare them with
our regionally applicable results. At the moment, these comparisons provide promising indica-
tions that our results depict the main trends when comparing regions and LU classes. This pro-
vides some qualitative assurance for the main insights obtained from the results in this paper. In
most cases where the literature reports measured SOC stocks that are significantly different
from our results, the relative comparisons between LU classes and locations are very rarely dif-
ferent. If a given LU class in a given region has a higher measured SOC stock than another LUC
class in another region, our model results typically also report a higher SOC stock.

Besides individual studies, we also compared our results with pre-established SOC data-

bases. This enables a systematic assessment of the agreement between results. We used the
geo-referenced 19,000-measurement LUCAS database for the EU [

17

]. A geospatial analysis of

LUCAS depicts the same general conclusions of our work. Considering only EU countries, the
regions with higher attainable SOC stocks are Northern Europe, Southeast France and the
British Isles, which are the regions where higher SOC field measurements are located in
LUCAS. This fact can be explained by the combination of three factors: temperatures are rela-
tively low temperate, precipitation is relatively moderate (in the case of British Isles) or low
(Southeast France) and organic soils, which are typically high in SOM [

], are present.

Results of geospatial correlation analysis shows that for 13 LU classes out of 21 evaluated SOC
stocks in our work and LUCAS measurements are significantly correlated at 5% (mean Pear-
son’s r is 0.212). For sugar beet and rapeseed no spatial correlation was found. This result is
due to the fact that our results depict attainable SOC stocks, while LUCAS measurements
depict SOC stocks at an unknown time after transformation, with no assurance that SOC at
each location has stabilized. The fact that a positive correlation was found is reassuring, as a
stronger correlation could not be expected.

Comparison of assumptions with previous uses of the model

Regarding the choice of parameters made in previously published works using RothC at point/
plot scales, we observed a widespread use of default DPM/RPM ratios and decomposition

Global modelling of soil organic carbon

PLOS ONE |

September 19, 2019

16 / 27

rates for the different organic matter pools (

Table 1

). Most model applications (197) use only

defaults (i.e. DPM/RPM ratio of 1.44 for cropland and improved grasslands, 0.67 for unim-
proved grassland and scrub and 0.25 for forest; and decomposition rate of 10, 0.30, 0.66, 0.02
1/year for DPM, RPM, HUM and BIO, respectively). The decomposition rate of the HUM
pool was the parameter most often re-parameterized, followed by the DPM/RPM ratio of the
soil C inputs. The HUM pool is most often re-parameterized (as for example [

86

–

88

]) as it has

the most sensitive kinetic constant (due to the longer residence time). Regarding land classes
(data shown in

S2 File

), all three main LU types (cropland, grassland and forest) are covered

by the sampled papers, and for each of them multiple sub-classes are considered—e.g. in the
cropland class, the studies assessed involve 20 different crops/rotations (among others, maize,
wheat and rice). Forest is the class that required more re-parametrizations (23 out of 52 appli-
cations), mainly in the HUM decomposition rate. Nevertheless, twenty-two of those re-param-
etrizations are from only one study (Shirato et al. [

86

]). Croplands required the least changes

(only 5), which was expected considering that RothC was originally developed for croplands.
This low number of re-parameterizations in studies that applied the model at very fine scales
provides some assurance that the default configuration of the model used in this work works
well for most regions of the World.

The resilience of our results is also demonstrated when comparing our results to results

obtained using different parametrizations for specific regions. For example, RothC was origi-
nally developed for humid regions and is known to overestimate C mineralization in dry
regions [

89

]. A small change in the soil moisture rate modifier was proposed (where the mini-

mum value of the rate modifier was changed to 0.1 –originally it was 0.2) [

89

]. Comparing the

attainable SOC stocks in a dry region (Alentejo, Portugal) that were obtained in a previous
regional study [

] using the modified version of RothC to the attainable SOC stocks obtained

in this paper for the same region (using the “default” RothC parameter set), there are no statis-
tically significant differences in results (considering the confidence intervals for both applica-
tions of the model). For example, the average attainable SOC stocks for irrigated maize with
residues left on the field obtained in this paper was about 38

±

7 t C/ha, while in the regionally

modified application it was 41

±

4 t C/ha. The average obtained in each model application for

each land class always falls within the confidence interval obtained in the other. Therefore,
despite of the reported overestimation of C mineralization under dry conditions, there is no
statistically significant difference in results for our application at UHTU scale.

Comparison of results with large-scale RothC applications

We only found 2 continental/global applications [

https://doi.org/10.1371/journal.pone.0222604.t001

,

36

] of RothC in the first 200 references

(out of the 521). As in this paper, all of them used the model without re-parametrization.
Those studies include the main three LU classes, cropland, grassland and forest, without dis-
criminating between sub-classes (e.g. specific crops). For grasslands, the model application
performed by Smith et al. [

36

] at European scale considered a DPM/RPM ratio of 1.44 (i.e, the

Table 1. Accounting of the RothC model parametrization. BIO microbial biomass; DPM—easily decomposable plant material; HUM—humified organic matter;
RPM—resistant plant material.

Was

it changed?

DPM/RPM ratio

Decomposition rate

DPM

RPM

HUM

BIO

Yes

10

(4%)

2

(1%)

4

(2%)

36

(15%)

1

(

<

1%)

No

222

(96%)

231

(99%)

229

(98%)

197

(85%)

232

(

>

99%)

Global modelling of soil organic carbon

PLOS ONE |

September 19, 2019

17 / 27

default for improved grassland), while the global application performed by Gottschalk et al.
[

] considered 0.67 (i.e. the default for unimproved grasslands). Our study is therefore

aligned with past practices in terms of parameter selection but adds significant regional and
LU detail when compared to those prior applications. Both papers focus on the effects of cli-
mate change and estimate a SOC loss in cropland and grassland areas in Europe due to warmer
temperatures. They also note that SOC loss may be slowed by decreased soil moisture. Those
results cannot be directly compared to the results obtained in this paper because our study was
not forward-looking and did not take the effects of future climate change (temperature and
precipitation) into account. However, our results have implications for studies aimed to esti-
mating the effect of climate change on SOC loss. We show that SOC loss/gain is highly related
with the specific crop type. In the same UHTU, land occupation with one crop type could lead
to SOC gains and with another crop to SOC loss due to different yields and their response to
changing climate conditions. These results therefore suggest that future similar studies should
consider at minimum the same level of differentiation between crop types and regions used in
this paper.

The use of a Monte Carlo approach in this work also had additional value as it avoided

assuming that the number chosen for each input variable was representative of each entire
region. Input data varied per iteration according to a pre-set probability distribution depicting
the regional variability of the data within each UHTU. The average of the 100 iterations per-
formed for each UHTU and land class, which is indicated in this work as the most likely result
(with accompanying uncertainty intervals), is therefore representative of the UHTUs they
depict.

Limitations, extensions and future work

Comparison of RothC with other models. We selected RothC for this paper, but RothC

is one of several soil process-based models available that could have been used instead. Among
those, two notable models that were considered as alternatives for the research conducted here
were the widely used CENTURY [

90

] and DNDC [

91

] models. Originally, all models simulate

the first 20 cm depth, but all were already used at 30 cm depth. All models run using a monthly
time step. CENTURY has a broader scope than RothC, including not only soil C, but also N,
phosphorous and sulphur dynamics. DNDC includes the C and N cycle in the soil. Having a
larger scope for modelled processes, a limiting factor for using CENTURY and DNDC in full
is the higher number of variables/parameters required. Some variables are difficult to obtain
for a global application with the characteristics of the one in this paper. We focused specifically
on processes influencing C dynamics, and, as RothC includes only the C cycle, the number of
parameters involved and data required is lower. The lower number of required input data
makes a spatially differentiated and crop-specific global assessment feasible. Higher computa-
tional time requirement to run CENTURY and DNDC are additional limiting factors to use
these models at global scale.

For C, RothC considers a total of five pools (two litter pools, DPM and RPM, and three soil

pools, HUM, BIO and IOM), while DNDC has four pools (plant residue, microbial biomass,
active humus and passive humus) and CENTURY includes three SOM pools (active, slow and
passive) plus two litter pools (structural and metabolic) that are not included in the soil. The
simplicity of RothC, however, comes at a cost as it does not model some processes. It is the
only of the three models that omits leaching of organic matter and that does not include a
management practices module (other models can simulate the effects of fertilization, harvest-
ing, fire, grazing, irrigation and erosion–in some cases through extensions or more recent
modules built on top of the original model). The lack of depictions of these processes are

Global modelling of soil organic carbon

PLOS ONE |

September 19, 2019

18 / 27

limitations that may be overcome through extensions to RothC and are suggested for further
revisions of work similar to what we present in this paper.

Model parameterization. Despite several limitations inherent to the task of performing

global SOC modelling, such as the simplified depiction of soil biogeochemical processes by the
model and uncertainty of the input data (e.g. plant residues), our results qualitatively depict
well-established regional and global trends, such as the fact that SOC stocks are lower for agri-
cultural LU classes and that SOC loss in transitions to croplands is a faster process than SOC
recovery. Our results are also robust to differences within each LU class, e.g. attainable SOC
stock for wheat production is higher in northern Europe than southern Asia [

76

78

]. Here, we

used a “generic” parametrization of the RothC model, which together with other factors deter-
mines the obtained attainable SOC stock and mineralization rates. It is likely that general esti-
mation accuracy could be improved using local parametrization of the model, as performed by
some prior RothC applications [

86

88

], for those regions and classes where the estimated

results deviate from field-measured values. We expect locally specific parameterizations to
improve results because the agreement of the model with other published studies seems to be
better for temperate regions closer to the ones used for the original calibration of the model.
The quantitative comparison of results with data presented here is also promising despite the
fact that more work is needed to assess the regions and classes for which the model works
properly and the ones where the model fails to depict realistic SOC stock dynamics.

Attainable SOC stock and mineralization rates are also influenced by other methodological

choices done in this work. The Monte Carlo determination of intra-UHTU uncertainty did
not take into consideration the uncertainty in RothC parameters. Only input variables had an
assigned probability distribution within each UHTU, but even for those we assumed as a sim-
plification that all input variables were normally distributed. When we obtained the data from
a geospatial dataset (i.e. we estimated the parameters for a normal distribution starting from
pixel-level data within each UHTU), we verified statistically that the assumption was plausible.
For some data, however, we were unable to do so as we had to assume a distribution based on
one likely value plus a confidence interval (e.g. IPCC-based data on yields). It is unlikely that
all input variables are in fact distributed normally. For example, clay content is frequently used
in modelling according to a normal distribution (e.g. [

]), but empirically that is rarely the

case [

92

]. This information influences the depiction of intra-regional variability through the

uncertainty quantified using Monte Carlo analysis.

Calculation of soil carbon inputs. Soil C input is a critical factor in explaining differences

between attainable SOC, as demonstrated by complementary analysis reported in

Yield data (FAOSTAT [

44

]) and the factors used to calculate above and belowground residues

(IPCC [

42

]) contribute the most to the soil C inputs and their uncertainty. In this work, for

all agricultural LU classes, we considered constant plant yields, as well as constant above-
ground and belowground residue C inputs and C from manure application.

This method introduces error in cases where belowground plant growth does not scale with

aboveground yield. There is some evidence that belowground C varies less between crops than
aboveground C. For example, Cagnarini et al. [

35

] conclude that the use of constant factors

overestimates belowground biomass when compared with field-measured belowground bio-
mass, for maize and wheat in Switzerland. Taghizadeh-Toosi et al. [

93

] used long-term experi-

mental data to conclude that belowground biomass growth was independent of the crop. In
general, the literature overall shows that the type of crop influences belowground biomass (e.g.
[

,

]), but belowground C is also highly influenced by the year and farming system (e.g.

conventional or organic) [

95

]. The method used in this paper prevented us from simulating

intra-annual changes in crop residue C inputs from belowground as we use constant near-past
average monthly temperature and precipitation. Studies that, unlike this one, are interested in

Global modelling of soil organic carbon

PLOS ONE |

September 19, 2019

19 / 27

understanding the dynamics of SOC gains and losses due to climate change or other dynamic
processes should take into account that belowground biomass can be nearly uncorrelated with
the response of the aboveground biomass, which is typically easier to estimate.

Additionally, in order to simplify the calculation of C inputs, we also considered the same C

content of the crop, forest and grassland residues (i.e. 0.40 kg C/kg DM) and constant C/N
ratio of manure applied (i.e. 14 kg C/kg N). While there is good evidence for the C content of
living plant biomass [

], there are many types of manure with varying chemical compositions

[

96

]. Due to lack of data, we were unable to consider specific types of manure applied in each

UHTU, or their quantities. As future work, a sensitivity analysis should be performed for those
critical parameters.

Depiction of agricultural management practices. A remaining significant limitation is

that the model was unable to assess a wider variety of management practices. For example, this
work does not yet include an assessment of SOC changes in no-till systems or crop rotations,
as well as the influence of forest fires [

97

]. For woody crops, an evaluation of the differences

between bare soil and use of cover crops is also missing. All these effects are crucial and deserve
inclusion in future models or iterations of this work. Another limitation of this work is that
only one grassland class was included, despite the multitude of grazing systems and respective
management practices in the World. Follow-up work should involve multiple grassland clas-
ses, as well as the explicit modelling of one or more shrubland classes. Having multiple grass-
land types (and also grassland management) and expanding the scope of cropland
management practices should also enable more comprehensive comparisons of attainable
SOC stocks within and between LU classes. Within classes, it will produce results regarding
the most SOC-increasing practices in each UHTU. Between classes, it will produce updates
results regarding where in the world and for which combinations of practices particular crop-
land can accumulate more SOC than grassland or forest.

At the moment the IPCC [

] only provides SOC stock estimates (under native vegetation)

for a combination of 9 climate regions and 6 soil types. To calculate SOC stock for other classes
(e.g. long term cultivated land), the IPCC method applies simple management factors to the
estimated SOC stock under native vegetation depending on climate region, moisture regime
and tillage type. However, this procedure is not crop-specific and limited in biogeographical
differentiation. We thus consider that the use of the results provided in this study is a signifi-
cant step forward that can provide an alternative, more accurate quick assessment tool. These
results can be used for large-scale assessments of SOC change in the absence of local informa-
tion, such as Tier 2 IPCC defaults [

] involved in the calculation of C sequestration and C

emissions due to LUC.

Potential extensions of the work. The data in this paper also enables the historical assess-

ment of SOC gains and losses due to past land use change, as well as prospective studies for
scenario assessment of future land conversions. For example, studies that assess the effects of
land use in the recent past on SOC can use the data presented here to estimate the sign and
magnitude of the difference before and after changes. Average mineralization rates per LU
class and region can also be used for additional modelling exercises.

A different application of the work presented here is in Life Cycle Assessment (LCA) stud-

ies. There has been a push in recent years to produce accurate and meaningful LCA impact
assessment models determining the impacts of land use [

98

102

]. For the role of land use on

the loss of the soil’s biotic production potential, SOC depletion is commonly used as a proxy
indicator. After the original models based on IPCC data [

103

104

], more advanced models

using statistical interpretation of global SOC maps have been proposed [

84

105

]. These models

are very limited by data availability and end up depicting a very limited number of land classes
or miss critical biogeographical differences [

105

]. Process-based modelling of SOC can

Global modelling of soil organic carbon

PLOS ONE |

September 19, 2019

20 / 27

overcome those limitations and provide accurate and highly detailed LCA characterization fac-
tors. A demonstration of the potential of using RothC for this end was published for the region
of Alentejo, Portugal [

]. The data presented in this paper enables the extension of that work

for the World.

Conclusion

This work provides a new global SOC assessment enabled by process-based modelling rather
than conventional statistical upscaling of local measurements. Global modelling is not new,
but here we introduce important innovations in SOC modelling regarding the depth (number
of land classes) and reach (spatial representation) of the analysis. We calculated annual and
attainable SOC stocks for 80 combinations of LU and management practices in approximately
17,000 regions in the World. This study therefore significantly increases simultaneously the
number of LU classes and the spatial resolution of modelling, compared to other global estima-
tions of potential global SOC stocks.

We conclude that using aggregated “cropland”, “agriculture”, “arable land” or other similar

land use or land cover classes in global assessment is too coarse to produce meaningful results,
as demonstrated by the extreme variability of results for each cropland sub-class within the
same region and between regions. We showed that this approach to modelling can produced
novel results absent from past modelling exercises. For example, we also showed that, by per-
forming this differentiated, crop-specific analysis, some crops in particular regions of the
World can accumulate more SOC than forests and grasslands. This effect exists mostly when
crop residues are left on the field, which has a double effect: the amount of C entering the soil
is higher, and the soil remains covered for a longer period of time, reducing aeration and min-
eralization. These are particular cases that do not invalidate the generic observation that con-
versions to cropland overall reduce SOC, but they should be acknowledged as those exceptions
are relatively frequent and exist even for highly important and representative crops such as
maize.

Supporting information

S1 File. Supplementary methods, results and future applications.
(DOCX)

S2 File. Comparison between model results and long-term experiment data.
(XLSX)

S3 File. Unique Homogeneous Territorial Unit (UHTU) map.
(DOCX)

S4 File. Curve fitting results.
(DOCX)

Acknowledgments

This work was supported by FCT/MCTES (PIDDAC) through project UID/EEA/50009/2019
and by project “Animal Future—Steering Animal Production Systems towards Sustainable
Future", funded by the Fundac¸

ã

o para a Ci

ê

ncia e Tecnologia (SusAn/0001/2016). T.G. Morais

was supported by grant SFRH/BD/115407/2016 and R.F.M. Teixeira by grant SFRH/BPD/
111730/2015 from Fundac¸

ã

o para a Ci

ê

ncia e Tecnologia.

Global modelling of soil organic carbon

PLOS ONE |

https://doi.org/10.1111/gcbb.12044

September 19, 2019

21 / 27

Author Contributions

Conceptualization: Tiago G. Morais, Ricardo F.M. Teixeira, Tiago Domingos.

Data curation: Tiago G. Morais.

Formal analysis: Tiago G. Morais.

Investigation: Tiago G. Morais, Ricardo F.M. Teixeira, Tiago Domingos.

Methodology: Tiago G. Morais, Ricardo F.M. Teixeira, Tiago Domingos.

Resources: Ricardo F.M. Teixeira.

Supervision: Ricardo F.M. Teixeira, Tiago Domingos.

Validation: Tiago G. Morais, Ricardo F.M. Teixeira.

Visualization: Tiago G. Morais.

Writing – original draft: Tiago G. Morais, Ricardo F.M. Teixeira, Tiago Domingos.

Writing – review & editing: Ricardo F.M. Teixeira, Tiago Domingos.

References

1.

Buchholz T, Friedland AJ, Hornig CE, Keeton WS, Zanchi G, Nunery J. Mineral soil carbon fluxes in
forests and implications for carbon balance assessments [Internet]. GCB Bioenergy. 2014. pp. 305–
311.

2.

Lehmann J, Kleber M. The contentious nature of soil organic matter. Nature. Nature Publishing Group;
2015; 528: 60.

https://doi.org/10.1038/nature16069

PMID:

26595271

3.

Arnell NW, Lowe JA, Brown S, Gosling SN, Gottschalk P, Hinkel J, et al. A global assessment of the
effects of climate policy on the impacts of climate change. Nat Clim Chang. Nature Research; 2013; 3:
512–519.

https://doi.org/10.1038/nclimate1793

4.

Solomon S, Qin D, Manning M, Chen Z, Marquis M, Averyt KB, et al. Climate Change 2007 Synthesis
Report. Contribution of Working Groups I, II and III to the Fourth Assessment Report of the Intergov-
ernmental Panel on Climate Change. Climate change 2007: Synthesis Report. Contribution of Working
Group I, II and III to the Fourth Assessment Report of the Intergovernmental Panel on Climate Change.
Summary for Policymakers. 2007.

https://doi.org/10.1002/ep.670210305

5.

Poeplau C, Don A, Vesterdal L, Leifeld J, Van Wesemael B, Schumacher J, et al. Temporal dynamics
of soil organic carbon after land-use change in the temperate zone—carbon response functions as a
model approach [Internet]. Global Change Biology. Blackwell Publishing Ltd; 2011. pp. 2415–2427.

https://doi.org/10.1111/j.1365-2486.2011.02408.x

6.

Popp A, Humpeno¨der F, Weindl I, Bodirsky BL, Bonsch M, Lotze-Campen H, et al. Land-use protec-
tion for climate change mitigation. Nat Clim Chang. Nature Publishing Group; 2014; 4: 1095–1098.

https://doi.org/10.1038/nclimate2444

7.

Guo LB, Gifford RM. Soil carbon stocks and land use change: A meta analysis. Glob Chang Biol.
Blackwell Science Ltd; 2002; 8: 345–360.

https://doi.org/10.1046/j.1354-1013.2002.00486.x

8.

Lal R. Soil carbon sequestration impacts on global climate change and food security. Science (80-).
2004; 304: 1623–1627 ST-Soil carbon sequestration impacts.

https://doi.org/10.1126/science.

1097396

PMID:

15192216

9.

Doetterl S, Stevens A, Six J, Merckx R, Van Oost K, Casanova Pinto M, et al. Soil carbon storage con-
trolled by interactions between geochemistry and climate. Nat Geosci. Nature Research; 2015; 8:
780–783.

https://doi.org/10.1038/ngeo2516

10.

Reichstein M, Bahn M, Ciais P, Frank D, Mahecha MD, Seneviratne SI, et al. Climate extremes and
the carbon cycle. Nature. Nature Publishing Group, a division of Macmillan Publishers Limited. All
Rights Reserved.; 2013; 500: 287–295. Available:

https://doi.org/10.1038/nature12350

PMID:

23955228

11.

Koven CD, Hugelius G, Lawrence DM, Wieder WR. Higher climatological temperature sensitivity of
soil carbon in cold than warm climates. Nat Clim Chang. Nature Publishing Group; 2017; 7: 817–822.

https://doi.org/10.1038/nclimate3421

Global modelling of soil organic carbon

PLOS ONE |

http://dx.doi.org/10.1038/nclimate2322

September 19, 2019

22 / 27

12.

Giardina CP, Litton CM, Crow SE, Asner GP. Warming-related increases in soil CO2 efflux are
explained by increased below-ground carbon flux. Nat Clim Chang. Nature Publishing Group; 2014; 4:
822–827. Available:

13.

Harris NL, Brown S, Hagen SC, Saatchi SS, Petrova S, Salas W, et al. Baseline Map of Carbon Emis-
sions from Deforestation in Tropical Regions. Science (80-). 2012; 336: 1573–1576.

PMID:

Schrumpf M, Schulze ED, Kaiser K, Schumacher J. How accurately can soil organic carbon stocks
and stock changes be quantified by soil inventories? Biogeosciences. 2011; 8: 1193–1212.

org/10.5194/bg-8-1193-2011

15.

Smith P. How long before a change in soil organic carbon can be detected? Glob Chang Biol. Blackwell
Science Ltd; 2004; 10: 1878–1883.

https://doi.org/10.1111/j.1365-2486.2004.00854.x

16.

Va´clavı´k T, Lautenbach S, Kuemmerle T, Seppelt R. Mapping global land system archetypes. Glob
Environ Chang. Pergamon; 2013; 23: 1637–1647.

https://doi.org/10.1016/J.GLOENVCHA.2013.09.

004

17.

To´th G, Jones A, Montanarella L. The LUCAS topsoil database and derived information on the regional
variability of cropland topsoil properties in the European Union. Environ Monit Assess. 2013; 185:
7409–25.

https://doi.org/10.1007/s10661-013-3109-3

PMID:

23371251

18.

de Brogniez D, Ballabio C, Stevens A, Jones RJA, Montanarella L, van Wesemael B. A map of the top-
soil organic carbon content of Europe generated by a generalized additive model. Eur J Soil Sci. Black-
well Publishing Ltd; 2015; 66: 121–134.

https://doi.org/10.1111/ejss.12193

19.

IPCC. 2006 IPCC Guidelines for National Greenhouse Gas Inventories. Institute for Global Environ-
mental Strategies (IGES) for the Intergovernmental Panel on Climate Change. Kanagawa, Japan:
The Intergovernmental Panel on Climate Change (IPCC); 2006.

20.

Del Grosso SJ, Gollany HT, Reyes-Fox M, Del Grosso S, Ahuja L, Parton W. Simulating Soil Organic
Carbon Stock Changes in Agroecosystems using CQESTR, DayCent, and IPCC Tier 1 Methods. Syn-
thesis and Modeling of Greenhouse Gas Emissions and Carbon Storage in Agricultural and Forest
Systems to Guide Mitigation and Adaptation. American Society of Agronomy, Inc., Crop Science Soci-
ety of America, Inc., and Soil Science Society of America, Inc.; 2016. pp. 89–110.

2134/advagricsystmodel6.2013.0001.5

21.

Borgen SK, Dalsgaard L, Arnoldussen A. CO2 emissions from Norwegian cropland: influence of IPCC
tier level. Soil carbon sequestration Clim food Secur Ecosyst Serv. 2013; 27: 101.

22.

Berhongaray G, Alvarez R. The IPCC Tool for predicting soil organic carbon changes evaluated for the
Pampas, Argentina. Agric Ecosyst Environ. 2013; 181: 241–245.

https://doi.org/10.1016/j.agee.2013.

10.002

23.

Batjes NH. Soil organic carbon stocks under native vegetation—Revised estimates for use with the
simple assessment option of the Carbon Benefits Project system. Agric Ecosyst Environ. 2011; 142:
365–373.

https://doi.org/10.1016/j.agee.2011.06.007

24.

ECCC. National Inventory Report 1990–2015: Greenhouse Gas Sources and Sinks in Canada. Envi-
ronment and Climate Change Canada (ECCC); 2017.

25.

CITEPA. Rapport National d’Inventaire pour la France au titre de la Convention cadre des Nations
Unies sur les Changements Climatiques et du Protocole de Kyoto. Centre Interprofessionnel Tech-
nique d’Etudes de la Pollution Atmosphe

´ rique; 2017.

26.

APA. Portuguese National Inventory Report on Greenhouse Gases, 1990–2018. Amadora, Portugal:
Portuguese Environmental Agency; 2018.

27.

Monforti F, Lugato E, Motola V, Bodis K, Scarlat N, Dallemand J-F. Optimal energy use of agricultural
crop residues preserving soil organic carbon stocks in Europe. Renew Sustain Energy Rev. Perga-
mon; 2015; 44: 519–529.

https://doi.org/10.1016/j.rser.2014.12.033

28.

Luo Y, Keenan TF, Smith M. Predictability of the terrestrial carbon cycle. Glob Chang Biol. 2015; 21:
1737–1751.

https://doi.org/10.1111/gcb.12766

PMID:

25327167

29.

Campbell EE, Paustian K. Current developments in soil organic matter modeling and the expansion of
model applications: a review. Environ Res Lett. IOP Publishing; 2015; 10: 123004.

1088/1748-9326/10/12/123004

30.

Coleman K, Jenkinson DS, Crocker GJ, Grace PR, Klı´r J, Ko¨rschens M, et al. Simulating trends in soil
organic carbon in long-term experiments using RothC-26.3. Geoderma. 1997; 81: 29–44.

org/10.1016/S0016-7061(97)00079-7

31.

Liu DL, Chan KY, Conyers MK, Li G, Poile GJ. Simulation of soil organic carbon dynamics under differ-
ent pasture managements using the RothC carbon model. Geoderma. 2011; 165: 69–77.

org/10.1016/j.geoderma.2011.07.005

Global modelling of soil organic carbon

PLOS ONE |

September 19, 2019

23 / 27

32.

Gottschalk P, Smith JU, Wattenbach M, Bellarby J, Stehfest E, Arnell N, et al. How will organic carbon
stocks in mineral soils evolve under future climate? Global projections using RothC for a range of cli-
mate change scenarios. Biogeosciences. Copernicus GmbH; 2012; 9: 3151–3171.

https://doi.org/10.1016/j.jclepro.2018.05.061

5194/bg-9-3151-2012

33.

Morais TG, Silva C, Jebari A, A

´ lvaro-Fuentes J, Domingos T, Teixeira RFM. A proposal for using pro-

cess-based soil models for land use Life cycle impact assessment: Application to Alentejo, Portugal. J
Clean Prod. 2018; 192: 864–876.

34.

Morais TG, Teixeira RFM, Rodrigues NR, Domingos T. Characterizing livestock production in Portu-
guese sown rainfed grasslands: Applying the inverse approach to a process-based model. Sustainabil-
ity. Multidisciplinary Digital Publishing Institute; 2018; 10: 4437.

https://doi.org/10.3390/su10124437

35.

Cagnarini C, Renella G, Mayer J, Hirte J, Schulin R, Costerousse B, et al. Multi-objective calibration of
RothC using measured carbon stocks and auxiliary data of a long-term experiment in Switzerland. Eur
J Soil Sci. John Wiley & Sons, Ltd (10.1111); 2019; ejss.12802.

https://doi.org/10.1111/ejss.12743

36.

Smith J, Smith P, Wattenbach M, Zaehle S, Hiederer R, Jones RJA, et al. Projected changes in mineral
soil carbon of European croplands and grasslands, 1990–2080. Glob Chang Biol. Wiley/Blackwell
(10.1111); 2005; 11: 2141–2152.

https://doi.org/10.1111/j.1365-2486.2005.001075.x

37.

Ingram JSI, Fernandes ECM. Managing carbon sequestration in soils: concepts and terminology.
Agric Ecosyst Environ. Elsevier; 2001; 87: 111–117.

https://doi.org/10.1016/S0167-8809(01)00145-1

38.

Kutsch WL, Bahn M, Heinemeyer A. Soil Carbon Dynamics: An Integrated Methodology [Internet].
Cambridge: Cambridge University Press; 2009. Available:

https://www.amazon.com/Soil-Carbon-

Dynamics-Integrated-Methodology/dp/0521865611

39.

Weihermu¨ller L, Graf A, Herbst M, Vereecken H. Simple pedotransfer functions to initialize reactive
carbon pools of the RothC model. Eur J Soil Sci. Blackwell Publishing Ltd; 2013; 64: 567–575.

https://

doi.org/10.1111/ejss.12036

40.

Chapagain AK, Hoekstra AY. Water footprint of nations. Volume 1: Main report. Value of Water
Research Report Series—Volume 1: Main Report, Value of Water Research Report Series No. 16.
Delft. Delft, The Netherlands; 2004.

41.

FAO/IIASA. Global Agro-ecological Zones (GAEZ v3.0). IIASA, Laxenburg, Austria FAO, Rome, Italy
iv. 2012;

https://doi.org/10.1029/97GB03657

42.

IPCC. Revised 1996 IPCC Guidelines for National Greenhouse Gas Inventories. Volume 2—Work-
book. The Intergovernmental Panel on Climate Change (IPCC), the Organization for Economic Co-
operation and Development (OECD) and the International Energy Agency (IEA); 1997.

43.

IPCC. Good Practice Guidance for Land Use, Land-Use Change and Forestry. Institute for Global
Environmental Strategies (IGES) for the Intergovernmental Panel on Climate Change. Kanagawa:
The Intergovernmental Panel on Climate Change (IPCC); 2003.

44.

FAO. Food and Agriculture Organization of the United Nations—Statistics Division [Internet]. 2018
[cited 25 Mar 2015]. Available:

http://faostat.fao.org/

45.

Jebari A, del Prado A, Pardo G, Rodrı´guez Martı´n JA, A

´ lvaro-Fuentes J. Modeling Regional Effects of

Climate Change on Soil Organic Carbon in Spain. J Environ Qual. The American Society of Agronomy,
Crop Science Society of America, and Soil Science Society of America, Inc.; 2018; 47: 644.

org/10.2134/jeq2017.07.0294

PMID:

30025039

46.

Jebari A. Estimacio

´ n de los cambios en los stocks de carbono del suelo agrı´cola a escala regional:

Impacto de los usos del suelo y del manejo en la Comunidad Auto´noma de Arago

´ n. Master Thesis.

2016.

47.

Robinson TP, William Wint GR, Conchedda G, Van Boeckel TP, Ercoli V, Palamara E, et al. Mapping
the global distribution of livestock. Baylis M, editor. PLoS One. Public Library of Science; 2014; 9:
e96084.

https://doi.org/10.1371/journal.pone.0096084

PMID:

24875496

48.

Mueller ND, Gerber JS, Johnston M, Ray DK, Ramankutty N, Foley JA. Closing yield gaps through
nutrient and water management. Nature. Nature Publishing Group; 2012; 490: 254–257.

https://doi.org/10.1021/es1041755

org/10.1038/nature11420

PMID:

22932270

49.

Roma

´ n P, Martı´nez MM, Pantoja A. Farmer’s Compost Handbook—Experiences in Latin America.

Santiago, Chile: Food and Agriculture Organization of the United Nations (FAO); 2015.

50.

Pfister S, Bayer P, Koehler A, Hellweg S. Environmental impacts of water use in global crop produc-
tion: hotspots and trade-offs with land use. Environ Sci Technol. American Chemical Society; 2011;
45: 5761–8.

PMID:

21644578

51.

Thornthwaite CW. An Approach toward a Rational Classification of Climate. Geogr Rev. 1948; 38: 55.

https://doi.org/10.2307/210739

Global modelling of soil organic carbon

PLOS ONE |

https://lpdaac.usgs.gov/dataset_discovery/modis/modis_products_table/

September 19, 2019

24 / 27

52.

DAAC L. MODIS/Terra Land Surface Temperature and Emissivity Monthly L3 Global 0.05Deg CMG
[Internet]. 2016. Available:

mod11c3

53.

NASA. Global Precipitation Analysis [Internet]. 2016. Available:

http://precip.gsfc.nasa.gov/

54.

Joshi MM, Lambert FH, Webb MJ. An explanation for the difference between twentieth and twenty-first
century land–sea warming ratio in climate models. Clim Dyn. Springer Berlin Heidelberg; 2013; 41:
1853–1869.

https://doi.org/10.1007/s00382-013-1664-5

55.

IPCC. Climate Change 2013—The Physical Science Basis. Climate Change 2013—The Physical Sci-
ence Basis. 2014.

https://doi.org/10.1017/cbo9781107415324

56.

Donat MG, Lowry AL, Alexander L V., O’Gorman PA, Maher N. More extreme precipitation in the
world’s dry and wet regions. Nat Clim Chang. Nature Publishing Group; 2016; 6: 508–513.

http://esdac.jrc.ec.europa.

org/10.1038/nclimate2941

57.

ESDAC. Global Soil Organic Carbon Estimates [Internet]. 2012. Available:

eu/content/global-soil-organic-carbon-estimates

58.

FAO, IIASA, ISRIC, ISSCAS, JRC. Harmonized World Soil Database (version 1.2). FAO, Rome, Italy
and IIASA, Laxenburg, Austria; 2012. 3123

59.

NASA LP DAAC. Land Cover Type Yearly L3 Global 0.05Deg CMG (MCD12C1). In: NASA EOSDIS
Land Processes DAAC, USGS Earth Resources Observation and Science (EROS) Center [Internet].
2017 [cited 20 Feb 2017]. Available:

https://lpdaac.usgs.gov/dataset_discovery/modis/modis_

products_table/mcd12c1

60.

Fischer G, Nachtergaele F, Prieler S, van Velthuizen H, Verelst L, Wiberg D. Global Agro-ecological
Zones Assessment for Agriculture [Internet]. 2008. Available:

http://www.iiasa.ac.at/Research/LUC/

luc07/External-World-soil-database/HTML/index.html?sb=1

61.

Thematic Mapping. World Borders Dataset [Internet]. 2009 [cited 20 Feb 2017]. Available:

http://

thematicmapping.org/downloads/world_borders.php

62.

Rothamsted Research. Rothamsted Carbon Model (RothC) [Internet]. 2017.

63.

Metropolis N, Ulam S. The Monte Carlo method. J Am Stat Assoc. 1949; 44: 335–341.

10.1080/01621459.1949.10483310

PMID:

18139350

64.

Clarivate Analytics. InCites Journal Citation Report [Internet]. 2017 [cited 29 Dec 2017]. Available:

https://jcr.incites.thomsonreuters.com/

65.

ESDAC. LUCAS 2009 TOPSOIL data [Internet]. 2017 [cited 8 Jan 2017]. Available:

http://eusoils.jrc.

ec.europa.eu/content/lucas-2009-topsoil-data

66.

Nilsson M-C, Wardle DA. Understory vegetation as a forest ecosystem driver: evidence from the north-
ern Swedish boreal forest. Front Ecol Environ. Ecological Society of America; 2005; 3: 421–428.

https://doi.org/10.1890/1540-9295(2005)003[0421:UVAAFE]2.0.CO;2

67.

Mokany K, Raison RJ, Prokushkin AS. Critical analysis of root: Shoot ratios in terrestrial biomes. Glob
Chang Biol. Wiley/Blackwell (10.1111); 2006; 12: 84–96.

https://doi.org/10.1111/j.1365-2486.2005.

001043.x

68.

Frank AB, Liebig MA, Tanaka DL. Management effects on soil CO2 efflux in northern semiarid grass-
land and cropland. Soil Tillage Res. Elsevier; 2006; 89: 78–85.

https://doi.org/10.1016/J.STILL.2005.

06.009

69.

Lee X, Huang Y, Huang D, Hu L, Feng Z, Cheng J, et al. Variation of Soil Organic Carbon and Its Major
Constraints in East Central Asia. Hui D, editor. PLoS One. Public Library of Science; 2016; 11:
e0150709.

https://doi.org/10.1371/journal.pone.0150709

PMID:

26934707

70.

Wieder WR, Bonan GB, Allison SD. Global soil carbon projections are improved by modelling microbial
processes. Nat Clim Chang. Nature Publishing Group; 2013; 3: 909–912.

https://doi.org/10.1038/

nclimate1951

71.

Bationo A, Kihara J, Vanlauwe B, Waswa B, Kimetu J. Soil organic carbon dynamics, functions and
management in West African agro-ecosystems. Agric Syst. Elsevier; 2007; 94: 13–25.

10.1016/j.agsy.2005.08.011

72.

Wang G, Luo Z, Han P, Chen H, Xu J. Critical carbon input to maintain current soil organic carbon
stocks in global wheat systems. Sci Rep. Nature Publishing Group; 2016; 6: 19327.

PMID:

Luo Z, Wang E, Sun OJ. Can no-tillage stimulate carbon sequestration in agricultural soils? A meta-
analysis of paired experiments. Agric Ecosyst Environ. Elsevier; 2010; 139: 224–231.

10.1016/J.AGEE.2010.08.006

Global modelling of soil organic carbon

PLOS ONE |

September 19, 2019

25 / 27

74.

Cai ZC, Qin SW. Dynamics of crop yields and soil organic carbon in a long-term fertilization experiment
in the Huang-Huai-Hai Plain of China. Geoderma. Elsevier; 2006; 136: 708–715.

1016/J.GEODERMA.2006.05.008

75.

Zhang W, Liu K, Wang J, Shao X, Xu M, Li J, et al. Relative contribution of maize and external manure
amendment to soil carbon sequestration in a long-term intensive maize cropping system. Sci Rep.
Nature Publishing Group; 2015; 5: 10791.

https://doi.org/10.1038/srep10791

PMID:

26039186

76.

Lemke RL, VandenBygaart AJ, Campbell CA, Lafond GP, Grant B. Crop residue removal and fertilizer
N: Effects on soil organic carbon in a long-term crop rotation experiment on a Udic Boroll. Agric Eco-
syst Environ. Elsevier; 2010; 135: 42–51.

https://doi.org/10.1016/J.AGEE.2009.08.010

77.

Ludwig B, John B, Ellerbrock R, Kaiser M, Flessa H. Stabilization of carbon from maize in a sandy soil
in a long-term experiment. Eur J Soil Sci. John Wiley & Sons, Ltd; 2003; 54: 117–126.

10.1046/j.1365-2389.2003.00496.x

78.

Manna MC, Swarup A, Wanjari RH, Ravankar HN, Mishra B, Saha MN, et al. Long-term effect of fertil-
izer and manure application on soil organic carbon storage, soil quality and yield sustainability under
sub-humid and semi-arid tropical India. F Crop Res. Elsevier; 2005; 93: 264–280.

1016/j.fcr.2004.10.006

79.

Bationo A, Buerkert A. Soil organic carbon management for sustainable land use in Sudano-Sahelian
West Africa. Nutr Cycl Agroecosystems. Kluwer Academic Publishers; 2001; 61: 131–142.

org/10.1023/A:1013355822946

80.

Jarvis SC, Stockdale EA, Shepherd MA, Powlson DS. Nitrogen Mineralization in Temperate Agricul-
tural Soils: Processes and Measurement. Adv Agron. Academic Press; 1996; 57: 187–235.

org/10.1016/S0065-2113(08)60925-6

81.

Rey A, Petsikos C, Jarvis PG, Grace J. Effect of temperature and moisture on rates of carbon minerali-
zation in a Mediterranean oak forest soil under controlled and field conditions. Eur J Soil Sci. Blackwell
Science Ltd; 2005; 56: 589–599.

https://doi.org/10.1111/j.1365-2389.2004.00699.x

82.

Ding S, Xue S, Liu G. Effects of long-term fertilization on oxidizable organic carbon fractions on the
Loess Plateau, China. J Arid Land. Science Press; 2016; 8: 579–590.

https://doi.org/10.1007/s40333-

016-0007-x

83.

Duiker SW, Lal R. Crop residue and tillage effects on carbon sequestration in a Luvisol in central Ohio.
Soil Tillage Res. Elsevier; 1999; 52: 73–81.

https://doi.org/10.1016/S0167-1987(99)00059-8

84.

Morais TG, Domingos T, Teixeira RFM. A spatially explicit life cycle assessment midpoint indicator for
soil quality in the European Union using soil organic carbon. Int J Life Cycle Assess. Springer Berlin
Heidelberg; 2016; 21: 1076–1091.

https://doi.org/10.1007/s11367-016-1077-x

85.

Jones C, McConnell C, Coleman K, Cox P, Falloon P, Jenkinson D, et al. Global climate change and
soil carbon stocks; predictions from two contrasting models for the turnover of organic carbon in soil.
Glob Chang Biol. Blackwell Science Ltd; 2005; 11: 154–166.

https://doi.org/10.1111/j.1365-2486.

2004.00885.x

86.

Shirato Y, Hakamata T, Taniyama I. Modified rothamsted carbon model for andosols and its validation:
changing humus decomposition rate constant with pyrophosphate-extractable Al. Soil Sci Plant Nutr.
Taylor & Francis Group; 2004; 50: 149–158.

https://doi.org/10.1080/00380768.2004.10408463

87.

Luo Z, Wang E, Zheng H, Baldock JA, Sun OJ, Shao Q. Convergent modelling of past soil organic car-
bon stocks but divergent projections. Biogeosciences. 2015; 12: 4373–4383.

https://doi.org/10.5194/

bg-12-4373-2015

88.

A

´ lvaro-Fuentes J, Morell FJ, Plaza-Bonilla D, Arru´e JL, Cantero-Martı´nez C. Modelling tillage and

nitrogen fertilization effects on soil organic carbon dynamics. Soil Tillage Res. Elsevier; 2012; 120: 32–
39.

https://doi.org/10.1016/J.STILL.2012.01.009

89.

Farina R, Coleman K, Whitmore AP. Modification of the RothC model for simulations of soil organic C
dynamics in dryland regions. Geoderma. Elsevier; 2013;200–201: 18–30.

https://doi.org/10.1016/J.

GEODERMA.2013.01.021

90.

Parton WJ, Schimel DS, Cole C V., Ojima DS. Analysis of Factors Controlling Soil Organic Matter Lev-
els in Great Plains Grasslands. Soil Sci Soc Am J. Soil Science Society of America; 1987; 51: 1173–
1179.

https://doi.org/10.2136/SSSAJ1987.03615995005100050015X

91.

Li C, Frolking S, Frolking TA. A model of nitrous oxide evolution from soil driven by rainfall events: 1.
Model structure and sensitivity. J Geophys Res Atmos. Wiley-Blackwell; 1992; 97: 9759–9776.

https://

doi.org/10.1029/92JD00509

92.

Brejda JJ, Moorman TB, Smith JL, Karlen DL, Allan DL, Dao TH. Distribution and Variability of Surface
Soil Properties at a Regional Scale. Soil Sci Soc Am J. 2010;

https://doi.org/10.2136/sssaj2000.

643974x

Global modelling of soil organic carbon

PLOS ONE |

https://doi.org/10.1038/srep32568

September 19, 2019

26 / 27

93.

Taghizadeh-Toosi A, Christensen BT, Glendining M, Olesen JE. Consolidating soil carbon turnover
models by improved estimates of belowground carbon input. Sci Rep. Nature Publishing Group; 2016;
6: 32568.

PMID:

27580672

94.

Smith J, Smith P, Wattenbach M, Gottschalk P, Romanenkov VA, Shevtsova LK, et al. Projected
changes in the organic carbon stocks of cropland mineral soils of European Russia and the Ukraine,
1990–2070. Glob Chang Biol. John Wiley & Sons, Ltd; 2007; 13: 342–356.

https://doi.org/10.1111/j.

1365-2486.2006.01297.x

95.

Hu T, S

ø

rensen P, Wahlstro

¨ m EM, Chirinda N, Sharif B, Li X, et al. Root biomass in cereals, catch

crops and weeds can be reliably estimated without considering aboveground biomass. Agric Ecosyst
Environ. Elsevier; 2018; 251: 141–148.

https://doi.org/10.1016/J.AGEE.2017.09.024

96.

Aguilera E, Guzma

´ n GI, A

´ lvaro-Fuentes J, Infante-Amate J, Garcı´a-Ruiz R, Carranza-Gallego G, et al.

A historical perspective on soil organic carbon in Mediterranean cropland (Spain, 1900–2008). Sci
Total Environ. 2018; 621: 634–648.

https://doi.org/10.1016/j.scitotenv.2017.11.243

PMID:

29202285

97.

Powlson DS, Stirling CM, Jat ML, Gerard BG, Palm CA, Sanchez PA, et al. Limited potential of no-till
agriculture for climate change mitigation. Nat Clim Chang. Nature Publishing Group; 2014; 4: 678–
683.

https://doi.org/10.1038/nclimate2292

98.

Souza DM, Teixeira RFM, Ostermann OP. Assessing biodiversity loss due to land use with Life Cycle
Assessment: are we there yet? Glob Chang Biol. 2015; 21: 32–47.

https://doi.org/10.1111/gcb.12709

PMID:

25143302

99.

Teixeira RFM, Maia de Souza D, Curran MP, Anto´n A, Michelsen O, Mil

à

i Canals L. Towards consen-

sus on land use impacts on biodiversity in LCA: UNEP/SETAC Life Cycle Initiative preliminary recom-
mendations based on expert contributions. J Clean Prod. 2016; 112: 4283–4287.

1016/j.jclepro.2015.07.118

100.

Curran M, de Souza DM, Anto´n A, Teixeira RFM, Michelsen O, Vidal-Legaz B, et al. How Well Does
LCA Model Land Use Impacts on Biodiversity?-A Comparison with Approaches from Ecology and
Conservation. Environ Sci Technol. American Chemical Society; 2016; 50: 2782–95.

https://doi.org/10.1016/j.jclepro.2016.05.077

10.1021/acs.est.5b04681

PMID:

26830787

101.

Vidal Legaz B, Maia De Souza D, Teixeira RFM, Anto´n A, Putman B, Sala S. Soil quality, properties,
and functions in life cycle assessment: an evaluation of models. J Clean Prod. 2017; 140: 502–515.

102.

Teixeira RFM, Morais TG, Domingos T. A Practical Comparison of Regionalized Land Use and Biodi-
versity Life Cycle Impact Assessment Models Using Livestock Production as a Case Study. Sustain-
ability. Multidisciplinary Digital Publishing Institute; 2018; 10: 4089.

https://doi.org/10.3390/

SU10114089

103.

Mil

à

i Canals L, Mu

ñ

oz I, McLaren S, Brand

ã

o M. LCA Methodology and Modelling Considerations for

Vegetable Production and Consumption [Internet]. CES Working Papers 02/07; 2007. Available:

http://

www.ces-surrey.org.uk/

104.

Brand

ã

o M, Mil

à

i Canals L. Global characterisation factors to assess land use impacts on biotic pro-

duction. Int J Life Cycle Assess. 2013; 18: 1243–1252.

https://doi.org/10.1007/s11367-012-0381-3

105.

Teixeira RFM, Morais TG, Domingos T. Consolidating Regionalized Global Characterization Factors
for Soil Organic Carbon Depletion Due to Land Occupation and Transformation. Environ Sci Technol.
American Chemical Society; 2018; 52: 12436–12444.

https://doi.org/10.1021/acs.est.8b00721

PMID:

30253100

Global modelling of soil organic carbon

PLOS ONE |