Biogeochemical processes are of tremendous importance for determining the fate of many organic and inorganic compounds in the subsurface. Most global elemental cycles involve biogeochemical transformation, and the recycling of carbon and nutrients relies almost exclusively on biogeochemical processes. In particular, the majority of natural organic compounds are biogeochemically reactive, but also a large number of anthropogenic organic carbon compounds can be biogeochemically transformed, for instance, during the biodegradation of organic contaminants. Furthermore, inorganic compounds such as e.g., many nitrogen, phosphorus or sulfur compounds, metal compounds or minerals are directly or indirectly affected by biogeochemical reactions. To which extent and at which conditions a biogeochemical reaction takes place depends not only on the properties of the involved chemical reactants and products but also on the behavior of the microbial community (or communities) catalyzing the biogeochemical transformation. Porous media—in particular natural porous media—are complex and often heterogeneous structures, which imposes severe challenges in determining the exact physical, chemical and ecological conditions the microbial community is exposed to and to which extent it is able to provide any ecosystem service, such as the catalysis of a biogeochemical reaction.

Reactive transport models have become an established mean for the investigation and quantification of countless chemical transformations in porous media (e.g., Lichtner et al. (1996), Xiao et al. (2018), other chapters of this issue), and modeling approaches addressing the biogeochemical dynamics in natural porous media such as soil, aquifers and aquatic sediments exist already for decades (Berner 1980; Lichtner et al. 1996; Boudreau 1997; Murphy and Ginn 2000; Barry et al. 2002; Brun and Engesgaard 2002; Meysman et al. 2003; Thullner et al. 2007; McGuire and Treseder 2010; Meile and Scheibe 2018, 2019). Over the last years the improvement of experimental techniques has led to a better knowledge on the behavior of microorganisms that drive biogeochemical reactions and, more broadly, to an increased understanding of coupled reactive-transport processes in the subsurface. This has led to increasingly sophisticated geomicrobial modeling approaches which have been used already in a reaction-transport framework or which provide the potential for being used to simulate biogeochemical processes and their dependency on the abundance and activity of microbial communities in subsurface environments.

The aim of this chapter is to build upon established modeling approaches for the simulation of biogeochemical processes in soils, aquifers and sediments across scales and to provide an overview of recent extensions of these established approaches towards increasingly complex representations of microbial dynamics in the subsurface. We thus start with section ‘Traditional Approaches for Simulating Biogeochemical Processes in the Subsurface’ giving a brief summary of these established modeling approaches. This will be followed in section ‘Thermodynamically Informed Models’ by a discussion of thermodynamic constraints on microbial activity and how to include such constraints into the simulation of biogeochemical reactions. Thereafter two sections will focus on the explicit representation of microbial dynamics. In recent years the increasing computational power as well as the increasing amount of experimental information on microbial systems has led to numerous new approaches for modeling the behavior of individual microbial species as well as of microbial communities (Song et al. 2014) and only a limited overview can be provided here. We separate these approaches into two (certainly overlapping) aspects: section ‘Integration of Microbial Sequencing Data in Geomicrobial Models’ describes approaches of high complexity which consider a maximum of available data, mainly from various sequencing techniques, and section ‘Considering the Ecological Behavior of Microorganisms’ describes approaches which put more emphasis on the ecological behavior of the microorganisms represented by certain traits and which are (implicitly) based on the assumption that the complexity of an ecological model should be limited (May 1974). Stable isotope approaches have become increasingly established for the investigation of biogeochemical processes and stable isotope data can be used to constrain the elemental fluxes between reactants and the microbial biomass or reaction products, and to identify which microbial pathway(s) has or have been involved into a degradation reaction. The section ‘Considering Stable Isotope Signatures’ provides an overview on how such information can be embedded in reactive transport modeling approaches. Modeling microbially driven processes in a porous medium implicitly involves various spatial scales ranging from the micro scale on which microbial cells act to the observation scale which ranges from the cm scale in high resolution laboratory experiments up to the 100 km scale in global land and ocean models. In section ‘Scale effects’ some consequences of linking microbial and transport scales are discussed and approaches on how to consider them within a modeling framework are presented.


Reaction networks for biogeochemistry

Microorganisms generate the energy required for growth and maintenance by catalyzing the transfer of electrons from electron donor substrates (ED) to terminal electron acceptors (EA). As a result, they exert a major control on the redox chemistry of their immediate surroundings (Thullner et al. 2007). In surface and subsurface systems such as rivers and lakes, sediments, soils and groundwater, the vastly dominant source of electrons is most often reduced carbon bound in complex organic molecules of natural or anthropogenic origin. Their progressive degradation (typically oxidation) shapes the redox conditions of the surrounding environment and classically leads to a redox stratification. This zonation results from the sequential utilization of EA's (O2, NO3, Mn(VI), Fe(III), SO4=) followed by methanogenesis, an order which is consistent with a progressive decrease in Gibbs energy yields released through these metabolic pathways (Claypool and Kaplan 1974; Froelich et al. 1979; Stumm and Morgan 1996).

In Reactive Transport Model (RTM) applications, redox reactions have traditionally been simulated using the well-known Michaelis–Menten model for enzymatically catalyzed reactions (e.g., Van Cappellen and Gaillard (1996); Boudreau (1997); Barry et al. (2002)) but other approaches have also been proposed, e.g., the ‘reverse Michaelis–Menten’ model (Schimel and Weintraub 2003) or the ‘equilibrium chemistry approximation’ (Tang and Riley 2013). The Michaelis–Menten model derives from theory of enzyme kinetics and is consistent with observations that show saturation behavior with increasing availability of ED substrate (Thullner et al. 2007). Field and laboratory observations also reveal similar saturation behavior with respect to the concentrations of EA's, and the reaction rate for a redox metabolic pathway is classically represented as:


where k is the maximum rate and KED and KEA denote half-saturation constants with respect to the electron donor and acceptor, respectively. Inhibition terms can be introduced in the rate law to account for the suppression of a specific metabolic pathway by the presence of higher energy-yielding EA. This classical approach to simulate the redox zonation has been extensively used in RTM applications of natural and engineered environments, see, e.g., the reviews by Thullner et al. (2007), Arndt et al. (2013) and Paraska et al. (2014) for further details.

(Sub)-surface environmental conditions can often be traced back directly or indirectly to the oxidation of reduced carbon molecules. As a result, an accurate representation of these processes, including their temporal evolution, is often of the upmost importance for characterization of the biogeochemical dynamics in these systems. Redox processes control, among others, the recycling of inorganic carbon and nutrients, the precipitation/dissolution of carbonate and sulfide minerals, the pH of (pore)-waters, the sorption/desorption of toxic elements as well as the carbon sequestration in soils and sediments, all of which have been addressed in RTM applications (see, e.g., Paraska et al. 2014; Li et al. 2017 for recent reviews). Models relying on kinetic rate laws similar to Equation (1) or variations of it have been applied across all spatial scales, from pore networks to the assumed horizontally homogeneous plot or core scale (e.g., Gharasoo et al. 2012; De Biase et al. 2013; Druhan et al. 2014), and from catchment scales (e.g., Bao et al. 2014) to regional and global scales (e.g., Thullner et al. 2009; Wania et al. 2010; Krumins et al. 2013; Raivonen et al. 2017; Hülse et al. 2018). They remain popular, although microbial rate laws explicitly accounting for the dependence of the microbial reaction rate on the abundance of the microorganisms are increasingly applied in RTMs, and begin to be used in large scale Earth System Models (e.g., Wieder et al. 2013, 2014; Sulman et al. 2014; Wang et al. 2016; Huang et al. 2018).

Kinetics of geomicrobial reactions

The Monod model for microbial growth.

Models in which microbial biomasses are explicitly included are typically applied to study the response of microbial communities to fluctuations in environmental conditions and the competition of different microbial groups for a common substrate (Boudreau 1999; Wirtz 2003; Thullner et al. 2005, 2007; Dale et al. 2008a; Regnier et al. 2011; Arndt et al. 2013). These biomass-explicit models have been coupled to RTMs to address environmental questions in a wide range of natural and engineered systems including water bodies (e.g., Vanderborght et al. 2002; Reed et al. 2014), groundwaters (e.g., Thullner and Schäfer 1999; Barry et al. 2002; Yabusaki et al. 2007; Li et al. 2009; Tartakovsky et al. 2009; Bao et al. 2014), lake sediments (e.g., Jin et al. 2013), marine sediments, (e.g., Dale et al. 2010; Regnier et al. 2011), and soils (e.g., Neill and Gignoux 2006; Wieder et al. 2013; Sierra et al. 2015; Huang et al. 2018).

Geomicrobial models rely on rate laws that explicitly account for the dependence of the microbial reaction rate on the biomass of the microorganisms, according to the Monod microbial growth equation (Monod 1949) or its derivatives (e.g., Soetaert and Herman 2009):

dBdt=μB=μmax BSKS+S

