Microbial models with minimal mineral protection can explain long-term soil organic carbon persistence

Soil organic carbon (SOC) models currently in widespread use omit known microbial processes, and assume the existence of a SOC pool whose intrinsic properties confer persistence for centuries to millennia, despite evidence from priming and aggregate turnover that cast doubt on the existence of SOC with profound intrinsic stability. Here we show that by including microbial interactions in a SOC model, persistence can be explained as a feedback between substrate availability, mineral protection and microbial population size, without invoking an unproven pool that is intrinsically stable for centuries. The microbial SOC model based on this concept reproduces long-term data (r2 = 0.92; n = 90), global SOC distribution (rmse = 4.7 +/− 0.6 kg C m−2), and total global SOC in the top 0.3 m (822 Pg C) accurately. SOC dynamics based on a microbial feedback without stable pools are thus consistent with global SOC distribution. This has important implications for carbon management, suggesting that relatively fast cycling, rather than recalcitrant, SOC must form the primary target of efforts to build SOC stocks.


Calculation of rate modifying constants
With the exception of the microbial biomass rate modifier (µ), all rate modifiers are based on established models (RothC and CENTURY) in order that differences between SOMic and traditional 1st order models can be attributed to the microbial dynamics rather than to different soil temperature and moisture responses.

Temperature
The rate modifying factor for temperature (τ) is given by a generalised Poisson function, which is the temperature response function of the CENTURY model 8 multiplied by a normalising coefficient ( f t ) to give τ = 1.0 at the average mean annual temperature of the calibration sites (10.8 • C).  (3) in which T is the soil temperature ( • C), f t = 4.99, T max = 45 • C, and T opt = 35 • C.

Moisture
The rate modifying factor for soil moisture (θ) is derived in one of two ways, depending on whether soil moisture data are available. If soil moisture is not known, as is often the case in field experiments (including the results presented in the main manuscript for long-term agricultural experiments), it is estimated from precipitation and potential evapotranspiration (PET) according to the RothC algorithm 1 : where max_md is the maximum possible soil moisture deficit, and amd is the accumulated moisture deficit (see ref. 1 for further detail on definitions and how these are calculated).
If soil moisture data are available, as for example in the calculation of global SOC stocks using soil moisture and temperature values from the Community Land Model (Section 2.4), the calculation of θ is based on a modified version of the RothC Equation, which uses the ratio of soil moisture to field capacity in the place of AMD: where θ m is the soil moisture and θ f is the field capacity.

Soil cover
The soil cover rate modifying factor (ν) represents the faster mineralization of topsoil organic carbon following tillage that has removed vegetation. It is calculated according to the RothC method 1 as

Microbial biomass
SOC decomposition shows a non-linear response to enzyme concentration, which is believed to be due to competition for enzyme binding sites on substrates, as predicted by Langmuir adsorption isotherm theory 9 . This non-linear relationship can be approximated by a reverse Michaelis-Menten (MM) equation 2 , which, assuming the approximation that exoenzyme production is proportional to microbial biomass 10 , can be expressed as where, µ is the rate modifying factor for microbial biomass, µ max is its maximum (enzyme saturated) value, [MB] is the concentration of microbial biomass, and the MM constant (K M ) is the [MB] at which the reaction rate is at half-maximum. Whereas, in practice, different substrates could have different values of µ max and K M , the simplifying assumption was made in SOMic version 1.0 to use common values for all microbially-mediated processes to keep the number of model parameters small and thus avoid overfitting.

