Reactive Interfaces in Direct Numerical Simulation of Pore-Scale Processes

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 …


INTRODUCTION
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 defi ned as the largest spatial scale at which it is possible to distinguish the different fl uid 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).Specifi cally, 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 fl uid 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 infi ltration 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 fl uid-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 fl ow, 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.

Flow
Given the low compressibility of water, the incompressible Navier-Stokes equations can be used to accurately describe the fl ow 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 fl uid velocity,  the fl uid density, and μ the dynamic viscosity.Inertial forces have their origin in the convective acceleration of the fl uid as it fl ows 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 simplifi ed 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 fl uid phase and delineated by the fl uid-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 c i is the molal concentration of species i in solution M is the mass

Molins
fraction of water (kg H 2 O kg -1 ), i D is the diffusion coeffi cient (m 2 s -1 ) and i R is the reaction term (mol m 3 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 fl uid 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 signifi cant under acidic conditions (Molins et al. 2012;Ovaysi and Piri 2013).The reaction term R i 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 dissolutionprecipitation take place at solid-fl uid boundaries and are not included in Equation ( 5).

Surface reactions
Interfacial reactions can instead be expressed as a boundary condition at the solid-fl uid boundaries where m r is the surface reaction rate (expressed in units of mass per unit time per unit surface) and im  is the stoichiometric coeffi cient 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 k m is the intrinsic rate constant (mol m -2 reactive surface s -1 ), a E is the activation energy (kcal mol -1 ), n i a  is a product representing the inhibition or catalysis of the reaction by various ions in solution raised to the power n, with i a being the activity of species I, G is the Gibbs free energy (kcal), with 1 m , 2 m , and 3 m being three parameters that affect the affi nity 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 affi nity or G term can often be neglected, as in the studies of chlorite dissolution under high-pCO 2 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 fl uid-solid interface, assuming microbial cells were uniformly distributed on the surface (Fig. 1 where k m 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 K m is the half-saturation constant associated with the redox species i that defi nes 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 (s i mol m -2 surface) as a function of species concentrations: Reactive Interfaces in Direct Numerical Simulation of Pore-Scale Processes 465 max , 1 where s max and  represent the maximum surface concentration and the equilibrium constant, respectively (Zaretskiy et al. 2012).

DIRECT NUMERICAL SIMULATION
Direct numerical simulation involves the use of conventional discretization methods to solve the fl ow, transport, and geochemical equations.These include Eulerian, meshbased methods such as fi nite-difference, fi nite-element, or fi nite-volume methods, in which the differential equations are discretized by defi ning the value of the unknowns at fi xed points in space.Finite differences have been applied to the simulation of fl ow 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 fi nite-element method has been used for the simulation of fl ow and transport in sphere packs (Cardenas 2008(Cardenas , 2009) ) and as a mixed fi nite-element-fi nite-volume method for fl ow and reactive transport in Fontainebleau sandstone (Zaretskiy et al. 2012).A fi nite-volume method has been applied for fl ow and reactive transport in computer-generated and experimentally derived simulation domains (Molins et al. 2012(Molins et al. , 2014)).Direct numerical simulation has also been used in combination with other approaches.For example, Yoon et al. ( 2012) used fi nite-volume methods to simulate reactive transport with a fl ow 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 fi nite-difference solver was employed for Navier-Stokes fl ow.Both structured (e.g., Molins et al. 2014) andunstructured (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 classifi ed 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, fl uid-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 fl ow 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 fi nite-volume unstructured methods are not affected by this issue because they can capture arbitrary surfaces by appropriate meshing strategies (e.g., Cardenas 2008Cardenas , 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.

Molins
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 specifi cally for surface reactions while still solving for fl ow and transport within a binary, staircase domain.Another approach, the embedded boundary method, uses cut cells that result from intersecting the irregular fl uid-solid interfaces with the structured Cartesian grid (Fig. 3a) to discretize the equations using the fi nite-volume method (Trebotich et al. 2014;Trebotich and Graves 2015).The embedded boundary method is thus a more rigorous approach in that fl ow, 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 surfi cial 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 fl uid-solid interface,   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 suffi cient 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 Reactive Interfaces in Direct Numerical Simulation of Pore-Scale Processes 467 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 Cranfi eld 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 fl ow 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 fl ow 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 fl ow 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 fl ow in multiscale domains, such as fractured or vuggy media, can instead be simulated using the Stokes-Darcy or the Stokes-Brinkman equations (Golfi er 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-fl ow subdomains.To close the model, mass-and momentum-conservation equations are specifi ed 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 (Golfi er et al. 2002;Popov et al. 2009;Gulbransen et al. 2010) are modifi ed 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 fl uid 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 infi nity, 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 fl uid fl ow (Popov et al. 2009;Yang et al. 2014).A very similar approach has been recently used by Scheibe et al. (2015) to simulate fl ow 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 fl uxes 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 fl ow.

SURFACE AREA ACCESSIBILITY AND EVOLUTION IN MINERAL REACTIONS
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 fl ow 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.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-pCO 2 conditions.The laboratory experiment consisted of the injection of a solution at 4-bar pCO 2 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 Darcyscale and pore-scale effl uent 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-fl owing pore spaces that exchanged mass by diffusion with fast fl ow 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 signifi cant 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 infi ltrated by CO 2 -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 infi ltration of CO 2 -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 infi ltrating solution) precedes dissolution of the individual grains of sparite (Fig. 7).Noiriel et al. ( 2009) described mathematically the time-dependent RSA, A r , with the expression where 0 r A is the initial surface area, rm A 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 n 1 , n 2 , and n 3 are empirical coeffi cients 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 (CaSiO 3 ) 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 3 CaCO ( ),  with 0 a being the specifi c surface area in units of mol m -2 , ' p k is a proportionality factor that can be obtained assuming that when t→∞, A eff (∞) = 0, and 2 / 3 p  for monodisperse spherical particles, and 2 / 3 p  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).

PORE-SPACE EVOLUTION
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 fl uid-mineral interface, the resulting changes have impacts at larger scales.For example, (bio)precipitates can signifi cantly 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 CO 2 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-fl uid interface ( ) n u  can be described by, e.g., Li et al. (2010): where V m 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-fl uid (Hirt andNichols 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 fi xed Cartesian grids without having to parameterize these surfaces.Diffuse interface models, such as the phase-fi eld 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-fi eld method are described below in the context of recent applications to modeling pore-scale dissolution and precipitation processes.Reactive Interfaces in Direct Numerical Simulation of Pore-Scale Processes 473

Level set method
In the level-set method, the fl uid-solid interface   ( )  x is represented by a contour level of the level-set function   ,t  x , Equation (10).The evolution of the level set is governed by the following advection equation, which describes the motion of the fl uid-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 affi nity term to account for nearequilibrium conditions.The three-dimensional effects of fl ow conditions and reaction rates were explored quantitatively.The simulation showed that the evolution of the permeabilityporosity relationship depended on particular parameter values used (i.e., fl ow 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-fi eld method
Phase-fi eld methods are based on the idea that the free energy depend on an order parameter (the phase-fi eld 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 fi eld.The fi eld changes smoothly in a diffuse zone around the interface that has a fi nite width.
A sharp-interface asymptotic analysis of the phase-fi eld 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 fi eld is a function of the evolution of the phase fi eld: where  is proportional to the molar volume of the mineral.The evolution of the phase fi eld is described with where is the phase-fi eld characteristic time parameter,  is closely related to the interface thickness, and  controls the strength of the coupling between the phase fi eld  and the concentration c.The relationship between   and  is the n obtained from an asymptotic analysis, which ensures that the phase-fi eld equations converge to the sharp interface solution: Two-dimensional dendritic growth due to solute precipitation was simulated using this phase-fi eld model (Xu and Meakin 2011).The simulations were performed under diffusionlimited conditions by setting the chemical reaction rate much larger than the rate of diffusioni.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 snowfl ake-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-fi eld 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-fl uid 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 infi ltration 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 infi ltration 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-fl ow-velocity pathways (wormholes) over slower pathways.
Direct numerical simulation at the Darcy scale has shown that the wormholes become more highly ramifi ed and ultimately diffuse when the Damköhler number is small (surfacereaction-controlled) for a given system (Steefel and Lasaga 1990;Steefel and Maher 2009).More recently, Golfi er et al. (2002) used the multiscale  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 Golfi er 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 coeffi cient (Golfi er et al. 2002), thus avoiding the direct representation of the surface as in the level set or phase-fi eld method examples presented above.
The model was able to reproduce the dissolution regimes observed experimentally with increasing fl ow rates (Fig. 11

UPSCALING OF SURFACE REACTIONS BY VOLUME AVERAGING
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;Golfi er 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 fi rst 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.Golfi er et al. ( 2002) and  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 i c  is the spatially averaged concentration of the reactant in the fl uid portion of the pore space ( ),  V is an effective pore water velocity, defi ned as the intrinsic average of the fl uid velocity vector in the averaging volume ( ),  u * D is the effective hydrodynamic dispersion tensor,  is the porosity, v A is the surface area per unit volume of porous medium, and ,eff i R is the macroscale effective reaction rate (specifi cally, 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 fi rst-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, c i  .This would have signifi cantly increased the complexity of the problem.Rather, they proposed the following semi-empirical nonlinear form for the macroscale effective reaction rate term: where K m,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 fi eld 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 specifi cally because the order of averaging operations cannot be interchanged for nonlinear processes such as Michaelis-Menten kinetics.

SUMMARY AND OUTLOOK
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 fl uid 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 fl uid-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 signifi cant challenge for direct numerical simulation approaches.The level-set and the phase-fi eld

F igure 3 .
Example of an irregular geometry on a Cartesian grid (left), in which shaded areas represent volume of cells excluded from domain: (a) Embedded boundary representation with interfaces "cutting" regular cells (left) and single-cut cell showing boundary fl uxes (right); (b) Binary "staircase" representation of the interface [Adapted from Trebotich et al 2014].

Figure 5 .
Figure 5. Schematic representation of the transport phenomena in an idealized fracture

Figure 6 .
Figure 6.Results from reactive transport simulations of a crushed calcite capillary experiment (Molins et al. 2014).pH values are shown on calcite grain surfaces while vectors (black arrows) indicate direction and magnitude of fl uid velocities in the pore space.Values of pH along fast fl ow paths are lower than in slow-fl ow paths, resulting in heterogeneous dissolution rates on minerals surfaces located in close proximity.The diameter of the capillary tube is 524 μm.[Modifi ed from Molins S, Trebotich D, Yang L, Ajo-Franklin JB, Ligocki TJ, Shen C, Steefel CI (2014) Pore-scale controls on calcite dissolution rates from fl ow-through laboratory and numerical experiments.Environmental Science and Technology, Vol.48, p. 7453-7460 Abstract Art with permission.© 2014 American Chemical Society.]

Figure 7 .
Figure 7. (top) Dissolution process according to the sugar-lump model.The matrix is composed of (a) spherical grains of surface area 0 r A (b) which dissociate into smaller grains, thus increasing the waterexposed surface area.(c) Subsequently, the individual particles dissolve, which reduces the surface area.(bottom) Thin section of the rock observed with optical microscopy in transmitted light.(b) Scanning electron microscopy (SEM) observations of the sparite and the micrite.[Reprinted from Noiriel C, Madé B, Gouze P (2007) Impact of coating development on the hydraulic and transport properties in argillaceous limestone fracture.Water Resources Research, Vol.43, W09406, Fig 1, 2 with permissions.]

Figure 8 .
Figure 8.(a) Backscattered SEM images of a cross-section of a carbonated wollastonite grain after 2 days of reaction in circumneutral pH conditions.Note the succession of the inner intact core of reacting wollastonite, a fractured layer composed with calcite and silica, and the compact continuous and poorly permeable coating of calcite.(b) Normalized extents of wollastonite carbonation as a function of time in circumneutral pH conditions.The lower two curves correspond to fi ts of the data assuming the progressive formation of a passivating coating of calcites.[Reprinted from Daval D, Martinez I, Corvisier J, Findling N, Goffé B, Guyot F (2009) Carbonation of Ca-bearing silicates, the case of wollastonite: Experimental investigations and kinetic modeling.Chemical Geology, Vol.265, Figs.3c, 6, with permission.]

Figure 9 .
Figure 9. Evolution of the three-dimensional solid surface together with the velocity and concentration fi elds for precipitation in a pore throat with an initially sinusoidal-shaped aperture under different conditions expressed in the form of the dimensionless Damköhler (Da) and Péclet (Pe) numbers.[Reprinted from Li X, Huang H, Meakin P (2010) A three-dimensional level set simulation of coupled reactive transport and precipitation/dissolution. International Journal of Heat and Mass Transfer, Vol.53, p. 2908-2923, Fig, 6 with permission.] ): (a) face dissolution (diffusion dominates over advection and instabilities do not develop), (b) conical wormholes (instabilities are present but still strongly infl uenced by diffusion), (c) dominant wormholes (advection-dominated instabilities), (d) ramifi ed wormholes, and (e) uniform dissolution.Golfi er 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 fl ow 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 fl ow in the multiscale domain, with Stokes fl ow being prominent in the free-fl ow regions (i.e., wormholes) and Darcy fl ow dominant in the partially dissolved porous medium.

Figure 10 .
Figure 10.Snapshots of interfacial patterns of unstable dendritic growth obtained with the phase fi eld method of Xu and Meakin (2008) for precipitation under diffusion-limited conditions, at various simulation times (dimensionless) : a) t = 0.002 b) t = 0.004 c) t = 0.006 d) t = 0.008.[Reprinted from Xu Z and Meakin P (2011) Phase-fi eld modeling of two-dimensional solute precipitation/dissolution: Solid fi ngers and diffusion-limited precipitation.Journal of Chemical Physics, Vol.134, 044137, Fig 3. with permission.]