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.


Darcy-scale simulation of geochemical reactive transport has proven to be a useful tool to gain mechanistic understanding of the evolution of the subsurface environment under natural or human-induced conditions. At this scale, however, the porous medium is typically conceptualized as a continuum with bulk parameters that characterize its physical and chemical properties based on the assumption that all phases coexist in each point in space. In contrast, the pore scale can be defined as the largest spatial scale at which it is possible to distinguish the different fluid and solid phases that make up natural subsurface materials. Because the pore scale directly accounts for the pore-space architecture within which mineral reactions, microbial interactions and multi-component transport play out, it can help explain biogeochemical behavior that is not understood or predicted by considering smaller or larger scales (Fig. 1). Specifically, the nonlinear interaction between the coupled physical and geochemical processes may result in emergent behavior, including changes in permeability, diffusivity, and reactivity that is not captured easily by a Darcy-scale continuum description.

Reactive processes in porous media such as microbially mediated reduction–oxidation (Fig. 1) or mineral dissolution–precipitation (Fig. 2) take place at interfaces between fluid and solid phases. Because the different phases are distinguishable at the pore scale, experimental and modeling studies need to consider these interfaces so as to accurately determine reaction rates. An interface is the surface between two phases that differ in their physical state or chemical composition. Depending on the scale of observation, the appearance of the interface can vary. Sharp interfaces are those in which the physical and chemical characteristics change abruptly across the interface. Diffuse interfaces are those in which the characteristics change smoothly over a layer of varying thickness. Reactive processes themselves can change the appearance of the interface. For example, mineral heterogeneity can result in the creation of degraded zones (Fig. 2), where dissolution of a faster-dissolving mineral (e.g., calcite) from within a matrix of relatively insoluble minerals (e.g., dolomite and silicates) leaves behind a porous continuum Deng et al. (2013).

As advances in experimental and imaging techniques allow for improved characterization of pore-scale processes, modeling approaches are being challenged to incorporate the textural and mineralogical heterogeneity of natural porous media, in particular with regard to their treatment of interfaces. Simulation of the evolution of reactive interfaces is critical to capture the processes that lead to emergent behavior, including the reactive infiltration instability (Ortoleva et al. 1987; Hoefner and Fogler 1988; Steefel and Lasaga 1990) or reactivity evolution (Luquot and Gouze 2009; Noiriel et al. 2009). The focus of this chapter is on the direct numerical simulation of reactive processes at the pore scale, with an emphasis on the role of fluid–solid interfaces (Kang et al. 2007; Tartakovsky et al. 2007). Direct numerical simulation (DNS) employs mesh-based methods that often require an explicit representation of these reactive interfaces. Many concepts discussed in this chapter, however, are applicable more generally, including to other modeling approaches such as the lattice Boltzmann or particle methods (Kang et al. 2007; Tartakovsky et al. 2007).

In this chapter, the equations for pore-scale processes of flow, transport and geochemical reactions are succinctly presented, followed by a description of the methods for their direct numerical simulation within pore-scale domains. Next, the representation of reactive surfaces in DNS applications is reviewed, with an emphasis on surface reactivity evolution and transport limitations to reactive surfaces. This is followed by a review of methods to capture the physical evolution of the pore space with a focus on reactive instabilities. To end the chapter, an approach for upscaling interfacial reactions from the pore scale to the Darcy scale is presented.



Given the low compressibility of water, the incompressible Navier–Stokes equations can be used to accurately describe the flow of water in the pore space via the conservation of momentum and mass, respectively, for most subsurface conditions:





where the left-hand side of Equation (1) describes the inertial forces, and the right-hand side includes the pressure gradient (∇p) and the viscous forces, with u being the fluid velocity, ρ the fluid density, and μ the dynamic viscosity. Inertial forces have their origin in the convective acceleration of the fluid as it flows through the tortuous pore space. Viscous forces originate in the friction between water molecules and are responsible for the dissipation of energy. When viscous forces dominate (i.e., at very small Reynolds number) the above equations can be simplified to the steady-state Stokes equations (Steefel et al. 2013):





Equations (1–2) and (3–4) are to be solved within the pore space occupied by the fluid phase and delineated by the fluid–solid interfaces (Fig. 1).

Multicomponent reactive transport

