Long residence times of rapidly decomposable soil organic matter : application of a multi-phase , multi-component , and vertically resolved model ( BAMS 1 ) to soil carbon dynamics

Accurate representation of soil organic matter (SOM) dynamics in Earth system models is critical for future climate prediction, yet large uncertainties exist regarding how, and to what extent, the suite of proposed relevant mechanisms should be included. To investigate how various mechanisms interact to influence SOM storage and dynamics, we developed an SOM reaction network integrated in a one-dimensional, multi-phase, and multi-component reactive transport solver. The model includes representations of bacterial and fungal activity, multiple archetypal polymeric and monomeric carbon substrate groups, aqueous chemistry, aqueous advection and diffusion, gaseous diffusion, and adsorption (and protection) and desorption from the soil mineral phase. The model predictions reasonably matched observed depth-resolved SOM and dissolved organic matter (DOM) stocks and fluxes, lignin content, and fungi to aerobic bacteria ratios. We performed a suite of sensitivity analyses under equilibrium and dynamic conditions to examine the role of dynamic sorption, microbial assimilation rates, and carbon inputs. To our knowledge, observations do not exist to fully test such a complicated model structure or to test the hypotheses used to explain observations of substantial storage of very old SOM below the rooting depth. Nevertheless, we demonstrated that a reasonable combination of sorption parameters, microbial biomass and necromass dynamics, and advective transport can match observations without resorting to an arbitrary depth-dependent decline in SOM turnover rates, as is often done. We conclude that, contrary to assertions derived from existing turnover time based model formulations, observed carbon content and 114C vertical profiles are consistent with a representation of SOM consisting of carbon compounds with relatively fast reaction rates, vertical aqueous transport, and dynamic protection on mineral surfaces.


