Modeling Reactive Transport Processes in Fractures

Fractures are ubiquitous and important features in the Earth subsurface (Berkowitz 2002; Pyrak-Nolte et al. 2015). They are created as a result of rock failure when the critical stress (i.e., fracture toughness) is exceeded, or in the case of subcritical crack growth, when cracks propagate under stress conditions below fracture toughness, facilitated by chemical reactions (Atkinson 1984). The necessary conditions for fracture growth can be created by natural unloading from land erosion (Engelder 1987), tectonic events (Molnar et al. 2007), and crystal growth in presence of fluids supersaturated with respect to solid phases (Royne and Jamtveit 2015). Fractures can also be artificially created for enhanced energy recovery, through excess fluid pressure and change of thermal stress, as in the case of geothermal energy extraction and unconventional oil and gas production (McClure and Horne 2014; Lampe and Stolz 2015). Fractures can be observed by surveying rock outcrops, inferred by fluid flow and geochemical measurements, and detected using geophysical techniques (Berkowitz 2002; St Clair et al. 2015; Walton et al. 2015). When open, fractures act as preferential flow pathways because of their high permeability, and thus typically control fluid migration and solute transport in fractured rocks. For this reason, fractures are avoided when siting and designing geologic isolation systems, such as for nuclear waste and CO₂ storage, in order to prevent undesired fluid and chemical migration (Kovscek 2002; Lewicki et al. 2007; Birkholzer et al. 2012). In the Earth's critical zone, fractures control the availability of water for rock weathering and hence the development of the regolith layer (Brantley et al. 2017). It has also become accepted that weathering itself typically controls fracture permeability in hard rock aquifers (Lachassagne et al. 2011).

Simulations of flow and transport in fractured porous media are of prime importance in studies related to oil and geothermal exploration, waste disposal, and environmental stewardship of contaminated sites. With that aim, three broad categories of conceptual and numerical models have been developed, to capture the unique flow and transport properties of fractures: equivalent continuum models, dual or multiple continua models, and discrete fracture networks. The equivalent continuum model treats the fractures and rock matrix as a single continuum with a permeability corrected for the effect of fractures (Long et al. 1982). In the dual-continua (dual-porosity and dual-permeability) models (Warren and Root 1963), a separate fracture continuum is added, co-located with the rock matrix continuum, with its own hydraulic properties including porosity-permeability relationships that can differ from those applied to the matrix. The conceptualization of dual continuum models can vary, depending on the types of connections implemented between fracture and matrix grid blocks, with advantages and disadvantages, as discussed by Lichtner (2000). Multiple-continua models build upon the dual-continua concept by adding more co-located interacting media, thus allowing consideration of more hydrogeological components (Doughty 1999). For example, multiple interacting continua (MINC) models can be used to capture transient processes resulting from mass transfer between fractures and rock matrices at different distances from the fracture surface (Pruess 1991). Similarly, a triple-continuum approach has been proposed to investigate the effect of small fractures on flow through fractured bedrock (Wu et al. 2004).

In contrast to the dual- and multiple-continua models, discrete fracture networks represent the individual fractures explicitly, although with reduced dimensionality, i.e., (n − 1)-dimensional fractures in an n-dimensional domain (Hyman et al. 2015; Lei et al. 2017). Because the discrete fracture network is generated based on the geometrical properties of each individual fracture (e.g., aperture, length, orientation), and can account for internal aperture variability (Makedonska et al. 2016), the flow and transport results are arguably more accurate than with the dual/multiple-continua approach. However, the application of discrete fracture networks is typically not suitable for large-scale systems with numerous fractures, or when information on fracture geometrical properties is unavailable.

The morphology and thus hydrophysical properties of fractures evolve dynamically as a function of mineral dissolution and precipitation driven by fluid composition, temperature, and pressure variations. This is the case when advective flow in fractures introduces fluids that are out of chemical equilibrium with the wallrock, or displaces fluids to zones of different temperature and/or pressure. For example, solubility changes by cooling, and boiling through depressurization are key processes leading to fracture alteration in hydrothermal systems (e.g., Browne 1978). Mixing of working fluids in geothermal systems with shallower and cooler groundwater can also lead to significant mineral deposition (Griffiths et al. 2016). In the case of geologic carbon storage systems, the introduction of CO₂ into deep aquifers creates an acidic and carbonated fluid that reacts with minerals (Rochelle et al. 2004), especially carbonates, with the potential of impacting fracture permeability and caprock integrity (Fitts and Peters 2013). Acidic fluids are also introduced in hydraulic fracturing processes to dissolve minerals and thus create and maintain open fractures (Kalfayan 2008). In all these cases, a detailed understanding of reactive transport mechanisms leading to mineral precipitation and dissolution in fractures is crucial to the development of predictive tools that can be applied with confidence in studies such as engineered system design and optimization, energy recovery feasibility studies, and environmental impact assessments.

Much work has been done, for decades, towards the understanding and modeling of fluid flow in fractured systems (National Research 1996). Many studies have also been performed since the 1960's towards understanding the geochemical processes responsible for the mineralogical alteration of fractures (veins) observed in hydrothermal systems and ore deposits (e.g., Hemley and Jones 1964; Meyer 1967; Browne 1978; Barnes 1997; Parry 1998), including the concept of water/rock ratio to characterize mass transfer between fractures and rock matrix (Giggenbach 1984). Comparatively much fewer investigations have focused on understanding and modeling the actual mechanisms of coupled chemical reaction and transport within fractures and between fractures and rock matrix. Models of rock matrix–fracture interactions, including multicomponent geochemical effects, have been developed using both discrete fracture networks and dual/multi-continua approaches in the context of hydrothermal systems, geologic nuclear waste disposal and geologic carbon sequestration, both at the field scale (e.g., Steefel and Lichtner 1998; Sonnenthal et al. 2005; Xu et al. 2006; Gherardi et al. 2007; Liu et al. 2017) and smaller scales (e.g., Steefel and Lasaga 1994; Xu and Pruess 2001; Dobson et al. 2003). However these studies relied on simplified conceptualizations of fracture properties (Snow 1970; Witherspoon et al. 1980), reactive surface areas, and porosity-permeability relationships (Sonnenthal et al. 2005). Advances in computational power and microscopic imaging have since allowed the development of more sophisticated and mechanistic models that can more accurately capture the intricate interplay between fluid flow and reaction in fractures, including effects of local fluid dynamics, matrix diffusion, fluid chemistry variability and mixing, and the impacts of mineral precipitation/dissolution on reactive surface areas, as further discussed here.

