# Effect of Pore-Scale Mineral Spatial Heterogeneity on Chemically Induced Alterations of Fractured Rock: A Lattice Boltzmann Study

Authors: Hossein Fazeli; Ravi Patel; Helge Hellevang
Journal: Geofluids
 (2018)

---

## Abstract

Fractures are the main flow path in rocks with very low permeability, and their hydrodynamic properties might change due to interaction with the pore fluid or injected fluid. Existence of minerals with different reactivities and along with their spatial distribution can affect the fracture geometry evolution and correspondingly its physical and hydrodynamic properties such as porosity and permeability. In this work, evolution of a fracture with two different initial spatial mineral heterogeneities is studied using a pore-scale reactive transport lattice Boltzmann method- (LBM-) based model. The previously developed LBM transport solver coupled with IPHREEQC in open-source Yantra has been extended for simulating advective-diffusive reactive transport. Results show that in case of initially mixed structures for mineral assemblage, a degraded zone will form after dissolution of fast-dissolving minerals which creates a resistance to flow in this region. This causes the permeability-porosity relationship to deviate from a power-law behavior. Results show that permeability will reach a steady-state condition which also depends on transport and reaction conditions. In case of initially banded structures, a comb-tooth zone will form and the same behavior as above is observed; however, in this case, permeability is usually less than that of mixed structures.


---

## Body

## 1. Introduction

During reactive transport processes in tight rocks, fractures play an important role as it is the main flow path for species transport. Existence of fractured seals in geological CO2 sequestration and fractures in the host rock of nuclear waste disposal sites are some practical examples in which hydrodynamic properties of fractures can help to better understand the long-term evolution of the system [1]. These hydrodynamic properties will change due to chemical disequilibrium resulting from interactions between pore fluid and rock minerals. There are experiments which have shown that the presence of reactive and nonreactive (or less reactive) minerals and their spatial distribution with respect to each other can affect the way a fracture evolves as a result of reactions with reactive fluids [2]. For instance, when slow- and fast-reacting minerals are mixed in the rock matrix around a fracture, dissolution of fast-reacting minerals will form a degraded zone around the fracture [2, 3]. In another example, when fractures crossing zones of slow- and fast-reacting minerals form banded structures, the dissolution of fast-reacting minerals will generate a comb-tooth structure around the fracture [2, 3]. Therefore, fracture permeability evolution might not follow the same behavior when minerals form initially banded and mixed structures, and in addition to the experiments, numerical studies are needed to address the effect of spatial mineral distributions on fracture evolution. There have been different continuum-scale numerical studies that have investigated the fracture evolution [4–7]. In these models, mineral spatial heterogeneities are generally represented in a single discretized cell, but pore-scale models resolving spatial heterogeneities can provide better description of processes and help to understand the role of spatial heterogeneities on fracture evolution [8, 9]. Different pore-scale reactive transport modeling approaches, involving geometry changes due to chemical interactions, have been developed [10–17]. Lattice Boltzmann method- (LBM-) based models have also been developed which are mostly focusing on single mineral dissolution in nonfractured media [5, 18–33]. Recently, Chen et al. [26] investigated the effect of mineral heterogeneity on evolution of a single fracture surrounded by a rock matrix composed of binary minerals, but the model is limited to a simple mineralogy. In the present work, we utilize a similar setup used in Chen et al.’s work [26], a single fracture surrounded by a rock matrix composed of two different minerals, and focus more on the effect of mineral spatial heterogeneity on fracture geometry evolution by constructing initially mixed and banded mineral structures that are conceptualizations of the structures observed by Ellis et al. [2]. The pore-scale simulations are carried out using a multicomponent reactive transport LBM model proposed by Patel et al. [30] which has been implemented in open-source code Yantra [34]. This model couples the LBM transport solver with a well-known geochemical solver IPHREEQC which allows accounting for complex and realistic reaction networks. In this study, we have further extended the model to include the advective flow component and the kinetic reactions. We use mixed and banded structures in two different parts. In the first part of the paper, we perform single-species simulations to investigate fracture geometry evolutions and also study the effect of different flow and reaction conditions by choosing different dimensionless numbers. The second part of the paper presents results related to a realistic case where a single fracture (representing a leakage pathway) in a caprock, above a formation used for CO2 storage, is considered where the rock matrix is assumed to be composed of carbonate (calcite) and clay (kaolinite) minerals. This rock is representative of a carbonate-rich caprock such as some intervals in the Draupne shale [35].


## 2. Numerical Model Description

### 2.1. Lattice Boltzmann Method for Fluid Flow

The flow of incompressible fluids can be described by continuity and Navier-Stokes (NS) equations [36]
(1)∇·u=0,ρ∂u∂t+ρu·∇u=−∇P+μ∇2u+F,where ρ (kg/m3) is density, u (m/s) is velocity, P (Pa) is pressure, μ (Pa·s) is dynamic viscosity, and F (N/m3) is body force.LBM has been used in this study to solve the Navier-Stokes equation. LBM describes the behavior of a collection of particles, and instead of macroscopic equation of fluid dynamics, it is based on the Boltzmann equation which describes dynamics of a gas on a mesoscopic scale [37]:
(2)∂f∂t+ξα∂f∂xα+Fαρ∂f∂ξα=Ωf,where f is density of particles with velocity ξ in the α-direction at position x and time t. Also, F is the force term and Ω is a collision operator which represents local distribution of f because of collisions of particles. It can be shown that (2) can recover the NS equation meaning that if one solves (2) and obtains f, then one is able to calculate the macroscopic quantities such as velocity and pressure using f. To solve for f, (2) must be discretized in time, velocity space, and physical space [37]:
(3)fix+ciΔt,t+Δt=fix,t+Ωif+FiΔt,where
(4)Ωif=−Δtτfi−fieq,fieq=wiρ1+u·cics2+u·ci22cs4−u·u2cs2,Fi=wiciαcs2Fα,cs2=13·ΔxΔt2.In the above equations,Δx and Δt are space and time resolutions. fi is the discrete velocity distribution function, fieq is the equilibrium distribution function, and τ is the relaxation time. Discrete velocities ci and their corresponding set of weighing coefficients wi will form velocity sets. There are different velocity sets used in LBM, and in this study and for fluid flow LBM, we use D2Q9 where 2 is the number of spatial dimensions (x and y) and 9 is the set’s number of velocities [37]:
(5)wi=49,i=0,19,i=1,2,3,4,136,i=5,6,7,8,ci=0,0,i=0,cosπi−12,sinπi−12,i=1,2,3,4,2cosπi−52+π4,sinπi−52+π4,i=5,6,7,8.Equation (3) consists of two parts named collision and streaming. Particles move with velocity ci to a neighboring node at position x+ciΔt during time step t+Δt, and at the same time, they are affected by the collision operator. So, (3) can be broken into two parts [37]:
(6)Collisionplussourceterm:fi∗x,t=fix,t+Ωif+FiΔt,Streaming:fix+ciΔt,t+Δt=fi∗x,t.Once these two steps are performed, distribution functionsfi in the current timestep can be calculated and then one is able to compute the macroscopic fluid density and velocity [37]:
(7)ρ=∑i=08fi,u=1ρ∑i=08fici+FΔt2ρ.As was mentioned, it can be shown that the LB equation can recover the NS equation using Chapman-Enskog analysis [38], and it can be observed that the kinematic viscosity ν takes the following form in terms of Lattice parameters:
(8)ν=cs2τ−Δt2.


### 2.2. Lattice Boltzmann Method for Advection-Diffusion-Reaction Equation

Reactive transport of chemical species can be described using the advection-diffusion-reaction (ADR) equation [39]:
(9)∂Cj∂t+∇·−Dj∇Cj+uCj=Rj,where Cj is the concentration of species j, Dj is the diffusion coefficient which is considered constant for all species, and Rj represents the source or sink term due to reactions for species j. Rj is related to the homogeneous reactions while heterogeneous reactions happening at the solid-fluid interface are usually described as boundary conditions.LBM can also be applied to solve the ADR equation. So, the LB equation for each speciesj [30]
(10)gijx+ciΔt,t+Δt=gijx,t+Ωijg+wiRtotaljx,tΔt,with
(11)Ωijg=−Δtτggij−gieq,j,(12)gieq,j=wiCj1+u·cics2,solves the ADR equation for concentration of species j
(13)Cj=∑i=04gij.In the above equations,g is the distribution function used for species concentration. It should be pointed out that in the LBM model above, the D2Q5 lattice was used for the ADR equation. Rtotalj is a source-sink term due to both heterogeneous and homogeneous reactions. While in most of previous LBM models [2, 10, 18, 20, 21, 23–27, 29, 40–42], heterogeneous reactions are treated as a boundary flux, in this study, LBM treats heterogeneous reactions as pseudohomogeneous reactions meaning that Rtotalj at the bulk only includes the homogeneous reactions, but at the fluid nodes adjacent to a solid node, Rtotalj also includes heterogeneous reactions (now treated as pseudohomogeneous reactions) [30]. Using this approach, both homogeneous and heterogeneous reactions can be treated similarly, and it gives the possibility to couple the LBM transport part with an external geochemical solver such as IPHREEQC which has been used in this work. The validation of the LBM model has been presented in the Appendix.


### 2.3. Numerical Implementation

