INTRODUCTION

Reactive transport in the Earth and Environmental Sciences is at a crossroads today. The discipline has reached a level of maturity well beyond what could be demonstrated even 15 years ago. This is shown now by the successes with which complex and in many cases coupled behavior have been described in a number of natural Earth environments, ranging from corroding storage tanks leaking radioactive Cs into the vadose zone (Zachara et al. 2002; Steefel et al. 2003; Lichtner et al. 2004), to field scale sorption behavior of uranium (Davis et al. 2004; Li et al. 2011; Yabusaki et al. 2017) to the successful prediction of mineral and pore solution profiles in a 226 ka chemical weathering profile (Maher et al. 2009), to the prediction of ion transport in compacted bentonite and clay rocks (Tournassat and Steefel 2015; Soler et al. 2019; Tournassat and Steefel 2019, this volume). Yet for those thinking deeply about Earth and Environmental Science problems impacted by reactive transport processes, it is clear that many challenges remain.

A common theme here is that improved scientific understanding of complex Earth and environmental processes become possible as new conceptual and numerical models for reactive transport are developed. In some examples, this success has hinged on demonstrating that reactive transport models were computationally feasible. In others, it was also critical to employ the flexible framework of these multi-component, multi-dimensional models to show that it was possible to parameterize them for application to complex field systems. Certainly, the computational burden increased with the use of these models, as well as the considerable effort and knowledge required of researchers to develop and apply them. The benefits that have resulted include a much improved and generalized set of model capabilities, as well as a deeper level of understanding of the underlying coupled physical and biogeochemical processes of Earth systems. The maturity of the discipline has also been demonstrated by recent benchmarking activities that have shown that as many as ten different software packages can simulate complex natural reactive transport problems and achieve essentially the same results (Steefel et al. 2015).

So, is the “reactive transport problem” solved? Where do we go from here? This is the crossroads we are at now as we decide what are the challenges that need to be faced so as to continue advancing the field. Arguably the most significant challenges we now face are associated with the huge range of length scales that need to be addressed, since these extend from the molecular to nanoscale to pore scale all the way up to the watershed and continental scale (Molins and Knabner 2019, this volume). Across this extreme range of scales, the constitutive equations and parameters that are used to describe reactive transport processes often change as well, thus requiring mathematical and numerical models to become “scale aware”. Charged porous media offers special challenges, since ion mobility can be strongly affected by electrostatic interactions, and this can lead to effects such as anion exclusion that are not captured by Fick's Law (Appelo and Wersin 2007; Appelo et al. 2010; Tournassat and Steefel 2015). Where the charged porous media involves nanoscale porosity, off-diagonal coupling effects on transport between such master variables as fluid pressure, electrical current and chemical composition may become important (Tournassat and Steefel 2019, this volume). At the watershed and continental scales, reactive transport is further complicated by the coupling with diverse Earth surface processes, including subsurface and surface water, vegetation, and the atmosphere, all played out typically in highly heterogeneous and transient settings (Li 2019, this volume). This can lead to “hot spots” and “hot moments” that may have an outsized effect on system function even though they represent a limited percentage of the total land surface area (Dwivedi et al. 2018a).

HISTORICAL DEVELOPMENT

Geochemical modeling of Earth and environmental systems began with equilibrium descriptions of the thermodynamic state of a particular multi-component solution, and such approaches are still in wide use today for the purposes of interpreting the chemistry of natural waters. An important next step was the development of reaction path models that captured the sequence of chemical/mineralogical states resulting from such processes as chemical weathering and hydrothermal alteration. It could be argued, in fact, that these reaction path models represent the first “reactive transport” models insofar as they address irreversible geochemical/mineralogical processes for the first time. Early versions of the reaction path models, pioneered by Helgeson and co-workers, did not include an explicit treatment of real-time kinetics, but rather quantified the system evolution instead as a function of reaction progress (Helgeson 1968; Helgeson et al. 1969). These models provided a way of interpreting quantitatively the sequence of minerals observed in nature as the natural consequence of the dissolution of some primary phase (e.g., feldspar), which is itself out of equilibrium due to either the initial state of the system, or more commonly, due to the flux of reactive constituents.

An explicit treatment of transport processes is not factored into these reaction path approaches. The method can be used to describe chemical processes in a batch or closed system (e.g., a laboratory beaker), or for exceedingly simplified transport. However, such conditions are of limited interest in the geosciences where the driving force for most reactions is often transport. Lichtner (1988) clarified the application of the reaction path models to water–rock interaction by demonstrating that they could be used to describe pure advective transport through porous media. By adopting a reference frame which followed the fluid packet as it moved through the medium, the reaction progress variable could be thought of as travel time instead.

Multi-component reactive transport models that could treat any combination of transport and biogeochemical processes date back to the mid-1980s. Lichtner (1985) outlined much of the basic theory of a continuum model for reactive transport. Yeh and Tripathi (1989) also presented the theoretical and numerical basis for the treatment of reactive contaminant transport, demonstrating for the first time the inapplicability of classical linear distribution, or Kd, approaches in predicting contaminant dynamics. Steefel and Lasaga (1994) presented a reactive flow and transport model for non-isothermal, kinetically controlled water–rock interaction and fracture sealing in hydrothermal systems based on simultaneous numerical solution of both reaction and transport. This study was the first to consider multicomponent reactive transport in the context of non-isothermal flow fields, an important subject for geothermal and hydrothermal systems. It was apparently also the first to consider how reaction-induced permeability change (clogging) could alter the behavior of a hydrothermal or other flow system. Along with earlier studies of mineral dissolution related feedbacks to the flow system through permeability (Ortoleva et al. 1987; Steefel and Lasaga 1990), this work introduced the important theme of coupled processes in reactive transport analysis (Seigneur et al. 2019, this volume).

In what follows, we briefly review the state of the science in terms of conceptual and mathematical models for reactive transport, considering how the equations change as we proceed down scale from continuum to pore to molecular models. This is followed with an overview of the numerical methods that are used to solve these potentially complex problems. We then proceed to discussions of individual topics where reactive transport has had a significant impact since the publication of Reviews in Mineralogy volume 34 in 1996 (Lichtner et al. 1996) including an update to the persistent issue of contaminants.