This chapter focuses on the geochemical alteration of single fractures, which is critical for improving the conceptualization and prediction of flow and transport in fractured rocks at the reservoir scale. Here, we concentrate our discussions on single-phase saturated systems, whereas a detailed discussion of multiphase processes can be found in Sin and Corvisier (2019, this volume). First, experimental observations of fracture alteration caused by fluid–rock interactions are summarized. This is followed by a review of the modeling efforts dedicated to simulating reactive transport processes in individual fractures. We also discuss the major factors controlling fracture alteration as elucidated by the modeling efforts, and conclude the chapter with an outlook on directions where further research is needed.

The detailed observation of fracture alteration driven by geochemical reactions is made possible by the advancement of non-invasive imaging techniques to characterize fracture surfaces and geometry, such as profilometry (Ameli et al. 2013), Nuclear Magnetic Resonance Imaging (Dijk et al. 1999), Positron Emission Tomography (Loggia et al. 2004; Tenchine and Gouze 2005), and Computed Tomography (Gouze et al. 2003). In particular, the combination of fracture-flow experimental apparatus with computed tomography imaging has enabled in situ visualization of fracture morphology changes due to reactive fluid flow (Deng et al. 2015; Ajo-Franklin 2017). Improvement of imaging detectors and resolutions, and development of new image processing algorithms, have also further improved our capability to quantify fracture geometries; this can be challenging because fracture aperture is orders of magnitude smaller than the dimensions of the fracture plane. For example, an inverse point-spread function was developed to accurately determine fracture apertures that are one-tenth of the voxel size (Ketcham et al. 2010), and filters that detect the unique characteristics of fracture morphology, i.e., the planar geometry, were used to distinguish fracture from large pores/vugs in the rocks (Deng et al. 2016a). Segmentation algorithms have also been developed and applied to differentiate the fracture void, the rock matrix immediately adjacent to the fracture that has gone through geochemical alteration, and the intact rock matrix (Noiriel et al. 2007; Deng et al. 2013). Furthermore, coupling of tomography imaging with microscopic techniques such as Scanning Electron Microscope (SEM), Energy Dispersive X-ray Spectroscopy (EDS) and X-Ray Fluorescence (XRF), and the application of techniques such as dual-energy microtomography provide spatially resolved geochemical data in addition to bulk information derived from fluid chemistry (Noiriel 2015). By correlating the grayscale values of the 3D tomography images with high resolution 2D SEM/EDS images, algorithms can be trained to extract mineralogical information from tomography images, and to map mineral composition on fracture surfaces (Ellis and Peters 2016). Such information is critical for investigating the controls of mineralogy on the changes in fracture flow properties.

The continuous replenishment of fresh reactive fluid can accelerate mineral reactions in fractures. When the fluid is under-saturated with respect to the minerals in contact, minerals dissolve, resulting in fracture aperture enlargement (Dijk et al. 2002; Detwiler et al. 2003). Fracture opening caused by the dissolution of carbonate minerals when exposed to acidic fluid is of particular interest because of its prevalence in the subsurface and its relevance to various natural processes (e.g., karst formation) and human activities (e.g., geologic carbon storage). Given the fast kinetics of carbonate minerals (Chou et al. 1989), especially calcite, considerable fracture opening within days or hours has been observed in experimental studies that use aggressive synthetic acidic fluids (Singurindy and Berkowitz 2004; Ellis et al. 2011; Elkhoury et al. 2013; Noiriel et al. 2013; Emmanuel and Levenson 2014; Deng et al. 2015; Garcia-Rios et al. 2015; Ajo-Franklin 2017). Figure 1a shows 3D reconstructions of a limestone fracture at different time intervals of a flow-through experiment, demonstrating a uniform fracture opening over time caused by calcite dissolution (Noiriel et al. 2013). The average aperture increased from 48 µm to 346 µm after reacting with a CO₂-charged fluid (P_CO2= 0.1 MPa) for 55 hours. Figure 1b shows a series of fracture aperture maps (2D projections on the fracture plane of the 3D fracture volumes reconstructed from in situ tomography images) and histograms that highlight localized fracture enlargement (Deng et al. 2015). After reacting the limestone fracture with CO₂-rich brine (P_CO2 = 77 bar) for 42 hours, the average aperture increased from 244 µm to 344 µm, while local apertures exceeded 1000 µm.

The intrinsic fracture permeability increases in response to fracture enlargement caused by mineral dissolution. For a parallel wall fracture, the volumetric flux flowing through the fracture (Q, [m³/s]) for a given pressure gradient (∇P, [Pa/m]) is proportional to the cube of the fracture aperture (b, [m]) (i.e., the cubic law, Eqn. 1), and thus fracture permeability (k_f, [m²]) is proportional to the square of fracture aperture (Eqn. 2) (Snow 1969):

where W [m] is the width of the fracture, and μ [Pa s] is the dynamic viscosity of the fluid.

For non-parallel fractures, a hydraulic aperture (b_h, [m]) can be defined such that the cubic law still holds. The increase of the hydraulic aperture does not necessarily follow the geometric aperture. Fracture permeability increases more substantially than would be expected based on the cubic law if preferential flow channels develop (Deng et al. 2015; Garcia-Rios et al. 2017), while the increase in fracture permeability is mitigated when fracture enlargement is accompanied by an increase in fracture roughness and hydraulic tortuosity (Ellis et al. 2011; Deng et al. 2013; Noiriel et al. 2013). For example, unreacted mineral bands can make fracture surfaces rougher and serve as flow strictures (Fig. 2).

Fracture aperture enlargement caused by mineral dissolution may not be a linear function of time because of decreasing subsequent mineral dissolution. One mechanism that causes mineral dissolution to slow down is the development of preferential channels (Deng et al. 2015). As more reactive fluid is directed into these channels, the surface area in the nonchannelized regions becomes less accessible for reactions. Furthermore, as preferential channels grow, transverse diffusion (i.e., perpendicular to the fracture plane) in these channels becomes important and the reactions are increasingly limited by the transport process. Another mechanism that limits subsequent mineral dissolution is the retreat of the mineral front into the rock matrix (Noiriel et al. 2007; Abdoulghafour et al. 2013; Davila et al. 2016a; Ajo-Franklin 2017). As such, reactants and products of the dissolution reactions have to diffuse through the remaining rock matrix, creating a transport limitation on chemical reactions.

