Abramoff et al. (2022). Improved global-scale predictions of soil carbon stocks with Millennial Version 2

rose.abramoff@lsce.ipsl.fr

Improved global-scale predictions of soil carbon stocks with Millennial

Version 2

Rose Z. Abramoff

Raphael A. Viscarra Rossel

Laboratoire des Sciences du Climat et de l

’

Environnement, LSCE/IPSL, CEA-CNRS-UVSQ, Universit

´

e Paris-Saclay, F-91191, Gif-sur-Yvette, France

b

Laboratoire de G

´

eologie, UMR 8538, Ecole Normale Sup

´

erieure, IPSL, PSL Research University, CNRS, Paris, France

c

Universit

´

e Libre des Bruxelles, 1050, Bruxelles, Belgium

d

Department of Earth System Science, Stanford University, Yang

&

Yamazaki Building, 473 Via Ortega, Stanford, CA, 94305-4216, USA

e

Physical and Life Sciences Directorate, Lawrence Livermore National Laboratory, Livermore, CA, 94551, USA

f

Biology Department, San Diego State University, San Diego, CA, 92130, USA

g

Soil and Landscape Science, School of Molecular and Life Sciences, Curtin University, GPO Box U1987, Perth WA, 6845, Australia

h

School of Atmospheric Sciences, Sun Yat-Sen University, Zhuhai, 510245, Guangdong, China

A R T I C L E I N F O

Keywords:

Soil carbon modeling

Microbial decomposition

Mineral association

Soil organic carbon stocks

A B S T R A C T

Soil carbon (C) models are used to predict C sequestration responses to climate and land use change. Yet, the soil

models embedded in Earth system models typically do not represent processes that reflect our current under-

standing of soil C cycling, such as microbial decomposition, mineral association, and aggregation. Rather, they

rely on conceptual pools with turnover times that are fit to bulk C stocks and/or fluxes. As measurements of soil

fractions become increasingly available, it is necessary for soil C models to represent these measurable quantities

so that model processes can be evaluated more accurately. Here we present Version 2 (V2) of the Millennial

model, a soil model developed to simulate C pools that can be measured by extraction or fractionation, including

particulate organic C, mineral-associated organic C, aggregate C, microbial biomass, and low molecular weight C.

Model processes have been updated to reflect the current understanding of mineral-association, temperature

sensitivity and reaction kinetics, and different model structures were tested within an open-source framework.

We evaluated the ability of Millennial V2 to simulate total soil organic C (SOC), as well as the mineral-associated

and particulate fractions, using three independent data sets of soil fractionation measurements spanning a range

of climate and geochemistry in Australia (N

=

495), Europe (N

=

175), and across the globe (N

=

659). When

using all the data together (N

=

1329), the Millennial V2 model predicted SOC (RMSE

=

3.3 kg C m

−

2

, AIC

=

675,

R

2

in

=

0.31,

R

2

out

=

0.26) better than the widely-used first-order decomposition model Century (RMSE

=

3.4

kg C m

−

2

, AIC

=

696,

R

2

in

=

0.21,

R

2

out

=

0.18) across sites, despite the fact that Millennial V2 has an increase in

process complexity and number of parameters compared to Century. Millennial V2 also reproduced the observed

fraction of C in MAOM and larger particle size fractions for most latitudes and biomes, and allows for a more

detailed understanding of the pools and processes that affect model performance. It is important to note that this

study evaluates the spatial variation in C stock only, and that the temporal dynamics of Millennial V2 remain to

be tested. The Millennial V2 model updates the conceptual Century model pools and processes and represents our

current understanding of the roles that microbial activity, mineral association and aggregation play in soil C

sequestration.

* Corresponding author.

E-mail addresses:

(R.Z. Abramoff),

betrand.guenet@ens.fr

(B. Guenet),

Haicheng.Zhang@ulb.ac.be

(H. Zhang),

georgiou@stanford.edu

r.viscarra-rossel@curtin.edu.au

georgiou1@llnl.gov

(K. Georgiou),

xxu@sdsu.edu

(X. Xu),

(R.A. Viscarra Rossel),

yuanwp3@mail.sysu.edu.cn

(W. Yuan),

philippe.

ciais@lsce.ipsl.fr

(P. Ciais).

Contents lists available at

ScienceDirect

Soil Biology and Biochemistry

journal homepage:

www.elsevier.com/locate/soilbio

https://doi.org/10.1016/j.soilbio.2021.108466

Received 21 May 2021; Received in revised form 14 October 2021; Accepted 18 October 2021

Soil Biology and Biochemistry 164 (2022) 108466

2

1. Introduction

Soils are not only a vast storage pool of carbon (C) but also a

potentially important feedback to climate change, as they store water,

mineral nutrients and organic matter, and exchange materials with local

waterways and the atmosphere. Exchanges of greenhouse gases between

soils and the atmosphere are particularly relevant to C cycle-climate

feedbacks, and are controlled jointly by biological activity including

plant input and decomposition of organic matter by soil microorgan-

isms, and by physical processes controlling soil water and temperature

(

Lajtha et al., 2018

). Soil decomposition models are used to assess cur-

rent soil organic matter (SOM) stocks and make predictions under future

climate conditions. Most models project soil C loss in response to future

warming, implying a positive feedback to climate change (

Sulman et al.,

2018

). Yet, large uncertainties arise from soil carbon-climate feedbacks

due to structural and parametric uncertainties of soil C models (

). These uncertainties can be partially constrained with data, but

only if measurements can be directly related to modeled quantities.

When models are constrained or evaluated using data, it is often

necessary to make assumptions about what modeled quantities repre-

sent and often the measurements themselves are proxies for underlying

biogeochemical processes (

Bailey et al., 2018

). However, the most

commonly-used SOM models define pools based on turnover time that

are not easily or consistently related to physical measurements (

Parton

et al., 1987

) or they define physically-based pools that use older defi-

nitions (

Jenkinson and Rayner, 1977

). Our best understanding of bio-

logical, physical, and chemical processes in soil has advanced in the

decades since these models were popularized. In the past decade alone,

many models have incorporated microbial decomposition (

) and stabilization of SOM

via aggregation (

Segoli et al., 2013

;

Stamati et al., 2013

) or mineral

association (

;

;

), reflecting field and laboratory studies emphasizing the impor-

tance of these processes for SOM stocks and persistence (

Tisdall and

Oades, 1982

;

Torn et al., 1997

;

Kallenbach et al., 2016

There is a general consensus about which processes are important for

controlling soil organic C (SOC) cycling (

), yet there is less consensus about the best mathematical

formulation of different processes (

). The kinetics of chemical reactions can be approximated in

different ways depending on assumptions. For example,

Michaelis-Menten kinetics assumes that the concentration of the sub-

strate greatly exceeds the concentration of the enzyme (forward

Michaelis-Menten) or vice versa (reverse Michaelis-Menten) (

). Equilibrium Chemistry Approxima-

tion (ECA) makes no assumption about the relative concentrations of

substrate versus enzyme (

Tang and Riley, 2013

Tang, 2015

), and may

therefore be the most generally applicable approximation of reaction

kinetics. However, a recent study comparing multiple approximation

methods suggested that because depolymerization is effectively limited

by enzyme binding sites, reverse Michaelis-Menten may be more

appropriate than ECA for modeling decomposition (

Tang and Riley,

). Similar to kinetics approximations, equations that relate soil

temperature and moisture to soil processes vary widely in their forms,

causing divergent model predictions (

Rodrigo et al., 1997

). The earliest

models often used empirically-derived functions for temperature and

moisture effects, but more recent models make use of relationships that

approximate thermodynamic or diffusive principles, respectively

(

Arrhenius, 1889

Davidson et al., 2012

Van Veen and Kuikman, 1990

Ghezzehei et al., 2018

). Lastly,

model representation of bond formation between organic matter and

mineral surfaces varies widely across models (

Sulman et al., 2018

),

reflecting the still unknown contribution of dissolved organic C sorption

(

Abramoff et al., 2021

) versus other forms of complexation that immo-

bilize OM (

Masiello et al., 2004

),

including aggregation (

The original framework of the Millennial model (

Abramoff et al.,

2018

) sought to update early models such as Century and Roth-C with

two goals: 1) to define C pools which would be related more directly to

field measurements, and 2) to reflect current understanding of soil mi-

crobial and physicochemical processes. To that end, the Millennial

model included explicit representation of microbial activity, association

with minerals via sorption, and aggregation of organic matter. The

original version (

Abramoff et al., 2018

) was compared with the SOM

model currently used in many Earth System Models (e.g., E3SM, CESM,

ORCHIDEE), Century (

Parton et al., 1987

), but was not tested against an

independent dataset of measurements. Therefore, it remained an open

question whether the Millennial model could indeed predict SOC stocks

better than a first-order decomposition model

–

an important gap that

we address herein. In this study, we update the equations of the

Millennial model, test alternate model structures, and evaluate the

ability of the model to predict spatial variation in SOC stocks and un-

derlying soil fractions across multiple biomes based on first principles.

2. Methods

2.1. Model development from Version 1

The original Millennial model (hence, Version 1 or V1) equations

were developed to facilitate comparisons with the Century model

(

Abramoff et al., 2018

). As such, many of the equations of this original

version follow a similar structure to, or are borrowed directly from, the

Century model. For example, the temperature and moisture scalars used

in Millennial V1 are taken from the daily-time step version of the Cen-

tury model (

Del Grosso et al., 2005

Parton et al., 2010

). In this paper,

we maintain the conceptual model as presented in

, but update the governing equations to reflect recent de-

velopments in representations of temperature and moisture dynamics,

microbial mortality, and protection of organic matter by association

with minerals. The full equations of this new version, Millennial Version

2 (V2), are presented in the next section, but the main differences be-

tween Millennial V1 and V2 are summarized in

Table 1

. Specifically, the

temperature function was updated from the Century temperature scalar

used in

to the Arrhenius equation (

Davidson

et al., 2012

Abramoff et al., 2017

). The moisture function was updated

from the Century moisture scalar used in

Parton et al. (2010)

to a

relationship representing the effects of matric potential, oxygen limita-

tion, and diffusion on reaction rates (

Ghezzehei et al., 2018

). The

maximum capacity of minerals to sorb organic matter was estimated in

Millennial V1 using a relationship derived from 72 incubations of dis-

solved organic C with different soil types (

Mayes et al., 2012

). In

Millennial V2 we use a broader empirical relationship based on the C

content of

>

1200 measurements of mineral-associated organic matter

(

Georgiou et al., 2021

). We updated the microbial mortality equation

from one based on a fixed rate constant to one that includes

density-dependence of microbial biomass, sensu

Georgiou et al. (2017)

.

We also tested several alternate approximations of reaction kinetics,

described in detail in Section

2.3

. Millennial V1 has 23 fittable param-

eters, and Millennial V2 has 24 fittable parameters.

2.2. Model description of Millennial V2

The system of equations below follows the conceptual figure (

Fig. 1

),

tracking the size of and transfers between five C pools: particulate

organic matter (POM; denoted

P

in equations), low molecular weight

carbon (LMWC,

L

), aggregate C (AGG,

A

), mineral-associated organic

matter (MAOM,

M

), and microbial biomass (MIC,

B

). For all Millennial

V2 equations below, descriptions, units, and default values for variables

can be found in

. The change in POM (

P

) stock with time is

governed by the balance between plant litter C input and aggregate C

breakdown, aggregate C formation, and decomposition,

R.Z. Abramoff et al.

Soil Biology and Biochemistry 164 (2022) 108466

3

dP

/

dt

=

p

i

F

i

+

p

a

F

a

−

F

pa

−

F

pl

[1]

where

p

i

is the proportion of C input allocated to POM;

F

i

is the flux of

aboveground plant litter, root litter and root exudates;

p

a

is the pro-

portion of C in aggregate breakdown allocated to POM;

F

a

is aggregate C

breakdown;

F

pa

is the aggregate carbon formation from POM; and

F

pl

is

depolymerization of POM into LMWC. Depolymerization of POM is

governed by the reverse Michaelis-Menten equation (

Tang and Riley,

),

F

pl

=

V

pl

S

w

,

D

P

B

K

pl

+

B

[2]

V

pl

=

α

pl

e

−

Ea

pl

/

[

R

(

T

+

273

.

15

)]