MATHEMATICAL FORMULATION

Presentations of the basic reactive transport equations can take different directions depending on the interests of the author, but for this Rev Mineral Geochem volume, it seems appropriate to emphasize the geochemical and mineralogical aspects. The primary variables in the reactive transport problem are the concentrations of component i (Ca, Na, …) in phase γ (aqueous, gas, and mineral). Typically, reactive transport equations are developed for each of the phases, or they are combined with an overall mass balance that covers all of the phases present, with differing transport properties for each of the phases. The concentrations, along with temperature and pressure, are the primary variables. Various derived or secondary variables follow from the primary variables, including the reaction rates (for minerals, dependent on the aqueous and/or surface concentrations and the mineral surface area), the mineral volume fractions (essentially linear combinations of the component concentrations in the mineral phase), and the mineral surface areas (related to the volume fractions through the specific surface area). Whether the mineral volume fractions are formally included within a single nonlinear solve depends on the software considered. Often the time scale separation resulting from slow kinetics of common mineral reactions in many cases justifies the update of these at the end of the time step. In the cases where all of the components are in equilibrium in the multiple phases present, we regain the Gibbs Phase Rule, with degrees of freedom equal to the number of components plus 2 (for T and P).

The mineral volumes are further related to the porosity of the medium through the relation:

 
ϕ=1Σm=1Nmϕm
(1)

The permeability of the medium, which determines the flow rate for a given pressure or hydraulic head gradient (thus completing the feedback circle), is determined by a porosity–permeability relationship applicable to the medium in question. Typically, the Kozeny–Carman equation is used for porous media (Bear 1972):

 
k=ϕ3(1ϕ)215M2
(2)

where M is the fluid–solid interfacial area. The applicability of this relation has been questioned in many cases, particularly for reaction-induced porosity–permeability change. For fractured media, the cubic law is used:

 
k=b312d
(3)

where b is the fracture aperture and d is the separation of the fractures perpendicular to the fracture plane (Witherspoon et al. 1980; Steefel and Lasaga 1994; Deng et al. 2016). The effective diffusivity in porous media, which includes the effect of tortuosity, also depends on the porosity, often through Archie's Law (Seigneur et al. 2019, this volume).

The geochemical relations in Figure 1 form the basis for the reactive transport equations. For a fully kinetic formulation in which each aqueous species is treated separately (i.e., they are not assumed to be in equilibrium), the equation is given by (Steefel et al. 2015):

 
(ϕSLCi)tAccumulation Term=(ϕSLDi*Ci)Dispersion(qCi)Advectionr=1NrvirRrAqueous Reactionsm=1NmvimRmMineral Reactionsg=1NgvigRgGas Reactions
(4)

In Equation (4), SL is the liquid saturation in the porous media, D* is the dispersion tensor, q is the Darcy velocity, and the νir, νim and νig refer to the stoichiometric coefficients for component i in aqueous reaction r, mineral reaction m, or gas reaction g respectively (Steefel et al. 2015).

For the partial equilibrium case, where some number of aqueous and/or surface complexes are assumed to be in equilibrium and thus described by algebraic mass action equations, we have:

 
(ϕSLΨi)t=(ϕSLDi*Ψi)(ΦSLνϕi)Σr=1NrνirRrΣm=1NmνimRmΣg=1NsνigRg
(5)

where the total concentration, Ψi is defined as a linear combination of the concentrations of the primary and secondary species:

 
ΨiTotalConcentration=CiPrimarySpecies+l=1Nsvi,lClSecondarySpecies=Ci+l=1Nsvi,l[γpCpvp,lKl1]
(6)

where γp is the activity coefficient for the primary species, and Kl is the equilibrium constant for the secondary species reaction. Surface complexes can be added similarly. The right most term in Equation (6) is the total concentration rewritten entirely in terms of primary species—the so-called “Direct Substitution Approach”. The Direct Substitution Approach (of DSA) involves making use of the laws of mass action for secondary species assumed to be in equilibrium with the primary species.

Most of the software in use today for continuum reactive transport makes use of some form of Equation (6), with or without the Direct Substitution Approach, although alternative formulations have been proposed (Yeh et al. 2014). The total concentration approach (or “canonical formulation”, Lichtner 1985) lends itself to treatment of reactions with a mixed equilibrium and kinetic approach (see Steefel and MacQuarrie 1996 for a more detailed discussion). Other approaches, particular involving the use of Gibbs free energy minimimization (Steefel and MacQuarrie 1996; Leal et al. 2014) are now in common use, although they require the use of separate kinetic routines if such chemical processes are to be included.

It is also possible to treat fully equilibrium reaction networks with the canonical formulation, as described in Steefel et al. (2015). In this case, all of the reaction rates are eliminated, replaced typically with the mass action equations in logarithmic form:

 
1Δt[Ci+Σl=1NsνilCl+Σm=1NmνimCm]n+1[Ci+Σl=1NsνilCl+Σm=1NmνimCm]n=0i=1,Nc
(7)
 
logCl=Σl=1Nsνillog(γiCin+1)logγln+1logKl1=1,NslogCmn+1=Σm=1Nmνimlog(γin+1Cin+1)logKmm=1,Nm
(8)

where n and n + 1 refer to the present and future time level respectively, and the γi's and γj's are the activity coefficients for the secondary and primary species respectively, and Kl and Km are the equilibrium constants for the secondary species and minerals, respectively. A similar approach can be taken to include equilibrium surface complexes and gases.

Details of the formulations for the transport terms (Darcy's Law, Fick's Law) are provided in other chapters in this volume, as well as Steefel et al. (2015).

NUMERICAL FORMULATION

The numerical treatment of reactive transport was described in some detail in Steefel and MacQuarrie (1996) and in Steefel et al. (2015). In what follows, we first consider the treatment of the typically nonlinear reaction term first without transport (i.e., a set of nonlinear ordinary differential equations). Then we proceed to a brief discussion of how reactions and transport can be coupled over the spatial domain (a set of nonlinear partial differential equations).

Reaction terms

