# Improved global-scale predictions of soil carbon stocks with Millennial Version 2

Rose Z Abramoff, Bertrand Guenet, Haicheng Zhang, Katerina Georgiou, Xiaofeng Xu, Raphael Viscarra Rossel, Wenping Yuan, Philippe Ciais

<figure><img src="image_1.png" /></figure>

To cite this version:

Rose Z Abramoff, Bertrand Guenet, Haicheng Zhang, Katerina Georgiou, Xiaofeng Xu, et al. Improved global-scale predictions of soil carbon stocks with Millennial Version 2. Soil Biology and Biochemistry, 2022, 164, pp.108466. 10.1016/j.soilbio.2021.108466. hal-03419469

HAL Id: hal-03419469

https://hal.science/hal-03419469

Submitted on 15 Nov 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Improved global-scale predictions of soil carbon stocks with Millennial
Version 2

Rose Z. Abramoff<sup>a, *</sup>, Bertrand Guenet<sup>b</sup>, Haicheng Zhang<sup>c</sup>, Katerina Georgiou<sup>d,e</sup>, Xiaofeng Xu<sup>f</sup>,
Raphael A. Viscarra Rossel<sup>g</sup>, Wenping Yuan<sup>h</sup>, Philippe Ciais<sup>a</sup>

ARTICLE INFO

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 understanding 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⁻², AIC = 675, R²<sub>in</sub> = 0.31, R²<sub>out</sub> = 0.26) better than the widely-used first-order decomposition model Century (RMSE = 3.4 kg C m⁻², AIC = 696, R²<sub>in</sub> = 0.21, R²<sub>out</sub> = 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.

<footnote>ª Laboratoire des Sciences du Climat et de l'Environnement, LSCE/IPSL, CEA-CNRS-UVSQ, Université Paris-Saclay, F-91191, Gif-sur-Yvette, France</footnote>

<footnote>ᵇ Laboratoire de Géologie, UMR 8538, Ecole Normale Supérieure, IPSL, PSL Research University, CNRS, Paris, France</footnote>

<footnote>ᶜ Université Libre des Bruxelles, 1050, Bruxelles, Belgium</footnote>

<footnote n="4">⁴ Department of Earth System Science, Stanford University, Yang & Yamazaki Building, 473 Via Ortega, Stanford, CA, 94305-4216, USA</footnote>

<footnote n="5">⁵ Physical and Life Sciences Directorate, Lawrence Livermore National Laboratory, Livermore, CA, 94551, USA</footnote>

<footnote>f Biology Department, San Diego State University, San Diego, CA, 92130, USA</footnote>

<footnote n="8">⁸ Soil and Landscape Science, School of Molecular and Life Sciences, Curtin University, GPO Box U1987, Perth WA, 6845, Australia</footnote>

<footnote>h School of Atmospheric Sciences, Sun Yat-Sen University, Zhuhai, 510245, Guangdong, China</footnote>

<footnote>* Corresponding author.</footnote>

<footnote>E-mail addresses: rose.abramoff@lsce.ipsl.fr (R.Z. Abramoff), betrand.guenet@ens.fr (B. Guenet), Haicheng.Zhang@ulb.ac.be (H. Zhang), georgiou@stanford.edu, georgiou1@llnl.gov (K. Georgiou), xxu@sdsu.edu (X. Xu), r.viscarra-rossel@curtin.edu.au (R.A. Viscarra Rossel), yuanwp3@mail.sysu.edu.cn (W. Yuan), philippe.ciais@lsce.ipsl.fr (P. Ciais).</footnote>
# 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 microorganisms, and by physical processes controlling soil water and temperature (Lajtha et al., 2018). Soil decomposition models are used to assess current 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 (Luo et al., 2016; Wieder et al., 2017; Shi et al., 2018; Ito et al., 2020; Xu et al., 2020). 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 represent 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 definitions (Jenkinson and Rayner, 1977). Our best understanding of biological, 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 (Allison et al., 2010; German et al., 2012; Wieder et al., 2014) and stabilization of SOM via aggregation (Segoli et al., 2013; Stamati et al., 2013) or mineral association (Wang et al., 2013; Ahrens et al., 2015; Tang and Riley, 2015), reflecting field and laboratory studies emphasizing the importance 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 (Schmidt et al., 2011; Lehmann et al., 2020), yet there is less consensus about the best mathematical formulation of different processes (Wieder et al., 2015a; Sulman et al., 2018). 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 substrate greatly exceeds the concentration of the enzyme (forward Michaelis-Menten) or vice versa (reverse Michaelis-Menten) (Michaelis and Menten, 1913; Bailey, 1989). Equilibrium Chemistry Approximation (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, 2019). 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; 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 immobilize OM (Masiello et al., 2004; Mikutta et al., 2011; Weng et al., 2017), including aggregation (Van Veen and Kuikman, 1990; von Lützow et al.,

2007).

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 microbial 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 underlying 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 Century model (Del Grosso et al., 2005; Parton et al., 2010). In this paper, we maintain the conceptual model as presented in Abramoff et al. (2018), but update the governing equations to reflect recent developments 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 between Millennial V1 and V2 are summarized in Table 1. Specifically, the temperature function was updated from the Century temperature scalar used in Del Grosso et al. (2005) 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 limitation, 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 dissolved 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 parameters, 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 Table A1. 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,
<caption>**Table 1**</caption>

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 Table A1 for definitions of variables, their units, and values.