Figure1 shows how the LBM models for fluid flow and mass transport and the reaction solver have been coupled. To initialize the model, first, initial, and boundary concentrations of different species are computed using IPHREEQC. In each time step, first we solve for the velocity profile using the LBM flow solver (as described in Section 2.1). In this study, since the time scale of the fluid flow is much shorter than the other processes, velocity reaches steady-state condition during the chosen time step. Equilibrium distribution functions are calculated using (12), and the collision step is then executed according to the following equation:
(14)gij∗x,t=gijx,t+Ωijg.Figure 1
Flowchart of the LBM reactive transport model.We do not include the reaction source-sink term in this step. In the next step, molar concentrationsCtrnpj are transferred to IPHREEQC. IPHREEQC is a class of the PHREEQC geochemical solver that facilitates the usage of PHREEQC along with other software such as transport solvers [43]. After importing concentrations Ctrnpj from the transport step into the geochemical solver, new concentrations Cphrqj are calculated by the PHREEQC. Finally, Rtotalj is calculated using
(15)Rtotaljx,t=Cphrqj−Ctrnpj∆t.All these steps are done using a PHREEQC wrapper file which is part of the LBM code used in this study and has been detailed in Patel’s work [34] and in open-source code Yantra 1.0.0 which is available online (https://bitbucket.org/yantralbm/yantra/overview).The geometry is updated using static update rules detailed in Patel et al. [30]. Before performing the streaming step, to apply the effect of chemical reactions, we update the distribution functions using
(16)gij∗∗x,t=gij∗x,t+wiRtotaljx,tΔt,and then the streaming step is executed:
(17)gijx+ciΔt,t+Δt=gij∗∗x,t.Once the streaming step is performed, species concentrations can be obtained using (13), and we go to the next time step. It should be pointed out that in the next time step, we solve for the velocity when there is at least 1% change in the porosity.


## 2.1. Lattice Boltzmann Method for Fluid Flow

The flow of incompressible fluids can be described by continuity and Navier-Stokes (NS) equations [36]
(1)∇·u=0,ρ∂u∂t+ρu·∇u=−∇P+μ∇2u+F,where ρ (kg/m3) is density, u (m/s) is velocity, P (Pa) is pressure, μ (Pa·s) is dynamic viscosity, and F (N/m3) is body force.LBM has been used in this study to solve the Navier-Stokes equation. LBM describes the behavior of a collection of particles, and instead of macroscopic equation of fluid dynamics, it is based on the Boltzmann equation which describes dynamics of a gas on a mesoscopic scale [37]:
(2)∂f∂t+ξα∂f∂xα+Fαρ∂f∂ξα=Ωf,where f is density of particles with velocity ξ in the α-direction at position x and time t. Also, F is the force term and Ω is a collision operator which represents local distribution of f because of collisions of particles. It can be shown that (2) can recover the NS equation meaning that if one solves (2) and obtains f, then one is able to calculate the macroscopic quantities such as velocity and pressure using f. To solve for f, (2) must be discretized in time, velocity space, and physical space [37]:
(3)fix+ciΔt,t+Δt=fix,t+Ωif+FiΔt,where
(4)Ωif=−Δtτfi−fieq,fieq=wiρ1+u·cics2+u·ci22cs4−u·u2cs2,Fi=wiciαcs2Fα,cs2=13·ΔxΔt2.In the above equations,Δx and Δt are space and time resolutions. fi is the discrete velocity distribution function, fieq is the equilibrium distribution function, and τ is the relaxation time. Discrete velocities ci and their corresponding set of weighing coefficients wi will form velocity sets. There are different velocity sets used in LBM, and in this study and for fluid flow LBM, we use D2Q9 where 2 is the number of spatial dimensions (x and y) and 9 is the set’s number of velocities [37]:
(5)wi=49,i=0,19,i=1,2,3,4,136,i=5,6,7,8,ci=0,0,i=0,cosπi−12,sinπi−12,i=1,2,3,4,2cosπi−52+π4,sinπi−52+π4,i=5,6,7,8.Equation (3) consists of two parts named collision and streaming. Particles move with velocity ci to a neighboring node at position x+ciΔt during time step t+Δt, and at the same time, they are affected by the collision operator. So, (3) can be broken into two parts [37]:
(6)Collisionplussourceterm:fi∗x,t=fix,t+Ωif+FiΔt,Streaming:fix+ciΔt,t+Δt=fi∗x,t.Once these two steps are performed, distribution functionsfi in the current timestep can be calculated and then one is able to compute the macroscopic fluid density and velocity [37]:
(7)ρ=∑i=08fi,u=1ρ∑i=08fici+FΔt2ρ.As was mentioned, it can be shown that the LB equation can recover the NS equation using Chapman-Enskog analysis [38], and it can be observed that the kinematic viscosity ν takes the following form in terms of Lattice parameters:
(8)ν=cs2τ−Δt2.


## 2.2. Lattice Boltzmann Method for Advection-Diffusion-Reaction Equation

Reactive transport of chemical species can be described using the advection-diffusion-reaction (ADR) equation [39]:
(9)∂Cj∂t+∇·−Dj∇Cj+uCj=Rj,where Cj is the concentration of species j, Dj is the diffusion coefficient which is considered constant for all species, and Rj represents the source or sink term due to reactions for species j. Rj is related to the homogeneous reactions while heterogeneous reactions happening at the solid-fluid interface are usually described as boundary conditions.LBM can also be applied to solve the ADR equation. So, the LB equation for each speciesj [30]
(10)gijx+ciΔt,t+Δt=gijx,t+Ωijg+wiRtotaljx,tΔt,with
(11)Ωijg=−Δtτggij−gieq,j,(12)gieq,j=wiCj1+u·cics2,solves the ADR equation for concentration of species j
(13)Cj=∑i=04gij.In the above equations,g is the distribution function used for species concentration. It should be pointed out that in the LBM model above, the D2Q5 lattice was used for the ADR equation. Rtotalj is a source-sink term due to both heterogeneous and homogeneous reactions. While in most of previous LBM models [2, 10, 18, 20, 21, 23–27, 29, 40–42], heterogeneous reactions are treated as a boundary flux, in this study, LBM treats heterogeneous reactions as pseudohomogeneous reactions meaning that Rtotalj at the bulk only includes the homogeneous reactions, but at the fluid nodes adjacent to a solid node, Rtotalj also includes heterogeneous reactions (now treated as pseudohomogeneous reactions) [30]. Using this approach, both homogeneous and heterogeneous reactions can be treated similarly, and it gives the possibility to couple the LBM transport part with an external geochemical solver such as IPHREEQC which has been used in this work. The validation of the LBM model has been presented in the Appendix.


## 2.3. Numerical Implementation

Figure1 shows how the LBM models for fluid flow and mass transport and the reaction solver have been coupled. To initialize the model, first, initial, and boundary concentrations of different species are computed using IPHREEQC. In each time step, first we solve for the velocity profile using the LBM flow solver (as described in Section 2.1). In this study, since the time scale of the fluid flow is much shorter than the other processes, velocity reaches steady-state condition during the chosen time step. Equilibrium distribution functions are calculated using (12), and the collision step is then executed according to the following equation:
(14)gij∗x,t=gijx,t+Ωijg.Figure 1
Flowchart of the LBM reactive transport model.We do not include the reaction source-sink term in this step. In the next step, molar concentrationsCtrnpj are transferred to IPHREEQC. IPHREEQC is a class of the PHREEQC geochemical solver that facilitates the usage of PHREEQC along with other software such as transport solvers [43]. After importing concentrations Ctrnpj from the transport step into the geochemical solver, new concentrations Cphrqj are calculated by the PHREEQC. Finally, Rtotalj is calculated using
(15)Rtotaljx,t=Cphrqj−Ctrnpj∆t.All these steps are done using a PHREEQC wrapper file which is part of the LBM code used in this study and has been detailed in Patel’s work [34] and in open-source code Yantra 1.0.0 which is available online (https://bitbucket.org/yantralbm/yantra/overview).The geometry is updated using static update rules detailed in Patel et al. [30]. Before performing the streaming step, to apply the effect of chemical reactions, we update the distribution functions using
(16)gij∗∗x,t=gij∗x,t+wiRtotaljx,tΔt,and then the streaming step is executed:
(17)gijx+ciΔt,t+Δt=gij∗∗x,t.Once the streaming step is performed, species concentrations can be obtained using (13), and we go to the next time step. It should be pointed out that in the next time step, we solve for the velocity when there is at least 1% change in the porosity.


## 3. Results and Discussion

### 3.1. Generation of Initially Mixed and Banded Mineral Assemblages

To investigate the effect of initial mineral distributions on fracture geometry evolution, the quartet structure generation set (QSGS) algorithm [44] was used to generate mixed and banded mineral distributions around a single fracture (Figures 2 and 3). Both of these two geometries are used for both single- and multispecies simulations. In single-species simulations, a synthetic reactive mineral and an inert mineral are present while in the multispecies simulations, the rock matrix is composed of calcite (fast-reacting) and kaolinite (slow-reacting) minerals. The QSGS algorithm is a reconstruction method which can generate the microstructure of materials composed of different constituents such as rocks with different minerals. There are three tunable parameters (cd, Di, and Φ) in the QSGS algorithm to control how different phases (minerals) are growing inside a nongrowing phase. cd represents the probability of a location in the nongrowing phase to become a nucleus for the first growing mineral. In this study, the nonreactive (in a single-species case) and kaolinite minerals (in a multispecies case) were considered as the growing phases which grow inside the reactive (in the single-species case) and calcite minerals (in the multispecies case), respectively. Di is the probability for a nongrowing phase cell in direction i, to merge into a neighboring growing cell. In this study, eight different growth directions are considered for the growing phases. Also, Φ indicates the volume fraction of the growing phase in the nongrowing phase and is considered to be equal to 0.3 in this study. In both mixed and banded structure cases, cd is equal to 0.015. To generate mixed mineral structures (Figure 2), we used Di=1,2,3,4=0.04 and Di=5,6,7,8=0.01. To generate banded structures, values of Di=2,4=0.04 and Di=1,3,5,6,7,8=0 were chosen to allow the vertical growth of growing minerals so that it can generate initial banded structures for the mineral spatial distribution (Figure 3).Figure 2
The geometry of a single fracture (in light gray) surrounded by reactive (in black) and nonreactive (in dark gray) minerals. In multispecies simulations, dark gray is kaolinite and black is calcite. The QSGS algorithm was used to grow nonreactive minerals inside a reactive mineral (in black) to form a mixed structure.Figure 3
The geometry of a single fracture (in light gray) surrounded by reactive (in black) and nonreactive (in dark gray) minerals. In multispecies simulations, dark gray is kaolinite and black is calcite. The QSGS algorithm was used to grow nonreactive minerals inside a reactive mineral (in black) to form a banded structure.


### 3.2. Single-Species Fracture Dissolution

In this section, it is assumed that the rock matrix around the fracture consists of two different minerals. For the single-species case, it is considered that one mineral is reactive (with molar volume of 2 × 10−3 m3/mol) and the other is inert. A single synthetic species A(aq) is injected into the fracture and reacts with the reactive mineral A(s) according to a first-order kinetic reaction with dissolution rate given by the following equation [45]:
(18)RA=krscA,where cA is the concentration of Aaq and s is the reactive surface area.


#### 3.2.1. Initial Mineral Distribution: Mixed Structures

Figure2 is used as the geometry for the mixed-structure case. The dimension of the domain is 1 mm × 0.5 mm with a resolution of 5 μm. Initially, concentration of species A(aq) is zero in the fracture. The flow with kinematic viscosity of 10−6 m2/s (representing water) is driven using a constant pressure gradient in the x-direction and with no flow top and bottom boundaries. Species A(aq) is injected into the fracture from the inlet where it is a constant concentration boundary with a concentration equal to 0.5 mol/L and all other boundaries are zero-gradient concentration boundaries. To see how hydrodynamic properties (permeability and porosity) of the fracture change during dissolution under different transport and reaction conditions, we run simulations at different Peclet (Pe) and Damkohler (Da) numbers defined as follows:
(19)Pe=ULD,Da=krLD,where L is the initial aperture of the fracture equal to 200 μm, D is the diffusion coefficient with a value of 10−11 m2/s, kr is the reaction rate constant obtained from kr=5×10−8·Da m/s, and U is the average velocity calculated by U=5×10−8·Pe m/s. Normalized permeability (permeability divided by initial permeability) and their corresponding normalized porosity (porosity divided by initial porosity) values have been obtained for different Pe and Da numbers. As shown in Figure 4, for a constant Da, higher Pe will lead to a higher permeability increase for a given porosity. As the flow regime will change from a diffusive to advective regime, it results in faster flow of the reactant through the fracture consequently leading to more dissolution of minerals along the fracture walls and quicker increase in fracture aperture. Figure 5 compares the fracture evolution for two Pe numbers of 26 and 0.026 while Da is 20 for both cases. These two cases have the same porosity (after dissolution) of 0.59, but permeability values (after dissolution) in Figures 5(a) and 5(b) are 1.60 × 10−9 and 1.48 × 10−9 m2, respectively. In Figure 5(a), the dissolution front quickly moves along the fracture increasing the aperture along the fracture wall which translates into higher fracture permeability. However, in Figure 5(b), most of the dissolution happens near to the inlet, and it does not increase fracture aperture as much as Figure 5(a) along the flow direction.Figure 4
Normalized permeability versus normalized porosity at different Pe and Da numbers for initially mixed mineral assemblage.Figure 5
Concentration (mol/L) of speciesA(aq) when (a) Pe = 26 and Da = 20 and (b) Pe = 0.026 and Da = 20. In both (a) and (b), porosity (after dissolution) is near 0.59 but permeability (after dissolution) is 1.60 × 10−9 m2 for (a) and 1.48 × 10−9 m2 for (b). The nonreactive minerals are illustrated in light gray. The areas showing a zero concentration of species A(aq) are representing reactive minerals.
(a)
(b)Figure4 also shows how increasing the reaction rate (Da number) results in lower permeability values. Figures 6(a) and 6(b) are related to the cases with Da = 20 and Da = 10, respectively, while both have the same porosity of 0.50 (after dissolution) and Pe equal to 26. As it can be observed, in Figure 6(b), the fracture walls have been dissolved more than Figure 6(a) along the flow direction; hence, dissolution in Figure 6(b) has more effect on increasing the fracture permeability, and permeability values confirm this difference where permeability in Figure 6(a) is 1.47 × 10−9 m2 and in Figure 6(b) is 1.56 × 10−9 m2. In all previous cases, either when Pe is constant and Da differs or when Da is constant and Pe changes, dissolution of reactive minerals in regions far away from the main channel has less impact on increasing permeability since the generated porous zone does not play an important role in increasing permeability and flow regime is more diffusive in this porous layer.Figure 6
Concentration (mol/L) of speciesA(aq) when (a) Pe = 26 and Da = 20 and (b) Pe = 26 and Da = 10. In both (a) and (b), porosity (after dissolution) is near 0.50 but permeability (after dissolution) is 1.47 × 10−9 m2 for (a) and 1.56 × 10−9 m2 for (b). The nonreactive minerals are illustrated in light gray. The areas showing a zero concentration of species A(aq) are representing reactive minerals.
(a)
(b)We also ran a case in which distribution of inert minerals inside the reactive minerals form finer mixed structures. Figure7 shows that if the initial mixed distribution of nonreactive minerals inside the reactive minerals contains finer nonreactive minerals, the porous region formed after dissolution has more negative effect on permeability enhancement. Figures 8(a) and 8(b) are snapshots of species concentration having a porosity of 0.56 (after dissolution) and permeabilities (after dissolution) equal to 1.56 × 10−9 and 1.41 × 10−9 m2, respectively. Although the extent of dissolution (along the fracture wall and in the flow direction) in Figure 8(b) is even a bit more than Figure 8(a), the generated porous zone in Figure 8(b) creates more resistance to flow, and this region would be more tortuous than a porous region containing coarser inert minerals and that is why we observe lower permeability for the case with finer nonreactive minerals. This effect can be observed in Figures 9(a) and 9(b) where the velocity magnitude (in the degraded zone) for the case with finer minerals is less than that of the case with coarser minerals.Figure 7
Normalized permeability versus normalized porosity at Pe = 26 and Da = 20 for initially mixed mineral assemblage where two different mineral sizes are used for nonreactive minerals.Figure 8
Concentration (mol/L) of speciesA(aq) when Pe = 26 and Da = 20, for (a) coarser and (b) finer nonreactive minerals. In both (a) and (b), porosity (after dissolution) is near 0.56 but permeability (after dissolution) is 1.56 × 10−9 m2 for (a) and 1.41 × 10−9 m2 for (b). The nonreactive minerals are illustrated in light gray. The areas showing a zero concentration of species A(aq) are representing reactive minerals.
(a)
(b)Figure 9
Velocity profile (m/s) for cases with (a) coarser and (b) finer nonreactive minerals when Pe = 26 and Da = 20. In both (a) and (b), porosity (after dissolution) is near 0.56 but permeability (after dissolution) is 1.56 × 10−9 m2 for (a) and 1.41 × 10−9 m2 for (b).
(a)
(b)


#### 3.2.2. Initial Mineral Distribution: Banded Structures

For the banded structures, the geometry illustrated in Figure3 is used. All the initial and boundary conditions for flow and mass transfer and also the type of reaction are the same as the ones used in the mixed-structure case. Simulations were performed at different Pe and Da numbers. As with Figure 4 (for the mixed structures), Figure 10 also displays that at a constant Da, higher permeability values can be reached when an advective flow regime is more dominant (higher Pe).Figure 10
Normalized permeability versus normalized porosity at different Pe and Da numbers for initially banded mineral assemblage.The contour maps in Figures11(a) and 11(b) compare the fracture aperture along fracture walls (when Da number is 20 and porosity (after dissolution) is 0.58) showing that in Figure 11(a), which corresponds to Pe = 26, the aperture increase is more uniform compared to that in Figure 11(b), with Pe of 0.026, in which the aperture is increasing more locally and near the inlet. Therefore, permeability is higher in Figure 11(a). Also, in the case of constant Pe and different Da numbers, it can be observed in Figure 10 that there is a negative correlation between the Da number and permeability increase. Figure 12(a) which corresponds to a larger Da number indicates a face-dissolution behavior and less permeability while the dissolution in Figure 12(b) is happening uniformly along the fracture wall leading to a higher permeability.Figure 11
Concentration (mol/L) of speciesA(aq) when (a) Pe = 26 and Da = 20 and (b) Pe = 0.026 and Da = 20. In both (a) and (b), porosity (after dissolution) is near 0.58 but permeability (after dissolution) is 1.42 × 10−9 m2 for (a) and 1.37 × 10−9 m2 for (b). The nonreactive minerals are illustrated in light gray. The areas showing a zero concentration of species A(aq) are representing reactive minerals.
(a)
(b)Figure 12
Concentration (mol/L) of speciesA(aq) when (a) Pe = 26 and Da = 20 and (b) Pe = 26 and Da = 10. In both (a) and (b), porosity (after dissolution) is near 0.52 but permeability (after dissolution) is 1.39 × 10−9 m2 for (a) and 1.43 × 10−9 m2 for (b). The nonreactive minerals are illustrated in light gray. The areas showing a zero concentration of species A(aq) are representing reactive minerals.
(a)
(b)In the banded structure case as well as the mixed structures, it can be observed that the permeability increases until it reaches a constant value. Pe and Da numbers can affect the time it takes for the permeability to reach this steady-state condition, but it is the existence of nonreactive minerals that causes this plateau in permeability-porosity curve. This behavior in normalized permeability-porosity relationship indicates that the relationk=k0φ/φ03 (which is normally used in many continuum modeling approaches to relate fracture permeability changes to porosity [39]) is not always valid when minerals with different reactivities exist in the fractured media. This cubic relation is valid for single mineral cases where fracture planes are parallel and smooth. Figure 13 displays normalized permeability versus normalized porosity for a fracture geometry such as Figure 2 with only reactive minerals present in the rock matrix (no nonreactive mineral was considered). Figure 13 shows that in absence of nonreactive minerals, there will be no steady-state behavior in permeability values and normalized permeability-porosity relationship takes a power-law form, and the exponent of this power-law equation can be equal to 3 for a specific Da and Pe number.Figure 13
Normalized permeability versus normalized porosity at different Pe and Da numbers when a single fracture is surrounded only by a single reactive mineral and no nonreactive mineral exists in the rock matrix.


### 3.3. Multispecies Fracture Dissolution

In the previous section, effects of mineral spatial heterogeneity and different reactive transport conditions on fracture evolution were investigated when synthetic single-species reactive and nonreactive minerals were present in the model. In reality, however, the minerals are composed of different species and have more complicated kinetic rate equations which might affect the fracture textural evolution in ways that are not the same as evolutions related to the synthetic minerals. Therefore, in this section, simulations are performed to show the ability of the developed pore-scale reactive transport model to handle more complicated cases where different minerals with more sophisticated kinetic rate equations exist in the model. As discussed before, in this section, we simulate preferential dissolution of calcite minerals while kaolinite minerals are also present in a fractured caprock, under two different initial mineral distributions.


#### 3.3.1. Initial Mineral Distribution: Mixed Structures

In this section, the geometry illustrated in Figure2 was used with the black and dark gray areas representing calcite and kaolinite minerals, respectively. The domain dimensions, grid resolution, type of flow, and mass transfer boundary conditions are similar to the synthetic mineral case. A constant pressure gradient is applied on the domain such that the average velocity in the entire domain equals to U=5×10−7·Pe m/s. The diffusion coefficient for all species is considered to be 10−10 m2/s. The initial and inlet species concentrations are listed in Table 1. The inlet solution is assumed to be a CO2-saturated brine which has a pH around 4 and is undersaturated with respect to calcite with saturation index Ω (Ω = log(Q/K), where Q is the ion activity product and K is the equilibrium constant) of −3. When solving our chemical reaction system in IPHREEQC, there are different speciation and kinetic reactions. Among those reactions, calcite and kaolinite reactions are considered as kinetic reactions. For the calcite reaction, the following reaction rate equation based on the transition state theory (TST) is used [46]:
(20)RCaCO3=k1aH++k2aH2CO3∗+k31−aCa2+aCO32−Keq,CaCO3,which corresponds to the following reaction:
(21)CaCO3s⇌Ca2++CO32−.Table 1
Concentrations used for the initial and inlet solutions.
Parameter
Initial value in the fracture
Value at the inlet
pH
7.4
3.9
Ca
3.401 × 10−3 mol/kgw
5.510 × 10−3 mol/kgw
C
6.975 × 10−3 mol/kgw
1. 334 mol/kgw
Na
1 mol/kgw
1 mol/kgw
Cl
1 mol/kgw
1 mol/kgw
Al
1.777 × 10−6 mol/kgw
—
Si
1.777 × 10−6 mol/kgw
—In (20), ai is the activity of species i; RCaCO3 is the reaction rate; k1, k2, and k3 are reaction rate constants which are equal to 0.89, 5.01 × 10−4, and 6.6 × 10−7 mol/(m2·s), respectively. The equilibrium constant of the reaction (Keq,CaCO3) is 10–8.49.The kaolinite reaction is described as(22)Al2Si2O5OH4s+6H+⇌2Al3++2H4SiO4+H2O,and its overall reaction rate equation is given by [47]:
(23)RAl2Si2O5OH4=kHaH+0.4+kOHaOH−0.31−aAl3+2aH4SiO42aH+6Keq,Al2Si2O5OH40.9,where reaction rate constants kH and kOH are 10–10.8 and 10–15.7 mol/(m2·s), respectively, and reaction equilibrium constant Keq,Al2Si2O5OH4 is 107.435. It should be pointed out that all the reaction rate and equilibrium constants correspond to 25°C.We performed simulations at two different Pe numbers of 2.6 and 0.026. Ca and Al concentrations at three different times are shown in Figures14–17 for these two Pe numbers. As it can be observed in these figures, for the lower Pe number, the dissolution is limited to the inlet (face dissolution) while for the higher Pe number, since the carbonated water can react with more reactive minerals along the fracture length, the dissolution is uniform along the fracture. These two different dissolution patterns (face dissolution and uniform dissolution) cause different permeability values for the fracture. For the face-dissolution case, the aperture increase happens closer to the inlet and a degraded zone forms in this region where it does not contribute to the permeability increase due to the high flow resistance in this layer. However, in the case of uniform dissolution, the aperture increases uniformly throughout the fracture and it has more positive effect on the permeability enhancement. Due to this behavior, it can be seen in Figure 18 that when the flow regime is more diffusive, the permeability is lower than that of a more advective flow regime.Figure 14
Ca concentrations (mol/L) at (a)t=1.5 hr, (b) t=4.9 hr, and (c) t=8.7 hr. In all three cases, Pe = 2.6 and kaolinite is illustrated in light gray.
(a)
(b)
(c)Figure 15
Ca concentrations (mol/L) at (a)t=6.2 hr, (b) t=20.5 hr, and (c) t=42 hr. In all three cases, Pe = 0.026 and kaolinite is illustrated in light gray.
(a)
(b)
(c)Figure 16
Al concentrations (mol/L) at (a)t=1.5 hr, (b) t=4.9 hr, and (c) t=8.7 hr. In all three cases, Pe = 2.6 and kaolinite is illustrated in light gray.
(a)
(b)
(c)Figure 17
Al concentrations (mol/L) at (a)t=6.2 hr, (b) t=20.5 hr, and (c) t=42 hr. In all three cases, Pe = 0.026 and kaolinite is illustrated in light gray.
(a)
(b)
(c)Figure 18
Normalized permeability versus normalized porosity at two different Pe numbers for initially mixed mineral assemblage (multispecies case).Also, Ca and Al concentrations at the outlet are plotted in Figures19 and 20, respectively. For the lower Pe number, Ca concentration increases and then it becomes constant whereas for the higher Pe number, it increases to a maximum and then starts decreasing. This behavior can be attributed to a decrease in calcite dissolution rate. When calcite surfaces are near to the main flow channel, the fluid velocity is higher and the fluid reactivity can be maintained along the fracture channel. Hence, the dissolution rate is higher which leads to the increase in Ca concentration. However, when calcite surfaces dissolve and become further away from the main flow channel, the reactive fluid (here, the low-pH solution) cannot easily flow into the degraded zone because the flow barrier is high in this region. Therefore, the solution in those regions becomes buffered, and the saturation index is closer to zero meaning that the rate of calcite dissolution would be lower and limited by diffusion. So, we observe a decreasing trend in the Ca concentration at the outlet. Also, Figure 20 shows that Al concentration decreases quickly when the inlet solution is introduced to the domain. This is because the inlet solution does not contain Al and also kaolinite dissolution rate is very slow, so the Al concentration will decrease in the domain. But a closer look at Figure 20 indicates that after a sharp decrease at the beginning, Al concentration starts to increase very slowly in which this increase confirms the kaolinite dissolution reaction which is of course very slow at the time scale considered here.Figure 19
Effluent Ca concentrations versus time at two different Pe numbers.Figure 20
Effluent Al concentrations versus time at two different Pe numbers.


#### 3.3.2. Initial Mineral Distribution: Banded Structures

As in Section3.3.1, we simulate two different cases where Pe = 2.6 and Pe = 0.026. Similar dissolution patterns as in the mixed structures can be seen (Figures 21–24); a smaller Pe number results in face dissolution and higher Pe number leads to uniform dissolution. Also, a smaller Pe number translates into lower permeability values compared to the higher permeability values obtained for the higher Pe number (Figure 25). As in the mixed structures, in the banded structure case, when fluid velocity is higher and uniform aperture increase is visible along the fracture walls, dissolution has a more positive role in increasing fracture permeability than the case with lower velocity leading to more localized dissolution near the inlet.Figure 21
Ca concentrations (mol/L) at (a)t=1.3 hr, (b) t=2.8 hr, and (c) t=6.4 hr. In all three cases, Pe = 2.6 and kaolinite is illustrated in light gray.
(a)
(b)
(c)Figure 22
Ca concentrations (mol/L) at (a)t=5.1 hr, (b) t=12.6 hr, and (c) t=29.3 hr. In all three cases, Pe = 0.026 and kaolinite is illustrated in light gray.
(a)
(b)
(c)Figure 23
Al concentrations (mol/L) at (a)t=1.3 hr, (b) t=2.8 hr, and (c) t=6.4 hr. In all three cases, Pe = 2.6 and kaolinite is illustrated in light gray.
(a)
(b)
(c)Figure 24
Al concentrations (mol/L) at (a)t=5.1 hr, (b) t=12.6 hr, and (c) t=29.3 hr. In all three cases, Pe = 0.026 and kaolinite is illustrated in light gray.
(a)
(b)
(c)Figure 25
Normalized permeability versus normalized porosity at two different Pe numbers for an initially banded mineral assemblage (multispecies case).It should be pointed out that in the banded structure, compared to the mixed-structure case, the permeability values are smaller indicating that when nonreactive minerals are distributed as banded structures, they act as a flow barrier and, as can be observed in Figures26(a) and 26(b), generate more resistance to flow compared to the porous zone formed after dissolution in the case of mixed structures.Figure 26
Velocity profile (m/s) for cases with (a) mixed and (b) banded structures when Pe = 2.6 for a multispecies case. In both (a) and (b), porosity (after dissolution) is near 0.47 but permeability (after dissolution) is 1.58 × 10−9 m2 for (a) and 1.43 × 10−9 m2 for (b).
(a)
(b)Also, similar to the mixed-structure case, while calcite dissolves and recedes from the main flow channel, the rate of calcite dissolution decreases due to buffering effect and lower velocity in the areas between two vertical kaolinite bands (comb-tooth zones). This will cause a decrease in Ca concentration after it reaches a maximum at the beginning (Figure27). However, Figure 28 displays that effluent Ca concentrations for the case of banded structures are higher than those of mixed structures. Although in both cases of mixed and banded structures, after calcite dissolution, the calcite surfaces will be further away from the main flow channel, but as Figure 26 showed, in the case of a banded structure, it would be more difficult for the carbonated water to enter the comb-tooth zone compared to the degraded zone. Therefore, in the case of mixed structures, the extent of buffering would be more than the banded structure case, and that is why lower Ca concentrations are observed at the outlet for mixed-structure cases. Also, higher concentrations of Al between thinner kaolinite bands near to the inlet (Figures 23(c) and 24(c)), also confirm the increasing trend in Al concentration indicated in Figure 29.Figure 27
Effluent Ca concentrations versus time at two different Pe numbers.Figure 28
Effluent Ca concentrations versus time for mixed and banded structures when Pe = 2.6.Figure 29
Effluent Al concentrations versus time at two different Pe numbers.


## 3.1. Generation of Initially Mixed and Banded Mineral Assemblages

To investigate the effect of initial mineral distributions on fracture geometry evolution, the quartet structure generation set (QSGS) algorithm [44] was used to generate mixed and banded mineral distributions around a single fracture (Figures 2 and 3). Both of these two geometries are used for both single- and multispecies simulations. In single-species simulations, a synthetic reactive mineral and an inert mineral are present while in the multispecies simulations, the rock matrix is composed of calcite (fast-reacting) and kaolinite (slow-reacting) minerals. The QSGS algorithm is a reconstruction method which can generate the microstructure of materials composed of different constituents such as rocks with different minerals. There are three tunable parameters (cd, Di, and Φ) in the QSGS algorithm to control how different phases (minerals) are growing inside a nongrowing phase. cd represents the probability of a location in the nongrowing phase to become a nucleus for the first growing mineral. In this study, the nonreactive (in a single-species case) and kaolinite minerals (in a multispecies case) were considered as the growing phases which grow inside the reactive (in the single-species case) and calcite minerals (in the multispecies case), respectively. Di is the probability for a nongrowing phase cell in direction i, to merge into a neighboring growing cell. In this study, eight different growth directions are considered for the growing phases. Also, Φ indicates the volume fraction of the growing phase in the nongrowing phase and is considered to be equal to 0.3 in this study. In both mixed and banded structure cases, cd is equal to 0.015. To generate mixed mineral structures (Figure 2), we used Di=1,2,3,4=0.04 and Di=5,6,7,8=0.01. To generate banded structures, values of Di=2,4=0.04 and Di=1,3,5,6,7,8=0 were chosen to allow the vertical growth of growing minerals so that it can generate initial banded structures for the mineral spatial distribution (Figure 3).Figure 2
The geometry of a single fracture (in light gray) surrounded by reactive (in black) and nonreactive (in dark gray) minerals. In multispecies simulations, dark gray is kaolinite and black is calcite. The QSGS algorithm was used to grow nonreactive minerals inside a reactive mineral (in black) to form a mixed structure.Figure 3
The geometry of a single fracture (in light gray) surrounded by reactive (in black) and nonreactive (in dark gray) minerals. In multispecies simulations, dark gray is kaolinite and black is calcite. The QSGS algorithm was used to grow nonreactive minerals inside a reactive mineral (in black) to form a banded structure.


## 3.2. Single-Species Fracture Dissolution

In this section, it is assumed that the rock matrix around the fracture consists of two different minerals. For the single-species case, it is considered that one mineral is reactive (with molar volume of 2 × 10−3 m3/mol) and the other is inert. A single synthetic species A(aq) is injected into the fracture and reacts with the reactive mineral A(s) according to a first-order kinetic reaction with dissolution rate given by the following equation [45]:
(18)RA=krscA,where cA is the concentration of Aaq and s is the reactive surface area.


### 3.2.1. Initial Mineral Distribution: Mixed Structures

Figure2 is used as the geometry for the mixed-structure case. The dimension of the domain is 1 mm × 0.5 mm with a resolution of 5 μm. Initially, concentration of species A(aq) is zero in the fracture. The flow with kinematic viscosity of 10−6 m2/s (representing water) is driven using a constant pressure gradient in the x-direction and with no flow top and bottom boundaries. Species A(aq) is injected into the fracture from the inlet where it is a constant concentration boundary with a concentration equal to 0.5 mol/L and all other boundaries are zero-gradient concentration boundaries. To see how hydrodynamic properties (permeability and porosity) of the fracture change during dissolution under different transport and reaction conditions, we run simulations at different Peclet (Pe) and Damkohler (Da) numbers defined as follows:
(19)Pe=ULD,Da=krLD,where L is the initial aperture of the fracture equal to 200 μm, D is the diffusion coefficient with a value of 10−11 m2/s, kr is the reaction rate constant obtained from kr=5×10−8·Da m/s, and U is the average velocity calculated by U=5×10−8·Pe m/s. Normalized permeability (permeability divided by initial permeability) and their corresponding normalized porosity (porosity divided by initial porosity) values have been obtained for different Pe and Da numbers. As shown in Figure 4, for a constant Da, higher Pe will lead to a higher permeability increase for a given porosity. As the flow regime will change from a diffusive to advective regime, it results in faster flow of the reactant through the fracture consequently leading to more dissolution of minerals along the fracture walls and quicker increase in fracture aperture. Figure 5 compares the fracture evolution for two Pe numbers of 26 and 0.026 while Da is 20 for both cases. These two cases have the same porosity (after dissolution) of 0.59, but permeability values (after dissolution) in Figures 5(a) and 5(b) are 1.60 × 10−9 and 1.48 × 10−9 m2, respectively. In Figure 5(a), the dissolution front quickly moves along the fracture increasing the aperture along the fracture wall which translates into higher fracture permeability. However, in Figure 5(b), most of the dissolution happens near to the inlet, and it does not increase fracture aperture as much as Figure 5(a) along the flow direction.Figure 4
Normalized permeability versus normalized porosity at different Pe and Da numbers for initially mixed mineral assemblage.Figure 5
Concentration (mol/L) of speciesA(aq) when (a) Pe = 26 and Da = 20 and (b) Pe = 0.026 and Da = 20. In both (a) and (b), porosity (after dissolution) is near 0.59 but permeability (after dissolution) is 1.60 × 10−9 m2 for (a) and 1.48 × 10−9 m2 for (b). The nonreactive minerals are illustrated in light gray. The areas showing a zero concentration of species A(aq) are representing reactive minerals.
(a)
(b)Figure4 also shows how increasing the reaction rate (Da number) results in lower permeability values. Figures 6(a) and 6(b) are related to the cases with Da = 20 and Da = 10, respectively, while both have the same porosity of 0.50 (after dissolution) and Pe equal to 26. As it can be observed, in Figure 6(b), the fracture walls have been dissolved more than Figure 6(a) along the flow direction; hence, dissolution in Figure 6(b) has more effect on increasing the fracture permeability, and permeability values confirm this difference where permeability in Figure 6(a) is 1.47 × 10−9 m2 and in Figure 6(b) is 1.56 × 10−9 m2. In all previous cases, either when Pe is constant and Da differs or when Da is constant and Pe changes, dissolution of reactive minerals in regions far away from the main channel has less impact on increasing permeability since the generated porous zone does not play an important role in increasing permeability and flow regime is more diffusive in this porous layer.Figure 6
Concentration (mol/L) of speciesA(aq) when (a) Pe = 26 and Da = 20 and (b) Pe = 26 and Da = 10. In both (a) and (b), porosity (after dissolution) is near 0.50 but permeability (after dissolution) is 1.47 × 10−9 m2 for (a) and 1.56 × 10−9 m2 for (b). The nonreactive minerals are illustrated in light gray. The areas showing a zero concentration of species A(aq) are representing reactive minerals.
(a)
(b)We also ran a case in which distribution of inert minerals inside the reactive minerals form finer mixed structures. Figure7 shows that if the initial mixed distribution of nonreactive minerals inside the reactive minerals contains finer nonreactive minerals, the porous region formed after dissolution has more negative effect on permeability enhancement. Figures 8(a) and 8(b) are snapshots of species concentration having a porosity of 0.56 (after dissolution) and permeabilities (after dissolution) equal to 1.56 × 10−9 and 1.41 × 10−9 m2, respectively. Although the extent of dissolution (along the fracture wall and in the flow direction) in Figure 8(b) is even a bit more than Figure 8(a), the generated porous zone in Figure 8(b) creates more resistance to flow, and this region would be more tortuous than a porous region containing coarser inert minerals and that is why we observe lower permeability for the case with finer nonreactive minerals. This effect can be observed in Figures 9(a) and 9(b) where the velocity magnitude (in the degraded zone) for the case with finer minerals is less than that of the case with coarser minerals.Figure 7
Normalized permeability versus normalized porosity at Pe = 26 and Da = 20 for initially mixed mineral assemblage where two different mineral sizes are used for nonreactive minerals.Figure 8
Concentration (mol/L) of speciesA(aq) when Pe = 26 and Da = 20, for (a) coarser and (b) finer nonreactive minerals. In both (a) and (b), porosity (after dissolution) is near 0.56 but permeability (after dissolution) is 1.56 × 10−9 m2 for (a) and 1.41 × 10−9 m2 for (b). The nonreactive minerals are illustrated in light gray. The areas showing a zero concentration of species A(aq) are representing reactive minerals.
(a)
(b)Figure 9
Velocity profile (m/s) for cases with (a) coarser and (b) finer nonreactive minerals when Pe = 26 and Da = 20. In both (a) and (b), porosity (after dissolution) is near 0.56 but permeability (after dissolution) is 1.56 × 10−9 m2 for (a) and 1.41 × 10−9 m2 for (b).
(a)
(b)


### 3.2.2. Initial Mineral Distribution: Banded Structures

For the banded structures, the geometry illustrated in Figure3 is used. All the initial and boundary conditions for flow and mass transfer and also the type of reaction are the same as the ones used in the mixed-structure case. Simulations were performed at different Pe and Da numbers. As with Figure 4 (for the mixed structures), Figure 10 also displays that at a constant Da, higher permeability values can be reached when an advective flow regime is more dominant (higher Pe).Figure 10
Normalized permeability versus normalized porosity at different Pe and Da numbers for initially banded mineral assemblage.The contour maps in Figures11(a) and 11(b) compare the fracture aperture along fracture walls (when Da number is 20 and porosity (after dissolution) is 0.58) showing that in Figure 11(a), which corresponds to Pe = 26, the aperture increase is more uniform compared to that in Figure 11(b), with Pe of 0.026, in which the aperture is increasing more locally and near the inlet. Therefore, permeability is higher in Figure 11(a). Also, in the case of constant Pe and different Da numbers, it can be observed in Figure 10 that there is a negative correlation between the Da number and permeability increase. Figure 12(a) which corresponds to a larger Da number indicates a face-dissolution behavior and less permeability while the dissolution in Figure 12(b) is happening uniformly along the fracture wall leading to a higher permeability.Figure 11
Concentration (mol/L) of speciesA(aq) when (a) Pe = 26 and Da = 20 and (b) Pe = 0.026 and Da = 20. In both (a) and (b), porosity (after dissolution) is near 0.58 but permeability (after dissolution) is 1.42 × 10−9 m2 for (a) and 1.37 × 10−9 m2 for (b). The nonreactive minerals are illustrated in light gray. The areas showing a zero concentration of species A(aq) are representing reactive minerals.
(a)
(b)Figure 12
Concentration (mol/L) of speciesA(aq) when (a) Pe = 26 and Da = 20 and (b) Pe = 26 and Da = 10. In both (a) and (b), porosity (after dissolution) is near 0.52 but permeability (after dissolution) is 1.39 × 10−9 m2 for (a) and 1.43 × 10−9 m2 for (b). The nonreactive minerals are illustrated in light gray. The areas showing a zero concentration of species A(aq) are representing reactive minerals.
(a)
(b)In the banded structure case as well as the mixed structures, it can be observed that the permeability increases until it reaches a constant value. Pe and Da numbers can affect the time it takes for the permeability to reach this steady-state condition, but it is the existence of nonreactive minerals that causes this plateau in permeability-porosity curve. This behavior in normalized permeability-porosity relationship indicates that the relationk=k0φ/φ03 (which is normally used in many continuum modeling approaches to relate fracture permeability changes to porosity [39]) is not always valid when minerals with different reactivities exist in the fractured media. This cubic relation is valid for single mineral cases where fracture planes are parallel and smooth. Figure 13 displays normalized permeability versus normalized porosity for a fracture geometry such as Figure 2 with only reactive minerals present in the rock matrix (no nonreactive mineral was considered). Figure 13 shows that in absence of nonreactive minerals, there will be no steady-state behavior in permeability values and normalized permeability-porosity relationship takes a power-law form, and the exponent of this power-law equation can be equal to 3 for a specific Da and Pe number.Figure 13
Normalized permeability versus normalized porosity at different Pe and Da numbers when a single fracture is surrounded only by a single reactive mineral and no nonreactive mineral exists in the rock matrix.


## 3.2.1. Initial Mineral Distribution: Mixed Structures

Figure2 is used as the geometry for the mixed-structure case. The dimension of the domain is 1 mm × 0.5 mm with a resolution of 5 μm. Initially, concentration of species A(aq) is zero in the fracture. The flow with kinematic viscosity of 10−6 m2/s (representing water) is driven using a constant pressure gradient in the x-direction and with no flow top and bottom boundaries. Species A(aq) is injected into the fracture from the inlet where it is a constant concentration boundary with a concentration equal to 0.5 mol/L and all other boundaries are zero-gradient concentration boundaries. To see how hydrodynamic properties (permeability and porosity) of the fracture change during dissolution under different transport and reaction conditions, we run simulations at different Peclet (Pe) and Damkohler (Da) numbers defined as follows:
(19)Pe=ULD,Da=krLD,where L is the initial aperture of the fracture equal to 200 μm, D is the diffusion coefficient with a value of 10−11 m2/s, kr is the reaction rate constant obtained from kr=5×10−8·Da m/s, and U is the average velocity calculated by U=5×10−8·Pe m/s. Normalized permeability (permeability divided by initial permeability) and their corresponding normalized porosity (porosity divided by initial porosity) values have been obtained for different Pe and Da numbers. As shown in Figure 4, for a constant Da, higher Pe will lead to a higher permeability increase for a given porosity. As the flow regime will change from a diffusive to advective regime, it results in faster flow of the reactant through the fracture consequently leading to more dissolution of minerals along the fracture walls and quicker increase in fracture aperture. Figure 5 compares the fracture evolution for two Pe numbers of 26 and 0.026 while Da is 20 for both cases. These two cases have the same porosity (after dissolution) of 0.59, but permeability values (after dissolution) in Figures 5(a) and 5(b) are 1.60 × 10−9 and 1.48 × 10−9 m2, respectively. In Figure 5(a), the dissolution front quickly moves along the fracture increasing the aperture along the fracture wall which translates into higher fracture permeability. However, in Figure 5(b), most of the dissolution happens near to the inlet, and it does not increase fracture aperture as much as Figure 5(a) along the flow direction.Figure 4
Normalized permeability versus normalized porosity at different Pe and Da numbers for initially mixed mineral assemblage.Figure 5
Concentration (mol/L) of speciesA(aq) when (a) Pe = 26 and Da = 20 and (b) Pe = 0.026 and Da = 20. In both (a) and (b), porosity (after dissolution) is near 0.59 but permeability (after dissolution) is 1.60 × 10−9 m2 for (a) and 1.48 × 10−9 m2 for (b). The nonreactive minerals are illustrated in light gray. The areas showing a zero concentration of species A(aq) are representing reactive minerals.
(a)
(b)Figure4 also shows how increasing the reaction rate (Da number) results in lower permeability values. Figures 6(a) and 6(b) are related to the cases with Da = 20 and Da = 10, respectively, while both have the same porosity of 0.50 (after dissolution) and Pe equal to 26. As it can be observed, in Figure 6(b), the fracture walls have been dissolved more than Figure 6(a) along the flow direction; hence, dissolution in Figure 6(b) has more effect on increasing the fracture permeability, and permeability values confirm this difference where permeability in Figure 6(a) is 1.47 × 10−9 m2 and in Figure 6(b) is 1.56 × 10−9 m2. In all previous cases, either when Pe is constant and Da differs or when Da is constant and Pe changes, dissolution of reactive minerals in regions far away from the main channel has less impact on increasing permeability since the generated porous zone does not play an important role in increasing permeability and flow regime is more diffusive in this porous layer.Figure 6
Concentration (mol/L) of speciesA(aq) when (a) Pe = 26 and Da = 20 and (b) Pe = 26 and Da = 10. In both (a) and (b), porosity (after dissolution) is near 0.50 but permeability (after dissolution) is 1.47 × 10−9 m2 for (a) and 1.56 × 10−9 m2 for (b). The nonreactive minerals are illustrated in light gray. The areas showing a zero concentration of species A(aq) are representing reactive minerals.
(a)
(b)We also ran a case in which distribution of inert minerals inside the reactive minerals form finer mixed structures. Figure7 shows that if the initial mixed distribution of nonreactive minerals inside the reactive minerals contains finer nonreactive minerals, the porous region formed after dissolution has more negative effect on permeability enhancement. Figures 8(a) and 8(b) are snapshots of species concentration having a porosity of 0.56 (after dissolution) and permeabilities (after dissolution) equal to 1.56 × 10−9 and 1.41 × 10−9 m2, respectively. Although the extent of dissolution (along the fracture wall and in the flow direction) in Figure 8(b) is even a bit more than Figure 8(a), the generated porous zone in Figure 8(b) creates more resistance to flow, and this region would be more tortuous than a porous region containing coarser inert minerals and that is why we observe lower permeability for the case with finer nonreactive minerals. This effect can be observed in Figures 9(a) and 9(b) where the velocity magnitude (in the degraded zone) for the case with finer minerals is less than that of the case with coarser minerals.Figure 7
Normalized permeability versus normalized porosity at Pe = 26 and Da = 20 for initially mixed mineral assemblage where two different mineral sizes are used for nonreactive minerals.Figure 8
Concentration (mol/L) of speciesA(aq) when Pe = 26 and Da = 20, for (a) coarser and (b) finer nonreactive minerals. In both (a) and (b), porosity (after dissolution) is near 0.56 but permeability (after dissolution) is 1.56 × 10−9 m2 for (a) and 1.41 × 10−9 m2 for (b). The nonreactive minerals are illustrated in light gray. The areas showing a zero concentration of species A(aq) are representing reactive minerals.
(a)
(b)Figure 9
Velocity profile (m/s) for cases with (a) coarser and (b) finer nonreactive minerals when Pe = 26 and Da = 20. In both (a) and (b), porosity (after dissolution) is near 0.56 but permeability (after dissolution) is 1.56 × 10−9 m2 for (a) and 1.41 × 10−9 m2 for (b).
(a)
(b)


## 3.2.2. Initial Mineral Distribution: Banded Structures

For the banded structures, the geometry illustrated in Figure3 is used. All the initial and boundary conditions for flow and mass transfer and also the type of reaction are the same as the ones used in the mixed-structure case. Simulations were performed at different Pe and Da numbers. As with Figure 4 (for the mixed structures), Figure 10 also displays that at a constant Da, higher permeability values can be reached when an advective flow regime is more dominant (higher Pe).Figure 10
Normalized permeability versus normalized porosity at different Pe and Da numbers for initially banded mineral assemblage.The contour maps in Figures11(a) and 11(b) compare the fracture aperture along fracture walls (when Da number is 20 and porosity (after dissolution) is 0.58) showing that in Figure 11(a), which corresponds to Pe = 26, the aperture increase is more uniform compared to that in Figure 11(b), with Pe of 0.026, in which the aperture is increasing more locally and near the inlet. Therefore, permeability is higher in Figure 11(a). Also, in the case of constant Pe and different Da numbers, it can be observed in Figure 10 that there is a negative correlation between the Da number and permeability increase. Figure 12(a) which corresponds to a larger Da number indicates a face-dissolution behavior and less permeability while the dissolution in Figure 12(b) is happening uniformly along the fracture wall leading to a higher permeability.Figure 11
Concentration (mol/L) of speciesA(aq) when (a) Pe = 26 and Da = 20 and (b) Pe = 0.026 and Da = 20. In both (a) and (b), porosity (after dissolution) is near 0.58 but permeability (after dissolution) is 1.42 × 10−9 m2 for (a) and 1.37 × 10−9 m2 for (b). The nonreactive minerals are illustrated in light gray. The areas showing a zero concentration of species A(aq) are representing reactive minerals.
(a)
(b)Figure 12
Concentration (mol/L) of speciesA(aq) when (a) Pe = 26 and Da = 20 and (b) Pe = 26 and Da = 10. In both (a) and (b), porosity (after dissolution) is near 0.52 but permeability (after dissolution) is 1.39 × 10−9 m2 for (a) and 1.43 × 10−9 m2 for (b). The nonreactive minerals are illustrated in light gray. The areas showing a zero concentration of species A(aq) are representing reactive minerals.
(a)
(b)In the banded structure case as well as the mixed structures, it can be observed that the permeability increases until it reaches a constant value. Pe and Da numbers can affect the time it takes for the permeability to reach this steady-state condition, but it is the existence of nonreactive minerals that causes this plateau in permeability-porosity curve. This behavior in normalized permeability-porosity relationship indicates that the relationk=k0φ/φ03 (which is normally used in many continuum modeling approaches to relate fracture permeability changes to porosity [39]) is not always valid when minerals with different reactivities exist in the fractured media. This cubic relation is valid for single mineral cases where fracture planes are parallel and smooth. Figure 13 displays normalized permeability versus normalized porosity for a fracture geometry such as Figure 2 with only reactive minerals present in the rock matrix (no nonreactive mineral was considered). Figure 13 shows that in absence of nonreactive minerals, there will be no steady-state behavior in permeability values and normalized permeability-porosity relationship takes a power-law form, and the exponent of this power-law equation can be equal to 3 for a specific Da and Pe number.Figure 13
Normalized permeability versus normalized porosity at different Pe and Da numbers when a single fracture is surrounded only by a single reactive mineral and no nonreactive mineral exists in the rock matrix.


## 3.3. Multispecies Fracture Dissolution

In the previous section, effects of mineral spatial heterogeneity and different reactive transport conditions on fracture evolution were investigated when synthetic single-species reactive and nonreactive minerals were present in the model. In reality, however, the minerals are composed of different species and have more complicated kinetic rate equations which might affect the fracture textural evolution in ways that are not the same as evolutions related to the synthetic minerals. Therefore, in this section, simulations are performed to show the ability of the developed pore-scale reactive transport model to handle more complicated cases where different minerals with more sophisticated kinetic rate equations exist in the model. As discussed before, in this section, we simulate preferential dissolution of calcite minerals while kaolinite minerals are also present in a fractured caprock, under two different initial mineral distributions.


### 3.3.1. Initial Mineral Distribution: Mixed Structures

In this section, the geometry illustrated in Figure2 was used with the black and dark gray areas representing calcite and kaolinite minerals, respectively. The domain dimensions, grid resolution, type of flow, and mass transfer boundary conditions are similar to the synthetic mineral case. A constant pressure gradient is applied on the domain such that the average velocity in the entire domain equals to U=5×10−7·Pe m/s. The diffusion coefficient for all species is considered to be 10−10 m2/s. The initial and inlet species concentrations are listed in Table 1. The inlet solution is assumed to be a CO2-saturated brine which has a pH around 4 and is undersaturated with respect to calcite with saturation index Ω (Ω = log(Q/K), where Q is the ion activity product and K is the equilibrium constant) of −3. When solving our chemical reaction system in IPHREEQC, there are different speciation and kinetic reactions. Among those reactions, calcite and kaolinite reactions are considered as kinetic reactions. For the calcite reaction, the following reaction rate equation based on the transition state theory (TST) is used [46]:
(20)RCaCO3=k1aH++k2aH2CO3∗+k31−aCa2+aCO32−Keq,CaCO3,which corresponds to the following reaction:
(21)CaCO3s⇌Ca2++CO32−.Table 1
Concentrations used for the initial and inlet solutions.
Parameter
Initial value in the fracture
Value at the inlet
pH
7.4
3.9
Ca
3.401 × 10−3 mol/kgw
5.510 × 10−3 mol/kgw
C
6.975 × 10−3 mol/kgw
1. 334 mol/kgw
Na
1 mol/kgw
1 mol/kgw
Cl
1 mol/kgw
1 mol/kgw
Al
1.777 × 10−6 mol/kgw
—
Si
1.777 × 10−6 mol/kgw
—In (20), ai is the activity of species i; RCaCO3 is the reaction rate; k1, k2, and k3 are reaction rate constants which are equal to 0.89, 5.01 × 10−4, and 6.6 × 10−7 mol/(m2·s), respectively. The equilibrium constant of the reaction (Keq,CaCO3) is 10–8.49.The kaolinite reaction is described as(22)Al2Si2O5OH4s+6H+⇌2Al3++2H4SiO4+H2O,and its overall reaction rate equation is given by [47]:
(23)RAl2Si2O5OH4=kHaH+0.4+kOHaOH−0.31−aAl3+2aH4SiO42aH+6Keq,Al2Si2O5OH40.9,where reaction rate constants kH and kOH are 10–10.8 and 10–15.7 mol/(m2·s), respectively, and reaction equilibrium constant Keq,Al2Si2O5OH4 is 107.435. It should be pointed out that all the reaction rate and equilibrium constants correspond to 25°C.We performed simulations at two different Pe numbers of 2.6 and 0.026. Ca and Al concentrations at three different times are shown in Figures14–17 for these two Pe numbers. As it can be observed in these figures, for the lower Pe number, the dissolution is limited to the inlet (face dissolution) while for the higher Pe number, since the carbonated water can react with more reactive minerals along the fracture length, the dissolution is uniform along the fracture. These two different dissolution patterns (face dissolution and uniform dissolution) cause different permeability values for the fracture. For the face-dissolution case, the aperture increase happens closer to the inlet and a degraded zone forms in this region where it does not contribute to the permeability increase due to the high flow resistance in this layer. However, in the case of uniform dissolution, the aperture increases uniformly throughout the fracture and it has more positive effect on the permeability enhancement. Due to this behavior, it can be seen in Figure 18 that when the flow regime is more diffusive, the permeability is lower than that of a more advective flow regime.Figure 14
Ca concentrations (mol/L) at (a)t=1.5 hr, (b) t=4.9 hr, and (c) t=8.7 hr. In all three cases, Pe = 2.6 and kaolinite is illustrated in light gray.
(a)
(b)
(c)Figure 15
Ca concentrations (mol/L) at (a)t=6.2 hr, (b) t=20.5 hr, and (c) t=42 hr. In all three cases, Pe = 0.026 and kaolinite is illustrated in light gray.
(a)
(b)
(c)Figure 16
Al concentrations (mol/L) at (a)t=1.5 hr, (b) t=4.9 hr, and (c) t=8.7 hr. In all three cases, Pe = 2.6 and kaolinite is illustrated in light gray.
(a)
(b)
(c)Figure 17
Al concentrations (mol/L) at (a)t=6.2 hr, (b) t=20.5 hr, and (c) t=42 hr. In all three cases, Pe = 0.026 and kaolinite is illustrated in light gray.
(a)
(b)
(c)Figure 18
Normalized permeability versus normalized porosity at two different Pe numbers for initially mixed mineral assemblage (multispecies case).Also, Ca and Al concentrations at the outlet are plotted in Figures19 and 20, respectively. For the lower Pe number, Ca concentration increases and then it becomes constant whereas for the higher Pe number, it increases to a maximum and then starts decreasing. This behavior can be attributed to a decrease in calcite dissolution rate. When calcite surfaces are near to the main flow channel, the fluid velocity is higher and the fluid reactivity can be maintained along the fracture channel. Hence, the dissolution rate is higher which leads to the increase in Ca concentration. However, when calcite surfaces dissolve and become further away from the main flow channel, the reactive fluid (here, the low-pH solution) cannot easily flow into the degraded zone because the flow barrier is high in this region. Therefore, the solution in those regions becomes buffered, and the saturation index is closer to zero meaning that the rate of calcite dissolution would be lower and limited by diffusion. So, we observe a decreasing trend in the Ca concentration at the outlet. Also, Figure 20 shows that Al concentration decreases quickly when the inlet solution is introduced to the domain. This is because the inlet solution does not contain Al and also kaolinite dissolution rate is very slow, so the Al concentration will decrease in the domain. But a closer look at Figure 20 indicates that after a sharp decrease at the beginning, Al concentration starts to increase very slowly in which this increase confirms the kaolinite dissolution reaction which is of course very slow at the time scale considered here.Figure 19
Effluent Ca concentrations versus time at two different Pe numbers.Figure 20
Effluent Al concentrations versus time at two different Pe numbers.


### 3.3.2. Initial Mineral Distribution: Banded Structures

As in Section3.3.1, we simulate two different cases where Pe = 2.6 and Pe = 0.026. Similar dissolution patterns as in the mixed structures can be seen (Figures 21–24); a smaller Pe number results in face dissolution and higher Pe number leads to uniform dissolution. Also, a smaller Pe number translates into lower permeability values compared to the higher permeability values obtained for the higher Pe number (Figure 25). As in the mixed structures, in the banded structure case, when fluid velocity is higher and uniform aperture increase is visible along the fracture walls, dissolution has a more positive role in increasing fracture permeability than the case with lower velocity leading to more localized dissolution near the inlet.Figure 21
Ca concentrations (mol/L) at (a)t=1.3 hr, (b) t=2.8 hr, and (c) t=6.4 hr. In all three cases, Pe = 2.6 and kaolinite is illustrated in light gray.
(a)
(b)
(c)Figure 22
Ca concentrations (mol/L) at (a)t=5.1 hr, (b) t=12.6 hr, and (c) t=29.3 hr. In all three cases, Pe = 0.026 and kaolinite is illustrated in light gray.
(a)
(b)
(c)Figure 23
Al concentrations (mol/L) at (a)t=1.3 hr, (b) t=2.8 hr, and (c) t=6.4 hr. In all three cases, Pe = 2.6 and kaolinite is illustrated in light gray.
(a)
(b)
(c)Figure 24
Al concentrations (mol/L) at (a)t=5.1 hr, (b) t=12.6 hr, and (c) t=29.3 hr. In all three cases, Pe = 0.026 and kaolinite is illustrated in light gray.
(a)
(b)
(c)Figure 25
Normalized permeability versus normalized porosity at two different Pe numbers for an initially banded mineral assemblage (multispecies case).It should be pointed out that in the banded structure, compared to the mixed-structure case, the permeability values are smaller indicating that when nonreactive minerals are distributed as banded structures, they act as a flow barrier and, as can be observed in Figures26(a) and 26(b), generate more resistance to flow compared to the porous zone formed after dissolution in the case of mixed structures.Figure 26
Velocity profile (m/s) for cases with (a) mixed and (b) banded structures when Pe = 2.6 for a multispecies case. In both (a) and (b), porosity (after dissolution) is near 0.47 but permeability (after dissolution) is 1.58 × 10−9 m2 for (a) and 1.43 × 10−9 m2 for (b).
(a)
(b)Also, similar to the mixed-structure case, while calcite dissolves and recedes from the main flow channel, the rate of calcite dissolution decreases due to buffering effect and lower velocity in the areas between two vertical kaolinite bands (comb-tooth zones). This will cause a decrease in Ca concentration after it reaches a maximum at the beginning (Figure27). However, Figure 28 displays that effluent Ca concentrations for the case of banded structures are higher than those of mixed structures. Although in both cases of mixed and banded structures, after calcite dissolution, the calcite surfaces will be further away from the main flow channel, but as Figure 26 showed, in the case of a banded structure, it would be more difficult for the carbonated water to enter the comb-tooth zone compared to the degraded zone. Therefore, in the case of mixed structures, the extent of buffering would be more than the banded structure case, and that is why lower Ca concentrations are observed at the outlet for mixed-structure cases. Also, higher concentrations of Al between thinner kaolinite bands near to the inlet (Figures 23(c) and 24(c)), also confirm the increasing trend in Al concentration indicated in Figure 29.Figure 27
Effluent Ca concentrations versus time at two different Pe numbers.Figure 28
Effluent Ca concentrations versus time for mixed and banded structures when Pe = 2.6.Figure 29
Effluent Al concentrations versus time at two different Pe numbers.


## 3.3.1. Initial Mineral Distribution: Mixed Structures

In this section, the geometry illustrated in Figure2 was used with the black and dark gray areas representing calcite and kaolinite minerals, respectively. The domain dimensions, grid resolution, type of flow, and mass transfer boundary conditions are similar to the synthetic mineral case. A constant pressure gradient is applied on the domain such that the average velocity in the entire domain equals to U=5×10−7·Pe m/s. The diffusion coefficient for all species is considered to be 10−10 m2/s. The initial and inlet species concentrations are listed in Table 1. The inlet solution is assumed to be a CO2-saturated brine which has a pH around 4 and is undersaturated with respect to calcite with saturation index Ω (Ω = log(Q/K), where Q is the ion activity product and K is the equilibrium constant) of −3. When solving our chemical reaction system in IPHREEQC, there are different speciation and kinetic reactions. Among those reactions, calcite and kaolinite reactions are considered as kinetic reactions. For the calcite reaction, the following reaction rate equation based on the transition state theory (TST) is used [46]:
(20)RCaCO3=k1aH++k2aH2CO3∗+k31−aCa2+aCO32−Keq,CaCO3,which corresponds to the following reaction:
(21)CaCO3s⇌Ca2++CO32−.Table 1
Concentrations used for the initial and inlet solutions.
Parameter
Initial value in the fracture
Value at the inlet
pH
7.4
3.9
Ca
3.401 × 10−3 mol/kgw
5.510 × 10−3 mol/kgw
C
6.975 × 10−3 mol/kgw
1. 334 mol/kgw
Na
1 mol/kgw
1 mol/kgw
Cl
1 mol/kgw
1 mol/kgw
Al
1.777 × 10−6 mol/kgw
—
Si
1.777 × 10−6 mol/kgw
—In (20), ai is the activity of species i; RCaCO3 is the reaction rate; k1, k2, and k3 are reaction rate constants which are equal to 0.89, 5.01 × 10−4, and 6.6 × 10−7 mol/(m2·s), respectively. The equilibrium constant of the reaction (Keq,CaCO3) is 10–8.49.The kaolinite reaction is described as(22)Al2Si2O5OH4s+6H+⇌2Al3++2H4SiO4+H2O,and its overall reaction rate equation is given by [47]:
(23)RAl2Si2O5OH4=kHaH+0.4+kOHaOH−0.31−aAl3+2aH4SiO42aH+6Keq,Al2Si2O5OH40.9,where reaction rate constants kH and kOH are 10–10.8 and 10–15.7 mol/(m2·s), respectively, and reaction equilibrium constant Keq,Al2Si2O5OH4 is 107.435. It should be pointed out that all the reaction rate and equilibrium constants correspond to 25°C.We performed simulations at two different Pe numbers of 2.6 and 0.026. Ca and Al concentrations at three different times are shown in Figures14–17 for these two Pe numbers. As it can be observed in these figures, for the lower Pe number, the dissolution is limited to the inlet (face dissolution) while for the higher Pe number, since the carbonated water can react with more reactive minerals along the fracture length, the dissolution is uniform along the fracture. These two different dissolution patterns (face dissolution and uniform dissolution) cause different permeability values for the fracture. For the face-dissolution case, the aperture increase happens closer to the inlet and a degraded zone forms in this region where it does not contribute to the permeability increase due to the high flow resistance in this layer. However, in the case of uniform dissolution, the aperture increases uniformly throughout the fracture and it has more positive effect on the permeability enhancement. Due to this behavior, it can be seen in Figure 18 that when the flow regime is more diffusive, the permeability is lower than that of a more advective flow regime.Figure 14
Ca concentrations (mol/L) at (a)t=1.5 hr, (b) t=4.9 hr, and (c) t=8.7 hr. In all three cases, Pe = 2.6 and kaolinite is illustrated in light gray.
(a)
(b)
(c)Figure 15
Ca concentrations (mol/L) at (a)t=6.2 hr, (b) t=20.5 hr, and (c) t=42 hr. In all three cases, Pe = 0.026 and kaolinite is illustrated in light gray.
(a)
(b)
(c)Figure 16
Al concentrations (mol/L) at (a)t=1.5 hr, (b) t=4.9 hr, and (c) t=8.7 hr. In all three cases, Pe = 2.6 and kaolinite is illustrated in light gray.
(a)
(b)
(c)Figure 17
Al concentrations (mol/L) at (a)t=6.2 hr, (b) t=20.5 hr, and (c) t=42 hr. In all three cases, Pe = 0.026 and kaolinite is illustrated in light gray.
(a)
(b)
(c)Figure 18
Normalized permeability versus normalized porosity at two different Pe numbers for initially mixed mineral assemblage (multispecies case).Also, Ca and Al concentrations at the outlet are plotted in Figures19 and 20, respectively. For the lower Pe number, Ca concentration increases and then it becomes constant whereas for the higher Pe number, it increases to a maximum and then starts decreasing. This behavior can be attributed to a decrease in calcite dissolution rate. When calcite surfaces are near to the main flow channel, the fluid velocity is higher and the fluid reactivity can be maintained along the fracture channel. Hence, the dissolution rate is higher which leads to the increase in Ca concentration. However, when calcite surfaces dissolve and become further away from the main flow channel, the reactive fluid (here, the low-pH solution) cannot easily flow into the degraded zone because the flow barrier is high in this region. Therefore, the solution in those regions becomes buffered, and the saturation index is closer to zero meaning that the rate of calcite dissolution would be lower and limited by diffusion. So, we observe a decreasing trend in the Ca concentration at the outlet. Also, Figure 20 shows that Al concentration decreases quickly when the inlet solution is introduced to the domain. This is because the inlet solution does not contain Al and also kaolinite dissolution rate is very slow, so the Al concentration will decrease in the domain. But a closer look at Figure 20 indicates that after a sharp decrease at the beginning, Al concentration starts to increase very slowly in which this increase confirms the kaolinite dissolution reaction which is of course very slow at the time scale considered here.Figure 19
Effluent Ca concentrations versus time at two different Pe numbers.Figure 20
Effluent Al concentrations versus time at two different Pe numbers.


## 3.3.2. Initial Mineral Distribution: Banded Structures

As in Section3.3.1, we simulate two different cases where Pe = 2.6 and Pe = 0.026. Similar dissolution patterns as in the mixed structures can be seen (Figures 21–24); a smaller Pe number results in face dissolution and higher Pe number leads to uniform dissolution. Also, a smaller Pe number translates into lower permeability values compared to the higher permeability values obtained for the higher Pe number (Figure 25). As in the mixed structures, in the banded structure case, when fluid velocity is higher and uniform aperture increase is visible along the fracture walls, dissolution has a more positive role in increasing fracture permeability than the case with lower velocity leading to more localized dissolution near the inlet.Figure 21
Ca concentrations (mol/L) at (a)t=1.3 hr, (b) t=2.8 hr, and (c) t=6.4 hr. In all three cases, Pe = 2.6 and kaolinite is illustrated in light gray.
(a)
(b)
(c)Figure 22
Ca concentrations (mol/L) at (a)t=5.1 hr, (b) t=12.6 hr, and (c) t=29.3 hr. In all three cases, Pe = 0.026 and kaolinite is illustrated in light gray.
(a)
(b)
(c)Figure 23
Al concentrations (mol/L) at (a)t=1.3 hr, (b) t=2.8 hr, and (c) t=6.4 hr. In all three cases, Pe = 2.6 and kaolinite is illustrated in light gray.
(a)
(b)
(c)Figure 24
Al concentrations (mol/L) at (a)t=5.1 hr, (b) t=12.6 hr, and (c) t=29.3 hr. In all three cases, Pe = 0.026 and kaolinite is illustrated in light gray.
(a)
(b)
(c)Figure 25
Normalized permeability versus normalized porosity at two different Pe numbers for an initially banded mineral assemblage (multispecies case).It should be pointed out that in the banded structure, compared to the mixed-structure case, the permeability values are smaller indicating that when nonreactive minerals are distributed as banded structures, they act as a flow barrier and, as can be observed in Figures26(a) and 26(b), generate more resistance to flow compared to the porous zone formed after dissolution in the case of mixed structures.Figure 26
Velocity profile (m/s) for cases with (a) mixed and (b) banded structures when Pe = 2.6 for a multispecies case. In both (a) and (b), porosity (after dissolution) is near 0.47 but permeability (after dissolution) is 1.58 × 10−9 m2 for (a) and 1.43 × 10−9 m2 for (b).
(a)
(b)Also, similar to the mixed-structure case, while calcite dissolves and recedes from the main flow channel, the rate of calcite dissolution decreases due to buffering effect and lower velocity in the areas between two vertical kaolinite bands (comb-tooth zones). This will cause a decrease in Ca concentration after it reaches a maximum at the beginning (Figure27). However, Figure 28 displays that effluent Ca concentrations for the case of banded structures are higher than those of mixed structures. Although in both cases of mixed and banded structures, after calcite dissolution, the calcite surfaces will be further away from the main flow channel, but as Figure 26 showed, in the case of a banded structure, it would be more difficult for the carbonated water to enter the comb-tooth zone compared to the degraded zone. Therefore, in the case of mixed structures, the extent of buffering would be more than the banded structure case, and that is why lower Ca concentrations are observed at the outlet for mixed-structure cases. Also, higher concentrations of Al between thinner kaolinite bands near to the inlet (Figures 23(c) and 24(c)), also confirm the increasing trend in Al concentration indicated in Figure 29.Figure 27
Effluent Ca concentrations versus time at two different Pe numbers.Figure 28
Effluent Ca concentrations versus time for mixed and banded structures when Pe = 2.6.Figure 29
Effluent Al concentrations versus time at two different Pe numbers.


## 4. Conclusions

In this study, a previously developed LBM transport solver coupled with IPHREEQC has been extended for simulating advective-diffusive reactive transport. The LBM was used to investigate the effects of initial mineral spatial heterogeneity on textural alteration of a fracture and consequently on its hydrodynamic properties. The LBM coupled with IPHREEQC geochemical solver enables us to model a broad range of different geochemical reactions. To simulate the effects of different mineral spatial distributions, synthetic models were constructed. Results showed that the generated degraded and comb-tooth zones [3] as observed in experiments done by Ellis et al. [2] will cause the permeability of the fracture to increase until a steady-state condition while in most of continuum-scale models, a power-law relation is used for permeability-porosity relationship which can overestimate the permeability values when multiminerals exist in the system. In addition to the existence of different minerals, this permeability-porosity relation is a function of transport and reaction conditions. At a constant Da number, increasing Pe causes higher permeability values, and at constant Pe, increasing the Da number will lead to a permeability decrease.Moreover, to show the ability of the LBM model for realistic cases, we simulated a case with calcite and kaolinite minerals which were assumed to represent two main minerals in a rock matrix around a single fracture inside a caprock. It was observed that the tortuosity in the degraded zone (in the case of initially mixed distributions) and distance from the main flow channel (in the case of initially banded structures) lowers the permeability of the fracture. Also, a decreasing trend in the effluent Ca concentration confirmed the effect of degraded and comb-tooth zones on decreasing the dissolution rate of calcite in those regions.3D simulations need to be done in order to better understand the effects of initial mineral distributions on fracture hydrodynamic properties. The developed LBM can also be used to simulate the 3D geometries as well, and this will be published in a future work.



---

*Source: 6046182-2018-07-18.xml*