Steefel and MacQuarrie (1996) reviewed the range of fully kinetic formulations and equilibrium or mixed equilibrium–kinetic formulations for the reaction terms. Steefel et al. (1996) presented the mixed equilibrium–kinetic approach based on the canonical formulation (Lichtner 1985), essentially the numerical approach commonly employed by many modern reactive transport software packages such as CrunchFlow, MIN3P, and PFLOTRAN. Other mathematical and numerical treatments, however, have been considered. For example, Yeh and co-workers (Yeh and Tsai 2014; Yeh et al. 2014) have proposed and made use of an alternative to the canonical formulation.

The accumulation and reaction terms are discretized at the present and future time step, giving rise to a set of nonlinear ordinary differential equations that must be solved numerically in most cases. Various methods can be used to treat the time derivatives, such as backwards Euler (the simplest) and Runge–Kutta (the most common). The nonlinear ordinary differential equations can be solved with Newton's method given by:

 
Σk=1NcfiCkCk=fi
(9)

where fi are the function residuals (mass balance equations written typically in terms of the total concentrations in which the sum of the terms equals zero), and fiCk are elements of the Jacobian matrix (the derivatives of the function residuals, or mass balance equations, with respect to the primary k unknowns). For the case where only the accumulation and reaction terms are considered for the sake of simplicity (no transport terms), we have the ordinary differential equations given by:

 
ϕSL[ψi,jxn+1ψi,jxn]Δt+Σm=1Nm[Ri,mn+1]=0
(10)

Using the Direct Substitution Approach (rewriting the secondary equilibrium species appearing in the total concentration in terms of primary species), Equation (10) can be written for a single Newton iteration as:

 
ϕΔt[δi,jCi+l=1NsνklCkνil[p=1NcCpνplKeq1]+AmkmνkmCkp=1NcaiνpmKm1]n+1Jacobian Matrix=ϕΔt[(Ci+l=1Nsνil[p=1NcCpνplKeq1]in+1ψin)Amkm([1p=1NcaiνpmKm1])]Function Residuals
(11)

Note that even software employing an operator splitting approach (discussed below) to reaction and transport typically carries out a nonlinear solve of the accumulation and reactions terms. As discussed, the Gibbs free energy routines are based on minimization of free energy rather than root finding and thus do not follow this method.

Coupling of reaction and transport terms

To handle the coupling of reaction and transport processes over the spatial domain (now partial differential equations), the normal choices are between operator splitting approaches (either with or without iteration) and the global implicit or one-step approach.

The operator splitting approach, whether in continuum or pore scale models, is the most widely used and in non-iterative form begins with a solve of the conservative transport equations followed by the reaction terms using the transported concentrations in the accumulation term. To simplify the presentation, we introduce the differential operator

 
L(Ψi)=[(qD*)]Ψi,
(12)

and write the operator splitting approach as

 
ϕSL(ψitranspψin)Δt=L(ψi)n+1,(i=1,,Nc),
(13)

followed by a reaction step using the transported total concentrations:

 
ϕSL(ψin+1ψitransp)Δt=Σi=1NmRin+1(i=1,,Nc)
(14)

where ψin and ψin+1 are the concentrations at the current (n) and future (n + 1) time levels, respectively, and ψitransp is the transported total concentration calculated with Equation (13). The right hand side of Equation (13) could be written in terms of concentrations at the current time step (n), in which case the transport is said to be treated with an explicit in time approach. Alternatively, one can solve Equation (13) implicitly in time, in which case a speciation of the total concentration in the reaction step in Equation (14) is required, either using the Direct Substitution Approach as in Equation (6), or by including the mass action equations for the secondary species in logarithmic form, as in Equation (8).

In the global implicit or one–step approach, the transport and reaction terms are solved simultaneously. In the Direct Substitution Approach, we solve simultaneously for the complexation (which are assumed to be at equilibrium), the heterogeneous reactions, and transport terms. This means that the primary species (Cj's) rather than the total concentrations are the unknowns in the transport equations:

 
ϕSLΔt[(Ci+Σl=1Nsνil[Πp=1NcCpνplKeq1]in+1ψin)]=L[Cj+Σi=1Nrνijγi1Ki1Πj=1Nc(γjCj)νij]+Amkm([1Πp=1NcaiνpmKm1])
(15)

In the case of a global implicit treatment, the size of the Jacobian matrix which must be constructed and solved becomes larger, since each function will include contributions from the concentrations both in the grid cell itself and from neighboring grid cells that are used in the discretization. For example, in the case of one-dimensional transport and Nc unknown concentrations at each nodal point, the form of the Newton equations to be solved is:

 
Σk=1Ncfi,jxlnCk,jxδlnCk,j+Σk=1Ncfi,jxlnCk,jx+1δlnCk,j+1+Σk=1Ncfi,jxlnCk,jx1δlnCk,jx1=fi,jx
(16)

where i refers to the component number, k is the unknown component species number, and jx, jx + 1, and jx − 1 are the nodal points, and the logarithms of the concentrations are solved for because of the improved numerical stability this provides. The Jacobian matrix in the case of one-dimensional transport takes a block tridiagonal form

 
[A1,1A1,20A2,1A2,2A2,3AN1,N2AN1,N1AN1,N0AN,N1AN,N][δlnC1δlnC2δlnCN1δlnCN]=[f1f2fN1fN]
(17)

where the entries in the Jacobian matrix (the A) are submatrices of dimension Nc by Nc in the case where there are Nc unknowns per nodal point. The entries d ln Ci refer here to the entire vector of unknown concentration corrections in logarithmic form at any particular nodal point, while the functions fi include the entire vector of equations for the unknown concentrations at each nodal point.

The chief advantage of the operator splitting approach is the ability to use modular reaction routines without implementing a potentially complicated global implicit solve. There are no Jacobian entries for the transport terms in the operator splitting treatment of reactive transport, so less computation and less code development is involved. The approach works reasonably well where there is scale separation of processes, thus less coupling between reaction and transport. The disadvantage of the approach is that operator splitting error occurs when the time step is greater than the Courant number (Valocchi and Malmstead 1992). Applying the Courant condition would require that mass not be transported more than a single grid cell in any single time step—with a Courant number greater than 1, the transport skips over entire grid cells without reaction. This can be an issue in heterogeneous domains where locally very fast transport rates occur, driving the time step based on the local Courant condition to a very small value.

The chief advantage of the global implicit approach is that, unlike the operator splitting method for the reactive transport problem, one is not restricted to time steps less than the Courant condition. The other advantage of the global implicit is in systems where reaction and transport are strongly coupled. In such systems, an operator splitting approach may fail to converge to a fully coupled solution quickly. The quadratic convergence of the Newton method thus may provide substantially improved numerical stability and even faster execution time. The diffusion-dominant electrostatic calculations discussed later in the chapter and in Tournassat and Steefel (2019, this volume) are prime examples where the global implicit method dramatically outperforms the operator splitting approaches.

REACTIVE TRANSPORT ADVANCES

In the following sections, the author's view as to the advances in the development and use of reactive transport models since the publication of Reviews in Mineralogy Volume 34 in 1996 (Lichtner et al. 1996) is given. The list is far from exhaustive, however, and a more comprehensive discussion of advances in the application of reactive transport modeling are given in the other chapters in this volume.

Pore scale and continuum models

Historically, most of the developments and applications in reactive transport have centered on continuum scale models. Continuum models are those that treat the medium as having continuous properties, such as reactive surface area, mineral volume fractions, porosity permeability, and diffusivity (Fig. 2). This approach forms the basis of most of the reactive transport codes in use today (Steefel et al. 2015). However, an important recent development has been that of pore scale reactive transport models (Molins et al. 2012, 2014; Emmanuel et al. 2015). The pore scale models are distinguished from continuum models by the fact that interfaces between fluid, minerals, and gases are explicitly resolved.

The constitutive equations for continuum and pore scale reactive transport are different. For pore scale, flow is described with the Navier–Stokes (or Stokes) equation, while Darcy's Law is used in the case of a continuum formulation. For transport, the pore scale will consider molecular diffusion in the pore fluid while the continuum model includes dispersion as an upscaled parameter (Fig. 3). Reactive surface area is used in continuum treatments, while diffusion to and from explicit reacting mineral surfaces are normally considered in the pore scale models.

Contaminant transport

Modern reactive transport methods have made significant contributions to the topic of contaminant transport in recent years. The contributions can be broadly divided into two important themes: 1) the more comprehensive and rigorous treatment of geochemistry in contaminant transport models, and 2) the role of aquifer heterogeneity on transport and mixing rates. The second of these topics is not considered further here, since it is primarily the purview of hydrology and transport disciplines—we focus here on the geochemical and mineralogical aspects of the problem.