where dB/dt is the rate of biomass production, B is the biomass, µmax is the maximum specific growth rate, S is the concentration of the substrate, and Ks is the Monod half-saturation constant, that is, the substrate concentration when µ = 0.5 µmax. In a related model approach—the Contois model—KS is considered to depend linearly on B (Contois 1959). About a decade ago, the functional dependency on substrate and biomass of the classical microbial growth model was challenged by Schimel and Weintraub (2003) and the reverse Michaelis–Menten kinetics was proposed as an alternative to simulate soil C decomposition. This model assumes that the rate is a nonlinear function of biomass while depending linearly on substrate concentration (Wang et al. 2016). Several later studies using microbial models showed that the classical Michaelis–Menten model may indeed fail to resolve the biogeochemical dynamics in both laboratory and natural settings (see Tang (2015) for an overview). This debate revolving around alternative nonlinear microbial models has received considerable attention in global scale applications of soil C decomposition and its response to global change. Comparative analysis of the model's emerging properties has also recently been carried out (Wang et al. 2016) and it has been demonstrated that both models are two special cases of the ‘equilibrium chemistry approximation’ (ECA) kinetics proposed by Tang (2015).

A particular biomass group i may rely on more than one substrate to sustain its metabolic needs. This dependency includes catabolic reactions involving the transfer of electrons from an external electron donor (ED) to a terminal electron acceptor (EA) as well as the possibility for a microbial group to rely on j electron donors to sustain its metabolic needs (Dale et al. 2006). The generalized form of Equation (2) thus reads:

dBidt=Σj=1nμmax,ijBiEDjKij+EDj KiEA+EA=Σj=1nμmax,ijBiFK,ij

where n is the number of ED's sustaining the growth of biomass group i, and FK,ij are dimensionless kinetic limitation terms (0-1) for each of the j catabolic reaction that are function of the uptake efficiencies of the ED's and EA's by the microorganisms, as well as the local availabilities of the ED and EA. Note that Equation (3) can further be modified by including inhibition terms to account for the suppression of a specific metabolic pathway by the presence of higher energy-yielding EA's or for the presence of toxic compounds. This is usually realized by introducing a series of hyperbolic functions, each characterized by an inhibition constant Kin (see, e.g., Van Cappellen et al. 1993; Boudreau 1997). In case where microorganisms are exposed to a new substrate an initial delay of their growth and degradation activity may be observed as the enzymatic apparatus of the cells needs to adapt to the new compound and a lag effect needs to be considered in the modeling approach (Wood et al. 1995; Ginn 1999; Nilsen et al. 2012).

The rates of biomass growth and substrate utilization are directly coupled through the growth yield, Y, expressed as the amount of biomass carbon produced per unit of mass electron donor consumed (Regnier et al. 2011). Values of Y depend on the Gibbs energy generated by the catabolic reaction j, the Gibbs energy needed for the formation of a new biomass i (anabolism), and the efficiency with which the organisms utilize energy (VanBriesen 2002). In other words, the growth yield partitions the flow of electrons between the catabolic pathway of energy generation and the anabolic pathway of microbial growth (Rittmann and McCarty 2001); accurately predicting Y is thus key to properly represent geomicrobial processes in RTMs (see next section for further details). The change in electron donor concentration due to microbial processes can then be represented by:

dEDjdt=Σi=1m1Yijμmax,ijBi FK,ij

where m is the number of microbial groups relying on the substrate (electron donor) j to ustain their metabolic energy needs and Yij represents an observed growth yield, which is the efficiency of converting carbon into microbial products (Sinsabaugh et al. 2013; Bradley et al. 2018b). A similar equation can be written for the electron acceptor consumption, which can be catalyzed by multiple microbial groups relying on several electron donor substrates. The coupled Equations (3, 4) are at the foundation of geomicrobial models. An alternative approach to simulate microbially driven reactions is to constrain the substrate consumption through Michaelis–Menten kinetics and then quantify the microbial growth rate by multiplying the consumption rate by the yield (e.g., Thullner et al. (2007)). Both approaches are identical if Yij is constant but if it varies, the choice defines if changes in yield imply unchanged growth rate but changing substrate consumption (Monod) or unchanged substrate consumption but changing growth rate (Michaelis–Menten).

As further discussed in section ‘Thermodynamically Informed Models’ of this chapter, the computation of Y values depends on the Gibbs energies generated and consumed by the catabolic and anabolic (new biomass synthesis) reactions, respectively. The underlying principles to combine the catabolic and anabolic pathways in so-called macrochemical equations are described in detail in Dale et al. (2006) and Smeaton and Van Cappellen (2018). In short, electrons supplied by an ED are used partially to generate energy and partially for biomass synthesis, the partitioning of this electron flow depending of the stoichiometry of the half-redox reactions involved in catabolism and anabolism. As an example, the macrochemical equation for anaerobic methane oxidation catalyzed by methane oxidizing archea (MOA),


is obtained from the combination of the following two redox reactions :


for anabolism and catabolism, respectively. See Dale et al. (2006) for the determination of the stoichiometric coefficients in Equations (57). We thus note that only a fraction of the carbon source is fully oxidized to dissolved inorganic carbon (HCO3), while the other fraction is only party oxidized into biomass (C5H7O2N). Furthermore, in this reaction the MOAs produce a reactive intermediate (H2) that can be used by other anaerobic microorganisms to sustain their metabolic needs, most notably sulfate reducing bacteria (SRB). The latter group uses H2 as energy source according to the following macrochemical equation:

9.7 SO42+40.7 H2+10.5 H++0.2 NH4++CO30.2 C5H7O2N+9.7 HS+41.3 H2O

Reaction 5 is endergonic under standard state conditions (see below for a definition) and thermodynamic calculations have in fact shown that CH4 can only be oxidized to H2 under a narrow range of in-situ pressure, temperature and solution composition (LaRowe et al. 2008). This thermodynamic constraint requires maintaining low H2 concentrations, which is typically achieved through H2 consumption by the SRB (Hoehler et al. 1994; Dale et al. 2008b). The syntrophic association of MOA and SRB thus catalyzes the net reaction of anaerobic oxidation of methane coupled to sulfate reduction, which can be obtained by combining Equations (5,8). If MOA and SRB biomasses are not explicitly simulated, the combination of the two catabolic pathways CH4+3H2O4H2+HCO3+H+ and 4H2+SO42+H+HS+4H2O only leads to the much simpler net reaction stoichiometry for anaerobic oxidation of methane (assuming that the reactive intermediate (H2) produced by the MOA is entirely consumed by the SRB):


In RTMs, the classical approach is to implement the net reaction Equation (9) for simulating the coupled methane-sulfur cycles (Regnier et al. 2011), but several authors have also explicitly modeled the dynamics of reactive intermediates (e.g., Dale et al. 2008b; Orcutt and Meile 2008; Alperin and Hoehler 2009b), thereby providing useful insights into the mechanisms and environmental controls of microbially mediated reactions. This metabolic modeling approach is straightforward to implement in RTMs as, for the above example, it only requires the addition of a new mass conservation equation for the intermediate species, H2, while the rate expression for the net oxidation of methane by sulfate is replaced by two rate expressions, one for each of the individual reaction steps (Regnier et al. 2011).

Mortality and maintenance terms in geomicrobial models.

The death of microorganisms and the contribution of dead biomass to substrate pools are included in geomicrobial models by assuming that mortality scales to the active biomass pool via a first-order decay rate constant µe (e.g., Dale et al. 2006). If the electron donor is an organic carbon substrate fed by the death of microbes, the conservation equations for biomass and substrates are modified according to:

dEDjdt=Σi=1m(1Yijμmax,ijBi FK,ijzij μe,ijBi)

with zji as fraction of the dead biomass returning to the pool of substrate j.

The above two equations form the cornerstone of biomass-explicit modeling approaches in RTMs (e.g., Thullner et al. 2007). However, the observed growth yields in this model does not distinguish the energetic costs associated to growth from those due to maintenance (Bradley et al. 2018b). The classical models of Herbert (1958) and Pirt (1965) allow to separate the energetic costs required to generate new biomass (Lipson 2015) from the energetic costs to perform all maintenance functions. These classical models differ by the provenance of maintenance energy, Herbert considering maintenance costs as endogenous catabolism, i.e., the consumption of biomass while Pirt assumes that maintenance leads to an additional consumption of substrate. Both models have their inherent limitations (see, e.g., Wang and Post (2012); and references therein for details). In addition, observations show that the specific maintenance rate, and the provenance of maintenance energy, can vary under different environmental conditions (van Bodegom 2007). Together, these elements call for the development of an hybrid model that allows for alternative supplies of maintenance energy from biomass and/or substrate, as proposed by Wang and Post (2012) and Bradley et al. (2018b) for soils and sediments, respectively:

dBidt=Σj=1nμmax,ij Bi FK,ijmq,ij Bi(1FK,ij)μe,ij Bi

where mq,ij represents the specific maintenance rate, and the term in parenthesis tends to zero when substrates are abundant, consistent with the idea that microorganisms do not consume biomass for their maintenance requirements in this case. Electron donor substrates are consumed instead and calculated according to:

dEDjdt=Σi=1m1YG,ijμmax,ij Bi FK,ij+1YG,ijmq,ijBiFK,ijzijμe,ijBi

where YG,ij is now the “true” growth yield, that is a parameter that only reflects the expenditure of energy to generate new biomass (Lipson 2015). Conversely, when substrates are scarce, the second term in Equation (13) tends to zero while the term in parenthesis in Equation (12) approaches 1.