Transport and reaction of dissolved species in the aqueous phase can be described by the following conservation equation:



where ci is the molal concentration of species i in solution (mol kg−1H2O), MH2O is the mass fraction of water (kg H2O kg−1), Di is the diffusion coefficient (m2 s−1) and Ri is the reaction term (mol m3 s−1). In Equation (5), transport takes place by advection (the translation in space of dissolved or suspended material at the rate of movement of a bulk fluid phase) and diffusion (mixing of solutes in the multicomponent mixture driven by concentration gradients). For simplicity, electrochemical migration associated with diffusion of charged species (Steefel et al. 2013) has not been included in Equation (5), although its contribution may be significant under acidic conditions (Molins et al. 2012; Ovaysi and Piri 2013). The reaction term Ri is expressed in volumetric terms and includes homogeneous reactions such as aqueous complexation. Aqueous complexation reactions are typically fast enough that they can be considered as locally at equilibrium, and Equation (5) can be written in the canonical form in terms of components (Steefel et al. 2014). Interfacial reactions such as mineral dissolution–precipitation take place at solid–fluid boundaries and are not included in Equation (5).

Surface reactions

Interfacial reactions can instead be expressed as a boundary condition at the solid–fluid boundaries



where rm is the surface reaction rate (expressed in units of mass per unit time per unit surface) and ξim is the stoichiometric coefficient of the i-th component in each surface reaction m.

A wide range of rate expressions can be considered. For mineral dissolution–precipitation reactions, the transition state theory law is often employed (e.g., Molins et al. 2014)



where km is the intrinsic rate constant (mol m−2 reactive surface s−1), Ea is the activation energy (kcal mol−1), Πain is a product representing the inhibition or catalysis of the reaction by various ions in solution raised to the power n, with ai being the activity of species I, ΔG is the Gibbs free energy (kcal), with m1, m2, and m3 being three parameters that affect the affinity dependence, R is the ideal gas constant (kcal K−1 mol−1), and T is the temperature (K). In the case of far-from-equilibrium dissolution, the affinity or ΔG term can often be neglected, as in the studies of chlorite dissolution under high-pCO2 conditions (Smith et al. 2013).

Other rate expressions have been used for interfacial reactions. For example, Wood et al. (2007) used Michaelis–Menten kinetics for the enzyme-mediated reaction rate applicable at the fluid–solid interface, assuming microbial cells were uniformly distributed on the surface (Fig. 1.I):



where km is the surface reaction rate, which includes the effect of the enzyme or biomass present on the mineral surface (mol m−2 reactive surface s−1), and Km is the half-saturation constant associated with the redox species i that defines the kinetic rate. In addition to kinetic models, equilibrium models, which do not require an explicit calculation of the rate, have also been employed at the pore scale. For example, Zaretskiy et al. (2012) used an equilibrium Langmuir model to calculate surface concentrations of adsorption species (si mol m−2 surface) as a function of species concentrations:



where smax and β represent the maximum surface concentration and the equilibrium constant, respectively (Zaretskiy et al. 2012).


Direct numerical simulation involves the use of conventional discretization methods to solve the flow, transport, and geochemical equations. These include Eulerian, mesh-based methods such as finite-difference, finite-element, or finite-volume methods, in which the differential equations are discretized by defining the value of the unknowns at fixed points in space. Finite differences have been applied to the simulation of flow and transport in sandstones (Sallès et al. 1993; Bijeljic et al. 2011a,b; Blunt et al. 2013) and carbonates (Bijeljic et al. 2013). The finite-element method has been used for the simulation of flow and transport in sphere packs (Cardenas 2008, 2009) and as a mixed finite-element–finite-volume method for flow and reactive transport in Fontainebleau sandstone (Zaretskiy et al. 2012). A finite-volume method has been applied for flow and reactive transport in computer-generated and experimentally derived simulation domains (Molins et al. 2012, 2014). Direct numerical simulation has also been used in combination with other approaches. For example, Yoon et al. (2012) used finite-volume methods to simulate reactive transport with a flow solution obtained with a lattice Boltzmann method, and a random-walk method was used for reactive transport by Sadhukhan et al. (2012) while a direct finite-difference solver was employed for Navier–Stokes flow. Both structured (e.g., Molins et al. 2014) and unstructured (e.g., Zaretskiy et al. 2012) meshes have been used in these applications.

