Beneath our feet is a fascinating world of flowing water, cosmopolitan microbes, and complex mineral assemblages. Yet we see none of it from above. Our quest to investigate these complex subsurface interactions has led to the development of reactive transport models. These are computer algorithms that allow us to explore, in a virtual way, the natural dynamics of Earth's systems and our anthropogenic impact on those systems. Here, we explain the concepts behind reactive transport models—which include the transport of aqueous species and the descriptions of biogeochemical reactions involving solutes, surfaces and microorganisms—and introduce to reactive transport applications in terrestrial and marine environments.

INTRODUCTION

Despite the vast technological and manufacturing innovations of the Second Industrial Revolution (~ 1870 to 1914), the chemical engineer of the 1930s still relied heavily on trial and error to optimize the design of industrial processes, from chemical separations and petroleum refining to combustion engines (Inger 2001). The question these engineers sought to answer was seemingly simple: how long did a fluid or gas need to stay in the reactor to achieve maximum conversion? The answer was also simple, in principle: the time the fluid or gas spends in the reactor should be longer than the timescale of the chemical reaction(s). Quantitative understanding, however, required the development of mathematical constructs. To address this problem, a young assistant professor at Göttingen University (Germany), Gerhard Damköhler, explored an extension of similarity theory to “nondimensionalize” the governing equations of reactive flows. What emerged was a set of dimensionless numbers that bears his name: the Damköhler (Da) numbers. The Da compares the timescale for solute or gas transport (i.e., how long it takes for a hypothetical particle of fluid or gas to traverse the reactor) to the timescales for the reactions to go to completion, i.e., how long it takes to complete a set of reactions as governed by the kinetics and thermodynamics that set the rates of reactions and final concentrations (Damköhler 1936). In this sense, the ratio reflected in the Da number is useful in that it defines the competition between transport and reaction (Fig. 1). Although developed in the context of chemical engineering, this concept is equally valid for the interactions between the reactions and transport processes below our feet, i.e., in the Earth's upper crust.

Today, the Da number has come to define the overarching principle of reactive transport. All open systems involve a competition between transport and reaction. Thus, the advantage of the Da number is that, although it is based on a mathematical construct, it does not require a computer model to evaluate. It is, in effect, a very practical rule-of-thumb. However, in practice, most open systems are much more complex than a singular Da number can describe, particularly if your reactor is the Earth. Increasingly, reactive transport problems are tackled numerically, allowing a continual proliferation of application areas and scientific questions. As the reactive transport concepts inspired by early chemical engineers infused into the physical sciences, including the geosciences, they became the basis for the modern field of “reactive transport modeling”. Today, advanced algorithms rapidly translate our conceptual models into numerical representations, pulling in carefully compiled data on the thermodynamics and kinetics of relevant species, as well as information about the porosity, permeability, and chemical composition of the porous media. As such, these algorithms provide a multidimensional tracking system for mineral and microbial reactions, as well as migrating chemicals or contaminants.

Given the power to peer virtually into the subsurface, almost all fields of the geosciences have been impacted by reactive transport approaches (Steefel et al. 2005). This issue of Elements focuses on the most common applications of reactive transport models (RTMs): those that involve low-temperature fluids and gases (~ <150 °C) in the subsurface and the minerals and microorganisms that these fluids interact with. In particular, the articles in this issue reflect on how RTMs have uniquely enabled us to obtain new insight into complex natural and engineered processes. In addition to the topics contained in this issue of Elements (Fig. 2), geoscientists over the decades have also used reactive transport principles to understand processes occurring in deeper environments and at higher temperatures, including magma chambers (e.g., Richter 1986) and metamorphic reaction fronts (e.g., Baumgartner and Ferry 1991).

FUNDAMENTAL PRINCIPLES OF REACTIVE TRANSPORT