When the rock matrix is composed of minerals of varying reactivity, preferential dissolution of the fast reacting minerals can result in the development of an altered porous layer on the fracture surface (Noiriel et al. 2007; Ellis et al. 2011; Elkhoury et al. 2015; Davila et al. 2016a; Ajo-Franklin 2017). For example, Figure 3 shows an altered layer composed of quartz and other silicates left behind after calcite dissolution (Davila et al. 2016a). Within the altered layer, removal of the fast reacting minerals also exposes more surface area of the remaining minerals (Garcia-Rios et al. 2015).

The impacts of the alteration in the near-fracture region on fracture permeability is complicated. The permeability of the altered layer is negligible compared to that of the open fracture (Noiriel et al. 2007; Chen et al. 2014). If the altered layer remains near the fracture surface, the initial flow path is mostly maintained and the fracture permeability change is limited (Davila et al. 2016a). Mineral grains in the altered layer, however, may reorganize and lead to local aperture decreases (Noiriel et al. 2007). The altered layer may also detach from the fracture surface, causing fracture aperture to increase locally (Noiriel et al. 2007; Andreani et al. 2008). Detachment of the altered layer may result from the decrease in cohesion relative to the shear stress imposed by the fluid flow (Noiriel et al. 2007), or due to the change of fluid pH at the grain boundaries (Pepe et al. 2010). The fracture permeability can increase following the detachment of the altered layer if the fracture flow flushes out the released particles (Deng et al. 2017); whereas the fracture permeability is decreased if the released particles re-deposit in the fractures and cause substantial fracture aperture reduction (Fig. 3d) (Noiriel et al. 2007; Ellis et al. 2013; Davila et al. 2016a).

The alteration of the near-fracture region is also common in cements. Typically, three altered layers are observed (Fig. 4), including a Ca-depleted layer close to the unaltered cement, an amorphous Si-rich layer bordering the bulk fluid, and a calcite-enriched layer in between (Abdoulghafour et al. 2013; Luquot et al. 2013; Walsh et al. 2014a; Abdoulghafour et al. 2016). These altered layers are products of the reaction sequence of portlandite dissolution, the precipitation of amorphous silica, and calcite precipitation due to local enrichment of Ca²⁺ at the portlandite dissolution front, respectively. More details about acid–cement interactions is discussed in Cama et al. (2019, this volume.) While the net result of the reactions is the removal of materials, the degree of fracture permeability change is widely variable. This change depends on the relative thickness of the three layers, which is partly controlled by flow rate, and the dissolution pattern in the fracture, which is affected by the initial fracture geometry. Fracture permeability increases when channelization (Abdoulghafour et al. 2016) or net precipitation of a low porosity calcite layer (Luquot et al. 2013) dominates, and remains unchanged when the boundary of the secondary Si-rich layer tracks the initial geometry. In some cases, fracture permeability decreases. Because the amorphous silica layer has a large porosity, the thickness of this layer can be considerable even with a small amount of amorphous silica, leading to a decrease in fracture aperture and permeability (Abdoulghafour et al. 2013). The decrease in fracture permeability is also partly attributed to the weakened mechanical properties (Walsh et al. 2013, 2014a), i.e., deformation of the altered layer, or failure of the asperities under confining pressure.

Fracture permeability can also decrease as a result of mineral precipitation when the injected fluids are supersaturated with respect to solid phases. The supersaturation state can be achieved by changing the temperature of the fluid, or by mixing two fluids of different compositions (Emmanuel and Berkowitz 2005). By injecting a well-mixed CaCl₂–NaHCO₃ solution supersaturated with respect to calcite into a transparent fracture analogue, Jones and Detwiler (2016) observed calcite precipitation in the artificial fracture, and showed that the precipitation patterns are sensitive to the initial mineral heterogeneity, i.e., seeding of calcite on the initial fracture surface.

Mineral precipitation is also possible when the necessary aqueous species are provided by simultaneous mineral dissolution, such as in the cases of calcite dissolution with gypsum precipitation (Garcia-Rios et al. 2015; Davila et al. 2016a) (Fig. 3b,c), and dolomite dissolution with calcite precipitation (Singurindy and Berkowitz 2004). Precipitation of carbonate minerals in basaltic fractures has also been observed, because the dissolution of minerals such as pyroxene and olivine releases a considerable amount of divalent cations in solution (Menefee et al. 2017).

In these cases, the direction and extent of fracture permeability change depends on the interplay between dissolution and precipitation. While fast dissolution supplies more reactants for the precipitation reaction, it also causes substantial fracture opening that can outcompete the impact of mineral precipitation and result in fracture permeability increase (Singurindy and Berkowitz 2005; Garcia-Rios et al. 2015; Davila et al. 2016a). It was also shown that transport limitation, such as experienced by dead ends of fractures and less connected fractures (Singurindy and Berkowitz 2005; Menefee et al. 2017), or due to slow flow velocity (Davila et al. 2016a), favors precipitation and fracture closing. Figure 5 summarizes experimental observations in a smooth fracture going through coupled calcite dissolution and gypsum precipitation, highlighting the important roles of fluid chemistry and flow velocity (Singurindy and Berkowitz 2005). In addition to the factors discussed above, the molar volume and the porosity of the precipitates can also impact fracture permeability change caused by precipitation. For the same amount of precipitates, a larger porosity in the precipitate corresponds to a sharper reduction in the fracture aperture and permeability, as in the case of the secondary Si-rich layer in cement fractures (Abdoulghafour et al. 2013).

As highlighted by experimental observations, the key processes involved in geochemical fracture alteration include: advective fluid flow, diffusion along the fracture plane, diffusion across the fracture walls and into/from the altered layer adjacent to the fracture, and accompanying mineral reactions (Fig. 6). In this section, we summarize reactive transport models of different complexities developed to examine one or more of these processes in single fractures. While theoretically, all these models can be applied to individual fractures of any length scale (typically between ~100 µm to ~100 m), in practice, the complexity of the model needs to be reduced as the length scale of the individual fracture being investigated increases. Another common practice is typically to apply pore-scale models to cases when the processes across fracture aperture (with length scales on the order of µm–mm) are of interest, whereas continuum-scale models are used when processes along the flow direction or in the fracture plane are the focus of the investigation.

Mesh-based direct numerical simulations of fracture alteration involves numerically solving a system of equations. The mathematical problem involves solving three governing equations (for continuity, momentum conservation, and mass conservation) and several constitutive relations (e.g., reaction rate laws, mass laws). The modeling approach can be categorized as either pore-scale or continuum-scale, depending on whether the fluid–rock interface (e.g., fracture surface) is explicitly resolved, or whether the solid phase and fluid phase co-exist in each model (numerical) grid cell, respectively. Here, we provide a summary of the key features of the pore-scale and continuum-scale modeling approach. More details can be found in Molins (2015) and Steefel et al. (2015), respectively.

With this type of modeling, the void space itself is discretized and the fluid–solid interface is tracked (Molins et al. 2012, 2014). Therefore, geometry evolution due to mineral precipitation and dissolution is explicitly modeled. The generic formulation of the governing equations of pore-scale models are given as follows:

where u is the velocity vector [m s⁻¹], ρ [kg m⁻³] and μ [Pa s] are the density and dynamic viscosity of the fluid, respectively, P [Pa] is pressure, and F [N m⁻³] is the body force such as gravity. C_i [mol m⁻³] and D_mi [m² s⁻¹] are respectively the concentration and the molecular diffusion coefficient of aqueous species i. The mass flux of species i from reaction of mineral m (r_m [mol m⁻² s⁻¹]) that is in contact with the fluid is treated as a boundary condition:

where υ_i,m is the stoichiometric coefficient of species i in the reaction of mineral m, and n is the unit normal vector.

The continuum-scale governing equations are based on variables averaged over each grid cell where the void space and the solid phase co-exist, with its void space modeled using porosity. Assuming steady state flow, the continuity equation (Eqn. 7) and Darcy's law (Eqn. 8) are solved for the flow field:

where v [m s⁻¹] is the pore velocity, q [m³ m⁻²_{porous medium} s⁻¹] is the Darcy flux, k [m²] is the permeability, and g [m s⁻²] is the gravitational acceleration. In the mass balance equation (Eqn. 9), ϕ [m³_{pore space} m⁻³_{porous medium}] is the porosity, D is the dispersion tensor, which can be used to account for dispersion in addition to molecular diffusion in porous media. Concentration changes due to mineral reactions are treated as source/sink terms (R_i, [mol m⁻³_{porous medium} s⁻¹]).

A common practice is to write the mass balance equation (Eqn. 9) with respect to the total concentration of component i (ψ_i), which is the summation of the concentration of the corresponding primary species p_i, and the concentrations of all secondary species (s_i) that are formulated using primary species p_i:

The concentrations of the secondary species (C_si) is related to the concentrations of primary species by the law of mass action (Eqn. 11), which assumes instantaneous equilibrium for aqueous reactions.

where K_eq,si and υ_{si, pi} are the equilibrium constant and the stoichiometric coefficient of the aqueous reaction that relates the secondary species s_i to the primary species p_i, respectively, and γ is the activity coefficient.

Mineral reactions are typically treated kinetically. The reaction rate of a given mineral, r_m [mol m⁻² s⁻¹], can be calculated as follows using the transition state theory rate laws (Lasaga 1984):

where k_rxn,m [mol m⁻² s⁻¹] is the kinetic coefficient (pre-exponential factor), E_a [kcal mol⁻¹] is the activation energy, R [kcal K⁻¹ mol⁻¹] is the ideal gas constant, T [K] is the absolute temperature, a_i is the activity of any catalytic or inhibitory species, and n_i is the corresponding exponent. The last term is the chemical affinity that depends on the ionic activity product (IAP_m), the solubility constant (K_eq,m), and the empirical exponents m₁ and m₂. Note that theoretically, 1/m₁ is the Temkin coefficient, which corresponds to the stoichiometric coefficient of the elementary reaction step (Lichtner 2016).

Unlike the pore-scale models, in which Equation (12) is used as a boundary condition (Eqn. 6) at the surface of each mineral, in the continuum-scale models, the reaction term (R_i, Eqn. 9) is treated as a source/sink term, accounting for the contribution from all minerals that are present in the grid cell and include species i. It relates to Equation (12) as:

where υ_i,m is the stoichiometric coefficient of species i in the reaction of mineral m, and A_m [m² m⁻³_{porous medium}] is the reactive surface area, which is related to specific surface areas [m²/g_mineral or m²/m³_mineral] through porosity and mineral molar volume.

In the continuum-scale models, porosity is updated based on the amount of mineral precipitation and dissolution taking place within each simulation time step:

where V_m [m³/mol] is the molar volume of the mineral, and superscripts t and t + Dt represent time before and after a time interval Δt, respectively. Permeability is updated from porosity based on user-specified constitutive relationships.

In continuum models, surface areas are updated based on simplified conceptualizations of geometric evolution, and typically follow a power law relationship with mineral volume fraction and/or porosity (Noiriel et al. 2009). In contrast, in pore-scale models, the fluid–solid interface is tracked as the reaction progresses and the reactive surface area is updated based on the evolution of the actual geometry. In addition, because the hydrodynamic processes are resolved at the pore scale, any transport limitation is accounted for explicitly instead of being lumped into an effective reactive surface area, as is common practice in continuum-scale models (Noiriel and Daval 2017).

One-dimensional continuum models have been typically used to investigate the alteration of rock matrix immediately adjacent to a fracture. Because this type of model focuses on reaction front propagation into an essentially impermeable rock matrix, a continuum approach is adopted. The only transport mechanism is diffusion, and fluid chemistry in the fracture defines the Dirichlet (fixed) boundary condition. The governing equation is therefore simplified to

where D_eff is the effective diffusion coefficient in the porous media, and is related to porosity (f) by some correlation, such as Archie's Law:

where a and c, which is typically known as the cementation exponent, are empirical coefficients depending on the textures of the rock.