In what follows, we review several quantitative approaches aiming at the integration of thermodynamic constraints in RTMs that explicitly simulate geomicrobial population dynamics. We restrict our analysis to concepts derived from equilibrium thermodynamics, which have been extensively applied within the framework of RTMs over the last decade. Appealing extensions to such classic, equilibrium thermodynamic approaches are models relying on the Maximum Entropy Production (MEP) concept (e.g., Vallino 2010). In these models, MEP is used either as an optimization goal or as a governing principle to understand and model biogeochemical processes (Meysman and Bruers 2010; Vallino 2010; Song et al. 2014). To our knowledge, however, application of these non-equilibrium thermodynamic approaches, including those relying on statistical thermodynamics (Song et al. 2014), have remained largely theoretical and conceptual. Despite their significant potential, the modeling of microbial processes based on MEP is thus not covered in this chapter.

Incorporating bioenergetics in models of microbial dynamics

Microbial dynamics is constrained by classical thermodynamics in two ways: (1) Growth yields (Y) are dependent on catabolic energy gains (Roden and Jin 2011) and (2) redox reactions can only proceed when the energy yield of the catabolic reaction exceeds a metabolic threshold. This bioenergetic limitation can be included in the Monod geomicrobial model for redox processes (section ‘Traditional Approaches for Simulating Biogeochemical Processes in the Subsurface’) using a functional dependency on the thermodynamic driving force for the reaction, which depends on the Gibbs energy yield (Regnier et al. 2011):

dEDjdt=Σi=1m1Yijμmax,ij Bi FK,ijFT,ij

where the thermodynamic limiting term, FT, represents the limitation that energy yields have on catabolic reactions rates (Jin and Bethke 2002; Dale et al. 2008c; LaRowe et al. 2014). Therefore, the implementation of thermodynamic constraints in geomicrobial models require calculation of the energetics of catabolic reactions and subsequent quantification of the FT term as well as derivation of a relation between Y and the energetics of cellular metabolism.

Energetics of catabolic reactions

Modeling studies on the energetics of redox reactions catalyzed by microorganisms have mostly concentrated on the relative Gibbs energy yields associated with different TEAs (Arndt et al. 2013), often to predict their classical sequential utilization (O2, NO3, MnO2, Fe(OH)3, SO4) in subsurface environments (e.g., Thullner et al. 2007 and references therein). This well documented redox sequence may however no longer hold under non-standard state conditions prevailing in natural environments (Amend and Teske 2005; Bethke et al. 2011; LaRowe and Van Cappellen 2011) and sulfate may, for instance, become an energetically more favorable TEA than iron oxides (LaRowe and Van Cappellen 2011).

Recently, LaRowe and Van Cappellen (2011) have extended such classical thermodynamic analysis to a wide variety of natural electron donors, mainly organic matter compounds. The approach has the major advantage of not requiring structural information in order to estimate the energetic potential of complex, natural organic matter. Instead, it uses the average Nominal Oxidation State of the Carbon (NOSC) as a proxy to scale the bonding in organic compounds to their energetic content (Arndt et al. 2013). The NOSC is readily determined for a wide array of organic compounds as it only depends on the net charge Z and the stoichiometric numbers of the elements C, H, N, O, P and S in a given organic compound, respectively denoted by a, b, c, d, e and f:


Thermodynamic calculations for the oxidation half reactions of organic compounds show that Gibbs energies of reactions and NOSC values follow an inverse relationship (Fig. 1). This remarkable linear correlation can thus be used to estimate the energetics of catabolic reactions (ΔGcat) coupling organic matter oxidation with any TEA as long as the average NOSC, which is directly calculated from element ratios, is known. These findings are potentially of high value for reactive-transport modeling and, more specifically, geomicrobial modeling. Indeed, while much attention has been given to the simulation of the redox sequence in subsurface environments, much remains to be done regarding the quality and diversity of substrates available for microbial growth and how they shape the biogeochemistry of the subsurface. A better representation of the energetics of substrates in RTMs would thus help predict important patterns of organic matter degradation in natural environments. For instance, the degradation of labile organic carbon compounds (e.g., planktonic biomass, polysaccharides) leads to similarly high Gibbs energy yields in both oxic and anoxic settings (e.g., Henrichs (2005)). In stark contrast, the degradation rates of more refractory compounds such as lignins or lipids are much more sensitive to redox conditions (Canfield 1993; Canuel and Martens 1996; Hartnett et al. 1998; Henrichs 2005; Jin and Bethke 2009). In oxic environments, the degradation rates of refractory compounds remain high because of the very high oxidative potential and the resulting weak sensitivity towards the depletion of energy-rich organic compounds while in anoxic environments deprived of energy rich-organics and powerful TEAs, the degradation rate becomes thermodynamically limited (Arndt et al. 2013).

Bioenergetic theory: the FT term

Microorganisms channel part of the catabolic energy released by redox reactions into metabolism and growth via intracellular synthesis of adenosine triphosphate, ATP (Regnier et al. 2011). In natural settings, the fraction of energy an organism uses to make ATP is essentially unknown (Harold 1986; Russell and Cook 1995). Therefore, the energetics of catabolic redox reactions only provides a quantitative estimate of the maximum amount of ATP that can be synthesized (LaRowe and Helgeson 2007; LaRowe et al. 2008; LaRowe et al. 2012). ATP synthesis implies that catabolic reactions must generate useable energy under the non-standard state conditions of natural environments (Smeaton and Van Cappellen 2018), that is, ΔGcat < 0:


where ΔGcat and Δcat are the Gibbs energy change of the catabolic reaction under actual and standard-state conditions, respectively, and Qcat refers to the reaction quotient of the reaction. The notion of standard state results from the impossibility to define absolute values of some thermodynamic quantities. Only changes can be determined, which require definition of a baseline, or a standard state, for substances. Precise definitions for gases, pure liquids, solids and admixtures are provided in Cox (1982). For application of the concept of standard state to substances in solutions, which is particularly relevant to this chapter, the composition of the system (as well as the pressure) must be defined and is customarily taken as the standard state molality of 1 mol kg−1. Note that temperature does not intervene in the definition of the standard states, but most thermodynamic tables report values at the recommended temperature of 298.15 K.

Values of Qcat, are calculated using:


where ai designates the activity of the ith species and vi corresponds to the stoichiometric coefficient of the ith species in the given catabolic reaction (LaRowe et al. 2014).

It is thought that redox reactions can however only proceed when the Gibbs energy yield for the catabolic reaction ΔGcat exceeds a minimum metabolic threshold. Such requirement can be added to the thermodynamic limiting term, FT, in the Monod model. The standard formulation (Boudart 1976; Jin and Bethke 2002; Jin and Bethke 2007) relates the minimum excess catabolic energy production to the energy required to synthesize one mole of ATP from ADP and monophosphate, ΔGATP, in such a way that the FT term reads:


In Equation (18), m is the number of ATP molecules produced per catabolic formula reaction, R is the universal gas constant, T is the absolute temperature, and χ is the average stoichiometric number. The latter is equivalent to the number of times the rate limiting step occurs per mole of ATP made multiplied by the number of electrons transferred in the jth reaction (Jin and Bethke 2002; Jin and Bethke 2003; Jin and Bethke 2005). Thermodynamics therefore imposes that ΔGcat+mΔGATP must be negative for a microbially mediated reaction to proceed. Because ΔGATP is defined as a positive (energy-requiring) value, ΔGcat must be sufficiently negative to exceed mΔGATP in absolute magnitude (Regnier et al. 2011). Equations (16, 17) show that the accumulation of reaction products (such as H2 in Eqn. 5) may limit the thermodynamic drive of the catabolic reaction and thus reduce rates in the Monod model via the FT term. For microbial reaction processes operating close to their thermodynamic limit, the reaction rates may become very sensitive to FT (Thullner et al. 2007). This has been shown, for example, in the case of anaerobic oxidation of methane in marine sediments (Regnier et al. 2005; Dale et al. 2006, 2008a) or for methanogenesis and sulfate reduction with H2 in hydrothermal systems (LaRowe et al. 2014). When geomicrobial reactions become thermodynamically limited, most energy generated by the catabolic process is then diverted to maintenance functions, with little energy left to invest in growth.

Alternative expressions accounting for thermodynamic limitations of microbial reactions have been proposed in the literature (Cupples et al. 2004; Thullner et al. 2007). Recently, LaRowe et al. (2012) proposed a formulation for FT that (1) relies on only one adjustable parameter rather than the three required by Equation (18)GATP, m and χ) and (2) circumvent the need to set a priori values of mΔGATP. Here, the approach relies on the amount of energy required to maintain a membrane potential (ΔGmp= FΔΨ) as a proxy for the minimum amount of energy that a microbe needs to be considered active:

FT={1exp(ΔGcat+FΔψRT)+1for ΔGcat00for ΔGcat0

where ΔGcat denotes the Gibbs energy of the catabolic reaction (per electron transferred), F is the Faraday constant and ΔΨ is the electric potential (in volts) across an energy transducing membrane. Such expression has been used to simulate, e.g., microbial dynamics in hydrothermal vent systems (LaRowe et al. 2014) and coupled CH4–SO4 cycles in sediments on an active continental margin offshore New Zealand (Dale et al. 2010) where methane-rich fluids migrate upwards through the sediment (LaRowe et al. 2012). The sensitivity of a reaction rate to the FT term in the above equation is shown for methanogenesis (MET) occurring in the wall of a hydrothermal vent chimney, as an example. Figure 2 shows that the computed MET rates are significantly reduced by the low energy yields in the inner part of the chimney wall, and become 0 towards the outer part of the chimney when the values of ΔGcat are positive. Without the FT term, the active zone of methanogenesis within the chimney wall would be larger.

The dependence of microbial reaction kinetics on bioenergetics is accounted for in the FT term. However, catabolic energy gains also constrain the value of Y (Roden and Jin 2011) and thermodynamics is also needed to relate Y to the energetics of cellular metabolism (Smeaton and Van Cappellen 2018). Therefore, the growth yield, Y, provides another link between the kinetics and thermodynamics of microbially driven redox reactions.

Quantifying growth yields

Growth yields (Y) are generally incorporated into biomass-explicit kinetic models using the Monod formulation, as summarized in section ‘Traditional Approaches for Simulating Biogeochemical Processes in the Subsurface’ (e.g., Thullner et al. 2007). In simple terms, Y values depend on the Gibbs energy generated by the catabolic reaction, the Gibbs energy needed for the formation of new biomass, and the efficiency with which organisms utilize energy (VanBriesen 2002). Two distinct lines of approach that relate Y to the catabolic energy yield either through empirical relationships (e.g., Rittmann and McCarty 2001; Roden and Jin 2011) or through bioenergetically based models (Heijnen and Van Dijken 1992; Heijnen et al. 1992; McCarty 2007) such as the Gibbs Energy Dissipation model, have then been followed. In short, the latter approach accounts for the loss of energy due to entropy production and heat dissipation by cells (see Dale et al. 2006 for further details). As shown by VanBriesen (2002) and for a variety of organic carbon compounds and metabolic pathways (aerobic degradation, denitrification and methanogenesis), energy dissipation approaches are consistent with empirically based methods and predict quite similar yields with laboratory-determined values of Y. A slightly less complex theoretical approach has also recently been developed (the so-called Microbial Turnover to Biomass (MTB) method of Trapp et al. (2018)) and applied to the quantification of pesticide degradation and formation of biogenic non-extractable residues (Brock et al. 2017). Overall, the major drawback of these theoretical calculations is that they provide only maximum Y values, and do not consider energy changes due to variations in chemical composition of the system (Thullner et al. 2007). That is, Y values are calculated under biogeochemical standard state conditions that may deviate significantly from those encountered in natural settings (LaRowe and Amend 2015b).

To circumvent this important limitation, Smeaton and Van Cappellen (2018) recently developed a quantitative method for Y that accounts for changes in physical (e.g., temperature) and chemical conditions under which the metabolic processes are occurring. Their resulting semi-empirical model, the Gibbs Energy Dynamic Yield Method (GEDYM), is an extension of the one by Heijnen et al. (1992) and relies upon a much larger database of experimental Y values, most of which are relevant to low energy yielding catabolic processes such as methanogenesis and sulfate reduction. Briefly, the method computes for a given microorganism, Gibbs energy changes of metabolic reactions by linking Gibbs energy changes of their corresponding catabolic (ΔGcat) and anabolic (ΔGsyn, see below) reactions through their growth yield. In their approach, Gibbs energy changes of metabolic, catabolic and anabolic reactions all depend on Q as shown by, e.g., Equation (16). Therefore, the resulting Y values account for changes in the chemical environment surrounding the cells (Smeaton and Van Cappellen 2018). As an example of application of the GEDYM model, Figure 3 reports Y values for sulfate reduction coupled to acetate oxidation for a broad range of environmental conditions, which in turn impact values of ΔGcat and ΔGsyn. Y values applied in geochemical models generally fall within their predicted range but are assumed constant (Smeaton and Van Cappellen (2018) and references therein).

At this stage, it is important to note that the metabolic energy balance reported above only accounts for the catabolic energy invested in anabolism, that is, biomass synthesis. It therefore differs from the total cellular energy balance, which also includes the energy needed to sustain all maintenance functions. In a biomass-explicit kinetic model, maintenance energy requirements are implemented separately through an additional term, as shown in Equations (12, 13). The following section elaborates on the quantification of anabolic and maintenance energy requirements through classical thermodynamics.

Energetics of anabolism and maintenance

The calculation of the Gibbs energy required to synthesize biomolecules, ΔGsyn, from simple material such as CO2 or acetate, involves three important steps: 1) definition of cellular biomass composition; 2) definition of the anabolic reaction stoichiometry; and 3) computation of the resulting Gibbs energy of the anabolic reaction.

In RTM applications, biomass synthesis generally assumes a unique biomass composition such as C5H7O2N (Rittmann and McCarty 2001; Dale et al. 2006, 2008a) or C5H9O2.5N which are based on the measured elemental concentrations of microbial species (Roels 1980; Smeaton and Van Cappellen 2018). A more complex approach consists in explicitly accounting for the variety of molecules making up the cell, that is, individual amino acids, nucleotides, lipids, saccharides, amines and other compounds (McCollom and Amend 2005) and subsequently calculate the energetics associated to each of these biomolecules. The polymerization of these various compounds in their respective biomacromolecules can also be accounted for in the energetics of anabolism, as performed in, e.g., Amend et al. (2013) and LaRowe and Amend (2016). To our knowledge, such approach as not yet been applied in a reactive-transport framework, despite being in principle compatible with RTM simulations.

Methods to constrain the anabolic reaction stoichiometry are presented in detail and illustrated in Rittmann and McCarty (2001), Dale et al. (2006), and Smeaton and Van Cappellen (2018). Briefly, establishing the full anabolic reaction requires combination of the half-redox reaction for biomass synthesis and the half-redox reaction providing the carbon source and electrons. If the carbon in microbial biomass is less oxidized than the carbon source, then the latter can also serve as electron donor in the full anabolic reaction. If not, the carbon source cannot be used to balance electrons during anabolism and an electron donor distinct from the carbon source is required in this case. Following LaRowe and Van Cappellen (2011), Smeaton and Van Cappellen (2018) have shown that the relative oxidation states of biomass (ROSB= 4 − NOSCB) and carbon source (ROSS= 4NOSCS) can directly be compared using the NOSC of both compounds. Thus, only when ROSS< ROSB the carbon source can also serve as electron donor.

As an example, consider biomass synthesis using ethanol as carbon substrate:


Here, NOSCS and ROSS of ethanol are −2 and 6, respectively while NOSCB and ROSB of biomass are 0 and 4, respectively. Therefore, since ROSS > ROSB, an electron acceptor other than the C source is needed to consume the electrons released by anabolism. For instance, the two electrons produced in Equation (20) can be consumed during sulfate reduction:

0.25 SO42+2.25 H++ 2.0e0.25 HS+H2O

which leads to the full anabolic reaction:

0.5 C2H6O+0.25 SO42+0.2 NH4++0.05 H+0.2 C5H7O2N+0.25 HS+1.1 H2O

Importantly, the oxidation states of other elementary building blocks (e.g., N, S…) required for biomass synthesis also influence the stoichiometry of the full anabolic reaction (LaRowe and Amend 2016). For instance, in well-oxidized conditions NO3 will be used as nitrogen source during biosynthesis while under reducing condition, NH4+ will be used instead. In the latter case, the microorganisms will not need to invest a stream of electrons to reduce the nitrogen compounds during biosynthesis (LaRowe and Amend 2016). This example highlights that the prevailing environmental conditions not only affect the energetics of catabolism but may also influence the energetic of anabolism, an aspect not yet fully accounted for in RTM applications.

The full anabolic reaction formula derived above can then be used to compute the energetics of biomass synthesis. This step is relatively straightforward if a single biomass composition and standard state conditions are assumed, as shown by the RTM simulations of, e.g., Dale et al. (2008a). However, as already pointed out, similar computation using a range of biomacromolecular compounds synthesized under variable environmental conditions that depart from standard states has, to our knowledge, not yet been performed in the context of reactive-transport modeling. In principle, such coupling would nevertheless be achievable and ΔGsyn values would then be calculated according to LaRowe and Amend (2016):


where ΔGr,i stands for the Gibbs energy of the reaction describing the synthesis of the ith biomolecule and ΔGpoly denotes the Gibbs energy required to polymerize biomolecules into their respective biomacromolecules, a term neglected when simple non-polymerized biomass composition are assumed. Values of ΔGr,i and ΔGpoly can then be calculated as a function of temperature and pressure using equations similar to Equations (16, 17).

The results by LaRowe and Amend (2016) clearly demonstrate that values of ΔGsyn depend strongly on the combination of the redox state of the precursor compounds and that of the environment, with values varying by up to about 40 kJ (g cell)−1 depending on the different combinations. LaRowe and Amend also show that the contribution of ΔGpoly to ΔGsyn is significant when biomass synthesis occurs in the most conducive environmental conditions. Overall, the important control of environmental conditions on the energetics of anabolism calls for more direct linkages between spatio-temporal gradients in redox conditions simulated by RTMs and their effects on biomass synthesis.

In addition to the energy required to synthesize biomolecules, the total cellular energy balance also includes the maintenance energy ΔGmain which is needed to sustain all other functions in support of viability and that do not result in new biomass (LaRowe and Amend 2016). The total cellular energy requirement, ΔGcell, thus reads:


where we note that the above formulation is certainly a simplification because energies required for biosynthesis and maintenance are not completely decoupled (van Bodegom 2007).