<table><thead><tr><th>Process</th><th>Millennial V1 [Equation Number in Abramoff et al., 2018]</th><th>Millennial V2 [Equation Number in this paper (Section 2.2)]</th></tr></thead><tbody><tr><td>Temperature function</td><td>S<sub>t</sub> = <div style="display: inline-block; vertical-align: middle; text-align: center;"><div style="border-bottom: 1px solid black; padding: 0 5px;">t<sub>2</sub> + (<div style="display: inline-block; vertical-align: middle; text-align: center;"><div style="border-bottom: 1px solid black; padding: 0 5px;">t<sub>3</sub></div><div style="padding: 0 5px;">π</div></div>atan[π(T - t<sub>1</sub>)]</div> [3]</td><td>V<sub>x</sub> = α<sub>x</sub>e<sup>-E<sub>α,x</sub>/[R(T+273.15)]</sup> [3,14]</td></tr><tr><td>Moisture function</td><td>S<sub>w</sub> = <div style="display: inline-block; vertical-align: middle; text-align: center;"><div style="border-bottom: 1px solid black; padding: 0 5px;">t<sub>2</sub> + (<div style="display: inline-block; vertical-align: middle; text-align: center;"><div style="border-bottom: 1px solid black; padding: 0 5px;">t<sub>3</sub></div><div style="padding: 0 5px;">π</div></div>atan[π t<sub>4</sub> (T<sub>ref</sub> - t<sub>1</sub>)]</div> 1</td><td>S<sub>w</sub> = e<sup>J<sub>pl</sub></sup> k<sub>a,min</sub> + (1 - k<sub>a,min</sub>) (<div style="display: inline-block; vertical-align: middle; text-align: center;"><div style="border-bottom: 1px solid black; padding: 0 5px;">θ</div><div style="padding: 0 5px;">φ</div></div>)<sup>0.5</sup> (<div style="display: inline-block; vertical-align: middle; text-align: center;"><div style="border-bottom: 1px solid black; padding: 0 5px;">θ</div><div style="padding: 0 5px;">φ</div></div>)<sup>0.5</sup> [4,15]</td></tr><tr><td>Sorption to minerals</td><td>Q<sub>max</sub> = BD 10<sup>cl</sup> log<sub>(%clay)</sub> + c<sub>2</sub> [11]</td><td>Q<sub>max</sub> = depth BD %claysilt p<sub>c</sub> [11]</td></tr><tr><td>Microbial mortality</td><td>F<sub>bm</sub> = k<sub>mm</sub>S<sub>w</sub>B [17]</td><td>F<sub>bm</sub> = k<sub>bd</sub>B<sup>2</sup> [16]</td></tr></tbody></table>

<figure><img src="image_3.png" /><figcaption>**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 Table A1 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.</figcaption></figure>

$$dP/dt = p_i F_i + p_a F_a - F_{pa} - F_{pl} \quad [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 proportion 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, 2019),

$$F_{pl} = V_{pl} S_{w,D} P \frac{B}{K_{pl} + B} \quad [2]$$

$$V_{pl} = \alpha_{pl} e^{-E_{a,pl} / [R(T+273.15)]} \quad [3]$$

$$S_{w,D} = \left(\frac{\theta}{\varphi}\right)^{0.5} \quad [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 $\alpha_{pl}$, an activation energy $E_{a,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 limitation of substrates defined in Ghezzehei et al. (2018) as the square root of the ratio of the volumetric water content ($\theta$) and the total porosity ($\varphi$).

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 \quad [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 \quad [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,

$$\frac{dL}{dt} = F_i(1 - p_i) - F_l + F_{pl} - F_{lm} - F_{lb} + (1 - p_b)F_{bm} + F_{ld} \quad [7]$$

where $F_l$ is the LMWC leaching loss,

$$F_l = k_l S_{w,D} L \quad [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. (2013),

$$F_{lm} = S_{w,D} K_{lm} L \left(1 - \frac{M}{Q_{max}}\right) \quad [9]$$

where $K_{lm}$ is the binding affinity that is adjustable based on the pH. $Q_{max}$ is the maximum sorption capacity (gC/m²), assuming a 1 m soil profile. $K_{lm}$ is further defined by the parameters $p_1$, $p_2$, and the desorption coefficient ($K_{ld}$),

$$K_{lm} = e^{-p_1 p H - p_2} K_{ld} \quad [10]$$

where parameters $p_1$ and $p_2$ are the coefficients for computing the first term ($e^{-p_1 p H - 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
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 \quad [11]$$

Desorption of LMWC is a function of MAOM, $Q_{max}$, and the desorption coefficient ($K_{ld}$) sensu Wang et al. (2013).

$$F_{ld} = K_{ld} \frac{M}{Q_{max}} \quad [12]$$

Microbial uptake of LMWC is a function of microbial biomass and LMWC concentration, temperature, and water,

$$F_{lb} = V_{lb}S_{w,B}B \frac{L}{K_{lb} + L} \quad [13]$$

$$V_{lb} = \alpha_{lb}e^{-Ea_{lb}/[R(T+273.15)]} \quad [14]$$

$$S_{w,B} = e^{\lambda\phi} \left[ k_{a,min} + (1 - k_{a,min}) \left( \frac{\phi - \theta}{\phi} \right)^{0.5} \right] S_{w,D} \quad [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 (Ghezzeheiri et al., 2018), where $\lambda$ is the dependence of the rate on the matric potential ($\phi$) and $k_{a,min}$ is the minimum relative rate in saturated soil. This relationship incorporates the effects of matric potential ($e^{\lambda\phi}$), oxygen limitation $\left[k_{a,min} + (1 - k_{a,min}) \left(\frac{\phi - \theta}{\phi}\right)^{0.5}\right]$ 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 \quad [16]$$

Both MAOM and POM can enter the aggregate C pool (A),

$$dA / dt = F_{ma} + F_{pa} - F_a \quad [17]$$

$$F_{ma} = k_{ma}S_{w,D}M \quad [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) \quad [19]$$

where $F_{bm}$ is the carbon flow from microbial biomass to MAOM, namely buildup of microbial necromass. The partitioning of microbial necromass 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 microbial necromass that forms mineral associations by other means, accounting 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} \quad [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})]\} \quad [21]$$

$$F_{bg} = F_{lb} [CUE_{ref} - CUE_{T}(T - T_{ae-ref})] \quad [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} \quad [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 reaction kinetics approximations in soil microbial decomposition models (Allison et al., 2010; Davidson et al., 2012; German et al., 2012; Wang et al., 2013; Wieder et al., 2014). 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 \frac{B}{K_{pl} + B + P} \quad [2a]$$

$$F_{lb} = V_{lb}S_{w,B}B \frac{L}{K_{lb} + L + B} \quad [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 \frac{B}{K_{pl}} \quad [2b]$$

$$F_{lb} = V_{lb}S_{w,B}B \frac{L}{K_{lb}} \quad [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,
$$S_{ij} = \frac{dy_i}{d\theta_j} \frac{\theta_j}{y_i} \qquad [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 Sij with the collin function of R's FME Version 1.3.6.1 package, and is more fully described in Soetaert (2016), Abramoff et al. (2017), and Marschmann et al. (2019). Very high values can be common with multiple parameters; for example, in a 5-parameter model of bacterial growth described in Soetaert (2016), the collinearity index was 2.4 × 10⁶ 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 g C m⁻² d⁻¹], soil temperature [-3.47 – 29.9 °C], volumetric water content [0.10–0.48 m³m⁻³], 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 includes three soil pools (active, slow and passive) and two litter pools (structural and metabolic). The equations governing interactions between the five Century model pools are defined in Parton et al. (1987), 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. Performance metrics for this model include the out-of-sample R² 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).

## 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 purposes, 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 Viscarra Rossel and Hicks (2015) 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⁻². SOC was fractionated into three fractions, originally reported as particulate 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 ¹³C nuclear magnetic resonance spectroscopy (Baldock et al., 2013). 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 productivity (NPP) derived from the land surface model BIOS2 (Haverd et al., 2013). Soil temperature (°C) and soil moisture (mm³ mm⁻³) 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 (Beaudoin and Rodell, 2016). 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 concentrations were reported in g C kg⁻¹ soil and converted to g C m⁻² 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
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 production from Li et al. (2019) and the root production estimated from Li et al. (2019) and spawn et al. (2020) 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óth 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., 2019). These measurements were used to evaluate the MEMS model, another model proposing measurable soil C pools (Robertson et al., 2019). 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 Table A1 for parameters and their fitted values). Parameters were optimized by minimizing the sum of squared residuals between the model and observations 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²) 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

<caption>**Table 2**</caption>

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

<table><thead><tr><th>Dataset</th><th>Region</th><th>Particle Size Fractions</th><th>NPP</th><th>Soil temperature &amp; moisture</th><th>Clay &amp; Silt</th><th>pH</th><th>Bulk Density</th></tr></thead><tbody><tr><td>VR</td><td>Australia</td><td>(Viscarra Rossel and Hicks, 2015; Viscarra Rossel et al., 2019)</td><td>BIOS2 (Haverd et al., 2013)</td><td>GLDAS (Beaudoing and Rodell, 2016)</td><td>Viscarra Rossel et al. (2015)</td><td>Viscarra Rossel et al. (2015)</td><td>Viscarra Rossel et al. (2015)</td></tr><tr><td>KG</td><td>Global</td><td>Georgiou et al. (2021)</td><td>(Li et al., 2019; Spawn et al., 2020)</td><td>GLDAS (Beaudoing and Rodell, 2016)</td><td>Georgiou et al. (2021)</td><td>SoilGrids (Hengl et al., 2017)</td><td>SoilGrids (Hengl et al., 2017)</td></tr><tr><td>LUCAS</td><td>EU</td><td>(Panagos et al., 2012; Cotrufo et al., 2019)</td><td>(Li et al., 2019; Spawn et al., 2020)</td><td>GLDAS (Beaudoing and Rodell, 2016)</td><td>Tóth et al. (2013)</td><td>Tóth et al. (2013)</td><td>SoilGrids (Hengl et al., 2017)</td></tr></tbody></table>

$$RMSE = \sqrt{\frac{\sum_{i=1}^{n} (SOC_{obs,i} - SOC_{mod,i})^2}{n}} \quad [25]$$

Where SOC<sub>obs,i</sub> is the observed SOC at each site i, SOC<sub>mod,i</sub> 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 = \frac{\sum_{i=1}^{n} |SOC_{obs,i} - SOC_{mod,i}|}{n} \quad [26]$$

$$MBE = \frac{\sum_{i=1}^{n} (SOC_{obs,i} - SOC_{mod,i})}{n} \quad [27]$$

We calculated the Akaike information criterion (AIC) as a performance metric, which considers not only the error but also the number of model parameters,

$$AIC = n \times \ln \left( \sqrt{\frac{\sum_{i=1}^{n} (SOC_{obs,i} - SOC_{mod,i})^2}{n}} \right) + 2p \quad [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 variables, we fit the model to the entirety of each dataset. This is also how we calculate the in-sample R², denoted R²<sub>in</sub> 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²<sub>out</sub> (out-of-sample R²), 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

<caption>**Table 3**</caption>

Model fit of SOC to measured values. R²<sub>in</sub> = coefficient of determination of the training dataset (i.e., in-sample), R²<sub>out</sub> = 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.

<table><thead><tr><th>Dataset</th><th>Model</th><th>R<sup>2</sup><sub>in</sub></th><th>R<sup>2</sup><sub>out</sub></th><th>RMSE (kg C m<sup>-2</sup>)</th><th>AIC</th><th>MAE (kg C m<sup>-2</sup>)</th><th>MBE (kg C m<sup>-2</sup>)</th></tr></thead><tbody><tr><td rowspan="2">All (N = 1329)</td><td>Millennial V2</td><td>0.31</td><td>0.26</td><td>3.29</td><td>675</td><td>2.14</td><td>-1.01</td></tr><tr><td>Century</td><td>0.21</td><td>0.18</td><td>3.42</td><td>696</td><td>2.30</td><td>-0.72</td></tr><tr><td rowspan="2">VR (N = 495)</td><td>Millennial V2</td><td>0.46</td><td>0.37</td><td>1.83</td><td>163</td><td>1.25</td><td>-0.63</td></tr><tr><td>Century</td><td>0.40</td><td>0.32</td><td>2.11</td><td>189</td><td>1.46</td><td>-1.17</td></tr><tr><td rowspan="2">KG (N = 659)</td><td>Millennial V2</td><td>0.32</td><td>0.21</td><td>3.46</td><td>372</td><td>2.28</td><td>-1.18</td></tr><tr><td>Century</td><td>0.20</td><td>0.02</td><td>3.81</td><td>398</td><td>2.63</td><td>-0.35</td></tr><tr><td rowspan="2">LUCAS (N = 175)</td><td>Millennial V2</td><td>0.04</td><td>0.04</td><td>5.03</td><td>150</td><td>3.71</td><td>-1.10</td></tr><tr><td>Century</td><td>0.04</td><td>0.01</td><td>5.21</td><td>152</td><td>3.77</td><td>-1.19</td></tr></tbody></table>
<caption>**Table 4**</caption>

Millennial V2 model fit of MAOM and non-MAOM (POM + AGG + LMWC + MIC) pools to measured values. $R_{in}^2$ = coefficient of determination of the training dataset (i.e., in-sample), $R_{out}^2$ = 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.

<table><thead><tr><th>Dataset</th><th>Pool</th><th>N</th><th>R<sub>in</sub><sup>2</sup></th><th>R<sub>out</sub><sup>2</sup></th><th>MAE (kg C m<sup>-2</sup>)</th><th>MBE (kg C m<sup>-2</sup>)</th></tr></thead><tbody><tr><td rowspan="2">All</td><td>MAOM</td><td>1329</td><td>0.40</td><td>0.27</td><td>1.58</td><td>-1.14</td></tr><tr><td>non-MAOM</td><td>670</td><td>0.17</td><td>0.24</td><td>0.82</td><td>-0.077</td></tr><tr><td rowspan="2">VR</td><td>MAOM</td><td>495</td><td>0.42</td><td>0.37</td><td>0.89</td><td>-0.39</td></tr><tr><td>non-MAOM</td><td>495</td><td>0.29</td><td>0.16</td><td>0.52</td><td>-0.24</td></tr><tr><td>KG</td><td>MAOM</td><td>659</td><td>0.40</td><td>0.31</td><td>1.38</td><td>-0.62</td></tr><tr><td rowspan="2">LUCAS</td><td>MAOM</td><td>175</td><td>0.15</td><td>0.09</td><td>1.99</td><td>-0.77</td></tr><tr><td>non-MAOM</td><td>175</td><td>NS</td><td>NS</td><td>2.20</td><td>-0.33</td></tr></tbody></table>

used the World Wildlife Fund Terrestrial Ecoregions Map which classifies 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 dataset, 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/respiration 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 parameters, 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). However, 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 identifiable 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 (Fig. 3a). 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 (Fig. 3b). 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 (Fig. 3c, 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⁻², AIC = 675, $R_{in}^2$ = 0.31, $R_{out}^2$ = 0.26) better than the Century model (RMSE = 3.4 kg C m⁻², AIC = 696, $R_{in}^2$ = 0.21, $R_{out}^2$ = 0.18; Table 3) across sites. The Millennial V2 model also outperformed the Century model when fit to each dataset separately (Table 3, 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 3, 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 (Table A1). 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

<figure><img src="image_4.png" /><figcaption>**Fig. 2.** Mean collinearity of parameters for four model variants, Millennial V1 (Original; 23 parameters), V2 with Linear kinetics (V2 LIN; 24 parameters), V2 with Michaelis-Menten kinetics (V2 MM; 24 parameters), and V2 with Equilibrium Chemistry Approximation kinetics (V2 ECA; 24 parameters). Collinearity is measured using the index $\gamma$ with a lower bound of 1 and an upper bound of infinity, where a change in one parameter can be compensated by 1-1/$\gamma$ by changing other parameters (Section 2.4). Values >10⁶ are common for models with many parameters (Soetaert, 2016). AGG = aggregates, LMWC = low molecular weight C, MAOM = mineral-associated organic matter, MIC = microbial biomass, POM = particular organic matter. The error bars represent the standard error of different collinearity indices estimated for 100 randomly-selected subsets of parameters.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>**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.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**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.)</figcaption></figure>

differently to the observed environment. Fitted Century model parameters were also within the same order of magnitude as the default parameters, which were taken from Parton et al. (1987) 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 between 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

<figure><img src="image_7.png" /><figcaption>**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.)</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>**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.)</figcaption></figure>

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 concentrations 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. 7, Fig. S6) as well as the breakdown between 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 (Fig. 7), 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

<figure><img src="image_9.png" /><figcaption>**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.)</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**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.)</figcaption></figure>

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 performance 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).
This may be due to the difference in the way NPP and C stock are represented 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 (Fig. 3c, 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 (Table 5, Fig. S11, Fig. S12). 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$ (Table 3). Millennial V2 also predicts the C contained within the MAOM pool, which generally contains older organic material than other fractions, especially at depth (Conen et al., 2008; Hicks Pries, 2017; Poeplau et al., 2018). 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 distribution of C stocks across latitudes (Fig. 5), within different biomes (Fig. 7) and across gradients of plant productivity (Fig. S10). We show in our sensitivity analyses (Fig. 3c) 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 relationship 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 (Fig. 7). 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 accurately 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 observations, 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 effects, 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 determination, 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

<figure><img src="image_11.png" /><figcaption>**Fig. 9.** Density plots showing the distribution of turnover times in Chen et al. (2013) (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]).</figcaption></figure>
<caption>**Table 5**</caption>

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_{in}^2$ = coefficient of determination of the training dataset (i.e., in-sample), $R_{out}^2$ = 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.

<table><thead><tr><th rowspan="2">Dataset</th><th rowspan="2">R<sup>2</sup><sub>in</sub></th><th rowspan="2">R<sup>2</sup><sub>out</sub></th><th rowspan="2">RMSE (kg C m<sup>-2</sup>)</th><th rowspan="2">MAE (kg C m<sup>-2</sup>)</th><th colspan="5">Relative Influence (%)</th></tr><tr><th>Soil Temp (°C)</th><th>pH</th><th>Clay & Silt (%)</th><th>Soil Moisture (mm<sup>3</sup> mm<sup>-3</sup>)</th><th>NPP (g C m<sup>-2</sup> d<sup>-1</sup>)</th></tr></thead><tbody><tr><td>All</td><td>0.66</td><td>0.48</td><td>2.74</td><td>1.79</td><td>38.8</td><td>26.5</td><td>13.5</td><td>11.8</td><td>9.48</td></tr><tr><td>VR</td><td>0.83</td><td>0.68</td><td>1.21</td><td>0.87</td><td>40.0</td><td>21.2</td><td>16.4</td><td>7.33</td><td>15.1</td></tr><tr><td>KG</td><td>0.78</td><td>0.61</td><td>2.34</td><td>1.59</td><td>31.0</td><td>16.9</td><td>14.7</td><td>25.0</td><td>12.3</td></tr><tr><td>LUCAS</td><td>0.34</td><td>0.15</td><td>4.76</td><td>3.72</td><td>7.03</td><td>40.2</td><td>19.6</td><td>11.1</td><td>22.1</td></tr></tbody></table>

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 predict 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 parameters 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 standardized 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. Measurements 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 potential 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 storage (Davies et al., 2020; Spohn, 2020), and depth-dependence of microbial 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 as 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 (Baldock et al., 2013; Ramírez et al., 2021; Sanderman et al., 2021), creating datasets of repeat measurements for model evaluation. 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 repository: https://github.com/rabramoff/Millennial. 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ébault for helpful discussion about model formulation. R.Z.A. was supported by the French government grant "Make Our Planet Great Again" and by a Marie Skłodowska-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.
# 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

<caption>**Table A1**</caption>

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.

<table><thead><tr><th>Equation Number</th><th>Variable</th><th>Definition</th><th>Default values</th><th>Fit Values: All</th><th>Fit Values: VR</th><th>Fit Values: KG</th><th>Fit Values: LUCAS</th><th>Units</th><th>Sources</th></tr></thead><tbody><tr><td>1, 7</td><td>p<sub>l</sub></td><td>Proportion of C input allocated to POM</td><td>0.66</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>Oleson et al. (2013)</td></tr><tr><td>1, 19</td><td>p<sub>a</sub></td><td>Proportion of aggregate-C breakdown allocated to POM</td><td>0.33</td><td>0.33</td><td>0.40</td><td>0.34</td><td>0.33</td><td>-</td><td>Oleson et al. (2013)</td></tr><tr><td>2</td><td>K<sub>pl</sub></td><td>Half-saturation constant of POM decomposition to LMWC</td><td>10,000</td><td>6,443</td><td>8,617</td><td>7,735</td><td>12,094</td><td>g C m<sup>-2</sup></td><td>Wang et al. (2014)</td></tr><tr><td>3</td><td>a<sub>pl</sub></td><td>Pre-exponential constant for V<sub>pl</sub></td><td>2.5 &times; 10<sup>12</sup></td><td>1.8 &times; 10<sup>12</sup></td><td>2.6 &times; 10<sup>12</sup></td><td>2.3 &times; 10<sup>12</sup></td><td>1.8 &times; 10<sup>12</sup></td><td>g C m<sup>-2</sup> (g C m<sup>-2</sup>)<sup>-1</sup> d<sup>-1</sup></td><td>Abramoff et al. (2017)</td></tr><tr><td>3</td><td>E<sub>pl</sub></td><td>Activation energy for V<sub>pl</sub></td><td>64,320</td><td>63,909</td><td>63,339</td><td>64,064</td><td>64,646</td><td>J mol<sup>-1</sup></td><td>Abramoff et al. (2017)</td></tr><tr><td>3</td><td>R</td><td>Gas constant</td><td>8.31446</td><td>-</td><td>-</td><td>-</td><td>-</td><td>J K<sup>-1</sup> mol<sup>-1</sup></td><td>Abramoff et al. (2017)</td></tr><tr><td>4, 15</td><td>&phi;</td><td>total porosity</td><td>0.60</td><td>0.60</td><td>0.62</td><td>0.61</td><td>0.60</td><td>mm<sup>3</sup> mm<sup>-3</sup></td><td>site-specific</td></tr><tr><td>5</td><td>k<sub>pa</sub></td><td>Rate of aggregate formation from POM</td><td>0.020</td><td>0.018</td><td>0.012</td><td>0.014</td><td>0.026</td><td>d<sup>-1</sup></td><td>Segoli et al. (2013)</td></tr><tr><td>6</td><td>k<sub>b</sub></td><td>Breakdown rate of soil aggregate carbon</td><td>0.019</td><td>0.020</td><td>0.026</td><td>0.023</td><td>0.015</td><td>d<sup>-1</sup></td><td>Segoli et al. (2013)</td></tr><tr><td>7, 19, 20</td><td>p<sub>b</sub></td><td>Partitioning of necromass to MAOM and LMWC</td><td>0.50</td><td>0.50</td><td>0.52</td><td>0.50</td><td>0.61</td><td>-</td><td>-</td></tr><tr><td>8</td><td>k<sub>l</sub></td><td>Leaching rate of LMWC</td><td>0.0015</td><td>-</td><td>-</td><td>-</td><td>-</td><td>d<sup>-1</sup></td><td>Abramoff et al. (2018)</td></tr><tr><td>10, 12</td><td>K<sub>ld</sub></td><td>Specific desorption rate for LMWC and MIC</td><td>1</td><td>-</td><td>-</td><td>-</td><td>-</td><td>mg C L<sup>-1</sup> d<sup>-1</sup></td><td>Wang et al. (2014)</td></tr><tr><td>10</td><td>p<sub>1</sub></td><td>Coefficient for estimating the binding affinity for LMWC sorption</td><td>0.186</td><td>0.12</td><td>0.078</td><td>0.21</td><td>0.38</td><td>-</td><td>Mayes et al. (2012)</td></tr><tr><td>10</td><td>p<sub>2</sub></td><td>Coefficient for estimating the binding affinity for LMWC sorption</td><td>0.216</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>Mayes et al. (2012)</td></tr><tr><td>10</td><td>pH</td><td>pH</td><td>7</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>site-specific</td></tr><tr><td>11</td><td>BD</td><td>Bulk density</td><td>1000</td><td>-</td><td>-</td><td>-</td><td>-</td><td>kg soil m<sup>-3</sup></td><td>site-specific</td></tr><tr><td>11</td><td>p<sub>c</sub></td><td>Coefficient for estimating the maximum sorption capacity</td><td>0.86</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>Georgiou et al. (2021)</td></tr><tr><td>11</td><td>% claysilt</td><td>Percent of soil that is in the clay and silt fractions</td><td>80</td><td>-</td><td>-</td><td>-</td><td>-</td><td>%</td><td>site-specific</td></tr><tr><td>13</td><td>K<sub>lb</sub></td><td>Half saturation constant for microbial uptake</td><td>290</td><td>774.6</td><td>710.8</td><td>654.8</td><td>100.3</td><td>g C m<sup>-2</sup></td><td>Abramoff et al. (2017)</td></tr><tr><td>14</td><td>a<sub>lb</sub></td><td>Pre-exponential constant for V<sub>lb</sub></td><td>2.6 &times; 10<sup>12</sup></td><td>2.3 &times; 10<sup>12</sup></td><td>1.2 &times; 10<sup>12</sup></td><td>2.2 &times; 10<sup>12</sup></td><td>7.2 &times; 10<sup>12</sup></td><td>g C m<sup>-2</sup> (g C m<sup>-2</sup>)<sup>-1</sup> d<sup>-1</sup></td><td>Abramoff et al. (2017)</td></tr><tr><td>14</td><td>E<sub>lb</sub></td><td>Activation energy for V<sub>lb</sub></td><td>60,260</td><td>57,865</td><td>60,428</td><td>60,058</td><td>57,795</td><td>J mol<sup>-1</sup></td><td>Abramoff et al. (2017)</td></tr><tr><td>15</td><td>&lambda;</td><td>Dependence of rate on matric potential</td><td>2.1 &times; 10<sup>-4</sup></td><td>-</td><td>-</td><td>-</td><td>-</td><td>kPa<sup>-1</sup></td><td>Ghezzhei et al. (2018)</td></tr><tr><td>15</td><td>k<sub>a,min</sub></td><td>Minimum relative rate in saturated soil</td><td>0.2</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>Ghezzhei et al. (2018)</td></tr><tr><td>15</td><td>&phi;</td><td>matric potential</td><td>-15</td><td>-</td><td>-</td><td>-</td><td>-</td><td>kPa</td><td>site-specific</td></tr><tr><td>16</td><td>k<sub>bd</sub></td><td>Microbial death rate</td><td>0.0036</td><td>0.0045</td><td>0.0044</td><td>0.0040</td><td>0.0036</td><td>m<sup>2</sup> gC<sup>-1</sup> d<sup>-1</sup></td><td>Abramoff et al. (2017)</td></tr><tr><td>18</td><td>k<sub>ma</sub></td><td>Rate of aggregate formation from MAOM</td><td>0.020</td><td>0.0048</td><td>0.0052</td><td>0.0038</td><td>0.0052</td><td>d<sup>-1</sup></td><td>Segoli et al. (2013)</td></tr><tr><td>21, 22</td><td>CUE<sub>ref</sub></td><td>Reference CUE</td><td>0.60</td><td>0.19</td><td>0.53</td><td>0.40</td><td>0.52</td><td>-</td><td></td></tr><tr><td>21, 22</td><td>CUE<sub>T</sub></td><td>CUE dependence on temperature</td><td>0.012</td><td>-</td><td>-</td><td>-</td><td>-</td><td>&deg;C<sup>-1</sup></td><td></td></tr><tr><td>21, 22</td><td>T<sub>ae-ref</sub></td><td>Reference temperature for temperature control on CUE</td><td>15</td><td>-</td><td>-</td><td>-</td><td>-</td><td>&deg;C</td><td></td></tr><tr><td>1, 2, 5</td><td>P</td><td>POM</td><td></td><td></td><td></td><td></td><td></td><td>g C m<sup>-2</sup></td><td>Calculated</td></tr><tr><td>7, 8, 9, 13</td><td>L</td><td>LMWC</td><td></td><td></td><td></td><td></td><td></td><td>g C m<sup>-2</sup></td><td>Calculated</td></tr><tr><td>2, 13, 16, 20</td><td>B</td><td>Microbial biomass</td><td></td><td></td><td></td><td></td><td></td><td>g C m<sup>-2</sup></td><td>Calculated</td></tr><tr><td>9, 12, 18, 19</td><td>M</td><td>MAOM</td><td></td><td></td><td></td><td></td><td></td><td>g C m<sup>-2</sup></td><td>Calculated</td></tr><tr><td>6, 17</td><td>A</td><td>Aggregate C</td><td></td><td></td><td></td><td></td><td></td><td>g C m<sup>-2</sup></td><td>Calculated</td></tr><tr><td>1, 6, 17, 19</td><td>F<sub>a</sub></td><td>Aggregate breakdown</td><td></td><td></td><td></td><td></td><td></td><td>g C m<sup>-2</sup> d<sup>-1</sup></td><td>Calculated</td></tr><tr><td>1, 5, 17</td><td>F<sub>pa</sub></td><td>Aggregate formation from POM</td><td></td><td></td><td></td><td></td><td></td><td>g C m<sup>-2</sup> d<sup>-1</sup></td><td>Calculated</td></tr><tr><td>1, 2, 7</td><td>F<sub>pl</sub></td><td>Decomposition of POM into LMWC</td><td></td><td></td><td></td><td></td><td></td><td>g C m<sup>-2</sup> d<sup>-1</sup></td><td>Calculated</td></tr><tr><td>7, 8</td><td>F<sub>l</sub></td><td>LMWC leaching loss</td><td></td><td></td><td></td><td></td><td></td><td>g C m<sup>-2</sup> d<sup>-1</sup></td><td>Calculated</td></tr></tbody></table>

(continued on next page)
<table><caption>**Table A1 (continued)**</caption><thead><tr><th>Equation Number</th><th>Variable</th><th>Definition</th><th>Default values</th><th>Fit Values: All</th><th>Fit Values: VR</th><th>Fit Values: KG</th><th>Fit Values: LUCAS</th><th>Units</th><th>Sources</th></tr></thead><tbody><tr><td>7, 9, 19</td><td>F<sub>lm</sub></td><td>Adsorption of LMWC to minerals</td><td></td><td></td><td></td><td></td><td></td><td>g C m<sup>-2</sup> d<sup>-1</sup></td><td>Calculated</td></tr><tr><td>7, 13, 20, 21, 22</td><td>F<sub>lb</sub></td><td>Uptake of LMWC by microbial biomass</td><td></td><td></td><td></td><td></td><td></td><td>g C m<sup>-2</sup> d<sup>-1</sup></td><td>Calculated</td></tr><tr><td>7, 16, 19, 20</td><td>F<sub>bm</sub></td><td>Microbial mortality</td><td></td><td></td><td></td><td></td><td></td><td>g C m<sup>-2</sup> d<sup>-1</sup></td><td>Calculated</td></tr><tr><td>12, 19</td><td>F<sub>ld</sub></td><td>Desorption</td><td></td><td></td><td></td><td></td><td></td><td>g C m<sup>-2</sup> d<sup>-1</sup></td><td>Calculated</td></tr><tr><td>22</td><td>F<sub>bg</sub></td><td>Microbial growth</td><td></td><td></td><td></td><td></td><td></td><td>g C m<sup>-2</sup> d<sup>-1</sup></td><td>Calculated</td></tr><tr><td>17, 18, 19</td><td>F<sub>ma</sub></td><td>Aggregate formation from MAOM</td><td></td><td></td><td></td><td></td><td></td><td>g C m<sup>-2</sup> d<sup>-1</sup></td><td>Calculated</td></tr><tr><td>20, 22, 23</td><td>F<sub>mr</sub></td><td>Microbial respiration</td><td></td><td></td><td></td><td></td><td></td><td>g C m<sup>-2</sup> d<sup>-1</sup></td><td>Calculated</td></tr><tr><td>2, 3</td><td>V<sub>pl</sub></td><td>Maximum rate of POM decomposition to LMWC</td><td></td><td></td><td></td><td></td><td></td><td>d<sup>-1</sup></td><td>Calculated</td></tr><tr><td>13, 15</td><td>S<sub>w,B</sub></td><td>Water scalar</td><td></td><td></td><td></td><td></td><td></td><td>-</td><td>Calculated</td></tr><tr><td>2, 4, 5, 6, 8, 9, 15, 18</td><td>S<sub>w,D</sub></td><td>Diffusion limitation of substrates</td><td></td><td></td><td></td><td></td><td></td><td>-</td><td>Calculated</td></tr><tr><td>9, 10</td><td>K<sub>lm</sub></td><td>LMWC and microbial necromass adsorption rate</td><td></td><td></td><td></td><td></td><td></td><td>d<sup>-1</sup></td><td>Calculated</td></tr><tr><td>9, 11</td><td>Q<sub>max</sub></td><td>Maximum sorption capacity</td><td></td><td></td><td></td><td></td><td></td><td>g C m<sup>-2</sup></td><td>Calculated</td></tr><tr><td>13, 14</td><td>V<sub>lb</sub></td><td>Potential LMWC uptake rate</td><td></td><td></td><td></td><td></td><td></td><td>d<sup>-1</sup></td><td>Calculated</td></tr><tr><td>1, 7</td><td>F<sub>i</sub></td><td>C input</td><td></td><td></td><td></td><td></td><td></td><td>g C m<sup>-2</sup> d<sup>-1</sup></td><td>Input</td></tr><tr><td>3, 16, 21, 22</td><td>T</td><td>Soil temperature</td><td></td><td></td><td></td><td></td><td></td><td>°C</td><td>Input</td></tr><tr><td>4, 15</td><td>θ</td><td>volumetric water content</td><td></td><td></td><td></td><td></td><td></td><td>mm<sup>3</sup> mm<sup>-3</sup></td><td>Input</td></tr><tr><td>23</td><td>CO<sub>2</sub></td><td>carbon dioxide production</td><td></td><td></td><td></td><td></td><td></td><td>g C m<sup>-2</sup> d<sup>-1</sup></td><td>Output</td></tr></tbody></table>

<caption>**Table A2**</caption>

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.

<table><thead><tr><th>Variable</th><th>Definition</th><th>Default values</th><th>Fit Values: All</th><th>Fit Values: VR</th><th>Fit Values: KG</th><th>Fit Values: LUCAS</th><th>Units</th><th>Sources</th></tr></thead><tbody><tr><td>w1</td><td>Water scalar parameter</td><td>30.0</td><td>33.61</td><td>17.38</td><td>36.97</td><td>22.61</td><td>-</td><td>Parton et al. (2010)</td></tr><tr><td>w2</td><td>Water scalar parameter</td><td>9.00</td><td>8.42</td><td>10.55</td><td>7.77</td><td>7.84</td><td>-</td><td>Parton et al. (2010)</td></tr><tr><td>t1</td><td>x-axis location of inflection point</td><td>15.4</td><td>18.08</td><td>16.83</td><td>15.77</td><td>18.15</td><td>°C</td><td>Del Grosso et al. (2005)</td></tr><tr><td>t2</td><td>y-axis location of inflection point</td><td>11.75</td><td>12.99</td><td>10.39</td><td>9.32</td><td>13.78</td><td>-</td><td>Del Grosso et al. (2005)</td></tr><tr><td>t3</td><td>Distance from the maximum point to the minimum point (step size)</td><td>29.7</td><td>25.92</td><td>30.10</td><td>20.98</td><td>28.91</td><td>-</td><td>Del Grosso et al. (2005)</td></tr><tr><td>t4</td><td>Slope of line at inflection point</td><td>0.031</td><td>0.038</td><td>0.025</td><td>0.032</td><td>0.035</td><td>-</td><td>Del Grosso et al. (2005)</td></tr><tr><td>c1</td><td>Intercept of clay fraction relationship</td><td>0.85</td><td>0.83</td><td>0.85</td><td>0.85</td><td>0.86</td><td>-</td><td>Parton et al. (1987)</td></tr><tr><td>c2</td><td>Slope of clay fraction relationship</td><td>0.68</td><td>0.67</td><td>0.71</td><td>0.74</td><td>0.60</td><td>-</td><td>Parton et al. (1987)</td></tr><tr><td>k<sub>ls</sub></td><td>Turnover rate of structural litter pool</td><td>0.01</td><td>-</td><td>-</td><td>-</td><td>-</td><td>d<sup>-1</sup></td><td>Parton et al. (1987)</td></tr><tr><td>k<sub>lm</sub></td><td>Turnover rate of metabolic litter pool</td><td>0.045</td><td>-</td><td>-</td><td>-</td><td>-</td><td>d<sup>-1</sup></td><td>Parton et al. (1987)</td></tr><tr><td>k<sub>a</sub></td><td>Turnover rate of active pool</td><td>0.020</td><td>-</td><td>-</td><td>-</td><td>-</td><td>d<sup>-1</sup></td><td>Parton et al. (1987)</td></tr><tr><td>k<sub>s</sub></td><td>Turnover rate of slow pool</td><td>5.0 × 10<sup>-4</sup></td><td>6.7 × 10<sup>-4</sup></td><td>1.9 × 10<sup>-4</sup></td><td>5.2 × 10<sup>-4</sup></td><td>4.0 × 10<sup>-4</sup></td><td>d<sup>-1</sup></td><td>Parton et al. (1987)</td></tr><tr><td>k<sub>p</sub></td><td>Turnover rate of passive pool</td><td>2.0 × 10<sup>-5</sup></td><td>9.1 × 10<sup>-6</sup></td><td>8.2 × 10<sup>-6</sup></td><td>1.9 × 10<sup>-5</sup></td><td>1.3 × 10<sup>-5</sup></td><td>d<sup>-1</sup></td><td>Parton et al. (1987)</td></tr><tr><td>p<sub>li</sub></td><td>Proportion of plant residue to structural litter pool</td><td>0.66</td><td>0.62</td><td>0.67</td><td>0.64</td><td>0.67</td><td>-</td><td>Analogous to Millennial litter partitioning assumption</td></tr><tr><td>p<sub>lma</sub></td><td>Fraction of metabolic litter to active pool</td><td>0.45</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>Parton et al. (1987)</td></tr><tr><td>p<sub>lsa</sub></td><td>Fraction of structural litter to active pool</td><td>0.5</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>Parton et al. (1987)</td></tr><tr><td>p<sub>lss</sub></td><td>Fraction of structural litter to slow pool</td><td>0.7</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>Parton et al. (1987)</td></tr><tr><td>p<sub>sa</sub></td><td>Fraction of slow pool to active pool</td><td>0.42</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>Parton et al. (1987)</td></tr><tr><td>p<sub>sp</sub></td><td>Fraction of slow pool to passive pool</td><td>0.03</td><td>0.020</td><td>0.029</td><td>0.022</td><td>0.041</td><td>-</td><td>Parton et al. (1987)</td></tr><tr><td>p<sub>pa</sub></td><td>Fraction of passive pool to active pool</td><td>0.45</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>Parton et al. (1987)</td></tr><tr><td>p<sub>ap</sub></td><td>Fraction of active pool to passive pool</td><td>0.004</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>Parton et al. (1987)</td></tr><tr><td>LigFrac</td><td>Fraction of litter that is lignin</td><td>0.2</td><td>0.20</td><td>0.20</td><td>0.21</td><td>0.20</td><td>-</td><td>Zhang et al. (2018)</td></tr><tr><td>%claysilt</td><td>Percent of soil that is in the clay and silt fractions</td><td>80</td><td>-</td><td>-</td><td>-</td><td>-</td><td>%</td><td>site-specific</td></tr><tr><td>F<sub>c</sub></td><td>Field capacity</td><td>0.39</td><td>-</td><td>-</td><td>-</td><td>-</td><td>mm<sup>3</sup></td><td>Input: observed maximum</td></tr><tr><td>F<sub>i</sub></td><td>C input</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>mm<sup>-3</sup></td><td>volumetric water content</td></tr><tr><td>T</td><td>Soil temperature</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>°C</td><td>Input</td></tr><tr><td>θ</td><td>volumetric water content</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>mm<sup>3</sup></td><td>Input</td></tr><tr><td>StrLitter</td><td>Structural litter pool</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>g C m<sup>-2</sup></td><td>Calculated</td></tr><tr><td>MetLitter</td><td>Metabolic litter pool</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>g C m<sup>-2</sup></td><td>Calculated</td></tr><tr><td>ACTIVE</td><td>Active pool</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>g C m<sup>-2</sup></td><td>Calculated</td></tr><tr><td>SLOW</td><td>Slow pool</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>g C m<sup>-2</sup></td><td>Calculated</td></tr><tr><td>PASSIVE</td><td>Passive pool</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>g C m<sup>-2</sup></td><td>Calculated</td></tr><tr><td>S<sub>t</sub></td><td>Temperature scalar</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>Calculated</td></tr><tr><td>S<sub>w</sub></td><td>Water scalar</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>Calculated</td></tr><tr><td>F<sub>tex</sub></td><td>Soil texture effect on decomposition of active pool</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>Calculated</td></tr><tr><td>F<sub>ls</sub></td><td>Structural litter decomposition flux</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>Calculated</td></tr></tbody></table>

(continued on next page)
<table><caption>**Table A2 (continued)**</caption><thead><tr><th>Variable</th><th>Definition</th><th>Default values</th><th>Fit Values: All</th><th>Fit Values: VR</th><th>Fit Values: KG</th><th>Fit Values: LUCAS</th><th>Units</th><th>Sources</th></tr></thead><tbody><tr><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td>g C m<sup>-2</sup><br>d<sup>-1</sup></td><td></td></tr><tr><td>F<sub>lm</sub></td><td>Metabolic litter decomposition flux</td><td></td><td></td><td></td><td></td><td></td><td>g C m<sup>-2</sup><br>d<sup>-1</sup></td><td>Calculated</td></tr><tr><td>F<sub>a</sub></td><td>Active pool decomposition flux</td><td></td><td></td><td></td><td></td><td></td><td>g C m<sup>-2</sup><br>d<sup>-1</sup></td><td>Calculated</td></tr><tr><td>F<sub>s</sub></td><td>Slow pool decomposition flux</td><td></td><td></td><td></td><td></td><td></td><td>g C m<sup>-2</sup><br>d<sup>-1</sup></td><td>Calculated</td></tr><tr><td>F<sub>p</sub></td><td>Passive pool decomposition flux</td><td></td><td></td><td></td><td></td><td></td><td>g C m<sup>-2</sup><br>d<sup>-1</sup></td><td>Calculated</td></tr></tbody></table>

## Appendix B

### Century model equations

Model inputs are identical in description and units to those used in the Millennial model (Table A1). All other parameters are described in Table A2, including their default and fitted values for each dataset.

$$S_t = \frac{t2 + \frac{t_3}{\pi} \operatorname{atan}[t4 * \pi (T - t1)]}{t2 + \frac{t_3}{\pi} \operatorname{atan}[t4 * \pi (30 - t1)]} \qquad [B1]$$

$$S_w = \frac{1}{1 + w1 * e^{-w^2 * \frac{a}{K}}} \qquad [B2]$$

$$F_{tex} = c1 - c2 \%claysilt 0.01 \qquad [B3]$$

$$F_{ls} = \text{StrLitter} k_{ls} S_t S_w e^{-3 \text{LigFrac}} \qquad [B4]$$

$$F_{lm} = \text{MetLitter} k_{lm} S_t S_w \qquad [B5]$$

$$F_a = \text{ACTIVE} k_a S_t S_w F_{tex} \qquad [B6]$$

$$F_s = \text{SLOW} k_s S_t S_w \qquad [B7]$$

$$F_p = \text{PASSIVE} k_p S_t S_w \qquad [B8]$$

$$\frac{d\text{StrLitter}}{dt} = p_{li}F_i - F_{ls} \qquad [B9]$$

$$\frac{d\text{MetLitter}}{dt} = (1 - p_{li})F_i - F_{lm} \qquad [B10]$$

$$\frac{d\text{ACTIVE}}{dt} = (1 - \text{LigFrac})p_{isa}F_{ls} + p_{ima}F_{lm} + F_s p_{sa} + F_p p_{pa} - F_a \qquad [B11]$$

$$\frac{d\text{SLOW}}{dt} = \text{LigFrac} p_{isa} F_{ls} + F_a (1 - F_{tex} - p_{ap}) - F_s \qquad [B12]$$

$$\frac{d\text{PASSIVE}}{dt} = F_a p_{ap} + F_s p_{sp} - F_p \qquad [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. https://doi.org/10.1007/s10533-017-0409-7.

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

* Allison, S.D., Wallenstein, M.D., Bradford, M a, 2010. Soil-carbon response to warming dependent on microbial physiology. Nature Geoscience 3, 336–340. https://doi.org/10.1038/ngeo846.

* Arrhenius, S., 1889. Über die Reaktionsgeschwindigkeit bei der Inversion von Rohrzucker durch Säuren. Zeitschrift für Physikalische Chemie 4, 226–248.

* Bailey, C.J., 1989. Enzyme kinetics of cellulose hydrolysis. Biochemical Journal 262, 1001–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.1111/gcb.13926.

* 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. https://doi.org/10.1071/SR12374.

* Beaudoing, H., Rodell, M., 2016. GLDAS Noah Land Surface Model L4 Monthly 1 X 1 Degree V2, p. 1.
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. https://doi.org/10.1016/j.gloplacha.2012.10.006.

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

Cotrugo, 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. https://doi.org/10.1038/s41561-019-0484-6.

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–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.1111/gcb.15381.

Del Grosso, S.J., Parton, W.J., Mosier, A.R., et al., 2005. Modeling soil CO2 emissions from ecosystems. *Biogeochemistry* 73, 71–91. https://doi.org/10.1007/s10533-004-0898-z.

Dove, N.C., Torn, M.S., Hart, S.C., Taş, 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., Siccamo, 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.

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*. https://doi.org/10.1038/s41467-017-01116-z.

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.

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. https://doi.org/10.1111/j.1365-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.

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.1111/geb.13186.

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. https://doi.org/10.5194/bg-10-2011-2013.

Hengl, T., Jesus, JM De, Huvelink, 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* 1313–139.

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–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. https://doi.org/10.1038/ncomms13630.

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. https://doi.org/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*. https://doi.org/10.1088/1748-9326/ab2ee4.

Luo, Y., Ahlström, 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.

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. https://doi.org/10.1016/j.envsoft.2019.104518.

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.1029/2004GB002219.

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. https://doi.org/10.2136/sssaj.

Michaelis, L., Menten, M.O., 1913. The kinetics of the inversion effect. *Biochemische Zeitschrift* 49, 333–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. https://doi.org/10.1016/j.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.

Panagos, P., Van Liedekerke, M., Jonesa, A., Montanarella, L., 2012. European Soil Data Centre: response to European policy support and public data requirements. Land Use Policy 29, 329–338. https://doi.org/10.1016/j.landusepol.2011.07.003.

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.

Parton, W.J., Schimel, D.S., Cole, C.V., et al., 1987. Analysis of factors controlling soil organic matter levels in Great plains grasslands. *Soil Science Society of America Journal* 51, 1173–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. https://doi.org/10.1016/j.soilbio.2018.06.025.

Poeplau, C., Kätterer, 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.1111/gcb.13491.

Ramírez, P.B., Calderón, 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.3389/fenvs.2021.634472.

Richter, D. deB., Hofmockel, M., Callaham, M.A., et al., 2007. Long-term soil experiments: keys to managing earth's rapidly changing carbon fractions. *Soil Science Society of America Journal* 71, 266–279. https://doi.org/10.2136/sssaj2006.0181.

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.5194/bg-2018-430.

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. https://doi.org/10.1016/S0304-3800(97)00067-7.

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.

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.1016/j.ecolmodel.2013.04.010.

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.1038/s41467-018-04526-9.

Sierra, C.A., 2012. Temperature sensitivity of organic matter decomposition in the Arrhenius equation: some theoretical considerations. *Biogeochemistry* 108, 1–15. https://doi.org/10.1007/s10533-011-9596-9.

Sierra, C., Muller, M., 2015. A general mathematical framework for representing soil organic matter dynamics. *Ecological Monographs* 16, 16881. https://doi.org/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. https://doi.org/10.1111/gcb.13556.

Søtaert, K., 2009. rootSolve: Nonlinear Root Finding, Equilibrium and Steady-State Analysis of Ordinary Differential Equations.

Søtaert, K., 2016. R Package FME: Inverse Modelling, Sensitivity, Monte Carlo-Applied to a Dynamic Simulation Model.

Søtaert, K., Herman, P.M.J., 2009. A Practical Guide to Ecological Modelling. Using R as a Simulation Platform. Springer.

Spann, 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, 1–22.

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.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. https://doi.org/10.1016/j.geoderma.2013.06.014.

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.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. https://doi.org/10.1016/j.soilbio.2019.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. https://doi.org/10.5194/gmd-8-3823-2015.

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
decomposition. Biogeosciences 10, 8329–8351. https://doi.org/10.5194/bg-10-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.1038/38260.

Tóth, 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. https://doi.org/10.1007/s10661-013-3109-3.

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.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.1111/ejs.12237.

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. https://doi.org/10.1071/SR14366.

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.

von Lützow, M., Kögel-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. https://doi.org/10.1016/j.soilbio.2007.03.007.

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.

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.1038/ismej.2014.120.

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. https://doi.org/10.1890/12-0681.1.

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–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. https://doi.org/10.1002/2015GB005188. 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. https://doi.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. https://doi.org/10.1111/gcb.13979.

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. https://doi.org/10.1029/2018JG004777.

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.

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. https://doi.org/10.1111/geb.12029.

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. https://doi.org/10.5194/gmd-11-4779-2018.