Partitioning of DOC between sorption and microbial uptake
Only DOC is assumed to be taken up by microbes, because compounds must be in solution to cross the cell membrane. Microbes must also compete for DOC with abiotic sorption of DOC to minerals which also removes carbon from the DOC pool. The aggregate amount of C removed in a given time period from the DOC pool (D doc ) by both microbial uptake and sorption is given by where k doc is the modified rate constant for the combined processes of microbial uptake and sorption. Note that, although the form of equation 8 appears similar to a first-order reaction, the dependence of the rate factor k doc on reverse MM dynamics and on competition between microbes and sorption (see equations 9-13 below) ensures that it is not a simple 1st order process. k doc is calculated as the weighted mean of the modified rate constants for sorption (k sorb ) and microbial uptake (k mu ) according to in which the weighting applied to the sorption coefficient is the fraction ( f sorb ) of D doc that is sorbed, and the weighting for the microbial uptake is the remainder (1 − f sorb ) of D doc that is taken up by microbes. The modified rate constant k mu is the product of the base rate constant for microbial uptake (k mu ) and the rate modifying factors itemised above in Section 1.4: The modified rate constant k sorb is calculated similarly, but without any dependence on the microbial biomass rate-modifier (µ), instead having a rate modifying factor for soil clay content. For the clay rate modifier (c), we used a linear function of clay content (Equation 11), as is common practice in most current models (although future models may be improved by considering also clay mineralogy). 11, 12 where, m c is the slope of the clay response function, and c 0 is the clay content at which c is equal to one. In order that c represents the change in rate constant relative to the average clay content in the calibration sites, c 0 was defined as the mean clay content at the calibration sites (23%). m c was then calculated by numerical optimisation as one of the calibration parameters (see Section 1.8).
The modified rate constant k sorb is then given by The fraction of carbon removed from the DOC pool due sorption can now be calculated from the modified sorption and microbial uptake rate constants as whereby the carbon sorbed (S) during a time period is and the microbial uptake (M) is

Partitioning of microbial uptake between growth and respiration
Microbial uptake (M) is partitioned between growth (G) and respiration (R) by the microbial carbon use efficiency (CUE).
The balance between microbial growth and respiration, and thereby CUE, is known to vary with environmental parameters, and in particular temperature [13][14][15][16][17] . The temperature dependence of CUE was modelled as a linear function according to where T is the soil temperature ( • C), CUE 0 is the CUE at 15 • C, and m cue is the rate of change of CUE with temperature. CUE 0 and m cue were estimated by model calibration as described in Section 1.8 below. 3/20

Cycling of carbon between pools
Using the rate-modifying factors and partition functions derived above, the evolution of the MB, MAC, and DOC pools and respired CO 2 -carbon over time is given by the following equations, in which the quantity of carbon in a pool in time period i is the carbon remaining from time i − 1 plus the new carbon added to that pool: where, k desor b is the base rate constant for desorption of MAC, and k mt is the base rate constant for microbial biomass turnover. The dependence of MAC desorption rate on reverse MM dynamics in equation 19 is in accordance with observations that desorption of organic matter from mineral surfaces is accelerated by microbial enzyme activity, which in turn is limited by enzyme activity site availability. [18][19][20][21] Mechanisms by which microbial activity can destabilise and increase desorption of MAC include, for example, the destabilisation of soil aggregates by decomposition of organic binding agents 22 ; the utilisation of sorbed organic compounds by microorganisms that adhere to the mineral surfaces [22][23][24] ; microbial transformation of sorbed organic compounds into more readily desorbed compounds 25 ; and displacement by microbial exudates 21 .

Differential form of model equations
The rate of change of carbon in each pool can also be written in differential form (equations 23-27 below) to make it easier to see structure of the carbon fluxes between pools. For brevity, the model carbon pools in Eq. 23-27 are denoted as C 1 (SPM), C 2 (IPM), C 3 (DOC), C 4 (MAC), and C 5 (MB). Note that, in these differential equations, k 1 , k 2 , k 3 , k 4 , and k 5 are the decomposition rate factors for the pools C 1 to C 5 , respectively, after the rate modifying factors have been applied. Therefore, although these differential equations appear similar in form to first order decay, the dependence of the rate factors on the size of other pools ensures that they are not, in fact, first order. In particular, the dependence of k 1 , k 2 , k 3 , k 4 , and k 5 on microbial biomass abundance through MM means that these factors are not constants but are themselves dynamic quantities.
where, dCO 2 dt is the rate of respired CO 2 evolution; and d L dt is the rate of fresh plant litter input to the soil.  27 . Crop rotation is a 2-year winter wheat/fallow system, with a 15-month fallow and a 9-month cropping season. The climate is semi-arid. Temperature, precipitation, and PET were calculated using the Food and Agriculture Organisation of the United Nations' (FAO's) New_LocClim climate-station database and interpolation software, by inverse distance weighted average (IDWA), with vertical and horizontal regression correction, from a maximum of 50 climate stations within a maximum of 1000 km (Table 2).