The starting point for contaminant hydrogeology has always been the linear distribution coefficient (or Kd) models for sorption, which have the advantage that they can be easily incorporated into the transport equations without rendering them nonlinear. It has also been asserted that they have the advantage that the data required for the implementation is minimal, or at least less than that required by the multicomponent models discussed below. However, it should be pointed out that the Kd models, while potentially being as simple as the determination of the sorbed and aqueous concentration under a set of environmental conditions, typically require unique constraints for each and every site to which they are applied. In addition, as will be apparent from the discussion below, such Kd models might not even apply to a single site where conditions (temperature, salinity, competing ion concentrations, sorption site density) change over time. In other words, there is no real generality in the case of linear distribution coefficients, in contrast to the more rigorous ion exchange and surface complexation models (like that of metals on iron hydroxides) that can be applied quite widely.

The primary issue with the Kd models, and even somewhat more advanced sorption models like the Langmuir or Freundlich formulations, is that they do not consider competitive sorption. In contrast, surface complexation or multicomponent ion exchange models explicitly consider competitive effects. This is of course well known to the geochemical community, so perhaps the failure to embrace the more rigorous multicomponent sorption models is hard to understand. Multicomponent ion exchange and surface complexation models are briefly reviewed below, and then considered in the context of reactive transport.

Ion exchange and transport

Ion exchange reactions can be described via a mass action expression with an associated equilibrium constant. The exchange reaction can be written in generic form, assuming here a chloride solution, as:

 
νAClu(aq)+uBXν(s)uBClν(aq)+νAXu(s)
(18)

where X refers to the exchange site occupied by the cations Au+ and Bv+. The equilibrium constant or selectivity coefficient, Keq, for this reaction can be written as

 
Keq=[BClν]u[AXu]ν[AClu]ν[BXν]u
(19)

where the parentheses [] refer to the thermodynamic activities. Here it is clear that the sorption of any particular contaminant (for example, 137Cs+ or 90Sr2+) will be affected by the competing cation concentrations, for example, sodium via the reaction:

 
Na++CsXNaX+Cs+
(20)

A strongly saline solution like that associated with the tank leaks at the Hanford site, arguably one of the most significant environmental problems facing the United States, will thus have an important role in determining the effective 137Cs linear distribution coefficient. The real-world problem can be made more complicated (and potentially worse) when multiple exchange sites occur in the sediment or soil, as is the case at the Hanford 200 tank farm (Zachara et al. 2002). In this case, the effective Kd depends also on the Cs+ concentration (Fig. 4), as noted by Zachara et al. (2002) and by Steefel et al. (2003).

An example 2D reactive transport simulation was presented in Steefel et al. (2005) to demonstrate the effect of the competing NaNO3 (a major component of the leaking tank fluids) concentration. This simulation was carried out using ion exchange selectivities and site concentrations determined from Hanford bulk sediment (Steefel et al. 2003). Figure 5 shows the transport of the non-reactive nitrate, with the plume extending over the entire domain size considered. The fractional migration of 137Cs at 1 M NaNO3 is significantly less, in keeping with the expectation that some retardation of this contaminant will occur (although still more than would be the case for 137Cs in a typically dilute soil or vadose zone water). At 5 M NaNO3, the migration of the 137Cs is significantly farther, even if still retarded relative to the nitrate. The enhanced 137Cs observed below many of the Hanford 200 tanks can thus be explained largely on the basis of a classical ion exchange model that accounts for the competing NaNO3 concentrations in the plume, and on the elevated concentrations of Cs+ that exceed the number of high affinity sites that are available for strong sorption. The model also suggests that as leaking of highly concentrated NaNO3 tank fluids ceases, further migration of the 137Cs is unlikely.

Surface complexation and transport

