This is an Open Access article distributed under the terms of the Creative Commons Attribution License (, which permits unrestricted use, distribution and reproduction in any medium provided that the original work is properly attributed.


Heterogeneity in pore structure and reaction properties including grain size and mineralogy, pore size and connectivity, and sediment surface area and reactivity is a common phenomenon in subsurface materials. Heterogeneity affects transport, mixing, and the interactions of reactants that affect local and overall geochemical and biogeochemical reactions. Effective reaction rates can be orders of magnitude lower in heterogeneous porous media than those observed in well-mixed, homogeneous systems as a result of the pore-scale variability in physical, chemical, and biological properties, and the coupling of pore-scale surface reactions with mass-transport processes in heterogeneous materials.

Extensive research has been performed on surface reactions at the pore-scale to provide physicochemical insights on factors that control macroscopic reaction kinetics in porous media. Mineral dissolution and precipitation reactions have been frequently investigated to evaluate how intrinsic reaction rates and mass transfer control macroscopic reaction rates. Examples include the dissolution and/or precipitation of calcite (Bernard 2005; Li et al. 2008; Tartakovsky et al. 2008a; Flukiger and Bernard 2009; Luquot and Gouze 2009; Kang et al. 2010; Zhang et al. 2010a; Molins et al. 2012, 2014; Yoon et al. 2012; Steefel et al. 2013; Luquot et al. 2014), anorthite and kaolinite (Li et al. 2006, 2007), iron oxides (Pallud et al. 2010a,b; Raoof et al. 2013; Zhang et al. 2013a), and uranyl silicate and uraninite (Liu et al. 2006; Pearce et al. 2012). Adsorption and desorption at the pore-scale have been investigated to understand the effect of pore structure heterogeneity on reaction rates and rate scaling from the pore to macroscopic scales (Acharya et al. 2005; Zhang et al. 2008, 2010c, 2013b; Zhang and Lv 2009; Liu et al. 2013a). Microbially mediated reactions have also been studied at the pore-scale including denitrification (Raoof et al. 2013; Kessler et al. 2014), sulfate reduction (Raoof et al. 2013), organic matter and nutrient transformation (Knutson et al. 2007; Gharasoo et al. 2012; Raoof et al. 2013), and biomass growth (Dupin and McCarty 2000; Dupin et al. 2001; Nambi et al. 2003; Knutson et al. 2005; Zhang et al. 2010b; Tartakovsky et al. 2013). These studies indicate that pore-scale heterogeneity and coupling of reaction and transport have major impacts on the macroscopic manifestation of reactions and reaction rates.

Various experimental and numerical approaches have been developed and applied for pore-scale investigation of surface reactions in porous media. A micromodel is one of the experimental systems most widely used for studying pore-scale reactions under flow conditions. It is a 2-dimensional (2-D) flow cell system typically fabricated with silicon materials, into which a desired pore structure can be imprinted to form a pore network (Zhang et al. 2010a). The interfaces between pores and solids in the pore network can be coated with certain redox-sensitive materials such as hematite (Zhang et al. 2013a) to provide reactive sites for surface reactions. The micromodel allows real time observation of surface reactions at the pore-scale using a suite of spectroscopic and microscopic techniques including optical microscopy (Dupin and McCarty 2000; Nambi et al. 2003; Knutson et al. 2005; Zhang et al. 2010a; Yoon et al. 2012), synchrotron X-ray microprobe (XMP) and X-ray absorption spectroscopy (XAS) (Pearce et al. 2012), and Raman spectroscopy (Zhang et al. 2013a). Other experimental systems for pore-scale research include flow-cell reactors (Fredd and Fogler 1998; Panga et al. 2005; Pallud et al. 2010a,b), capillary tubes (Molins et al. 2014), and columns (Bernard 2005; Luquot and Gouze 2009; Fridjonsson et al. 2011; Noiriel et al. 2012; Luquot et al. 2014). These experimental apparatus complemented with non-destructive imaging techniques, such as X-ray microtomography (XMT) (Bernard 2005; Luquot and Gouze 2009; Navarre-Sitchler et al. 2009; Noiriel et al. 2012; Luquot et al. 2014; Molins et al. 2014) and nuclear magnetic resonance (NMR) (Fridjonsson et al. 2011), provide powerful tools for studying pore-scale surface reactions, reaction-induced pore structure and permeability changes, and feedback between surface reactivity and reactive transport in porous media (Dupin and McCarty 2000; Nambi et al. 2003; Knutson et al. 2005; Willingham et al. 2008; Zhang et al. 2010b; Edery et al. 2013).

Numerical approaches are an important component of pore-scale geochemical research. Pore-scale flow and reactive transport models are used to simulate the temporal and spatial evolution of surface reactions in coupling with transport, to analyze and interpret experimental results, and to complement experimental approaches by extending experimental conditions. Commonly used numerical approaches include pore-network models (Acharya 2005; Li et al. 2006, 2007; Mehmani et al. 2012; Raoof et al. 2013; Varloteaux et al. 2013), lattice Boltzmann models (Kang et al. 2006, 2010, 2014; Knutson et al. 2007; Lichtner and Kang 2007; Zhang et al. 2008; Huber et al. 2014), smoothed particle hydrodynamics models (Tartakovsky et al. 2008b), and conventional finite-difference or finite-volume methods (van Duijn and Pop 2004, 2005; Willingham et al. 2008; Orgogozo et al. 2010; Porta et al. 2012, 2013; Yoon et al. 2012; Liu et al. 2013a; Steefel et al. 2013; Molins et al. 2014; Trebotich et al. 2014).

These models can explicitly incorporate pore geometry, reactive surface area, and distribution, and even molecular reaction rates to simulate coupled transport and reactions at the local scale and to provide insights into effective rates at the macroscopic scale.

Both experimental and numerical studies indicate that effective reaction rates can decrease by orders of magnitude from the molecular to macroscopic scale, with the magnitude decrease dependent on specific reactions and the structures of porous media (Swoboda-Colberg and Drever 1993; Malmström et al. 2000; Acharya et al. 2005; Li et al. 2006; Meile and Tuncay 2006; Navarre-Sitchler and Brantley 2007; Zhang et al. 2008; Bi et al. 2009; Flukiger and Bernard 2009; Zhang and Lv 2009; Miller et al. 2010; Liu et al. 2013a; Zhang et al. 2013b; Raoof et al. 2013; Salehikhoo et al. 2013). Various factors have been identified that may individually or collectively contribute to the scale-dependent behavior of geochemical reaction rates (Alekseyev et al. 1997; White and Brantley 2003; Meile and Tuncay 2006; Lichtner and Kang 2007; Steefel and Maher 2009; Maher 2010; Miller et al. 2010; Dentz et al. 2011; Liu et al. 2013a, 2014a). These factors can be grouped into two categories:

  1. scale-dependent variations in mineral surface properties, such as surface roughness that controls reactive surface area and secondary mineralization that alters surface reactivity (White 1995; Alekseyev et al. 1997; Benner et al. 2002; White and Brantley 2003), and

  2. heterogeneous distributions of flow, transport, and chemical properties in subsurface environments (Malmström et al. 2004; Li et al. 2006; Meile and Tuncay 2006; Lichtner and Kang 2007; Maher 2010; Shang et al. 2011; Liu et al. 2013a, 2014a; Zhang et al. 2013a).

Previous research into the existence and cause of scale-dependent reaction rates has targeted specific examples with relatively constrained physical and chemical conditions. General theories and/or approaches for systematically describing and conceptually predicting the scaling effect are under-developed. In this chapter, we begin with the development of a theoretical base for predicting the scale-dependent behavior of effective surface reaction rates in heterogeneous porous media, including explicit linkage of the effective reaction rate constant with its intrinsic rate constant and the important influence of pore-scale variations in reactant concentrations (see the Section Theoretical Consideration of Effective Reaction Rates). Approaches to estimate intrinsic reaction rates and rate constants are then discussed with emphasis on molecular simulations (see Intrinsic Rates and Rate Constants). We follow with examples of how the coupling of pore-scale subsurface reaction and transport processes creates variations in reactant concentrations, and how these, in turn, cause the effective rate constants to deviate from the intrinsic reaction rate constant as a function of time and space in: i) grain-scale reactive diffusion systems (see Grain-scale reactions, sub-grain process coupling, and effective rates) and ii) in reactive transport systems with fluid flow (see Subgrid variations in reactant concentrations and effective rate constants under flow conditions). We conclude with a summary of the results and discussion in this chapter, a concept to minimize the scale-dependent behavior of effective surface reaction rates in macroscopic models through consideration of subgrid heterogeneity, and a brief discussion of other potential important factors affecting reaction rate scaling.