To perform direct numerical simulation at the pore scale, experimental images of porous media resolved at the micrometer scale (e.g., from X-ray computed microtomography or XCMT) need to be converted to computational domains. Typically, segmentation is used to identify the discrete materials in an image (Wildenschild and Sheppard 2013). A number of techniques are used to segment images to obtain the morphology and topology of the pore space (Wildenschild and Sheppard 2013), frequently by binarization (e.g., Noiriel et al. 2009) but more recently using ternary segmentation methods (e.g., Deng et al. 2013; Scheibe et al. 2015). Binary reconstructed domains consist of voxels that are classified as either pore space or solid, while ternary reconstructed domains may allow for assigning a porosity value to areas of the domain that are under resolved by the characterization method (Scheibe et al. 2015) or may have been subject to reaction (Deng et al. 2013).

Interface representation

Binary domains can be directly incorporated into numerical models. In this scenario, fluid–solid interfaces can take the shape of a staircase in the case of structured meshes, where complex surfaces are represented as perpendicular walls locally at the grid cell level. Most structured-grid applications to simulation of flow and transport use this approach (e.g., Sallès et al. 1993; Bijeljic et al. 2011a,b; Blunt et al. 2013). However, for reactive applications the approach can lead to an overestimation of the actual interfacial area available for surface reaction (Fig. 3b). Finite-element and finite-volume unstructured methods are not affected by this issue because they can capture arbitrary surfaces by appropriate meshing strategies (e.g., Cardenas 2008, 2009; Zaretskiy et al. 2012). However, this advantage may be lost when the mineral surface evolves due to dissolution and precipitation, as re-meshing may be required.

To remedy the issue of reactive surface overestimation for a structured mesh, Molins et al. (2011) used the marching cubes algorithm to evaluate the interfacial area specifically for surface reactions while still solving for flow and transport within a binary, staircase domain. Another approach, the embedded boundary method, uses cut cells that result from intersecting the irregular fluid–solid interfaces with the structured Cartesian grid (Fig. 3a) to discretize the equations using the finite-volume method (Trebotich et al. 2014; Trebotich and Graves 2015). The embedded boundary method is thus a more rigorous approach in that flow, transport, and reactions are solved within the same domain. This method also makes it possible to leverage the advantages of structured methods while capturing complex surficial geometries (Molins et al. 2014). In Molins et al. (2014), a level set was used on the segmented microtomographic image to obtain the implicit function representing the calcite surface on the Cartesian grid. With a level set, the fluid–solid interface, Γ(x), is represented by a contour of a function φ(x,t) such that:



where c is a constant, and the level-set function φ is greater than c for one phase, and less than c for the other phase. The so-called level-set method builds on this description of the interface to track its motion (in particular, as a result of dissolution–precipitation reactions) and is discussed in the section devoted to pore-space evolution.

Micro-continuum and multiscale approaches

Binary representations of porous media as discussed above may not be sufficient where reactive processes result in degraded porous zones due to mineral heterogeneity (Deng et al. 2013) (see Fig. 2) or in porous precipitates that allow for diffusion (Yoon et al. 2012). The porosity in these reacted porous regions may be under-resolved at a given pore-scale model resolution; thus, the interface can be better described as discontinuous and diffuse rather than continuous and sharp. In fact, one could argue that binary representations do not represent natural porous materials where porosity is distributed—as a fractal or otherwise—over a very large range of spatial scales (Anovitz et al. 2013), or where reactivity may be preferentially found in connected nanoporosity, e.g., within chlorite in the Cranfield sandstone (Landrot et al. 2012).

The micro-continuum approach assumes the existence of a porous medium continuum at very small spatial scales. This assumption is valid over porous volumes in which the properties of the medium are continuously distributed (region II in Fig. 4). Darcy’s law is the governing equation for flow in the micro-continuum, while reactive transport is described by the mass balance equation of each species. Micro-continuum equations are parameterized using bulk parameters such as porosity and reactive surface area. For more details, the reader is kindly referred to the chapter of this volume devoted to the micro-continuum approach (Steefel et al. 2015, this volume).