Introduction
Soil organic matter (SOM) represents a large stock of carbon that can be exchanged with the atmosphere over short (seconds to weeks; Trumbore, 2000), intermediate (months to annual; Baldocchi et al., 2001), and long time frames (decades to centuries; Baisden et al., 2002;Ciais et al., 2012;Hsieh, 1993).Current estimates are that more than three times as much carbon is stored in terrestrial soils (2344 Pg C) than is in the atmosphere (Jobbagy and Jackson, 2000), although recent estimates for high-latitude systems indicate that value for soil carbon may be an underestimate (Schuur et al., 2009;Tarnocai et al., 2009).Soil temperature and moisture are expected to change as climate changes over the 21st century and are important controllers of the rate at which soil carbon is retained or released to the atmosphere (Parton et al., 1987; Published by Copernicus Publications on behalf of the European Geosciences Union.Todd-Brown et al., 2013).Therefore, SOM dynamics could generate important feedbacks to climate change.However, recent reviews have indicated that current site and global land biogeochemical (BGC) models do not represent the suite of mechanisms necessary for prediction of future SOM stocks (Conant et al., 2011;Dungait et al., 2012;Schmidt et al., 2011), and that their SOM predictions have large uncertainties (Friedlingstein et al., 2006;Todd-Brown et al., 2013).
At a particular place in the soil, temporal variation in organic carbon in its various phases may depend on decomposing surface litter, root mortality and exudation, dissolved organic matter (DOM) and inorganic carbon diffusion and advection, colloidal transport, gross carbon fluxes adsorbing and desorbing from soil mineral surfaces, biomineral aggregate formation and destruction, carbon assimilation and release from microbes, transport into plants, and others (Dungait et al., 2012).Each of these fluxes has specific, and mostly poorly characterized, controlling mechanisms with different environmental sensitivities, and predicting their net impact at any particular location in the soil is therefore difficult.Relevant environmental controllers on these fluxes include soil temperature and moisture, pH, redox state, alternative electron acceptor, chemical and physical inhibitors, soil structure, and others (e.g., Conant et al., 2011;Schimel and Schaeffer, 2012).These controllers are also affected by the wide range of relevant temporal scales (minutes to millennia) and spatial heterogeneity, making prediction of the carbon feedbacks between soil and atmosphere even more difficult.
Although the mechanisms impacting SOM stocks mentioned above have been hypothesized, and sometimes demonstrated, to impact SOM dynamics, studying their effects using observations or experiments within the complex soil environment is difficult, since synergistic and competing processes may mask the effect of any particular mechanism.A complimentary approach is to develop mechanistic numerical models that specifically represent processes hypothesized to be relevant.However, biogeochemical models always require some process simplification, and those simplifications may not represent all relevant interactions.We discuss some examples of this conundrum below, but note that the level to which mechanistic detail needs to, or can, be included in land models remains unclear.One goal of the current work is to begin to develop a modeling structure where hypotheses regarding model complexity can be addressed.
Given the importance of SOM dynamics in affecting atmospheric greenhouse gas (i.e., CO 2 , N 2 O, and CH 4 ) concentrations and soil nutrient availability for plant growth (Denman et al., 2007), developing mechanistic and reliable numerical models that can be integrated in general circulation or Earth system models (GCMs or ESMs, respectively) is crucial.However, the current suite of land models integrated in ESMs represent at most only a few of the mechanisms described above, making it impossible to analyze the impact of individual or interacting mechanisms using this suite of models.Current ESM belowground biogeochemistry submodels range from very simple, single-pool, turnover-time models to those that broadly follow the single soil layer structure of the Century Model or its progeny DayCent and ForCent (Parton et al., 1987(Parton et al., , 1998(Parton et al., , 2010)).Recently, a few models applied, or designed to be applied, at climate-model resolution have considered vertically resolved SOM processes (e.g., Braakhekke et al., 2011;Jenkinson et al., 2008;Koven et al., 2013;Tang et al., 2013) imposed on a Century-like carbon pool structure.However, none of these models explicitly represent the processes discussed above (such as adsorption and protection, desorption, and microbial activity).
The Century (Parton et al., 1987(Parton et al., , 1998(Parton et al., , 2010) ) and RothC models (Jenkinson and Coleman, 2008;Jenkinson et al., 2008) are archetypes of the most common belowground carbon biogeochemistry models.This type of model represents carbon moving from aboveground and root litter into several carbon pools, with the distribution of inputs determined by a metric of litter quality (e.g., lignin content, C : N ratio).SOM decomposition is predicted from pre-determined carbon pool turnover times that are modified based on temperature, soil moisture, and clay content.These models do not explicitly represent microbial activity, although several recent publications have highlighted its importance in terrestrial BGC models (Singh et al., 2010;Tang and Riley, 2013;Wang et al., 2013;Wieder et al., 2013;C. G. Xu et al., 2011).A growing body of evidence also suggests that soil mineralogy plays a large role in determining SOM stocks and turnover times, especially below the active root zone (Masiello, 2004;Masiello et al., 2004;Mikutta et al., 2006;Torn et al., 1997).The inclusion of soil texture or clay content in some models has been proposed as a proxy for organo-mineral interactions.The field studies cited above document large differences in soil carbon stock and turnover times that could be explained by mineralogy but not by texture or clay content.
There are a group of models that have not applied the intrinsic turnover time concepts underlying Century, and these models do include some of the important processes described above.Perhaps the most robust example is the ecosys model (Grant et al., 2003(Grant et al., , 2012)), which explicitly represents multiple microbial functional groups, coupled reaction kinetics constrained by oxidation-reduction energy yields, transport, impacts of plants, spatial heterogeneity, and interactions with the atmosphere.These models form an important proof of concept for the type of model developed here, which attempts to resolve complex simultaneous reaction-diffusionadvection networks, yet has a design compatible with integration into a global climate model.
It has been suggested that the next generation of land BGC models to be used for climate prediction should include, explicitly or implicitly, representation of vertical carbon transport, microbial activity, mineral surface interactions, temperature and moisture controls on independent processes, and nutrient dynamics (Schmidt et al., 2011).However, because of the complex nature of the expected interactions when this complex group of mechanisms is integrated into a reactive transport framework, the high computational cost of performing these types of simulations, and the issues associated with vertical and horizontal heterogeneity, a promising first step would be to develop and apply one-dimensional mechanistic modeling frameworks that can be applied at the site scale.Model reduction techniques (e.g., Riley, 2013), parameter inversion (e.g., (Luo et al., 2009;Vrugt et al., 2008)), and spatial scaling methods will likely be required in translating these approaches into representations appropriate for ESMs.
Toward this end, we develop here a model of SOM dynamics (called BAMS1; Biotic and Abiotic Model of SOM version 1) by integrating into a subsurface reactive transport model (TOUGHREACT; Gu et al., 2009;Maggi et al., 2008;Xu, 2008;Xu et al., 2011) a microbially mediated reaction network with two functional groups (heterotrophic aerobic bacteria and fungi), dynamic adsorption and desorption, multiple SOM species with varying parameters representing different aspects of reactivity, and mineral surface exchanges.To make this first attempt at developing a tractable model structure, we have omitted several of the mechanisms described above, including soil aggregation and its potential for protecting SOM, colloidal transport, and nutrient interactions.After describing the model structure, we present comparisons against observations from a compilation of grassland sites in several western US states.We then present a series of analyses designed to test the relative importance of the mechanisms described above in affecting long-term SOM dynamics.

Reactive transport solver
The numerical code we used to solve the BAMS1 reaction network and transport is TOUGHREACT (Xu, 2008), a three-dimensional reactive transport simulator based on the TOUGH model (Pruess et al., 1999) into which an arbitrary reaction network can be integrated.The model can represent a large suite of processes, including aqueous, gaseous, and sorbed phases; advection and diffusion; multiple competing microbial populations; adsorption and desorption; hydrology; the soil energy budget; and equilibrium and nonequilibrium chemical reactions.TOUGHREACT has been applied to analyze many aspects of soil systems, including near-surface carbon and N cycling (e.g., Gu et al., 2009;Gu and Riley, 2010;Maggi et al., 2008;Spycher et al., 2009;Xu et al., 2004).Details of the numerical methods and process representations used in TOUGHREACT can be found in the model's technical guide (Xu et al., 2013).
For BAMS1, we applied a one-dimensional version of TOUGHREACT, which solves simultaneous mass balance equations for aqueous and gaseous phase tracer concentrations (C i (mol-C L −1 )): where t (s) is time, z is depth (m), D (m 2 s −1 ) is the effective aqueous or gaseous diffusivity, and v (m s −1 ) is bulk aqueous or gaseous velocity.The first and second right hand side terms of Eq. ( 1) account for diffusive and advective transport and the third right hand side term accounts for an arbitrary number (m) of sources, sinks, and exchanges between phases.

Carbon decomposition reaction network
For this first version of the model, we assumed the active decomposers in soils consist of heterotrophic aerobic bacteria and fungi.We chose these two microbial groups because they are known to have specific affinities to decompose plant litter and other SOM compounds (DeAngelis et al., 2013;Neely et al., 1991;Romani et al., 2006;Thevenot et al., 2010) and to produce different necromass (Frostegard and Baath, 1996).
Although including only two functional groups of microbes substantially under-represents observed functional diversity in soils (Goldfarb et al., 2011), several factors influenced our choice: (1) observations sufficient to develop models of multiple SOM decomposing microbial functional groups are not available, nor has there been an attempt to synthesize observations in a manner amenable for inclusion in the type of model applied here (although see Bouskill et al. (2012) for an example of a trait-based approach applied to nitrification) and (2) our goal is to make a first attempt at analyzing interactions between the processes hypothesized to be important for SOM dynamics, not to fully constrain every possible complexity in the system.The numerical relationships and baseline parameters (and sensitivity analyses around these baseline values) needed to simulate microbial dynamics are discussed in the next section.Soil organic matter consists of thousands of plant and microbe synthesized and degraded compounds, and numerical models cannot represent this complexity in its entirety.Further, partitioning these compounds into groups according to properties relevant to SOM dynamics is difficult because the many relevant processes are controlled by different compound traits.For example, the set of rules used to group compounds into classes of intrinsic decomposition rates may be different than the classification needed to characterize sorptive interactions with mineral surfaces.Other compound classifications of interest for SOM cycling would reflect the influence of the compound chemistry on carbonuse efficiency (Brant et al., 2006;Frey et al., 2013;Sinsabaugh et al., 2013), high-affinity interactions with soil minerals (Gordon and Millero, 1985;Gu et al., 1994;Kleber et al., 2005;Mikutta et al., 2007), susceptibility to exoenzyme decomposition (Allison et al., 2010;Schimel and Weintraub, 2003;Sinsabaugh and Shah, 2012) temperature sensitivity (Conant et al., 2011;Davidson and Janssens, 2006).The impacts of these mechanisms depend strongly on the environment in which they operate.Carbon use efficiency, for example, depends on the metabolic disposition of the decomposer biomass, and varies depending on whether the microbial community at a given point in time is disposed to grow (to incorporate C) or in need of a critical nutrient that requires an energy consuming, dissimilatory metabolic mode.It is thus not always defensible to characterize a given substrate as having a certain carbon use efficiency, as suggested by some authors (e.g., Frey et al., 2013).
In this first attempt at integrating some of these processes into a reactive transport model, we implemented an approach that is more mechanistic while capturing variations in substrate chemistry.
Our approach is to group compounds by quantifiable properties that can be considered to be relevant for metabolic processing, i.e., oxygen : carbon (O / C) ratio, positive or negative charge, and degree of polarity (Table 1).The rationale behind this strategy is that organic compounds that are already oxidized (as expressed by the O / C ratio) will require less oxygen and lower energy inputs to be fully oxidized to CO 2 than more reduced compounds.Polar compounds will be more soluble and thus more mobile in a largely aqueous soil system than non-polar compounds.Ionization state and the sign of an electrical charge are determinants for protective sorptive and other interactions with both the mineral matrix and charged organic colloids.Our approach is designed to be flexible enough to incorporate other hypotheses regarding controls on decomposition and transport such as accessibility and aggregation.
The model's decomposition reaction network consists of four groups of organic compounds: (1) above and belowground litter and root exudates; (2) organic polymers; (3) simpler organic monomers; and (4) fungal and bacterial biomass (Fig. 1).Carbon is brought into the system through aboveground litter, belowground litter, and root exudates.Through implicit (i.e., not explicitly represented as a separate pool) exoenzyme activity, the litter can be degraded into several simpler organics (i.e., monosaccharides, amino sugars, organic acids, amino acids, phenols, lipids, and nucleotides) in proportions that can be manipulated in the model for sensitivity analyses (baseline values shown in Fig. 2).Root exudates are considered to be simple monomeric organics.
The model assumes fungi and bacteria can assimilate the polymeric and monomeric organics directly, resulting in biomass increases and CO 2 production.The proportions of compounds that can be assimilated by either microbial population are set as baseline values (Fig. 2) and can be manipulated in sensitivity analyses.Carbon in microbial necromass is returned to solution (Fig. 1) and can participate in further reactions, be transported vertically, or interact with soil grain surfaces.We parameterized a fraction (f m ) of the microbial necromass as peptidoglycan, which is intended to represent polymeric cell wall material (Fig. 1).The remaining necromass carbon is returned to the monomer pools, in proportions that can be manipulated for a particular simulation (Fig. 2).Microbial exudation or carbon overflow (Tempest and Neijssel, 1992) from stoichiometric imbalance can be an important contributor to overall microbial carbon turnover in some cases.However, we have not included them in the current version of the model.

Decomposition reactions of organic monomers
We assume the decomposition of monomers (monosaccharide, amino acid, amino sugar, organic acid, lipid, and nucleotide) is carried out by aerobic heterotrophic bacteria, which utilize monomers as a source of energy and carbon for biomass growth and maintenance, thus leading to biomass yield and CO 2 production.The reaction describing decomposition of 1 mole of generic monomer k C i , with k carbon atoms can be written as follows: where B A (mg C-wet-biomass L solution −1 ) is the aerobic heterotrophic bacteria biomass concentration, a is the assimilation-to-respiration ratio describing the fraction of carbon atoms being assimilated into the cell as compared to the fraction released as CO 2 , and R O,i is the oxidative ratio of the compound, which describes the contribution of aqueous O 2 uptake as compared to the O 2 available from the substrate.The biomass yield coefficient Y i = ak for B A describes the microbial biomass gain per mole of consumed substrate ((g wet-biomass) (mol substrate) −1 ).The effect of temperature on the decomposition rate was not included in Eq. (2).We leave the determination of an appropriate parameterization for the temperature effect and related analyses to future work.

Stoichiometry for fungal depolymerization of organic compounds
Large polymers (represented in this model as cellulose (C 6 HO 5 ) n , hemicellulose (C 9 H 15 O 7.5 ) n , and lignin (C 10 H 12 O 3 ) n ) can be decomposed (depolymerized) into monomers by fungi.The fungi may then assimilate part of the carbon from the substrate polymer, consume O 2 , produce energy and CO 2 , and release free monomers available to other microorganisms such as aerobic bacteria.The specific monomer products vary with the source polymer (Fig. 2).A general depolymerization reaction of polymer m j C i , with k carbon atoms made of n monomers m j C j with m j carbon atoms, is derived from stoichiometric constraints: where B F (mg wet-biomass L −1 ) is the fungal biomass concentration, W is the depolymerization efficiency describing  the fraction of carbon released in the monomers, and x j are the stoichiometric coefficients in Eq. ( 2) for each monomer.In Eq. ( 3), the biomass yield coefficient Y i is Note that Eq. ( 3) can also describe depolymerization reactions by exoenzymes (i.e., no direct assimilation and biomass growth) when W = 1 (i.e., Y i = 0).

Biological reaction rates
The microbial decomposition rate of substrate m j C i is generically described using Michaelis-Menten kinetics (e.g., Maggi et al., 2008): where is the rate of change of C i for decomposition (mol m −3 s −1 ); subscript d represents decomposition; µ i (s −1 ) is the maximum specific consumption rate of substrate i (s −1 ) (i.e., the intrinsic decomposition rate); O 2 (mol O 2 m −3 ) is aqueous O 2 concentration; B (mg wetbiomass L −1 ) is biomass of bacteria (B A ) or fungi (B F ); Y i is the biomass yield on substrate i; and K i and K O 2 (mol L −1 ) are half-saturation coefficients for C i and O 2 , respectively (Table 2).The functions f (θ ), g (pH), and h (T ) describe the environmental effects on microbial activity of soil moisture, pH, and temperature T , respectively (Maggi et al., 2008).In the current suite of experiments, we did not include the effects of pH or soil water stresses (i.e., g (pH) and f (θ ) are set to 1) on decomposition rates, although these factors can be important (e.g., Schimel et al., 2011).We note that recent work has shown that Michaelis-Menten kinetics may be inaccurate when simulating complex consumer-substrate interactions (Tang and Riley, 2013), but leave incorporation of these effects into the model for later work.O 2 is consumed during aerobic microbial respiration on monomers and fungal depolymerization (Fig. 1) with rate: CO 2 is produced in both reactions with rates that depend on the carbon content in the substrate and the reaction structure.The reaction rates for monomer oxidation and depolymerization are respectively, where W is the depolymerization efficiency, k is the number of carbon atoms in C i , n is the number of monomers in the polymer, and m j is the number of carbon atoms in the product monomer C j .
The decomposer biomass grows as substrate is assimilated, and declines as cells die, according to the following: where δ (s −1 ) is microbial death rate.
We followed the hypothesis of Kleber ( 2010) that the O / C ratio (R O/C,i ) of relatively simple organic compounds provides a proxy for the maximum specific consumption rate (Table 2).We fit a simple exponential curve to the observations shown in Fig. 2 of Kleber (2010), and normalized the result to the rate associated with monosaccharides: For peptidoglycan, Eq. ( 10) was not used but a relatively slow intrinsic rate and higher adsorption rate were imposed to represent the tendency of these larger molecules to decompose more slowly (Miltner et al., 2012).

Abiotic processes
CO 2 gas dissolution is calculated assuming local equilibrium between aqueous and gaseous phases (Maggi et al., 2008).Diffusivities are estimated based on the molecular weight of each component and the Millington and Quirk (1961) approach to account for tortuosity.Water flow is modeled using the Darcy-Richards equation (Pruess et al., 1999) and the van Genuchten (1980) relationships for matric potential and hydraulic conductivity with soil water saturation.Advective gaseous transport as a result of pressure gradients was not considered in the present study, although it can be included using the equations of states available in the TOUGHREACT standard features.
Adsorption and desorption are complex processes that depend on characteristics of the organic molecules, soil mineral surfaces and properties, and aqueous chemistry (Dudal and Gerard, 2004).Further, there are large differences in temperature sensitivities for the various sorption mechanisms (Conant et al., 2011), motivating explicit representation of these processes in models.To examine the impact of sorption on predicted SOM stocks, we imposed forward (adsorption; k f (s −1 )) and reverse (desorption; k r (s −1 )) rates.In the absence of competing sources and sinks of a particular species, this formulation would result in an effective equilibrium linear sorption relationship: K d = k f /k r .In our formulation, sorption reactions are subsumed in the m ∂C i ∂t m terms of Eq. (1), k f is taken to be 6.6 × 10 −8 s −1 , and sorbed species are protected from decomposition.
Linear isotherms have been commonly applied in soil sorption studies (e.g., Neff and Asner, 2001).However, other studies (Gu et al., 1994;Kothawala et al., 2008Kothawala et al., , 2009;;Mayes et al., 2012) concluded that DOM sorption in soils more closely follows a Langmuir relationship than a linear relationship.We note that both of these approaches assume equilibrium between the sorbed and aqueous phases, although the Langmuir relationship implies a partitioning dependent on soil properties and the sorbed concentration.We intend to examine the impact of the functional form of the sorption isotherm on SOM dynamics in future work.Soil aggregation (Six et al., 2001(Six et al., , 2004) is another important process that can lead to carbon persistence in soil.The density of aggregates and rates of aggregation and disaggregation in soil depend on many ecosystem properties (Conant et al., 2011), so properly representing their impacts on protected SOM may be important in climate change studies.Although important for SOM residence times, the density of aggregates and rates of aggregation and disaggregation have not yet been represented in the model.

Climate forcing, boundary conditions, and initial conditions
Climate forcing (precipitation, evapotranspiration, temperature) and carbon inputs used in the BAMS1 simulations were taken from a CLM4 (Lawrence et al., 2011) simulation of a US Great Plains grassland.CLM4 is the land model integrated in the Earth system model called CESM1 (http: //www.cesm.ucar.edu/models/cesm1.0/), and represents surface and subsurface hydrology, ecosystem energy exchanges with the atmosphere, and plant and soil biogeochemistry.
The equilibrium simulations for comparison to SOM observations were forced with repeated 1948-1972 cycles from Qian et al. (2006).For the BAMS1 simulations, we imposed the CLM4 predicted infiltration and carbon inputs from leaf, wood, and root litter partitioned into the soil using exponentially decaying depth profiles with length scales of 1, 7, and 12 cm, respectively.Baseline simulations were run for 10 000 years from initial conditions of no soil carbon and for constant climate forcing.This spin-up period was sufficient to ensure that overall changes in SOM stocks were less than 0.01 % yr −1 by the end of the simulation.
Having a fully coupled model that includes reactive transport capabilities and the aboveground process representations in CLM4 would be the preferred model structure to perform the analyses we are presenting here because it would allow for the coupling between above and belowground processes.However, to our knowledge such a capability does not currently exist (although efforts are underway, e.g., Tang et al., 2013).

Comparison to SOM and 14 C observations
In this work we were primarily interested in exploring hypotheses that could explain variations in SOM residence times; therefore, our comparisons to SOM observations were designed to engender confidence that the suite of processes we have included were internally consistent and produced realistic SOM predictions, but not that they describe precisely the conditions and SOM profiles of a specific site.We therefore compared the predicted SOM profiles using the baseline set of parameters to observations from grasslands located between 39 • N and 43 • N latitude across Nebraska and Colorado.Compiled vertically resolved observations of SOM content from the National Soil Carbon Database (NSCN; available at http://www.fluxdata.org/nscn/SitePages/Home.aspx) included 618 observations in that region; Fig. 3).Along with SOM profiles, root profiles from grasslands in the same region were selected from the Global Root Distribution Profile database (Schenk and Jackson, 2005).
As another check on model simulations, we compiled observations of the fraction of total biomass that are fungi and bacteria.To this end, a database of microbial soil vertical profiles from existing literature was compiled with data from Alvarez et al. (1995), Anderson and Domsch (1989), Ekelund et al. (2001), Fierer et al. (2003), Zhou et al. (2008), Selvam et al. (2010), Tsai et al. (2007), Agnelli et al. (2004), and Yeates et al. (1997).SOM radiocarbon ( 14 C) profiles are potentially valuable tracers of the integrated suite of processes affecting vertical SOM profiles and have recently been applied to evaluate climate-scale land surface model SOM predictions (Jenkinson and Coleman, 2008;Koven et al., 2013).We predicted SOM 14 C values in the year 2003 by running a parallel simulation to the bulk SOM simulation, but imposing a first order decay term (with turnover time of 8267 year corresponding to the 14 C radioactive decay rate) for each modeled component in both aqueous and protected phases.To impose the "bomb" atmospheric 14 C content, we applied the northern hemispheric values from Levin and Hesshaimer (2000) for 1950-1976and Levin and Kromer (2004) for 1977-2003.

Model analyses
We performed ten sensitivity tests (Table 3) and a number of model experiments to explore how system properties and changes in soil microclimate may impact SOM dynamics.To investigate the relative importance of adsorption and desorption rates from soil mineral surfaces on equilibrium SOM stocks, we performed simulation scenarios varying the sorption rates k f and k r (Table 3).First, we modified k f and k r concurrently by factors of 0.8 (S1) and 1.2 (S2) (resulting in an unchanged value of K d ).Next, we investigated the impacts of the different temperature sensitivities of desorption and adsorption.As noted in Conant et al. (2011), warmer temperatures favor desorption for high-affinity soil OM-mineral interactions (e.g., non-covalent bonds) and favor adsorption for low-affinity soil OM-mineral interactions.For sensitivity analyses S3 and S4, we assumed a mean temperature increase throughout the profile that led to a relative increase in desorption (k r ) of a factor of 1.2 (S3; high affinity) and a relative increase in adsorption (k f ) of a factor of 1.2 (S4; low affinity).Although this sensitivity analysis substantially under-represents the complexity of the temperature dependencies of sorption and protection mechanisms, it does allow us to qualitatively investigate the impact of temperaturedependent sorption on vertically resolved SOM content.We note that protection associated with occlusion and aggregation will also have temperature sensitivities (Conant et al., 2011); we did not analyze those affects here.
We also manipulated the microbially mediated transformation rates by multiplying µ i by factors of 0.8 (S5) and 1.2 (S6).Because interactions between sorption and the fraction of each compound accessible to decomposition (i.e., in the aqueous phase) will impact SOM content, we performed perturbations to µ i in combination with perturbations to k f and k r (scenarios S7, S8, S9, and S10).Total SOM profiles for these scenarios were again compared to the baseline scenarios.For all the sensitivity analyses we also estimated the 14 C SOM vertical profile.
To further examine carbon and microbial dynamics in the system, we performed an experiment with two 500-year simulations, both starting from the end of the 10 000-year simulation.The first simply continued the 10 000-year simulation for an additional 500 years.In the second simulation, we doubled all chemical species initial concentrations in the top 20 cm from those at the end of the 10 000-year simulation, and performed a 500-year simulation.Differences between these simulations over the 500-year period were used to investigate coupled C, microbial, and transport dynamics.
Finally, we performed pulse carbon input experiments by doubling the steady-state concentration of all compounds in seven depth intervals (0-10, 10-20, 20-30, 30-40, 40-75, 75-125, and 125-200 cm) independently and tracking carbon flows throughout the column.These simulations were performed to analyze the effective turnover and transit times of the various organic carbon inputs as they move between pools and with depth in the soil column.From the 10 000-year spin-up initialization state, the model was run for 5000 years, and comparisons were made to the baseline simulation, which was also run for 5000-year post initialization state.Anomalies of total SOM were calculated against the baseline simulation, and normalized to the amount of carbon injected for each experiment.Although the responses of individual compounds and of total SOM are complex and do not generally follow first-order dynamics, we calculated a simple effective first-order turnover time for total SOM and total DOM as a function of depth interval by treating them as single pools with constant turnover times.We also discuss the complexity of interpreting turnover times in a complex reaction network such as the one proposed here.

Equilibrium predictions for the grassland sites
Using the baseline set of parameters (Table 2), climate, carbon inputs, and a 10 000-year simulation, the model estimated total SOM contents roughly in agreement with the grassland observations we compiled (Fig. 3).At steady state, most of the input polymers (cellulose, hemicellulose, lignin) were predicted to be in the top meter of soil (Fig. 4a), consistent with the depth-distribution of the source carbon and the assumed inability of protected-phase compounds to move vertically.Lignin had the largest concentration of the input polymers, with a peak at about 10 cm depth.The fraction of total SOM predicted to be lignin was highest (∼ 1-2 %, as compared to more than 20 % for leaf and root C inputs) in the zone where direct carbon inputs occurred and decreased to ∼ 0 % by 1 m depth (Fig. 5b).The mean predicted proportion was about 1 % between 0 and 50 cm depth and compared well with the values from the synthesis by Thevenot et al. (2010), who reported a range for grasslands of ∼ 0.7-2.7 %.Cellulose and hemicellulose concentrations had flatter vertical concentration gradients and levels about an order of magnitude lower than that of lignin.Monomers protected (sorbed) on surfaces (Fig. 4b) had a deeper concentration profile than the input polymers, consistent with their ability to partition into the aqueous phase and move vertically with water advection and through diffusion.Of the sorption-protected monomers, the predicted relative concentration ranking, from highest to lowest, was monosaccharides, amino acids, lipids, phenols, and organic acids.Very little amino sugars or nucleotides were predicted to persist in the protected phase.
The vertical profile shape of the predicted total dissolved monomer content was similar to that for protected-phase monomers, while the relative abundances were 2-3 orders of magnitude lower (Fig. 4c).Dissolved monosaccharides had the highest predicted concentrations, followed by amino acids, lipids, phenols, and organic acids.The other dissolved monomers were predicted to have much lower concentrations.
The proportion of total microbial biomass predicted to be fungal was largest at the surface and declined rapidly to ∼ 0 % at 40 cm depth (Fig. 5a), consistent with the expectation that fungi inhabit the portion of the soil column dominated by recent plant carbon inputs and the few observations available for comparison (Ekelund et al., 2001;Zhou et al., 2008).

Sensitivity analysis for total SOM
To better understand the interactions among the represented controls on SOM dynamics, we analyzed sensitivity to various model parameters (S1-S10; Fig. 6; Table 3).Scenarios S1-S4 examined changes to sorption parameters only.First, we modified k f and k r concurrently by factors of 0.8 (S1) and 1.2 (S2) (resulting in unchanged partitioning coefficients compared to the baseline simulations).These perturbations, which were applied throughout the soil column, had small impacts on the SOM profile or total stocks (Fig. 6a; Table 3).We then investigated the impacts of asymmetric changes in k r and k f consistent with a temperature perturbation leading to  a factor of 1.2 increase in desorption rates (S3) and factor of 1.2 increase in adsorption rates (S4) (Fig. 6a).These changes also had very small impacts on predicted SOM content.
Modifying the bacterial SOM consumption rates (µ i ) on all aqueous compounds by factors of 0.8 (S5) and 1.2 (S6) had relatively larger impacts on SOM content (Fig. 6b).As expected, decreasing µ i led to increased SOM contents, while increasing µ i led to lower SOM contents throughout the soil column.As discussed below, these changes highlight the potential impacts of a temperature-dependent microbial growth rate.
Finally, modifying k f , k r , and µ i concurrently also had large impacts on SOM (Fig. 6c).Decreasing k f and µ i simultaneously (S7) led to increases in both 0-50 cm and 50-200 cm SOM content that were comparable (but slightly smaller) to those from S5 (which only decreased µ i ).Decreasing k f and increasing µ i simultaneously (S8) led to about a 30-50 % reduction in SOM, consistent with the impact of increasing µ i alone and the very small impact of decreasing k f alone.Decreasing k r and µ i simultaneously (S9) resulted in large increases in SOM throughout the column.Finally, decreasing k r and increasing µ i (S10), which alone have opposite impacts, led to a moderate decrease in SOM content.In these sensitivity scenarios, the impacts of changes in microbial growth rate dominated the impacts from changes in sorption parameters.

Sensitivity analysis for 14 C of total SOM
The baseline simulation resulted in 14 C profiles in soils that had a qualitative shape (i.e., monotonic depletion with depth) and values at depth consistent with observations (Fig. 6d), although the near-surface values were too depleted compared to the few observational data sets published for grasslands.We discuss potential reasons for this mismatch below, but note here that we did not attempt to tune a separate non-equilibrium pool near the soil surface that would allow for a better match with this expectation.We expected the 14 C value at about 1 m depth to be in the range [−600 ‰, −400 ‰], based on pasture and grassland observations at Paragominas, Brazil (Trumbore et al., 1995) and Riverbank, California and Judgeford, New Zealand (Baisden and Parfitt, 2007).
For scenarios S1, S2, S3, and S4, the predicted changes in the 14 C value of total SOM at 1 m depth were between ∼ 0 and −80 ‰.Using a simple one-box donor-controlled turnover model (i.e., first order loss) implies that the 14 Cinferred changes in turnover times between, e.g., S1 and the baseline simulations at 50 cm depth, were relatively larger than those predicted for the total SOM content changes (33 % versus 7 %, respectively).This pattern was consistent across scenarios S1, S2, S3, and S4 for all depths.Decreasing the microbial growth rate (µ i ; S5) and increasing µ i (S6) while holding all other parameters constant led to enhanced depletions of about 50 ‰ in 14 C below about 0.5 m depth (Fig. 6e).Scenarios S7-S10 all led to enhanced 14 C depletion compared to the baseline (Fig. 6f), with the largest enhanced depletion (∼ 100 ‰) corresponding to S9 (decreasing k r and decreasing µ i ) and S10 (decreasing k r and increasing  3).
µ i ) and the smallest (∼ 30 ‰) to S7 (decreasing k f and µ i ) and S8 (decreasing k f and increasing µ i ).

Transient pulse experiments
We performed several transient pulse experiments to investigate SOM turnover dynamics, vertical transport, and effective turnover times of the various compounds in the system.First, we applied a pulse of all compounds by doubling the initial concentrations in the top 20 cm of the soil column and performing a 500-year simulation starting with conditions from the end of the 10 000-year baseline simulation.Concentration differences compared to a 500-year simulation with baseline parameters (i.e., effectively continuing the 10 000year spin-up simulation another 500 years) varied widely between the various simulated carbon compounds (Fig. 7).The anomalies are shown as log 10 transformations of the absolute % differences from the baseline 500-year simulation; negative anomalies are indicated by white contours.Anomalies in cellulose and hemicellulose were 0.01-0.1 % in the top 5 cm of soil and persisted for the entire simulation.Negative anomalies (i.e., priming) were predicted between about 10 and 20 cm out to about 300 years.Lignin anomalies in the top 5 cm of soil had similar patterns to cellulose, but in contrast the positive perturbation of about 1 % below about 5 cm persisted for decades.Both aerobic bacteria and fungi were predicted to have ∼ 10 % increases in the top 20 cm of soil shortly after the perturbation, and these increases were predicted to persist for many decades.It appears that the microbial biomass has moved into another relatively persistent state.
The largest anomalies in protected monomers were predicted in the amino acids, amino sugars, phenols, and lipids pools, all of which showed perturbations of about 10 % in the top 20 cm of soil for hundreds of years.Beyond the polymer species, only the monosaccharides and organic acids were predicted to have negative anomalies, and the magnitudes of these anomalies were relatively small (less than ∼ 1 %).Anomalies in total SOM were predicted to be ∼ 1 %, 0.1 %, and 0.01 % in the 0-10, 10-100, and 100-200 cm depth ranges, respectively, by the end of the 500-year simulation.
We also performed a series of pulse experiments by doubling the initial total SOM in explicit depth intervals (0-10, 10-20, 20-30, 30-40, 40-75, 75-125, and 125-200 cm) to characterize the dynamic response of the system as it reequilibrates.We analyzed the predicted perturbed concentrations over 5000 years as scaled anomalies from the baseline simulation by normalizing total SOM and total DOM responses to the predictions immediately following the pulse input (Fig. 8a and b).For all depth intervals there were large predicted decreases in the SOM and DOM anomalies (∼ 30-60 % and ∼ 40-90 %, respectively)  To highlight important features, we: (1) used different ordinate axes for cellulose, hemicellulose, and lignin; (2) excluded anomalies less than 0.001 %; and (3) used solid filled contours where the anomalies were positive and white contour lines surrounding regions where the anomalies were negative (e.g., for the monosaccharides below about 20 cm after about 200 years).
years associated primarily with rapid microbial responses and advective transport.We applied a first-order approximation to estimate the turnover time (τ ) for total SOM and DOM anomalies to re-equilibrate as a function of the depth interval of the perturbation (Fig. 8c).
The SOM and DOM turnover times followed the same qualitative pattern with depth: a slight increase from an intermediate value down to 30 cm, an increase to about 1 m, and then a relatively constant τ below 1 m.However, the DOM τ was always lower than the total SOM τ , with the ratio between them being ∼ 100 down to ∼ 70 cm and a ratio of ∼ 3 below ∼ 70 cm.Turnover times estimated from the radiocarbon profiles show modest coherence with the pulse-based estimates for total SOM, but much higher turnover times for DOM dynamics.

Model predictions
Our model reproduced observed profiles of SOM 14 C values by representing several mechanisms that empirical studies have postulated to be important: sorption, advection, microbial dynamics (fungal and bacteria groups), and carbon compound specific properties.For example, the vertical DOM advection predicted by the model was consistent with that inferred in observational (e.g., Rumpel et al., 2002;Sanderman and Amundson, 2008) and other modeling analyses (Braakhekke et al., 2013).The model predicts that there are significant interactions between SOM and DOM beyond those expected from simple non-reactive transport of DOM from surficial layers (Sanderman et al., 2008), and that these interactions are responsible for the predicted older SOM and DOM 14 C values (Fig. 6) and long turnover times inferred from discrete pulse experiments and (c) Inferred first-order turnover times and 14 C ages for total SOM and total DOM as a function of depth interval following the perturbation.
the 14 C values (Fig. 8c).The prediction that total DOM concentrations were three orders of magnitude lower than protected-phase SOM concentrations (Fig. 4) is consistent with observations of a vertisol grassland (Don and Schulze, 2008), mollisol grassland, and ultisol forest sites (Sanderman et al., 2008).However, DOM concentrations are much more variable than SOM, so this ratio is not in general a good metric for comparison with model predictions.Combining all the dissolved monomer concentrations with the water flux through the system resulted in a predicted annual DOM leaching flux of 6.2 g m −2 year −1 , which matches very well the estimate by Kindler et al. (2011) for grasslands of 5.3 ± 2 g m −2 year −1 .
The role of microbial cell wall components as stable carbon compounds in soils has been discussed in several recent articles (e.g., Clemente et al., 2012).In the reaction network developed here, the compound we assigned to represent cell wall material (peptidoglycan) accounted for 5-10 % of SOM at steady state; this is a large enough fraction of SOM that resolving uncertainty in the partitioning of new microbial necromass into different species would be very useful.
Different model structures have reproduced observed decreases in SOM content with depth in a variety of ecosystems (Jenkinson et al., 2008;Koven et al., 2011Koven et al., , 2013;;Tang et al., 2013).However, because there are many ways that multiplepool SOM models can be calibrated against bulk SOM observations, we do not take the match of our vertical profile SOM content predictions with observations (Fig. 3) to be a particularly strong validation of our model structure.Since so many different model structures and parameterizations can match bulk SOM observations, the development of clear falsifiable hypotheses for these observed patterns remains an important goal for SOM model development.
Vertical profiles of 14 C values of total and dissolved organic carbon are a valuable constraint on process representation in SOM dynamics models.Although we did not simulate a particular site with 14 C observations, our predictions are consistent with the general structure of a monotonic decrease in 14 C value with depth to values less than −400 ‰ by ∼ 1 m depth.The decreasing 14 C values with depth occur because of a decrease in the decomposition rate with depth associated with decreasing substrate concentrations and therefore microbial biomass, rather than transport times from the surface, as also seen in sensitivity analyses in Koven et al. (2013).Over time, the accumulation of material on the mineral surfaces and 14 C decay leads to a relatively 14 C depleted SOM profile at depth.SOM pools were more 14 C depleted in the top ∼ 30 cm of the soil column than most observations suggest.We did not try to tune or optimize an arbitrary vertically variable mechanism to account for this deficit, as has been done for discrepancies with 14 C observations in other modeling exercises (e.g., Jenkinson and Coleman, 2008;Koven et al., 2013), because our main goal was to explore the effects of different mechanisms rather than to match a particular data set.We leave it to future work to propose and test specific alternative mechanisms that might lead to a larger proportion of younger carbon in the top portion of the soil column.However, for context, we estimated that an increase of ∼ 30 % in SOM concentrations resulting from plant inputs over the past several decades would lead to a close agreement with observations.Thus, an additional young and protected carbon pool of modest size (Fig. 3), and effectively not in equilibrium with the aqueous phase, can explain the bias between our predicted and commonly observed 14 C values near the soil surface.Our model allows for an additional non-equilibrium carbon pool that could be tuned to match these 14 C and SOM profiles, but we have avoided that type of tuning here.Processes that may be good candidates for this level of protection include aggregation and formation of colloids, which have been shown to substantially affect chemical mobility and carbon decomposition rates in soils (Daynes et al., 2013;Kausch and Pallud, 2013;Six et al., 2000).Development of process representations that will improve these comparisons is a task we leave for future research.
Other mechanisms not currently included in our model have also been hypothesized, and in some cases demonstrated, to play a role in the dynamics of deep SOM and its relatively depleted 14 C values.First, because the relative abundance of substrates at depth is low, physical separation between microbes and substrates could limit decomposition rates (Chabbi et al., 2009;Ekschmitt et al., 2008).The inaccessibility of substrates is expected to increase nonlinearly with soil drying.Models for this particular mechanism have been proposed by considering the effective diffusivity of substrate material as a function of pore geometry and moisture content (Or et al., 2007).Second, the density of mineral sorption sites increases with depth, which can lead to an effective inhibition of exoenzyme activity as these enzymes interact with mineral surfaces instead of substrate (Quiquampoix et al., 2002;Tietjen and Wetzel, 2003).Third, on the long time scales that we have simulated here, bioturbation may be an important vertical mixing mechanism, although it has been shown to be negligible in a forest soil (Braakhekke et al., 2013).As a final example, formation of aggregates that effectively shield SOM from decomposition could also play a role in some soils (see the review by Six et al., 2004).
Characterizing model structural uncertainty requires a more complete set of observational metrics that are specific to the model structure being tested.For a Century-like model without observationally defined pools, even with carefully constrained carbon inputs, one can only compare modeled and observed total SOM content and 14 C values.But, since specific combinations of SOM content and 14 C values in these types of models can be achieved with various combinations of parameters (Jenkinson and Coleman, 2008;Koven et al., 2013), and since the carbon content and 14 C value of modeled individual pools cannot be measured separately, it is difficult to tell if the partitioning of carbon between pools, vertical transport of material, and protection mechanisms were properly represented.This particular problem is partly alleviated in a model like the one described here; the problem then becomes finding observations of the specific carbon compounds and their partitioning between phases.More generally, modeling mechanisms or pools that can be approximated by measured quantities should enable much more informative testing and model development.
An additional technique to assess model structural uncertainty would be that of building a set of increasingly simplified models based on the one presented and analyzed here, and then comparing their goodness of fits against existing metrics after parameter estimation.We did not apply such an approach here since we did not have sufficient observations from one specific site, but could only assess model performance against trends from heterogeneous observations.However, testing model structural uncertainty will be performed in future work to gain confidence in the model structure and to inform the level of complexity required for a particular application.

Sensitivity analyses
Our sensitivity analyses were designed to better understand system dynamics and to highlight areas of uncertainty that could motivate future observations.We found that uncertainty in sorption rates (S3, S4; Fig. 6), designed to mimic sorption sensitivity to warming, had only a small effect on total-column SOM after 10 000 years.Of course, SOM dynamics under decadal-scale climate change will be influenced by many different temperature dependencies and their interactions.
Microbial growth rates had asymmetric impacts on SOM profiles.In particular, increased microbial growth rates (S6) decreased SOM content to a lesser degree than decreased growth rates (S5) increased SOM content.These scenarios both led to more 14 C depleted total SOM compared to the baseline, but with a larger decline associated with increased growth rates.Even though the high-growth scenario (S6) had less stabilized SOM at depth than scenario S5, the increased microbial growth rates led to about a doubling of peptidoglycan (i.e., microbial cell wall residuals after death) levels compared to the baseline.These compounds were strongly adsorbed and had low reactivity, thus they had relatively depleted 14 C values compared to the other protected species and resulted in 14 C depletion of the total SOM pool.
In the experiments manipulating sorption rates (k f , k r ) and microbial growth rate (µ i ), all four combinations led to SOM contents mostly outside the standard deviations of the observations.The largest biases were associated with both desorption and intrinsic microbial growth rates being simultaneously decreased (S9).Both scenarios that increased the microbial growth rate (S8, S10) led to low SOM content and all four scenarios led to more depleted SOM 14 C values as compared to the baseline simulation.

Pulse experiments
The pulse injection experiments demonstrated the temporal and vertical complexity of SOM responses.For example, some quantities peaked shortly after the addition of substrate and then declined and many had multiple increases and decreases in stock over long periods of time (Fig. 7), indicating Geosci.Model Dev., 7, 1335-1355, 2014 www.geosci-model-dev.net/7/1335/2014/ a relatively long memory in this dynamic system.In a real system with continuously or seasonally pulsing inputs, such a signal would be indistinguishable after a few carbon-input pulses.Some of the chemical species' had negative anomalies from the baseline that only manifested several decades after the perturbation.Although these anomalies were relatively small compared to the perturbation, they show that the overall system responses can have interesting transients and other complexities.Critically, the intrinsic rates affecting carbon turnover that are integrated in the model are faster than the overall system response characteristic time.Thus, the prediction that the total SOM characteristic response time is longer than these intrinsic rates indicates a substantial amount of recycling of carbon species between the aqueous, microbial, and protected pools.
The pulse injection experiment, where SOM was injected and followed for 5000 years, also showed relatively longer total SOM (∼ 100-300 years) and shorter DOM (0.5-10 years) simplified first-order turnover times down to about 30 cm (Fig. 8c).The decrease, and then increase, in turnover times with depth for this transient case is different than that inferred from our predictions of a monotonically decreasing 14 C profile of DOM and total SOM at steady state (Figs.6d  and 8c).This discrepancy highlights the differences in interpretation of SOM turnover rates based on the metric used to characterize turnover.Further, there appear to be at least two characteristic turnover times each for DOM and SOM in the pulse experiments: ∼ 1 year (Fig. 8a) and > 100 years (Fig. 8b).

Development of new SOM models
The model structures introduced by the Century (Parton et al., 1987) and RothC (Jenkinson and Coleman, 2008;Jenkinson et al., 2008) models include several SOM pools with pre-assigned turnover times that are modified by soil temperature, soil moisture, and soil texture.In that class of models, the input path of carbon to the various SOM pools may depend on the input lignin content, soil texture, or other ecosystem properties.Those models depend on the wide range in imposed SOM turnover times (months to centuries) to make equilibrium predictions consistent with observations, yet it is not possible to explicitly measure or identify the groups of compounds existing in any of the assumed pools.Further, several important aspects of that model structure have been challenged based on recent analyses using fine-resolution visualization (e.g., Lehman et al., 2008), molecular characterization (e.g., Kleber et al., 2010), field studies (e.g., Torn et al., 1997), and isotopic probes (e.g., Gleixner, 2013).
Because existing models do not represent many of the ecosystem processes that have strong climate dependencies, predictions of global terrestrial carbon cycle responses to climate change over the next century may be substantially in error.This contention is supported by recent analyses of Earth system model predictions (Todd-Brown et al., 2013).
We believe that the (1) recognition that important climaterelevant processes are missing from land models; (2) importance of terrestrial carbon dynamics on atmospheric greenhouse gas levels and climate; and (3) poor performance of these ESMs with respect to SOM dynamics, all argue for a re-evaluation of the methods used to predict carbon dynamics in land models.
Over the past decade, there have been major changes in the conceptual framework of SOM dynamics, evolving from one of carbon pools defined by assumed characteristic turnover times, to one in which organic material is cycled among pools of different physical-chemical state (such as adsorbed to a mineral surface or occluded in an aggregate) with turnover times primarily determined by the interaction of microbial and physical-chemical factors.This move away from a reliance on intrinsic recalcitrance to explain dynamics follows a parallel new understanding that soil organic matter comprises only a small amount of selectively preserved plant inputs and is composed more of microbial necromass and relatively simpler organic molecules with higher intrinsic decomposability (Miltner et al., 2012).The observed longevity, quantity, and vertical distribution of SOM profiles are inconsistent with a high level of decomposability, and are considered to be the result of various protection mechanisms (e.g., mineral interactions, aggregation), physical separation of substrates and their consumers, microbial population dynamics and activity (Allison et al., 2010;Lawrence et al., 2009;Tang and Riley, 2013;Wang et al., 2013), and transport mechanisms (Sanderman and Amundson, 2008).Since these mechanisms have different soil temperature and moisture sensitivities, have different characteristic temporal and spatial scales, and affect the various components of SOM differently, characterizing the impact of any single mechanism on emergent, vertically resolved SOM dynamics is difficult without a model capable of representing this coupled system.
The model structure proposed here is a start toward developing that type of process representation.Many other processes could also be included for a robust assessment of the emergent system behavior of soils under changing climate, changing vegetation type and status, and spatially varying soils, hydrology, and vegetation properties.

Observations needed to improve SOM models
The model structure described in this paper includes several of the processes and concepts recently identified as important in soil organic matter dynamics, including transport; mineral surface association and protection; explicit representation of fungi and bacteria; and explicit representation of input polymers, monomers, and recycling of microbial necromass.While improvements in representing these processes are needed, for global change predictions it would also be valuable to integrate and test the importance of mechanistic representation of at least seven additional processes: (a) aggregation (Segoli et al., 2013;Six et al., 2001) transport of colloids (Flury and Qiu, 2008;Thompson et al., 2006); (c) surface interactions (Conant et al., 2011); (d) enzyme dynamics (Allison et al., 2010); (e) nutrient-microbe interactions; (f) microbe-plant interactions; and (g) representation of subgrid-scale heterogeneity in soil properties and climate.
We propose that a new modeling methodology is needed that uses observations to improve explicit representation of the relevant processes while characterizing uncertainty and avoiding parameter equifinality (i.e., multiple parameter value-combinations resulting in a comparable match to the emergent-scale observations; Tang and Zhuang, 2008).In our study, moderate parameter perturbations caused large changes in predicted profile SOM content (Fig. 6), highlighting the value of coherent observational data sets that can be used to better constrain model parameter estimation and test model predictions (Moorhead et al., 2013).Although we argue for a more complex model structure to represent the many processes affecting SOM dynamics, the resulting system complexity makes determining individual parameter values difficult, especially if only observations of the emergent properties and behaviors of the system are used.For example, it has been common practice to use surface CO 2 flux measurements to infer temperature sensitivity of heterotrophic respiration.However, a careful analysis of the vertically resolved temperature and CO 2 production indicates that such an approach may be overly simplistic (Gu et al., 2004;Phillips et al., 2013;Tang and Riley, 2013).
For parameter calibration in a complex SOM model, it may be useful to use observations that target a specific mechanism under controlled conditions.For example, laboratory sorption experiments could be carried out in sterilized soil with known aqueous composition and temperature.Controlling observations in this manner reduces the influence of other processes (e.g., plant inputs, microbial decomposition, bulk transport), but we also recognize that laboratory experiments have artifacts and rates tend to be much faster in the laboratory than in situ.Uncertainty in the parameter values should be propagated through the remainder of the full model (Luo et al., 2009) and tested under relatively controlled field conditions.With careful calibration and uncertainty characterization, confidence in the coupled model predictions is enhanced or, if the total-system uncertainty is large, submodels that significantly reduce predictability could be identified for further work.Finally, such an uncertainty quantification and parameter estimation framework could also inform the level of complexity required to answer a particular question.

Application to regional-and climate-scale models
As with almost all soil biogeochemistry models, we have only applied the model introduced here in a one-dimensional column, but terrestrial ecosystems have spatial heterogeneity that affect SOM dynamics in edaphic properties (Bird et al., 2002), soil moisture (Riley and Shen, 2014;Sivapalan, 2005) and temperature (Davidson andJanssens, 2006), vegetation (Turner et al., 2004), and so on.While spatial variability in SOM dynamics occurs at scales from pore (King et al., 2010;Molins et al., 2012) to meters (Frei et al., 2012;King et al., 2010;Mishra and Riley, 2012) to km (Li et al., 2008), in climate models (typical resolution ∼ 100 km) the spatial heterogeneity in ecosystem properties is typically represented by spatially non-specific tiling of, e.g., plant functional types (Lawrence et al., 2011).In the absence of alternatives, it is often taken as an article of faith that singlecolumn representations of SOM dynamics can be spatially scaled by simple area averaging, but whether this approach gives a reasonable approximation to existing or future SOM stocks is not well tested.Addressing this question will require a suite of modeling approaches combined with finely resolved observations.In the modeling context, we contend that it is important to represent the underlying mechanisms of SOM dynamics at the spatial resolution at which they can be observed, at least in benchmark or exploratory models.It is unlikely, however, that a robust three-dimensional reactive transport solver will be developed in the near-term that can operate regionally or globally yet also simulate the very long time scales (> 1 K year) and fine spatial heterogeneity (µm -10 km) known to affect SOM dynamics.Thus, two promising strategies in the near term are to investigate soil system behavior on relatively smaller domains commensurate with the scale and availability of observations and to apply model reduction techniques (e.g., Olson et al., 2012;Robinson et al., 2012;Riley, 2013) to the full model.

Conclusions
The field of soil biogeochemistry faces pressing questions about how the many mechanisms known to be important for SOM cycling interact, which processes are most important for explaining current patterns, and which will be important for predicting future dynamics and transient responses.The model presented here takes a step towards addressing these questions by including sufficient complexity to explore competing hypotheses for explaining observed patterns in SOM stocks and turnover times.The model includes representations of two microbial functional groups, an equilibrium protection mechanism (sorption), and vertical transport.With these mechanisms, relatively rapid microbial transformation rates yielded SOM turnover times up to several thousand years at depth because of the transformation of plant material to microbial necromass that tended to be protected on mineral surfaces and low concentrations of dissolved, assimilable substrates.In other words, the persistence of theoretically rapidly decomposable substrates and increase in SOM turnover times with depth could be explained by microbial activity and transformation, sorption kinetics, and, to a lesser extent, vertical transport.The model also predicted complex transient responses to a perturbation in plant inputs, which would be difficult to determine experimentally.By showing that simple extrapolations to the future may not be accurate, these results confirm the need for models that are testable yet mechanistic enough to explore system responses to novel conditions such as those under climate change.Performing more sophisticated sensitivity analyses, parameter inversions, and perhaps developing reduced order models should allow a determination of the trade offs between increasing model complexity, parameter uncertainty, and model structural uncertainty.

Figure 1 .
Figure 1.Reaction Network Diagram.FUN and AER stand for fungi and heterotrophic aerobic bacteria, respectively.

Figure 2 .
Figure 2. Partitioning of the various carbon source pools into the eleven SOM destination pools.The values in each column sum to 1.0; each value indicates the fraction of the source pool allocated to a particular destination pool.

Figure 3 .
Figure 3. Predicted and observed SOM content.Observed values are represented as average and standard deviation of 618 observations from grasslands extracted from the NSCN database.

Figure 4 .
Figure 4. Predicted steady state (a) protected-phase polymers; (b) protected-phase monomers; and (c) dissolved monomers as a function of depth after 10 k year of simulations in the grassland.

Figure 5 .
Figure 5. Proportion of biomass that is fungal (a) and proportion of total SOM that is lignin (b) as functions of depth.

Figure 6 .
Figure6.Total SOM content (top row) and 14 C of total SOM (bottom row) at end of 10 000-year simulation for the ten sensitivity scenarios (a and d: S1-S4; b and e: S5-S6; c and f: S7-S10) that varied combinations of k f , k r , and µ i (Table3).

Figure 7 .
Figure7.Time and depth profiles of microbial and chemical species concentration anomalies (log10( % |anomaly|) compared to the baseline simulation) for the perturbation experiment that doubled all chemical species concentrations in the top 20 cm of soil as an initial condition.To highlight important features, we: (1) used different ordinate axes for cellulose, hemicellulose, and lignin; (2) excluded anomalies less than 0.001 %; and (3) used solid filled contours where the anomalies were positive and white contour lines surrounding regions where the anomalies were negative (e.g., for the monosaccharides below about 20 cm after about 200 years).

Figure 8 .
Figure 8. (a, b) Predicted 5000 year time histories of the scaled total SOM and total DOM anomalies resulting from a doubling of all compounds in the specified depth intervals (a shows first 3 years on a linear time scale and b shows 3-5000 years on a log time scale).(c)Inferred first-order turnover times and 14 C ages for total SOM and total DOM as a function of depth interval following the perturbation.

Table 1 .
Simplified overview of possible molecular functionalities (often dependent on solution pH, monomers only), number of repeating units for polymers, O / C ratios (R (O/C) ), oxidative ratio (R O ), and the sorption partitioning coefficient (k f k r ) for relevant compound classes.Y indicates that compound class possesses this particular option to interact with other molecules or surfaces (YY = possesses the option in extraordinary fashion).

Table 2 .
Parameters defining decomposition of organic matter pools by bacteria and fungi for maximum specific consumption rate (µ), assimilation-to-respiration ratio (α), and yield (Y).

Table 3 .
Description of parameters and mechanisms used in the sensitivity analysis.The three multipliers indicate the factors applied to the parameter for each sensitivity scenario.The last two columns indicate the impacts of each sensitivity scenario parameter changes on 0-50 cm and 50-200 cm total SOM contents.