Abiotic and Biotic Controls on Soil Organo–Mineral Interactions: Developing Model Structures to Analyze Why Soil Organic Matter Persists

Soil organic matter (SOM) represents the single largest actively cycling reservoir of terrestrial organic carbon, accounting for more than three times as much carbon as that present in the atmosphere or terrestrial vegetation (Schmidt et al. 2011; Lehmann and Kleber 2015). SOM is vulnerable to decomposition to either CO₂ or CH₄, which can increase atmospheric greenhouse gas concentrations (GHGs) and serve as a positive feedback to climate change. Conversely, the formation and stabilization of SOM within aggregates or associated with soil minerals can lead to carbon sequestration, representing a negative feedback to climate change. However, the conundrum as to why some SOM decomposes rapidly, while other thermodynamically unstable SOM can persist on centennial time scales (Hedges et al. 2000), leads to substantial uncertainty in model structures, as well as uncertainty in the predictability of the land carbon sink trajectory.

This chapter aims to tackle a part of this problem through generating recommendations for improved model structure in the representation of SOM cycling within global Earth System Models. The chapter is organized as follows: We first review the recent work contributing to the understanding of SOM stability. Then, we discuss current model structures and how these model structures deal with various processes, including data-model integration, microbial modeling, organo–mineral interactions, and SOM persistence. Finally, we conclude with recommendations for the future development of soil carbon models and the required relevant processes that determine the stability of SOM.

Traditional conceptual models of SOM formation and stability have focused on the chemical structure of SOM compounds. These models consider the recalcitrance of certain SOM compounds to biological decomposition as the primary determinant of persistence and long-term stability (Schmidt et al. 2011; Lehmann and Kleber 2015). Within this perspective, plant-derived polymers are the primary source of organic matter (Fig. 1a). These polymers were thought to form through a suite of reactions, via biotic (Fig. 1a.i) or abiotic mechanisms, leaving large, dark-colored (humic) compounds rich in carbon and nitrogen (N) (Lehmann and Kleber 2015). It was thought that the formation of these compounds renders SOM resistant to further decomposition. However, the hypothesis that chemical recalcitrance (and humification) explains SOM stability has gradually fallen out of favor, due in part to a failure to demonstrate that the secondary synthesis of humic compounds is relevant in natural systems (Keiluweit et al. 2015; Kleber et al. 2015), but also because chemically complex compounds (e.g., lignin) that should theoretically persist in soil can have shorter residence times than supposedly labile molecules (Schmidt et al. 2011). In response, new theories of SOM formation and stability have begun to emerge (Schmidt et al. 2011; Riley et al. 2014; Cotrufo et al. 2015; Lehmann and Kleber 2015; Castellano et al. 2015; Dwivedi et al. 2017a).

Many of these emerging theories put interactions between physical and geochemical conditions and biology at the forefront of SOM stability. Schmidt et al. (2011), for example, concluded that SOM persistence depended on ecosystem properties; they described an array of characteristic environmental traits that contribute to the stability of SOM and sequestration of carbon in soils. These properties include climate, physical properties (e.g., the formation of preferential flow paths within soils that determine the degree of connectivity between microbes and substrates—Hunt and Manzoni 2016), chemical properties (e.g., soil pH or the chemical quality of plant exudates and litter—Castellano et al. 2015), physical protection mechanisms (e.g., the sorption of organic matter onto mineral surfaces, Fig. 1b.iii—Kleber et al. 2015), soil redox status, the ecology of belowground biota (including bacteria, archaea, fungi and soil macrofauna, Fig. 1b.v—Bardgett and Van Der Putten 2014), and the metabolic constraints on OM mineralization. In such a network of factors, the diverse biotic components interact chemically and physically with each other, and with the surrounding heterogeneous matrix of minerals and soil aggregates, resulting in the emergence of SOM.

The majority of SOM is decomposed by microbes and respired into the atmosphere over short (< 50 year) time scales. However, a significant portion of this SOM can become chemically or physically protected—before or after microbial processing—and persist for much longer (100 to 1,000 year) time scales. Several studies have recognized the importance of diverse microorganisms in the formation of SOM (Allison et al. 2010; Knicker 2011; Cotrufo et al. 2013). The Microbial Efficiency-Matrix Stabilization framework, for example, describes SOM dynamics resulting from individual physiologies of diverse microorganisms (Cotrufo et al. 2013). The concept behind the “framework” asserts that plant exudates and litter, while being significant components of SOM themselves (Castellano et al. 2015), are decomposed (Fig. 1a/b.ii), assimilated (Fig. 1b.v), and processed by the microbial community.

Plants play an active role in soil carbon cycling through litterfall and root exudation, thereby influencing microbial growth (Fig. 1b.iv). For instance, roots excrete small-molecular-weight compounds to aid plant acquisition of nutrients and water. In addition, when plants supply small molecular-weight compounds such as sugars or reducible compounds, they prime microbial growth in soil. In general, the ‘priming’ of the soil microbial community can stimulate the production of exoenzymes near the root zone, increasing depolymerization of large biopolymers to monomers, which can become accessible to plant roots (Fig. 1b.i) (Kuzyakov et al. 2000; Gunina and Kuzyakov 2015). Moreover, roots also play a significant role in SOM stability through activities aimed at liberating and mobilizing nutrients (Fig. 1b.vi). Organic acids and chelators can disrupt the ligand bonds between mineral surfaces and organic matter or nutrients. Such actions can subsequently release previously immobilized phosphorus and other inorganic nutrients (Hocking 2001; Wu et al. 2018), and destabilize mineral-associated organic matter (MAOM) (Keiluweit et al. 2015; Jilling et al. 2018).

