WATER RESOURCES RESEARCH, VOL. 50, 440–452, doi:10.1002/2013WR014289, 2014
Pore-space alteration induced by brine acidification in subsurface
geologic formations
Saeed Ovaysi1 and Mohammad Piri1
Received 18 June 2013; revised 6 December 2013; accepted 20 December 2013; published 17 January 2014.
[1] A new Lagrangian particle-based method is presented to simulate reactive transport in
natural porous media. This technique is based on Modified Moving Particle Semi-implicit
(MMPS) and takes as input high-resolution voxel images of natural porous media. The flow
field in the medium is computed by solving the incompressible Navier-Stokes equations.
Moreover, a multicomponent ion transport model is coupled with a homogeneous and
heterogeneous reactions module to handle pore-space alteration (i.e., pore-wall dissolution).
The model is first successfully validated against the experimental data available in the
literature. Subsequently, X-ray microtomographic images of two naturally occurring porous
media are used to investigate the impact of reaction kinetics and pore-space topology on
pore-space alteration induced by brine acidification in subsurface conditions. We observed
that at the normal rates of reactions no significant change in porosity and permeability takes
place in the short term. Whereas, higher reaction rates caused major changes in the
macroscopic properties (e.g., porosity and permeability) of the rocks. We also show that
these changes are strongly affected by the rocks’ pore-scale topologies.
Citation: Ovaysi, S., and M. Piri (2014), Pore-space alteration induced by brine acidification in subsurface geologic formations,Water
Resour. Res., 50, 440–452, doi:10.1002/2013WR014289.
1. Introduction acid injection into carbonate rocks [Hoefner and Fogler,
1988; Fredd and Fogler, 1998; Golfier et al., 2002;
[2] Subsequent to CO2 injection into underground geo- Luquot and Gouze, 2009] have provided valuable insights
logical formations, a series of heterogeneous rock/fluid on the optimum conditions under which permeability
chemical reactions take place. These reactions, which also enhancement is best achieved. In particular, Fredd and
account for the permanent mineral trapping of carbon, Miller [2000] present a classification of dissolution patterns
gradually alter the pore-space topology of the rock forma- based on injection flow rate. This subject is later expanded
tion. The resulting pore-space topology alterations translate by Golfier et al. [2002], where a dissolution pattern dia-
into changes in the properties of the medium such as poros- gram is presented to predict a priori the formation of worm-
ity and permeability. These changes affect fluids flow in holes using Peclet (Pe) and Damkohler (Da) dimensionless
the system and may have implications for the process. numbers. Similar classification of both dissolution and pre-
Clearly, it is highly desirable to avoid conditions whereby cipitation patterns is also reported by Daccord et al.
the rock structure undergoes changes favoring CO2 leakage [1993]. Furthermore, several experimental studies have
to the surface. Wellbore stimulation via acid injection in focused on calcite precipitation in glass micromodels
hydrocarbon reservoirs is another example where rock/fluid [Algive et al., 2007; Yoon et al., 2012] and flow cells [Katz
interactions play a critical role. In this case, unlike carbon et al., 2011] to evaluate the impact of calcite deposition
sequestration, the objective is to enhance the permeability pattern on the porosity and permeability of the medium.
of the medium in order to improve production. [4] Application of the above-mentioned findings in field-
[3] These processes, despite their opposite objectives, scale models requires an additional step to develop consti-
are closely related as similar pore-scale physics dictate the tutive relations for Darcy-scale parameters that convey the
outcomes. Therefore, a fundamental examination of the rel- correct physics. This may be achieved by further under-
evant phenomena at the pore scale can shed light on how to standing the phenomena responsible for rock/fluid interac-
successfully implement the above-mentioned processes. tions at the pore scale. For example, Gouze and Luquot
Experimental investigation of pore-space dissolution due to [2011] and Luquot and Gouze [2009] used X-ray microto-
mography to characterize changes in reactive surface area,
porosity, and permeability at the pore scale induced by
1Department of Chemical and Petroleum Engineering, University of CO2 injection into pure calcite rock samples. Pore-scaleWyoming, Laramie, Wyoming, USA.
modeling of reactive transport is an alternative way to
Corresponding author: S. Ovaysi, Institute for Computational Engineer- obtain relevant information to upscale [Kechagia et al.,
ing and Sciences, University of Texas at Austin, 1 University Station 2002; Battiato and Tartakovsky, 2011] these processes.
C0200, ACES, Austin, TX 78712-0027, USA. (s.ovaysi@gmail.com)
Various numerical methods such as lattice Boltzmann
©2013. American Geophysical Union. All Rights Reserved. [Sullivan et al., 2005; Kang et al., 2010; Chen et al.,
0043-1397/14/10.1002/2013WR014289 2013], pore-network modeling [Varloteaux et al., 2013],
440
OVAYSI AND PIRI: PORE-SPACE ALTERATION INDUCED BY BRINE ACIDIFICATION
and level set [Li et al., 2008a] have been used to study min- [7] In this paper, we present an extension of MMPS
eral dissolution and precipitation in two-dimensional and model to simulate pore-level reactive transport of fluids
three-dimensional porous media. Numerical schemes based through naturally occurring porous media in order to enable
on random walks and finite difference formulations have the study of both of the above-mentioned areas. We first
also been employed to study dissolution [Bekri et al., 1995] discuss the porous mediums studied in this work and the
and precipitation [Salles et al., 1993] in artificially gener- numerical discretization that would best suit MMPS. We
ated two-dimensional and three-dimensional porous media then elaborate on the extended MMPS. This is followed by
as well as fractures [Dijk and Berkowitz, 1998]. Conven- model validation against experimental results available in
tional mesh-based methods [Yang et al., 2013] such as the literature. Finally, we present and analyze the simula-
finite volume have been applied to study pore-scale reac- tion results followed by concluding remarks.
tive transport with varied degrees of success [Molins et al.,
2012; Steefel and Maher, 2009].
[5] Despite these efforts, numerical simulation of reactive 2. Pore-Space Representation
transport at the pore scale requires more research in two [8] In accordance with the objectives of this study, we
main areas: (a) inclusion of electrochemical migration in the selected two samples taken from rocks with different pore-
transport equation and (b) development of robust algorithms space topologies, i.e., Gambier limestone (sample A) and
to account for changes in pore-space topology induced by Berea sandstone (sample B). Experimental measurements
rock/fluid chemical reactions. As evidenced by Ovaysi and revealed 7800 mD and 980 mD permeabilities for the
Piri [2013], electrochemical transport of ionic species is a Gambier limestone and Berea sandstone, respectively. The
critically important mechanism that cannot be ignored in a measured porosities of these rocks are 49% and 22.7%,
pore-scale model of reactive transport. Similar observations respectively. Since carbonates are more susceptible to rock/
were made by Molins et al. [2012], where deviation of the fluid reactions under subsurface conditions, we assumed
results from experimental data at low pH was attributed to that both of these rocks are completely composed of cal-
ignoring electrochemical migration in their model. Reliable cium carbonate. This assumption is not consistent with the
simulation of flow and transport in systems with evolving mineralogy of Berea sandstone; nonetheless, a visual
solid boundaries is another area where more research is inspection of the X-ray microtomography images of these
required. In this work, we focus on both of these areas. rocks, i.e., Figures 1 and 2, reveals vastly different pore-
[6] To this end, one can use mesh-based numerical tech- space topologies, see Blumenfeld et al. [2013] and Øren
niques [Berger et al., 2003; LeVeque and Li, 1994; Pember and Bakke [2003] for pore-size distributions. It is the
et al., 1995; Miller and Trebotich, 2012]. These methods, impact of these different pore-space topologies on dissolu-
however, require relatively significant improvements to tion/precipitation patterns that we investigate in this work.
handle changing boundary conditions expected at the pore Therefore, the original mineralogy of these rocks were not
scale in reactive systems. On the other hand, particle-based the main focus of the study. Furthermore, since the changes
methods are potentially, due to their fully Lagrangian in porosity and permeability appear after long periods of
nature, capable of handling evolving rock/fluid boundaries time, we selected smaller samples to significantly reduce
more reliably. Among these techniques, Modified Moving the time span of the simulations. We cut 0.42 mm 3 0.42
Particle Semi-implicit (MMPS) [Ovaysi and Piri, 2010; mm 3 1.06 mm samples from the original images. We then
Ovaysi, 2010; Ovaysi and Piri, 2011, 2012] and Smoothed
Particle Hydrodynamics (SPH) [Monaghan, 2005; Tarta-
kovsky et al., 2008] have been applied to reactive transport
problems. However, there are major differences between
MMPS and SPH. First, the particle-based summations in
SPH impose several assumptions on the quality of the ker-
nel as well as the existence of a full support for the fluid
particles. These restrictions significantly limit the ability of
SPH to reliably deal with open flow boundaries. As a con-
sequence, only periodic boundary conditions have been
handled by SPH. The particle-based summations in MMPS,
on the other hand, are less restrictive insofar as various
boundary conditions on the open flow boundaries can be
easily implemented [Ovaysi and Piri, 2010]. Another major
difference is in the solution strategy: SPH solves the time-
dependent governing equations in a fully explicit manner
and relies on an artificial equation of state to simulate a
quasi incompressible fluid; MMPS solves the governing
equations semiimplicitly. Here, pressure is computed
through the solution of an elliptic equation which is derived
from the continuity equation for incompressible fluids. The
above factors have enabled MMPS to be validated success-
fully against various experimental data [Ovaysi and Piri,
2010, 2011, 2013], whereas no such quantitative valida- Figure 1. A (3.6 mm)3 isosurface of Gambier limestone
tions is available in the literature for SPH. at 5.406 mm resolution.
441
OVAYSI AND PIRI: PORE-SPACE ALTERATION INDUCED BY BRINE ACIDIFICATION
Figure 4. Distribution of static pressure (in Pa) in a 0.42
mm 3 0.42 mm 3 1.06 mm sample of Berea sandstone
(sample B). Only the fluid particles are visualized.
where Ni is the total number of particles in the neighbor-
hood of i and W, kernel, simply gives a higher weight to
the particles close to i than those in distance. To reduce
computational time and limit the number of neighboring
particles, the kernel is usually chosen so that it adopts zero
beyond a certain boundary. This boundary is called kernel
size and denoted by h. One example of such a kernel that is
Figure 2. A (2.138 mm)3 isosurface of Berea sandstone used in the present work reads
at 5.345 mm resolution. 8 2
placed one particle in the center of each voxel of the discre- >> 2r2 12 0 hr <
tized image, which provided a map to mark the solid/fluid <> h 2
particles. Figures 3 and 4 illustrate the void space occupied W 5> > 2r 2ij h (2)by fluid particles in the Gambier and Berea samples, 22 r < h
respectively. >: h 2
0 r h
3. Particle-Based Summations
where Wij5W rij and r5rij5jri2rjj.
[9] To resolve the governing equations of flow and trans- [11] Furthermore, in equations (3)–(5), we present the
port in a system represented by a set of i51; :::;Np particles particle-based approximations of the gradient, divergence,
identified by their distance vector ri, one needs to define and Laplacian operators, respectively, acting on particle i.
fairly accurate approximations of the differential operators
used in those equations. The particle-based approximations d XNi A 2A
presented in this section are those given in the original r j iiA5 2 rj2ri Wij (3)
work by Koshizuka et al. [1995]. N i j rij
[10] Consider a continuous function A(r). A particle- d XNi vj2vi : rj2ri
based approximation of this function can be written as ri:v5N Wij (4)i j r2ij
X 2d XNiNi 2
A Wðjr 2r jÞ ri A5 kN Aj2Ai Wij (5)j i j i j
j
Ai5AðriÞ5 X (1)Ni
Wðjr 2r jÞ where d is the number of spatial dimensions and the param-i j
j eter k is
ð
ðWðrÞr
2dv
k5 V (6)
WðrÞdv
V
[12] Also, N i, particle number density, is defined as
XNi
N i5 Wij (7)
j
Figure 3. Distribution of static pressure (in Pa) in a 0.42 4. Model Description
mm 3 0.42 mm 3 1.06 mm sample of Gambier limestone [13] The transport of a reactive fluid through a porous
(sample A). Only the fluid particles are visualized. medium is a complex process encompassing several
442
OVAYSI AND PIRI: PORE-SPACE ALTERATION INDUCED BY BRINE ACIDIFICATION
!
physicochemical processes. As such, a realistic pore-level
2 q 1 Xn k k
reactive transport model must be able to address fluid flow, r Pd5 r:vE1 r:v Dt (14)Dt Dt
transport of charged and neutral species, and homogeneous k51
and heterogeneous chemical reactions that result in pore- [21] Further details regarding the above-mentioned equa-
space alteration. In the previous studies, we have discussed tions are discussed elsewhere [Ovaysi and Piri, 2010]. The
the first three aspects of the model which we briefly explain solution to these equations will resolve the pore-level flow
here. We then complete the model description by discus- of a single-component nonreactive fluid in porous media.
sing the heterogeneous fluid/rock reactions that lead to the
alteration of the pore space. 4.2. Multicomponent Reactive Transport
[22] In addition to the equations described in section 4.1,
4.1. Fluid Flow transport of a multicomponent reactive fluid is also gov-
[14] The incompressible Navier-Stokes equations govern erned by
the pore-level laminar flow of a Newtonian incompressible
fluid in porous media. Given the Lagrangian nature of MMPS, DCi DCi DCi DCi5 1 1
here we introduce the Lagrangian version of these equations. Dt Dt diffusion Dt electrical coupling Dt chemical reactions
We begin by presenting the equation of motion, i.e., (15)
Dv 1r 1r l52 P 2 P 1 r2v1g (8) where Ci is the molarity of species i.
Dt q d q s q [23] Equation (15) is solved explicitly at each time step
for all species present in the system. Below, we discuss the
where v is the velocity vector, q is density, l is viscosity, g three terms on the right-hand side of this equation, namely
is the gravity vector, Ps is static pressure, and Pd is dynamic diffusion, electrical coupling, and chemical reactions.
pressure. We consider the total pressure to be the summation 4.2.1. Diffusion
of static and dynamic pressures, i.e., P5Ps1Pd . [24] The first term on the right-hand side of this equation
[15] MMPS follows a semiimplicit approach known as refers to the contribution of molecular diffusion to the total
pressure projection method where the total velocity is rep- transport and is computed by
resented by
DCi 2
v5vE1vI (9) 5Di r Ci1r:ðCir ln ciÞ (16)Dt diffusion
where subscripts E and I denote explicit and implicit, where Di is the molecular diffusion coefficient and crespectively. i
is the
activity coefficient of species i. In this work, the activity
[16] vE is computed using
coefficients of the neutral species are assumed to be unityn whereas for the ions we use the Davies activity model
v 5vn
1
1 2 r lP 1 r2v1g Dt (10) [Samson et al., 1999], i.e.,E q s q
pffiffi
Az2 I
where superscript n and Dt denote the previous time step ln ci52
i pffiffi1Cz2i I (17)
and the size of time step, respectively. 11B I
[17] Static pressure in equation (10) is directly calculated
using the external pressure gradient applied on the medium where A, B, and C are constants. zi is the valence of the spe-
while taking into account its pore-space topology as cies i and I, ionic strength, is defined by
implied by XN
2
r2Ps50 (11) I5 Cizi (18)
i
[18] The remainder of equation (8) is solved implicitly to
compute vI where N is the total number of species.
D 4.2.2. Electrical Couplingt
vI52 rPd (12) [25] The different ions present in the system diffuse atq different rates which can lead to an artificial state where the
[19] To do that, the unknown field of Pd must be calcu- solution is electrically charged. To prevent this, the electri-
lated using the continuity equation. However, to compensate cal coupling term must ensure that the electroneutrality
for the errors resulting from the numerical operations, we condition always prevails at every point in the system, i.e.,
present a modified continuity equation that takes into account
minor deviations occurred in previous time steps, i.e., XN
1 X ziCi50 (19)n r:vkDtk i1r:v50 (13)
Dt
k51 [26] The standard approach [Marchand et al., 2001] to
[20] Combining equations (13), (9), and (12), Pd is assure the electroneutrality condition is to solve the
obtained from Poisson-Nernst-Planck equations. However, that approach
443
OVAYSI AND PIRI: PORE-SPACE ALTERATION INDUCED BY BRINE ACIDIFICATION
Table 1. The Homogeneous Reactions Used in This Studya (Re << 1), we believe this assumption is plausible. Fur-
thermore, we assume the pore walls are completely made
Index Reaction log Keq of calcite. Therefore, the reactive surface can easily be cal-
1 H O () H12 1OH2 213.49 culated by first marking the solid particles that interface
2 CO 21H2O () H2CO 3 22.77 with the fluid phase. A subset of these surface particles,
3 H2CO () HCO 21H13 3 26.27
2 () 22 1 i.e., S, can then be identified at the neighborhood of4 HCO 3 CO 3 1H 210.16
any fluid particle through the summations discussed in
aEquilibrium constants are taken from Ellis et al. [2010]. section 3. Also, since each particle represents a cubic grid
of size Dp, original particle spacing, the volume occupied
by each particle is calculated by V5aD3p, where a 0 is a
is computationally expensive and requires restrictions on shrinkage factor describing the size of particle compared to
the size of the time step. In a previous study [Ovaysi and its original size. Also, the reactive surface of particle i, Ai,
Piri, 2013], we presented a new approach that delivers the is calculated through the assumption that particles shrink/
same results at a lower computational cost. We suggest grow homogeneously, i.e.,
computing the electrical coupling term using
XX WijDC N 2i jF
52kziC z
3
i jCj (20) Ai56aiD
2 (23)
Dt p Nelectrical coupling j i
where k is calculated using where F is a subset of particles in the neighborhood of i
that includes only fluid particles. The summations in this
X X equation are introduced to account for only the exposedN N
z C 1Dt z DCi surface of the solid particles.i i i Dt diffusion
i X i X [29] Also, the contribution of kinetic reactions in equa-k5 (21)N N tion (22) for fluid particle i can be calculated using
Dt ziCi z
2
i Ci
i i
DC 23XXNri;j 10
5 Rkj;lAkWik (24)
which is obtained by combining equations (19), (20), and Dt k ViN i kS l
first-order time integration of Ci.
where R is the production rate of species j through the
4.2.3. Chemical Reactions kj;l
kinetic reaction l listed in Table 2 and Nr is the total num-[27] Depending on their rate, chemical reactions can be
ber of kinetic reactions. All the variables of this equation
categorized as equilibrium controlled or kinetics controlled.
are in SI system of units.
In the present work, we study the interaction of brine and
[30] To account for shrinkage/growth of the surface solid
calcite. Therefore, in Tables 2 and 1, we list seven chemical
particles, we first calculate the amount of mass gained/lost
reactions that are believed to take place [Chou et al., 1989].
as a result of the kinetic reactions, i.e.,
We assume the four homogeneous reactions that are listed
in Table 1 take place instantaneously and hence their equi-
librium must be enforced at every point inside the fluid A DtXXNri
M5M 1 R A W (25)
phase. On the other hand, the remaining three reactions, i i N kCaCO 3;l k iki kS l
listed in Table 2, are heterogeneous surface reactions that
proceed slowly according to their respective rate equations.
First, we need to account for the equilibrium and kinetic where M is in mol. Then a which is a measure of the extent
reactions separately. To do that we write of shrinkage/growth can be calculated by
M
DCi DCi DCi a
i
5 (26)
5 1 (22) q D3
Dt CaCO 3 pchemical reactions Dt k Dt eq
4.2.3.1. Kinetic Reactions where qCaCO 527110 mol m
23.
3
[28] The kinetic reactions listed in Table 2 occur at the [31] The shrinkage/growth of the solid particles continu-
surface of the calcite crystals. In the present work, we ously alters the pore space. This phenomenon is accounted
assume that the flow field does not entrain the partially dis- for by introducing a modified kernel used in the particle-
solved solid particles. Since the studied flow rates are low based summations discussed in section 3. Using a similar
Table 2. The Kinetic Reactions Used in This Studya
Index Reaction 2 dCaCO 3dt ðmol cm22s21Þ
1 CaCO 31H
1 ) Ca 211HCO23 8:931025aH1
2 CaCO 31H 21 2
28
2CO 3 ) Ca 12HCO 3 5310 aH CO 2 3
3 CaCO 3 () Ca 211CO2 6:53102113 aH O21:931022aCa 212 aCO 223
aReaction rates are taken from Chou et al. [1989].
444
OVAYSI AND PIRI: PORE-SPACE ALTERATION INDUCED BY BRINE ACIDIFICATION
Table 3. Canonical Form of the Equilibrium Reactions DCi DCi
5 2FiðCÞ (31)
0 Dt Dt
Index Reaction K eqeq
1 H2O2H
1 () OH2 Keq;1 where
2 H2O1CO 2 () H2CO 3 Keq;1Keq;2
3 H2O1CO 22H
1 () HCO23 Keq;1Keq;2Keq;3
4 H2O1CO 222H
1 () CO 223 K DC DC DCeq;1Keq;2Keq;3Keq;4 ðCÞ i i iFi 5 1 1 (32)
Dt diffusion Dt electrical coupling Dt k
concept presented earlier [Ovaysi and Piri, 2010], the can be known explicitly using the concentration values in
modified kernel reads the previous time step.
[36] Combining equations (30) and (31) yields
Wij;modified5aiajWij (27)
DðC51C11C21C31C4Þ
5F5ðCÞ1F1ðCÞ1F2ðCÞ1F3ðCÞ
Dt
4.2.3.2. Equilibrium Reactions
[32] The equilibrium reactions are treated according to 1F4ðCÞ
the method described by Lichtner et al. [1996]. As the first DðC62C12C322C4Þ
step, the equilibrium reactions must be written in their 5F6ðCÞ2F1ðCÞ2F3ðCÞ22F4ðCÞDt
canonical forms shown in Table 3. This leaves the primary
0
and secondary species as listed in Table 4. Consider R DðC71C21C31C Þej 4 5F7ðCÞ1F2ðCÞ1F3ðCÞ1F4ðCÞ
the progress of the canonical reaction j, then for every sec- Dt
ondary species i (33)
DC 0 [37] Using mass action equations for the reactions ini
5R ei; i51; 2; 3; 4 (28)
Dt Table 3, we obtaineq
0
K C c C
[33] Performing mass balance for primary species on the eq;1 5 5 tC15
reactions listed in Table 3 reveals C6c6c1
0Keq;2C5c5C7c7C 5
DC 25 0 0 0 0
52R e12R
c2Ct
e22R e32R e4
Dt eq 0
(34)
Keq;3C5c5C7c7
DC C356 0 0 0
5R 1R C6c6c3e1 e312R e4 (29)
Dt eq 0Keq;4C5c5C7c7Ct
DC C 57 0 0 0 4
52R 2R 2R C
2 2
6c6ce2 e3 e4 4
Dt eq
where Ct is the total molar concentration.
[34] Substituting equation (28) in equation (29), we [38] Using equation (34) in equation (33), we reach at a
derive set of three nonlinear equations. In the present work, we
use Newton-Raphson method to solve for concentration of
DC DC DC DC DC the primary species. Following that, concentration of the5 1 2 3 4
1 1 1 1 50 secondary species can also be known using equation (34).
Dt eq Dt eq Dt eq Dt eq Dt eq [39] All the results presented in the following sections
DC DC are generated using our parallel code on a hybrid multi-6 1 DC3 DC4
2 2 22 50 GPU platform with 54 GPUs.
Dt eq Dt eq Dt eq Dt eq
DC7 DC2 DC3 DC4
1 1 1 50 5. Validation
Dt eq Dt eq Dt eq Dt eq
[40] We validate the above model against a microfluidic
(30)
reactive flow experiment performed by Li et al. [2008b].
For this purpose, we first generate a 4200 mm-long medium
[35] Also, from equations (15) and (22), we know of a cylindrical tube with 500 mm diameter and 20 mm reso-
lution (see Figure 5). To replicate the experimental condi-
Table 4. Primary and Secondary Species tions, we assume the bordering solid particles to the tube in
the middle 4000 mm are completely made of calcite. The
Secondary Index Species Primary Index Species
remaining length of the tube in the inlet and outlet is
1 OH2 5 H2O assumed to be nonreactive in compliance with the experi-
2 H2CO3 6 H
1
mental setup that indicates ceramics. Next, water with
3 HCO23 7 CO2
4 CO 22 10 mM ionic strength (NaCl) and pH5 4 and 5 is injected3 to the tube at two different flow rates, i.e., 4.72 and
445
OVAYSI AND PIRI: PORE-SPACE ALTERATION INDUCED BY BRINE ACIDIFICATION
Table 5. Brine Compositions in mol/L
Component Initial Brine Acidic Brine
H2O 55.56 55.56
H1 1e27 0.1
Na1 0.9 0.9
Ca21 0.1 0.1
Cl2 1.1 1.2
Figure 5. Visualization of Ca21 molar concentration in
the 4000 mm-long reactive portion of a microtube at t5 15
min. The image is obtained from our simulation results in 6.1. Kinetics
this study when injecting a pH5 4 solution at 9.39 lL/min [42] To study the impact of kinetics on pore-space altera-
flow rate. Only the fluid particles are visualized. tion, the acidic brine is injected to sample A under two dif-
ferent conditions. In case I, the kinetic reactions take place
9.39 lL/min. To achieve these flow rates, we apply 0.19 under the normal rates reported in Table 2. Whereas, the
and 0.38 Pa pressure difference across the length of the kinetic reaction rates for case II were 1003 faster than the
tube, respectively. Concentration of Ca21 in the outlet normal rates. Other variables such as the equilibrium con-
stream is then recorded in each case after 15 min and plot- stants and hydrodynamic conditions were identical in both
ted in Figure 6. As shown in this figure, our simulations cases. The applied pressure difference across the length
produced comparable results with the experimental data of the sample is 1 Pa. Under these conditions, we calculate
and within the same range of accuracy of the modeling Pe5 ulD 53:4 and Re5
qul
l 50:03, where u is the mean1
results presented in the same work, i.e., Li et al. [2008b]. It Hinterstitial velocity and l is the specific length calculated
should be stressed that, in simulating the above experiment, using the following analysis [Ovaysi and Piri, 2011].
only the reactions listed in Tables 1 and 2 were considered. Knowing both the permeability and grain size of Berea
sandstone [Øren and Bakke, 2003], a virtual grain size of
6. Results and Discussion 800 mm is estimated for Gambier limestone through
1:4
[41] The samples presented in section 2 are first saturated K / a , where K and a denote permeability and grain
with water at a 1.2 M ionic strength (NaCl and CaCl ) and size, respectively [Shepherd, 1989]. Although grain size is2
pH5 7. Then, at t5 0 a highly acidic brine with pH5 1 is of no physical meaning in a limestone, we use this number
continuously injected to the inlet of the samples. Brine con- as a measure for specific length to determine the hydrody-
centrations are given in Table 5. Since the samples are namics and dispersion regimes under which the simulations
entirely made of calcium carbonate, the heterogeneous are performed. Also, to have a rough estimate of dissolu-
reactions listed in Table 2 are expected to take place at the tion pattern inside the sample, we define
rock/fluid interfaces. We assume the process takes place
under subsurface conditions of 40C and 12.7 MPa [Ellis rate of dissolution wDa5 5 CaCO 3 (35)
et al., 2010], where CO2 exists as a super-critical phase. rate of convection qMQ
Molecular diffusion coefficients of the ions involved are
given in Table 6. where w dMCaCO 5 dt is the rate of rock dissolution and M is3
the molar mass of the dry sample. qM and Q denote molar
density and volume flow rate of the injected acidic brine,
100 respectively. Although, due to variation of permeability
Exp-pH=4 and reactive surface area, Da changes continuously over
90
time, we obtain an average Da5 2.47 and 9.24 for cases I
80 Li et al-pH=4 and II, respectively. Note the above Pe indicates a transi-
70 Present work-pH=4 tional dispersion zone where molecular diffusion coeffi-
cient plays a significant role. This is evidenced by a
60 Exp-pH=5 uniform distribution and hence the absence of concentra-
1
50 Li et al-pH=5 tion gradient of H ions across the width of the pores in
40 Figure 7. Furthermore, a Da5 2.47 implies that even underPresent work-pH=5 the normal reaction rates, transport of acid to the rock
30
20 Table 6. Molecular Diffusion Coefficients of Different Ions in
10 Water at 40C Taken From Ellis et al. [2010] and Newman and
Thomas-Alyea [2004]
0
0 2 4 6 8 10
Species Dm (m
2/s)3 109 Species D 2m (m /s)3 10
9
Flow rate (µL/min)
H1 9.312 CO2 1.94
Figure 6. Comparison of Ca21 concentration in the outlet OH2 5.26 H2CO3 1.5
1 2
after 15 min as measured in the experiments, modeling Na 1.334 HCO 3 1.2421 22
results from Li et al. [2008b], and our simulation results in Ca 0.792 CO 3 0.968
Cl2 2.032
the present work.
446
Ca2+ Concentration (µM)
OVAYSI AND PIRI: PORE-SPACE ALTERATION INDUCED BY BRINE ACIDIFICATION
case I where rock dissolution takes place in regions far
away from the inlet, see Figure 15. This localized dissolu-
tion pattern is in agreement with previous studies where a
similar behavior was observed both at the core scale [Gol-
fier et al., 2002; Gouze and Luquot, 2011] and field scale
[Kalia and Balakotaiah, 2007; Cohen et al., 2008].
Although normal reaction rates yield a uniform dissolution
pattern, the magnitude of this dissolution during the time
span of our simulations is small and its impact on porosity
and permeability of the sample is insignificant. The faster
kinetics, on the other hand, cause significant change in both
porosity and permeability as shown in Figure 9. In this fig-
ure, we plot the normalized porosities and permeabilities of
sample A versus time. The normalized porosities and per-
meabilities are obtained based on their respective values
computed in the absence of chemical reactions. It is shown
that gradual dissolution of the pore walls leads to an imbal-
ance between the inlet and outlet flow rates where the inlet
flow rate has to be greater in order to fill the newly created
pore space. To highlight this, in Figure 9 we have plotted
both the inlet and outlet permeabilities. Note that calcium
carbonate can both precipitate and dissolve. However, the
highly acidic environment inside the medium forces disso-
lution to be dominant, i.e., we did not observe a > 1. This
leads to a net CO2 production in the medium. However, the
CO2 content of the brine remains below its solubility limit
reported by Ellis et al. [2010].
Figure 7. pH distribution at t5 20 s in sample A with (a)
normal kinetics and (b) 1003 heightened kinetics. Only the
fluid particles with pH< 2 are shown.
surface is inferior to the rate by which the acid is consumed
on the rock surface. In other words, the actual rate of
kinetic reactions is controlled by the transport of reactants
and reaction products to and from the rock-fluid surface.
Therefore, a 1003 increase in the rate of kinetic reactions
will not be proportionally translated into higher dissolution
of the rock matrix. Instead, only a 3.73 increase in the rate
of dissolution is observed. Also, notable in this figure is the
slower progress of the H1 plume in case II. This is due to
the higher rate of H1 consumption in the 1003 heightened
kinetic reactions which leaves less H1 to be transported
downstream. Therefore, the acidic front dissolves more and
travels less into the sample. In Figure 8, we illustrate how
this phenomenon affects the dissolution pattern. As clearly
seen, higher reaction rates have lead to more dissolution of
the inlet pore walls in Figure 8. However, this fast rate of
dissolution has consumed most of the H1 ions and, conse-
quently, the pore walls in the downstream are less affected
when compared to Figure 8 where the kinetic reactions are
slower. Comparing the variation of rock surface areas for
the above cases in Figure 15, it becomes clear that the dis-
solution pattern in case II is localized at the inlet of the
sample. In case I, on the other hand, the dissolution process
takes place more uniformly throughout the medium. There-
fore, the rate of rock dissolution in case I does not increase Figure 8. Shrinkage/growth factor a for the solid par-
sharply as in case II. Given that rock dissolution in both ticles in sample A at t5 20 s under (a) normal kinetics and
cases is transport limited, this can be directly attributed to (b) 1003 heightened kinetics. Only the solid particles with
the continuous supply of H1 ions throughout the sample in a < 0:95 are shown.
447
OVAYSI AND PIRI: PORE-SPACE ALTERATION INDUCED BY BRINE ACIDIFICATION
1.25 1.06
K /K walls of the exposed narrow channels creating new activeinlet 0
Koutlet /K0
φ/φ conduits and hence improving the permeability of the0
1.05 medium. The above phenomenon is demonstrated in Figure
1.2
11 where we illustrate the active flow channels in sample B
1.04 at different times during the simulation. Note that dissolu-
1.15 tion of the solid walls in the marked area gives birth to a
1.03 new active flow conduit which connects to the main flow
channel in the medium. Once the acidic brine diffuses
1.1 through the sample and reaches the narrow conduit that
1.02
connects the large pore on top to the main flow channel, it
1.05 gradually widens the conduit. After some time, this wid-
1.01 ened conduit acts as a reliable connection between the large
pore on top and the main flow channel. In Figure 11, a sim-
1 1
0 5 10 15 20 25 ilar phenomenon has created another active flow channel
t (s) just below the marked area. The collective impact of these
new active flow channels translates to a significant
Figure 9. Variation of porosity and permeability in sam- improvement in permeability of sample B as shown in Fig-
ple A during the injection of a highly acidic brine with ure 12. Noteworthy, creation of new active flow channels
1003 faster kinetics. in sample B has lead to a situation where a minor change in
porosity yields a significant change in permeability. This is
6.2. Pore-Space Topology also evidenced by the relatively insignificant amount of
pore-wall dissolution seen in Figure 13 when compared
[43] Sample B presents a medium with drastically differ-
with that in sample A, cf. Figure 8, which is also demon-
ent pore-scale features than that of sample A and hence a
strated in Figures 14 and 15. Also noteworthy in Figure 12,
suitable candidate to investigate the impact of these fea-
the inlet permeability grows much more quickly than the
tures on pore-space alteration. Identical to section 6.1, we
outlet permeability. This is to pull more fluid particles into
inject a highly acidic brine into sample B that is initially
the medium to supplant the dissolved solid walls. However,
saturated with a neutral brine by applying a 1 Pa pressure
since sample B has a much lower permeability than that of
difference across the length of the sample. Furthermore, to
sample A, the gap between the inlet and outlet permeabil-
observe a noticeable alteration of the pore space, the kinetic
ities are wider in Figure 12 when compared with Figure 9.
reactions were sped up 1003. This compares to case II in
[44] Using equation (35) we calculate Da5 62.6 which
the previous section. As expected, the low permeability of
compared to case II in section 6.1 is a much larger number
sample B hinders progress of the acidic brine through the
m and implies a faster rate of rock matrix dissolution. Thissample, see Figure 10. Using a 200 m grain size for Berea
can be mainly attributed to susceptibility of the pore space
sandstone [Øren and Bakke, 2003], we calculated a
in sample B to the creation of new flow channels and,
Re5 0.009 and Pe5 0.1 which places the dispersion
hence, a greater exposure of the rock surface to the acidic
regime well within the region dominated by molecular dif-
brine which further speeds up the dissolution process. The
fusion. Under this regime, the acidic brine diffuses through
existence of large well-connected pores in sample A makes
the fluid channels it comes in contact with regardless of
the creation of new flow channels less likely. Instead, the
their conductivity. This implies that the less conductive
main mechanism by which pore-wall dissolution contrib-
(passive) flow channels have an equal chance of coming in
utes to porosity/permeability alteration is by widening of
contact with the acidic brine as the more conductive
the already existing channels, see Figure 16. Note that this
(active) channels. The acidic brine can then dissolve the
process creates a flow pattern where more fluid has to be
injected into the medium. Mass conservation implies that
more fluid has to pass faster through the same unaltered
conduits downstream of the medium. It is this increase in
velocity that translates to a more actively connected porous
medium in Figure 16b. Futhermore, it is important to note
that the dissolution process gradually smooths the pore
walls, leading to a gradual decrease in the surface area
available for further dissolution. This in turn leads to a
gradual decline in the rate of dissolution. Formation of new
channels can enhance this process by introducing sharp
decreases in the rock surface area. As shown in Figure 15,
the surface area available to the kinetic reactions in sample
B decreases rapidly after 5 s into the simulation which is
caused by the formation of new channels. This process
brings some of the inactive pores in contact with a supply
of acidic brine. Therefore, as evidenced in the above figure,
Figure 10. pH distribution at t5 18 s in sample B with rate of dissolution increases rather sharply after 5 s.
1003 heightened kinetics. Only the fluid particles with This process continues until t5 10 s when the newly cre-
pH< 2 are shown. ated channels have been widened enough to be visible in
448
K/K0
φ/φ0
OVAYSI AND PIRI: PORE-SPACE ALTERATION INDUCED BY BRINE ACIDIFICATION
Figure 11. The active flow channels in sample B at (a) 2 s, (b) 6 s, (c) 10 s, (d) 14 s, and (e) 18 s. Only
the particles faster than 73 1025 m/s are shown. The color bar indicates velocity in m/s.
Figure 11. These channels are then connected to the main After this point, the dissolution process slows down due to
flow channel. Upon establishing the connections, perme- the decreased surface area available for reactions. Note that
ability of the medium increases sharply, see Figure 12. permeability enhancement in sample B after t5 10 s is
mainly due to widening of the pore channels which is the
1.25 1.06
K /K sole mechanism responsible for permeability enhancementinlet 0
Koutlet /K0
φ/φ0 in sample A.
1.05
1.2 6.3. Constitutive Relations for Permeability and
1.04 Surface Area
1.15 [45] The quantities computed in our pore-scale simula-
1.03 tions can be used to propose constitutive relations describ-
ing the evolution of permeability and rock surface area as a
1.1
1.02 function of variation in porosity in the samples studied.a
Luquot and Gouze [2009] propose rr 5
/
/ to account for 1.05 0 0
1.01 the variation of specific surface area r5 Að12/ÞV versus
porosity. In this equation, V is the bulk volume of the sam-
1 1
0 2 4 6 8 10 12 14 16 18 20 ple and subscript 0 denotes the initial conditions. Given
t (s) that V is constant, to find the powerlaw exponent a, in Fig-a
A 12/ /
Figure 12. Variation of porosity and permeability in sam- ure 17 we fit our data to A 50 12/0 / to get a527.05,0
ple B during the injection of a highly acidic brine at 1003 23.48, and 22.77 for case I sample A, case II sample A,
faster kinetics. and case II sample B, respectively. Clearly, reaction surface
449
K/K0
φ/φ0
OVAYSI AND PIRI: PORE-SPACE ALTERATION INDUCED BY BRINE ACIDIFICATION
0
-5e-08
-1e-07
-1.5e-07
Sample A 1X
-2e-07 Sample A 100X
Sample B 100X
-2.5e-07
-3e-07
-3.5e-07
-4e-07
0 0.0002 0.0004 0.0006 0.0008 0.001 0.0012
L (m)
Figure 13. Shrinkage/growth factor a for the solid par- Figure 14. Variation of rock surface area along the length
ticles in sample B at t5 18 s under 1003 heightened
of samples A (cases I and II) and B (case II) at t5 18 s.
kinetics. Only the solid particles with a < 0:95 are shown.
The acidic brine is injected at L5 0.
area in sample A is more sensitive to changes in porosity The results revealed that pore-scale features play a signifi-
for case I. However, this sensitivity is not proportionally cant role on the type of pore-space alterations that are
translated to permeability enhancement as noted earlier. On induced by rock/fluid chemical reactions. Those in turn
the contrary, in sample B where reduction in surface area determine the changes in macroscopic properties such as
versus porosity is much milder, we observe the highest porosity and permeability. We observed that the dissolution
enhancement in permeability, cf. Figure 18. Furthermore, caused by the exposure of pore walls to acidic brine leads
since permeability enhancement in this sample is nonlinear, to the widening of the exposed pores. In the medium where
we found that the power law exponent in K / /b [Carman, the most conductive flow channels are connected to the rest
1937] is a function of porosity variation. Given our earlier of the system through narrow pores, pore wall dissolution
discussion regarding the birth of new flow channels in sam- can significantly enhance permeability. This significant
ple B after 10 s into the simulation, we calculated b5 1.13 enhancement in permeability, however, is not necessarily
for t< 10 s and b5 14.1 for t> 10 s which is consistent accompanied by a significant increase in porosity. On the
with previous studies on this subject [Noiriel et al., 2004]. other hand, if the porous medium is well-connected through
Permeability variation in case II sample A, on the other a network of conductive flow channels, widening of the
hand, is more uniform and we obtained b5 5.76 for the already conductive pores is the main mechanism by which
entire time span of the simulation. One should note that permeability is enhanced. This implies that a considerable
applicability of these relations to larger samples of these change in porosity has to occur before we note a major
rock types is a subject for further investigation. improvement in permeability.
[47] The simulations carried out in this study covered a
dispersion regime where molecular diffusion is dominant.
7. Conclusions
[46] In this study, we presented a technique to model 1 7e-09
pore-space alteration caused by rock/fluid chemical interac- Sample A 1XSample A 100X
tions. The new model is built on the MMPS platform which Sample B 100X 0.98 6e-09
is a Lagrangian, direct pore-scale modeling technique. Two
rock samples, namely Berea sandstone and Gambier lime- 5e-09
0.96
stone, possessing vastly different pore-scale topologies,
were imaged using X-ray microtomography. Furthermore, 4e-09
0.94
since carbonates are more reactive in acidic environments,
3e-09
it was assumed that both of these rocks are made of cal-
cium carbonate only. We then injected a highly acidic brine 0.92 2e-09
into these samples and using normal reaction rates, we did
not observe any noticeable alteration of the pore space in 0.9 1e-09
20 s (real time). Since the main objective of this study was
to present a Lagrangian pore-scale modeling technique that 0.88 0 2 4 6 8 10 12 14 16 18 20
is capable of handling changing rock-fluid boundaries that t (s)
occur during acid injection into carbonate rocks, we artifi-
cially created circumstances where rock-fluid boundaries Figure 15. Evolution of rock surface area (dashed lines
did change; namely, the kinetic rock/fluid reactions were with symbols) and dissolution rate (solid lines) over time in
sped up 1003. This enabled us to reach our goal of investi- samples A (cases I and II) and B (case II). The initial reac-
gating the impact of chemical pore-space alterations on the tion surface areas in samples A and B are 1.28 3 1025 m2
rocks macroscopic properties in a practical amount of time. and 7.18 3 1026 m2, respectively.
450
A/A 20 ΔA (m )
ψCaCO (mol/s)3
OVAYSI AND PIRI: PORE-SPACE ALTERATION INDUCED BY BRINE ACIDIFICATION
Figure 16. The active flow channels in sample A at (a) 2 s and (b) 20 s. Only the particles faster than 3
3 1024 m/s are shown. The color bar indicates velocity in m/s.
This means that all the pore-scale channels that lie on an 1.25
Sample A 100X
equal distance in the pore space from the injection site are Sample B 100X
equally exposed to the acidic brine. At high Peclet num-
1.2
bers, however, the injected acidic brine prefers the more
conductive channels and hence the chances for pore-wall
dissolution in those channels increase. At the same time, 1.15
the less conductive channels are rarely exposed to the
acidic brine. This can significantly modify the dissolution 1.1
pattern in porous media at high Peclet numbers. A future
work should shed light on this aspect of pore space altera-
tion at high Peclet numbers. 1.05
[48] We also presented constitutive relations describing
the variations in permeability and reaction surface area 1
with changes in porosity for the samples studied here. 1 1.005 1.01 1.015 1.02 1.025 1.03
However, we note that one would need to investigate their φ/φ0
applicability for larger samples. Even though, this question
has been partly answered for single component transport in Figure 18. Variation of the rate of rock dissolution over
small pore-scale samples, see Ovaysi and Piri [2011], the time in sample A and B (both case II). Dots represent the
subject still requires more investigation due to complex simulation results whereas the line is obtained by fitting
processes involved in reactive transport. those results to K5c1ð / Þb/ . For sample A, we obtain (c,0
b)5 (0.033, 5.76) whereas for sample A we calculate (c,
b)5 (0.047, 1.13) at t< 10 s and (c, b)5 (2.27, 14.1) at
1 t> 10 s.
0.95 [49] Acknowledgments. We gratefully acknowledge financial support
of Encana and the School of Energy Resources and the Enhanced Oil
0.9 Recovery Institute at the University of Wyoming.
0.85 References
Algive, L., S. Bekri, M. Robin, and O. Vizika (2007), Reactive transport:
Experiments and modeling, paper presented at the International Sympo-
0.8
sium Society of Core Analysts, Institut Français du Petrole, Calgary,
Canada.
0.75 Battiato, I., and D. M. Tartakovsky (2011), Applicability regimes for mac-
Sample A 1X
Sample A 100X roscopic models of reactive transport in porous media, J. Contam.
Sample B 100X
0.7 Hydrol., 120–121, 18–26.
1 1.005 1.01 1.015 1.02 1.025 1.03 Bekri, S., J. F. Thovert, and P. M. Adler (1995), Dissolution of porous
media, Chem. Eng. Sci., 50(17), 2765–2791.
φ/φ0 Berger, M. J., C. Helzel, and R. J. Leveque (2003), h-box methods for the
approximation of hyperbolic conservation laws on irregular grids, SIAM
Figure 17. Variation of reaction surface area versus J. Numer. Anal., 41(3), 893–918.
porosity in samples A (cases I and II) and B (case II). Dots Blumenfeld, R., M. Blunt, B. Bijeljic, E. Boek, O. Gharbi, P. Mostaghimi,
represent the simulation results whereas the lines are and A. Raeini (2013), Imperial College consortium on pore-scale model-
ing, Yearly Progress Report, 3 Jan. Imperial College, London, U. K.
obtained by fitting those results to A 5 12/ / aA 12/ ð/ Þ . Exponent Carman, P. (1937), Fluid flow through granular beds, Trans. Inst. Chem.0 0 0
a is 27.05 (sample A case I), 23.48 (sample A case II), Eng., 15, 150–166.
Chen, L., Q. Kang, B. A. Robinson, Y. L. He, and W. Q. Tao (2013), Pore-
and22.77 (sample B case II). scale modeling of multiphase reactive transport with phase transitions
451
K/K0
K/K0
OVAYSI AND PIRI: PORE-SPACE ALTERATION INDUCED BY BRINE ACIDIFICATION
and dissolution-precipitation processes in closed systems, Phys. Rev. E, Marchand, J., E. Samson, and J. J. Beaudoin (2001), Modeling ion transport
87, 043306. mechanisms in unsaturated porous media, National Research Council
Chou, L., R. M. Garrels, and R. Wollast (1989), Comparative study of the Canada, Report No. NRCC-44698, p. 1–18.
kinetics and mechanisms of dissolution of carbonate minerals, Chem. Miller, G. H., and D. Trebotich (2012), An embedded boundary method for
Geol., 78, 269–282. the Navier-Stokes equations on a time-dependent domain, Commun.
Cohen, C. E., D. Ding, M. Quintard, and B. Bazin (2008), From pore scale Appl. Math. Comput. Sci., 7(1), 1–31.
to wellbore scale: Impact of geometry on wormhole growth in carbonate Monaghan, J. J. (2005), Smoothed particle hydrodynamics, Rep. Prog.
acidization, Chem. Eng. Sci., 63, 3088–3099. Phys., 68, 1703–1759.
Daccord, G., O. Lietard, and R. Lenormand (1993), Chemical dissolution of Molins, S., D. Trebotich, C. I. Steefel, and C. Shen (2012), An investigation
a porous medium by a reactive fluid—II. Convection vs reaction, behav- of the effect of pore scale flow on average geochemical reaction rates
ior diagram, Chem. Eng. Sci., 48(1), 179–186. using direct numerical simulation,Water Resour. Res., 48, W03527, doi:
Dijk, P., and B. Berkowitz (1998), Precipitation and dissolution of reactive 10.1029/2011WR011404.
solutes in fractures,Water Resour. Res., 34(3), 457–470. Newman, J., and K. E. Thomas-Alyea (2004), Electrochemical Systems,
Dullien, F. A. L. (1992), Porous Media: Fluid Transport and Pore Struc- JohnWiley, Hoboken, N. J.
ture, Academic, San Diego, Calif. Øren, P.-E., and S. Bakke (2003), Reconstruction on Berea sandstone
Ellis, B. R., L. E. Crandell, and C. A. Peters (2010), Limitations for brine and pore-scale modeling of wettability effects, J. Pet. Sci. Eng., 39, 177–
acidification due to SO2 co-injection in geologic carbon sequestration, 199.
Int. J. Greenhouse Gas Control, 4, 575–582. Noiriel, C., P. Gouze, and D. Bernard (2004), Investigation of porosity and
Fredd, C. N., and H. S. Fogler (1998), Influence of transport and reaction on permeability effects from microstructure changes during limestone
wormhole formation in porous media, AIChE J., 44(9), 1933–1949. dissolution, Geophys. Res. Lett., 31, L24603, doi:10.1029/
Fredd, C. N., and M. J. Miller (2000), Validation of carbonate matrix stimu- 2004GL021572.
lation models, paper presented at the SPE International Symposium on Ovaysi, S., and M. Piri (2010), Direct pore-level modeling of incompressi-
Formation Damage Control, Soc. of Pet. Eng., Lafayette, La. ble fluid flow in porous media, J. Comput. Phys., 229, 7456–7476.
Golfier, F., C. Zarcone, B. Bazin R. Lenormand, D. Lasseux, and M. Ovaysi, S. (2010), Direct pore-level modeling of fluid flow in porous media,
Quintard (2002), On the ability of a Darcy-scale model to capture worm- PhD dissertation, Dep. Chem. Pet. Eng., Univ. of Wyoming, Wyo.
hole formation during the dissolution of a porous medium, J. Fluid Ovaysi, S., and M. Piri (2011), Pore-scale modeling of dispersion in disor-
Mech., 457, 213–254. dered porous media, J. Contam. Hydrol., 124, 68–81.
Gouze, P., and L. Luquot (2011), X-ray microtomography characterization Ovaysi, S., and M. Piri (2012), Multi-GPU acceleration of direct pore-scale
of porosity, permeability and reactive surface changes during dissolution, modeling of fluid flow in natural porous media, Comput. Phys. Commun.,
J. Contam. Hydrol., 120–121, 45–55. 183, 1890–1898.
Hoefner, M. L., and H. S. Fogler (1988), Pore evolution and channel forma- Ovaysi, S., and M. Piri (2013), Pore-scale dissolution of CO21SO2 in deep
tion during flow and reaction in porous media, AIChE J., 34(1), 45–54. saline aquifers, Int. J. Greenhouse Gas Control, 15, 119–133.
Kalia, N., and V. Balakotaiah (2007), Modeling and analysis of wormhole Pember, R. B., J. B. Bell, P. Colella, W. Y. Crutchfield, and M. L. Welcome
formation in reactive dissolution of carbonate rocks, Chem. Eng. Sci., 62, (1995), An adaptive Cartesian grid method for unsteady compressible
919–928. flow in irregular regions, J. Comput. Phys., 120(2), 278–304.
Kang, Q., P. C. Lichtner, and D. R. Janecky (2010), Lattice Boltzmann Salles, J., J. F. Thovert, and P. M. Adler (1993), Deposition in porous media
method for reacting flows in porous media, Adv. Appl. Math. Mech., 2(5), and clogging, Chem. Eng. Sci., 48(16), 2839–2858.
545–563. Samson, E., G. Lemaire, J. Marchand, and J. J. Beaudoin (1999), Modeling
Katz, G. E., B. Berkowitz, A. Guadagnini, and M. W. Saaltink (2011), chemical activity effects in strong ionic solutions, Comput. Mater. Sci.,
Experimental and modeling investigation of multicomponent reactive 15, 285–294.
transport in porous media, J. Contam. Hydrol., 120–121, 27–44. Shepherd, R. G. (1989), Correlations of permeability and grain size,
Kechagia, P. E., I. N. Tsimpanogiannis, Y. C. Yortsos, and P. Lichtner Ground Water, 27(5), 633–638.
(2002), On the upscaling of reaction-transport processes in porous media Steefel, C. I., and K. Maher (2009), Fluid-rock interaction: A reactive trans-
with fast or finite kinetics, Chem. Eng. Sci., 57, 2565–2577. port approach, Rev. Mineral. Geochem., 70, 485–532.
Koshizuka, S., H. Tamako, and Y. Oka (1995), A particle method for Sullivan, S. P., F. M. Sani, M. L. Johns, and L. F. Gladden (2005), Simula-
incompressible viscous flow with fluid fragmentation, Comput. Fluid tion of packed bed reactors using lattice Boltzmann methods, Chem.
Dyn. J., 4(1), 29–46. Eng. Sci., 60(12), 3405–3418.
LeVeque, R. J., and Z. L. Li (1994), The immersed interface method for Tartakovsky, A. M., G. Redden, P. C. Lichtner, T. D. Scheibe, and P.
elliptic equations with discontinuous coefficients and singular sources, Meakin (2008), Mixing-induced precipitation: Experimental study and
SIAM J. Numer. Anal., 31(4), 1019–1044. multiscale numerical analysis, Water Resour. Res., 44, W06S04, doi:
Li, X., H. Huang, and P. Meakin (2008a), Level-set simulation of coupled 10.1029/2006WR005725.
advection-diffusion and pore structure evolution due to mineral precipi- Varloteaux, C., M. T. Vu, S. Bekri, and P. M. Adler (2013), Reactive trans-
tation in porous media, Water Resour. Res., 44, W12407, doi:10.1029/ port in porous media: Pore-network model approach compared to pore-
2007WR006742. scale model, Phys. Rev. E, 87, 023010.
Li, L., C. I. Steefel, and L. Yang (2008b), Scale dependence of mineral dis- Yang, X. F., T. D. Scheibe, M. C. Richmond, W. A. Perkins, S. J. Vogt, S.
solution rates within single pores and fractures, Geochim. Cosmochim. L. Codd, J. D. Seymour, and M. I. McKinley (2013), Direct numerical
Acta, 72, 360–377. simulation of pore-scale flow in a bead pack: Comparison with magnetic
Lichtner, P. C., C. I. Steefel, and E. H. Oelkers (1996), Reactive Transport resonance imaging observations, Adv. Water Resour., 54, 228–241.
in Porous Media, Mineral. Soc. of Am., Chantilly VA, USA. Yoon, H., A. J. Valocchi, C. J. Werth, and T. Dewers (2012), Pore-scale
Luquot, L., and P. Gouze (2009), Experimental determination of porosity simulation of mixing-induced calcium carbonate precipitation and disso-
and permeability changes induced by injection of CO2 into carbonate lution in a microfluidic pore network, Water Resour. Res., 48, W02524,
rocks, Chem. Geol., 265, 148–159. doi:10.1029/2011WR011192.
452