This model has been used to investigate the alteration of rock matrix bordering a hyper-alkaline fluid-filled fracture under conditions relevant to nuclear waste repositories (Steefel and Lichtner 1994). The study predicted calcite precipitation in the rock matrix immediately adjacent to the fracture; the precipitated calcite was shown to serve as an armoring layer preventing the rock matrix from neutralizing the hyper-alkaline fluid in the fracture and from sorbing and retarding the transport of radionuclides in the long run. This type of 1D model has also been shown to qualitatively capture the three altered layers created by portlandite dissolution, calcite precipitation and calcite re-dissolution, and amorphous silica precipitation (Luquot et al. 2013).

An alternative approach for tracking the reaction front into the rock matrix neighboring a fracture treats the reaction front migration as a moving boundary problem (Ulm et al. 2003; Walsh et al. 2014a). In this approach, the transport between fronts is still governed by diffusion (Eqn. 17), and the front movement is described by mass conservation of the controlling species at the front (Eqn. 18). The concentrations of the aqueous species at each reaction front are controlled by local chemical equilibrium with mineral phases.

In these equations ⟦ ⟧ calculates the difference across the front, ρ_i [kg m⁻³] is the density of species i in the solid, Z_f is the location of the front, and M_i [kg mol⁻¹] is the molar mass of species i.

If the fluid chemistry boundary condition at the fracture surface is fixed, the 1D model implicitly assumes that the local fracture fluid chemistry does not vary with the mineral reactions and is independent of flow and transport in the fracture. This assumption is valid if the fluid velocity in the fracture is sufficiently fast or if the region of interest is close to the inlet/recharge zone (Steefel and Lichtner 1994). Alternatively, the 1D model can be coupled with modules that simulate the flow and transport within the fracture separately, and update the boundary conditions accordingly (Walsh et al. 2013, 2014a).

In order to capture the changes of fluid chemistry along the flow direction, as a result of upstream reactions in addition to rock matrix alteration in the near-fracture region, a 2D domain that is perpendicular to the fracture plane needs to be simulated. Equation (9) is applicable to the entire domain if the fracture is treated as a continuum. However, properties of the fracture grid cells need to be defined in different ways: the porosity is calculated from the local fracture aperture (Eqn. 19), and the permeability is defined and updated following the cubic-law formulation (Eqn. 20),

where Δz denotes the size of the fracture grid cells, other variables as defined previously, and the subscript 0 is used for the initial value (at time t = 0).

This type of 2D model has been used to successfully reproduce the fluid chemistry evolution observed in a series of flow-through experiments reacting fractured marl cores with carbonated water, as well as the thickness of the altered layer on the fracture wall due to calcite dissolution (Davila et al. 2016b). This type of 2D model would be the equivalent of a 1D dual-/multi-continua model (Dobson et al. 2003), in which the computational domain was discretized along the fracture and the interface area was configured to account for the mass transfer between the fracture and rock matrix continua.

In the continuum approach discussed above, because fracture aperture is not resolved explicitly, the underlying assumption is that fluid within the fracture is well-mixed and the concentrations of all aqueous species are homogeneous. This assumption does not necessarily hold because the fluid velocity is not uniform, and the velocity gradient can influence solute transport and create concentration gradients between the fracture walls. For example, in parallel fractures, fluid velocity follows a parabolic profile and approaches zero at the fracture surface (Poiseuille flow). Mineral reaction products tend to accumulate within this hydrodynamic boundary layer, where solute transport is dominated by diffusion.

To resolve the velocity gradient and the resulting concentration gradient, a pore-scale approach is needed. For parallel fractures, the analytical solution of Equation (4), i.e., the parabolic velocity profile, can be directly substituted into Equation (5):

where u₀ is the average velocity. The reaction is given by the boundary condition, i.e., the diffusive flux at the fracture surface (at z = ±b / 2) is equal to the mass flux due to reactions. For simple first-order surface reactions and equilibrium reactions, analytical solutions can be derived for the average concentration across the fracture aperture (Berkowitz and Zhou 1996). These analyses provide insights regarding the breakthrough of reactive solutes in the fractures and the effective dispersion coefficients that can be used for upscaling.

Numerical solutions of this model can be used to simulate systems with non-linear reaction rates. For example, Li et al. (2008) performed numerical simulations for several minerals under different kinetic dissolution constraints in parallel fractures, and investigated the flow and reaction conditions that would lead to the breakdown of the well-mixed assumption, i.e., the continuum treatment (Fig. 7a). The authors found two necessary conditions for the development of concentration gradient across the fracture aperture: (i) comparable advection and reaction rates, and (ii) reaction rates that are faster than diffusion across the aperture. The model can also be applied in rough fractures. Instead of assuming a parabolic velocity profile, the full Navier–Stokes equation is solved to capture the hydrodynamics that arise from surface roughness (Fig. 7b). As such, the impacts of surface roughness on local reaction at the fracture wall and overall reaction in the fracture can be investigated (Deng et al. 2018b). Fracture roughness has been observed to increase following mineral dissolution (Gouze et al. 2003; Ellis et al. 2011), therefore considering the compounded impacts and feedback of surface roughness on mineral reactions is of great importance.

Multiscale models can be used to simultaneously capture pore-scale fluid flow in the fracture and diffusion-controlled reactive transport in the rock matrix (Molins et al. 2019). In this approach, the computational domain is divided into the fracture domain, in which the pore-scale governing equations are solved, and the porous media domain where the continuum-scale equations are used. At the boundary between the two domains, diffusive fluxes are balanced. The study of (Wen et al. 2016) represents a simplified implementation of the hybrid approach, in which fracture aperture is discretized and assigned velocity values according to a modified parabolic profile that corrects for roughness. Using this model, the authors investigated the evolution of the diffusivity of the rock matrix bordering a rough fracture as a result of preferential dissolution of calcite in the rock matrix, and its dependence on calcite abundance. More details of the multi-scale modeling approach can be found in Molins and Knaber (2019, this volume).

Pore-scale models can also be applied in both the fracture and rock matrix domains, which has the advantage of explicitly accounting for hydrodynamics and surface area change in the altered layer. This is, however, computationally expensive for finite-volume based reactive transport models. Alternative approaches, such as the Lattice Boltzmann method, have been explored (Chen et al. 2014; Fazeli et al. 2018). Chen et al. (2014) investigated fracture evolution in a binary mineral system using a pore-scale model, and showed the development of an altered layer, in which the flow velocity is negligible. This confirms that the dominant transport mechanism in the altered near-fracture region is diffusion, as has been assumed in the 1D continuum-scale modeling discussed above.