Geochemical reaction rates on mineral surfaces can be generally described using the mass action law:



where ri is the reaction rate, kins is the intrinsic rate constant, Ai is the reactive surface area, cj and γj are the concentration and activity coefficient of chemical species j, vij is the reaction order with respect to chemical species j. Rigorously, the mass action law (Eqn. 1) is only applicable at the pore scale, and concentrations and reaction rates in Equation (1) are defined locally at specific reactive surface locations or sites. As described in the following sections, the reactive surface locations and chemical species concentrations are generally heterogeneously distributed in subsurface materials. The observed and simulated chemical reaction rates in subsurface materials are effective reaction rates that are either explicitly or implicitly defined within a numerical grid cell volume that contains many pores and mineral surfaces. The grid cell volume can be the representative elemental volume (REV) to define Darcy-scale physical properties (Bear 1979), but most often it is a volume of porous media with its size constrained by computational resources and resolution need (Mo and Friedly 2000; Chiogna and Bellin 2013). Grid cells may be selected at different scales depending on applications (Fig. 1) including:

  1. Earth system models (Harrison et al. 2009; Poudel et al. 2013) where an entire watershed may represent one numerical grid cell;

  2. field scale numerical problems with numerical grid sizes ranging from meters to hundred meters (Scheibe et al. 2006; Maher 2010; Molins et al. 2010; Li et al. 2011; Bao et al. 2014);

  3. core scale numerical simulations with numerical grid size from centimeter to meters (Zinn and Harvey 2003; Shang et al. 2011, 2014; Liu et al. 2014a); and

  4. pore-scale research at the scales of μm to mm (Willingham et al. 2008; Zhang et al. 2013a).

Subgrid heterogeneity ranging from pore to large scales is a common feature of all these applications, indicating that Equation (1) is not directly applicable to describe effective reaction rates in heterogeneous subsurface materials.

To derive the effective reaction rate, Equation (1) is averaged over a porous medium volume (ΔV) to yield:



where the variables with a bar are the average variables over the pore space within ΔV. Using the average variables, pore-scale variables in Equation (1) can be written as the follows:







where Ai′, cj′, and γj′ are the deviations from their corresponding average values within ΔV. Note that the average of Ai′, cj′, or γj′ within ΔV is zero. Replacing Equations (3–5) into Equation (2) yields:



In subsurface systems, the relative deviation of the activity coefficient from its average value is typically small (Lasaga 1979). However, the relative deviations of pore-scale reactive surface area and chemical species concentrations from their average values within ΔV can be significant in heterogeneous materials. Equation (6) reveals that the effective reaction rate defined for ΔV is not only a function of average surface area and chemical species concentrations, but also depends on the pore-scale variations of these variables within ΔV. The calculation of the pore-scale variations in Equation (6) requires detailed knowledge of pore structure, pore-scale chemical species concentrations, and reactive surface area or site concentration distributions. Such information is difficult to obtain in natural subsurface environments. Consequently, the pore-scale variations in Equation (6) are lumped into an effective rate constant:



and the effective reaction rate in ΔV can be expressed using the average variables:



Equation (8) has the same form as Equation (1). However, the variables and rate parameters are different. Equation (8) equals Equation (1) if the spatial variability of all variables is minimal in ΔV, such as in a well-mixed system. Well-mixed systems are rare, and consequently Equation (8) is generally different from Equation (1). However, Equation (8) is used interchangeably with Equation (1) in most subsurface reactive transport studies and simulations. In these cases, an effective rate constant must be used instead of the intrinsic rate constant.

The effective rate constant is influenced by the relative deviations of pore-scale concentrations from their average values within ΔV (Eqn. 7). The relative deviation of reactive surface area depends on reactive mineral distribution and accessibility, as well as geochemical/ biogeochemical reactions that act on the surfaces and their potential to alter surface reactivity as reactions proceed. The concentrations of dissolved species will depend on not only surface reactions, but also transport processes that supply reactants at the pore scale and remove reaction products. Consequently, the effective rate constant (keff) will be strongly affected by the coupling of reaction and transport, and will generally vary with space and time. Next we present several generic scenarios to demonstrate potential changes of keff from its intrinsic value. In sections that follow, we discuss reactive transport studies to derive effective rate constants in heterogeneous porous media.

Well-mixed conditions

Under well-mixed conditions, the spatial variations in dissolved species concentrations are negligible. Although reactive surface area or reactive surface site concentration can be spatially variable, the effective reaction rate constant is the same as the intrinsic reaction rate constant:



Mass transport limited conditions

For the convenience of demonstration, we first consider a case with only one dissolved species limiting the rate of a surface reaction (i.e., j = 1 and νij = 1 in Eqn. 1). Equation (7) can then be simplified after ignoring the pore-scale variation in activity coefficient:



The pore-scale distribution of reactive area (or site concentration) is often negatively correlated with that of dissolved reactant concentrations under mass transport limited conditions in porous media. For example, finer-grained materials typically have a higher surface area and more reactive sites. These finer materials can form lenses, aggregates, coatings, or spots, which are less accessible to advective fluid due to the limitation of local hydraulic conductivity or pore connectivity. Consequently these regions may have lower concentrations of dissolved reactants supplied through fluid flow. Under this condition, Equation (10) predicts a lower keff than kins, because the second term on the right hand site of Equation (10) is a negative value. With increasing size of ΔV, large scale hydraulic heterogeneity, such as macroscopic flow channels, can affect the spatial distributions of reactants. This may further increase the range of dissolved reactant concentrations within ΔV, and thus increase the negative magnitude of the second term on the right hand side of Equation (10). The net result is that the effective rate constant decreases with increasing scale (domain I in Fig. 2). Equation (10), however, also predicts an alternative scenario when the variations in reactive surface site and dissolved reactant concentration are positively correlated. When this happens, keff will be larger than kins (domain II in Fig. 2). This scenario requires that more reactive sites reside at locations where dissolved reactants are rapidly supplied.

More complex scenarios can develop when two or more chemical reactants limit the rate of a surface reaction, or when a reaction is not first order with respect to a reactant. For example, when a second-order rate expression is considered with respect to dissolved species j (i.e., j = 1, νi1 =2 in Eqn. 1), then



The bracketed value in Equation (11) not only depends on the pore-scale correlation between the variations in reactive surface site and dissolved reactant concentrations, but also on the relative strength of the self-correlation for variations in reactant concentration. A large, self-correlation in the variations of dissolved reactant concentration can decrease the effect of the negative correlation between dissolved species and reactive surface site concentrations, leading to less scale-dependent behavior. Similar complex scenarios can exist for reactions involving multiple reactants. For these cases, the effective rate constants in Figure 2 may not uniformly decrease or increase as scale increases.


The effective rate constant is proportional to the intrinsic rate constant (Eqn. 7). Well-mixed experimental systems such as batch and stirred flow-cell reactors are often used to characterize reaction mechanisms and to estimate intrinsic rate constants that are independent of transport processes. Recently, molecular-scale simulations have been used to estimate kinetic rates and intrinsic rate constants. In this section, the most commonly used molecular simulation methods to calculate reaction rates are briefly described. The application of these methods for determining molecular-scale rates and rate constants of surface reactions are then discussed using several surface reaction examples.

Approaches to calculate molecular-scale reaction rates

Molecular dynamics (MD) is a powerful technique to study the molecular-scale behavior of geochemical systems (Stack et al. 2013). However, because the integration time step in MD simulations is constrained by the timescale for molecular vibrations that is usually on the order of 1 femtosecond, the simulated time that can be achieved routinely with MD simulations is limited to tens to hundreds of nanoseconds. Therefore, a range of techniques, referred to as rare-event techniques hereafter, have emerged with the goal of extending the timescale of MD simulations to investigate kinetic progress with rates that are lower than those manageable with standard MD algorithms. Examples of rare-event methods include transition path sampling (Dellago et al. 1998), metadynamics (Laio and Parrinello 2002), and reactive flux (Chandler 1987). A frequent application of the reactive flux approach to geochemically relevant systems is the study of water exchange reactions around aqueous ions (Rey and Guardia 1992; Rey and Hynes 1996a,b; Spångberg et al. 2003; Kerisit and Parker 2004b; Loeffler et al. 2006; Rustad and Stack 2006; Stack and Rustad 2007; Kerisit and Rosso 2009; Dang and Annapureddy 2013; Kerisit and Liu 2013). This same approach was subsequently applied to quantify the rates of surface complexation of aqueous ions on mineral surfaces (Kerisit and Parker 2004b) and eventually to look at attachment and detachment reactions at kink sites to simulate the elementary steps of mineral dissolution and growth (Stack et al. 2012).