A well-established example of the use of the continuum approach for interfacial reactive processes is that of reactive transport in fractures (Steefel and Lichtner 1998a,b; Noiriel et al. 2007). In fractured media, there is a sharp contrast—with a clear separation of scales—between the porosity in the fracture and in the rock matrix. As a result, the fracture is modeled as a fast flow path, where transport is dominated by advection, and the rock matrix, where transport is dominated by diffusion (Fig. 5). Heterogeneous reactions take place within the porous rock matrix or the layer of precipitate that coats the surface (Fig. 5). To simulate flow in discrete fractures in Darcy-scale models, the cubic law is often used to obtain a fracture permeability under the assumption of parallel smooth walls, e.g., Steefel and Lichtner (1998a,b) and Noiriel et al. (2007). However, this model does not typically provide an accurate estimation of fracture hydrodynamic properties at realistic fracture roughness, e.g., Deng et al. (2013).

Pore-scale flow in multiscale domains, such as fractured or vuggy media, can instead be simulated using the Stokes–Darcy or the Stokes–Brinkman equations (Golfier et al. 2002; Popov et al. 2009; Gulbransen et al. 2010; Yang et al. 2014). In the Stokes–Darcy approach, Darcy’s law and mass conservation are applied in the porous subdomains, and the Stokes equations in the free-flow subdomains. To close the model, mass- and momentum-conservation equations are specified at the interfaces between domains (Popov et al. 2009). When the geometry of these interfaces is very complex, however, the Stokes–Brinkman model is advantageous, in that a single set of equations is used over the entire domain. In the Stokes–Brinkman model, the Navier–Stokes equations (Yang et al. 2014) or the Stokes equations (Golfier et al. 2002; Popov et al. 2009; Gulbransen et al. 2010) are modified to add a Darcy term, e.g.:





where V is the effective velocity, K is the permeability tensor, P is the pressure, and μ̃ is an effective viscosity. The uppercase notation indicates here that these quantities apply to a porous medium continuum, where both fluid and solid phases are assumed to coexist at each point in space, rather to each discrete phase separately as in a pore-scale description. The process of obtaining these quantities from the microscale (pore-scale) description, formally referred to as upscaling, is discussed in a separate section below. In the pore space, permeability tends to infinity, the Darcy term becomes negligible and Equations (11–12) simplify to the Stokes equations (Eqn. 3–4). In the porous regions, on the other hand, both viscous and inertial terms become negligible due to slow velocities and Darcy’s law governs fluid flow (Popov et al. 2009; Yang et al. 2014). A very similar approach has been recently used by Scheibe et al. (2015) to simulate flow and tracer transport in a soil column imaged by XCMT and segmented using a ternary approach making it possible to explain the breakthrough curves observed experimentally.

For reactive transport problems, multiscale hybrid (Battiato et al. 2011; Roubinet and Tartakovsky 2013) and mortar methods (Mehmani et al. 2012) have been proposed. In these, Darcy-scale and pore-scale models are used in different subdomains, with the appropriate continuity of mass fluxes being enforced at the interfaces between subdomains (Battiato et al. 2011; Mehmani et al. 2012; Roubinet and Tartakovsky 2013). While these studies have simulated the macroscale treating the porous medium as a continuum (i.e., at the Darcy-scale, region II in Fig. 4) and the microscale with a pore-scale representation (region I in Fig. 4), the methods are general enough that they could be applied in multiscale problems in which a pore-scale representation were used for the macroscale (inhomogeneous region III in Fig. 4) and a micro-continuum representation for unresolved nanoscale porous regions (region II in Fig. 4). In this sense, they would be conceptually equivalent to the Stokes–Darcy approach for multiscale flow.


Mineral rates measured in laboratory experiments are often several orders of magnitude faster than those estimated from natural systems, e.g., Malmström et al. (2000). These differences in rates have been attributed to a variety of factors including, among others, reactive surface area accessibility in natural porous media (Peters 2009; Landrot et al. 2012), limitations on flow and transport in heterogeneous material (Drever and Clow 1995; Malmström et al. 2000; Salehikhoo et al. 2013; Li et al. 2014), or transport, rather than interface control of rates (Drever and Clow 1995; Steefel and Lichtner 1998b). Pore-scale modeling can be used to address some of these hypotheses by explicitly accounting for the rate-limiting effect of transport, and by incorporating mechanistic descriptions for the evolution of reactive surface area.

Transport control on rates