Flow variations in the fracture plane (i.e., along its surface) can arise from perturbations such as initial fracture aperture, and influence the diffusive transport of solutes in fractures (Nowamooz et al. 2013). Such flow instabilities and their positive feedback with geochemical reactions, as demonstrated by experimental observations, can result in self-organization phenomena, e.g., fracture channelization (Ortoleva et al. 1987; Deng et al. 2015). Models that discretize the two dimensions in the fracture plane have been used to explicitly capture these flow and transport processes (Hanna and Rajaram 1998; Detwiler and Rajaram 2007; Szymczak and Ladd 2012; Deng and Peters 2018). In this type of 2D model, the third dimension is considered by integrating over the fracture aperture (Eqn. 22):

where C_i is the depth-averaged concentration, R_i is the effective reaction rate [mol m⁻² s⁻¹], and q [m² s⁻¹] is the local flux, solved using the Reynolds equation, i.e., depth-averaged Stokes equation:

For single mineral systems, assuming the reactive surface area is the geometric fracture surface area (i.e., 2xy), the aperture is updated using

It is noted that in reality the reactive surface area is affected by the topography of the fracture surfaces. Typically, it is corrected by multiplying the geometric fracture surface area by a surface roughness factor. The surface roughness factor can be determined independently if the surface geometry is mapped or BET surface area is measured. However, BET measurements may result in over-estimation because they capture, and are typically dominated by, pore space surface area that is not in direct contact with the fluid flow in the fracture. The surface roughness factor can also be inferred from the extent of reactions (Deng et al. 2016b). Overall, in the continuum approach, for which the geometry is not explicitly resolved, the determination of reactive surface area has been challenging and requires guidance based on knowledge of the system or pore-scale investigations as illustrated by Deng et al. (2018b).

Dividing Equation (22) by the thickness of the fracture domain considered, which is essentially the grid size in the z direction (z, [m]) in the 2D model and can be an arbitrary value larger than the local aperture, we get the mass balance equation of the continuum model (Eqn. 9) where the porosity follows the definition in Equation (19). The same pore velocity as calculated from Equation (23) can also be recovered from Darcy's law (Eqn. 8) by using a permeability that is the intrinsic fracture permeability (defined in Eqn. 2) weighted by the porosity (defined in Eqn. 19).

Although the 2D fracture plane model does not explicitly resolve the non-uniform velocity profile across the fracture aperture, i.e., the hydrodynamic boundary layer close to the fracture wall, modifications can be made to account for the resulting impacts on mineral reactions. Because the major transport mechanism at the fracture surface is diffusion, for fast reacting minerals, the reaction rates are limited by solute diffusion to and from the fracture surface. The diffusive flux (i.e., diffusion controlled reaction rate R_diff, [mol m⁻² s⁻¹]) can be expressed as

where Sh is the Sherwood number, which is a function of reaction rate and the velocity profile across the aperture and is bounded by two asymptotic limits 7.54 and 8.24 (Gupta and Balakotaiah 2001), and C_is is the concentration at the fracture surface. Assuming the concentration at the fracture surface is at equilibrium with respect to the mineral, R_diff can be evaluated independently from known variables. Hanna and Rajaram (1998) used the smaller of R_diff and the kinetically controlled reaction rate to account for the diffusion limitation caused by a strong hydrodynamic boundary layer.

Alternatively, an effective reaction rate coefficient (k^’_{rxn, eff} [m s⁻¹]) can be derived analytically for first order reactions with a stoichiometric coefficient of one, given that the diffusive flux has to be balanced by the mass flux from kinetic reaction (R_rxn):

where C_ieq is the equilibrium concentration, and k^’_rxn is the kinetic rate coefficient for the first order reactions written in unit of [m s⁻¹], which can be estimated from k_rxn used in Equation (12). This formulation allows a smooth transition between the kinetically controlled and transport controlled reaction rate (Detwiler and Rajaram 2007; Szymczak and Ladd 2012).

The so-called 2.5D model builds upon the 2D fracture plane model and accounts for the alteration of the near fracture region and the diffusion process through the altered layer (Deng et al. 2016b, 2017, 2018a). In the 2.5D model, each grid cell is a continuum where the fracture and rock matrix on the fracture wall co-exist with a total thickness of Δz. This model solves the same governing equations as the 2D fracture plane model, but with two important modifications.

The fracture aperture is defined for each mineral (m) and measures the distance between the reaction fronts of this mineral in the two fracture halves. It is calculated from f_i,j,m, the volume fraction of the mineral measured for the intact rock matrix and V_i,j,m, the volume fraction of the mineral in the grid (i, j):

Therefore, the reaction front for each mineral is tracked and the distance between the reaction front and the fracture surface (constrained by the mineral that retreats the slowest), d_AL, can be calculated.

The diffusion controlled reaction rate due to the presence of the altered layer (R_ALdiff [mol m⁻² s⁻¹]) is calculated as:

An effective reaction rate (R_eff) is then defined as the harmonic mean of the kinetically (R_rxn) and diffusion (R_ALdiff) controlled reaction rates, such that the lower one of the two dictates the reaction. This is similar to how the diffusion limitation caused by the hydrodynamic boundary layer is accounted for (Eqn. 28):

The 2.5D model successfully reproduced channeling in the fracture plane, the development of the altered layer, and the resulting diffusion limitation and decreasing overall dissolution of calcite observed in a dolomite fracture reacting with carbonated water (Deng et al. 2016b). The conceptualization embodied in the 2.5D model can also be implemented using a 2D dual-contniua model, in which the mass transfer between the fracture continuum and the rock matrix continuum is explicitly defined by Equations (30) and (31).

By capturing the development of the altered layer, the 2.5D model also enables considering the erosion of the altered layer, which is proportional to the thickness of the altered layer (d_AL) as observed in previous experiments (Andreani et al. 2008). A phenomenological law (Eqn. 32) has been implemented in the 2.5D model to successfully simulate fracture opening due to detachment of the altered layer observed in a carbonate-rich shale (Deng et al. 2017):

