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

Fazeli, Hossein; Patel, Ravi; Hellevang, Helge

doi:https://doi.org/10.1155/2018/6046182

Geofluids

On this page

Abstract Introduction Results and Discussion Conclusions Appendix Conflicts of Interest Acknowledgments References Copyright Related Articles

Special Issue

Contribution of Pore-Scale Approach to Macroscale Geofluids Modelling in Porous Media

View this Special Issue

Research Article | Open Access

Volume 2018 | Article ID 6046182 | https://doi.org/10.1155/2018/6046182

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

Hossein Fazeli,¹Ravi Patel,²and Helge Hellevang¹

Academic Editor: Xiangzhao Kong

Received29 Dec 2017

Revised01 May 2018

Accepted15 May 2018

Published18 Jul 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.

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 CO₂ 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 CO₂ 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] where (kg/m³) is density, (m/s) is velocity, (Pa) is pressure, (Pa·s) is dynamic viscosity, and (N/m³) 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]: where is density of particles with velocity in the -direction at position and time . Also, is the force term and is a collision operator which represents local distribution of because of collisions of particles. It can be shown that (2) can recover the NS equation meaning that if one solves (2) and obtains , then one is able to calculate the macroscopic quantities such as velocity and pressure using . To solve for , (2) must be discretized in time, velocity space, and physical space [37]: where

In the above equations, and are space and time resolutions. is the discrete velocity distribution function, is the equilibrium distribution function, and is the relaxation time. Discrete velocities and their corresponding set of weighing coefficients 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]:

Equation (3) consists of two parts named collision and streaming. Particles move with velocity to a neighboring node at position during time step , and at the same time, they are affected by the collision operator. So, (3) can be broken into two parts [37]:

Once these two steps are performed, distribution functions in the current timestep can be calculated and then one is able to compute the macroscopic fluid density and velocity [37]:

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:

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]: where is the concentration of species , is the diffusion coefficient which is considered constant for all species, and represents the source or sink term due to reactions for species . 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 species j [30] with solves the ADR equation for concentration of species

In the above equations, 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. 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 at the bulk only includes the homogeneous reactions, but at the fluid nodes adjacent to a solid node, 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

Figure 1 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:

We do not include the reaction source-sink term in this step. In the next step, molar concentrations 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 from the transport step into the geochemical solver, new concentrations are calculated by the PHREEQC. Finally, is calculated using

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 and then the streaming step is executed:

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 (, , and ) in the QSGS algorithm to control how different phases (minerals) are growing inside a nongrowing phase. 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. is the probability for a nongrowing phase cell in direction , 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, is equal to 0.015. To generate mixed mineral structures (Figure 2), we used and . To generate banded structures, values of and 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).

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⁻³ m³/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]: where is the concentration of and is the reactive surface area.

3.2.1. Initial Mineral Distribution: Mixed Structures

Figure 2 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⁻⁶ m²/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: where L is the initial aperture of the fracture equal to 200 μm, D is the diffusion coefficient with a value of 10⁻¹¹ m²/s, is the reaction rate constant obtained from m/s, and is the average velocity calculated by 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⁻⁹ and 1.48 × 10⁻⁹ m², 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.

(a)

(b)

Figure 4 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⁻⁹ m² and in Figure 6(b) is 1.56 × 10⁻⁹ m². 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.

(a)

(b)

We also ran a case in which distribution of inert minerals inside the reactive minerals form finer mixed structures. Figure 7 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⁻⁹ and 1.41 × 10⁻⁹ m², 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.

(a)

(b)

(a)

(b)

3.2.2. Initial Mineral Distribution: Banded Structures

For the banded structures, the geometry illustrated in Figure 3 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).