4/20
The nine treatments modeled here varied in terms of mineral nitrogen fertilizer additions, residue management (fall burning, spring burning, or no burning), and organic additions (none, manure, or pea vine), as described in Table 3.
Above-ground carbon additions to the soil are from Rasmussen and Parton 28 . The experiment was revised in 1967 to change the wheat type from a medium-tall to a semi-dwarf variety. Below-ground carbon allocation was estimated using a root:shoot ratio of 0.625 for the long-straw variety, and 0.5325 for the short straw 29 . This provided the estimated total carbon inputs in Table 4.
Soil carbon measurements from Rasmussen et al. 27 are given in Table 5.

Sanborn Field
Sanborn Field, is located on the campus of the University of Missouri-Columbia (39.03 • N, 94.58 • W), and was established in 1888. The soil is a Mexico silt loam with 28.3% mean clay content in the maize and wheat plots modelled 30 . Mean annual temperature is 13 • C, mean annual precipitation is 973 mm, and mean annual potential evapotranspiration is 790 mm 31 . Monthly mean values of these climate parameters are given in Table 6. The experiment station comprises rotation and manure treatments on 39 plots. The treatments modelled were continuous maize with full mineral fertiliser and regular tillage (cmf), continuous maize with no fertiliser of manure (cmn), continuous wheat with full mineral fertiliser (cwf), continuous wheat with manure (cwm), and continuous wheat with no fertiliser of manure (cwn). Total carbon additions to the soil from crop residues and below ground production are given in Table 3 of Buyanovsky and Wagner 31 . In addition, the manure treatments received 13.4 Mg manure annually. The Phyllis2 Database for biomass and waste 32 gives mean water content of farmyard manure as 70%, mean ash content 33 %, and mean dry ash-free carbon content as 48 %. On this basis, the annual carbon addition in manure was estimated at 1.3 Mg C y −1 . Details on the management of the treatments and soil carbon measurements can be found in Buyanovsky and Wagner 31 .

Rothamsted Experimental Station
The Rothamsted Experimental Farm (51 • 49' N, 0 • 21' W) is the oldest long-term agricultural experiment site globally, established in 1843. The soil is a silty clay loam (FAO classification is Chromic Luvisol) with an average 23.4% clay content 33 . Complete experimental data including annual management, yields, meteorology and soil carbon are available from the Rothamsted electronic archive 34 , with a summary of the modelled sites provided here. The mean annual temperature is 9.2 • C, mean annual precipitation is 704 mm, and mean annual potential evapotranspiration is 450 mm 1 . Monthly mean values are given in Table 7.
Broadbalk: The first experimental winter wheat crop was sown on Broadbalk in autumn 1843, and every year since then 35 . All plots received annual returns of dead crop material (roots, stubble and part of the chaff), and were fallow only between Aug harvest and Nov sowing. Treatments varied according to the amounts and types of mineral and organic fertilisers that were applied (Table 8). Wheat straw was assumed to be 85% dry matter 34 , and 48.79% C 32 . In the plots that received 35 Mg ha −1 y −1 farmyard manure, this was assumed to provide 3 Mg C ha −1 y −133 .
Hoosfield: The Hoosfield experiment was started in 1852. Spring barley grown every year except 1912, 1933, 1943 and 1967 when the whole experiment was fallowed to control weeds. In contrast to Broadbalk, because it is spring sown, it has only been necessary to fallow these four times to control weeds 34 . The annual input of plant residues from the barley was estimated as 2.80 Mg C ha −1 y −1 in treatments and years in which farmyard manure was applied, and 1.6 Mg C ha −1 y −1 in unmanured treatments and years 1 . The modelled treatments are described in Table 9. In the plots that received 35 Mg ha −1 y −1 farmyard manure, this was assumed to provide 3 Mg C ha −1 y −133 .