where E [m s⁻¹] is the erosion rate, and d_ALc is the critical thickness, which has to be exceeded for the erosion to take place. Both η [m^1−e s⁻¹] and e are empirical coefficients that are expected to vary with mineral composition of the rock matrix and can be well constrained by experimental observations.

Full 3D pore-scale modeling of the geochemical alteration of fractures involves numerically solving the pore-scale governing equations (Eqns. 3–5) in the actual fracture geometry. It captures the detailed fluid flow, diffusion and reaction dynamics, including flow instabilities in the fracture plane and the velocity variations across the fracture aperture. However, this modeling approach is computationally intensive because of the fine mesh resolution needed to resolve the fracture aperture and the relatively large domain as determined by the size of the fracture plane.

There are only a few studies performed in simple fracture geometries (Starchenko et al. 2016; Starchenko and Ladd 2018). Their model was implemented in OpenFOAM with customized libraries for boundary conditions and mesh updates. It was shown that the 3D model results deviate noticeably from the 2D fracture plane model results after channels are established, even with the implementation of the diffusion limited reaction rate in the 2D model. This was primarily attributed to the breakdown of Reynolds equation, which is used in the 2D model for fluid flow simulations. It has been shown that the use of the Reynolds equation may introduce large errors when the flow velocity is high, i.e., Reynolds number is large, or when the fracture has a large roughness (Brown et al. 1995; Zimmerman and Bodvarsson 1996; Oron and Berkowitz 1998).

There are even fewer 3D simulations on the evolution of fracture alteration with more complex, yet realistic geometries (e.g., Szymczak and Ladd 2009, using the Lattice Boltzmann method), let alone with mineralogical heterogeneity in the near fracture region.

Numerical models that capture the important reactive transport processes in fractures enable numerical experiments for systematic exploration of a broad parameter space. Flow rate, as it directly relates to residence time (Glassley et al. 2002), is the most widely studied controlling factor. The role of mineralogy, including mineral composition and spatial distribution, has received increasing attention. The following discussion focuses on dissolution-driven fracture alteration, based on the studies that are available.

Discussions of the impacts of flow rate on fracture dissolution have been commonly framed using the dimensionless Peclet (Pe) and Damköhler (Da) numbers, which are defined for fractures as follows:

where u is the average velocity in fracture, L is the length of the fracture, and k_rxn^’ is the reaction rate coefficient in unit of m s^−1,Da_I captures the relative magnitude of reaction rate and advection rate, whereas Da_II measures the ratio between reaction rate and diffusive mass transfer rate.

These dimensionless numbers were originally derived from dimensional analysis of the governing equations (Berkowitz and Zhou 1996; Detwiler and Rajaram 2007), and involve simplifications of geometry and reaction kinetics. Therefore, it should be noted that the applicability of these dimensionless numbers is limited. As the impact of mineralogy will be discussed in the following subsection, we focus on a single mineral system for the discussion of flow rate, for which the traditional definitions of Pe and Da still provide a useful framework.

Consistent observations have been made on the impact of flow rate on fracture dissolution patterns (Detwiler and Rajaram 2007; Szymczak and Ladd 2009, 2011; Elkhoury et al. 2013; Deng et al. 2018a; Starchenko and Ladd 2018). Some of the data points from both numerical simulations and experimental observations are summarized in Figure 8 (Starchenko and Ladd 2018), where Da_eff is a variant of Da_I (Eqn. 34a). The choice of Da_I versus Da_II is somewhat arbitrary in constructing this type of figures (a similar alternative figure can be readily derived from Pe and Da_II), and the discussion of the dependence of dissolution patterns on either Da_I or Da_II can be convoluted when transport transitions from being advection-dominated to diffusion-dominated. At large Pe, when advection dominates, fracture dissolution patterns largely depend on Da_I, i.e., the relative magnitude of the reaction and advection rates. When the advection rate is large compared to the reaction rate (low Da_I), dissolution is uniform because the fast replenishment of reactive fluid homogenizes the concentration field and mineral reactions in the fracture. At large Da_I, the advection rate is slower than the reaction rate, and the reactants are consumed faster than they are transported. The resulting dissolution pattern is referred to as face dissolution (in flow-through experiments), as the dissolution is highly localized to the inlet, and does not result in much change in fracture permeability. At intermediate Da_I when the flow rate and the advection rate are comparable, the feedback between flow and reaction is most pronounced. Initial perturbations are magnified to create preferential channels in the fractures, leading to a dramatic fracture permeability increase. The morphology of the channels are also shown to be dependent on Da_I. The channels are more compact at slightly larger Da_I, and are more diffuse at smaller Da_I as fracture dissolution transitions into a uniform dissolution regime. At small Pe and small Da_II, diffusion is the dominant transport process and tends to smooth out reaction instabilities and, therefore, results in a relatively uniform dissolution pattern. If the reaction rate is much faster than the diffusion rate (large Da_II), the reaction will be localized at the inlet, leading to face dissolution. Overall, this is consistent with the fact that fewer data points fall into the wormholing regime at low Pe on Figure 8.

Both mineral composition and the spatial distribution of various, especially reactive, minerals influence fracture alteration. In rocks where reactive minerals are well mixed with non-reactive or slow-reacting mineral phases, fracture alteration typically involves the development and erosion of the altered layer depending on the volume percentage of fast-reacting minerals. Deng et al. (2018a) performed numerical investigations of two multi-mineral scenarios. Based on these simulation results, conceptual models (Fig. 9) and multi-mineral Damköhler number (mDa) were developed for analyses of fracture alteration patterns in these systems. If fast-reacting minerals only account for a small volume fraction of the rock matrix, the altered layer in the near fracture region subsists or slowly dissolves. In this case, the fracture dissolution pattern and its dependence on the flow rate are controlled by the slow reacting minerals (Eqn. 35). In contrast, if the fast-reacting minerals account for a large volume fraction of the rock matrix, detachment of the porous altered layer is expected. The fracture geometry change is then controlled by the thickness of the altered layer, i.e., the fast reacting minerals (Eqn. 36).

where V_frac is the volume of the fracture, Q is the volumetric flux, b is the average aperture, and f_r is the volume fraction of a given mineral.

The multi-mineral Da highlights the observation that the flow regimes corresponding to different dissolution patterns shift depending on the mineral composition.