[3]

S

w

,

D

=

(

θ

φ

)

0

.

5

[4]

where

V

pl

is the maximum rate of POM decomposition modified by an

Arrhenius temperature relationship (

Davidson et al., 2012

;

Sierra,

2012

), and is a function of a pre-exponential constant

α

pl

, an activation

energy

Ea

pl

, the universal gas constant

R

, and the soil temperature in

Celsius (

T

).

K

pl

is the half-saturation constant, and

B

is the microbial

biomass carbon. The moisture modifier

S

w,D

refers to the diffusion lim-

itation of substrates defined in

Ghezzehei et al. (2018)

as the square root

of the ratio of the volumetric water content (

θ

) and the total porosity (

φ

).

The formation (

F

pa

) of aggregate C (

A

) from POM is a function of the

rate of aggregate formation (

k

pa

) and

S

W

,

D

,

F

pa

=

k

pa

S

w

,

D

P

[5]

Similarly, soil aggregate C breakdown (

F

a

) into POM and MAOM is a

function of the rate of breakdown

(

k

b

) and

S

W

,

D

.

F

a

=

k

b

S

w

,

D

A

[6]

The change in LMWC (

L

) depends on LMWC input, the leaching rate

(

F

l

), decomposition of POM, adsorption to minerals (

F

lm

), and microbial

uptake (

F

lb

), and the proportion of microbial necromass that enters

LMWC (

p

b

), microbial mortality (

F

bm

), and desorption (

F

ld

). In a multi-

layer version of the Millennial model, LMWC would also depend on

leaching input from other soil layers, but in this single layer version we

assume that the leaching input is included in the LMWC input,

dL

dt

=

F

i

(

1

−

p

i

) −

F

l

+

F

pl

−

F

lm

−

F

lb

+ (

1

−

p

b

)

F

bm

+

F

ld

[7]

where

F

l

is the LMWC leaching loss,

F

l

=

k

l

S

w

,

D

L

[8]

and where

k

l

is the leaching rate. Adsorption of LMWC to MAOM (

M

)

follows a Langmuir-type saturating relationship

sensu

Wang et al.

,

F

lm

=

S

w

,

D

K

lm

L

(

1

−

M

Q

max

)

[9]

where

K

lm

is the binding affinity that is adjustable based on the pH.

Q

max

is the maximum sorption capacity (gC/m

2

), assuming a 1 m soil profile.

K

lm

is further defined by the parameters

p

1

,

p

2

, and the desorption co-

efficient (

K

ld

),

K

lm

=

e

−

p

1

pH

−

p

2

K

ld

[10]

where parameters

p

1

and

p

2

are the coefficients for computing the first

term

(

e

−

p

1

pH

−

p

2

in

L

/

mg

)

from the site-level pH. This term is derived

from

Mayes et al. (2012)

and is equivalent to the equilibrium constant.

By multiplying the equilibrium constant by the desorption coefficient,

we obtain the coefficient of adsorption, following the principle that the

equilibrium constant is equal to the ratio of the adsorption coefficient

and desorption coefficient (i.e.,

K

eq

=

K

ads

/K

des

) (

Wang et al., 2013

). The

Table 1

Main differences between Millennial V1 and Millennial V2. Please refer to

Abramoff et al. (2018)

for a full description of Millennial V1 equations, to Section

2.2

for a

full description of Millennial V2 equations, and to

for definitions of variables, their units, and values.

Process

Millennial V1 [Equation Number in

Abramoff et al., 2018

]

Millennial V2 [Equation Number in this paper (Section

2.2

)]

Temperature function

S

t

=

t

2

+

(

t

3

π

)

atan

[

π

(

T

−

t

1

)]

t

2

+

(

t

3

π

)

atan

[

π

t

4

(

T

ref

−

t

1

)]

[3]

V

x

=

α

x

e

−

Ea

x

/

[

R

(

T

+

273

.

15

)]

[3,14]

Moisture function

S

w

=

1

+

w

1

exp

(−

w

2

RWC

)

[4]

S

w

=

e

λφ

[

k

a

,

min

+(

1

−

k

a

,

min

)

(

φ

−

θ

φ

)

0

.

5

](

θ

φ

)

0

.

5

[4,15]

Sorption to minerals

Q

max

=

BD

10

c

1

log

(

%

clay

) +

c

2

[11]

Q

max

=

depth BD

%

claysilt p

c

[11]

Microbial mortality

F

bm

=

k

mm

S

t

S

w

B

[17]

F

bm

=

k

bd

B

2

[16]

Fig. 1.

Conceptual model of Millennial V2, following Millennial V1 of

Abramoff et al. (2018)

. Black boxes are soil C pools. Solid arrows indicate fluxes between pools

(See

for definitions). Colored boxes indicate modeled processes that have been updated or changed from Version 1. The dash line indicates that microbial

biomass controls the depolymerization rate. POM

=

particulate organic matter, LMWC

=

low molecular weight carbon, MAOM

=

mineral-associated organic matter.

R.Z. Abramoff et al.

Soil Biology and Biochemistry 164 (2022) 108466

4

maximum sorption capacity

Q

max

in g C kg

−

1

soil can be estimated using

the empirical equation from

Georgiou et al. (2021)

, which depends on

the clay and silt content in percent (

%claysilt

) and a coefficient (

p

c

). This

expression is then converted to model units of g C m

−

2

using the

site-level sampling depth (

depth

) in m and bulk density (

BD

) in kg soil

m

−

3

.

Q

max

=

depth BD

%

claysilt p

c

[11]

Desorption of LMWC is a function of MAOM,

Q

max

, and the desorp-

tion coefficient (

K

ld

)

sensu

Wang et al. (2013)

F

ld

=

K

ld

M

Q

max

[12]

Microbial uptake of LMWC is a function of microbial biomass and

LMWC concentration, temperature, and water,

F

lb

=

V

lb

S

w

,

B

L

K

lb

+

L

[13]

V

lb

=

α

lb

e

−

Ea

lb

/

[

R

(

T

+

273

.

15

)]

[14]

S

w

,

B

=

e

λφ

[

k

a

,

min

+

(

1

−

k

a

,

min

)

(

φ

−

θ

φ

)

0

.

5

]

S

w

,

D

[15]

where

V

lb

is the maximum uptake rate of LMWC, modified by an

Arrhenius temperature relationship as in Equation

(3)

, and

K

lb

is the

half-saturation constant for microbial activity. The term

S

w,B

refers to

the total moisture sensitivity of biological activity (

Ghezzehei et al.,

2018

), where

λ

is the dependence of the rate on the matric potential (

φ

)

and

k

a

,

min

is the minimum relative rate in saturated soil. This relationship

incorporates the effects of matric potential (

e

λφ

), oxygen limitation

[

k

a

,

min

+(

1

−

k

a

,

min

)

(

φ

−

θ

φ

)

0

.

5

]

and diffusion limitation of substrates

(

S

W

,

D

).

Microbial biomass mortality is calculated using a density-dependent

formulation derived from

Georgiou et al. (2017)

, where the mortality

(

F

bm

) is a function of a fixed microbial death rate (

k

bd

), and the square of

the microbial biomass pool

B

, which can be derived from the logistic

growth equation of population dynamics.

F

bm

=

k

bd

B

2

[16]

Both MAOM and POM can enter the aggregate C pool (

A

),

dA

/

dt

=

F

ma

+

F

pa

−

F

a

[17]

F

ma

=

k

ma

S

w

,

D

M

[18]

where

F

ma

is the carbon flow from MAOM to aggregate C, and

k

ma

is the

aggregate formation rate from MAOM. MAOM is formed by sorption of

LMWC and microbial necromass, and is affected by transfer into and out

of the aggregate C pool,

dM

/

dt

=

F

lm

−

F

ld

+

p

b

F

bm

−

F

ma

+

F

a

(

1

−

p

a

)

[19]

where

F

bm

is the carbon flow from microbial biomass to MAOM, namely

buildup of microbial necromass. The partitioning of microbial necro-

mass to MAOM versus LMWC (

p

b

) is not controlled by Langmuir

adsorption but rather assumes that MAOM is made up of both mineral

surfaces which respond to sorption (Equations

9 and 12

) as well as mi-

crobial necromass that forms mineral associations by other means, ac-

counting for the observed discrepancy between maximum sorption

capacities measured using DOC sorption experiments (

Abramoff et al.,

2021

) and the observed maximum capacity of the mineral fraction as a

whole (

Georgiou et al., 2021

).

Microbial biomass changes as a result of uptake, mortality, and loss

via respiration,

dB

/

dt

=

F

lb

−

F

bm

−

F

mr

[20]

Uptake is partitioned into respiration (

F

mr

) and growth (

F

bg

) based on

a temperature-dependent carbon use efficiency (CUE),

F

mr

=

F

lb

{

1

−

[

CUE

ref

−

CUE

T

(

T

−

T

ae

−

ref

)]}

[21]

F

bg

=

F

lb

[

CUE

ref

−

CUE

T

(

T

−

T

ae

−

ref

)]

[22]

where

CUE

ref

is the reference CUE, and

CUE

T

is the CUE dependence on

temperature.

T

ae-ref

and

T

are the reference and current soil temperature,

respectively. Therefore, the total CO

2

released through heterotrophic

respiration is

dCO

2

/

dt

=

F

mr

[23]

2.3. Alternate approximations of reaction kinetics

We tested three methods for approximating the reaction kinetics

governing depolymerization and microbial uptake. In the Millennial V1

model, we used the double Michaelis-Menten equation to approximate

depolymerization and forward Michaelis-Menten for uptake. However,

recent models use a variety of kinetics approximations ranging from

complex to simple. Here we test three of them: (1) a combination of

reverse and forward Michaelis-Menten (MM), (2) equilibrium chemistry

approximation (ECA), and (3) linear (LIN) kinetics.

Following arguments from

Tang and Riley (2019)

on the different

limiting factors for depolymerization and microbial uptake, we defined a

model using reverse Michaelis-Menten kinetics for depolymerization,

and forward Michaelis-Menten kinetics for microbial uptake.

Michaelis-Menten approximations are arguably the most popular reac-

tion kinetics approximations in soil microbial decomposition models

(

Allison et al., 2010

;

Davidson et al., 2012

;

). This model variant is the default for

Millennial V2, and depolymerization and uptake are described by

Equations

2 and 13

, respectively, from the previous section.

Second, we tested the equilibrium chemistry approximation, which

is a more complex approximation than Michaelis-Menten because it is

closer to the full reaction kinetics (

Tang, 2015

). The ECA is commonly

used in models with microbial and/or mineral interactions (

Tang and

Riley, 2015a

Abramoff et al., 2017

). For this model variant, we replaced

Equations

2 and 13

with 2a and 13a, respectively.

F

pl

=

V

pl

S

w

,

D

P

B

K

pl

+

B

+

P

[2a]

F

lb

=

V

lb

S

w

,

B

L

K

lb

+

L

+

B

[13a]

Lastly, linear kinetics is the simplest approximation, wherein the

reaction is linear with respect to each contributing pool. However, note

that because the reaction rate depends on two pools, the interaction is

second-order. For this model variant, we replaced Equations

2 and 13

with 2b and 13b, respectively.

F

pl

=

V

pl

S

w

,

D

P

B

K

pl

[2b]

F

lb

=

V

lb

S

w

,

B

L

K

lb

[13b]

2.4. Parametric sensitivity, collinearity, and sensitivity to inputs

For each model variant (V1, ECA, MM, LIN), we calculated a local

parameter sensitivity index (

S

ij

), which summarizes the effect of each

parameter

j

at timestep

i

for each model pool,

R.Z. Abramoff et al.

Soil Biology and Biochemistry 164 (2022) 108466

5

S

ij

=

dy

i

dθ

j

θ

j

y

i

[24]

Where

y

is the model output and

θ

is the parameter value (

Soetaert and

Herman, 2009

). We initialized the model using the steady-state C pools

at default parameter values, calculated using the

stode

function of R

Version 4.0.4

’

s rootSolve Version 1.8.2.1 package (

Soetaert, 2009

), as in

Wieder et al. (2015b)

. Then, we ran each instance of the model for 100

years using the same repeated global average year of forcing from

Abramoff et al. (2018)

, and evaluated the model output

y

at time step

i

=

100 years for each parameter value

θ