In the reactive flux approach, the reaction rate is written as the product of the transition state rate constant, kTST, and the transmission coefficient, κ. The former corresponds to the rate at which the system reaches the transition state and is defined as



where μ is the reduced mass of the atoms involved in the reaction, β is 1/kBT where kB is the Boltzmann constant and T is the temperature, δ is the reaction coordinate, δ* is the value of the reaction coordinate at the transition state, and W(δ) is the free energy of the system at δ, which can be obtained using a number of techniques, including the potential of mean force approach (PMF). In the PMF approach, separate MD simulations are performed with the system constrained at specific values of δ, and the mean force along the reaction coordinate is integrated to obtain W(δ). The transmission coefficient accounts for the departure from ideal transition state behavior due to barrier re-crossing (i.e., not all trajectories that reach the transition state are successful transitions between the reactants and products free energy wells) and is determined from the plateau value of the normalized reactive flux, which can be computed as



where θ[x] is the Heaviside function, which is 1 if x is larger than 0 and 0 otherwise, and δ(0) is the initial velocity along the reaction coordinate. The subscript ‘c’ denotes that the initial configurations are generated in the constrained reaction coordinate ensemble. The reactive flux is usually computed by running a large number of short trajectories forward and backward in time from initial configurations constrained at the transition state.

Molecular-scale rates of uranyl sorption reactions at mineral surface sites

Molecular simulation techniques have been used to study uranyl sorption at mineral surfaces (Kremleva et al. 2010; Geckeis et al. 2013). Mineral surfaces have included silica (Greathouse et al. 2002; Patsahan and Holovko 2007; Boily and Rosso 2011), alumina (Moskaleva et al. 2004; Glezakou and deJong 2011; Tan et al. 2013), iron and titanium oxides (Steele et al. 2002; Perron et al. 2006a,b, 2008; Drot et al. 2007; Sherman et al. 2008; Roques et al. 2009; Skomurski et al. 2011; Pan et al. 2012; Sebbari et al. 2012), phyllosilicates (Zaidan et al. 2003; Greathouse and Cygan 2005, 2006; Greathouse et al. 2005; Kremleva et al. 2008, 2011, 2012; Veilly et al. 2008; Hattori et al. 2009; Martorell et al. 2010; Lectez et al. 2012; Liu et al. 2013b; Yang and Zaoui 2013a,b; Teich-McGoldrick et al. 2014); feldspars (Kerisit and Liu 2012, 2014), and carbonates (Doudou et al. 2012). The majority of these studies focused on energy minimizations of adsorbed uranyl complexes at various mineral surfaces. As a whole, they have shown that significant variations exist in surface complex configurations and binding energies for different surfaces and sites of a given mineral. Although adsorption rates were not calculated in most cases, this body of work strongly suggests that molecular-scale heterogeneities should influence the adsorption/desorption rates of uranyl for a given mineral. For example, static density functional theory (DFT) calculations of uranyl surface complexation on the basal and edge surfaces of kaolinite (Al2Si2O5(OH)4) (Kremleva et al. 2008, 2011; Martorell et al. 2010) revealed significant differences between the tetrahedral Si layer and the octahedral Al layer, as well as a wide range of binding energies on different sites om the (010) edge surface. Similar conclusions were drawn from simulations of uranyl surface complexation on pyrophyllite (Al2Si4O10(OH)2) edge surfaces (Kremleva et al. 2012), where it was concluded that multiple adsorbed complexes likely exist on any given surface, and that different edge surfaces exhibit different preferred surface complex coordination environments.

A small subset of these studies investigated the kinetics of uranyl sorption. For example, MD simulations by Boily and Rosso (2011) using PMF calculations showed significantly different barrier heights for uranyl adsorption/desorption on three quartz (α-SiO2) surfaces, providing further evidence that sorption rates are heterogeneous at the molecular scale even for a given mineral. MD simulations of uranyl sorption on the most stable surfaces of orthoclase (KAlSi3O8 -tectosilicate) and kaolinite (phyllosilicate) (Kerisit and Liu 2012, 2014) revealed large quantitative differences in the free energy barriers for adsorption/desorption as well as in the affinity of uranyl for these two types of silicate minerals. The formation of uranyl carbonate complexes affected desorption free energy barriers on the two minerals to different extents, indicating that variations among mineral phases can be compounded by solution conditions that impact aqueous speciation. Desorption rates of adsorbed uranyl orthoclase were calculated to be on the order of 106 hr−1 (Kerisit and Liu 2012), a rate several orders of magnitude faster than the grain-scale rate constant (Liu et al. 2013a).

Molecular-scale rates of elementary mechanisms of mineral growth and dissolution

Initial applications of rare-event techniques such as reactive flux to calculate molecular-scale rates of mechanisms relevant to mineral growth and dissolution focused on adsorption of ions on ideal, planar surfaces (i.e., formation of adatoms) (Piana et al. 2006; Kerisit and Parker 2004a,b). For example, adsorption of alkaline-earth cations on the {101̄4} calcite surface was found to be very fast (~108 s−1) and desorption was calculated to be only marginally slower (~106–107 s1) (Kerisit and Parker 2004b), demonstrating the presence of labile species on flat surfaces. The exact values of these rates as well as the extent of the driving force for adsorption is dependent on the details of the force-field parameters (Raiteri and Gale 2010). Similar free-energy barriers for adsorption/desorption were obtained by Piana et al. (2006) for Ba2+ on a pure, flat {001} barite (BaSO4) surface. In contrast, detachment rates of Ba2+ from a step edge on the same surface were determined by Stack et al. (2012) to be much slower (~104 s−1). Therefore, orders of magnitude differences in adsorption and desorption rates can exist between different sites and surface structures on the same mineral surface.