In the subsurface, we encounter a variety of minerals of wildly different kinetic and thermodynamic properties, coated by surfaces decorated with a complex array of ions and molecules. Interacting with these “messy” surfaces are fluids, gases, and hungry microorganisms. To track all of these phases and components over space and time requires a robust approach. Hence, we rely on a menu of mass balance equations to provide a palette for not only visualizing the subsurface but for quantifying the emergent dynamics of systems. The mathematical and numerical approaches themselves can be gleaned from Steefel et al. (2015); here, we focus here on the conceptual and practical considerations of reactive transport formulations by using only simple one-dimensional representations.

Conceptualizing Fluid Flow and Chemical Reactions at the Earth's Surface

Systems with no physical transfer of reactants and products, or perfectly closed systems, can be conceptualized as “batch” systems where the chemical transformations are purely a function of time (Fig. 3, top panel and column A). However, most systems on Earth are characterized by some exchange of material, depending on the temporal and spatial scales over which we can observe them. The simplest mass balance description for an open system is usually the “well-mixed flow reactor” (WMFR). This is our early chemical engineer's well-mixed tank, which has fluid (in this case, water) flowing into and out of the reactor and facilitating the reactions occurring within (Fig. 3B). Oceans and lakes are often conceptualized this way. However, other systems, such as rivers and aquifers, are not well mixed (i.e., they display chemical gradients in space and time). A river or aquifer can instead be thought of as a series, or chain, of WMFRs linked together along a river's reach or aquifer flow path, with each reactor accomplishing a fraction of the overall reaction progress (Fig. 3C). Within an RTM, WMFRs make up each of the “grid cells” through which fluids, reactants, and products traverse, encountering a variety of stationary solid phases along the way, including surfaces covered in sorbed species and microbes. Hence, if a system can be conceptualized as a WMFR, you may not need an RTM: a multicomponent geochemical model may be sufficient. However, if the system of interest involves the evolution of chemical constituents not only in time but also in space, then an RTM approach is required to fully address the coupling between reaction and transport.

How do we set the fluxes of fluids and elements into and out of the grid cells? Conceptually, these require a mass balance description of flow and transport. In the case of aquifers, water flow is usually described by Darcy's law, where the water flux is proportional to the product of the hydraulic conductivity (i.e., the ease with which a fluid can move through pore spaces or fractures) and the pressure gradient. In addition, the transport of dissolved species requires a description of dispersivity and/or diffusion to fully describe the mass fluxes. In these formulations, hydraulic conductivity and transport parameters constitute upscaled and averaged conditions that are valid on the centimeter to kilometer scale. Models numerically capture fluid flow and transport processes in a variety of ways, and some include a capability for multiphase flow and sophisticated treatments of diffusive and dispersive processes (Steefel et al. 2015). Although powerful in many ways, these continuum modeling approaches are not well suited to deal with problems at the grain and/or micro scale. Recently, a new generation of pore-scale models have been developed to resolve fluid velocity gradients at the scale of individual pores (Molins 2015).

Modeling Systems of Reactants and Products

Once the flow and transport parameters have been determined, then the reactants and products defined by the stoichiometry of the relevant chemical reactions, including microbial reactions, form the reaction network. These reactions may include the formation of aqueous complexes that affect the sorption and redox potential of ions, surface complexation and the exchange of cations, the oxidation–reduction reactions, and the mineral dissolution and precipitation reactions (Bethke 2010). Importantly, the reactive transport approach honors the stoichiometric relationships between reactants and products, subject to thermodynamic constraints. Several different databases containing the reactions and the best available equilibrium constants have been developed to provide internally consistent data (e.g., Wolery and Colón 2017). Using such databases, geochemical speciation calculations, constrained by thermodynamic principles, can be used to determine the distribution of species in aqueous solutions, on mineral surfaces, and in solid phases. However, many reactions of interest are sufficiently slow to also be governed by the reaction rate, or kinetics, in addition to thermodynamics. Thus, to complete the reaction network, reactive transport models must be able to accommodate varied rate expressions and associated rate coefficients to adequately describe reaction progress and geochemical evolution. Most kinetic descriptions are derived from experimental studies. Thus, the integration between models and controlled experiments tends to strengthen the utility of the models. For a full description of common geochemical capabilities see Table 2 in Steefel et al. (2015), and for an overview of mineral–dissolution/precipitation rate laws see Palandri and Kharaka (2004).