Perhaps an even more convincing demonstration of the inadequacy of the classical Kd models is provided by reactive transport analyses of metal and radionuclide migration influenced by surface complexation on iron hydroxides. In fact, an entire generation of geochemists have investigated the surface complexation behavior of metals on ferric hydroxides over many years (Dzombak and Morel 1990; Davis et al. 1998, 2004), but only more recently has the use of such models been demonstrated convincingly in reactive transport frameworks at the field scale.

Perhaps the first successful demonstration of the use of surface complexation models to describe field-scale reactive transport was presented by Davis and co-workers for the Naturita, Colorado uranium-contaminated site (Curtis et al. 2004, 2006; Davis et al. 2004). Davis et al. (2004) used both electrostatic and non-electrostatic surface complexation models to describe uranium sorption on the Naturita sediment, but finally settled on the non-electrostatic models because of their flexibility in treating natural and complex multi-mineralic sediments. In order to match the pH dependence of sorption in particular, it was necessary to use a three site SCM consisting of weak, strong, and very strong sites. Figure 6 shows the match with the experimental data at a variety of CO2 partial pressures, an important variable because of the strong competition between uranium carbonate complexes in solution and the surface complexes developed on the Naturita sediment surfaces (Hsi and Langmuir 1985; Prikryl et al. 2001; Davis et al. 2004). Calcium uranium carbonate complexes in solution can further reduce the sorption of uranium, especially Ca2UO2(CO3)3(aq) (Bernhard et al. 2001; Brooks et al. 2003; Fox et al. 2006).

Curtis et al. (2006) used the non-electrostatic surface complexation model presented in Davis et al. (2004) to investigate field-scale transport at the Naturita, Colorado site (Fig. 7). They demonstrated that the surface complexation model, when combined with realistic flow and transport parameters for the aquifer, accurately reproduced the available data. They also demonstrated (again) the inadequacy of a constant Kd model (Fig. 7).

A number of other studies have been carried out demonstrating the ability of the surface complexation models to describe field scale behavior, including those at the Rifle CO uranium contaminated site (Yabusaki et al. 2007, 2017; Zachara et al. 2013), and at the contaminated Savannah River site (Bea et al. 2013; Arora et al. 2018).

The laboratory–field rate discrepancy

The laboratory–field rate discrepancy has been a longstanding topic of discussion in the geochemical literature (White and Brantley 2003). The suggestion has been that field rates are three to five orders of magnitude slower than rates constrained in laboratory settings for what is argued to be the same reactive pathway, but is this really the case? What (if anything) does reactive transport analysis have to contribute to this debate?

Although it is not straightforward to quantify, it is a reasonable conclusion that with the use of a modern reactive transport model that considers multiple reactive minerals and perhaps most importantly the approach to equilibrium with respect to various phases, one can reduce this apparent discrepancy by roughly one order of magnitude. Beyond this reconciliation, it appears that careful application of a rigorous reactive transport analysis that considers detailed reaction mechanisms and physical/chemical heterogeneity can eliminate most or all of the remaining disparity. The sources of the discrepancy (which in fact may or may not be real) can be attributed to at least three possible factors:

  1. Geochemical effects associated with inhibiting constituents in solution, or with incongruent reactions and a nonlinear dependence on the Gibbs free energy of the pore fluid in contact with the reacting phases (Zhu et al. 2004; Maher et al. 2006).

  2. Effects of mineralogical heterogeneity on reactivity (Li et al. 2006) and physical heterogeneity operating at the pore scale (Molins et al. 2012; Beckingham et al. 2017).

  3. Effects of flow field heterogeneity leading to non-uniform fluid travel times through structurally complex systems, thus leading to an apparent reaction rate that is the flux-weighted average of these distinct travel times (Steefel and Maher 2009; Maher 2010, 2011).

Maher et al. (2009) carried out a systematic reactive transport analysis of Terrace 5 (226 ka) in the Santa Cruz weathering chronosequence and demonstrated that there was essentially no lab-field discrepancy when several factors were accounted for in their 1D model (Fig. 8).

An advantage of the Santa Cruz field site was that the mostly granitic unconsolidated sediment allowed for truly 1D downward flow of pore fluids, so that the common issue of lateral diversion of flow in low-permeability crystalline rocks does not arise and confound the interpretation of the effective fluid–rock ratio along a flow path. To simulate the weathering in Terrace 5 (the oldest at 226 ka), Maher et al. (2009) set up a 1D reactive transport model that included K-feldspar and albite and allowed for kaolinite precipitation. Smectite was present as a primary phase, but did not have a large impact on the simulation results except for the role the mineral played in cation exchange. Much of the weathering took place in the variably saturated zone where it was important to include the diffusion of reactive species (mainly CO2) in the gas phase.

The key to resolving the “discrepancy”, however, was found in the suggestions from Zhu and co-workers (Zhu et al. 2004) that slow clay precipitation might explain why the system did not remain far from equilibrium where the rates of dissolution for the feldspars are at their maximum. Slow clay precipitation in the model was implemented so as to match observed supersaturation with respect to kaolinite, and then combined with the nonlinear dependence of feldspar dissolution rates on the reaction affinity (or ΔG) observed in experiments by Hellmann and Tisserand (2006). With these two coupled effects, the rates of dissolution of the feldspar were much slower than what would be expected with either fast (equilibrium) clay precipitation or a linear TST mineral dissolution model. Thus, it was possible to closely match the observed mineral profiles after 226 ka (Fig. 9, upper row). In addition, the model also matched present day mineral saturation states calculated based on pore water chemistry (Fig. 9, lower row).

Incorporating isotope fractionation

The incorporation of isotope fractionation into modern reactive transport models holds the potential for a large impact on many Earth systems applications, and such isotope-enabled simulations are beginning to gain recognition (Druhan and Winnick 2019). The reasons for this potential are illustrated by the significance of isotopes as an analytical component of modern geochemistry, namely their sensitivity to tracking reaction and transport processes far beyond the precision with which such processes can be monitored using other methods (e.g., conventional chemical analysis). This is particularly true for isochemical (or nearly so) systems like incongruent mineral dissolution and precipitation, where isotopes may be the only method available to unravel the underlying mechanisms of fluid–mineral interaction (DePaolo 2011).