If instead, the reactive minerals follow certain spatial patterns, the fracture alteration patterns are largely controlled by the spatial distribution of the reactive mineral phase (Deng and Peters 2018; Noiriel and Deng 2018). When the reactive minerals are discontinuous along the flow direction, fracture aperture enlargement is localized to where reactive minerals are present, but the change of hydraulic properties of the fracture is constrained by the non-reactive minerals. When the reactive minerals are continuous along the flow direction, reaction fronts tend to penetrate deeper into the fracture compared to the case when no mineral heterogeneity is present. The dissolution of the reactive minerals can thus result in considerable increase in fracture hydraulic properties (Fig. 10). Faster hydraulic opening is expected if the spatial distribution of the reactive minerals overlap with the preferential flow paths prescribed by the initial fracture geometry.

Overall, advancements in experimental techniques have provided valuable insights into important reactive transport processes in fractures. The experimental observations have also provided indispensable data sets for the development and testing of conceptual and numerical models, which have enabled systematic investigation of different controlling factors, and the development of more mechanistic approaches to improve our predictive capability. Enhanced understanding of the role of mineralogy on fracture alteration resulting from recent studies is of particular interest, as it is an important and necessary step in bridging fundamental understanding and studies of natural and engineered systems, where multiple minerals are typically present.

In this chapter, we summarized major reactive transport processes observed in single fractures. The morphology and hydraulic properties of fractures evolve as a result of the interplay between fluid flow in the fracture, diffusion in the fracture plane and across the fracture aperture, mineral dissolution/precipitation, and diffusion through altered zones adjacent to the fracture. Fracture permeability typically increases when the alteration is driven by mineral dissolution. The extent of the increase depends on the changes in fracture morphology, which is largely controlled by the relative magnitude of the time scale of transport and the time scale of reactions. While the increase of fracture permeability is limited when the dissolution is localized at the inlet, i.e., large Da, fracture permeability can increase sharply when the reaction rate and the transport rate are comparable and initial perturbations are amplified to produce self-organization features such as channels in the fractures. The altered layers bordering the fracture surface, which is created by the preferential dissolution of fast-reacting minerals in multi-mineral systems, can mitigate fracture permeability increase by increasing fracture surface roughness. The altered layers can also passivate subsequent mineral dissolution due to diffusion limitation. However, depending on the mechanical strength of the altered layers, these may detach from the fracture surfaces and result in fracture enlargement, while releasing fine particles that may eventually clog the fractures. Fracture permeability also decreases when precipitation dominates in the fractures.

In addition, we highlighted existing process-based models for the prediction of the alteration of an individual fracture driven by geochemical reactions in fully saturated rocks. This was not meant to be an exhaustive review of the literature but rather to provide some context for model selection. Model complexity, and, to a large extent, model accuracy, is inversely related to the computational requirement, and it is advised that model selection be based on the system and research question being investigated. Understanding the major reactive transport processes that each type of model addresses provides a basis for such consideration. Both the continuum-and pore-scale modeling approaches have been applied to simulate fracture alteration. 1D models along the flow direction or perpendicular to the fracture plane provide useful tools for investigating solute transport along a single direction through advection and diffusion, and their controls on mineral reactions. 2D models perpendicular to the fracture plane allow detailed examination of processes across the fracture aperture or into the bordering rock matrix, while considering flow and transport along the flow direction. When aperture and flow variations within the fracture plane cannot be ignored, 2D models within the fracture planes are necessary to capture the fracture morphological change and fracture permeability evolution from mineral reactions. The 2.5D models enables consideration of processes across the fracture aperture, such as diffusion limitation of the altered layers, while keeping the computational efficiency of a 2D model. This is especially useful in mineralogically complicated systems. 3D pore-scale models simulate explicitly all the physical and chemical processes that affect fracture alteration and thus are expected to produce more accurate results regarding fracture permeability change; but their computational costs are typically high. Therefore, their application in geochemically complex, yet realistic, systems where interactions between the fracture and the bordering rock matrix are important is unpractical. Hybrid models provide a promising alternative that balances the computational costs and the details of the processes involved.

By outlining the current landscape of this field, we also hope to help identify research gaps that remain. A few examples are discussed below.

The controlling factors of precipitation-driven fracture alteration are only partly understood. Precipitation is complicated by nature. It involves nucleation and crystal growth, both of which are highly sensitive to local heterogeneity and are dependent on a range of factors (Noiriel et al. 2016). The shape, size and spatial distribution of the neo-formed crystals are dependent on the material substrate (Noiriel et al. 2016), and the type of precipitate is affected by flow rate and local saturation state. For instance, aragonite was observed to precipitate in a fracture at faster flow rates, whereas calcite precipitated at slower flow rates (Singurindy and Berkowitz 2004). Therefore, while models described in this chapter may be applied directly to systems where mineral precipitation dominates, caution needs to be exercised. For example, it is expected that in some cases, the concentration gradient across the fracture aperture results in super-saturation at the fracture surface whereas the bulk fluid chemistry may indicate under-saturation.

Mineral reactions in fractures are typically accompanied by a redistribution of local stress and a change of the rock mechanical properties (Elkhoury et al. 2013, 2015). Therefore, a modeling capability that resolves coupled mechanical and geochemical processes in fractures is needed. Some groundbreaking efforts have been taken in this direction, and highlighted that geomechnical processes can mitigate fracture opening caused by dissolution reactions (Yasuhara et al. 2004; Ameli et al. 2014; Walsh et al. 2014b; Kim et al. 2015; Spokas et al. 2018). A mechanistic understanding of mechanical weakening of near fracture regions (such as the detachment of the altered layer) and modeling of microscopic processes are still lacking.

Multiphase flow dynamics in factures are highly dependent on fracture geometry (Yang et al. 2016). This geometry dependence on the fracture surface mineralogy and the resulting contact between the reactive fluid phase and minerals are also important factors that affect fracture evolution in multiphase systems and need to be further investigated.

Most importantly, findings derived from experimental and modeling studies at an individual fracture scale need to be integrated into reservoir-scale studies. Implementation of conceptualizations developed from studies discussed in this chapter into the dual-/multi-continua models discussed earlier provides a potential upscaling framework for reservoir scale applications.

The authors would like to thank the support by Laboratory Directed Research and Development (LDRD) funding from Berkeley Lab, provided by the Director, Office of Science, of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231. We also would like to thank reviewers and the editors for their constructive comments.