The work of Piana et al. (2006) on Ba2+ sorption on three barite surfaces also demonstrated that the free-energy barrier for adsorption is strongly surface specific and can be significantly reduced by the presence of anions adsorbed on the surface, indicating that the sequence of sorption events plays an important role in determining the overall reaction rate. MD simulations of Ca2+ and CO32− dissolution from calcite surfaces also resulted in free-energy barriers that were heavily surface specific (Spagnoli et al. 2006). In addition, detailed rare-event simulations of Ba2+ attachment and detachment to/from a stepped surface by Stack et al. (2012) highlighted that, even for this seemingly simple surface reaction, multiple intermediates were involved with a wide range of transition rates between intermediates (104 to 1010 s−1). Rare-event techniques have also improved understanding of mineral nucleation. For example, MD simulations demonstrated importance of precursor amorphous clusters in CaCO3 homogenous nucleation by demonstrating the presence of small free energy barriers and favorable reaction free energies for the addition of CaCO3 ion pairs to a growing nucleus (Tribello et al. 2009; Raiteri et al. 2010).

The studies described above made use of classical MD techniques due to the large size of the simulated systems. Another approach to computationally probe the rates of mineral growth and dissolution involves the use of small clusters that are treated quantum mechanically to represent surface sites and that can be immersed in a continuum dielectric model to simulate solvation effects. With this approach, potential energy surfaces can be obtained by progressively elongating bonds within the clusters to evaluate the energetics of mineral dissolution. Although the effects of the extended crystal structure are lost with this approach, it allows for the making/breaking of covalent bonds and investigating the effects of pH. In particular, this approach has been applied, for example, to study the dissolution of silicate minerals such as olivine (Morrow et al. 2010) and feldspar (Criscenti et al. 2005, 2006).


Subsurface materials contain microscopic pores and/or fractures within individual minerals, grains, and grain-aggregates (Anbeek 1993; Berkowitz and Scher 1998; Lee et al. 1998; White et al. 2001; Berkowitz 2002; Zachara et al. 2004; Liu et al. 2006; Luquot and Gouze 2009; Hay et al. 2011). These intra-granular porous domains have pore sizes ranging from nanometers to tens of micrometers, and have a high interfacial surface area to pore volume ratio (Anbeek 1993; Hurlimann et al. 1994). Various geochemical reactions can occur in the intra-granular pore domains that affect bulk solution chemistry in external pore regions (Barber Ii et al. 1992; McGreavy et al. 1992; Anbeek 1993; Lee et al. 1998; Liu et al. 2006; McKinley et al. 2006, 2007). The transport in these intra-granular domains is dominated by diffusion that is strongly affected by pore connectivity (Ewing et al. 2012). Slow diffusion process allows intra-granular solutions to evolve close to mineral equilibrium that can deviate from that in bulk solutions. Geochemical reactions such as sorption/desorption (Ball and Roberts 1991; Qafoku et al. 2005; Liu et al. 2008; Ginn 2009) and mineral precipitation and dissolution (McGreavy et al. 1992; Lee et al. 1998; Anderson et al. 2002; Liu et al. 2006; McKinley et al. 2006, 2007; Pallud et al. 2010a) can occur in defiance of thermodynamic conditions in bulk solutions. The slow diffusion, on the other hand, decreases mass exchange rate with external pore waters, leading to non-ideal behavior of chemical transport even for non-reactive species (Beven and Germann 1982; Ball and Roberts 1991; Ewing et al. 2010, 2012; Hay et al. 2011).

Micromodel experimental systems have been used to evaluate the relationship between the intrinsic and effective rate constants (Zhang et al. 2013a; Liu et al. 2015). Micromodels can be fabricated to simulate the complex physical structure of natural sediment grains and grain-aggregates. Figure 3 shows an example of micromodel systems. The four corner regions in the micromodel (Fig. 3) are composed of grain aggregates and intra-aggregate pores. Surrounded by the aggregates is a macropore domain where advection dominates. The micromodel system has been used to provide insights into the scale-dependent variation in the effective rate of uranyl desorption that is over 3–5 orders of magnitude slower than the intrinsic rate estimated from molecular simulations (Liu et al. 2013a). It has also been used to examine hematite reduction in heterogeneous media under flow conditions and to investigate the difference in effective rate constants between well-mixed and heterogeneous transport systems (Zhang et al. 2013a).

These examples display the importance of the coupling of transport and reaction at the pore-scale in controlling the grain-scale reaction rates. Accurate simulations of grain-scale reaction rates are, however, challenged by the difficulties in the estimation of effective reaction rates in both intra- and inter-granular domains (Liu et al. 2006; Zhang et al. 2013a). We now demonstrate how a pore-scale reactive diffusion model can be used to simulate the variability in reactant concentrations in intra-granular domains. The simulated concentrations are then used by Equation (7) to calculate effective rate constants for grain-scale reactive transport models. Grain-scale models based on mass transfer and diffusion will be compared to assess how different physical models influence the deviation between effective and intrinsic rate constants.

Pore-scale variability in reactant concentrations at the sub-grain scale

The evolution of reactant concentrations in intra-granular domains can be described using a reactive diffusion equation:


cjt=·(Djcj)+i=1Nrajiri,         j=1,2,N,

where Di is the molecular diffusion coefficient of aqueous chemical species j, aji is the stoichiometric coefficient of species j in reaction i, Nr and N are the total numbers of kinetic reactions and species in the system, respectively, and the other variables are defined before for Equation (1). Various numerical approaches exist for solving the multi-species diffusion problem described by Equation (14) using implicit or explicit numerical scheme, or their combination (Yeh and Tripathi 1991; Press et al. 1992; Liu et al. 2011).