Transport to mineral surfaces, especially in physically heterogeneous media, can lead to poorly mixed conditions at the pore scale. As a result, effective rates cannot be determined using an average concentration in the pore space. Instead concentrations need to be resolved within the pore space, and rates calculated at each mineral surface. The pore-scale description can thus be used to evaluate the impact on the effective rates of the departure from the assumption of well-mixed conditions.

Molins et al. (2014) used a combination of experimental, imaging, and modeling techniques to investigate the pore-scale transport and surface reaction controls on calcite dissolution under elevated-pCO2 conditions. The laboratory experiment consisted of the injection of a solution at 4-bar pCO2 into a capillary tube packed with crushed calcite. A high-resolution, pore-scale, direct numerical model was used to simulate the experiment based on a computational domain consisting of reactive calcite, pore space, and the capillary wall, constructed from volumetric X-ray microtomography images. Part of the reported discrepancy between simulated Darcy-scale and pore-scale effluent concentrations was apparently due to mass-transport limitations to and from reactive surfaces. These were most pronounced near the inlet where larger diffusive boundary layers formed around grains. Transport limitations resulted from the heterogeneity of the pore structure: in slow-flowing pore spaces that exchanged mass by diffusion with fast flow paths, the assumption of well-mixed conditions did not apply (Fig. 6). Although the difference between pore- and Darcy-scale results due to transport controls was discernible with the highly accurate methods employed, this difference is expected to be more significant in more heterogeneous porous media, such as in natural subsurface materials.

Surface area evolution

In the study of Molins et al. (2014) it was assumed that the surface and pore-space geometry did not evolve over the 16 s of simulated time. However, as a result of reactive processes, normally both the reactive surface area and pore structure evolves. Here, two studies are presented that considered evolution of reactive surface area using a micro-continuum approach. Methods to track pore-space evolution are introduced in the next section.

In a study by Noiriel et al. (2009) of limestone infiltrated by CO2-rich solution, Sr and Ca release rates were used to assess the relative dissolution rates of the sparitic and microcrystalline phases in the limestone subjected to infiltration of CO2-rich solution. The results demonstrated that the reactive surface area (RSA) of the sparite increased greatly, as recorded by the rate of dissolution of that phase over time. In contrast, its geometric surface area, as recorded by XCMT, decreased slightly. To describe the time-dependent behavior, Noiriel et al. (2009) proposed a “sugar cube” model in which disaggregation of the granular network (presumably resulting in the large increase in RSA of the sparitic phase, which is now exposed to more of the reactive infiltrating solution) precedes dissolution of the individual grains of sparite (Fig. 7). Noiriel et al. (2009) described mathematically the time-dependent RSA, Ar, with the expression



where Ar0 is the initial surface area, Arm is the maximum surface area given by the sum of all of the surface areas of the individual particles, φ is the concentration of calcite over time, φ0 is its initial concentration, and n1, n2, and n3 are empirical coefficients that depend on the geometry of the aggregate.

Dissolution rates can also be affected by the precipitation of secondary mineral or transformation products, as it can limit access to the surface of the dissolving mineral. This is the case of the carbonation of wollastonite (CaSiO3) studied by Daval et al. (2009):



Experimental results showed that under circumneutral pH conditions dissolution rates were inhibited by the formation of a dense calcite coating. The passivation effect was successfully accounted for by the use of an effective reactive surface area



This equation links the true reactive surface area of wollastonite to the amounts of wollastonite and neo-formed calcite (φCaCO3), with a0 being the specific surface area in units of mol m−2, kp is a proportionality factor that can be obtained assuming that when t→∞, Aeff(∞) = 0, and p = 2 / 3 for monodisperse spherical particles, and p ≥ 2 / 3 for particle size distributions with a non-spherical geometry. Although the model slightly overestimated the rate for intermediate values of the reaction extent, the overall behavior of the extent of carbonation as a function of time was reproduced closely (Fig. 8). This showed that secondary calcite precipitation can play an important role as coatings of reactive surfaces (Daval et al. 2009).