Thermodynamic models for the maintenance energy compute ΔGmain from typical doubling times of microbial populations and the power demand required to maintain cell integrity (LaRowe and Amend (2016) and references therein). The latter must nevertheless distinguish the power used while organisms are growing, from the so-called ‘basal maintenance power’ (Hoehler and Jørgensen 2013), which is the power required by microorganisms to remain viable only and is typically several orders of magnitude lower (LaRowe and Amend 2015a). Overall, maintenance energies can become important for populations that have long doubling/replacement times and large (cumulative) maintenance powers, ΔGmain exceeding ΔGsyn in this case, and may even become the dominant component of the total energy budget. Major avenues for future RTM research would thus be to calculate changes in power demands for maintenance as a function of environmental conditions and evolving states of the microorganisms as well as to better couple the energetics of biomolecule synthesis and maintenance into a unified modeling framework.

Carbon Use Efficiency as an alternative to Ys

Before proceeding, it is important to stress that thermodynamic constraints on microbial growth have been measured, reported and implemented in models using multiple currencies (Sinsabaugh et al. 2013). In particular, in the field of ecological modeling and especially in soil science, Y is typically not used and growth yields are constrained from rates of carbon transformation using the concepts of ‘Carbon Use Efficiency’ (CUE) and its equivalent Microbial Growth Efficiencies, MGE. The CUE is generally defined as the ratio of growth (µ) to assimilation, that is, CUE = µ /(µ + R), where R includes any C losses to respiration (Sterner and Elser 2002; Manzoni et al. 2012b). Reported CUE values determined in the laboratory and in the field have been synthesized by del Giorgio and Cole (1998) for aquatic systems and by Six et al. (2006) for soil environments. These empirically determined CUEs span a broad range, from typical values as low as 0.01 up to values close to 0.8, partly reflecting the diversity of processes captured by the applied empirical methods (physiological, but also characteristic of community or ecosystem dynamics) that influence C metabolism across varying spatio-temporal scales (Manzoni et al. 2012b; Sinsabaugh et al. 2013; Geyer et al. 2016). In order to better organize and properly interpret CUEs, Geyer et al. (2016) recently proposed a conceptual framework that structures its definition according to increasingly temporal and spatial scales of investigation. In short, Geyer and co-workers distinguish the CUE of populations, which is governed by species-specific metabolic and thermodynamic constraints as described by the growth yield Ys discussed above, from CUEs of communities and ecosystems which include many additional controls such as substrate stoichiometry, external physico-chemical conditions (temperature, moisture, pH,…) or substrate recycling. These distinctions are important because many large scale models of soil carbon dynamics, especially those embedded in Earth System Models of the coupled biogeochemical cycles and climate such as the CENTURY model (Parton et al. 1988) rely on the CUE (or the MGE) concept, and mostly address the larger ecosystem scale (e.g., Allison (2014), Wieder et al. (2014), Huang et al. (2018)). The reader is referred to e.g., Geyer et al. (2016) for further discussion on the use of the CUE and its relation to Y across scales.