The reactive diffusion equation (Eqn. 14) was used to simulate pore-scale variations in reactive surface site and dissolved species concentrations during hematite reduction in the intra-granular pore spaces as shown in the four corner regions in the micromodel (Fig. 3). In this example, hematite is assumed to be uniformly distributed on all the pore surfaces (bottom and side walls) in both intra- and inter-granular domains in the micromodel. The experimental design, micromodel fabrication, and hematite deposition on the pore surfaces for this example have been detailed elsewhere (Zhang et al. 2013a). To initiate hematite reduction in the intra-granular domain for the modeling purpose, reduced flavin mononucleotide (FMNH2) with a constant concentration (100 μM) was provided in the external pore space (center white region in Fig. 3). FMNH2 is a biogenic organic molecule capable of reducing iron oxides (Marsili et al. 2008; Shi et al. 2012). Driven by the concentration gradient, FMNH2 diffuses into the intra-granular domains where it reacts with hematite on pore surfaces. Hematite (Fe2O3) undergoes reductive dissolution to yield Fe(II) and FMN, which is an oxidized form of FMNH2. The kinetics of hematite reduction in a well-mixed system follows the rate expression (Liu et al. 2007; Shi et al. 2012; Zhang et al. 2013a):



where cFe2O3 is the pore-scale hematite concentration on the mineral surfaces (μM), cFMNH2 is the pore-scale aqueous concentration (μM), and kins is the rate constant at the pore-scale (μM−1h−1), which is equivalent to the intrinsic rate constant for this study (Eqn. 7). Equation (14) and Equation (15) were solved using the finite-volume method with a constant numerical node size of 10 × 10 μm2. The initial hematite concentration on each pore surface is 0.73 μmol/cm2 in the micromodel. The surface based hematite concentration was converted to volumetric concentration using a surface area of 10 × 28 μm2 for the side wall and 10 × 10 μm2 for the bottom surface wall.

The statistical distributions of FMNH2 and residual hematite concentrations in the intra-aggregate domains were calculated to vary with time (Fig. 4). Initially, FMNH2 concentrations were zero at all diffusion nodes (Fig. 4A). Initial hematite concentration was one of five possible values after normalizing to local pore node volume (5 red bars in Fig. 4B), a consequence of the fact that each pore node can contain 1 bottom pore wall, and additional zero to four side walls of hematite depending on the topology of water-grain interfaces in a grid node. The corresponding five possible hematite concentrations were 0.26, 0.99, 1.72, 2.45, and 3.18 M. The high molar concentration of hematite resulted from a high surface to pore volume ratio at the pore-scale. After 23 hours of FMNH2 diffusion and hematite reduction, the FMNH2 concentration distribution shifted toward higher concentrations (Fig. 4A), while the hematite concentration distribution shifted toward lower concentration (Fig. 4B). The hematite concentrations in nodes with four side pore walls did not change, because those nodes were not connected to external pore regions. Nodes with the larger decrease in hematite concentrations were those containing higher concentrations of FMNH2, which lead, in turn, to the negative correlation between hematite and FMNH2 concentrations at later times. This correlation affected the effective rate constant in the intra-granular domain as discussed in the next subsection.

Effective reaction rates and rate constants

The intra-granular diffusion domain is often termed as the immobile domain and external pore regions where water flows is often termed as the mobile domain in hydrology. The mass transfer between mobile and immobile domains is often described using the first-order mass exchange equation (Brusseau et al. 1992; Haggerty and Gorelick 1995; Šimunek et al. 2003):



where cjim is the average concentration of species j in the immobile domain, cjim is the average concentration in the mobile domain, and km is the mass transfer coefficient. In this model, the spatial change in species concentrations in the immobile domain is no longer considered. Using this model to include reactions in the immobile domain, Equation (16) becomes:



where riim is the average rate of reaction i in the immobile domain. Using terminology for Equation (7), the reaction rate in Equation (17) can be expressed as the follows for the case of hematite reduction in the micromodel (Fig. 4),



From Equation (7), we have



The ratio of the effective rate constant to intrinsic rate constant can be calculated using the pore-scale concentration results in Figure 4. The calculated effective rate constant (Fig. 5) was lower than the intrinsic rate constant and decreased with increasing time. This was because residual hematite concentration was increasingly negatively correlated with the FMNH2 concentration distribution in the diffusion domain as hematite reductively dissolved. Note that the simulation was only run for 23 hours in this case, and the effective rate constant is expected to further decrease with time, before it eventually increases back to 1 when all hematite is reduced and dissolves.

Diffusion models, such as one-dimensional slab or spherical diffusion models (Crank 1975; Ball and Roberts 1991) represent an alternative approach to describe intra-granular reactions and mass exchange with the mobile domain:



where Da is the effective diffusivity. The diffusion model is macroscopic without considering details in pore structure and connectivity in the intra-granular domain. However, Figure 6 shows that the effective reaction rate constants become closer to the intrinsic rate constant using the diffusion model. In this specific case, a slab reactive diffusion model was used, and the intra-granular macroscopic diffusion domain was divided into 3 numerical nodes (exterior, which is close to external pore space; interior, which is the farthest from external pore space; and middle, which is between exterior and interior nodes). The effective reaction rate constant in each node was calculated using the simulated pore-scale concentrations of FMNH2 and residual hematite within a specific numerical node. The calculated effective reaction rate constants still deviate from the intrinsic value in each numerical node with a largest deviation observed in the exterior node (Fig. 5). The deviation correlates with the concentration change that was the greatest in the exterior node. It is expected that the effective rate constant in each node will become closer to the intrinsic rate constant when more numerical nodes are used. However, complete elimination of the deviation between the effective and intrinsic rate constants cannot be achieved regardless of the number of numerical nodes used in discretizing Equation 20, because of the macroscopic nature of the grain-scale diffusion model.

The effective rate constant also changed with numerical grid size when a macroscopic model is used to describe the reaction in the intra-granular domain (Fig. 6). Figure 6 only shows the results for two types of grain-scale grid node sizes:

  1. the intra-granular domain was divided into 3 numerical nodes (grid node size = 700 mm, subgrain scale), and

  2. the intra-granular domain was treated as one numerical grid node (grid node size = 2000 mm, grain scale) in corresponding to the results in Figure 5.

The effective rate constant generally decreased with increasing grid size (Fig. 6); however, its value varied significantly at each grid size scale, reflecting spatial and temporal changes in pore-scale concentrations of reactants and their correlation (Fig. 4). The effective rate constants in Figure 6 were calculated using the pore-scale concentrations at time 5, 15, and 23 hours of FMNH2 diffusion and hematite reduction. It is expected to have a wider range of effective rate constant distribution when other time results are also included.


Subsurface materials are heterogeneous in many ways (Fig. 1) in addition to intra-granular complexity in pore structure and reactivity (Figs. 3 and 4). A numerical grid system used for simulating large scale reactive transport typically contains heterogeneity in both intra- and inter-granular domains within individual grid nodes. This heterogeneity is referred to as subgrid heterogeneity, which has been attributed to 5–7 orders of magnitude decrease in the effective rate of U(VI) desorption in sediments from the molecular to field scales (Liu et al. 2013a, 2014a; Ma et al. 2014), and four orders of magnitude decrease for microbial reduction of uranium from laboratory to field systems (Bao et al. 2014). Here the micromodel system (Fig. 3) was used again to demonstrate how to:

  1. calculate the variations in reactant concentrations in systems with subgrid heterogeneity, and

  2. use the calculated results to estimate effective rate constant as a function of space and time.

Macroscopic reactive transport models with and without considering subgrid heterogeneity are compared to evaluate their effect on the estimated effective rate constants.

Pore-scale concentration variations under flow conditions

Simulations of pore-scale reactive transport under flow conditions require calculating fluid flow and chemical reactive transport. Fluid flow at the pore-scale can be generally described using the Navier–Stokes (N-S) equation (Bird et al. 2007).



where u is the velocity vector, g is the gravitational constant, g is the gravitational force, h is the pressure or water head, and ν is the effective kinematic viscosity (ν = μ/ρ, where μ is the dynamic viscosity and ρ is the fluid density). The simulation of fluid flow using Equation (21) requires knowledge of the pore geometry. However, soil and sediment systems are heterogeneous, consisting of various sizes of pores, minerals, organic matter, and organisms, often exhibiting hierarchical structures spanning seven or more order of magnitude. Pore size distribution can be quantified to the nm scale by direct and indirect methods including mercury porosimetry and neutron scattering, but the important pore network is much more difficult to define, especially at scale of < 1μm (Ball et al. 1990; Hay et al. 2011). This creates a scenario that pore geometry for large pores can be spatially resolved, while that for small pores is below spatial resolution. The soil and sediment domains containing small pores will therefore have to be treated as porous media (mixed pores and solids with unknown pore geometry) in numerical discretization. This will prevent the direct application of Equation (21) to describe fluid flow in soils and sediments. An alternative model has recently been developed to describe fluid flow in mixed media containing both pores and porous domains (Yang et al. 2014). The model becomes N–S in the pore regions, and becomes a Darcy-law-based model for fluid flow in the porous domain.

When the fluid velocity field is specified, the reactive solute transport can be simulated using the pore-scale multi-species advection and diffusion equations coupled with reactions:



Equations (21) and (22) were used to simulate fluid flow, hematite reduction, and reactive transport in the micromodel (Fig. 3) to calculate reactant variations with time. FMNH2 (100 μM) was injected from the left side of the micromodel with a constant flow rate (270 μL/h). The injected FMNH2 reacts with hematite on pore surfaces in both inter- and intra-granular domains. Hematite was reduced and dissolved quickly in the advection domain (red region) as FMNH2 carried by fluid from the left side reacted with the hematite on the bottom pore wall in the advection domain (Fig. 7). There were regions near the grain–grain contacts in the inter-granular pore domain (center up and bottom, Fig. 7), where hematite reduction was much slower than that in the advection domain. These regions mimic the wedge spaces or dead end pores created by grain-grain contacts in porous media. This slow reduction region is termed as the macropore diffusion domain, as the advection bypasses these regions (Zhang et al. 2013a; Liu et al. 2015).

The hematite-reduction front migrated from the advection domain to the macropore and intra-granular diffusion domains through FMNH2 diffusion. Figure 7 only shows one snapshot. The hematite-reduction fronts moved further into the diffusion domains with continual FMNH2 injection. Eventually all hematite in the macropore diffusion domain was reduced. The reduction rate in the intra-granular domain was, however, much slower as discussed previously. Note that the initial FMNH2 concentration was zero and initial hematite concentration was one of 5 possible values depending on topology of water–grain interfaces in the grid node. After 23 hours of reduction, the FMNH2 concentration in most pores in the advection domain was equal to the FMNH2 concentration in the injection solution (= 100 μM); but its concentration was still zero in the interior of the intra-granular domain. The net result is a wide distribution of FMNH2 concentration (Fig. 8). The average concentration of FMNH2 increased in the micromodel, while the residual hematite concentration decreased with time (dashed lines in Fig. 8). Those nodes with a higher FMNH2 concentration have a lower concentration of hematite than their average concentrations, leading to a negative value for the second right-hand term in Equation (7).

Effective rate constants