.

Collinearity describes the linear dependence of model parameters,

and can be summarized in a collinearity index,

γ

, where a change in one

parameter can be compensated by 1-1/

γ

by changing other parameters.

For example, a model with

γ

=

20 can offset the effect of a parameter

change by 95% by changing other parameters.

γ

has a lower bound of 1

when all terms are orthogonal and an upper bound of infinity when all

terms are linearly dependent. It is calculated using the normalized

sensitivity matrix

S

ij

with the

collin

function of R

’

s FME Version 1.3.6.1

package, and is more fully described in

Marschmann et al. (2019)

. Very high values can be common

with multiple parameters; for example, in a 5-parameter model of bac-

terial growth described in

Soetaert (2016)

, the collinearity index was

2.4

×

10

6

when using all five parameters.

Using the steady state version of the model we conducted an analysis

of model sensitivity to inputs across the range of inputs represented in

the evaluation sites (see 2.6.1): net primary production [NPP;

0.007

–

1.99 gC m

−

2

d

−

1

], soil temperature [-3.47

–

29.9

◦

C], volumetric

water content [0.10

–

0.48 m

3

m

-3

], clay and silt percentage [1

–

98], and

pH [2.8

–

7.9]. Model steady state solutions were evaluated at 10 soil

temperatures, 10 volumetric water contents, and 10 values of NPP

equally spaced along the range in a full factorial design, for a total of

1000 model evaluations. The three V2 model variants were very similar

in their patterns, so we plotted only the model variant with the lowest

collinearity index, V2 MM. We performed the same sensitivity analysis

with the Century model, which is described more fully in the next

section.

2.5. Century and gradient-boosted models

Because the Century model is a commonly-used process-based model

that is included in many terrestrial biosphere models, we used it as a

point of comparison for the Millennial model. The Century model in-

cludes three soil pools (active, slow and passive) and two litter pools

(structural and metabolic). The equations governing interactions be-

tween the five Century model pools are defined in

,

and are also well-described in

Sierra and Muller (2015)

. We did not

include the nitrogen and plant submodels that are described in

Parton

et al. (1987)

, focusing only on carbon cycling in this study. The default

parameters are described in

Table A2

, and the equations reproduced in

Appendix B

.

Both the Millennial and Century models are process-based models

that use assumptions about soil processes to make SOC stock predictions.

In order to evaluate the relationship between measured SOC stocks and

the model forcing inputs (NPP, soil temperature, volumetric water

content, pH, and percentage of clay and silt) without making any

process-level assumptions, we created an empirical null model using a

gradient-boosted machine learning (GBM) algorithm (R package:

caret

V6.0-86

, functions:

createDataPartition, trainControl, train

with method

“

gbm

”

). GBM models were trained on 80% of the data and evaluated on

the final 20% with 10-fold cross validation, repeated 10 times. Perfor-

mance metrics for this model include the out-of-sample R

2

from cross-

validation, the root-mean-square error, and mean absolute error. The

relative influence of the forcing input predictors was calculated using

the decrease in error when the predictor was used to split regression