Capturing Connections Between Chemical Reactions and Transport Using Constitutive Equations

Solving the equations for flow, transport, and a host of bio- and geochemical reactions is useful in many contexts. However, other “constitutive” equations are often needed to link these various equations together. For example, over time, new minerals may precipitate, or microbial biomass may grow, which tends to cause a reduction in porosity and which can affect both the permeability and the diffusivity of ions and molecules. Thus, an additional set of relationships—such as the Kozeny–Carman equation for porosity–permeability feedbacks, or Archie's law for porosity–tortuosity interactions—relate the volume changes associated with geochemical reactions to the transport processes. This results in either negative or positive feedbacks between reactions and fluid flow. Positive feedbacks between chemical reactions and porosity can be particularly relevant for highly reactive systems (e.g., geologic CO2 storage), where increased porosity accelerates flow and, thus, enhances dissolution, leading to the phenomenon of wormholing. Negative feedbacks, or clogging, occur when a decrease in porosity limits fluid transport. Such negative feedbacks can be driven by biological processes (e.g., microbial biofilms) or net reactions that result in a positive volume change. Recent progress has been made using pore-scale models to understand these constitutive equations at a phenomenological level (Molins 2015).

EARTH SCIENCE APPLICATIONS OF REACTIVE TRANSPORT MODELS

We next provide a general overview of the myriad ways in which RTMs are used in the geosciences. In order to show the opportunities for integration across time and disciplines that RTMs present, we will take a somewhat coarse walk through geologic time and environments.

Decoding Earth's History

Most of Earth's history is recorded in chemical discontinuities in rocks and sediments. Long-term trends in elements or isotopes that are now considered as proxies for environmental conditions in the ocean/atmosphere, or on land, have led to pronounced signatures in the sedimentary rock record, often perceived as “jumps” or “wiggles”. To decipher the environmental and diagenetic processes that imparted these chemical signals, RTMs can be used to simulate conditions from the sediment–water interface to, potentially, hundreds of meters below it (Wang and Van Cappellen 1996). For example, in order to understand how the burial and diagenesis of phosphorous is altered by a shift from anoxic to oxic bottom waters, Dale et al. (2016) used an RTM to show how even shallow bioturbation, in a period such as the early Paleozoic, may have strongly enhanced phosphorous burial and helped to stabilize atmospheric O2. Although diagenesis—which alters the porosity, permeability, and compositions of sedimentary rocks—can complicate paleoenvironmental records, it can also be critical in understanding the evolution of porosity in petroleum reservoirs. For example, dolomitization affects, globally, nearly 50% of carbonate rocks, and most massive dolomites appear to form during early burial of platform carbonates where circulating seawater provided a source of magnesium through shallow circulation zones, reflux of evapoconcentrated seawater, or geothermal convection. Modeling studies have been employed to examine these mechanisms by accounting for nonisothermal fluid and heat transport [e.g., patterns of dolomite formation due to geothermal convection of cold seawater through a carbonate platform (Whitaker and Xiao 2010)].

In addition to directly simulating events in Earth history, RTMs elucidate key connections between the various parameters, including the removal of atmospheric CO2 via chemical weathering due to hydrologic and erosional patterns (Goddéris et al. 2019 this issue) or the relationships between organic matter decomposition in marine sediments and carbon source and age (Arndt et al. 2013). Using RTMs to understand the processes that sculpt our planet can integrate the disciplines of paleoceanography, sedimentology, geochemistry, hydrology, geophysics, and geothermal research. And this itself could bring new insights into Earth's wiggles, jumps, and trends.