Global SOC distribution
For the estimate of global SOC stocks and distribution, Community Land Model (CLM) spatial output data were used to estimate soil temperature, moisture and litter inputs. CLM version 4.5CN data were from the historic land post-processed monthly-average data for the years 1850-2010 36 available on request from the National Center for Atmospheric Research (NCAR) at https://www.earthsystemgrid.org.
Soil temperature in the top 0.3 m was calculated as the weighted mean of the top 5 layers from the CLM output file clm45cn_1deg4502_hist.clm2.h0.TSOI.185001-201012.nc (weighted by layer thickness). Soil moisture was the weighted mean of the top five layers of clm45cn_1deg4502_hist.clm2.h0.H2OSOI.185001-201012.nc. Saturation capacity (θ s ; m 3 m −3 ) was estimated using the following equation

5/20
where, D b is bulk density (Mg m −3 ), and clay is in percent 37 . Global gridded bulk density and soil texture estimates were obtained from the regridded harmonized world soils database v 1.2 38 .
Global ecoregions were based on the classification from Olson and Dinerstein 39 , with the associated spatial data available at http://maps.tnc.org/files/shp/terr-ecoregions-TNC.zip.

Radiocarbon age of SOC
The age of SOC as a function of depth was calculated by running the model on a sequence of soil layers from the top horizon down to below 1 m depth, with DOC leaching represented by the quantity of DOC removed from each soil layer in each time step being added to the next lower layer. Because C is exchanged between DOC, MAC, and MB as it percolates through the soil column, the mean DOC age in a soil layer is always (i) younger than the mean SOC age in that layer (because it contains some younger DOC added from fresh litter or from the layer above), and (ii) older than the fresh litter or DOC inputs, because it contains some C derived from desorbed MAC or MB turnover. This means that the DOC input into subsequent layers increases incrementally in age (the mean age of DOC in a layer is always older than the fresh C inputs into that layer), in a process that has been described as "cycling downwards" 25 .
The soil layer thicknesses were based on the CLM depth definitions and are given in Table 10. In each time period, the exchange of C between the DOC, MAC, and MB pools is first calculated according to the Equation 22. Then DOC advection from one layer to the next was calculated using an advection velocity estimated using CLM values for groundwater recharge rate for the years 1850-2010 36 . Time averaged water velocity for the Rothamsted site varies with depth from a mean surface infiltration rate of 310 mm yr −1 to an aquifer recharge rate of 45mm yr −1 -a value that averages saturated and unsaturated flow conditions over time. Daily values can be obtained from the CLM model output file clm45cn_1deg4502_hist.clm2.h0.QCHARGE. 185001-201012.nc available on request from the National Center for Atmospheric Research (NCAR). Soil temperature and moisture in each layer were also taken from the same CLM run, in files clm45cn_1deg4502_hist.clm2.h0.TSOI. 185001-201012.nc and clm45cn_1deg4502_hist.clm2.h0.H2OSOI.185001-201012.nc, respectively. To estimate the radiocarbon age profile at a specific location, the vertically layered model with DOC advection was spun up over 12,000 years using the litter input values and gridded soil characteristics as described in Section 2.4. For example, using the CLM forcing data at Rothamsted (51 • 49' N, 0 • 21' W) give predicted versus measured 40 radiocarbon ages with an adjusted R 2 = 0.67, p < 2.2x10-16, n=78) (Fig. 1). Note that the 78 data points from ref. 40 excluded measurements from the Park Grass surface soils, which had anomalously old 14 C ages, due to contamination with coal.   1931-1966 1967-1978 1979-1986 1931-1966 1967-1978 1979-1986 1931-1966 1967-1978 1979-1986 fB_N0  fB  fB  fB  0  0  0  0  0  0  sB_N0  sB  sB  sB  0  0  0  0  0  0  nB_N0  nB  nB  nB  0  0  0  0  0  0  sB_N45  fD  nB  sB  0  0  0  0  45  45  sB_N90  sD  nB  sB  0  0  0  0  90  90  nB_N45  fD  nB  nB  0  0  0  34  45  45  nB_N90  sD  nB  nB  0  0  0  34  90  90  nB_PV  nB  nB Table 4. Total organic carbon inputs to the soil at Pendleton OR, including above-ground residues, below-ground net primary production, and organic amendments (Mg C ha −1 crop −1 ). 1931-1941 1942-1951 1952-1966 1967-1976 1977-1986