trees in the model (

Friedman, 2001

Viscarra Rossel and Hicks (2015)

2.6. Model simulations

2.6.1. Data used to fit and evaluate models

We evaluated the Millennial V2 and Century models using three

datasets of soil fractionation measurements. For the purposes of this

evaluation, we used size fractionation measurements separating MAOM

from larger size fractions using a particle size upper bound that ranged

from 50 to 60

μ

m, depending on the dataset. It is important to note that

according to these size fractionation protocols, and some mixed density-

and size-fractionation protocols like that of

Poeplau et al. (2017)

, the

MAOM fraction is defined by size alone, but the aggregate and POM

fractions, though we do not isolate them in this evaluation, can be

defined in different ways. Traditional fractionation protocols disperse

the soil to measure POM and MAOM fractions, but measure micro- and

macroaggregate fractions on soil that is minimally dispersed (

Poeplau

et al., 2018

). As a result, aggregate fractions measured by sieving

minimally-dispersed soil often contain POM and MAOM. For our pur-

poses, we define the aggregate fraction in Millennial V2 to be the stable

microaggregates which remain after dispersion in the larger particle size

fraction (

>

50

–

60

μ

m), and therefore do not contain substantial MAOM.

This is analogous to the

“

heavy sand and stable aggregate fraction

”

in

Poeplau et al. (2017)

, or

“

coarse, heavy POM

”

within the framework of

Robertson et al. (2019)

, who used a combination of size and density

fractionation methods similar to

Poeplau et al. (2017)

. For our purposes,

and to be inclusive of multiple fractionation methods, we focus on

separating the MAOM fraction from larger particle size fractions (POM

and stable microaggregates) using a 10

μ

m range for the upper bound of

MAOM that encompasses most of the fractionation methods that are

commonly used (

Poeplau et al., 2018

). In most fractionation protocols,

microbial biomass and LMWC may be present in both MAOM and larger

size fractions, but because these pools were on average less than 1% of

the total modeled SOM, we grouped them into the larger size fractions,

both for simplicity and to conserve the accuracy of the MAOM fraction

measurement, as this was the only fraction explicitly measured by all of

the datasets.

The first dataset, Viscarra Rossel (VR), is derived from 495 Australian

sites described in

and

Viscarra Rossel

et al. (2019)

. Vegetation ranges from grass to forest with 155 grazed

sites and 345 sites under minimal use or nature conservation. SOC stocks

were reported in t/ha for the top 30 cm of soil and converted to g C m

−

2

.

SOC was fractionated into three fractions, originally reported as par-

ticulate organic C (POC), humic organic C (HOC), and resistant organic

C (ROC), corresponding to the 50

–

2000

μ

m particle size fraction, the

<

50

μ

m particle size fraction, and charcoal measured using solid-state

13

C nuclear magnetic resonance spectroscopy (

Baldock et al., 2013

Beaudoing and Rodell, 2016

For the purposes of this study, we define MAOM as the operational

measurement of soil particles

<

50

μ

m (HOC), while the remaining C

pools are considered to be comprised of soil particles 50

–

2000

μ

m

(POC). The VR dataset reported a number of soil and climate variables,

including bulk density, pH, percent clay and silt, and net primary pro-

ductivity (NPP) derived from the land surface model BIOS2 (

Haverd

et al., 2013

). Soil temperature (

◦

C) and soil moisture (mm

3

mm

−

3

) were

derived from the Global Land Data Assimilation System

’

s (GLDAS) Noah

1

◦

x1

◦

V2.1 land surface model output for the top 0

–

10 cm, 10

–

40 cm,

and 40

–

100 cm of soil (

). We calculated a

weighted average for the three soil layers, depending on the proportion

of soil in each layer which is determined by the depth of sampling at

each site.

The second dataset, Georgiou (KG), is derived from 659 globally-

distributed sites described in

Georgiou et al. (2021)

, including forest,

grass, and crops experiencing varying land management. SOC concen-

trations were reported in g C kg

−

1

soil and converted to g C m

−

2

using

the reported depth of sampling (median: 14 cm; 1st, 3rd quartile: 10 cm,

20 cm; range: 6 cm

–

62 cm), and bulk density for the top 1 m of soil

extracted from SoilGrids 250 m (

Hengl et al., 2017

) aggregated to 1 km.

SOC was fractionated into MAOM, defined as the soil particle size

R.Z. Abramoff et al.

Soil Biology and Biochemistry 164 (2022) 108466

6

fraction

<

60

μ

m. The KG dataset reported percent clay and silt, but not

pH which was derived for the top 1m of soil from SoilGrids in the same

way as bulk density. NPP was derived from a 0.5

◦

x0.5

◦

global data

product of aboveground litter production (

Li et al., 2019

). Because NPP

includes both above and belowground components, we estimated the

belowground (i.e., root) litter production from the aboveground litter

production using the ratio of aboveground-to-belowground biomass

carbon density (Mg/ha) from gridded global maps at 300 m spatial

resolution (

Spawn et al., 2020

). This makes the simplifying assumption

that the ratio of aboveground-to-belowground biomass is analogous to

the ratio of aboveground litter -to-fine root litter. However, the fraction

of belowground biomass in

Spawn et al. (2020)

(median: 0.44, range:

0

–

1) was similar to the fraction of NPP allocated to fine roots (median:

0.47, range: 0.08

–

0.94) in a globally-distributed dataset (N

=

112) (

Xia

et al., 2019

), so we accept this first order assumption for estimating

belowground litter production. The sum of the aboveground litter pro-

duction from

Li et al. (2019)

and the root production estimated from

was used as the NPP input to the

Millennial V2 model for the KG dataset. Soil temperature and moisture

were derived from GLDAS as described for the VR dataset.

The third dataset, LUCAS, is derived from 175 European sites

included in the Land Use/Cover Area Frame Survey (

T

´

oth et al., 2013

),

representing natural and re-vegetated forests and grassland. These sites

measured soil fractions by size fractionation on topsoil with vegetation

residues and litter removed, sampled from 0 to 20 cm (

Cotrufo et al.,

). These measurements were used to evaluate the MEMS model,

another model proposing measurable soil C pools (

Robertson et al.,

). The MAOM pool is defined as the

<

53

μ

m soil particle size

fraction, while the remaining C pools are considered to be comprised of

soil particles 53

–

2000

μ

m. The dataset also measures ancillary soil

properties such as pH and percent clay and silt. Bulk density, NPP, soil

temperature and soil moisture were not reported, and were derived in

the same way as described for the KG dataset. All data used in this study

are summarized in

Table 2

and a map of site locations plotted in

Fig. S1

.

2.6.2. Model fitting and evaluation

When fitting the Millennial V2 model and the Century model across

sites, we assumed that each site had reached a steady-state. As a result,

we could fit the model using the steady-state solution rather than

running the model dynamically. However, whenever the collinearity

index is greater than 20, it is not advisable to fit all of the parameters at

once (

Soetaert, 2016

). Therefore, we chose to optimize parameters with

a local parametric sensitivity index (Equation

(24)

) greater than |0.25|,

or 15/24 fittable parameters in the case of the Millennial V2 model and

13/22 fittable parameters in the case of the Century model (see

Viscarra Rossel and Hicks, 2015

for parameters and their fitted values). Parameters were optimized by

minimizing the sum of squared residuals between the model and ob-

servations using the

modFit

function of R Version 4.0.4

’

s FME Version

1.3.6.1 package with the Levenberg-Marquardt algorithm. Parameter

lower bounds were either -Inf or 0, and upper bounds were either 1 or

Inf. We calculated the coefficient of determination (R

2

) for observed and

predicted pools across sites using the optimal parameter set from fitting.

As a model performance metric, we calculated the root-mean-square

error (RMSE) as

RMSE

=

̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅

∑

n
i

=

1

(

SOC

obs

,

i

−

SOC

mod

,

i

)

2

n

√

[25]

Where

SOC

obs

,

i

is the observed SOC at each site

i

,

SOC

mod

,

i

is the modeled

SOC at each site, and

n

is the number of sites. We also consider the mean

absolute error (MAE) and mean bias error (MBE),

MAE

=

∑

n
i

=

1

⃒

SOC

obs

,

i

−

SOC

mod

,

i

⃒

n

[26]

MBE

=

∑

n
i

=

1

(

SOC

obs

,

i

−

SOC

mod

,

i

)

n

[27]

We calculated the Akaike information criterion (AIC) as a perfor-

mance metric, which considers not only the error but also the number of

model parameters,

AIC

=

n

×

ln

⎛

⎝

̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅

∑

n
i

=

1

(

SOC

obs

,

i

−

SOC

mod

,

i

)

2

n

√

⎞

⎠

+

2

p

[28]

where

p

is the total number of model parameters.

For the purposes of finding the best overall parameter set and

describing relationships between SOC stocks and environmental vari-

ables, we fit the model to the entirety of each dataset. This is also how we

calculate the in-sample R

2

, denoted

R

2

in

Tables 3 and 4

. However, for

the purposes of cross-validation we fit the model to a random draw of

80% of each dataset, holding back 20% for testing, repeated 10 times.

Model performance indices

R

2

out

(out-of-sample R

2

), RMSE, MAE, MBE,

and AIC are calculated by comparing modeled values to the testing

dataset and averaging across the 10 repeated fits.

To classify the range of predicted and observed C stocks by biome, we

Table 2

Sources of data used to fit the model (Particle Size Fractions), force the model (NPP, soil temperature, soil moisture), or define site-specific quantities (Clay, Silt, pH,

Bulk Density).

Dataset

Region

Particle Size Fractions

NPP

Soil temperature

&

moisture

Clay

&

Silt

pH

Bulk Density

VR

Australia

(

Viscarra Rossel et al., 2019

)

BIOS2 (

Haverd et al.,

2013

)

GLDAS (

)

KG

Global

Georgiou et al. (2021)

)

GLDAS (

)

SoilGrids (

SoilGrids (

)

LUCAS

EU

(

;

)

)

GLDAS (

)

SoilGrids (

)

Table 3

Model fit of SOC to measured values.

R

2

in

=

coefficient of determination of the

training dataset (i.e., in-sample),

R

2

out

=

coefficient of determination of the test

dataset (i.e., out-of-sample), RMSE

=

root mean square error, MAE

=

mean

absolute error, MBE

=

mean bias error. AIC

=

Akaike Information Criterion. N is

the number of SOC measurements included in each dataset.

Dataset

Model

R

2

in

R

2

out

RMSE

(kg C

m

−

2

)

AIC

MAE

(kg C

m

−

2

)

MBE

(kg C

m

−

2

)

All (N

=

1329)

Millennial

V2

0.31

0.26

3.29

675

2.14

−

1.01

Century

0.21

0.18

3.42

696

2.30

−

0.72

VR (N

=

495)

Millennial

V2

0.46

0.37

1.83

163

1.25

−

0.63

Century

0.40

0.32

2.11

189

1.46

−

1.17

KG (N

=

659)

Millennial

V2

0.32

0.21

3.46

372

2.28

−

1.18

Century

0.20

0.02

3.81

398

2.63

−

0.35

LUCAS

(N

=

175)

Millennial

V2

0.04

5.03

150

3.71

−

1.10

Century

0.04

0.01

5.21

152

3.77

−

1.19

R.Z. Abramoff et al.

Soil Biology and Biochemistry 164 (2022) 108466

7

used the World Wildlife Fund Terrestrial Ecoregions Map which clas-

sifies vegetated land into 14 biomes (

Olson et al., 2001

). For comparison

purposes, we made use of ancillary measurements collected by each of

the datasets, specifically land type classification from the LUCAS data-

set, grouped into broadleaved, mixed forest, coniferous, mixed grass,

pure grass and re-vegetated. Lastly, to evaluate the range of

model-predicted turnover time of soil C assuming that C stocks are at

equilibrium, we used a global dataset of N

=

470 estimates of turnover

time (

Chen et al., 2013

). Turnover time is defined here as SOC/respi-

ration rate, assuming that steady-state has been reached (

Sierra et al.,

2017

).

3. Results

We estimated the mean parameter sensitivity of C stocks in the five

model pools at equilibrium for Millennial V1 and for the V2 model

variants. Millennial V1 parameter sensitivity is lower than that of the V2

variants for all pools except for the microbial biomass pool and LMWC

(

Fig. S2

). The Millennial V2 variants have similar sensitivity to param-

eters, with the ECA kinetics variant having slightly greater sensitivity of

C stocks to parameter values across all pools except the microbial

biomass pool. The mean collinearity of parameters was highest for the

Millennial V2 variant with linear kinetics, followed by Millennial V1.

Note that Millennial V1 has one fewer fittable parameter than the V2

variants which may lower the collinearity index for V1 (

Fig. 2

). How-

ever, Millennial V2 with Michealis-Menten kinetics had the lowest

collinearity of all the model variants for all pools, indicating that its

parameters are more identifiable compared to the other model variants.

Using this index, we can select V2 MM as the model which is comparable

to the others in terms of parameter sensitivity, but has the most iden-

tifiable parameters, lending itself to easier parameter fitting.

Steady-state solutions of the model are sensitive not only to the

choice of parameters, as demonstrated above, but also the model inputs,

including soil temperature, volumetric water content (VWC), and net

primary production (NPP). Sensitivity analyses conducted across a range

of model inputs show that increasing soil temperature causes all model

pools to decrease in size, as a result of an increase in the reaction rate of

depolymerization and uptake, according to the Arrhenius equation

(

a). Except for the microbial biomass and LMWC pools (

Fig. S3

), all

model pools are smaller at high VWC, which tends to increase the rate of

C cycling and turnover in general (

b). Conversely, all model pools

have a positive relationship with NPP, with all pools increasing in

response to increased plant inputs. The relationship with NPP is not

linear, however, but rather begins to saturate at the upper end of the

observed NPP range because microbial biomass growth increases

decomposition that slows the C accumulation (

c,

Fig. S3

Increasing the percentage of clay and silt present in soil corresponds to

an increase in SOC by decreasing the C that is available to microbes via

mineral association until soil minerals become saturated. Conversely,

increasing pH has the opposite effect, affecting a parameter that inhibits

sorption to minerals which causes SOC to become available to microbes.

In comparison to the Century model, effects of model inputs on steady

state SOC are qualitatively similar for soil temperature and VWC

(

Fig. S4

). The Century model relationship between SOC and NPP is

linear, reflecting the first-order kinetic relationship between inputs and

decomposition rates. Further, the relationship between SOC and the

percentage of clay and silt increases without saturating at values of clay

and silt approaching 100% (

Fig. S4

).

When using all the data together (N

=

1329), the Millennial V2

model predicted SOC (RMSE

=

3.3 kg C m

−

2

, AIC

=

675,

R

2

in

=

0.31,

R

2

out

=

0.26) better than the Century model (RMSE

=

3.4 kg C m

−

2

, AIC

=

696,

R

2

in

=

0.21,

R

2

out

=

0.18;

) across sites. The Millennial V2

model also outperformed the Century model when fit to each dataset

separately (

Fig. 4

). Of the three datasets used, both the

Millennial V2 and Century models were able to explain the greatest

amount of site-level variation in SOC, MAOM and other fractions for the

VR dataset, and the least for the LUCAS dataset (

,

Table 4

). We

made very weak assumptions about the prior range of the parameters

(-Inf/0, 1/Inf), yet most parameters fit to Millennial V2 were within the

same order of magnitude as the default parameters and similar across

datasets (

). The half-saturation constant for microbial uptake

had a large absolute and relative change from the default parameters

after fitting. Other large absolute changes included the activation energy

of uptake and decomposition, as well as the half-saturation constant for

decomposition. On average, the LUCAS dataset

’

s fitted parameters were

less similar to the default parameters than those of the VR and KG

datasets, suggesting that SOC stocks in the LUCAS dataset respond

Table 4

Millennial V2 model fit of MAOM and non-MAOM (POM

+

AGG

+

LMWC

+

MIC) pools to measured values.

R

2

in

=

coefficient of determination of the training

dataset (i.e., in-sample),

R

2

out

=

coefficient of determination of the test dataset (i.

e., out-of-sample), RMSE

=

root mean square error, MAE

=

mean absolute error,

MBE

=

mean bias error. NS

=

not significant at the P

<

0.05 level. N is the

number of SOC measurements included in each dataset.

Dataset

Pool

N

R

2

in

R

2

out

MAE (kg C

m

−

2

)

MBE (kg C

m

−

2

)

All

MAOM

1329

0.40

0.27

1.58

−

1.14

non-

MAOM

670

0.17

0.24

0.82

−

0.077

VR

MAOM

495

0.42

0.37

0.89

−

0.39

non-

MAOM

495

0.29

0.16

0.52

−

0.24

KG

MAOM

659

0.40

0.31

1.38

−

0.62

LUCAS

MAOM

175

0.15

0.09

1.99

−

0.77

non-

MAOM

175

NS

2.20

−

0.33

Fig. 2.

Mean collinearity of parameters for four

model variants, Millennial V1 (Original; 23 parame-

ters), V2 with Linear kinetics (V2 LIN; 24 parameters),

V2 with Michaelis-Menten kinetics (V2 MM; 24 pa-

rameters), and V2 with Equilibrium Chemistry

Approximation kinetics (V2 ECA; 24 parameters).

Collinearity is measured using the index

γ

with a

lower bound of 1 and an upper bound of infinity,

where a change in one parameter can be compensated

by 1-1/

γ

by changing other parameters (Section

2.4

).

Values

>

10

6

are common for models with many pa-

rameters (

Soetaert, 2016

). AGG

=

aggregates, LMWC

=

low molecular weight C, MAOM

=

miner-

al-associated organic matter, MIC

=

microbial

biomass, POM

=

particular organic matter. The error

bars represent the standard error of different collin-

earity indices estimated for 100 randomly-selected

subsets of parameters.

R.Z. Abramoff et al.

Soil Biology and Biochemistry 164 (2022) 108466

8

differently to the observed environment. Fitted Century model param-

eters were also within the same order of magnitude as the default pa-

rameters, which were taken from

and

Del Grosso

et al. (2005)

. Similar to Millennial V2, the Century model had more

absolute and relative variation in parameters controlling responses to

the environment, such as temperature and soil moisture parameters

(

Table A2

).

Millennial V2 and Century generally capture the total soil C stock

across a latitudinal gradient (

Fig. 5

Fig. S5

), although both models tend

to underestimate the amount of MAOM at high latitudes (

Fig. 6

). It is

important to note however, that while the other latitude bins have be-

tween 40 and 762 sites, the

−

50

◦

S and 70

◦

N latitude bins have only 4

and 8 sites represented, respectively. In the observed data, MAOM

makes up 73% (median; range: 2

–

100%) of the soil C stock, while in

Millennial V2 MAOM is 68% (median; range: 35

–

69%) of SOC. The

Century model does not explicitly simulate a MAOM fraction, but its two

slowest turnover pools, PASSIVE and SLOW, account for 64% (median;

range: 59

–

66%) and 32% (median; range: 29

–

33%) of the soil C stock,

respectively. The Millennial model

’

s non-MAOM pools include stable

microaggregates (or heavy POM; 24% median; range: 23

–

34%), POM

(or light POM; 8% median; range: 7

–

28%), microbial biomass C (0.1%

median; range: 0.04%

–

2.2%), and LMWC (0.1% median; range:

0.04

–

1.3%;

Fig. 6

). Literature estimates of microbial biomass C vary

widely, from 0.4% of total SOC in

Fahey et al. (2005)

to 2% in

Walker

et al. (2018)

. However, our estimate is low compared to existing

global-scale estimates of ~2% for 0

–

30 cm (

Xu et al., 2013

). LMWC

Fig. 3.

Response of each SOC pool of Millennial V2 evaluated at steady state with Michaelis-Menten kinetics to (a) soil temperature, (b) volumetric water content

(VWC) (c) net primary production (NPP), (d) clay and silt percentage, and (e) pH.

All

Viscarra Rossel

Georgiou

LUCAS

Centur

y

Millennial V2

0

10 20 30 0

10 20 30

0

10

20

30

0

10

20

30

Modeled C stock (kgC

m

−

2

)

Obser

ved C stock (kgC

m

−

2

)

Pool

MAOM
SOC

Fig. 4.

The relationship between modeled C stock by Century and Millennial Model V2 compared with observed C stock in the MAOM pool (yellow symbols) as well

as total SOC stock (blue symbols). The Century model estimates SOC only. Model fit is shown for all three datasets combined (All) as well as for each dataset

separately. The dashed line is the 1:1 line. (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)

Fig. 5.

Boxplots showing the distribution of the observed (green), Century

model-predicted (pink), and Millennial V2 model-predicted (orange) total C

stock for each 10

◦

bin of latitude. (For interpretation of the references to colour

in this figure legend, the reader is referred to the Web version of this article.)

R.Z. Abramoff et al.

Soil Biology and Biochemistry 164 (2022) 108466

9

could be analogous to dissolved organic C (DOC) if there is enough

moisture in the soil matrix, and if we do not consider DOC molecules

that are too large to be taken up by microbes. Measured DOC concen-

trations from a recent global synthesis are roughly consistent with

Millennial V2 estimates, with a range of 0.1%

–

3% of SOC depending on

the biome, with a global average of 0.2% (

Guo et al., 2020

). Including

microbial biomass and DOC measurements as constraints could improve

global predictions by informing the uptake and turnover parameters that

control exchanges between these and other modeled pools.

Millennial V2 generally captures the distribution of total SOC stocks

across different biomes (

Fig. S6

) as well as the breakdown be-

tween MAOM and non-MAOM fractions (

Fig. 8

Fig. S7

). The Millennial

V2 model predicts total SOC stocks within the range of observed values

for all but the montane biome (

), although there is much more

variation when each dataset is considered individually (

Fig. S8

). The

Century model does not estimate the C stock in each soil fraction, but

Millennial V2 is generally within the range of the observed values of C

stocks in the MAOM and non-MAOM fractions (

Fig. 8

Fig. S9

). Here it is

possible to identify the fraction which contributes to poor model per-

formance in certain biomes. For example, Millennial V2 predicts MAOM

C stocks well in tropical grasslands, but overestimates C stocks in the

non-MAOM fraction. Conversely, Millennial V2

’

s poor performance in

the montane and desert biomes is due to the model systematically

underpredicting the MAOM fraction in these biomes.

The Millennial V2 model may capture the observed variation in C

stocks across a gradient of NPP better than the Century model (

Fig. S10

Fig. 6.

Barplots showing the median total C stock that is observed, Century model-predicted, and Millennial V2 model-predicted, for each 10

◦

bin of latitude. Colors

show the distribution of C pools in the observed data (gray

=

MAOM, yellow

=

non-MAOM), Century (light blue

=

PASSIVE, orange

=

SLOW, dark red

=

ACTIVE,

pink

=

Structural Litter, purple

=

Metabolic Litter), and Millennial V2 (gray

=

MAOM, dark blue

=

Aggregates, green

=

POM, light orange

=

LMWC, black

=

microbial biomass). Metabolic Litter, LMWC, and microbial biomass pools may be too small to distinguish. (For interpretation of the references to colour in this figure

legend, the reader is referred to the Web version of this article.)

Fig. 7.

Box plots showing the distribution of observed (green), Century model-

predicted (pink), and Millennial V2 model-predicted (orange) total C stock for

each biome. Biome descriptions can be found in

Table S1

. (For interpretation of

the references to colour in this figure legend, the reader is referred to the Web

version of this article.)

Fig. 8.

Box plots showing the distribution of observed (green) and Millennial

V2 model-predicted (orange) C stock for each biome in the MAOM pool and

non-MAOM pools. The Century model is not shown because it only predicts

total C stock. Biome descriptions can be found in

Table S1

. (For interpretation

of the references to colour in this figure legend, the reader is referred to the

Web version of this article.)

R.Z. Abramoff et al.

Soil Biology and Biochemistry 164 (2022) 108466

10

This may be due to the difference in the way NPP and C stock are rep-

resented in the two models. In the Century model, there is a constrained

linear relationship between NPP and C stock (

Fig. S4c

), whereas

Millennial V2 has a less linear relationship (

c,

Fig. S3c

). The model

also provides a reasonable distribution of turnover times compared to a

global distribution of turnover times measured using SOC/respiration

rate for the top 20 cm under an assumption of steady state (

Fig. 9

).

Though the represented depth varies for different sites in our dataset,

the median depth represented in the dataset was 20 cm (6 cm

–

62 cm).

The median turnover time of Millennial V2 (17.2 years) is higher than

that of

Chen et al. (2013)

(11.2 years), but it is much more similar to the

observations than the median turnover time of the Century model (37.0

years).

Finally, we used GBM models described in Section

2.5

trained on

each of the three datasets, to quantify the performance of a purely

empirical model in contrast to process-based models like Millennial V2

and Century, as well as the extent to which the environmental forcings

(e.g., NPP, soil temperature, pH) are related to the SOC stock. We found

that the GBM model was able to predict C stocks better than the process-

based models. Similar to Millennial V2 and Century, the GBM model was

most able to predict C stocks for the sites in the VR dataset, followed

closely by those in the KG dataset, and finally those in the LUCAS dataset

(

,

). For the VR and KG datasets, the model

forcing with the greatest relative influence on the C stock was the soil

temperature. Unlike the other two datasets, C stocks in the LUCAS

dataset were not influenced by soil temperature, rather relying more

heavily on the soil pH and a combination of the other forcings (

Table 5

).

4. Discussion

4.1. Millennial V2 and Century model performance

Millennial V2 simulates spatial variation in soil C stocks, including

measurable soil fractions, better than the widely-used Century model,

based on several model performance metrics, including RMSE, AIC,

MAE, MBE and R

2

(

). Millennial V2 also predicts the C contained

within the MAOM pool, which generally contains older organic material

than other fractions, especially at depth (

;

). MAOM was generally observed to be the

largest soil C pool across all dataset sites. Because Millennial V2 more

explicitly simulates processes relating C pools to mineral capacity, pH,

temperature and soil moisture, it better predicts not only the mean C

stock across different sites (

Fig. 4

), but also better predicts the distri-

bution of C stocks across latitudes (

Fig. 5

), within different biomes

(

) and across gradients of plant productivity (

Fig. S10

). We show in

our sensitivity analyses (

c) that the Century model has a first-order

relationship relating plant inputs to rates of C gain and loss. This is a

common feature of models without explicit representation of

microorganisms (

Georgiou et al., 2017

). Millennial V2 also considers

NPP to be positively related to C stock, all else equal, but this relation-

ship is both non-linear and strongly modified by other environmental

factors, such as soil temperature and pH. As a result, the emergent

relationship between NPP and C stock is less constrained to a positive

linear relationship, as in the Century model (

Fig. S4c

), and can better

match observations (

Fig. S10

Georgiou et al., 2017

).

Millennial V2 predicts C stocks across different biomes better than

the Century model, but still has significant biases in some biomes. For

example, Millennial V2 does not capture the observed distribution of C

stocks in montane biomes (

). However, while it remains difficult to

diagnose model failures in Century, the measurable pool framework of

the Millennial model allows us to better understand what processes

cause model failures. For example, we show that Millennial V2 accu-

rately predicts the larger particle size fractions (non-MAOM) in montane

biomes, and rather that it is the MAOM fraction that is consistently

underestimated by Millennial V2 (

Fig. 8

). This may have to do with the

weathering status and mineralogy of montane environments which may

have a different relationship with clay than do other biomes. It is

important to note that the biomes for which Millennial V2 (and Century)

do not predict C stocks well are also the biomes with the fewest obser-

vations, particularly the montane biome (N

=

9). Temperate coniferous

forests (N

=

30), deserts (N

=

63), boreal forest (N

=

111), and tropical-

moist broadleaved forests (N

=

171) are also underrepresented in

certain datasets (

Fig. S8

). Therefore, mismatches in distributions of

observed and predicted C stocks could be attributed to sample size ef-

fects, understanding that C stocks have a high spatial heterogeneity

within biomes.

The choice of dataset affects the accuracy of both Millennial V2 and

Century, not only due to sample size effects but due to the source and

quality of ancillary data used. We used a GBM model to empirically

quantify the relationship between the environmental factors used as

model forcing (i.e., soil temperature, clay and silt content, NPP, soil

moisture, pH) and C stocks (

Table 5

). Using this approach, it is clear that

there is less of a relationship between environmental factors at LUCAS

sites and the measured soil C stock, with only 15% of cross-validated

variance explained, compared to over 60% for the other two datasets.

The LUCAS sites also rank the influence of environmental factors on C

stock differently than the other datasets, with pH being by far the most

influential factor at LUCAS sites, compared to soil temperature at other

sites. This may explain why Millennial V2 only explains 4% of variation

in C stocks at LUCAS sites (compared to 31% and 46% for the other two

datasets), and why a different microbial-mineral model (

Robertson

et al., 2019

) fit to LUCAS site data has a similar coefficient of determi-

nation, although models are generally able to capture mean C stocks

across biome (

Figs. 7, 8

) and land type (

Fig. S13

) classifications. It is not

immediately clear what may cause greater heterogeneity in the LUCAS

dataset, but the fact that SOC stocks are more dependent on pH than soil

temperature suggests a potential role of land use history and past

management practices.

Past and current land management may be a source of uncertainty

across all the datasets, as well as unmet assumptions of steady state. For

example, the KG dataset contains 284 cropland sites, which are likely to

have changing C stocks. When we repeat our analysis excluding these

sites (

Table S2

,

Table S3

), some model performance metrics improved

while others worsened (

Table S4

,

Table S5

). In general, the coefficients

of determination improved while metrics based on error worsened.

Although we did not see clear evidence for model improvement by

excluding these sites, the difference in model performance supports the

idea that land use history and management may have important effects

on the spatial variation in C stocks that cannot be captured by models

that do not represent these processes.

4.2. Comparing process-based models to empirical models

The GBM model predicts steady-state C stocks with higher accuracy

Fig. 9.

Density plots showing the distribution of turnover times in

Chen et al.

https://github.com/rabramoff/Millennial

(denoted Observed; median [1st, 3rd quartiles]

=

11.2 [7.4, 17.5]), the

Century model (37.0 [31.6, 45.7]), and the Millennial V2 model (17.2

[11.6, 23.5]).

R.Z. Abramoff et al.

Soil Biology and Biochemistry 164 (2022) 108466

11

than both the Millennial V2 and Century models, raising the question,

why not replace process-based models entirely with empirical models?

One reason is that a process-based model can be more generalizable. For

example, when we used parameters that were fit to one dataset to pre-

dict stocks from a different dataset (

Table S6

), Millennial V2 was able to

explain 3

–

46% of variation in the other two datasets. Conversely, the

GBM model could only explain 0

–

22% of variation in the other datasets

(

Table S6

). Therefore, Millennial V2 predictions made by model pa-

rameters that have been fit to a subset of global-scale data are likely

more predictive than those made by an empirical model fit to the same

data. Because process-based models are more generalizable, they may be

more reliable under novel conditions, whereas empirical models are

most useful for interpolation. The same principle can explain why a

more process-rich model like Millennial V2 appears to outperform a

more empirical model such as Century (

Fig. 4

).

Another benefit of process-based models is the ability to predict

quantities other than the one that has been fitted. For example, by fitting

the model to C stocks, we also predicted reasonable turnover times,

which are related to both the stock size and the emission (respiration)

rate. Because C emissions are arguably the most important model output

for understanding climate feedbacks from soil, a model which uses

physical understanding to link stocks and emissions is inherently more

valuable than an empirical model of only C stocks. Though not tested

here, previous work on microbial soil models has also shown that

process-based models can capture seasonal hysteresis in C fluxes where

empirical models (with no explicit representation of time) cannot

(

Abramoff et al., 2017

). Yet, empirical models may be able to improve

process-based models by accelerating spin-up and parameterization, or

by replacing certain model components.

4.3. Future data and model needs

The large variation in model performance depending on the sample

size or the data source underscores the need for abundant and stan-

dardized data. Across these datasets, MAOM was measured using three

different particle size thresholds (

<

50, 53, and 60

μ

m), and POM and

aggregates were not separated in the larger particle size fraction. Mea-

surements of small C pools that are sensitive to environmental factors,

such as microbial biomass and low molecular weight C, were largely not

available. Although microbial biomass is generally less than 5% of total

SOC (

Xu et al., 2013

), it is the primary driver of C emissions from the

soil, and small changes to this pool may have large effects on total SOC

and C emissions. Recent work has also shown that DOC has a large po-

tential to contribute to C storage, especially in the mineral fraction

(

Abramoff et al., 2021

). Therefore, large, standardized datasets of

Millennial V2 pools (

Fig. 1

) and fluxes (e.g., respiration rate, enzyme

activity) will allow for better global-scale parameterization of the

model.

Although Millennial V2 is more complex than many soil models, it

has structural limitations that would benefit from future work, such as

consideration of nitrogen (N) and phosphorus (P) limitations to C stor-

age (

Davies et al., 2020

;

Spohn, 2020

), and depth-dependence of mi-

crobial activity and substrate availability (

Dove et al., 2021

). We further

acknowledge that this study evaluates the spatial variation in C stock

only, and that the temporal dynamics of Millennial V2 require further

testing. Several long-term experiments exist across the world (

Richter

et al., 2007

), but most have not made measurements of soil fractions at

different time periods. Future measurement campaigns, such the Joint

Research Center

’

s next planned sampling campaign for the LUCAS sites,

will add much needed temporal resolution to the existing dataset.

Recent and developing spectroscopy-based methods may also allow for

low-cost estimates of soil fractions from ongoing experiments as well as

archived soils (

;

), creating datasets of repeat measurements for model eval-

uation. C models are constantly evolving to represent the most updated

knowledge of the soil system. We now seek to associate that knowledge

more directly with measurable quantities, both to make models easier to

constrain and also to ensure that model predictions are accurate for the

right reasons, rather than the result of compensating biases. We hope

that Millennial V2 and similarly representative models will allow for

more realistic multiple constraints on the next generation of soil

biogeochemical models.

Code and data availability

The code for the models tested here are available on the github re-

pository:

. Millennial V1 and

V2 are available in both R and Fortran. We have also included the scripts

used to produce the figures and analyses presented in this manuscript as

R markdown notebooks within the same repository:

https://github.

com/rabramoff/Millennial/tree/master/R/analysis

. Some small data

files are also included in the same folder. Other datasets are described in

the manuscript with citations and when available their URLs are also

provided in the analysis code.

Declaration of competing interest

The authors declare that they have no known competing financial

interests or personal relationships that could have appeared to influence

the work reported in this paper.

Acknowledgements

R.Z.A. thanks Jinyun Tang and Elisa Th

´

ebault for helpful discussion

about model formulation. R.Z.A. was supported by the French govern-

ment grant

“

Make Our Planet Great Again

”

and by a Marie Sk

ł

o-

dowska

–

Curie Individual Fellowship (Grant no. 834-169) from the

European Union

’

s Horizon 2020 program. K.G. was supported as a

Lawrence Fellow at Lawrence Livermore National Lab (LLNL) by the

LLNL-LDRD Program under Project No. 21-ERD-045 and the US DOE

Office of Science, Office of Biological and Environmental Research,

Genomic Science Program as part of the LLNL Microbes Persist Scientific

Focus Area, SCW1632. Work at LLNL was conducted under the auspices

of US DOE Contract DE-AC52-07NA27344. X.X. has been supported by

the ORNL Terrestrial Ecosystem Science Scientific Focus Area (ORNL

TES-SFA) and NGEE Arctic projects and DE-SC0014416, which are

supported by the Office of Biological and Environmental Research in the

Department of Energy Office of Science.

Table 5

Performance metrics of the gradient-boosted machine learning algorithm prediction of soil C stock and the relative influence of different model forcings for the
different datasets.

R

2

in

=

coefficient of determination of the training dataset (i.e., in-sample),

R

2

out

=

coefficient of determination of the test dataset (i.e., out-of-sample),

RMSE

=

root mean square error, MAE

=

mean absolute error, Soil Temp

=

soil temperature, NPP

=

net primary production.

Relative Influence (%)

Dataset

R

2

in

R

2

out

RMSE (kg C m

−

2

)

MAE (kg C m

−

2

)

Soil Temp (

◦

C)

pH

Clay

&

Silt (%)

Soil Moisture (mm

3

mm

¡

3

)

NPP (g C m

¡

2

d

¡

1

)

All

0.66

0.48

2.74

1.79

38.8

26.5

13.5

11.8

9.48

VR

0.83

0.68

1.21

0.87

40.0

21.2

16.4

7.33

15.1

KG

0.78

0.61

2.34

1.59

31.0

16.9

14.7

25.0

12.3

LUCAS

0.34

0.15

4.76

3.72

7.03

40.2

19.6

11.1

22.1

R.Z. Abramoff et al.

Soil Biology and Biochemistry 164 (2022) 108466

12

Appendix A. Supplementary data

Supplementary data to this article can be found online at

https://doi.org/10.1016/j.soilbio.2021.108466

.

Appendices.

Appendix A

Table A1

Parameters, constants, pools, fluxes and input variables used in the Millennial model V2.

“

Calculated

”

variables are pools or fluxes whose values are a function of

parameters, initial values, and time.

Equation

Number

Variable

Definition

Default

values

Fit Values:

All

Fit Values:

VR

Fit Values:

KG

Fit Values:

LUCAS

Units

Sources

1, 7

p

i

Proportion of C input allocated to

POM

0.66

-

Oleson et al.

1, 19

p

a

Proportion of aggregate-C

breakdown allocated to POM

0.33

0.40

0.34

0.33

-

Oleson et al.

2

K

pl

Half-saturation constant of POM

decomposition to LMWC

10,000

6,443

8,617

7,735

12,094

g C m

−

2

Wang et al.

(2014)

3

a

pl

Pre-exponential constant for V

pl

2.5

×

10

12

1.8

×

10

12

2.6

×

10

12

2.3

×

10

12

1.8

×

10

12

g C m

−

2

(g C

m

−

2

)

−

1

d

−

1

3

Ea

pl

Activation energy for V

pl

64,320

63,909

63,339

64,064

64,646

J mol

−

1

3

R

Gas constant

8.31446

-

J K

−

1

mol

−

1

4, 15

φ

total porosity

0.60

0.62

0.61

0.60

mm

3

mm

−

3

site-specific

5

k

pa

Rate of aggregate formation from

POM

0.020

0.018

0.012

0.014

0.026

d

−

1

Segoli et al.

6

k

b

Breakdown rate of soil aggregate

carbon

0.019

0.020

0.026

0.023

0.015

d

−

1

Segoli et al.

7, 19, 20

p

b

Partitioning of necromass to MAOM

and LMWC

0.50

0.52

0.50

0.61

-

8

k

l

Leaching rate of LMWC

0.0015

-

d

−

1

10, 12

K

ld

Specific desorption rate for LMWC

and MIC

1

-

mg C L

−

1

d

−

1

Wang et al.

(2014)

10

p

1

Coefficient for estimating the

binding affinity for LMWC sorption

0.186

0.12

0.078

0.21

0.38

-

Mayes et al.

(2012)

10

p

2

Coefficient for estimating the

binding affinity for LMWC sorption

0.216

-

Mayes et al.

(2012)

10

pH

7

-

site-specific

11

BD

Bulk density

1000

-

kg soil m

−

3

site-specific

11

p

c

Coefficient for estimating the

maximum sorption capacity

0.86

-

Georgiou et al.

(2021)

11

%

claysilt

Percent of soil that is in the clay and

silt fractions

80

-

%

site-specific

13

K

lb

Half saturation constant for

microbial uptake

290

774.6

710.8

654.8

100.3

g C m

−

2

14

a

lb

Pre-exponential constant for V

lb

2.6

×

10

12

2.3

×

10

12

1.2

×

10

12

2.2

×

10

12

7.2

×

10

12

g C m

−

2

(g C

m

−

2

)

−

1

d

−

1

14

Ea

lb

Activation energy for V

lb

60,260

57,865

60,428

60,058

57,795

J mol

−

1

15

λ

Dependence of rate on matric

potential

2.1

×

10

−

4

-

kPa

−

1

Ghezzehei et al.

15

k

a,min

Minimum relative rate in saturated

soil

0.2

-

Ghezzehei et al.

15

φ

matric potential

−

15

-

kPa

site-specific

16

k

bd

Microbial death rate

0.0036

0.0045

0.0044

0.0040

0.0036

m

2

gC

−

1

d

−

1

18

k

ma

Rate of aggregate formation from

MAOM

0.020

0.0048

0.0052

0.0038

0.0052

d

−

1

Segoli et al.

21, 22

CUE

ref

Reference CUE

0.60

0.19

0.53

0.40

0.52

-

21, 22

CUE

T

CUE dependence on temperature

0.012

-

◦

C

−

1

21, 22

T

ae-ref

Reference temperature for

temperature control on CUE

15

-

◦

C

1, 2, 5

P

POM

g C m

−

2

Calculated

7, 8, 9, 13

L

LMWC

g C m

−

2

Calculated

2, 13, 16, 20

B

Microbial biomass

g C m

−

2

Calculated

9, 12, 18, 19

M

MAOM

g C m

−

2

Calculated

6, 17

A

Aggregate C

g C m

−

2

Calculated

1, 6, 17, 19

F

a

Aggregate breakdown

g C m

−

2

d

−

1

Calculated

1, 5, 17

F

pa

Aggregate formation from POM

g C m

−

2

d

−

1

Calculated

1, 2, 7

F

pl

Decomposition of POM into LMWC

g C m

−

2

d

−

1

Calculated

7, 8

F

l

LMWC leaching loss

g C m

−

2

d

−

1

Calculated

(

continued on next page

)

R.Z. Abramoff et al.

Soil Biology and Biochemistry 164 (2022) 108466

13

Table A1

(

continued

)

Equation

Number

Variable

Definition

Default

values

Fit Values:

All

Fit Values:

VR

Fit Values:

KG

Fit Values:

LUCAS

Units

Sources

7, 9, 19

F

lm

Adsorption of LMWC to minerals

g C m

−

2

d

−

1

Calculated

7, 13, 20, 21,

22

F

lb

Uptake of LMWC by microbial

biomass

g C m

−

2

d

−

1

Calculated

7, 16, 19, 20

F

bm

Microbial mortality

g C m

−

2

d

−

1

Calculated

12, 19

F

ld

Desorption

g C m

−

2

d

−

1

Calculated

22

F

bg

Microbial growth

g C m

−

2

d

−

1

Calculated

17, 18, 19

F

ma

Aggregate formation from MAOM

g C m

−

2

d

−

1

Calculated

20, 22, 23

F

mr

Microbial respiration

g C m

−

2

d

−

1

Calculated

2, 3

V

pl

Maximum rate of POM

decomposition to LMWC

d

−

1

Calculated

13, 15

S

w,B

Water scalar

-

Calculated

2, 4, 5, 6, 8, 9,

15,18

S

w,D

Diffusion limitation of substrates

-

Calculated

9, 10

K

lm

LMWC and microbial necromass

adsorption rate

d

−

1

Calculated

9, 11

Q

max

Maximum sorption capacity

g C m

−

2

Calculated

13, 14

V

lb

Potential LMWC uptake rate

d

−

1

Calculated

1, 7

F

i

C input

g C m

−

2

d

−

1

Input

3, 16, 21, 22

T

Soil temperature

◦

C

Input

4, 15

θ

volumetric water content

mm

3

mm

−

3

Input

23

CO

2

carbon dioxide production

g C m

−

2

d

−

1

Output

Table A2

Parameters, pools, fluxes and input variables used in the Century model.

“

Calculated

”

variables are pools or fluxes whose values are a function of parameters, initial

values, and time.

Variable

Definition

Default

values

Fit Values:

All

Fit Values:

VR

Fit Values:

KG

Fit Values:

LUCAS

Units

Sources

w1

Water scalar parameter

30.0

33.61

17.38

36.97

22.61

-

Parton et al. (2010)

w2

Water scalar parameter

9.00

8.42

10.55

7.77

7.84

-

Parton et al. (2010)

t1

x-axis location of inflection point

15.4

18.08

16.83

15.77

18.15

◦

C

t2

y-axis location of inflection point

11.75

12.99

10.39

9.32

13.78

-

t3

Distance from the maximum point to the

minimum point (step size)

29.7

25.92

30.10

20.98

28.91

-

t4

Slope of line at inflection point

0.031

0.038

0.025

0.032

0.035

-

c1

Intercept of clay fraction relationship

0.85

0.83

0.85

0.86

-

c2

Slope of clay fraction relationship

0.68

0.67

0.71

0.74

0.60

-

k

ls

Turnover rate of structural litter pool

0.01

-

d

−

1

k

lm

Turnover rate of metabolic litter pool

0.045

-

d

−

1

k

a

Turnover rate of active pool

0.020

-

d

−

1

k

s

Turnover rate of slow pool

5.0

×

10

−

4

6.7

×

10

−

4

1.9

×

10

−

4

5.2

×

10

−

4

4.0

×

10

−

4

d

−

1

k

p

Turnover rate of passive pool

2.0

×

10

−

5

9.1

×

10

−

6

8.2

×

10

−

6

1.9

×

10

−

5

1.3

×

10

−

5

d

−

1

p

li

Proportion of plant residue to structural

litter pool

0.66

0.62

0.67

0.64

0.67

-

Analogous to Millennial litter

partitioning assumption

p

lma

Fraction of metabolic litter to active pool

0.45

-

p

lsa

Fraction of structural litter to active pool

0.5

-

p

lss

Fraction of structural litter to slow pool

0.7

-

p

sa

Fraction of slow pool to active pool

0.42

-

p

sp

Fraction of slow pool to passive pool

0.03

0.020

0.029

0.022

0.041

-

p

pa

Fraction of passive pool to active pool

0.45

-

p

ap

Fraction of active pool to passive pool

0.004

-

LigFrac

Fraction of litter that is lignin

0.2

0.20

0.21

0.20

-

Zhang et al. (2018)

%claysilt

Percent of soil that is in the clay and silt

fractions

80

-

%

site-specific

Fc

Field capacity

0.39

mm

3

mm

−

3

Input: observed maximum

volumetric water content

F

i

C input

g C m

−

2

d

−

1

Input

T

Soil temperature

◦

C

Input

θ

volumetric water content

mm

3

mm

−

3

Input

StrLitter

Structural litter pool

g C m

−

2

Calculated

MetLitter

Metabolic litter pool

g C m

−

2

Calculated

ACTIVE

Active pool

g C m

−

2

Calculated

SLOW

Slow pool

g C m

−

2

Calculated

PASSIVE

Passive pool

g C m

−

2

Calculated

S

t

Temperature scalar

-

Calculated

S

w

Water scalar

-

Calculated

F

tex

Soil texture effect on decomposition of

active pool

-

Calculated

F

ls

Structural litter decomposition flux

Calculated

(

continued on next page

)

R.Z. Abramoff et al.

Soil Biology and Biochemistry 164 (2022) 108466

14

Table A2

(

continued

)

Variable

Definition

Default

values

Fit Values:

All

Fit Values:

VR

Fit Values:

KG

Fit Values:

LUCAS

Units

Sources

g C m

−

2

d

−

1

F

lm

Metabolic litter decomposition flux

g C m

−

2

d

−

1

Calculated

F

a

Active pool decomposition flux

g C m

−

2

d

−

1

Calculated

F

s

Slow pool decomposition flux

g C m

−

2

d

−

1

Calculated

F

p

Passive pool decomposition flux

g C m

−

2

d

−

1

Calculated

Appendix B

Century model equations

Model inputs are identical in description and units to those used in the Millennial model (

https://doi.org/10.1007/s10533-017-0409-7

). All other parameters are described in

Table A2

,

including their default and fitted values for each dataset.

S

t

=

t

2

+

t

3

π

atan

[

t

4

*

π

(

T

−

t

1

)]

t

2

+

t

3

π

atan

[

t

4

*

π

(

30

−

t

1

)]

[B1]

S

w

=

1

+

w

1

*

e

−

w

2

θ

fc

[B2]

F

tex

=

c

1

−

c

2 %

claysilt

0

.

01

[B3]

F

ls

=

StrLitter k

ls

S

t

S

w

e

−

3

LigFrac

[B4]

F

lm

=

MetLitter k

lm

S

t

S

w

[B5]

F

a

=

ACTIVE k

a

S

t

S

w

F

tex

[B6]

F

s

=

SLOW k

s

S

t

S

w

[B7]

F

p

=

PASSIVE k

p

S

t

S

w

[B8]

dStrLitter

dt

=

p

li

F

i

−

F

ls

[B9]

dMetLitter

dt

= (

1

−

p

li

)

F

i

−

F

lm

[B10]

dACTIVE

dt

= (

1

−

LigFrac

)

p

lsa

F

ls

+

p

lma

F

lm

+

F

s

p

sa

+

F

p

pa

−

F

a

[B11]

dSLOW

dt

=

LigFrac p

lss

F

ls

+

F

a

(

1

−

F

tex

−

p

ap

)

−

F

s

[B12]

dPASSIVE

dt

=

F

a

p

ap

+

F

s

p

sp

−

F

p

[B13]

References

Abramoff, R., Xu, X., Hartman, M., et al., 2018. The Millennial model: in search of

measurable pools and transformations for modeling soil carbon in the new century.

Biogeochemistry 137, 51

–

71.

.

Abramoff, R.Z., Davidson, E.A., Finzi, A.C., 2017. A parsimonious modular approach to

building a mechanistic belowground carbon and nitrogen model. Journal of

Geophysical Research Biogeosciences 1

–

17.

https://doi.org/10.1002/

2017JG003796

.

Abramoff, R.Z., Georgiou, K., Guenet, B., et al., 2021. How much carbon can be added to

soil by sorption? Biogeochemistry Letters.

https://doi.org/10.1007/s10533-021-

00759-x

. BIOGEOCHEMISTRY.

Ahrens, B., Braakhekke, M.C., Guggenberger, G., et al., 2015. Contribution of sorption,

DOC transport and microbial interactions to the 14C age of a soil organic carbon

profile: insights from a calibrated process model. Soil Biology and Biochemistry 88,

390

–

402.

https://doi.org/10.1016/j.soilbio.2015.06.008

Arrhenius, S., 1889. Über die Reaktionsgeschwindigkeit bei der Inversion von

Allison, S.D., Wallenstein, M.D., Bradford, M a, 2010. Soil-carbon response to warming

dependent on microbial physiology. Nature Geoscience 3, 336

–

Rohrzucker durch S

¨

auren. Zeitschrift für Physikalische Chemie 4, 226

248

Bailey, C.J., 1989. Enzyme kinetics of cellulose hydrolysis. Biochemical Journal 262,

1001

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

1002

.

Bailey, V.L., Bond-Lamberty, B., DeAngelis, K., et al., 2018. Soil carbon cycling proxies:

understanding their critical role in predicting climate change feedbacks. Global

Change Biology 24, 895

–

905.

https://doi.org/10.1071/SR12374

Baldock, J.A., Sanderman, J., MacDonald, L.M., et al., 2013. Quantifying the allocation

of soil organic carbon to biologically significant fractions. Soil Research 51,

561

–

576.

Beaudoing, H., Rodell, M., 2016. GLDAS Noah Land Surface Model L4 Monthly 1 X 1

Degree V2, p. 1

https://doi.org/10.1016/j.gloplacha.2012.10.006

R.Z. Abramoff et al.

Soil Biology and Biochemistry 164 (2022) 108466

15

Chen, S., Huang, Y., Zou, J., Shi, Y., 2013. Mean residence time of global topsoil organic

carbon depends on temperature, precipitation and soil nitrogen. Global and

Planetary Change 100, 99

–

108.

.

Conen, F., Zimmermann, M., Leifeld, J., et al., 2008. Relative stability of soil carbon

revealed by shifts in

δ

15N and C:N ratio. Biogeosciences 5, 123

https://doi.org/10.1038/s41561-019-0484-6

128

.

Cotrufo, M.F., Ranalli, M.G., Haddix, M.L., et al., 2019. Soil carbon storage informed by

particulate and mineral-associated organic matter. Nature Geoscience 12, 989

–

994.

Davidson, E.A., Samanta, S., Caramori, S.S., Savage, K., 2012. The Dual Arrhenius and

Michaelis

Menten kinetics model for decomposition of soil organic matter at hourly

to seasonal time scales. Global Change Biology 18, 371

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

384

.

Davies, C.A., Robertson, A.D., McNamara, N.P., 2020. The importance of nitrogen for net

carbon sequestration when considering natural climate solutions. Global Change

Biology 1

–

2.

https://doi.org/10.1007/s10533-004-

Del Grosso, S.J., Parton, W.J., Mosier, A.R., et al., 2005. Modeling soil CO2 emissions

from ecosystems. Biogeochemistry 73, 71

–

91.

0898-z

.

Dove, N.C., Torn, M.S., Hart, S.C., Tas¸, N., 2021. Metabolic capabilities mute positive

response to direct and indirect impacts of warming throughout the soil profile.

Nature Communications 12, 1

–

13.

https://doi.org/10.1038/s41467-021-22408-5

.

Fahey, T.J., Siccama, T.G., Driscoll, C.T., et al., 2005. The biogeochemistry of carbon at

Hubbard Brook. Biogeochemistry 75, 109

–

176.

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

004-6321-y

Friedman, J.H., 2001. Greedy function approximation: a gradient boosting machine.

Annals of Statistics 29, 1189

1232

https://doi.org/10.1038/s41467-017-01116-z

Georgiou, K., Abramoff, R.Z., Harte, J., et al., 2017. Microbial community-level

regulation explains soil carbon responses to long-term litter manipulations. Nature

Communications.

Georgiou, K., Jackson, R.B., Abramoff, R.Z., et al., 2021. Global capacity and controls of

mineral-associated carbon in soils. Nature Geoscience under review

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

German, D.P., Marcelo, K.R.B., Stone, M.M., Allison, S.D., 2012. The Michaelis-Menten

kinetics of soil extracellular enzymes in response to temperature: a cross-latitudinal

study. Global Change Biology 18, 1468

–

1479.

2486.2011.02615.x

.

Ghezzehei, T.A., Sulman, B., Arnold, C.L., et al., 2018. On the role of soil water retention

characteristic on aerobic microbial respiration. Biogeosciences Discussions 16,

1187

–

1209.

https://doi.org/10.5194/bg-2018-265

https://doi.org/10.1111/geb.13186

Guo, Z., Wang, Y., Wan, Z., et al., 2020. Soil dissolved organic carbon in terrestrial

ecosystems: global budget, spatial distribution and controls. Global Ecology and

Biogeography 29, 2159

–

2175.

https://doi.org/10.5194/bg-10-2011-2013

Haverd, V., Raupach, M.R., Briggs, P.R., et al., 2013. Multiple observation types reduce

uncertainty in Australia

’

s terrestrial carbon and water cycles. Biogeosciences 10,

2011

–

2040.

Hengl, T., Jesus, JM De, Heuvelink, G.B.M., et al., 2017. SoilGrids250m : Global Gridded

Soil Information Based on Machine Learning

Hicks Pries, C.E., 2017. The whole-soil carbon flux in response to warming. Science

Ito, A., Hajima, T., Lawrence, D.M., et al., 2020. Soil carbon sequestration simulated in

CMIP6-LUMIP models: implications for climatic mitigation. Environmental Research

Letters.

https://doi.org/10.1088/1748-9326/abc912

.

Jenkinson, D.S., Rayner, J.H., 1977. The turnover of soil organic matter in some of the

Rothamsted classical experiments. Soil Science 123, 298

https://doi.org/10.1038/ncomms13630

305

.

Kallenbach, C.M., Grandy, A., Frey, S.D., 2016. Direct evidence for microbial-derived soil

organic matter formation and its ecophysiological controls. Nature Communications

1

–

10.

Lajtha, K., Bailey, V.L., McFarlane, K., et al., 2018. Ch. 12: soils. In: Cavallaro, N.,

Shrestha, G., Birdsey, R., et al. (Eds.), Second State of the Carbon Cycle Report

(SOCCR2): A Sustained Assessment Report. U.S. Global Change Research Program,

Washington, DC, USA, pp. 469

506

Lehmann, J., Hansel, C.M., Kaiser, C., et al., 2020. Persistence of soil organic carbon

caused by functional complexity. Nature Geoscience 13, 529

–

534.

10.1038/s41561-020-0612-3

.

Li, S., Yuan, W., Ciais, P., et al., 2019. Benchmark estimates for aboveground litterfall

data derived from ecosystem models. Environmental Research Letters.

org/10.1088/1748-9326/ab2ee4

.

Luo, Y., Ahlstr

¨

om, A., Allison, S.D., et al., 2016. Toward more realistic projections of soil

carbon dynamics by Earth system models. Global Biogeochemical Cycles 30, 40

–

56.

https://doi.org/10.1002/2015GB005239

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

Marschmann, G.L., Pagel, H., Kügler, P., Streck, T., 2019. Equifinality , sloppiness , and

emergent structures of mechanistic soil biogeochemical models. Environmental

Modelling

&

Software 122, 104518.

envsoft.2019.104518

https://doi.org/10.1029/2004GB002219

Masiello, C.A., Chadwick, O.A., Southon, J., et al., 2004. Weathering controls on

mechanisms of carbon storage in grassland soils. Global Biogeochemical Cycles 18,

1

–

9.

https://doi.org/10.2136/sssaj

Mayes, M.A., Heal, K.R., Brandt, C.C., et al., 2012. Relation between soil order and

sorption of dissolved organic carbon in temperate soils. Soil Science Society of

America Journal 76, 61

–

69.

Michaelis, L., Menten, M.L., 1913. The kinetics of the inversion effect. Biochemische

Zeitschrift 49, 333

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

369

.

Mikutta, R., Zang, U., Chorover, J., et al., 2011. Stabilization of extracellular polymeric

substances (Bacillus subtilis) by adsorption to and coprecipitation with Al forms.

Geochimica et Cosmochimica Acta 75, 3135

–

3154.

gca.2011.03.006

Oleson, K.W., Lawrence, D.M., Bonan, G.B., et al., 2013. Technical Description of Version

4.5 of the Community Land Model (CLM). NCAR Tech., National Center for

Atmospheric Research, Bounder

Olson, D.M., Dinerstein, E., Wikramanayake, E.D., et al., 2001. Terrestrial ecoregions of

the world: a new map of life on Earth. BioScience 51, 933

938

https://doi.org/10.1016/j.landusepol.2011.07.003

Panagos, P., Van Liedekerke, M., Jonaes, A., Montanarella, L., 2012. European Soil Data

Centre: response to European policy support and public data requirements. Land Use

Policy 29, 329

–

338.

Parton, W.J., Schimel, D.S., Cole, C.V., et al., 1987. Analysis of factors controlling soil

Parton, W.J., Hanson, P.J., Swanston, C., et al., 2010. ForCent model development and

testing using the Enriched Background Isotope Study experiment. Journal of

Geophysical Research Biogeosciences 115, 1

–

15.

https://doi.org/10.1029/

2009JG001193

.

organic matter levels in Great plains grasslands. Soil Science Society of America

Journal 51, 1173

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

1179

.

Poeplau, C., Don, A., Six, J., et al., 2018. Isolating organic carbon fractions with varying

turnover rates in temperate agricultural soils

–

a comprehensive method comparison.

Soil Biology and Biochemistry 125, 10

–

26.

soilbio.2018.06.025

Poeplau, C., K

¨

atterer, T., Leblans, N.I.W., Sigurdsson, B.D., 2017. Sensitivity of soil

carbon fractions and their specific stabilization mechanisms to extreme soil warming

in a subarctic grassland. Global Change Biology 23, 1316

–

1327.

https://doi.org/10.3389/fenvs.2021.634472

10.1111/gcb.13491

.

Ramírez, P.B., Calder

´

on, F.J., Haddix, M., et al., 2021. Using diffuse reflectance

spectroscopy as a high throughput method for quantifying soil C and N and their

distribution in particulate and mineral-associated organic matter fractions. Frontiers

of Environmental Science 9, 1

–

13.

https://doi.org/10.2136/sssaj2006.0181

Richter, D. deB., Hofmockel, M., Callaham, M.A., et al., 2007. Long-term soil

experiments: keys to managing earth

’

s rapidly changing ecosystems. Soil Science

Society of America Journal 71, 266

–

279.

https://doi.org/10.5194/bg-2018-430

Robertson, A.D., Paustian, K., Ogle, S., et al., 2019. Unifying soil organic matter

formation and persistence frameworks: the MEMS model. Biogeosciences Discussions

16, 1225

–

124836.

https://doi.org/10.1016/S0304-3800(97)00067-7

Rodrigo, A., Recous, S., Neel, C., Mary, B., 1997. Modelling temperature and moisture

effects on C-N transformations in soils: comparison of nine models. Ecological

Modelling 102, 325

–

339.

.

Sanderman, J., Baldock, J.A., Dangal, S.R.S., et al., 2021. Soil organic carbon fractions in

the Great Plains of the United States: an application of mid-infrared spectroscopy.

Biogeochemistry.

https://doi.org/10.1007/s10533-021-00755-1

.

Schmidt, M., Torn, M.S., Abiven, S., et al., 2011. Persistence of soil organic matter as an

ecosystem property. Nature 478, 49

–

56.

https://doi.org/10.1038/nature10386

https://doi.org/10.1016/j.ecolmodel.2013.04.010

Segoli, M., De Gryze, S., Dou, F., et al., 2013. AggModel: a soil organic matter model with

measurable pools for use in incubation studies. Ecological Modelling 263, 1

–

9.

https://doi.org/10.1038/s41467-018-04526-9

Shi, Z., Crowell, S., Luo, Y., Moore, B., 2018. Model structures amplify uncertainty in

predicted soil carbon responses to climate change. Nature Communications 9, 1

–

11.

https://doi.org/10.1007/s10533-011-9596-9

Sierra, C.A., 2012. Temperature sensitivity of organic matter decomposition in the

Arrhenius equation: some theoretical considerations. Biogeochemistry 108, 1

–

15.

.

Sierra, C., Muller, M., 2015. A general mathematical framework for representing soil

organic matter dynamics. Ecological Monographs 16, 16881.

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

10.1890/15-0361.1

.

Sierra, C.A., Muller, M., Metzler, H., et al., 2017. The muddle of ages , turnover , transit ,

and residence times in the carbon cycle. Global Change Biology 23, 1763

–

1773.

Soetaert, K., 2009. rootSolve: Nonlinear Root Finding, Equilibrium and Steady-State

Analysis of Ordinary Differential Equations

Soetaert, K., 2016. R Package FME: Inverse Modelling, Sensitivity, Monte Carlo

to a Dynamic Simulation Model

Applied

Soetaert, K., Herman, P.M.J., 2009. A Practical Guide to Ecological Modelling. Using R as

a Simulation Platform. Springer

.

Spawn, S.A., Sullivan, C.C., Lark, T.J., Gibbs, H.K., 2020. Harmonized global maps of

above and belowground biomass carbon density in the year 2010. Scientific Data 7,

Spohn, M., 2020. Increasing the organic carbon stocks in mineral soils sequesters large

amounts of phosphorus. Global Change Biology 4169

–

4177.

https://doi.org/10.1016/j.geoderma.2013.06.014

10.1111/gcb.15154

.

Stamati, F.E., Nikolaidis, N.P., Banwart, S., Blum, W.E.H., 2013. A coupled carbon ,

aggregation , and structure turnover (CAST) model for topsoils. Geoderma 212,

51

–

64.

.

Sulman, B., Moore, J., Abramoff, R., et al., 2018. Multiple models and experiments

underscore large uncertainty in soil carbon dynamics. Biogeochemistry 141,

109

123

.

Tang, J., Riley, W.J., 2015. Weaker soil carbon-climate feedbacks resulting from

microbial and abiotic interactions. Nature Climate Change 5, 56

–

60.

https://doi.org/10.1016/j.soilbio.2019.107624

org/10.1038/nclimate2438

.

Tang, J., Riley, W.J., 2019. Competitor and substrate sizes and diffusion together define

enzymatic depolymerization and microbial substrate uptake rates. Soil Biology and

Biochemistry 139, 107624.

.

Tang, J.Y., 2015. On the relationships between the Michaelis

–

Menten kinetics , reverse

Michaelis

–

Menten kinetics , equilibrium chemistry approximation kinetics , and

quadratic kinetics. Geoscientific Model Development 8, 3823

–

3835.

10.5194/gmd-8-3823-2015

https://doi.org/10.5194/bg-10-

Tang, J.Y., Riley, W.J., 2013. A total quasi-steady-state formulation of substrate uptake

kinetics in complex networks and an example application to microbial litter

R.Z. Abramoff et al.

Soil Biology and Biochemistry 164 (2022) 108466

16

decomposition. Biogeosciences 10, 8329

–

8351.

8329-2013

.

Tisdall, J., Oades, J., 1982. Organic matter and water-stable aggregates in soils. Journal

of Soil Science 33, 141

–

163.

https://doi.org/10.1111/j.1365-2389.1982.tb01755.x

.

Torn, M.M.S., Trumbore, S.S.E., Chadwick, O.A.O., et al., 1997. Mineral control of soil

organic carbon storage and turnover. Nature 389, 3601

–

3603.

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

10.1038/38260

.

T

´

oth, G., Jones, A., Montanarella, L., 2013. The LUCAS topsoil database and derived

information on the regional variability of cropland topsoil properties in the

European Union. Environmental Monitoring and Assessment 185, 7409

–

7425.

Van Veen, J.A., Kuikman, P.J., 1990. Soil structural aspects of decomposition of organic

matter by micro-organisms. Biogeochemistry 11, 213

–

233.

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

10.1007/BF00004497

.

Viscarra Rossel, R., Hicks, W., 2015. Soil organic carbon and its fractions estimated by

visible

–

near infrared transfer functions. European Journal of Soil Science 438

–

450.

https://doi.org/10.1071/SR14366

Viscarra Rossel, R.A., Chen, C., Grundy, M.J., et al., 2015. The Australian three-

dimensional soil grid: Australia

’

s contribution to the GlobalSoilMap project. Soil

Research 53, 845

–

864.

.

Viscarra Rossel, R.A., Lee, J., Behrens, T., et al., 2019. Continental-scale soil carbon

composition and vulnerability modulated by regional environmental controls.

Nature Geoscience.

https://doi.org/10.1038/s41561-019-0373-z

https://doi.org/10.1016/j.soilbio.2007.03.007

von Lützow, M., K

¨

ogel-Knabner, I., Ekschmitt, K., et al., 2007. SOM fractionation

methods: relevance to functional pools and to stabilization mechanisms. Soil Biology

and Biochemistry 39, 2183

–

2207.

.

Walker, T.W.N., Kaiser, C., Strasser, F., et al., 2018. Microbial temperature sensitivity

and biomass change explain soil carbon loss with warming. Nature Climate Change

1.

https://doi.org/10.1038/s41558-018-0259-x

https://doi.org/10.1038/ismej.2014.120

Wang, G., Jagadamma, S., Mayes, M.A., et al., 2014. Microbial dormancy improves

development and experimental validation of ecosystem model. The ISME Journal 9,

226

–

237.

https://doi.org/10.1890/12-0681.1

Wang, G., Post, W.M., Mayes, M.A., 2013. Development of microbial-enzyme-mediated

decomposition model parameters through steady-state and dynamic analyses.

Ecological Applications 23, 255

–

272.

.

Weng, Z., Van Zwieten, L., Singh, B.P., et al., 2017. Biochar built soil carbon over a

decade by stabilizing rhizodeposits. Nature Climate Change 7, 371

https://doi.org/10.1002/2015GB005188

376

.

Wieder, W.R., Allison, S.D., Davidson, E.A., et al., 2015a. Explicitly representing soil

microbial processes in Earth system models. Global Biogeochemical Cycles

1782

–

1800.

. Received.

Wieder, W.R., Grandy, A.S., Kallenbach, C.M., et al., 2015b. Representing life in the

Earth system with soil microbial functional traits in the MIMICS model. Geoscientific

Model Development 8, 1789

–

1808.

https://doi.org/10.5194/gmd-8-1789-2015

Wieder, W.R., Grandy, A.S., Kallenbach, C.M., Bonan, G.B., 2014. Integrating microbial

physiology and physio-chemical principles in soils with the MIcrobial-MIneral

Carbon Stabilization (MIMICS) model. Biogeosciences 11, 3899

–

3917.

org/10.5194/bg-11-3899-2014

Wieder, W.R., Hartman, M.D., Sulman, B.N., et al., 2017. Carbon Cycle Confidence and

Uncertainty : Exploring Variation Among Soil Biogeochemical Models.

org/10.1111/gcb.13979

https://doi.org/10.1029/2018JG004777

Xia, J., Yuan, W., Lienert, S., et al., 2019. Global patterns in net primary production

allocation regulated by environmental conditions and forest stand age: a model-data

comparison. Journal of Geophysical Research Biogeosciences 124, 2039

–

2059.

.

Xu, W., Chang, J., Ciais, P., et al., 2020. Reducing uncertainties of future global soil

carbon responses to climate and land use change with emergent constraints. Global

Biogeochemical Cycles 34, 1

–

25.

https://doi.org/10.1029/2020GB006589

https://doi.org/10.1111/geb.12029

Xu, X., Thornton, P.E., Post, W.M., 2013. A global analysis of soil microbial biomass

carbon, nitrogen and phosphorus in terrestrial ecosystems. Global Ecology and

Biogeography 22, 737

–

749.

.

Zhang, H., Goll, D., Manzoni, S., et al., 2018. Modeling the effects of litter stoichiometry

and soil mineral N availability on soil organic matter formation using CENTURY-

CUE (v1.0). Geoscientific Model Development 11, 4779

–

4796.