What Shaped the Earth's Surface?

The surface of the Earth today carries the legacy of past events, from the geomorphic and climatic perturbations of the ice ages to the effects of human-induced land-cover changes. With the ability to simulate from a point in the past forward to today—a process called history matching—RTMs uniquely allow us to study the effect of changing boundary conditions. Examples include soil formation and weathering of the regolith (Li et al. 2017), the formation of ore deposits in sedimentary systems (e.g., Garven et al. 1993) and in hydrothermal systems (Raffensperger and Garven 1995), and the serpentinization of the upper mantle (Tutolo et al. 2018). Using models to understand which factors drive the distribution of elements we see today, in turn defines key environmental gradients, such as the lithological, climatic, and ecological variability that dictate the thickness of a weathered regolith and the importance of temperature as a control on the kinetics of serpentinization.

The Influence of Humans

Another application of RTMs is to quantify, and potentially ameliorate, the influence of humans—industrial activity not only provided us with the underlying mathematical construct for reactive transport, it has also left us with a legacy of soil and groundwater contamination. Reactive transport models have been useful as a critical tool for assessing soil and groundwater contamination. They have also been useful for developing and implementing methods for groundwater remediation, ranging from extraction of contaminated groundwater to in situ chemical treatment that changes the contaminant to a less mobile and, ideally, less toxic chemical form [e.g., uranium reduction by stimulated microorganisms (Li et al. (2009)]. Chemical spills have long presented a target for reactive transport modeling (Mayer et al. 2002), as has the need to safely dispose of energy by-products, namely, spent nuclear materials and carbon dioxide, which have been explored extensively using RTMs. As developed in other articles in this issue, models have allowed us to ask about the future of the different strategies for contaminant remediation or waste storage over timescales ranging from months to millions of years.

FUTURE DIRECTIONS

Existing modeling approaches are already advancing our understanding of the reactive transport processes that shape the Earth's surface and that control the fate of man-made contaminants. However, the exponential increase in Earth data resulting from airborne, satellite, and ground-based sensor technology is stimulating new approaches for data archiving, visualization, and assimilation, and, by extension, data-model integration and uncertainty quantification. Thus, new modes of data collection must be combined with robust frameworks for understanding uncertainty in the data and in the models themselves (Scheidt et al. 2018). The modelers toolkit (see Maher and Mayer 2019 this issue, page 117) will include a bevy of spatial and spatiotemporal statistical techniques, or collaborators with such knowledge.

Just as sensor technology is expanding the universe of Earth data, new three-dimensional and time-resolved nanoand microscale characterization approaches, including microcomputed X-ray tomography (~ 5–30 micrometer resolution) and full-field transmission X-ray spectroscopy (~ 2–40 nanometer resolution) are revealing the assembly of Earth materials at the pore scale. New modeling approaches enable us to leverage these techniques to evaluate pore-scale transport and reaction at the surfaces of individual grains (Molins 2015), as opposed to the continuum approach shown in Figure 3C. These new approaches currently hold promise for upscaling interfacial processes and may, ultimately, allow us to achieve multiscale simulation of environmental systems. Independent of scale, availability of high-resolution and high-quality data will allow us to better constrain RTMs, which, like all complex models, can be substantially affected by uncertainties and model nonuniqueness.

Although microorganisms have been included in models for decades, models are now starting to incorporate plant uptake and return of water and nutrients. For example, implementation of dynamic root-growth models into models to investigate feedbacks between growth and water or nutrient availability now extend from the single-plant (Gérard et al. 2017) to the plot-scale rhizosphere (Gérard et al. 2008), thereby expanding the scope of reactive transport modeling and further increasing interdisciplinary integration and reach.

Finally, research teams are becoming increasingly larger and inherently more interdisciplinary, a trend that also applies to the data streams. For example, field-based geoscience projects will often combine hydrologic and biogeochemical measurements with geophysical characterization, while also reaching across scales to conduct laboratory experiments and/or molecular-level characterization of field samples. The geosciences are already inherently interdisciplinary with respect to time—studies of modern systems may be aimed at reconstructing past intervals of Earth history or predicting the rock record of the future through natural analogue studies. Reactive transport models are bridging across disciplines, providing not only critical scientific tools for interrogating the Earth but also platforms for galvanizing interdisciplinary collaboration in the future. As we look forward, we see the potential for scientifically diverse teams to organize around RTMs not just as tools for understanding, but as tools for collaboration.

IN THIS ISSUE

Reactive transport concepts have become fundamental to a vast array of fields within the Earth sciences (Steefel et al. 2005). In this issue of Elements, we focus on topics extending throughout Earth's shallow crustal environments. Within this broad spatial domain, Christophe Meile and Timothy Scheibe (2019 this issue), in describing the assimilation of microbial processes into RTMs, outline our understanding of how the interactions between life and the environment unfold. Building on this contribution, Jennifer Druhan and Matthew Winnick (2019 this issue) reveal the importance of including isotopes into models of complex systems in order to constrain elemental transfers, ranging from those at the surfaces of minerals, to water vapor within the atmosphere. Yves Goddéris and coauthors (2019 this issue) focus on the evolution of the Critical Zone, past and present, using models to understand how water, erosion, and organisms shape the regolith. Because human activities generate distinct perturbations to natural systems, Henning Prommer and coauthors (2019 this issue) discuss how RTMs provide insight into water quality, including the impacts of both natural and human activity. Laurent De Windt and Nicolas Spycher (2019 this issue) detail the role that reactive transport modeling has played in understanding and ameliorating the risks associated with nuclear waste storage in deep geological repositories. Finally, humanity must come to terms with our dependence on fossil fuels and, in the short term, this may require engineered containment of carbon dioxide. Anna Harrison and coauthors (2019 this issue) outline the role that RTMs have played in developing strategies for both containing and optimizing the conversion of CO2 to carbonate minerals, reminding us of the important advances originally made by Gerhard Damköhler and his pioneering work on reactive transport.

Acknowledgments

KM acknowledges research support from the National Science Foundation (EAR-1254156) and the U.S. Department of Energy, Office of Biological and Environmental Research, Subsurface Biogeochemical Activity Program (DOE-DESC0014556; DE-AC02-76SF00515). KUM acknowledges research support from the Natural Sciences and Engineering Research Council of Canada in form of a Discovery Grant. Carissa Carter is thanked for her advice on visualization and Carl Steefel for his review. We are grateful to Principal Editor Friedhelm von Blanckenburg for his thoughtful for reviews and his support for this issue of Elements.

This is an open-access article distributed under the terms of the Creative Commons Attribution CC-BY-NC-ND 3.0 License, which permits users to copy and redistribute the work, provided this is not done for commercial purposes and further does not permit distribution of the work if it is changed or edited in any way, and provided that the user gives appropriate credit, provides a link to the license, and that the licensor is not represented as endorsing use of the work.Open access: Article available to all readers online. This article is CC BY-NC-ND.

Figures & Tables

Two different flow paths for a reactive chemical species result in different times spent in contact with a porous media. (top) The chemical species following the red flow path through relatively large grains traverses the entire reactor quickly with minimal time to react and, therefore, has a low Damköhler number (Da). (middle) The chemical species following the tortuous blue flow path through grains both large and small spends longer in the reactor, thereby allowing more time to convert from a reactant into a product resulting in a high Da. In reality, there may be a distribution of travel times within a given porous media. (bottom) Graph of the reactant concentrations, C [mol m-3] for the two systems relative to the initial concentration and equilibrium concentration. If the flow rate is very high, then the travel time is very short and conversion of the reactant is incomplete. If the flow rate is very low, most of the reaction occurs near the beginning of the reactor.

Two different flow paths for a reactive chemical species result in different times spent in contact with a porous media. (top) The chemical species following the red flow path through relatively large grains traverses the entire reactor quickly with minimal time to react and, therefore, has a low Damköhler number (Da). (middle) The chemical species following the tortuous blue flow path through grains both large and small spends longer in the reactor, thereby allowing more time to convert from a reactant into a product resulting in a high Da. In reality, there may be a distribution of travel times within a given porous media. (bottom) Graph of the reactant concentrations, C [mol m-3] for the two systems relative to the initial concentration and equilibrium concentration. If the flow rate is very high, then the travel time is very short and conversion of the reactant is incomplete. If the flow rate is very low, most of the reaction occurs near the beginning of the reactor.

Summary of the breadth of reactive transport modeling applications across the geosciences. These include marine and terrestrial systems, across time from past to future environments, modern natural systems, and human-impacted systems. Geoscience topics not shown include ore deposits, geothermal systems, and petroleum applications. Reactive transport models, including those that incorporate isotopic fractionation and microbial processes, therefore, have great applicability across and between a multitude of geological systems.

Summary of the breadth of reactive transport modeling applications across the geosciences. These include marine and terrestrial systems, across time from past to future environments, modern natural systems, and human-impacted systems. Geoscience topics not shown include ore deposits, geothermal systems, and petroleum applications. Reactive transport models, including those that incorporate isotopic fractionation and microbial processes, therefore, have great applicability across and between a multitude of geological systems.

(top) A graph of reactant concentration, C [mol m−3], versus time, t [s], (i.e., time-dependent chemical reaction) for three different mathematical conceptualizations of reactive transport. The reaction rate is described by a first-order rate constant, k [s−1]. (A) A batch reactor at a given time with reactant (red spheres) and product (yellow prisms and layers). (B) A well-mixed flow reactor where the change in concentration with time (∆C/∆t) is set by the volumetric flow rate into and out of the reactor, Q [m3 s−1], concentration at the inlet, and the reactions that occur within the control volume, V [m3] for a given residence time (V/Q). (C) The concentration along the flow path in a flow-through reactor (e.g.,Fig. 1) can be conceptualized in terms of a series of infinitely small reactors by defining Q in terms of the velocity, v [m s−1], across the face of the reactor, A [m2], at position x, and out of the reactor at position x + ∆x. Dividing by A∆x and taking the limit as ∆x goes to zero results in the advection–dispersion/diffusion–reaction equation, which is the fundamental partial differential equation solved numerically in reactive transport models. The second-order differential term represents the diffusive and dispersive fluxes, where D* is the sum of the molecular diffusion and mechanical dispersion coefficients.

(top) A graph of reactant concentration, C [mol m−3], versus time, t [s], (i.e., time-dependent chemical reaction) for three different mathematical conceptualizations of reactive transport. The reaction rate is described by a first-order rate constant, k [s−1]. (A) A batch reactor at a given time with reactant (red spheres) and product (yellow prisms and layers). (B) A well-mixed flow reactor where the change in concentration with time (∆C/∆t) is set by the volumetric flow rate into and out of the reactor, Q [m3 s−1], concentration at the inlet, and the reactions that occur within the control volume, V [m3] for a given residence time (V/Q). (C) The concentration along the flow path in a flow-through reactor (e.g.,Fig. 1) can be conceptualized in terms of a series of infinitely small reactors by defining Q in terms of the velocity, v [m s−1], across the face of the reactor, A [m2], at position x, and out of the reactor at position x + ∆x. Dividing by A∆x and taking the limit as ∆x goes to zero results in the advection–dispersion/diffusion–reaction equation, which is the fundamental partial differential equation solved numerically in reactive transport models. The second-order differential term represents the diffusive and dispersive fluxes, where D* is the sum of the molecular diffusion and mechanical dispersion coefficients.