The contour maps in Figures 11(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.

(a)

(b)

(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 relation (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.

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 Figure 2 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 m/s. The diffusion coefficient for all species is considered to be 10⁻¹⁰ m²/s. The initial and inlet species concentrations are listed in Table 1. The inlet solution is assumed to be a CO₂-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]: which corresponds to the following reaction:

In (20), is the activity of species i; is the reaction rate; , , and are reaction rate constants which are equal to 0.89, 5.01 × 10⁻⁴, and 6.6 × 10⁻⁷ mol/(m²·s), respectively. The equilibrium constant of the reaction () is 10^–8.49.

The kaolinite reaction is described as and its overall reaction rate equation is given by [47]: where reaction rate constants and are 10^–10.8 and 10^–15.7 mol/(m²·s), respectively, and reaction equilibrium constant is 10^7.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 Figures 14–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.

(a)

(b)

(c)

(a)

(b)

(c)

(a)

(b)

(c)

(a)

(b)

(c)

Also, Ca and Al concentrations at the outlet are plotted in Figures 19 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.

3.3.2. Initial Mineral Distribution: Banded Structures

As in Section 3.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.

(a)

(b)

(c)

(a)

(b)

(c)

(a)

(b)

(c)

(a)

(b)

(c)

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 Figures 26(a) and 26(b), generate more resistance to flow compared to the porous zone formed after dissolution in the case of mixed structures.

(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 (Figure 27). 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.

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.

Appendix

A. Validation of LBM Model

The applicability of the current LBM approach (without flow solver) used in this work has been previously demonstrated for ion exchange problems [48] or cases in which reactions were at local equilibrium [30]. In this study, we will perform two different benchmark tests to show the performance of our LBM model where kinetic reactions exist in the problem.

A.1 Reaction-Diffusion of a Single Species in a Square Domain

This benchmark shows how a single species A (with constant concentration) diffuses from left boundary into a square domain with size a × b (Figure 30) and reacts with a solid A_s (according to the reaction AA_s) at the top boundary. The governing equation for the stationary condition is [19]

with the boundary conditions

The analytical solution of (A.1) is [49] where and can be calculated by solving the following equation: In the above equations, is species concentration, is concentration at the inlet, is initial concentration at the domain, is reaction rate constant, and is diffusion coefficient. To compute the analytical and numerical solutions of the problem, a domain size of 1 cm × 1 cm was considered which was discretized into 50 × 50 grids for the LBM. Inlet and equilibrium (initial) concentrations are 10 × 10⁻⁵ mol/L and 1 × 10⁻⁵ mol/L, respectively. Also, the values used for reaction rate and diffusion constants are 4.8 × 10⁻³ m/s and 10⁻⁹ m²/s, respectively. As it was mentioned earlier, the reaction part is done using IPHREEQC. Figure 31 shows the steady-state solution obtained using analytical (the contour map and solid lines) and numerical (dashed lines) LBM. As it can be observed, analytical and numerical methods agree very well indicating the ability of the proposed LBM coupled with IPHREEQC.

A.2 Advection-Diffusion-Reaction of Multispecies Involved in a Kinetic Decay-Chain

To further verify the implementation of a coupled LBM-IPHREEQC reactive transport model, a 1D simulation of advection-diffusion-reaction (ADR) was considered in which four different species A, B, C, and D are reacting according to the following decay reactions [50, 51]. where (i = A, B, C, and D) refers to the reaction rate, is species concentration, and is the reaction rate constant. As illustrated in Figure 32, species A with constant concentration of 10⁻³ mol/L is introduced to the domain (with a length of 0.2 m) from inlet while other species concentrations are kept zero at this boundary. All other boundaries have a zero-gradient boundary condition for concentration. The initial concentration for all species is zero. Velocity is constant and equal to 4.8 × 10⁻⁶ m/s and diffusion coefficient is 10⁻⁸ m²/s. Once A enters the domain, a chain of decay reactions occurs according to (A.6)–(A.9). In these equations , , , and are 5 × 10⁻⁴, 2 × 10⁻⁴, 1 × 10⁻⁴, and 5 × 10⁻⁵ 1/s, respectively. The simulation was run until 10⁵ s, and comparisons of the results were made between the reactive transport modules of PHREEQC and LBM-IPHREEQC. The breakthrough curves for different species concentrations have been plotted in Figure 33. Results show a very good agreement between both solvers confirming the correctness of coupling of LBM with IPHREEQC.

Conflicts of Interest

The authors declare that there is no conflict of interest regarding the publication of this paper.

Acknowledgments

Funding for this study was provided by the PROTECT project (no. 233736) financed by the Research Council of Norway, Total E&P, and DEA which is gratefully acknowledged.

References

C.-F. Tsang, J. Birkholzer, and J. Rutqvist, “A comparative review of hydrologic issues involved in geologic storage of CO₂ and injection disposal of liquid waste,” Environmental Geology, vol. 54, no. 8, pp. 1723–1737, 2008.
View at: Publisher Site | Google Scholar
B. Ellis, C. Peters, J. Fitts et al., “Deterioration of a fractured carbonate caprock exposed to CO₂-acidified brine flow,” Greenhouse Gases: Science and Technology, vol. 1, no. 3, pp. 248–260, 2011.
View at: Publisher Site | Google Scholar
J. P. Fitts and C. A. Peters, “Caprock fracture dissolution and CO₂ leakage,” Reviews in Mineralogy and Geochemistry, vol. 77, no. 1, pp. 459–479, 2013.
View at: Publisher Site | Google Scholar
H. Deng, S. Molins, C. Steefel et al., “A 2.5D reactive transport model for fracture alteration simulation,” Environmental Science & Technology, vol. 50, no. 14, pp. 7564–7571, 2016.
View at: Publisher Site | Google Scholar
H. Deng, M. Voltolini, S. Molins et al., “Alteration and erosion of rock matrix bordering a carbonate-rich shale fracture,” Environmental Science & Technology, vol. 51, no. 15, pp. 8861–8868, 2017.
View at: Publisher Site | Google Scholar
F. Gherardi, T. Xu, and K. Pruess, “Numerical modeling of self-limiting and self-enhancing caprock alteration induced by CO₂ storage in a depleted gas reservoir,” Chemical Geology, vol. 244, no. 1-2, pp. 103–129, 2007.
View at: Publisher Site | Google Scholar
F. Liu, P. Lu, C. Zhu, and Y. Xiao, “Coupled reactive flow and transport modeling of CO₂ sequestration in the Mt. Simon sandstone formation, Midwest U.S.A,” International Journal of Greenhouse Gas Control, vol. 5, no. 2, pp. 294–307, 2011.
View at: Publisher Site | Google Scholar
J. P. Nogues, J. P. Fitts, M. A. Celia, and C. A. Peters, “Permeability evolution due to dissolution and precipitation of carbonates using reactive transport modeling in pore networks,” Water Resources Research, vol. 49, no. 9, pp. 6006–6021, 2013.
View at: Publisher Site | Google Scholar
C. I. Steefel, S. Molins, and D. Trebotich, “Pore scale processes associated with subsurface CO2 injection and sequestration,” Reviews in Mineralogy and Geochemistry, vol. 77, no. 1, pp. 259–303, 2013.
View at: Publisher Site | Google Scholar
M. Liu and P. Mostaghimi, “Characterisation of reactive transport in pore-scale correlated porous media,” Chemical Engineering Science, vol. 173, pp. 121–130, 2017.
View at: Publisher Site | Google Scholar
S. Molins, D. Trebotich, G. H. Miller, and C. I. Steefel, “Mineralogical and transport controls on the evolution of porous media texture using direct numerical simulation,” Water Resources Research, vol. 53, no. 5, pp. 3645–3661, 2017.
View at: Publisher Site | Google Scholar
L. Algive, S. Bekri, and O. Vizika, “Pore-network modeling dedicated to the determination of the petrophysical-property changes in the presence of reactive fluid,” SPE Journal, vol. 15, no. 3, pp. 618–633, 2010.
View at: Publisher Site | Google Scholar
C. Z. Qin, S. M. Hassanizadeh, and A. Ebigbo, “Pore-scale network modeling of microbially induced calcium carbonate precipitation: insight into scale dependence of biogeochemical reaction rates,” Water Resources Research, vol. 52, no. 11, pp. 8794–8810, 2016.
View at: Publisher Site | Google Scholar
J. Pereira Nunes, M. Blunt, and B. Bijeljic, “Pore-scale simulation of carbonate dissolution in micro-CT images,” Journal of Geophysical Research: Solid Earth, vol. 121, no. 2, pp. 558–576, 2016.
View at: Publisher Site | Google Scholar
S. Ovaysi and M. Piri, “Pore-space alteration induced by brine acidification in subsurface geologic formations,” Water Resources Research, vol. 50, no. 1, pp. 440–452, 2014.
View at: Publisher Site | Google Scholar
S. Kim and J. C. Santamarina, “Reactive fluid flow in CO₂ storage reservoirs: a 2-D pore network model study,” Greenhouse Gases: Science and Technology, vol. 5, no. 4, pp. 462–473, 2015.
View at: Publisher Site | Google Scholar
C. Varloteaux, M. T. Vu, S. Békri, and P. M. Adler, “Reactive transport in porous media: pore-network model approach compared to pore-scale model,” Physical Review E, vol. 87, no. 2, article 023010, 2013.
View at: Publisher Site | Google Scholar
M. Liu, M. Shabaninejad, and P. Mostaghimi, “Impact of mineralogical heterogeneity on reactive transport modelling,” Computers & Geosciences, vol. 104, pp. 12–19, 2017.
View at: Publisher Site | Google Scholar
Q. Kang, P. C. Lichtner, and D. Zhang, “An improved lattice Boltzmann model for multicomponent reactive transport in porous media at the pore scale,” Water Resources Research, vol. 43, no. 12, 2007.
View at: Publisher Site | Google Scholar
Q. Kang, P. C. Lichtner, and D. Zhang, “Lattice Boltzmann pore-scale model for multicomponent reactive transport in porous media,” Journal of Geophysical Research: Solid Earth, vol. 111, no. B5, 2006.
View at: Publisher Site | Google Scholar
Q. Kang, D. Zhang, S. Chen, and X. He, “Lattice Boltzmann simulation of chemical dissolution in porous media,” Physical Review E, vol. 65, no. 3, article 036318, 2002.
View at: Publisher Site | Google Scholar
E. S. Boek, I. Zacharoudiou, F. Gray, S. M. Shah, J. P. Crawshaw, and J. Yang, “Multiphase-flow and reactive-transport validation studies at the pore scale by use of lattice Boltzmann computer simulations,” SPE Journal, vol. 22, no. 3, 2017.
View at: Publisher Site | Google Scholar
Q. Kang, P. C. Lichtner, H. S. Viswanathan, and A. I. Abdel-Fattah, “Pore scale modeling of reactive transport involved in geologic CO₂ sequestration,” Transport in Porous Media, vol. 82, no. 1, pp. 197–213, 2010.
View at: Publisher Site | Google Scholar
M. Liu and P. Mostaghimi, “Pore-scale modelling of CO₂ storage in fractured coal,” International Journal of Greenhouse Gas Control, vol. 66, pp. 246–253, 2017.
View at: Publisher Site | Google Scholar
M. Liu and P. Mostaghimi, “Pore-scale simulation of dissolution-induced variations in rock mechanical properties,” International Journal of Heat and Mass Transfer, vol. 111, pp. 842–851, 2017.
View at: Publisher Site | Google Scholar
L. Chen, Q. Kang, H. S. Viswanathan, and W. Q. Tao, “Pore-scale study of dissolution-induced changes in hydrologic properties of rocks with binary minerals,” Water Resources Research, vol. 50, no. 12, pp. 9343–9365, 2014.
View at: Publisher Site | Google Scholar
A. Hiorth, E. Jettestuen, L. M. Cathles, and M. V. Madland, “Precipitation, dissolution, and ion exchange processes coupled with a lattice Boltzmann advection diffusion solver,” Geochimica et Cosmochimica Acta, vol. 104, pp. 99–110, 2013.
View at: Publisher Site | Google Scholar
J. Gao, H. Xing, Z. Tian et al., “Reactive transport in porous media for CO₂ sequestration: pore scale modeling using the lattice Boltzmann method,” Computers & Geosciences, vol. 98, pp. 9–20, 2017.
View at: Publisher Site | Google Scholar
Q. Kang, D. Zhang, and S. Chen, “Simulation of dissolution and precipitation in porous media,” Journal of Geophysical Research: Solid Earth, vol. 108, no. B10, 2003.
View at: Publisher Site | Google Scholar
R. A. Patel, J. Perko, D. Jacques, G. De Schutter, K. Van Breugel, and G. Ye, “A versatile pore-scale multicomponent reactive transport approach based on lattice Boltzmann method: application to portlandite dissolution,” Physics and Chemistry of the Earth, Parts A/B/C, vol. 70-71, pp. 127–137, 2014.
View at: Publisher Site | Google Scholar
Z. Tian and J. Wang, “Lattice Boltzmann simulation of CO₂ reactive transport in network fractured media,” Water Resources Research, vol. 53, no. 8, pp. 7366–7381, 2017.
View at: Publisher Site | Google Scholar
Z. Tian and J. Wang, “Lattice Boltzmann simulation of dissolution-induced changes in permeability and porosity in 3D CO₂ reactive transport,” Journal of Hydrology, vol. 557, pp. 276–290, 2018.
View at: Publisher Site | Google Scholar
A. Fakhari, Y. Li, D. Bolster, and K. T. Christensen, “A phase-field lattice Boltzmann model for simulating multiphase flows in porous media: application and comparison to experiments of CO₂ sequestration at pore scale,” Advances in Water Resources, vol. 114, pp. 119–134, 2018.
View at: Publisher Site | Google Scholar
R. Patel, Lattice Boltzmann Method Based Framework for Simulating Physico-Chemical Processes in Heterogeneous Porous Media and Its Application to Cement Paste, Ghent University, Gent, Belgium, 2016.
E. Skurtveit, E. Skurtveit, L. Grande et al., “Mechanical testing and sealing capacity of the Upper Jurassic Draupne Formation, North Sea,” in 49th U.S. Rock Mechanics/Geomechanics Symposium, San Francisco, CA, USA, June-July 2015.
View at: Google Scholar
W. Graebel, Advanced Fluid Mechanics, Academic Press, Cambridge, MA, USA, 2007.
T. Krüger, H. Kusumaatmaja, A. Kuzmin, O. Shardt, G. Silva, and E. M. Viggen, The Lattice Boltzmann Method, Springer, 2017.
View at: Publisher Site
S. Chapman and T. G. Cowling, “The mathematical theory of non-uniform gases: an account of the kinetic theory of viscosity, thermal conduction and diffusion in gases,” Tech. Rep., Cambridge University Press, Cambridge, UK, 1970.
View at: Google Scholar
C. Steefel, C. A. J. Appelo, B. Arora et al., “Reactive transport codes for subsurface environmental simulation,” Computational Geosciences, vol. 19, no. 3, pp. 445–478, 2015.
View at: Publisher Site | Google Scholar
F. Verhaeghe, J. Liu, M. Guo, S. Arnout, B. Blanpain, and P. Wollants, “Dissolution and diffusion behavior of Al₂O₃ in a Ca O–Al₂ O₃–SiO₂ liquid: an experimental-numerical approach,” Applied Physics Letters, vol. 91, no. 12, article 124104, 2007.
View at: Publisher Site | Google Scholar
X. He, N. Li, and B. Goldstein, “Lattice Boltzmann simulation of diffusion-convection systems with surface chemical reaction,” Molecular Simulation, vol. 25, no. 3-4, pp. 145–156, 2000.
View at: Publisher Site | Google Scholar
F. Verhaeghe, S. Arnout, B. Blanpain, and P. Wollants, “Lattice-Boltzmann modeling of dissolution phenomena,” Physical Review E, vol. 73, no. 3, article 036316, 2006.
View at: Publisher Site | Google Scholar
S. R. Charlton and D. L. Parkhurst, “Modules based on the geochemical model PHREEQC for use in scripting and programming languages,” Computers & Geosciences, vol. 37, no. 10, pp. 1653–1663, 2011.
View at: Publisher Site | Google Scholar
M. Wang, J. Wang, N. Pan, and S. Chen, “Mesoscopic predictions of the effective thermal conductivity for microscale random porous media,” Physical Review E, vol. 75, no. 3, article 036702, 2007.
View at: Publisher Site | Google Scholar
P. Szymczak and A. J. Ladd, “Instabilities in the dissolution of a porous matrix,” Geophysical Research Letters, vol. 38, no. 7, 2011.
View at: Publisher Site | Google Scholar
L. Chou, R. M. Garrels, and R. Wollast, “Comparative study of the kinetics and mechanisms of dissolution of carbonate minerals,” Chemical Geology, vol. 78, no. 3-4, pp. 269–282, 1989.
View at: Publisher Site | Google Scholar
L. Li, C. A. Peters, and M. A. Celia, “Upscaling geochemical reaction rates using pore-scale network modeling,” Advances in Water Resources, vol. 29, no. 9, pp. 1351–1370, 2006.
View at: Publisher Site | Google Scholar
R. A. Patel, J. Perko, D. Jaques, G. De Schutter, G. Ye, and K. Van Breugel, “Lattice Boltzmann based multicomponent reactive transport model coupled with geochemical solver for scale simulations,” in Coupled Problems 2013: Proceedings of the 5th International Conference on Computational Methods for Coupled Problems in Science and Engineering, Ibiza, pp. 806–817, Ibiza, Spain, June 2013.
View at: Google Scholar
H. Carlslaw and J. Jaeger, Conduction of Heat in Solids, Oxford Univ. Press, Oxford, 1959.
D. L. Parkhurst, K. L. Kipp, and S. R. Charlton, PHAST Version 2—A program for simulating groundwater flow, solute transport, and multicomponent geochemical reactions, US Geological Survey Techniques and Methods, 2010, 6, A35.
D. L. Parkhurst and L. Wissmeier, “PhreeqcRM: a reaction module for transport simulators based on the geochemical model PHREEQC,” Advances in Water Resources, vol. 83, pp. 176–189, 2015.
View at: Publisher Site | Google Scholar

Copyright

Copyright © 2018 Hossein Fazeli et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

PDF Download Citation

Download other formats

Order printed copies

Views

1150

Downloads

1339

Citations