Dissolution and precipitation of minerals and biogeochemical transformations such as bacterial growth change the geometry of the pore space. Although this evolution takes place at each fluid–mineral interface, the resulting changes have impacts at larger scales. For example, (bio)precipitates can significantly reduce the macroscopic porosity and permeability by plugging pore throats, e.g., Oelkers et al. (2008), Englert et al. (2009). Similarly, dissolution driven by injection of CO2 in the subsurface leads to an increase in porosity and permeability as well as changes in the mineral reactive surface area (Luquot and Gouze 2009).

Mathematically, solute dissolution and/or precipitation can be formulated as a moving boundary, or Stefan problem. Assuming uniform dissolution–precipitation of a single-mineral solid phase (m), the velocity of the moving solid–fluid interface (unΓ) can be described by, e.g.,



where Vm is the molar volume of the mineral. Equation (16) needs to be solved along with Equations (1, 2, 5, and 6) or Equations (3–6).

The front-tracking (Glimm et al. 1998), volume-of-fluid (Hirt and Nichols 1981), and level-set (Osher and Sethian 1988) methods have been used to solve moving boundary problems for tracking or capturing sharp interfaces. The advantage of these methods is that one can perform direct numerical simulations involving surfaces using the Eulerian approach on fixed Cartesian grids without having to parameterize these surfaces. Diffuse interface models, such as the phase-field method (Langer 1986), in which the interface is captured with a continuous variation of an order parameter, rather than explicitly represented as a sharp boundary, can also be used to solve the moving boundary problem. The level-set method and the phase-field method are described below in the context of recent applications to modeling pore-scale dissolution and precipitation processes.

Level set method

In the level-set method, the fluid–solid interface (Γ(x)) is represented by a contour level of the level-set function φ(x,t), Equation (10). The evolution of the level set is governed by the following advection equation, which describes the motion of the fluid–solid interface (Li et al. 2010):



where n is the norm of the level-set function (n = ∇φ/∇φ) at φ = c.

Li et al. (2010) used the level-set method for capturing interface evolution in simulations of dissolution and precipitation in a three-dimensional, single-pore throat (e.g., Fig. 9). The simulations were performed for a single-component system, where the reaction rate was described with both a zero-order kinetic term and an affinity term to account for near-equilibrium conditions. The three-dimensional effects of flow conditions and reaction rates were explored quantitatively. The simulation showed that the evolution of the permeability–porosity relationship depended on particular parameter values used (i.e., flow rates, rate constants). Further, the empirical Carman–Kozeny constitutive model did not capture evolution of permeability–porosity, especially for the conditions that led to non-uniform dissolution or precipitation patterns (Li et al. 2010).

Phase-field method

Phase-field methods are based on the idea that the free energy depend on an order parameter (φ, the phase-field variable) that acts as a function indicating what phase a point in space is in. The method replaces the boundary conditions at the interface (e.g., Eqn. 6) with a partial differential equation for the evolution of the phase field. The field changes smoothly in a diffuse zone around the interface that has a finite width.

A sharp-interface asymptotic analysis of the phase-field equations was developed by Xu and Meakin (2008) for the precipitation–dissolution problem considering diffusion as the only transport process. In this model, the evolution of the concentration field is a function of the evolution of the phase field:



where α is proportional to the molar volume of the mineral. The evolution of the phase field is described with



where τ is the phase-field characteristic time parameter, ɛ is closely related to the interface thickness, and λ controls the strength of the coupling between the phase field φ and the concentration c. The relationship between τ, ɛ, and λ is the n obtained from an asymptotic analysis, which ensures that the phase-field equations converge to the sharp interface solution:



Two-dimensional dendritic growth due to solute precipitation was simulated using this phase-field model (Xu and Meakin 2011). The simulations were performed under diffusion-limited conditions by setting the chemical reaction rate much larger than the rate of diffusion— i.e., the kinetics of the reaction at the interface was not relevant. The Mullins–Sekerka instability (Mullins and Sekerka 1964)—the same responsible for the dendritic growth of the snowflake—caused the growth process of a small circular nucleus placed in the center of a square domain to be very sensitive to perturbations (Xu and Meakin 2011). In the phase-field simulations of Xu and Meakin (2011), the perturbation was provided by the discretization of Equations (18–20) on a square lattice. As a result, diffusion-limited precipitation was observed to take place as an unstable dendritic growth (Fig. 10). The resulting solid–fluid interfacial pattern displayed a fractal geometry, whose fractal dimension agreed well with that estimated independently with a diffusion aggregation model (Meakin and Deutch 1984).