Overall, microbial processing of decomposed plant exudates and litter results in the microbial “exudation” of different metabolic byproducts, including enzymes, lipids, glycoproteins, extracellular polysaccharides (EPS), and the eventual accumulation of chemically heterogeneous necromass (Fig. 1b.iv) that are the major precursors to stable SOM (Miltner et al. 2012; Cotrufo et al. 2013, 2015; Kögel-Knabner 2017). Indeed, recent evidence suggests that the chemical composition of persistent SOM more closely resembles that of microbial cells and their byproducts than that of plant materials (Kögel-Knabner 2017). Decomposition-resistant SOM often exhibits a narrow carbon-to-nitrogen (carbon:nitrogen) ratio (between 8 to 10) that is more representative of microbial residues than the wide C:N ratio (upwards of 30) characteristic of plant materials (Baisden 2002; Conen et al. 2008). These microbial cells and excretions are often stabilized via organo–mineral and organo–metal–oxide interactions.

Physicochemical associations between soil minerals and plant or microbial-derived organic matter can result in the protection of theoretically thermodynamically unstable (or labile) compounds within soil for decades. Soil minerals often constitute a mixture of phyllosilicates, oxides, and hydroxides, as well as short-range order and noncrystalline minerals in some soil types. These particles are generally found in the clay (< 2 μm) size-fraction of soils. The accessibility of SOM to microorganisms can be hindered through mineral interactions by adsorbing or occluding organic compounds and immobilizing exoenzymes responsible for the depolymerization reactions. Recent studies have demonstrated that different classes of microbial biochemicals show variability in the degree to which they adsorb onto mineral surfaces (Jagadamma et al. 2012, 2014). For example, Swenson et al. (2015) and Newcomb et al. (2017) have demonstrated generalizable patterns in the competitive sorption of different functional groups to the surface of ferrihydrite, an amorphous iron oxide. Phosphate-, N-, and aromatic-containing metabolites and carboxylates showed strong adsorptive behavior onto mineral surfaces (Swenson et al. 2015). Conversely, differential biochemical adsorption shapes the composition, stoichiometry (e.g., carbon:nitrogen:phosphorous ratio), and accessibility of the dissolved organic matter (DOM) pool. For example, the immobilization of exoenzymes on mineral surfaces can alter catalytic function by blocking the active site or changing enzyme structure (Quiquampoix and Ratcliffe 1992; Baron et al. 1999; Norde and Giacomelli 2000; Servagent-Noinville et al. 2000). Although adsorbed exoenzymes generally decrease in catalytic capacity, their residence times are higher than free exoenzymes (Yan et al. 2010). These mechanisms, in turn, can shape the structure, function, and activity of the microbial community, which can promote the production of enzymes or organic acids that mineralize available polymers and desorb stable SOM compounds (Kaiser et al. 2014).

Finally, empirical observations and modeling efforts have typically concentrated on surface soils (< 20 cm) for the estimation of carbon storage and turnover rates (Trumbore et al. 1995; Hicks Pries et al. 2018). However, over half of global SOM is stored in subsoils (> 20 cm), and this SOM is significantly older than surface SOM (Jobbágy and Jackson 2000; Mathieu et al. 2015). These patterns suggest a distinct balance of decomposition and preservation processes at depth relative to that at the surface. In all, this emerging picture of MAOM places a premium on understanding the interdependent impacts of soil hydrology, geology, geochemistry, and ecology across the complete soil profile. In particular, it calls for an improved understanding of the feedback between soil properties and the metabolisms of microbes (fungi, bacteria, and archaea). Insight into these emergent properties of soil carbon cycling has in particular benefited from improvements in approaches to study fine-scale MAOM dynamics.

Emerging technologies are shaping the picture of MAOM dynamics as affected by biology, chemistry, and soil structure. Advances in analytical techniques are revealing the chemistry of minerals and the composition of soil organic matter. Extremely detailed information on the composition of extractable soil organic matter can be obtained by using Fourier Transform Ion Cyclotron Resonance Mass spectrometry (FT-ICR-MS) and quantified by ¹H and ¹³C solution-state Nuclear Magnetic Resonance (Simpson et al. 2001; Tfaily et al. 2015). For nonextractable organic matter, solid-state NMR can characterize ¹³C labeled organic matter sorbed on surfaces of non-paramagnetic minerals (Clemente et al. 2012). FT-Infrared spectrometry is compatible with living systems and can detect functional group characteristics of MAOM, even when associated with paramagnetic minerals such as iron oxides (Omoike and Chorover 2004; Lehmann et al. 2007). The capacity for imaging soil in nondestructive ways has grown. Advancements in X-ray computed tomography of intact soils yields information regarding the soil structure, pore network and even roots, which can be processed into three-dimensional images (Taina et al. 2011; Mairhofer et al. 2012). More importantly, when combining imaging techniques with chemical analysis, components of MAOM can be spatially resolved. Techniques such as synchrotron radiation scanning transmission X-ray microscopy (STXM) coupled with near edge X-ray absorption fine structure (NEXAFS) has been used to determine the spatial heterogeneity of MAOM in two dimensions (Lehmann et al. 2008). A more detailed literature review on the benefits and limitations of advanced analytical techniques used for MAOM has been presented previously (Kleber et al. 2015).