The effective overall rate constant for the micromodel (calculated using the results in Fig. 8 and Eqn. 7) decreased by two orders of magnitude during the first 23 hours of the hematite–FMNH2 reaction (Fig. 9A). Ratio of keff/kins decreased immediately from 1 to about 0.75 within a few seconds of FMNH2 injection because of a fast hematite reduction in the advection domain. After that, the effective rate decreased steadily with time. The effective rate constant for the micromodel was lower than the effective rate for the intra-granular domain (Fig. 6), as a result of larger concentration variations in the micromodel than those in the intra-granular domain only. The effluent concentration of Fe(II) produced from hematite reduction can be predicted using the effective rate for the micromodel under the assumption that the micromodel was a well-mixed system (Fig. 9B). By comparison, the model significantly over-predicted effluent Fe(II) if an intrinsic rate constant was used in the well-mixed model.

The above analysis of the effective rate constant treats the entire micromodel as a well-mixed reactor, or as one numerical grid volume in a reactive transport systems. The estimated effective rate constant reflects the pore-scale deviations in FMNH2 and hematite concentrations from their average values in the micromodel pore network (Fig. 6). Effective rate constants can also be calculated for each subgrid domain (advection domain, macropore diffusion domain, and intra-granular domain, Fig. 9A). The effective rate constants within individual subgrid domains deviate much less from the intrinsic value. The result implied that macroscopic models would improve the prediction by explicitly considering subgrid reactive domains if effective rate constants are unknown, and the intrinsic rate constant is used for model prediction in heterogeneous materials. However, Figure 9A demonstrates that the effective rate constants for individual domains still deviate from the intrinsic value. This deviation is expected to decrease with the application of more subgrid domains, but it cannot be completely eliminated. Model prediction using the intrinsic rate constant will always contain errors regardless how many subgrid domains are considered.


Geochemical and biogeochemical reactions occur fundamentally at the pore scale, and they are coupled and modified by transport processes in subsurface materials. These coupled pore-scale processes are, however, below the spatial resolution of most core and field scale observations. Consequently solute and solid phase concentrations and reaction rates observed at the macroscopic scales reflect the manifestation of coupled pore-scale reactions and transport that in most cases cannot be explicitly separated. Pore-scale reactive transport models can theoretically be used to simulate pore-scale processes and their coupling. These models, however, require pore-scale properties that are difficult to obtain for practical problems at the macroscopic scales. Consequently, medium-based reactive transport models with porous media reaction and transport properties defined on a numerical grid volume of porous media have to be adopted based on the mass conservation law. The theoretical analysis and results in this chapter demonstrate that reaction rates and rate constants defined on a numerical grid volume in heterogeneous materials are generally different from the intrinsic reaction rate and rate constants. The effective rate constants for the grid volume depend on the statistical distributions of reactants and their spatial correlation, which are expected to vary with both temporal and spatial scales, causing scale-dependent behavior of effective reaction rates and rate constants. Importantly, the deviation of the effective rate constants from the intrinsic value can be reduced if a macroscopic model considers reactive transport in subgrid domains. The subgrid domains may be derived from non-reactive tracer transport behavior that can be independently measured for a subsurface system at the appropriate scales of interest. However, subgrid domain models cannot completely eliminate the deviation of the effective rate constants from the intrinsic value, and thus contain uncertainty in model prediction.

Other factors beyond pore-scale process coupling can also affect the scaling of reaction rates from the molecular to macroscopic scales in heterogeneous subsurface materials. For example, the reaction rate at the molecular scale can be heterogeneous, even on a single mineral phase (Piana et al. 2006; Spagnoli et al. 2006; Boily and Rosso 2011; Lüttge et al. 2013); the rates of reactions occurring in nano-pores can be influenced by water ordering, surface curvature, and overlapping double layers (Kerisit and Liu 2009, 2012; Argyris et al. 2010; Bourg and Steefel 2012); and mineral surface reactivity can be modified as reactions proceed (Alekseyev et al. 1997; Benner et al. 2002; White and Brantley 2003; Bernard 2005; Noiriel et al. 2005; Luquot and Gouze 2009). The pore geometry and connectivity can also be changed by mineral precipitation (Steefel and Lichtner 1994; Davis et al. 2006; Zhang et al. 2010a; Yoon et al. 2012; Fanizza et al. 2013; Boyd et al. 2014) and dissolution (Noiriel et al. 2004, 2005, 2009; Bernard 2005; Flukiger and Bernard 2009; Luquot and Gouze 2009; Kang et al. 2010), microbial growth (Taylor and Jaffe 1990; Nambi et al. 2003; Thullner et al. 2004, 2007; Knutson et al. 2005; Seifert and Engesgaard 2007, 2012; Brovelli et al. 2009; Zhang et al. 2010b; Kirk et al. 2012; Tang et al. 2013) and biomineralization (Zhang and Klapper 2010; Cunningham et al. 2011; Pearce et al. 2012; Phillips et al. 2012; Hommel et al. 2013). These and other factors can be superimposed on each other. The effective rate equation (Eqn. 6), the effective rate constants (Eqn. 7), and the numerical results and data presented in this chapter represent a starting point for a systematic analysis of scale-dependent behavior of geochemical and biogeochemical reactions. Given the ubiquity of heterogeneity in subsurface materials, understanding, predicting, and minimizing the deviation of the effective reaction rate and rate constant from the intrinsic rate and rate constant is critically important in applying mechanism-based kinetic parameters for calibrating and predicting larger scale reactions and reactive transport.


This research is supported by the U.S. DOE, Office of Science, Biological and Environmental Research (BER) as part of the Subsurface Biogeochemical Research (SBR) Program through Pacific Northwest National Laboratory (PNNL) SBR Science Focus Area (SFA) Research Project. This research was performed using Environmental Molecular Science Laboratory (EMSL), a DOE Office of Science user facility sponsored by the DOE’s Office of BER and located at PNNL. PNNL is operated for DOE by Battelle Memorial Institute under contract DE-AC05-76RL01830.