Continuum and multiscale approaches

As apparent from the results presented previously, the application of interface tracking methods to direct numerical simulation of pore-scale processes has been limited to relatively small, ideal problems. Instead direct simulation of pore-space evolution has often relied on treating the porous medium as a continuum at the Darcy scale.

An example of pore-space evolution that has received wide attention in the literature is the well-known reactive infiltration instability, which results in the formation in wormholes (Ortoleva et al. 1987; Hoefner and Fogler 1988; Steefel and Lasaga 1990). Similar to the unstable dendritic-growth process simulated by Xu and Meakin (2011), in the reactive infiltration instability, mineral dissolution takes places under transport-controlled conditions. In this instability, however, dissolution is limited by the rate of advection, which leads to the more rapid dissolution in high-flow-velocity pathways (wormholes) over slower pathways.

Direct numerical simulation at the Darcy scale has shown that the wormholes become more highly ramified and ultimately diffuse when the Damköhler number is small (surface-reaction-controlled) for a given system (Steefel and Lasaga 1990; Steefel and Maher 2009). More recently, Golfier et al. (2002) used the multiscale Stokes–Brinkman equations (11–12) to simulate the formation of wormholes. The reactive transport equations were formulated by upscaling the pore-scale equations to the Darcy scale with the volume-averaging method (cf. Appendix A Golfier et al. 2002). This upscaling approach is discussed in the next section. As a result, the surface reaction (i.e., dissolution) was represented at the Darcy scale in its upscaled form, via a mass-transfer coefficient (Golfier et al. 2002), thus avoiding the direct representation of the surface as in the level set or phase-field method examples presented above.

The model was able to reproduce the dissolution regimes observed experimentally with increasing flow rates (Fig. 11): (a) face dissolution (diffusion dominates over advection and instabilities do not develop), (b) conical wormholes (instabilities are present but still strongly influenced by diffusion), (c) dominant wormholes (advection-dominated instabilities), (d) ramified wormholes, and (e) uniform dissolution. Golfier et al. (2002) pointed out that for (a)–(c) very rapid dissolution led to local equilibrium conditions at the sharp interfaces, while for (d)–(e) faster flow rates led to local non-equilibrium dissolution, with more diffuse interfaces, and even stable displacement observed for (e). The Stokes–Brinkman model was able to capture flow in the multiscale domain, with Stokes flow being prominent in the free-flow regions (i.e., wormholes) and Darcy flow dominant in the partially dissolved porous medium.


While numerical simulation of pore-scale processes can help explain biogeochemical behavior not easily predicted at larger scales, it can also be used to derive the parameters that apply to models at larger scales. For this purpose, the pore-scale problem, in which interfaces are resolved explicitly and thus the variation of medium properties is not continuous, needs to be upscaled to the Darcy-scale continuum, in which they vary continuously (Fig. 4). Homogenization by volume averaging is one of a number of upscaling methods. Volume averaging has been successfully used as an upscaling method for porous-media processes (Whitaker 1999; Golfier et al. 2002; Wood et al. 2007). In particular, upscaling by volume averaging can be used to generate a mechanistic description of the effective parameters from the microscale representation of the physical and biogeochemical properties as well as the geometry of the pore space (Wood et al. 2007).

In the approach of Wood et al. (2007), the pore-scale equations (5–6) are first averaged over a representative control volume and then the pore-scale concentrations and velocities are decomposed into volume averages and deviations from these averages. As a result of this process, a macroscale mass balance equation is obtained that provides a continuum representation of the pore-scale processes. However, in this form the solution still depends on the concentration deviations, which are pore-scale quantities, through the hydrodynamic dispersion, macrodiffusion and the macroscale reaction rate terms. A closure problem for the concentration deviations needs to be posed to model their behavior. Golfier et al. (2002) and Wood et al. (2007) assume that simple periodic unit cells (e.g., solid cubes or spheres) capture the essential features of the system for the computation of the effective rate. Substitution of the solution of the closure problem makes it possible to write the macroscopic equations as a function of averaged concentrations, now using effective parameters, e.g., Wood et al. (2007):