In recent years the improvement of biomolecular methods has led to a number of powerful approaches providing sequencing data from e.g., genomics, proteomics, transcriptomics or metabolomics and the combination of them. Such ‘omics’ approaches provide information on the metabolic potential of individual microbial species as well as of microbial communities (Müller and Hiller 2013; Franzosa et al. 2015). Together with biochemical data this has led to the set-up of extensive biochemical reaction networks describing the entire metabolism of a cell (Fig. 4). The resulting ‘genome-scale model’ provides the basis for the prediction of cellular functions such as growth or the formation of specific metabolites (Durot et al. 2009; Kim et al. 2012; Monk and Palsson 2014; O'Brien et al. 2015). For this purpose the stoichiometric matrix describing the reaction network needs to be combined with a ‘constrained-based modeling’ approach making use of additional knowledge and assumptions e.g., regarding the occurrence and magnitude of individual reaction fluxes and the steady-state assumption (Bordbar et al. 2014) and of the available experimental data e.g., on gene expression, protein expression or metabolite concentration (Reed 2012). Among a variety of such approaches the so-called ‘flux-balance analysis’ is one of the most commonly used (Orth et al. 2010). This methodology allows considering dynamic changes of the chemical environment of the cells (Mahadevan et al. 2002), but several other methods have also been proposed (Lewis et al. 2012; Ramkrishna and Song 2012; Song et al. 2014). While ‘omics’-based approaches have initially been used for modeling single microorganisms (or cell types) they have recently been expanded to describe metabolic fluxes in bulk microbial communities (Larsen et al. 2011; Biggs et al. 2015; Hanemaaijer et al. 2015; Gottstein et al. 2016; Perez-Garcia et al. 2016). The size of such communities is constrained by the required input data as well as the computational demands. So far, applications range from two-species systems (Klitgord and Segre 2010; Zomorrodi and Maranas 2012) to several hundreds of species (Magnúsdóttir et al. 2017).

Most metabolic flux approaches find their origin in the areas of medical and biotechnological research and have shown their potential for predicting the response of single cells to genetic modification or to drugs, or have been used to optimize the microbial production of specific compounds in industrial processes. However, they have also been used for the simulation of specific biochemical processes of relevance to the geosciences. Examples are methane oxidation by Methylococcus capsulatus (Lieven et al. 2018), methanogenesis by a community of Desulfovibrio vulgaris and Methanococcus maripaludis (Stolyar et al. 2007), and Fe(III) reduction by Geobacter metallireducens (Fang et al. 2012), by a community of Geobacter sulfurreducens and Rhodoferax ferrireducens (Zhuang et al. 2011) or by a community of the latter two species and Shewanella oneidensis (Zomorrodi et al. 2014). The Fe(III) reduction process is of particular interest in the context of uranium (U) reduction which is performed co-metabolically by some Fe(III) reducing microbes and leads to immobilization of uranium in contaminated aquifers.

The few applications of genome-scale descriptions of microorganisms in a RTM framework are also focusing on Fe-U cycling. King et al. (2009) coupled the dynamics of Geobacter sulfurreducens as described by Mahadevan et al. (2006) with a two-dimensional reactive transport solver and compared the resulting microbial growths and substrate consumptions with those obtained with more simplified description of microbial dynamics. Similarly Scheibe et al. (2009) coupled the same description of Geobacter sulfurreducens growth to the reactive transport model HYDROGEOCHEM to simulate acetate induced reduction of Fe(III) and U(VI) along a one-dimensional representation of an aquifer. In a follow-up study by Fang et al. (2011), a more integrated coupling of genome-scale data and reactive transport modeling was achieved. The same type of coupling was recently presented by Tartakovsky et al. (2013) using a “pore-scale smooth-particle hydrodynamics approach” (Tartakovsky et al. 2009). Their results were also compared to those obtained with a Monod-type model of bacterial growth and activity.

With the continuously increasing availability of sequencing data and computational power the inclusion of genome-scale models into reactive transport approaches has a clearly growing potential for applications focusing on the study of in situ processes occurring in porous media. There are however several important challenges limiting such applications. Natural microbial communities often consist of an unaccounted number of species and depending on the processes of interest the amount of data needed to set up a genome-scale model might be excessively large. An automated assembly of such data is supported by metabolic data bases (Karp et al. 2019) but results need to be checked for inconsistencies (Richter et al. 2015). In addition, a community-wide metabolic network implicitly considers all metabolites produced by a species to be fully available to all other species in the system, an approach which neglects limitations in availability due to transport into and out of a cell as well as to spatial separation of the cells. These limitations could be circumvented by considering the spatial distribution of the cells but existing approaches are limited to small number of species (Harcombe et al. 2014). Furthermore, substrate concentrations and environmental conditions may vary down to the pore-scale (Semple et al. 2004; Johnsen et al. 2005; Hesse et al. 2009; Schmidt et al. 2018), which challenges the genome-scale simulation of microbial dynamics in the same way as it does for other kinetic models of microbial activity.

Next to genome-scale models other less complex approaches that take available sequencing data into account exist. In these simpler methodologies, the data are used as proxies or biomarkers to constrain the dynamic behavior of a microbial community. While in the past, traditional microbial biomarkers such as cell counts, colony forming units, proteins or polysaccharides were used already as experimental references for the dynamics of microbial species or functional groups in the context of reactive transport simulations (Wirtz 2003; Thullner et al. 2004) more recent approaches combined these traditional approaches with marker genes to simulate the response of microorganism to fluctuating growth conditions (Stolpovsky et al. 2011) or carbon turnover and pesticide degradation in in a 1-D soil column (Pagel et al. 2016). Reed et al. (2014) and Louca et al. (2016) went further and applied a complex ‘gene-centric’ approach using marker genes as experimental references, the production of which was linked to biogeochemical reactions in the ocean.


The main reason for the simulation of microbial dynamics in the context of reactive transport modeling is their ability to catalyze biogeochemical reactions. This requires first of all modeling approaches describing the growth and metabolic activity of microbial species appropriately. However, the simulated systems represent an entire microbial ecosystem of varying complexity. In such an ecosystem the abundance of the microorganisms of interest can be affected not only by the concentration of substrates but also by a number of additional constraints, and the ability of the microorganisms to respond to these constraints also determines their abundance and their functional performance. There are numerous types of interactions between a (group of) microorganisms and their environment which for the sake of simplicity are separated here into two groups: microbe–microbe interactions, and the response of the microorganisms to abiotic constraints.

Microbe–microbe interactions

(Microbial) species may interact in different ways, with one species promoting or inhibiting the growth and activity of another species (Fig. 5, Lidicker 1979). While it can be challenging to determine such interactions in large communities (Faust and Raes 2012) specific interactions can be modeled using a generalized Lotka–Voltera approach describing the change of abundance Bi of a species i interacting with N other species according to:


with µi as specific growth rate and γi,j as strength of the interaction between species i and j. γ values can be either positive or negative depending on the type of interaction, a value of 0 indicating no influence of species j on species i (Fig. 5). Applications of the Lotka–Voltera model in the context of microbial community dynamics range from single predator-prey interactions (Mauclaire et al. 2003) for which this approach was originally developed (Lotka 1925; Voltera 1926) to interactions in large communities found in the human gut (Stein et al. 2013; Kuntal et al. 2019) or in cheese (Mounier et al. 2008). Next to direct interactions such as protists predation (Mauclaire et al. 2003) or phage infections (Jover et al. 2013), indirect interactions have also been addressed by the generalized Lotka–Voltera approach. Indirect interactions may be mediated by a chemical compound produced or provided by a microbial species, which then promotes or inhibits the growth and activity of another species. Such interactions can be modeled by including a dependency of growth rates on the abundance of other species as described above but can of course also be addressed explicitly, e.g by linking the rate expressions describing the two species to the concentration of such mediator compounds or via a multi-species metabolic network simulation (Freilich et al. 2011). Similarly, the competition of species for a chemical compound can be addressed by a generalized Lotka–Voltera approach but can also be explicitly implemented in reactive transport simulations for, e.g., describing redox stratification associated with organic carbon degradation (Thullner et al. 2007).

A different type of microbe-microbe interaction occurs in the terrestrial mycosphere where bacteria can benefit from the presence of extensive networks of fungal mycelia (Harms et al. 2011; Worrich et al. 2018). These networks can act as ‘fungal highways’ which promote the mobility of bacterial cells (Kohlmeier et al. 2005) (see below) and as ‘fungal pipelines’ that facilitate the transport of water and nutrients to bacteria (Furuno et al. 2012; Schamfuss et al. 2013; Worrich et al. 2017). These effects have been simulated by considering a highly increased diffusion coefficient for transport of bacteria and/or nutrients along the fungal networks (Banitz et al. 2013).

Interactions between microbes and their physical environment

Subsurface environments impose a multitude of physical constraints on the abundance and activity of microorganisms. The solid matrix and its properties define the space and surface area available for the microorganisms while water flow (and distribution in unsaturated media) controls not only the transport of dissolved compounds but also the relocation of microorganisms. Due to the heterogeneity of the subsurface, these physical constraints can exhibit spatial variations at scales ranging from the pore scale to the field scale. Furthermore, flow and transport conditions can vary in time due to surface processes such as precipitation events, seasonal effects or anthropogenic activities. As a result, the microbial communities that live in the shallow subsurface need to adapt to changing environmental conditions. In particular such changes can involve periods of unfavorable conditions or disturbances (Fig. 6) (Allison and Martiny 2008; Shade et al. 2012) which can occur either as extended stress periods (Manzoni et al. 2012a; Worrich et al. 2016; Rocca et al. 2019) or as (series of) short disturbance events.

Microbial transport of bacteria by the flowing water is generally described by advective–dispersive transport combined with rate expressions for attachment and detachment of cells (Ginn et al. 2002; Tufenkji 2007). Attachment and detachment of cells are commonly simulated using first-order rate expressions with values of the rate parameters depending on the flow velocity and the resulting shear forces, the size and surface properties of the bacterial cells as well as the structure and surface properties of the solid matrix. Considering cells as solid particles allows application of the filtration theory (Tien and Payatakes 1979) to predict the attachment rate coefficient katt:


with n as porosity, d as average grain size, v as transport velocity, η as collection efficiency (with values predicted from theory), and α as (empirically determined) collision efficiency. Limitations of the theoretical concepts of filtration theory to simulate the attachment process have been attributed e.g., to heterogeneities of the adhesion properties within the bacterial population (Simoni et al. 1998), to different microbial species affecting each other in their deposition behavior (Stumpp et al. 2011), or to repulsive interactions limiting the adhesion of cells on the matrix. While adhesion can alternatively be predicted by the (extended) DLVO theory (Hermansson 1999), there are also limitations to its application (Tufenkji 2007; Boks et al. 2008). Detachment rates depend on various factors (Peyton et al. 1995; Xavier et al. 2005) and are difficult to predict also as cells may increase their resistance to detachment producing extracellular polymeric substances (EPS) (Tay et al. 2001) or may in turn actively detach from surfaces as response to nutrient availability or other favorable/unfavorable conditions (see Ginn et al. 2002 and literature cited therein). A general theory (such as filtration theory for attachment) providing a quantitative description of detachment processes has not been introduced, yet.

Besides being passively transported by the water flow or by the movement of solid particles they are attached to (Thullner et al. 2005) some bacterial species are motile and are thus actively able to swim in water (Blair 1995) or slide along a surface (Kearns 2010). This active movement may occur as random or chemotactic movement, that is to say, by sensing and following chemical concentration gradients (Berg 2000; Alexandre et al. 2004). This movement allows bacterial cells to relocate to more favorable conditions e.g., by moving towards higher nutrient concentrations or by moving towards lower concentrations of hazardous compounds—an ability which is however at the expense of significant metabolic costs. Furthermore, bacteria can produce chemoattractants triggering the chemotactic self-attraction of their own species and promoting their aggregation (Mittal et al. 2003; Park et al. 2003). In addition, bacterial motility may affect the attachment rates of cells in porous media (Nelson and Ginn 2001; de Kerchove and Elimelech 2008), their resilience towards disturbances (König et al. 2017) and their macroscopic transport behavior (Ford and Harvey 2007; Bai et al. 2016; Creppy et al. 2019), and it may promote the formation of pronounced spatial distribution patterns (Budrene and Berg 1991; Keymer et al. 2006). To model bacterial motility, the time evolution of bacterial abundance B in space can be described as


with Db as random motility coefficient and χs,i as chemotactic sensitivities to the gradient of chemical compounds si. Sensitivities may be functions of the chemoattractant concentration (Keller and Segel 1971; Ford and Harvey 2007) or more simply considered as constants (Centler et al. 2011). In any case, positive χs,i values indicate attraction, while negative χs,i values indicate repulsion. Alternatively, the explicit link to the gradient of a chemical compound can be replaced by a dependency of the random motility coefficient Db on the concentration of a chemical compound (Banitz et al. 2011a). Typically only a single compound (e.g., a nutrient) is considered to drive the chemotactic motility of bacteria, but two compounds might be considered if, in addition to a substrate, a bacterially emitted chemoattractant is also present (Saragosti et al. 2010; Centler et al. 2011).

Implementation of these motility concepts into reactive transport models have allowed to match high resolution data on bacterial chemotaxis (Pedit et al. 2002; Banitz et al. 2012) and predict traveling bands of bacteria (Hilpert 2005; Saragosti et al. 2010) as well as to simulate the formation of aggregated population patterns (Centler et al. 2011; Gharasoo et al. 2014; Centler and Thullner 2015) or the spreading of chemotactic bacteria at larger scales (Valdés-Parada et al. 2009). Microorganisms can also promote the dispersal of other microbial species (Ben-Jacob et al. 2016). Assuming fungal hyphae as pathway of increased bacterial motility (‘fungal highways’) in the simulation of microbial systems allowed describing the benefit of such high motility networks for biodegradation (Banitz et al. 2011a,b), for horizontal gene transfer (Berthold et al. 2016) or for the resistance of bacterial populations to disturbances (König et al. 2018a,b). A link between such reactive transport approaches with models describing fungal growth in porous media (Cazelles et al. 2013) has not been introduced, yet.

Dormancy is another microbial strategy to endure stress periods caused by natural effects (Lennon and Jones 2011; Joergensen and Wichern 2018) or anthropogenic perturbations (Balaban et al. 2004). When facing such unfavorable conditions microorganisms can switch from an active state into a dormant state which is characterized by a reduced metabolic activity, lower maintenance requirements and thus a better survival of the cells. If environmental conditions become favorable again microorganisms can reactivate and grow again (Dworkin and Shah 2010). For modeling microbial dormancy microorganisms are classically subdivided into an active and a dormant fraction (Fig. 7), but approaches considering a single fraction with transient changes of its dormancy degree exist, too (Resat et al. 2012). The transition between these two fractions is typically linked to environmental conditions by using a kinetic rate expression (depending on environmental variables) for the deactivation rate (transition into dormancy) and a complementary expression for the reactivation or resuscitation (transition from dormancy to active state) rate (Bär et al. 2002; Jones and Lennon 2010; Wang et al. 2014). Alternatively, a switch function that depends on environmental conditions can be used to determine the direction of the transition (Stolpovsky et al. 2011; Mellage et al. 2015; Bradley et al. 2018a). Other concepts consider both, deactivation (Chihara et al. 2015) and reactivation (Epstein 2009; Buerger et al. 2012) to be random processes, or rely on a growth rate dependent dormancy index to determine the fraction of dormant bacteria (Wirtz 2003). Furthermore, different degrees of dormancy may be considered assuming either distinct subgroups of dormant microorganisms or considering the depth of dormancy to increase with the duration of the unfavorable conditions (Stolpovsky et al. 2011). Models using such concepts have been applied to match data from specific experiments (Stolpovsky et al. 2011; Wang et al. 2014, 2015) and to analyze environmental samples ranging from dynamic systems such as the surface layer of lakes (Jones and Lennon 2010) to low activity systems like the deep subsurface of marine sediments or to soil organic matter composition at arbitrarily large scales (Huang et al. 2018). Furthermore, such simulations highlighted the relevance of dormancy for microbial competition (Stolpovsky et al. 2016) and long-term survival (Bär et al. 2002) in periodically changing environments. When implemented into RTMs dormancy approaches have allowed to match observations from shallow marine sediments (Wirtz 2003) and to study the influence of dormancy on the coexistence of competing species in heterogeneous porous media (Stolpovsky et al. 2012).

Noteworthy, yet not addressed this chapter, are feedbacks of the microorganisms on their physical environments. This includes changes of the hydraulic properties of the solid matrix (Hommel et al. 2018) due to biomass aggregation (Baveye et al. 1998; Yarwood et al. 2006; Thullner 2010), or due to microbially induced mineral dissolution, precipitation (Barkouki et al. 2011; Ebigbo et al. 2012) or gas formation (Mahabadi et al. 2018).


Most elements, e.g., H, C, N, S, O associated in biogeochemically reactive compounds are found in more than one stable isotope form and the relative abundance of these isotopes can be determined with high accuracy. In recent years, isotopic methods have become increasingly used for the assessment of microbially driven processes in the subsurface. The most common approach is to use compound-specific stable isotope analysis (CSIA) to detect the stable isotope fractionation of reactants during (biogeo)chemical transformations (Brüchert 2004; Meckenstock et al. 2004; Elsner 2010; Thullner et al. 2012; Brunner et al. 2013; Blaser and Conrad 2016). The consideration of these fractionation effects in the context of reactive transport modeling is addressed in Druhan et al. (2019, this volume) and will thus not be discussed here. Instead, focus of this section is on the use of stable isotope signatures to trace the incorporation of inorganic compounds into microbial biomass, to identify the origin of microbial degradation products and to follow microbial degradation pathways.

Stable isotopes as tracers of matter fluxes into biomass

Since decades, the isotope signatures of animal biomass is known to be determined by the isotope signatures of their food sources (DeNiro and Epstein 1978). In the same way the isotopic composition of microbial biomass reflects the isotopic signatures of the used substrates (Coffin et al. 1989). Furthermore, stable isotope probing (SIP) approaches relying on the isotopic composition of specific biomarkers allow determining which (group of) microbial species pertaining to a complex community have taken up a specific substrate (Boschker and Middelburg 2002; Abraham 2014; Vogt et al. 2016b). Originally lipids were the commonly used biomarkers (Pancost and Sinninghe Damsté 2003; Yao et al. 2014; Wegener et al. 2016), which together with isotopically labeled substrates allowed for the detection of the substrate uptake into biomarkers (Pelz et al. 2001; Pombo et al. 2002; Stelzer et al. 2006). More recently, such SIP approaches were extended to the analysis of e.g., microbial DNA (Coyotzi et al. 2016), RNA (Lueders et al. 2016; Bradford et al. 2018) or proteins (Jehmlich et al. 2016) used as biomarkers for specific (groups of) microorganisms and the uptake of labeled substrate into their biomass.

The inclusion of stable isotope signatures into RTMs is well established. In short, heavy (h) and light isotope (l) fraction of a specific element in a reactive compound i are treated independently as two separate compounds with concentration hci and lci, respectively. The isotopic signature for each compound is then given by the concentration ratio Ri = hci / lci and commonly expressed as a relative deviation δi from a standard value (Rs):


Note that δ-values are typically expressed in ‰. Analogously, for each reactive transformation rate j expressions hrj and lrj have to be defined for each fraction. The rate expressions have to fulfill the relation hrj / lrj = αjRi if compound i is consumed by process j. The fractionation factor αj describes the magnitude of the stable isotope fractionation induced by the reactive transformation with αj = 1 indicating no fractionation. For further details—in particular regarding the simulation of stable isotope fractionation—see e.g., van Breukelen and Griffioen (2004), Thullner et al. (2008, 2012) and Druhan et al. (2019, this issue).

Although this RTM concept is well suited to describe the uptake of substrates and the resulting stable isotope signature into microbial biomass it has only rarely been applied for this purpose yet. Mauclaire et al. (2003) used stable isotope signatures to simulate carbon fluxes along a microbial food chain in a batch system (Fig. 8). 13C-labeled toluene was used as primary substrate, which was taken up by a bacterial species that served as a prey for a protist predator. Isotopic signatures of the bacterial and protist biomasses were determined using biomarker fatty acids while total bacterial biomass was quantified using cell counting combined with image analysis techniques. Rate expressions based on Monod-type growth of the bacteria and Lotka–Voltera type predation for the protist were used to simulate changes of concentration and carbon isotope signatures in the different carbon pools (toluene, bacterial biomass, protist biomass). Using this approach, a carbon conversion efficiency of approximately 30% for the bacterial growth on toluene, and of approximately 10% for the protist predation on the bacteria was determined. Relying on 13C-labeled E. coli cells added to a soil batch system Kindler et al. (2009) quantified the transformation of the labeled carbon during recycling of the biomass in the microbial food chain. Fatty acids were used to distinguish between living and dead biomass, the latter considered as part of the soil organic matter. Instead of using a classical kinetic model approach the process was modeled by applying a series of feeding cycles to quantify the flux of carbon between living biomass, soil organic matter and inorganic carbon (CO2). Results indicated that within each feeding cycle approximately 20% of the biomass food source was transferred into soil organic matter and that 33% of the added carbon source accumulated as soil organic matter. Carbon uptake by microbes and the associated changes of carbon isotope signatures in the biomass was simulated by Alperin and Hoehler (2009a) for anaerobic methane oxidation in marine sediments. To do so, they represented the network of biochemical reactions catalyzed by a consortium of methane oxidizing archea and sulfate reducing bacteria, and by organic matter oxidizing bacteria. In addition, various stable isotope fractionation effects were considered. Using this complex reaction network the authors showed that although the stable isotope signature of the biomass of the methane oxidizing consortia was similar to the stable isotope signature of methane, methane was not the only source of the biomass carbon, and the CO2 released by the organic carbon oxidizing bacteria was another likely contributor.

Stable isotope signatures have also been used as tracers for the formation of reaction products during (microbially induced) reactive transformations (Avery et al. 1999; Bolliger et al. 1999; Kindler et al. 2009). The interpretation of such signatures can be done by mass balance approaches not requiring more elaborate model simulations as long as fractionation effects can be neglected.

Stable isotopes as indicators of degradation pathways

Next to using stable isotope fractionation as qualitative and quantitative indicators of biodegradation processes (not discussed in this chapter), the combined analysis of the fractionation of two elements provides additional information on the dominant fractionation pathway (Elsner 2010; Vogt et al. 2016a)—also known as dual or 2D isotope analysis. The basic principle of this approach takes advantage of the fact that during a reaction a stable isotope fractionation of different elements in the reactant can occur, each of which described by a pathway- and element-specific fractionation factor α. When analyzing the change of the stable isotope signatures of two of these elements a linear relation between these changes (expressed in δ-notation) is often observed as long as fractionation effects are not too strong. The ratio between these changes, that is the slope of the linear correlation when plotting the changes against each other, is specific for each degradation pathway and given by the ratio between the isotope enrichment factors ε = 1 − α (Fig. 9). Therefore, such plots can be used to identify a dominating degradation pathway, e.g., to distinguish between aerobic and anaerobic pathways (Fischer et al. 2007) or between biotic and abiotic transformations (Badin et al. 2016), or to elucidate the specific reaction mechanism (Meyer et al. 2009; Dorer et al. 2014). If fractionation effects are very strong as e.g., for H fractionation or if two pathways take place simultaneously (Van Breukelen 2007; Centler et al. 2013) a more complex analysis is needed.


While the sections above focused on approaches describing the dynamics of microbial growth and activity two aspects have not been directly addressed: the measure by which the microbial abundance is quantified and the spatial representation of the microorganisms. Both of these interlinked aspects are explicitly or implicitly affected by the spatial resolution of the applied modeling approach.

Most reactive transport applications considering microbial abundance use a population-based concept in which biomass is quantified as concentration of microbial carbon mass, microbial cells or of specific parts (DNA, RNA, fatty acids or other biomarkers) of the microbial cells. Concentration may refer to individual species as sometimes implemented for simulating controlled laboratory experiments (Bauer et al. 2009; Dechesne et al. 2010; Banitz et al. 2011a; Monga et al. 2014) while for natural environments, the vast number of species in microbial communities often preclude the possibility to link the rate of a specific biogeochemical reaction to the abundance and activity of a restricted number of individual species. Instead the species are pooled into groups of specific traits, microbial guilds or functional groups e.g., catalysts of specific degradation pathways (Schäfer et al. 1998a; Brun et al. 2002; Wirtz 2003; Thullner et al. 2005; Yabusaki et al. 2011; Bouskill et al. 2012; Pagel et al. 2014; Yu and Zhuang 2019). In contrast, individual-based models consider the dynamic behavior of each single microbial cell in a system (Ferrer et al. 2008; Kreft et al. 2013; Hellweger et al. 2016; Jayathilake et al. 2017; Leveau et al. 2018), an approach that allows considering random variations of microbial behavior or intra-species variations of cell properties. Such individual-based models have been considered also within reactive transport approaches (Gras et al. 2011; Ebrahimi and Or 2014; Gharasoo et al. 2014; Centler and Thullner 2015; Ebrahimi and Or 2015; Kim and Or 2016) but applications are necessarily limited to small-scale (cm-scale or below) systems given typical abundances of 106 cells per cm3 or more in most porous medium environments.

In most aqueous systems microorganisms tend to be evenly distributed in the water without exhibiting significant biomass concentration gradients at small spatial scales. In porous media microorganisms are primarily associated with the solid matrix (Lehman et al. 2001; Griebler et al. 2002; Mellage et al. 2015), which favors the development of small-scale spatial heterogeneities. Even along the surface of the solid matrix large variations in microbial abundance may occur (Dechesne et al. 2003; Nunan et al. 2003; Or et al. 2007; Iltis et al. 2011). As a consequence not all microorganisms at a macroscopic location (defined e.g., by the experimental sampling volume) are exposed to the same environmental conditions. In particular, small scale variations in substrate concentration may expose microorganisms to environmental conditions that differ from the mean state, that is to say, substrate bioavailability is not the same everywhere (Semple et al. 2004). The bioavailability of a substrate is thus no only limited by its interaction with solid and liquid phases (Johnsen et al. 2005; Haws et al. 2006; Amos et al. 2007) but also by a mass transfer process that is required to provide the substrate at the microscopic location of the microbial cells (Bosma et al. 1997). In most cases a limited bioavailability of a substrate leads to reduced degradation rates but in cases where high substrate concentrations become toxic to the degrading microorganisms, limited bioavailability may exceptionally promote higher degradation rates (Hanzel et al. 2012; Gharasoo et al. 2015).

Mass transfer processes can be explicitly included into reactive transport models of high spatial resolution. Such implementation requires a theoretical framework for simulating the distribution of the microbial biomass. The simpler approach represents the interface between pore water and solid matrix as a reactive surface of no volume homogeneously covered by microorganisms (Hesse et al. 2009; Gharasoo et al. 2012). A slightly more sophisticated approach is to consider that microorganisms form a biofilm covering the surface of the solid matrix. Biofilms are aggregations of microbial cells and extracellular polymeric substances (Costerton 1999; Watnick and Kolter 2000; Branda et al. 2005; Carrel et al. 2018), either forming a homogeneous layer of biomass or heterogeneously distributed biomass aggregates on the solid matrix. Since the introduction of the first biofilm models (Rittmann and McCarty 1980) increasingly complex approaches have been implemented (Kreft et al. 2001; van Loosdrecht et al. 2002; Picioreanu et al. 2004; Alpkvist et al. 2006; Cumsille et al. 2014) and combined with different pore-scale reactive transport models (Suchomel et al. 1998; Dupin et al. 2001; Thullner et al. 2002; Kapellos et al. 2007; Thullner and Baveye 2008; von der Schulenburg et al. 2009; Bottero et al. 2013; Qin and Hassanizadeh 2015; Tang and Liu 2017). These simulations have allowed constraining the distribution and thickness of the biofilm as a result of the mass transfer of substrates to and into the biofilm and its subsequent in-situ degradation by microorganisms. Results from pore-scale simulations of mass-transfer limited microbial degradation in porous media (see Fig. 10 for an example) show that the extent of the mass transfer limitation depends (among several other factors) on the size of the pores, on the heterogeneity of the size distribution and on the pore-scale distribution of the biomass along the solid matrix surface (Hesse et al. 2009; Gharasoo et al. 2012; Schmidt et al. 2018).

At the continuum scale, i.e. the scale at which porous media are described by representative elementary volumes (REV), pore scale structures and the pore-scale distribution of the biomass can't be resolved anymore (Fig. 11), and although often still called biofilm, model approaches typically consider only a bulk biophase of unknown spatial distribution (Schäfer et al. 1998b; Schäfer 2001; Ebigbo et al. 2010; Bozorg et al. 2011). Regardless of the distribution type (colonies, aggregates, biofilm etc.), a linear exchange term


has been suggested to describe the mass flux Jtr between the bulk aqueous phase and the biophase (Baveye and Valocchi 1989). This flux term links the bulk substrate concentration cbulk (affected by the transport processes along the model domain) and the substrate concentration cbio inside the biophase, the latter describing the bioavailable concentration. All other factors controlling the magnitude of Jtr are subsumed into the value of the mass-transfer rate coefficient ktr. Predictions of the value of ktr can be obtained from formal upscaling approaches (Hesse et al. 2010; Orgogozo et al. 2013), which do however require a series of simplifying assumptions. As a result, in many reactive transport applications the transfer rate is not well constrained in quantitative terms. In cases of a single rate-limiting chemical compound mass transfer (Eqn. 29) and a biodegradation rate following Michaelis–Menten kinetics (Eqn. 1), both rates can be combined into a single expression known as the Best equation (Best 1955), which has been shown to provide adequate descriptions of mass transfer limited microbial degradation rates (Bosma et al. 1997; Simoni et al. 2001). Alternatively, formal upscaling approaches have been used to derive effective rate expressions where small-scale mass transfer effects are affecting the values of the degradation rate or transport parameters at larger scales (Wood et al. 2007; Golfier et al. 2009; Hesse et al. 2009).

Interestingly, we are not aware of applications of microbial explicit RTMs at scales larger than the plot or the sediment core scale. For instance, the coupled microbial-geochemical dynamics has not yet been addressed at the catchment scale, at least with regard to processes associated with subsurface flow. This contrasts with catchment-based RTM applications simulating surface flow hydrodynamics in, e.g., rivers, lakes or estuaries where heterotrophic and autotrophic microbial community dynamics has already been explicitly simulated (e.g., Vanderborght et al. (2002), Billen et al. (2013)). In these applications, the interaction between substrates and microbial biomasses ignore the potential effects of spatial heterogeneities at the sub-grid scale and the concentrations are thus assumed homogeneous at scales smaller than the model resolution (typically hundreds to thousands of meters). In stark contrast to subsurface catchment scale models, several biomass-explicit models have been developed for global-scale applications over the last 5–10 years, mostly in relation to soil C dynamics. In terms of biological process representation, these large-scale microbial models are not fundamentally different from those applied at (much) smaller spatial scales. A major alteration, however, lies on the assumption that soil physics and microbiology are homogenous within the soil grid to which it is applied (e.g., Huang et al. (2018)). Taken that the typical resolution of land surface schemes of Earth System Model is on the order of 0.5–2° (i.e. from tens to hundreds of km), numerous physical–biogeochemical–microbial interactions occurring at smaller spatial scales are unresolved in these models and their parametrization remains a grand challenge for the future.


The advancement in experimental methods has provided an enormous amount of information on microbial community dynamics in the subsurface and the ongoing development of new methods for the analysis of microbial processes will further increase our knowledge in this field. As outlined in this chapter, this has triggered the development of a large number of modeling concepts that support the exploitation of this new knowledge in the broader context of reactive-transport modeling. Not all of these concepts have yet been embedded in reactive transport modeling approaches, but the existing combinations indicate the strong potential of improving the representation of microbial processes in such models. Next to technical developments there is however the need to establish effective modeling concepts that can be applied for the simulation of microbial community dynamics and their biogeochemical function in natural porous media. The complex behavior of microbial communities still imposes severe challenges in their experimental assessment as observations obtained from any location in the subsurface may integrate the effects of several micro-environments each with distinct communities. Furthermore, many microbial species in a natural community are still not known or characterized. The required effective model approaches thus imperatively need to reduce the real world complexity of a natural community into a less complex community represented in the model framework. However, this reduced community should still reproduce the dynamics of the real community adequately. How to meet this challenge in the best possible way depends not only on the available data but also on the microbial functions of interest and thus on the purpose of the specific reactive transport modeling effort. For many subsurface compartments there is a growing need for reactive transport approaches to consider the functional dynamics of the microbial community more adequately, from the pore scale to the global scale. This is in particular the case for the Critical Zone and its compartments (Griebler et al. 2014; Vereecken et al. 2016; Li et al. 2017; Baveye et al. 2018; Vogel et al. 2018) partly because surface processes (e.g., weather, climate, anthropogenic activities) lead to a dynamic variation in the environmental conditions. However, this can also be the case for marine sediments (van de Velde et al. 2018). The increasing global population will lead to further anthropogenic disturbances of the subsurface and global change is expected to lead to more dynamic weather conditions. This further increases the need for model approaches that can predict the response of the microbial functions to such perturbations together with the adaptation to a changing climate.


MT was supported by funding from Helmholtz Centre for Environmental Research—UFZ via the integrated project Controlling Chemicals Fate (CCF) of the research topic Chemicals in the Environment (CITE) and by the German Research Foundation (DFG) CRC 1076 “AquaDiva”. PR received funding from the VERIFY project from the European Union's Horizon 2020 research and innovation programme under grant agreement No. 776810 and from the F.R.S.-FNRS (sabbatical leave). The authors thank D. LaRowe for commenting an earlier version of this chapter and N. Roevros for her technical support.

Open access: Article available to all readers online.