Reactive transport modeling in turn can also be used to evaluate some of the time-honored models of isotope geochemistry. For example, a distillation or Rayleigh model has been used to describe the change in isotopic ratio as a function of reaction progress in both open and closed systems (Druhan and Maher 2017) according to:

 
r=r0[f(α1)]
(21)

where r is the isotopic ratio, r0 is the original ratio, f is the fraction of the original reactant remaining, and α is the fractionation factor. As pointed out by van Breukelen and Prommer (2008), the Rayleigh model underestimates the fractionation in the case of heterogeneous flow systems. This conclusion was reinforced by Druhan and Maher (2017) in their systematic study. The Rayleigh model also overestimates the fractionation for systems involving sorption (van Breukelen and Prommer 2008). In fact, in any system where the fractionation factor is not constant, the Rayleigh model fails. This could be expected to be the case where kinetic isotope fractionation occurs during mineral precipitation along a flow path. If the precipitating mineral approaches equilibrium (where the kinetic isotope fractionation should disappear), then one expects the fractionation factor α to change continuously.

At a first level, modern reactive transport software (e.g., PHREEQc, CrunchTope, FLOTRAN, TOUGHREACT, see Steefel et al. 2015) are fully capable of treating isotopes as if they were separate chemical components (Druhan et al. this volume). For example, the rate laws for the reduction of total aqueous CrO5342 and CrO5242 can be written respectively as:

 
rtotal53=k53(CrO4253)total
(22)
 
rtotal52=k52(CrO4252)total
(23)

with the kinetic fractionation factor given by:

 
αk=k53k52
(24)

The formulation in Equations (22) and (23) can be adapted to microbially mediated reactions, as for the reduction of the isotopologues of SO3442 and SO3242 to 34HS and 32HS respectively (Druhan et al. 2012, 2014), but this requires a coupling of the rate expressions:

 
r34=BSRBμ34X34[CSO42CSO42+KSO42]FKacetateFT
(25)
 
r32=BSRBμ32X32[CSO42CSO42+KSO42]FKacetateFT
(26)

where BSRB is the biomass of sulfate-reducing bacteria, μ34 and μ32 are the maximum specific rates, and X34 and X32 are the mole fractions of SO3442 and SO3242, respectively. FKacetate is a Monod term for acetate and FT is a thermodynamic function (Jin and Bethke 2005; Dale et al. 2006). This formulation is simplified from the full expression given by Druhan et al. (2012) for the case where the half saturation constant, KSO42, for the two isotopologues is the same. The kinetic fractionation is captured by the use of slightly different maximum specific rate constants for the two isotopologues.

Using the formulations given in Equations (25) and (26), it is possible to track the isotopic evolution of the fluid and mineral phases as a function of space and time in sediment undergoing sulfate reduction due to electron donor injection (Williams et al. 2011). Druhan et al. (2014) simulated the time evolution of the effluent concentrations of iron and sulfate as well as the isotopologues of 34S and 32S as a result of sulfate reduction in a one meter column subject to a constant flow rate with high concentrations of the electron donor acetate (Fig. 10). The high Fe2+ concentrations are produced by the abiotic reaction of H2S with ferric hydroxide, which then acts also to titrate out the sulfide as precipitated FeS.

As discussed above, the treatment of kinetic fractionation of aqueous species is relatively straightforward. Handling precipitation of a mineral phase, as in the case of FeS precipitated in the Rifle example above, is more complex. Where mineral precipitation occurs, it is necessary to treat the resulting mineral phase as a solid solution, which has the effect of coupling the isotopologues through their activities in the solid. As an example, we consider the case of calcite precipitation with the isotopologues 44Ca and 40Ca (Druhan et al. 2013). The rates are given by:

 
R44=Acckb44X44([Ca44][CO32]KeqX441)
(27)
 
R40=Acckb40X40([Ca40][CO32]KeqX401)
(28)

Acc is the reactive surface area for the mineral calcite, 44kb and 40kb are the backward rate constants for precipitation, and 44X and 40X are the mole fractions (used to describe activity) for the two isotopologue mineral end-members, respectively. Keq is the equilibrium constant for the reaction, which if different in Equations (27) and (28), would result in equilibrium isotopic fractionation. A linear TST formulation has been assumed here for the sake of simplicity. Fractionation in more complex rate expressions is now being explored.

Note that the activities of the end-members effectively couple the mineral rate laws through the mole fraction:

 
X44=[CaCO344][CaCO344+CaCO340],X40=[CaCO340][CaCO344+CaCO340]
(29)

which is true even in the case where there is no kinetic fractionation associated with precipitation, as in the Rifle example above (Druhan et al. 2014) and in the Cr isotope benchmark problem (Wanner et al. 2015). However, when different rate constants are used for precipitation of the isotopologues, then kinetic fractionation occurs. In addition to resulting in kinetic fractionation during precipitation, the rate law proposed by Depaolo (2011) and implemented by Druhan et al. (2013) and Steefel et al. (2014) allows for back reaction of the isotopes over similar time scales. As shown in Figure 11, isotopic re-equilibration occurs after the system has reached chemical equilibrium. This behavior follows from the formulation of the TST rate law in DePaolo (2011) based on the principle of detailed balancing (Lasaga 1984), but whether this is realistic in most natural systems is another question. Any straightforward implementation of the DePaolo (2011) TST rate law results in kinetic fractionation during mineral dissolution, so one possibility is to require that no fractionation occurs during dissolution (an option now in the code CrunchTope).

In some cases back-reaction (or re-equilibration) may occur, but much more slowly than is predicted by the TST rate law. Such behavior is suggested by modeling of diagenetic processes at Site 984 in the Mid-Atlantic, where matching of the isotope profiles over the millions of years associated with sediment burial requires re-equilibration rates about 4 orders of magnitude slower than the bulk precipitation rates of calcite (Fig. 11, right) used to match the total Ca profile (Fig. 12, left) (Steefel et al. 2014).

Electrostatic effects on reactive transport

There are two developments over the last 20 years concerning electrostatic effects on reactive transport that are worth noting. The first is provided by the use of the Nernst–Planck equation, which deals rigorously with transport effects associated with the development of diffusion potentials as ions diffuse at different rates (Steefel et al. 2015). The second effect is associated with transport within the diffuse or electrical double layer (EDL) bordering charged surfaces along pores (Tournassat and Steefel 2019, this volume).