where ciβ is the spatially averaged concentration of the reactant in the fluid portion of the pore space (β), V is an effective pore water velocity, defined as the intrinsic average of the fluid velocity vector in the averaging volume (uβ), D* is the effective hydrodynamic dispersion tensor, θ is the porosity, Av is the surface area per unit volume of porous medium, and Ri,eff is the macroscale effective reaction rate (specifically, the contribution of the rate of the microbially mediated reaction to the mass of reactant i).

Essentially, the effective parameters contain the pore-scale information that has been averaged in the upscaling process. While Wood et al. (2007) obtained a closure problem solution for the zero- and first-order limits of the Michaelis–Menten reaction rate, it was not possible for the authors to use the nonlinear kinetics in Equation (8) for the solution of the closure problem without coupling the closure to the value of the macroscale concentration, ciβ. This would have significantly increased the complexity of the problem. Rather, they proposed the following semi-empirical nonlinear form for the macroscale effective reaction rate term:



where Km,eff is now an effective half-saturation constant for the reaction.

Comparison of results of a direct numerical simulation of pore-scale reactive transport (taken as ground truth) with the macroscopic, upscaled version suggested that the variance of the concentration field had a dramatic impact on the effective reaction rate (Fig. 12). As a result, the accuracy of Equation (22) as a predictor of the effective reaction rate decreased as the variance increased. These discrepancies occurred mostly because of the errors induced by the proposed average of the nonlinear term; more specifically because the order of averaging operations cannot be interchanged for nonlinear processes such as Michaelis–Menten kinetics.


Applications of direct numerical simulation of pore-scale processes in subsurface materials are growing in part due to the computational advantages of well-established computational fluid dynamics (CFD) methods, e.g., Bijeljic et al. (2011a,b), Zaretskiy et al. (2012), Blunt et al. (2013). Availability of high-level, open-source libraries such as OpenFOAM (OpenFOAM 2015) or Chombo (Adams et al. 2014; Colella et al. 2014) have also facilitated implementation of models into new simulation capabilities. In reactive systems, the ability to calculate the heterogeneous reaction rates at the fluid–solid interfaces makes DNS a suitable tool to mechanistically model processes that are not captured by models at larger scales, especially when well-mixed assumptions are made. In this chapter, selected applications to evaluate transport limitations to reactive surfaces and to capture the evolution of both the surface area and the pore space have been reviewed.

Explicit representation of reactive interfaces as they evolve remains, however, a significant challenge for direct numerical simulation approaches. The level-set and the phase-field method have been successfully applied to dissolution–precipitation reactions in rather simple computational domains (Li et al. 2010; Xu and Meakin 2011). In these simulations, a single mineralogy is assumed. However, natural reactive systems are characterized by physical and mineralogical heterogeneity at a variety of different scales. As a result, reactive interfaces are often discontinuous and diffuse, e.g., Noiriel et al. (2009) and Deng et al. (2013). In multiscale problems, direct numerical simulation has typically been performed using continuum methods rather than strictly using pore-scale methods. In these methods, however, surface reactions (as well as transport processes) are not explicitly calculated at the interface and are upscaled from the pore scale. Upscaling of strongly coupled non-linear processes is not straightforward and the solution of macro- and micro-scale problems is coupled. Applying the same conceptual model for the reaction rate expression (e.g., Eqn. 8) at the pore scale and at the Darcy scale with averaged quantities is in general not warranted.

In this context, pore-scale modeling has the potential to challenge the conceptual models that are routinely applied at the Darcy scale, in particular for heterogeneous reactions. To accomplish this objective, however, direct numerical simulation needs to be able to incorporate the physical and mineralogical heterogeneity of the subsurface environments at different scales that are captured by advanced characterization techniques, e.g., Landrot et al. (2012). Further, methods developed for the evolution of pore spaces, such as those of Li et al. (2010) and Xu and Meakin (2011), need to be applied to complex subsurface porous media, as characterized by these advanced techniques. It is ultimately the combination of new microscopic characterization, experimental, and modeling approaches that presents the opportunity to provide mechanistic explanations for many of the long-standing questions in geochemistry in porous media.


The author would like to thank Li Li and Carl Steefel for reviews of the manuscript. This material is based upon work supported as part of the Center for Nanoscale Control of Geologic CO2, an Energy Frontier Research Center funded by the U.S. Department of Energy, Office of Science, Office of Basic Energy Sciences under Award Number DE-AC02-05CH11231.