The emerging understanding of SOM dynamics also calls for greater efforts to identify and quantify the role of microbial functional traits in determining the fate of organic compounds within soils. Microbial communities are complex adaptive systems (Levin 1999), in which heterogeneous collections of individual organisms that coexist, co-operate, and compete across wide-ranging environmental gradients shape the evolution and emergence of community structure. The fitness of individuals within a community is governed by their inherent physiological traits and the energetic costs of intracellular trade-offs (Maharjan et al. 2013). Microbial traits can be defined as physiological or morphological characteristics closely linked with ecosystem function and performance (Viollea et al. 2014). These traits can be broadly separated into response (determining how organisms respond to environmental change, such as drought) and effect (determining an organism's effect on ecosystem processes, such as nitrification) traits (Allison 2012). The expression of, and resource allocation to, different traits (e.g., exoenzyme production) is regulated by intracellular energetic trade-offs (Beardmore et al. 2011). These factors determine an individual microbe's ability to compete for resources within a wider community, and to tolerate and thrive under extraneous selective pressures, including fluctuating environmental conditions. Traits and trade-offs that promote optimal relative fitness will subsequently be selected for and fixed (at least temporarily) within a population (Razinkov et al. 2013).

Many relevant traits (e.g., osmolyte production or exoenzyme production) show deep evolutionary divergence and nonrandom distributions (Schimel and Schaeffer 2012), meaning that environmental selection for specific ecological strategies can lead to concomitant changes in a suite of genome-linked traits related to SOM formation and degradation. Reorganization of microbial community composition can also change anabolic trait distribution and select for traits that convert plant-derived carbon to metabolites (e.g., osmolytes), or necromass that may be preferentially stabilized on soil minerals (Liang et al. 2017). On a broader scale, community changes that result in an increase in the soil fungal:bacterial ratio can also result in higher SOM due to the higher biomass associated with the proliferation of fungi (Strickland and Rousk 2010). The extent of bacterial secondary metabolism (e.g., antibiotics, non-structural carbohydrates, siderophores, etc.), which plays a large role in shaping the chemical composition of microbial exudates and necromass, is beginning to be quantified through advances in genomic and metabolomic technologies. For example, Crits-Christoph et al. (2018) recently used large-scale metagenomic sequencing to reconstruct the genomes of a number of common soil bacteria. They uncovered the potential for extensive secondary metabolite production, suggesting a largely underestimated capacity of these bacteria to produce and exude biosynthetic products and contribute to the composition of SOM. However, while “omic” techniques provide greater understanding of the complexity of the microbial community and its metabolic diversity, the need for large amounts of material for nucleic acid extraction (covering multiple niches) and the destructive nature of that extraction, limit understanding of temporal dynamics. These dynamics include how community-composition interactions necessitate the development and application of microbial centric models that permit testing of hypotheses related to cell-cell interactions, response to environmental fluctuations, and the microbial feedback to carbon dynamics (Riley et al. 2014; Georgiou et al. 2017).

While these small-scale processes undoubtedly interact, giving rise to larger-scale phenomena, a tractable balance of simplicity and complexity must be sought in the development of parsimonious, yet sufficiently accurate, macro-scale models of SOM cycling. To this end, it is critical to thoroughly understand key SOM decomposition and stabilization mechanisms, where and when they are important, and how nonlinear interactions between components may compound. Changes in climate and vegetation may also cause important shifts in the relative role of key mechanisms, with broad implications for the terrestrial carbon cycle.

There is a suite of SOM dynamics models incorporated into land models ranging from simple pool-based CENTURY-like representations (Parton et al. 1987, 1998, 2010; Jenkinson et al. 2008) to more complex models such as the Biotic and Abiotic Model of SOM (BAMS; Riley et al. 2014; Dwivedi et al. 2017; Tang et al. 2019) and the COntinuous representation of SOC in the organic layer and the mineral soil, Microbial Interactions and Sorptive StabilizatION (COMISSION; Ahrens et al. 2015), which explicitly represents physicochemical processes and microbially mediated reactions. Here we briefly describe the traditional and some of the emerging modeling approaches. More importantly, this section aims to describe differences in the traditional and emerging modeling approaches, and the need for a reactive transport modeling framework representing relevant mechanisms.

Traditionally, mean residence times (MRT) are used to categorize soil organic carbon into different “pools”. The decomposition reaction rate is the inverse of MRT, and therefore pools with shorter MRT have higher inherent decomposition reactivity (Parton et al. 1987; Davidson and Janssens 2006; Kleber 2010). These pools have been extensively used to investigate SOM decomposition in models such as CENTURY, ROTH-C, and their progeny models. These pools are intended to represent from the most labile to the most recalcitrant carbon in CENTURY-Like models (Parton et al. 1987, 2010), and from humified organic matter to inert in ROTH-C models (Jenkinson and Coleman 2008; Jenkinson et al. 2008). Similarly, SOM decomposition rates are also empirically referred to as fast, slow, and passive for pools ranging from the most labile, or humified, to the most recalcitrant, or inert. These traditional approaches have been very widely applied, but lack representation of the underlying microbial and abiotic processes described above. They have also been shown to have very high parameter equifinality, in spite of being calibratable in many systems (Luo et al. 2015, 2017).

The traditional CENTURY-like models are formulated by assuming sufficient microbial activity to support organic matter decomposition, because litter bag and soil incubation experiments often observe consistent patterns of exponential organic matter decomposition over time (Davidson et al. 1995; Berg and McClaugherty 2014). Soil minerals are considered to exert only a multiplicative effect on the turnover rates of the different components of SOM. Other multiplicative factors are similarly applied to represent how temperature and moisture regulate SOM decomposition. However, spectrum analysis of field measurements such as concentration-gradient-based soil surface CO₂ efflux time series often show that there are many dynamic signals that cannot be resolved with traditional models (e.g., Vargas et al. 2010; Manzoni et al. 2018). Notably, the traditional models fail to capture the priming effect in soil, which involves strong interactions between microbes and their bioproducts, substrates, and soil minerals (e.g., Fontaine and Barot 2005). Furthermore, when organo–mineral interactions are involved, substrate bio-availability usually decreases. This observation suggests that in soils where soil minerals are plentiful, oligotrophic microbes are favored. However, model representation of oligotrophic microbes may not be necessary, considering that microbes can live in a dormant state, making it difficult to determine whether oligotrophic microbes are the major contributors to observed microbial activity or simply whether a smaller population of more copiotrophic organisms are active (Joergensen and Wichern 2018).

Although these simplified pool-based models have been important for estimating SOM stocks, they are unable to capture SOM dynamics beyond first order effects. For example, most site- and global-scale terrestrial SOM models apply first-order relationships, and therefore cannot predict explicit turnover times for different compounds. Further, a major obstacle remains in matching the model with data in terms of understanding which components of SOM are decomposed and which are stabilized, because the pools in the traditional models are not structurally well defined (i.e., not observable; Abramoff et al. 2018). Most importantly, these traditional modeling approaches cannot answer why organic matter that is thermodynamically unstable persists in soils, sometimes for a millennium (Schmidt et al. 2011). In view of these modeling inconsistencies, new models have been developed and applied to understand the effects of microbial processes and organo–mineral interactions on SOM dynamics.

Interactions between microbes and soil minerals can be explicitly represented in more mechanistic models—some of the earlier attempts can be found in Smith (1979) and Grant et al. (1993). The recently published batch (i.e., no vertical resolution) models (e.g., MEND; Wang et al. 2013, RESOM; Tang and Riley 2014, CORPSE; Sulman et al. 2014, and MIMICS; Wieder et al. 2014) use diverse approaches to formulate organo–mineral interactions. There are also several new models (described below in Section ‘Current use of reactive transport modeling of SOM dynamics’) that explicitly represent microbes, soil mineral interactions, and vertical structure.

The MEND model incorporates a mineral-adsorbed soil carbon pool affected by separate adsorption and desorption kinetics that result in a Langmuir isotherm at equilibrium. However, Wang et al. (2013) argued that by incorporating microbial dormancy into the MEND modeling framework, only a single microbial functional group is required to reproduce laboratory aerobic incubation measurements of biomass and CO₂ respiration rates. A similar approach is adopted by the CORPSE model. RESOM uses the equilibrium form of the Langmuir isotherm, which is then integrated with the Equilibrium Chemistry Approximation kinetics (Tang and Riley 2013; Tang 2015; Zhu et al. 2017) to represent interactions between dissolvable carbon, mineral surfaces, exoenzymes, and microbes. The MIMICS model avoids an explicit formulation of organo–mineral interactions, and instead lumps these interactions by changing the half saturation constants used to drive Michaelis-Menten type microbial kinetics. These models were evaluated together in a recent model intercomparison (Sulman et al. 2018); that study concluded that the observational data were insufficient to invalidate any of the model structures. The authors suggested that long-term observations of soil carbon input (from plants) and perturbation experiments (such as warming) are needed to better evaluate and constrain these newly developed models.

Soil is a complex system and considering its physical and chemical evolution is critical for accurately capturing SOM dynamics (Riley et al. 2019). The following discussion is intended to review modeling needs for the following important topics: (1) microbial dynamics, (2) mineral associated organic matter (MAOM) interactions, and (3) the effects of SOM molecular structure. We subsequently make recommendations as to how these processes can be integrated into a reactive transport model.

The structure, function, and response to perturbation of the soil microbial community plays an important role in shaping soil biogeochemistry and SOM dynamics (see Section ‘Microbial and plant processes’ above), motivating the active representation of the belowground community in next-generation model structures. Allison (2012), using the DEMENT model to explore how microbial traits would covary with the evolution of decomposing litter, found that incorporating microbes with different growth traits enabled the model to explain more than 60% of the variance in decomposition rates of 15 Hawaii litter types. Using a similar modeling strategy, Kaiser et al. (2014) found that microbial diversity is essential in explaining the stoichiometric evolution in decomposing litters. However, relatively few studies explored the influence of microbial diversity on carbon formation and decomposition in the mineral soil. When this is done, microbial community structure is often aggregated into the r versus K strategists that are adopted from the MacArthur and Wilson (1967) classification scheme. For instance, the MIMICS model (Wieder et al. 2014) used a copiotrophic group to represent the r-strategists and an oligotrophic group to represent the K-strategists. This r vs K representation of microbial community structure has been argued to be critical to accurately simulating the priming effect (Fontaine et al. 2003). However, models like MEND, CORPSE, and RESOM are able to simulate priming with only one microbial group, as long as substrate diversity is included in the model. In addition, Lawrence et al. (2009) showed that short-term SOM dynamics resulting from interactions between microbes and environmental conditions (e.g., drying–rewetting events) could not be adequately represented using a first-order representation of decomposition (as in Century or RothC). They suggested that while first-order models can capture the bulk response of SOM dynamics under steady-state conditions, the inclusion of exoenzyme and microbial controls on decomposition is necessary to appropriately simulate pulsed rewetting dynamics.

Traditional models that employ a function of soil texture to represent soil mineral effects on soil carbon dynamics, as a multiplicative modifier to decomposition rates, are likely underestimating the roles played by soil minerals. For example, soil minerals shift between different oxidation status due to particular electron orbital configurations, which enable them to act as either reductants or oxidants of organic matter. A number of microorganisms are able to take advantage of this chemical flexibility and couple these soil minerals with organic carbon to sequester energy when the most favorable oxidant, i.e., oxygen, is in short supply. Likewise, none of the models mentioned above sufficiently accounts for organo–mineral interactions by only considering sorption–desorption dynamics.

Soil minerals are important factors in determining soil aqueous chemistry. Minerals, microbial metabolic byproducts (e.g., organic acids), and inorganic compounds exert strong controls on soil pH and redox potential. However, few contemporary soil carbon models have taken such organo–mineral interactions into account. Rather, when they are needed, such as in modeling methane dynamics, pH and redox potential are formulated as multiplier functions without accounting for their mechanistic evolutions (Zhang et al. 2002; Li et al. 2004; Zhuang et al. 2004; Riley et al. 2011). A slightly more mechanistically based model is attempted in Tang et al. (2016), where they merged the Windermere Humic Aqueous Model (WHAM; Tipping 1994) with the CLM-CN model to simulate the effect of iron cycling on CO₂ and CH₄ production in anoxic Arctic soil microcosms. Their model also simulated the evolution of pH using the aqueous chemistry capability provided by WHAM. However, the decomposition dynamics follow the traditional century-like formulation, so the model is therefore still not mechanistically representing organo–mineral interactions. Zheng et al. (2019) recently updated the framework from Tang et al (2016) with the PHREEQC 3.0 model (Charlton and Parkhurst 2011) and applied the resultant code to incubated organic soils from an Arctic tundra site. Although their simulations indicated some favorable comparisons with CO₂ and CH₄ fluxes, the representation of organo–mineral interactions was not improved.

The most comprehensive model integrating soil pH, redox, and sorption-desorption aspects of organo–mineral interaction is likely the ecosys model (e.g., Grant et al. 2003, 2017a, b), which includes a comprehensive parameterization of aqueous electrolyte dynamics and their relations with soil pH and ion-exchange capacity (Grant and Heaney 2010). However, this parameterization was primarily used to model inorganic phosphorus dynamics. As such, the role of soil minerals (e.g., ferrihydrite) as alternative electron acceptors for organic carbon decomposition is not modeled.

In summary, although much of the necessary chemistry knowledge exists, a modeling framework that systematically evaluates how organo–mineral interactions will play out in various aspects of soil biogeochemistry has not been developed (but see Riley et al. 2019 for a possible path forward in this regard). Further, more knowledge on how organo–mineral interactions affect biological aspects of microbial dynamics is needed. For instance, there is not a good method applicable to large-scale models to represent how organo–mineral interactions contribute to the formation of biofilms at different places within the soil, even though such structures are of great significance in determining how microbes deal with fluctuations in soil water (Or et al. 2007).

Most land models use CENTURY-Like SOM pools to represent SOM dynamics by applying pseudo-first-order approximations (Jenkinson et al. 2008; Parton et al. 2010). Although these simplified representations of SOM compounds are informative, some investigators have argued that the structures of diverse organic matter compounds may need to be represented to capture SOM dynamics, for several reasons. First, the Century-like pools are largely qualitative and vague, thereby leading to large uncertainty when confronting predictions of SOM dynamics with observations. In general, it is not possible to precisely categorize pools of SOM by turnover times that have such a broad range (e.g., a few days to thousands of years). In nature, the formation or decomposition of SOM compounds are properties of the ecosystem (Schmidt et al. 2011), and the traditional qualitative structural classification is a major obstacle toward providing a mechanistic basis for SOM dynamics. For instance, a group of SOM compounds may have a completely different response to soil moisture and temperature than another group of SOM compounds, even though these distinct groups may have comparable turnover times. Moreover, organo–mineral interactions may also show contrasting behaviors; for example, a group of compounds falling into comparable intrinsic decomposition rates may have different affinities to soil minerals or responses to environmental conditions (Gordon and Millero 1985; Gu et al. 1994; Mikutta et al. 2007; Kleber et al. 2011). Overall, we argue that there is a need to represent SOM compounds using their molecular structures in reactive transport models along with physical protection mechanisms (e.g., MAOM) and different functional groups of microbes (e.g., Riley et al. 2014; Dwivedi et al. 2017; Tang et al. 2019.)

As discussed in the Introduction and the sction ‘Current Use of Reactive Transport Modeling of SOM Dynamic’, SOM decomposition is governed by environmental conditions such as physical heterogeneity, physical disconnection, plant inputs, microbial diversity, and microbial activity, as well as the molecular structure of organic matter (Schmidt et al. 2011; Riley et al. 2014; Dwivedi et al. 2017a). To capture such complex physico-chemical and biological interactions, we propose that a vertically resolved reactive transport modeling framework is required.

Flow and reactive transport processes influence SOM formation and decomposition. In particular, geochemical processes can affect the mobility of chemical species through different mechanisms of dissolution-precipitation, sorption-desorption, ion-exchange, redox-reactions, complexation, and colloidal interactions (Arora et al. 2016, 2018; Yabusaki et al. 2017; Dwivedi et al. 2018b, a). There are a variety of RTMs and numerical codes—such as TOUGHReact (Xu et al. 2006; Maggi et al. 2008; Gu et al. 2009), PFLOTRAN (Hammond et al. 2014), CRUNCH (Steefel et al. 2015), ecosys (Grant 2013), BAMS (Riley et al. 2014; Dwivedi et al. 2017a; Tang et al. 2019), and BeTR (Tang et al. 2013; Tang and Riley 2018)—that are available to describe and can represent the interaction of various complex and competing biogeochemical processes across spatial and time scales. These numerical codes can represent aqueous, gaseous, and sorbed phases, along with a large suite of processes—including hydrology, soil energy dynamics, advection and diffusion, aqueous speciation, minerals precipitation and dissolution, microbial dynamics, adsorption and desorption, and equilibrium and nonequilibrium chemical reactions.

Several studies have demonstrated the use of reactive transport models to address SOM decomposition and stabilization (Grant et al. 2003; Ahrens et al. 2015; Dwivedi et al. 2017a; Georgiou et al. 2018). These models explicitly represent physical protection mechanisms and different functional groups of microbes. BAMS also represented a reduced number of SOM molecular structures (Riley et al. 2014; Dwivedi et al. 2017; Tang et al. 2019). Here, we provide a brief description of how these processes are represented in reactive transport models.

Organo–mineral interactions are represented in a fashion similar to the isotherm-based approach in the single layer batch model. For instance, BAMS1 (Riley et al. 2014) adopted linear kinetics to model sorption and desorption processes, whereas in BAMS2 (Dwivedi et al. 2017a) the formulation was replaced with the Surface Complexation Model (SCM), allowing variable biogeochemical conditions within a thermodynamic framework. Both BAMS1 and BAMS2 were able to reproduce the almost exponentially decreasing profile of soil carbon content and increasing age of soil radiocarbon with depth. The COMISSION model (Ahrens et al. 2015) similarly considers interactions between microbes, minerals, and substrates, using Langmuir kinetics and reactive transport modeling, and was able to obtain analogous results. In contrast, ecosys uses Freundlich kinetics to model the sorption dynamics of soluble carbon to soil mineral surfaces and unhydrolyzed organic matter (Grant et al. 2003). Of all these models, ecosys is the only one that also considers a full range of terrestrial ecosystem properties; it has been shown in many studies across various ecosystems to be able to capture much of the observed temporal variability in net ecosystem carbon exchanges (e.g., Grant et al. 2003, 2006, 2017a, b). Overall, these modeling exercises suggest transport combined with soil–mineral-regulated SOM stabilization and microbial decomposition help explain observed soil carbon content and age profiles. Furthermore, the specific mathematical form used to represent organo–mineral interaction plays a minor role, given the limited data available to evaluate the models.

Riley et al. (2014) made the first attempt to explicitly represent molecular structures of SOM compounds and integrate processes into a reactive transport model. Because SOM compounds consist of several plant-synthesized, degraded compounds, and microbial biomass as well as necromass, it is difficult to include a full representation of each species. Therefore, they grouped SOM compounds based on different metrics such as O:C ratio, charges (positive or negative), and degree of polarity. Riley et al. (2014) include above- and belowground litter, as well as root exudates for carbon inputs. The litter was degraded into several simpler organic compound groups such as cellulose, hemicellulose, lignin, monosaccharides, phenols, amino acids, lipids, and nucleotides. Root exudates were considered to be simpler organics, and were included among monomers (i.e., monosaccharides, amino acids, organic acids, and lipids pools). Subsequently, Dwivedi et al. (2017) modified this reaction network by collapsing all the monomers into a single pool. Although observations of individual SOM compounds were not available for comparison with predictions, individual SOM compounds highlighted the processes predicted to control SOM stocks and cycling.

To a certain extent, these reactive transport models demonstrated the value of including different microbial functional groups. BAMS included fungi and aerobic bacteria functional groups by assigning affinities to decompose plant litter and SOM (Neely et al. 1991; Romaní et al. 2006; Thevenot et al. 2010; DeAngelis et al. 2013). Subsequently, dissolved organic carbon (DOC) and peptidoglycan (i.e., necromass or dead cell wall material (DCWM); Frostegard and Baath 1996) were produced upon death of these microbial groups. Even though substantially underrepresenting functional diversity in soils, these models indicated key processes and pulse responses, and the value of examining SOM dynamics with this type of approach. It is important to realize that the governing ecological functions in the ecosystem are based upon biodiversity, and it is very data-intensive to decipher linkages between ecosystem function and biodiversity (Goldfarb et al. 2011).

There are some models that have gone beyond representing microbial diversity by considering soil fauna in the decomposition dynamics of SOM. For instance, Hunt et al. (1987) modeled the detrital food-web in a grassland by representing 15 groups of organisms, ranging from bacteria and fungi to predacious mites. The simulations show that trophic interactions are essential in maintaining the ecosystem wide nitrogen, echoing the message by Kaiser et al. (2014) that biodiversity is important for enabling nutrient-based regulation of the carbon cycle. Considering that SOM decomposition always involves a spectrum of bioreactions that include nitrification, denitrification, methanogenesis, and methanotrophy, it may be important to consider microbial diversity when all these different bioreaction pathways and different gas fluxes are to be modeled. For a carbon-only model, before more comprehensive data are available to benchmark model evaluations, it is hard to justify whether microbial community structure should be included in the model formulation.

Below, we describe the most important governing equations for simulating reactive flow and transport in the continuum subsurface environment.

Multiphase flow (e.g., liquid and gas) or the Richards equation can be used to describe variably saturated flow. The Richards equation assumes only one active phase (e.g., liquid), while the gas phase is assumed to be passive.

The Richards equation (Richards 1931) has been widely used in the literature for describing partially saturated flow assuming a passive air phase present at atmospheric pressure (e.g., Neuman 1973; Panday et al. 1993; Dwivedi et al. 2016, 2017b). The saturation-based formulation of the Richards equation is given as:

where ∇ denotes the divergence operator; h [m] is the hydraulic head; K [m s⁻¹] is the hydraulic conductivity tensor; t [s] is the time; ϕ [m³ void m⁻³ medium] is a dimensionless quantity, which represents the porosity of the porous media; S_a [m³ H₂O m⁻³ void] and dimensionless quantity k_ra [–] denote the saturation of the aqueous phase and the relative permeability, respectively; S_s [m⁻¹] is the specific storage coefficient; Q_a [m³ H₂O m⁻³ medium s⁻¹] is the volumetric source or sink term.

The corresponding pressure-based formulation is given as:

ρ_f [kg m⁻³ water] represents the fluid density; ∇P [kg m⁻¹ s⁻², or Pa m⁻¹] represents the gradient of the fluid pressure; g [m s⁻²] denotes the acceleration due to gravity; e_z is a unit vector in the vertical direction; k_sat [m²] is the permeability tensor for fully saturated conditions.

The occurrence of components like CO₂ or H₂O in various phases (gas, liquid, solid) can be described using multiphase flow. Multiphase flow is described by solving the mass conservation equation of multiphase transport of each component j in phase α:

where Y_ja is the mass fraction of component j in phase α; m [kg s⁻¹ m⁻¹] is the dynamic viscosity; and all other variables (e.g., k_ra, etc.) are same as defined earlier.

There are several models in the literature that can be used to describe the relationship between the aqueous phase saturation and capillary pressure. However, Brooks-Corey and van Genuchten formulations are the most well-known. In addition, the Mualem and Burdine formulations are typically used to describe relative permeability functions (Brooks and Corey 1964; Mualem 1976; van Genuchten 1980). We encourage readers to explore user's manuals of numerical codes such as TOUGHReact, PFLOTRAN, or Crunch for more details on these models.

To simulate reactive transport, the transport and biogeochemical reactions are combined in the continuity equation as follows:

The left-hand side of the Equation (4) denotes the accumulation term (mol m⁻³ medium s⁻¹), which is the product of the porosity (φ) and liquid saturation (S_L), and the concentration (C_i). On the right-hand side, the first and second term together describe advective-diffusive transport; the third, fourth, and fifth terms represent various reactions that are partitioned between aqueous-phase reactions (R_r), mineral reactions (R_m), and gas reactions (R_l), respectively. Furthermore, a stoichiometric relationship is used to compute transformation from reactant to product species. In Equation (4)$Di*$ denotes the diffusion coefficient, which is specific to chemical species considered as indicated by the subscript i; n_ij are the stoichiometric coefficients of reactant j in reaction i, whereas there are N_r, N_m, N_g reactants in aqueous, mineral, and gas phases, respectively.

These reactions can further be classified into equilibrium-based or kinetic. Equilibrium-based reactions typically require the assumption that there exists a local equilibrium between reactants and products. However, rate-based reactions (kinetic) are always the more general (Steefel and Lasaga 1994). SOM decomposition involves microbial reactions, which we describe below.

The microbial decomposition rate of substrate C_i can be described using Michaelis-Menten (Michaelis and Menten 1913) kinetics:

where μ_i (s⁻¹) is the maximum specific consumption rate of substrate i (s⁻¹); O₂ (mol O₂ m⁻³) is the aqueous O₂ concentration; and B (mg C-wet-biomass L⁻¹) is the wet biomass carbon of microbial functional groups considered in the reaction network (e.g., bacteria and fungi). The effects of pH, saturation levels of soil, and temperature (T) can be important for decomposition rates (e.g., Schimel et al. 2011). Equation (5) can be modified to using functions f₁(θ), f₂(pH), and f₃(T) to describe these environmental effects on microbial activity (Maggi et al. 2008), as follows:

While the Monod-Michaelis-Menten kinetics, including the related multi-substrate progeny, have enjoyed a wide range of applications, Tang and Riley (2013) showed that the Monod kinetics fail to account for the substrate limitation in their approximation to the law of mass action, and misplace such limitation as a linear competition by juxtaposing many Monod terms (e.g., for Eqn. (5), the Monod representation of competition may lead to simple addition of more Monod terms without accounting for one microbe's influence on another microbe's K parameter (Tang 2015; Tang and Riley 2017). This shortcoming will result in poor model performances under conditions of substrate limitation, a situation that is characteristic of soil. Further, the Monod kinetics will oversimplify the connections between substrates and consumers when describing a multi-substrate—multi-consumer network, resulting in the modeled fluxes being oversensitive to the kinetic parameters that would be formulated in the corresponding law of action model. The Equilibrium Chemistry Approximation kinetics (ECA) and its extended SUPECA (synthesizing unit plus ECA) kinetics are both able to adequately address these two shortcomings of the Monod kinetics, and also provide a very straightforward way to link the kinetic parameters with thermodynamics and biological traits (Tang and Riley 2013, 2017; Zhu et al. 2016a, b).

In most formulations of microbial growth in reactive transport models, substrate use efficiency is the only parameter to describe the partition of an assimilated substrate into either biomass or respired product. However, the catabolic and anabolic separation is more complex than this simple representation. When one takes into account building up metabolic reserves, cell maintenance and growth, and the enzyme and polymer exudation involved in microbial physiology, the emergent CUE can have strong nonlinear variations with respect to the environmental conditions, such as temperature or moisture (Tang and Riley 2014; Allison and Goulden 2017). CUE may be further modified by regulations of nutrient availability (Manzoni et al. 2018). All these conclusions suggest a more nuanced parameterization of substrate use is needed in future modeling of microbial-activity-regulated SOM dynamics.

SOM dynamics are influenced by minerals present in the geochemical system. For example, pH affects SOM reactions, and the presence or absence of calcite minerals may therefore affect SOM dynamics. Similarly, iron minerals (e.g., ferrihydrite) can be important for SOM dynamics, as a significant fraction of the aqueous phase of SOM can undergo sorption and stay protected from microbial reactions (Dzombak and Morel 1996; Dwivedi et al. 2017a). Mineral precipitation and dissolution reactions can be described using Transition State Theory (or TST) type rate laws:

where rate constants for neutral, acid or additional (j_th) reaction mechanisms are denoted by k_neutral, k_H⁺ and k_j, respectively; a_ij indicates the activity of the i_th aqueous species in the j_th reaction; finally, Q_m denotes the ion activity product of the m_th mineral phase and K_{eq, m} denotes its corresponding equilibrium constant.

Sorption is a complex process, one that depends upon organic molecule characteristics, surface area, and site density of minerals, as well as aqueous chemistry (Dudal and Gérard 2004). There are several approaches that can be used to describe sorption processes and MAOM. However, kinetic and equilibrium adsorption isotherms (e.g., linear, Langmuir, Freundlich; Davis and Kent 1990) have been widely used in the literature to describe SOM dynamics (Grant et al. 2011; Mayes et al. 2012; Riley et al. 2014; Ahrens et al. 2015). Typically, sorption isotherms describe a functional relationship between adsorbate and adsorbent for a constant-temperature condition. Correspondingly, empirical sorption isotherms have been developed to describe both kinetic (i.e., sorption site is slow and assumed to achieve adsorption equilibrium described using forward and reverse rates) and equilibrium (i.e., sorption site is assumed to reach adsorption equilibrium fast) systems. Because these empirical formulations cannot account for the out-of-sample effects of chemical conditions for which they are derived, thermodynamic surface complexation models (SCM) were also used (Dwivedi et al., 2017). A simple linear sorption relationship can be imposed using forward (adsorption; k (s⁻¹)) and reverse (desorption; k_r (s⁻¹)) rates. In the absence of any competing source or sink of a specific SOM species, an effective equilibrium linear sorption relationship will become as follows:

Nonlinear kinetic models require estimates of these rate constants (forward and reverse) and an additional exponential parameter (Goldfarb et al. 2011). On the other hand, SCMs offer a mechanistic approach for representing sorption processes and are generally more robust and applicable over variable geochemical conditions (Dzombak 1990; Goldfarb et al. 2011). The SCM models describe the sorption of solutes on solid surfaces as a chemical reaction between aqueous species and surface sites (surface complexation) (Dzombak 1990).

where A and S are the sorbate molecule and sorbent site (surface), respectively; and A − S is the sorbed complex. The sorption reactions are described by a mass action law equation assuming a local equilibrium:

where K_Eq denotes the equilibrium sorption constant; [A − S] is the concentration of surface complexes (mol kg⁻¹water); [A] denotes the concentration of sorbate molecule (mol kg⁻¹water), and [S] denotes the sorbent site (mol kg⁻¹water).

Although we are not aware of any study considering rate-limited and equilibrium-based models together to describe sorption processes, we argue that rate-limited and equilibrium-based sorption do co-occur in soils. We suggest future modeling approaches consider rate-limited and equilibrium-based sorption to determine the relative importance of these processes for SOM dynamics.

Soils store vast amounts of terrestrial organic carbon, more than the atmosphere and terrestrial vegetation combined. This carbon is vulnerable to release to the atmosphere under a changing climate. Although the mechanisms leading to the decomposition of SOM are well documented, uncertainties persist regarding the stability of SOM. To address this critical gap, a robust predictive understanding and modeling of SOM dynamics is essential for examining short- and long-term changes in soil carbon storage and feedbacks with climate. In this chapter, we reviewed recent research that improves emergent understanding of the important factors contributing to SOM stability. While there currently exists a suite of models representing SOM dynamics that span a range of complexity, some recent mechanistic models are more consistent with this emergent understanding of SOM persistence. Yet even those recent models do not represent several processes that can be important for SOM dynamics. We conclude that the next-generation models need to represent the full spectrum of quantitatively important mechanisms for determining SOM persistence—including rate-limited and equilibrium-based sorption, the formation of soil aggregates, representative soil minerals, microbial community dynamics, and vegetation interactions—to accurately predict short- and long-term SOM dynamics. Because this recommendation is obviously challenging, we have assembled an open-source, reactive-transport based SOM model that can be used to robustly integrate many of these processes (called BeTR-S; https://github.com/BeTR-biogeochemistry-modeling/sbetr) and invite the community to experiment with it (Riley et al. 2019).

Overall, it is important to understand SOM dynamics because SOM decomposition gives rise to potential greenhouse gases. Microbial dynamics, MAOM, and the molecular structure of SOM compounds have implications for understanding and modeling the short- and long-term responses of soil carbon stocks under local and regional climatic perturbations. Finally, this chapter illustrates the need to evaluate SOM dynamics in a reactive transport modeling framework that includes ecosystem properties.

This material is based upon work supported by the U.S. Department of Energy, Office of Science, Office of Biological and Environmental Research (as part of the Watershed Function Scientific Focus Area), and Office of Science, Office of Advanced Scientific Computing (as part of the project “Deduce: Distributed Dynamic Data Analytics Infrastructure for Collaborative Environments”) under Contract No. DE-AC02-05CH11231. Jinyun Tang acknowledges support from the Next Generation Ecosystem Experiment (NGEE-Arctic). The USDA NIFA Postdoctoral Fellowship program supported Katerina Georgiou. Nicholas Bouskill acknowledges support from a DOE Early Career Research Project (# FP00005182). William Riley acknowledges support from the Terrestrial Ecosystem Science Scientific Focus Area of Berkeley Lab. We thank Diana Swantek of Berkeley Lab for assistance with preparing Figure 1.