Nernst–Planck equation

The Nernst–Planck equation accounts for electrochemical migration associated with the development of a diffusion potential as ions diffuse at different rates (Steefel et al. 2015; Tournassat and Steefel 2015, 2019, this volume). The equation can be applied in both bulk (electrically neutral) solution, and as we shall see below and in Tournassat and Steefel (2019, this volume), within the electrical double layer. The Nernst–Planck equation can be derived from equation for the chemical potential mj (Steefel et al. 2015):

 
μj=μj0+RTln(γjCj)+zjFψ
(30)

where μj0 is the chemical potential at standard state, R is the gas constant, T is the absolute temperature, γj is the activity coefficient for the jth species, zj is its charge, F is Faraday's constant, and ψ is the electrical potential. Since the flux of an ion can be described in terms of the gradient of its chemical potential (Lasaga 1998):

 
Jj=DjCjRTμjx
(31)

we can derive the Nernst–Planck equation by substituting Equation (30) into Equation (31) to give:

 
Jj=DjCjziFRTDiCiψDjCjlnγj
(32)

Note that the first term on the right hand side of Equation (32) is the Fickian term used in most hydrogeological transport models. The last term on the right hand side only applies where a gradient in ionic strength is present (Tournassat and Steefel 2019, this volume). The second term on the right hand side is usually referred to as the diffusion potential since it is generated by the diffusion of charged ions at unequal rates. It goes to zero when all of the ion diffusion coefficients are the same (in which case Fick's Law is recovered if no ionic strength gradients are present).

The Nernst–Planck equation has now been implemented in a number of reactive transport codes, including FLOTRAN, CrunchFlow, MIN3P, and PHREEQc (Steefel et al. 2015). It was used to model ion profiles in marine sediments with upwelling hydrothermal fluids (Giambalvo et al. 2002), a geological system where matching solute profiles would be impossible with the use of a simple Fickian diffusion model.

A transient test problem to show the effects of the diffusion potential term in the Nernst–Planck equation was proposed by Lichtner (1998). In this problem, Na+ and Cl are the same concentration inside the domain and to the left of the Dirichlet boundary on the left side of the domain (the right hand boundary is treated as a no-flux boundary). So a simple application of Fick's Law would say there should be no flux of either Na+ or Cl (since ∂C/∂x = 0). But gradients in pH and NO3 develop because the concentrations outside (0.001 mM) and inside (0.1 mM) the domain are different for these ions. Through the effects of the Nernst–Planck equation (32), fluxes of both Na+ and Cl are created through the diffusion potential term (Fig. 13). The more rapid loss of H+ from the domain, which diffuses about 9 times faster than the NO3, results in a net influx of Na+ to compensate the charge. Cl is also removed from the internal domain to compensate for the higher out-flux of H+ compared to NO3.

Rolle et al. (2013) used the Nernst–Planck equation to model transverse diffusion/dispersion in a 2D flow field, and this problem became part of a benchmarking study of the Nernst–Planck equation described in Rasouli et al. (2015). Figure 14 shows the simulation results from CrunchFlow based on the problem described by Rolle et al. (2013), with injection of 0.29 mM KCl and MgCl2 from the left side (flow left to right) over an interval of 1 cm (the remainder of the injection ports, designed to produce 1D flow, contained deionized water). Note the greater spreading of the K+ plume resulting from its 2.8 times higher diffusivity than that of Mg2+.

Transport in the EDL

Electrostatic effects on surface complexation have been considered for many years, even if the non-electrostatic models like those employed by Davis and co-workers are now more commonly applied at the field scale. A new class of models has been developed that considers double layer or electrical double layer (EDL) screening of the charge of mineral surfaces bordering the pores. The EDL is often conceptually subdivided into two regions: the Stern layer located within the first water monolayers of the interface in which ions adsorb as inner-and outer-sphere surface complexes, and a diffuse layer (DL) located beyond the Stern layer in which a diffuse cloud of ions screens the remaining uncompensated surface charge (Tournassat and Steefel 2019, this volume). At infinite distance from the charged surface, the solution is electrically neutral and is described as bulk (or free) water, contrasting with the diffuse layer where electrical neutrality does not prevail. In the diffuse layer, water is considered to have properties similar to bulk water and the distribution of ions predicted by molecular dynamic simulations is often in close agreement with predictions made with the Poisson–Boltzmann (PB) equation (Tinnacher et al. 2016).

The Poisson–Boltzmann equation is derived by combining the Poisson equation that relates the local charge imbalance at a position y in the direction perpendicular to the charged surface to the second derivative of the electrostatic potential (ΨDL in V) at the same position:

 
d2ψDL(y)dy2=1εΣiziFCi,DL
(33)

where ε is the water dielectric constant (78.3 × 8.85419·10−12 F·m−1 at 298 K), F is the Faraday constant (96 485 C·mol−1), and zi is the charge of ion i. The concentration of species i in the DL, Ci,DL, is given by the Boltzmann equation:

 
ci,DL=ci,0exp(ziFψDL(y)RT)
(34)

where R is the gas constant (8.3145 J·mol−1·K−1), T is the temperature (K), and ci,0 is the concentration of species i in the bulk water at equilibrium with the diffuse layer (Tournassat and Steefel 2019, this volume).

There is currently interest in solving the Poisson–Boltzmann equation numerically, but to do so is challenging because of the fine discretization required (typically on the order of nanometers). In addition, the well-known analytical solutions to the Poisson–Boltzmann are inaccurate in the case of asymmetric electrolytes, or where the surface potential is high (Schoch et al. 2008; Kirby 2010). Since the Poisson–Boltzmann equation applies to equilibrium or static systems, a superior treatment is provided by the Poisson–Nernst–Planck (PNP) equation (Schoch et al. 2008; Kirby 2010), although the same issues with discretization arise. The PNP equation is implemented by combining the Poisson equation (33) with the Nernst–Planck equation (32). To the author's knowledge, there is no implementation of either the Poisson–Boltzmann or the Poisson–Nernst–Planck equation in any of the modern reactive transport simulators for natural waters (Steefel et al. 2015).

An alternative to the PNP that has been now widely applied is the mean electrostatic potential (MEP) model, essentially a dual continuum model in which electrically neutral bulk water constitutes one continuum and the EDL the second (Revil and Leroy 2004; Appelo and Wersin 2007; Appelo et al. 2010; Revil et al. 2011; Tournassat and Steefel 2015; Tournassat et al. 2016). This model averages ion concentrations in the diffuse layer Ci,DL¯ by scaling them to a mean electrostatic potential (MEP, ΨM) that applies to the diffuse layer volume (VDL):

 
Ci,DL¯=1VDLDLziCi,0exp(ziFψDL(x,y,z)RT)dxdydzCi,0exp(ziFψMRT)
(35)

The mean electrostatic potential value can be deduced from the charge balance in the diffuse layer:

 
ΣiziFCi,DL¯=ΣiziFCi,0exp(ziFψMRT)=QDL
(36)

where QDL (in C · m–3) is the charge that must be balanced in the diffuse layer (Tournassat and Steefel 2019, this volume).

The MEP model can be evaluated by comparing to an ingenious diffusion experiment conducted by Glaus et al. (2013). In this experiment, a gradient in concentration of NaClO4 from 0.5 M to 0.1 M was applied across a diffusion cell containing Na-montmorillonite. After equilibrating for two days, both reservoirs were spiked with the same concentration of the isotope 22Na. From Fick's Law alone, one would expect no change in the concentration of the 22Na in either reservoir, since no gradient in concentration is initially present. In fact, the experiments showed that the concentration of the 22Na increased continuously in the high concentration reservoir and decreased continuously in the low concentration reservoir. This was modeled with the MEP model implemented in CrunchClay (see Tournassat and Steefel 2019, this volume for figure showing the results).

Another recent demonstration of the MEP model is provided by a modeling study of the DR-A borehole diffusion experiment carried out in the Opalinus Clay in Switzerland (Soler et al. 2019). In this experiment, a cocktail of anions, cations, and non-reactive tracers were added to the borehole and their concentrations monitored over time as the cocktail solutes diffused out into the Opalinus Clay. Figure 15 shows the time evolution of two anions, iodide and bromide, as compared to the uncharged tritiated water (HTO). Note that over the same time period, the HTO drops to about 30% of its initial value, while the two anions decrease only to about 55% of their initial value. The MEP model captures the behavior of both the uncharged solute and the anions with a single model, which is not possible with a simple application of Fick's law. In addition, it captures accurately the phenomenon of anion exclusion (Soler et al. 2019).

Multiphase reactive systems

Arguably, one of the most significant achievements in reactive transport analysis in recent years has been the application to multiphase systems, including soils and the deep vadose zone where mostly water and air coexist (Dwivedi et al. 2019; Arora et al. 2019, both this volume), and the deep subsurface where other fluids and gases may be present (e.g., oil, gas, steam). It is well beyond the scope of this contribution to review these, but some discussion can be found in Sin and Corvisier (2019) and in Bildstein et al. (2019), both this volume. This is a challenging topic because of the presence of different phases, each of which has its own transport and chemical properties, and because of the complex, typically time-dependent interfaces between the phases.

Chemical–mechanical

Although some research has been initiated in recent years, particularly on the topic of pressure solution (Renard et al. 1997; Yasuhara et al. 2003, 2004, 2005; Yasuhara and Elsworth 2008; Taron and Elsworth 2010), and on coupled processes in nuclear waste respositories or geothermal systems (Taron et al. 2009; Gens et al. 2010; Taron and Elsworth 2010; Zheng et al. 2010; Rutqvist et al. 2014; Rutqvist 2015, 2017), much remains to be done to advance the topic together with the multicomponent chemistry we typically consider be a centerpiece of modern reactive transport analysis. Most of the analyses to date have been based on continuum models, but there is a new interest in incorporating chemistry into discontinuum mechanics model that can track the evolution of deforming and reacting interfaces explicitly (Hu et al. 2017).

Watersheds and beyond

Watershed reactive transport is another important frontier topic that is only now beginning to get the attention it deserves (Li 2019, this volume). Watersheds are situated in scale somewhere between the more local “Critical Zone” studies (e.g., Sullivan et al. 2019) and the larger regional and even global scale reactive transport models we can expect to see more of in the future (e.g., Goddéris et al. 2019). Progress is being made now in integrating reactive transport capabilities into the complex catchment scale hydrological models (Bao et al. 2017; Li et al. 2017; Zhi et al. 2019), which are complicated by the need to consider integrated surface and subsurface water flow, as well as the effects of both hydrological and biogeochemical hot spots (Dwivedi et al. 2018a,b).

CONCLUSIONS

The field of reactive transport has come a long way since the publication of Reviews in Mineralogy Vol. 34 in 1996. The modeling approach, which considers multicomponent chemical components and multiple reaction processes, has now been successfully applied in a large number of environments ranging from contaminated groundwater systems, hydrothermal and geothermal systems, nuclear waste repositories, chemical weathering, and now even watershed and regional scale river systems. Application of more sophisticated sorption models like those based on multicomponent ion exchange and surface complexation should have revolutionized the contaminant hydrogeology business, but so far the industry has resisted what appears to be the inevitable given the vastly superior performance of these models compared to the classical Kd model. Similarly, significant progress has been made in resolving the lab–field discrepancy using reactive transport modeling. The advent of pore scale models, which aim to capture explicitly the pore structure within which water–rock interaction occurs, is likely to provide further insights into the reactivity of natural soils, sediments, and rocks. The use of a new class of electrostatic models that consider the electrical double layer in charged media offers the promise of quantitatively explaining the migrations of ions in low permeability (especially clay) environments.

Acknowledgments

This work was supported by the Director, Office of Science, Basic Energy Sciences, Chemical Sciences, Geosciences, and Biosciences Division, of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231 to Lawrence Berkeley National Laboratory. This work was also supported as part of the Watershed Function Science Focus Area at Lawrence Berkeley National Laboratory funded by the U.S. Department of Energy, Office of Science, Biological and Environmental Research under the same contract. The author is grateful to the BRGM (Carnot project) and the University of Orleans (Labex Voltaire) in France for support for a research stay where at least some of the progress noted here was made, and where the view from the crossroads was first examined.

Open access: Article available to all readers online.