Geofluids

Volume 2018, Article ID 6046182, 28 pages

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

^{1}Department of Geosciences, University of Oslo, Postboks 1047, Blindern, Oslo, Norway^{2}Laboratory for Waste Management (LES), Paul Scherrer Institute, 5232 Villigen, Switzerland

Correspondence should be addressed to Hossein Fazeli; on.oiu.oeg@ilezaf.niessoh

Received 29 December 2017; Revised 1 May 2018; Accepted 15 May 2018; Published 18 July 2018

Academic Editor: Xiangzhao Kong

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.

#### 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_{2} 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_{2} 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^{3}) is density, (m/s) is velocity, (Pa) is pressure, (Pa·s) is dynamic viscosity, and (N/m^{3}) 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: