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

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.


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 fastreacting 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][5][6][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 porescale reactive transport modeling approaches, involving geometry changes due to chemical interactions, have been developed [10][11][12][13][14][15][16][17].Lattice Boltzmann method-(LBM-) based models have also been developed which are mostly focusing on single mineral dissolution in nonfractured media [5,[18][19][20][21][22][23][24][25][26][27][28][29][30][31][32][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].

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] ∇ • u = 0, where ρ (kg/m 3 ) is density, u (m/s) is velocity, P (Pa) is pressure, μ (Pa•s) is dynamic viscosity, and F (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 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]: where In the above equations, Δx and Δt are space and time resolutions.f i is the discrete velocity distribution function, f eq i is the equilibrium distribution function, and τ is the relaxation time.Discrete velocities c i and their corresponding set of weighing coefficients w i 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]: , i = 5, 6, 7, 8, Equation (3) consists of two parts named collision and streaming.Particles move with velocity c i to a neighboring node at position x + c i Δ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]: Once these two steps are performed, distribution functions f i 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 C j is the concentration of species j, D j is the diffusion coefficient which is considered constant for all species, and R j represents the source or sink term due to reactions for species j.R j 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 j 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.R j total 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 R j total at the bulk only includes the homogeneous reactions, but at the fluid nodes adjacent to a solid node, R j total 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.

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 C j trnp 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 C j trnp from the transport step into the geochemical solver, new concentrations C j phrq are calculated by the PHREEQC.Finally, R j total 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).

Geofluids
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.

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 singlespecies 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   Geofluids (c d , D i , and Φ) in the QSGS algorithm to control how different phases (minerals) are growing inside a nongrowing phase.c d 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 singlespecies 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.D i 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, c d is equal to 0.015.To generate mixed mineral structures (Figure 2), we used D i=1,2,3,4 = 0 04 and D i=5,6,7,8 = 0 01.To generate banded structures, values of D i=2,4 = 0 04 and D i=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).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 m 3 /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 firstorder kinetic reaction with dissolution rate given by the following equation [45]: where c A is the concentration of A aq and s is the reactive surface area.

Initial Mineral Distribution
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 m 2 /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   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.7 Geofluids 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 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 m 2 and in Figure 6(b) is 1.56 × 10 −9 m 2 .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.8 Geofluids 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 −9 and 1.41 × 10 −9 m 2 , 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     9 Geofluids with finer minerals is less than that of the case with coarser minerals.

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 10 Geofluids permeability while the dissolution in Figure 12(b) is happening uniformly along the fracture wall leading to a higher permeability.
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 k = k 0 φ/φ 0 3 (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 steadystate 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.11 Geofluids 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 U = 5 × 10 −7 • Pe m/s.The diffusion coefficient for all species is considered to be 10 −10 m 2 /s.The initial and inlet species concentrations are listed in Table 1.The inlet solution is assumed to be a CO 2 -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), a i is the activity of species i; R CaCO 3 is the reaction rate; k 1 , k 2 , and k 3 are reaction rate constants which are equal to 0.89, 5.01 × 10 −4 , and 6.6 × 10 −7 mol/ (m 2 •s), respectively.The equilibrium constant of the reaction (K eq,CaCO 3 ) is 10 -8. 49.
The kaolinite reaction is described as and its overall reaction rate equation is given by [47]:        where reaction rate constants k H and k OH are 10 -10.8 and 10 - 15.7 mol/(m 2 •s), respectively, and reaction equilibrium constant K eq,Al 2 Si 2 O 5 OH 4 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 facedissolution 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.
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      19 Geofluids    21 Geofluids 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.
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.
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.

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 combtooth 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.

Geofluids
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 A⇌A s )   where and β n can be calculated by solving the following equation: In the above equations, c is species concentration, c 0 is concentration at the inlet, c eq is initial concentration at the domain, k is reaction rate constant, and D 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 −5 mol/L and 1 × 10 −5 mol/L, respectively.Also, the values used for reaction rate and diffusion constants are 4.8 × 10 −3 m/s and 10 −9 m 2 /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.   25 Geofluids B, C, and D are reacting according to the following decay reactions [50,51].
where R i (i = A, B, C, and D) refers to the reaction rate, c i is species concentration, and k i is the reaction rate constant.As illustrated in Figure 32, species A with constant concentration of 10 −3 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.

Figure 1 :
Figure 1: Flowchart of the LBM reactive transport model.

Figure 2 :
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 :
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.

Figure 4 :
Figure 4: Normalized permeability versus normalized porosity at different Pe and Da numbers for initially mixed mineral assemblage.

Figure 5 :
Figure 5: Concentration (mol/L) of species A (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 m 2 for (a) and 1.48 × 10 −9 m 2 for (b).The nonreactive minerals are illustrated in light gray.The areas showing a zero concentration of species A (aq) are representing reactive minerals.
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 Figures5(a) and 5(b) are 1.60 × 10 −9 and 1.48 × 10 −9 m 2 , respectively.In Figure5(a), the dissolution front quickly moves along

Figure 6 :
Figure 6: Concentration (mol/L) of species A (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 m 2 for (a) and 1.56 × 10 −9 m 2 for (b).The nonreactive minerals are illustrated in light gray.The areas showing a zero concentration of species A (aq) are representing reactive minerals.

Figure 7 :
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 :
Figure 8: Concentration (mol/L) of species A (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 m 2 for (a) and 1.41 × 10 −9 m 2 for (b).The nonreactive minerals are illustrated in light gray.The areas showing a zero concentration of species A (aq) are representing reactive minerals.
) and 9(b) where the velocity magnitude (in the degraded zone) for the case

Figure 10 :
Figure 10: Normalized permeability versus normalized porosity at different Pe and Da numbers for initially banded mineral assemblage.

Figure 11 :
Figure 11: Concentration (mol/L) of species A (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 m 2 for (a) and 1.37 × 10 −9 m 2 for (b).The nonreactive minerals are illustrated in light gray.The areas showing a zero concentration of species A (aq) are representing reactive minerals.

Figure 18 :Figure 19 :
Figure 18: Normalized permeability versus normalized porosity at two different Pe numbers for initially mixed mineral assemblage (multispecies case).

Figure 20 :
Figure 20: Effluent Al concentrations versus time at two different Pe numbers.

Figure 25 :
Figure 25: Normalized permeability versus normalized porosity at two different Pe numbers for an initially banded mineral assemblage (multispecies case).

Figure 27 :
Figure 27: Effluent Ca concentrations versus time at two different Pe numbers.

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,

Figure 32 :
Figure 32: Sketch of domain used for 1D simulation of advection-diffusion-reaction of four species involved in a chain of decay reactions.The left boundary is maintained at constant concentration for different species, and the right boundary is a zero concentration gradient.

Figure 31 :
Figure 31: Steady-state concentrations, obtained by analytical (contour and solid lines) and LB (dashed lines) methods, for a single species diffusing through the left boundary of a square domain and reacting at the top boundary according to a first-order reaction while other boundaries are zero concentration gradients.

8 ×Figure 33 :
Figure33: Breakthrough curves for four different species at the outlet of a domain where species A is injected into the domain and a chain of decay reactions occurs between four species.Different lines show the results obtained with the LBM approach, and symbols indicate the results simulated using PHREEQC.

Table 1 :
Concentrations used for the initial and inlet solutions.