CO2-induced dissolution of low permeability carbonates. Part II: Numerical modeling of experiments

Advances in Water Resources 62 (2013) 388–408

Contents lists available at ScienceDirect

Advances in Water Resources journal homepage: www.elsevier.com/locate/advwatres

CO2-induced dissolution of low permeability carbonates. Part II: Numerical modeling of experiments Yue Hao ⇑, Megan Smith, Yelena Sholokhova, Susan Carroll Lawrence Livermore National Laboratory, 7000 East Avenue, Livermore, CA 94550, United States

a r t i c l e

i n f o

Article history: Available online 25 September 2013 Keywords: Enhanced oil recovery CO2 sequestration and storage Carbonate dissolution Core flood experiment Numerical simulation Reactive-transport model

a b s t r a c t We used the 3D continuum-scale reactive transport models to simulate eight core flood experiments for two different carbonate rocks. In these experiments the core samples were reacted with brines equilibrated with pCO2 = 3, 2, 1, 0.5 MPa (Smith et al., 2013 [27]). The carbonate rocks were from specific Marly dolostone and Vuggy limestone flow units at the IEAGHG Weyburn-Midale CO2 Monitoring and Storage Project in south-eastern Saskatchewan, Canada. Initial model porosity, permeability, mineral, and surface area distributions were constructed from micro tomography and microscopy characterization data. We constrained model reaction kinetics and porosity–permeability equations with the experimental data. The experimental data included time-dependent solution chemistry and differential pressure measured across the core, and the initial and final pore space and mineral distribution. Calibration of the model with the experimental data allowed investigation of effects of carbonate reactivity, flow velocity, effective permeability, and time on the development and consequences of stable and unstable dissolution fronts. The continuum scale model captured the evolution of distinct dissolution fronts that developed as a consequence of carbonate mineral dissolution and pore scale transport properties. The results show that initial heterogeneity and porosity contrast control the development of the dissolution fronts in these highly reactive systems. This finding is consistent with linear stability analysis and the known positive feedback between mineral dissolution and fluid flow in carbonate formations. Differences in the carbonate kinetic drivers resulting from the range of pCO2 used in the experiments and the different proportions of more reactive calcite and less reactive dolomite contributed to the development of new pore space, but not to the type of dissolution fronts observed for the two different rock types. The development of the dissolution front was much more dependent on the physical heterogeneity of the carbonate rock. The observed stable dissolution fronts with small but visible dissolution fingers were a consequence of the clustering of a small percentage of larger pores in an otherwise homogeneous Marly dolostone. The observed wormholes in the heterogeneous Vuggy limestone initiated and developed in areas of greater porosity and permeability contrast, following pre-existing preferential flow paths. Model calibration of core flood experiments is one way to specifically constrain parameter input used for specific sites for larger scale simulations. Calibration of the governing rate equations and constants for Vuggy limestones showed that dissolution rate constants reasonably agree with published values. However the calcite dissolution rate constants fitted to the Marly dolostone experiments are much lower than those suggested by literature. The differences in fitted calcite rate constants between the two rock types reflect uncertainty associated with measured reactive surface area and appropriately scaling heterogeneous distribution of less abundant reactive minerals. Calibration of the power-law based porosity–permeability equations was sensitive to the overall heterogeneity of the cores. Stable dissolution fronts of the more homogeneous Marly dolostone could be fit with the exponent n = 3 consistent with the traditional Kozeny–Carman equation developed for porous sandstones. More impermeable and heterogeneous cores required larger n values (n = 6–8). Ó 2013 Elsevier Ltd. All rights reserved.

1. Introduction Chemical interactions are known to have a major effect on rock porosity and permeability evolution and may alter the behavior or ⇑ Corresponding author. Tel.: +1 (925) 422 9657. E-mail address: [email protected] (Y. Hao). 0309-1708/$ - see front matter Ó 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.advwatres.2013.09.009

performance of natural and engineered reservoir systems. Such reaction-induced permeability evolution is of particular importance for the optimization of acid stimulation and CO2-enhanced oil recovery in carbonate reservoirs and has drawn considerable attention over the past decades (e.g., [1–27]). Predicting reactioninduced permeability evolution is also very important for geologic CO2 sequestration technologies designed to mitigate greenhouse

Y. Hao et al. / Advances in Water Resources 62 (2013) 388–408

gas emissions, because changes in rock porosity and permeability fundamentally affect the amount of CO2 that can be ultimately stored in the reservoir. This is particularly true for the transition of CO2-EOR fields to permanent CO2 storage, because the use of alternating water and gas floods to extract hydrocarbons will change the reservoir storage capacity. CO2 dissolution into reservoir fluids during EOR flooding mildly acidifies the resident brines and causes carbonate minerals to dissolve sometimes resulting in the formation of highly porous and conductive channels (often referred to in the petroleum acidizing literature as ‘‘wormholes,’’ e.g., [6]). Despite the large number of experimental and numerical studies devoted to understanding permeability evolution in carbonate rocks, only a few of them have made direct, comprehensive comparisons between predictive models and experimental data (e.g., [18–19]). The result is a lack of proper quantitative calibration, which limits the application of essential features of chemical fluid-rock interactions at laboratory or reservoir scales. Many research efforts (e.g., [7–11,13]) have employed the continuum modeling approach to explore key factors that affect mineral dissolution and related reactive transport processes, particularly focusing on investigation of the formation and development of unstable dissolution fronts or wormholes. Golfier et al. [7–9] presented a detailed study of dissolution channel development during acid dissolution of a porous medium. These authors outlined a three-dimensional Darcy-scale dissolution model that was derived from the pore-scale Stokes equation and volume-averaging theory, and summarized the governing equations and upscaling techniques used to determine macroscopic properties such as mass transfer coefficients, and dispersion properties. They further used the derived continuum model to examine the effects of advection, diffusion and reaction on wormhole formation, and developed parameter diagrams to describe different dissolution regimes and the transitions between them. More importantly, their results confirmed that the Darcy-scale continuum model was able to capture the general trends and patterns of wormhole formation. Panga et al. [10] have developed a more rigorous two-scale continuum model for simulations of wormholes in carbonate rocks by averaging pore scale physical and chemical behaviors into a continuum formation, so as to maintain the scale dependence of carbonate dissolution and transport processes. This model has also been used to analyze the influences of physical factors such as dispersion, heterogeneity magnitude, acid level, and pore-scale mass transfer on wormhole creation and pore volume required for breakthrough [10–11]. The continuum models that were used for simulation and analysis of carbonate dissolution processes are capable of predicting dynamic evolution of rather complex dissolution structures and wormhole geometry. However in most of these models mineral dissolution is often simplified as a one-component, first-order kinetics reaction that cannot provide a viable account of geochemical evolution due to mineral alteration. Furthermore, one of the largest challenges or uncertainties facing continuum-scale reactive-transport models is whether or not the so-called macroscopic properties (e.g., chemical rate kinetics, porosity–permeability relationships, reactive surface area expressions) can effectively describe the physical and chemical behaviors exhibited at the pore scale (e.g., [28–29]). For example the direct use of chemical rate kinetics measured from laboratory-based batch experiments to describe mineral alteration and associated reactive transport processes in continuum models, although quite common in the literature, has been questioned by recent researchers. Szecsody et al. [30] have found that directly implementing batch experiment-based chemical kinetics failed to reproduce reactive transport processes observed in column experiments. The fact that dissolution rates are several of orders of magnitude lower when measured in the field than the laboratory reflects a potential scale

389

dependence of reaction kinetics [31]. A number of numerical studies (e.g., [32–40]) have also recognized the scaling effects of macroscopic transport and geochemical reaction properties by comparison of continuum and pore-scale models, and attributed such scale dependence to pore-scale heterogeneities. Specifically, Li et al. [32] have developed a pore-scale network model to study the scale dependence of geochemical reaction kinetics controlling mineral dissolution and precipitation. Comparison of reaction rates averaged from pore-scale modeling results with those obtained by the continuum model using laboratory-measured chemical kinetics showed that scaling effects became important as pore-scale heterogeneity increased or acid levels were elevated within porous formations. They also used the same approach to examine the dependence of reaction rates on heterogeneous distribution of reactive minerals at a pore scale, finding larger scaling effects of reaction rates for a porous medium either with small mineral abundances or with large clusters of reactive minerals [33]. All of these findings emphasize the importance of calibrating macroscopic parameters used in the continuum models by either porescale simulations or experimental data in order to effectively account for scaling effects. In this study we applied a three-dimensional continuum reactive transport model to simulate CO2 core flood experiments performed by Smith et al. [27]. The numerical results were integrated with experimental data to analyze and interpret CO2 injection-induced carbonate dissolution phenomena from the pore to core scale. Specifically, we developed a comprehensive reactive transport modeling framework to interpret experimental results, and study the interplay between carbonate reactivity, flow velocity, effective permeability, and heterogeneity on the development of stable and unstable dissolution fronts. An empirical approach was developed and combined with advanced imaging techniques to capture the major pore structures and connectivity with a continuum representation, thereby leading to an improved understanding of effects of multi-scale heterogeneity on mineral alteration processes. This work adds to a growing body of experimental and numerical studies directed at understanding alteration processes for CO2 injection and storage in carbonate reservoirs (e.g., [18–21,24–26]). Our work differs from other studies in that the numerical model was constrained against detailed characterization data, and time-dependent solution chemistry and differential pressure data for two types of relatively low permeability rocks containing varying amounts of dolomite, calcite, and pore space. The rock formations studied here are specific to the reservoir flow units at the IEAGHG Weyburn-Midale CO2 Monitoring and Storage Project in south-eastern Saskatchewan, Canada (http://ptrc.ca/projects/weyburn-midale; [26]).

2. Methods 2.1. Overview of laboratory core-flood experiments The experiments that form the basis of the three-dimensional numerical reactive-transport modeling study are discussed briefly here. All references to experimental data and procedures are from Smith et al. [27]. Eight individual carbonate core samples collected from the Weyburn-Midale (Canada) CO2-enhanced oil recovery and storage site were tested to study the relationship between fluid flow, heterogeneity, and reaction by combining characterization, pressure/permeability, and solution chemistry data. These cores are from the Vuggy limestone V-6 unit (higher calcite abundance; lower initial porosity and permeability) or Marly dolostone M-3 unit (higher dolomite abundance; higher porosity and permeability). All experiments were conducted at 60 °C under constant 0.05 mL/min flow and 12.4 MPa fluid outlet pressure to prevent

390

Y. Hao et al. / Advances in Water Resources 62 (2013) 388–408

CO2 degassing. The Injected fluids contained pCO2 levels of 3, 2, 1, and 0.5 MPa and are representative of those calculated from site monitoring data. CO2-equilibrated brine was injected into a 15 mm-diameter, 30 mm-long core sample, and fluid samples were collected at the outlet and analyzed for calcium, magnesium, sulfate, total inorganic carbon (TIC), sodium, and chloride. Solution pH and mineral saturation indices were modeled with the EQ3/6 geochemical modeling tool [41] and the data.shv database [42]. Mineralogy and pore space distributions within each sample were characterized by X-ray computed microtomography (XCMT) prior to and after reaction at a final resolution of 14.8 lm/pixel. XCMT and scanning electron microscopy (SEM) analyses were combined to reconstruct three-dimensional pore space, mineral distribution, and initial accessible surface area to quantify the changes induced by carbonate dissolution. The initial sample porosity, permeability, and mineralogical data, and initial and injected brine composition are detailed in [27], and also reported in Tables 1 and 2, respectively. The total pore volumes for initial Marly and Vuggy samples are 1.6–1.9  106 and 7–8  107 m3. Under a constant 0.05 mL/ min injection flow the mean residence time of solutes ranged from 0.55 to 0.65 h for Marly cores, and from 0.23 to 0.28 h for Vyggy samples. We adopt the same naming convention as in the experiments to refer to specific core types and fluid pCO2 (e.g., ‘‘M-2’’ indicates a Marly dolostone core subjected to brine equilibrated with 2 MPa pCO2, ‘‘V-2’’ indicates a Vuggy limestone core subjected to brine equilibrated with 2 MPa pCO2). 2.2. Reactive transport modeling of core flood experiments 2.2.1. Mathematical formulation of reactive transport in porous media Coupled single-phase flow and multi-component reactive transport processes in porous media are mathematically constrained by the principle of mass conservation. Following Darcy’s law the transport of a single fluid phase in porous media is governed by

l

ðrp  qgÞ

ð2Þ

in which p is pressure, and K, l and g denote the intrinsic permeability, phase viscosity, and gravitational acceleration vector. In a multicomponent flow system with chemical reactions we assume homogeneous equilibrium in the aqueous phase and kinetically controlled mineral reactions, which are expressed as

Aj

X

Value

Pre-CO2 brine pH Brine composition [mole kg H2O1]

6.9 ± 0.2 [NaCl] = 1.01 [HCO3] = 0.00792 [SO4] = 0.0369 [CaCl2] = 0.0353 [MgCl2] = 0.0159 4.3, 0.33 m 4.5, 0.22 m 4.8, 0.11 m 4.9, 0.055 m

3 MPa pCO2 brine pH, [CO2(aq)] 2 MPa pCO2 brine pH, [CO2(aq)] 1 MPa pCO2 brine pH, [CO2(aq)] 0.5 MPa pCO2 brine pH, [CO2(aq)]

Am

X

mim Ai ;

mij Ai

ð3Þ

i

and

Table 1 Average mineral and pore space abundances (determined from pre-reaction tomography and SEM analyses) and bulk permeability (determined from differential fluid pressures recorded prior to CO2 introduction) [27]. Sample ID

Calcite (%)

Dolomite (%)

Porosity (%)

Permeability (mD)

M-0.5 core M-1 core M-2 core M-3 core V-0.5 core V-1 core V-2 core V-3 core

14 ± 9 20 ± 17 7.7 ± 4 19 ± 7 56 63 55 59

52 ± 7 49 ± 11 61 ± 4 45 ± 4 29 24 27 26

35 ± 5 31 ± 8 31 ± 3 36 ± 4 15 13 16 15

1.7 0.92 1.2 1.7 0.0091 0.024 0.032 0.032

ð4Þ

i

where Ai, Aj, and Am indicate primary and secondary aqueous species, and mineral phase, and mij and mim refer to stoichiometric coefficients of the ith primary species in jth aqueous and mth mineral reactions, respectively. Multicomponent transport processes are constrained by mass balance as,

X @ ð/cti Þ ¼ r  ðVcti Þ þ r  ð/Drcti Þ þ Ci  mim Rm : @t m

ð5Þ

While the left hand side of the equation denotes the rate of change in total mass of the ith primary species the terms on the right hand side account for contributions from advection, diffusion/dispersion, sources and sinks, and mineral dissolution/precipitation kinetic reactions. D, Ci and Rm represent combined diffusion and dispersion coefficients, source/sink terms of the ith species, and reaction rate for the mth mineral reaction, respectively. In this formulation cti is defined as the total concentration of the ith primary species,

ð1Þ

Here t, r, and C respectively represent time, spatial derivative operator, and source/sink terms. / is the porosity, q is the fluid density, and the Darcy flux vector V is expressed as

V¼

Experimental Parameter

cti ¼ ci þ

@ ð/qÞ ¼ r  ðqVÞ þ C: @t

K

Table 2 Initial and injected brine composition [27].

X

mij cj

ð6Þ

j

and the concentration for a secondary aqueous species is determined from the primary species concentrations through mass action equations,

Y cj ¼

mij

ðci ci Þ

i

K eq j cj

;

ð7Þ

in which ci and cj are the concentration of ith primary and jth secondary species, K eq j and cj are equilibrium constant of the aqueous reaction (Eq. (3)) and activity coefficient, respectively. Activity coefficients of charged aqueous species are calculated with an extended Debye–Hückel model (B-dot formulation) [43], and the activity coefficient for aqueous CO2 is determined using the Drummond formulation [44]. The changes in mineral phase volume and porosity due to mineral dissolution and precipitation are calculated by

dhm ¼ Rm V m dt

ð8Þ

and

X dhm d/ ¼ ; dt dt m

ð9Þ

with hm and Vm denoting the mth mineral volume fraction and molar volume. The permeability evolution resulting from porosity change is evaluated from a pre-defined porosity–permeability relationship, which will be discussed in a later section. A single continuum model was used to predict CO2-induced carbonate dissolution in Marly dolostone samples. We extended this

Y. Hao et al. / Advances in Water Resources 62 (2013) 388–408

Fig. 1. Dual-continuum (dual porosity/permeability) representation of the Vuggy limestones.

formulation to dual continuum (dual porosity/permeability) model for the case of more heterogeneous but less permeable Vuggy limestones. As will be further discussed later in Section 2.2.6, in this dual continuum model, the porous medium domain is conceptualized as two interacting continua in overlapping model grid regimes, one representing the interconnected macropore regions (e.g., fractures, vugs) and the other the low permeability rock matrix (Fig. 1). The mass balance equations for two continua, which are coupled through inter-continuum mass transfer functions, are separately defined in the form of Eqs. (1)–(9).

2.2.2. Numerical simulation tool The core flood experiments were simulated using reactive transport models that combine advective–diffusive transport with both aqueous chemical equilibrium and kinetic reactions. The reactive-transport simulations were performed by using the Nonisothermal Unsaturated–Saturated Flow and Transport code, NUFT [45–46]. The NUFT code, which is based on a Darcy’s flow approximation, represents a state-of-the-art simulator for modeling multiphase, multi-component heat and mass flow and reactive transport in unsaturated or saturated porous media. An integrated finite difference method was used for numerical discretization. The NUFT code is a highly flexible computer software package capable of running on PCs, workstations, and major parallel processing platforms. Some of its previous applications include: geothermal exploitation, nuclear waste disposal, CO2 sequestration, groundwater remediation, and subsurface hydrocarbon production (e.g., [47– 54]).

391

2.2.3. Model set-up Each core sample of 15 mm in diameter and 30 mm in length was discretized using a finite number of numerical grid blocks for the three-dimensional model, as shown in Fig. 2. A nominal grid size of 500  500  500 lm3 was used as the representative element volume, resulting in total 4.1  104–6.4  104 numerical grid blocks. The initial and boundary conditions in the numerical model mimicked those imposed experimentally, with impermeable lateral boundaries and constant pressure and flux conditions imposed at the outlet and inlet boundaries, respectively. The liquid phase density and viscosity were taken as 988.5 kg/m3 and 4.6  104 Pa s, respectively. The initial brine chemistry was calculated by using EQ3/6 [41]. The thermo-dynamic data were generated from the SUPCRT (data.shv) database [42]. Transport of chemical species is controlled by advection and diffusion, with a uniform diffusion coefficient of 109 m2/s for all species and no dispersion. Sensitivity analyses indicated minimal impacts of variation of diffusion and dispersion on chemical transport and mineral reactions. 2.2.4. Chemical equilibrium and reaction kinetics All aqueous reactions were assumed to proceed at equilibrium 2+ 2+ + with the primary species selected as H+, HCO 3 , Ca , Mg , Na , 2  Cl , and SO4 . Calcite and dolomite dissolution and precipitation reactions were modeled by using the following kinetic rate law that includes both acid and neutral mechanisms,

   EN EA 1 1 Qm m 1 m 1 N A ; Rm ¼ Sm km e R ðT 298:15Þ þ km e R ðT 298:15Þ anHþ 1  eq Km

ð10Þ

where Rm is the mineral reaction rate, Sm is the specific surface area for the mineral, T is the absolute temperature, R is the gas constant, K eq m is the equilibrium constant of mineral reaction and Qm is the corresponding activity product. aH+ and n denote the activity of N A H+ and reaction order, and km , ENm and km , EAm are the rate constants at 298.15 K and activation energies for neutral and acid mechanisms, respectively. Here Q = K implies equilibrium; Q < K implies undersaturation (and dissolution); and Q > K implies supersaturation (and precipitation). The magnitude and direction of the rate is determined by the relationship between Q and K. A negative value of Rm indicates dissolution reactions, whereas a positive value means precipitation reactions. The rate constants were initially based on neutral and acid mechanisms reported by Palandri and Kharaka [55], and fitted values were refined in the simulations to

Fig. 2. Model set-up. (a) core sample; (b) numerical discretization; (c) numerical model domain along with corresponding boundary conditions .

392

Y. Hao et al. / Advances in Water Resources 62 (2013) 388–408

Table 3 Major equilibrium and kinetic reactions considered in the simulations. log Keqa 60 °C

Reactions

Acid mechanism Em

+

2+

+ HCO 3 2+ 2+

Calcite + H = Ca Dolomite + 2H+ = Ca + Mg + 2 HCO 3 CO2(aq) + H2O = H+ + HCO 3 MgHCO3+ = Mg2+ + HCO 3 CaCO3(aq) + H+ = Ca2+ + HCO 3 þ  2+ CaHCO3 = Ca + HCO3  + MgCO3(aq) + H = HCO3 + Mg2+ MgCl+ = Cl + Mg2+ CaCl+ = Cl + Ca2+ CaCl2(aq) = 2Cl + Ca2+ CaSO4(aq) = SO42 + Ca2+ NaCl(aq) = Na+ + Cl a b c d e f

1.4011 1.4558 6.2059 1.1564 6.4028 1.1515 6.8816 0.0695 0.603 0.6674 2.2364 0.6693

b

14.4 36.1

n

c

1.0 0.5

Neutral mechanism acid

log k

d

0.30 3.19d

e

1.87 4.19e

f

0.43 4.19f

E mb

log kneutral

23.5 52.2

5.81d 7.53d

7.26e 7.53e

5.81f 7.53f

EQ3/6 [41] and SUPCRT (data.shv) database [42]. Activation energy (kJ/mol) from Palandri and Kharaka [55]. Reaction order with respect to H+ [55]. Reaction rate constants (mol/m2 s) at a temperature of 25 °C [55]. Reaction rate constants (mol/m2 s) at a temperature of 25 °C used for Marly dolostones. Reaction rate constants (mol/m2 s) at a temperature of 25 °C used for Vuggy limestones.

best match the experimental and characterization data. We chose to use the acid and neutral rate mechanism to capture changes in acidity resulting from variable pCO2 in the experiments and in carbon storage sites. Although it is possible to fit the data with a pH-independent constant [26], the fitting of the experiments to both the acid and neutral terms might allow the resulting equations to be extended to applications with a larger variation in pH such as acid stimulation operations. The major geochemical reactions and relevant kinetic rate constants and activation energies for calcite and dolomite reactions used in this study are listed in Table 3. Dolomite precipitation was explicitly disallowed in the model to reflect known difficulties for its precipitation under laboratory conditions [56–57]. 2.2.5. Porosity–permeability-surface area relationships The macroscopic porosity–permeability correlation is an important parameter for modeling flow and reactive transport processes at the continuum scale, because mineral reactions change transport properties of the rock. In this work, we employed a power-law equation to describe the dependence of permeability on porosity change due to mineral alteration (based on [58]):

Kt ¼ K0



/t /0

2.2.6. Initial core petrophysical properties It is of great importance to define and capture multi-scale heterogeneity observed in the carbonate rocks in the continuum model, because heterogeneity can influence chemical transport and mineral reaction [7,10,62–64]. Additionally, the representative elementary volume that forms the model grid blocks (REV = 500  500  500 lm3) was significantly larger than tomography voxel resolution (14.83 lm3), where each modeling grid block contained as many as 3.8  104 imaging voxels. In this section we describe the methods we used to upscale the high-resolution characterization data and estimate the effective grid-based porosity and permeability values for the continuum description.

2.2.6.1. Effective porosity and mineral phase composition. Several representative mixtures of pore space, calcite and dolomite abundance, referred to here as ‘‘matrices’’, were defined for specific gray-scale ranges from the XCMT characterization data of the

n ;

ð11Þ

where / and K are porosity and permeability; 0 and t indicate initial and evolving time values; and the exponent n is an empirical power index derived here from the match between simulations and experimental data. Although we recognize potential effects of changing mineral volumes on reactive surface area, our experimental datasets were not designed to provide the necessary parameters for use with the ‘‘sugar lump’’ model described by Noiriel et al. [59]. For this reason we used a simple representation of reactive surface changing in proportion to mineral phase abundances (based on the models of Emmanuel and Berkowitz [60] and Lichtner [61], also evaluated in Noiriel et al. [59]):

St ¼ S0

 2=3  2=3 ht /t ; h0 /0

ð12Þ

where S and h are the specific surface area and mineral volume fraction, and 0 and t indicate initial and evolving time values due to mineral alteration, respectively.

Fig. 3. Comparison of pore distribution in V-3 limestone prior to reaction between (a) XCMT image, and (b) porosity distribution in the numerical model.

Y. Hao et al. / Advances in Water Resources 62 (2013) 388–408

393

Fig. 4. Comparison of pore and mineral distributions in M-1 core prior to reaction between the characterization and model data. (a) XCMT image, and (b) porosity, (c) calcite, and (d) dolomite distributions in the numerical model.

unreacted cores. Higher-resolution SEM images of each matrix type were used to statistically estimate pore space, mineral composition and accessible surface area [27]. We employed a volumetric averaging approach on the calibrated XCMT gray-scale matrices to calculate effective porosity and mineral volume fraction for each numerical grid block as,

PN /¼

i

/i

N

ð13Þ

and

PN hm ¼

hm;i ; N

i

ð14Þ

where N is the total number of voxels contained in the grid block, /i and hm,i represents the microporosity and mineral volume fraction of voxel i, respectively.

As mentioned previously we modeled the Vuggy limestones as a heterogeneous dual-continuum system with different porous medium properties (e.g., porosity, mineral composition, and permeability). The representative elementary volume was further divided into two separate volumes (continua) based on their permeability characteristics (described in Section 2.2.6.2). Porosity, and calcite and dolomite volume fractions for each continuum were averaged in a similar way as expressed by Eqs. (13) and (14). Figs. 3 and 4 compare pore and mineral distributions between the characterization data and the continuum model representation for the V-3 limestone and M-1 dolostone cores. The spatial distribution of pore space, calcite, and dolomite imaged in the Vuggy limestone and Marly dolostone were well preserved at the continuum grid level. A comparison of the porosity cumulative distribution functions between model and XCMT data (Fig. 5) further illustrate that the adopted volumetric averaging approach appears

Fig. 5. Comparison of porosity cumulative distribution functions between XCMT tomography and model data for (a) M-0.5, (b) M-1, (c) M-2, and (d) M-3 dolostones, and (e) V-0.5, (f) V-1, (g) V-2, and (h) V-3 limestones.

394

Y. Hao et al. / Advances in Water Resources 62 (2013) 388–408

to sufficiently incorporate primary rock characteristics and heterogeneity into the numerical model.

2.2.6.2. Permeability distribution. Scaling the detailed characterization data into effective permeability tensors at the grid basis was not as intuitive and straightforward as the calculation of porosity and mineral volume fraction. The spatial variability of permeability across the samples is strongly enhanced in these carbonate units, particularly the Vuggy limestones, by pore space heterogeneity observed at multiple length scales. It was not possible to use a porescale model to calculate the initial permeability, because pore connectivity analysis showed that large vugs and macropores (voids occupying more than 10 voxels) remained largely disconnected throughout the cores, implying that macropores were connected, if at all, by micropore structures (of the size between 0.1 and 1 lm) at a scale far below XCMT resolution (see Fig. 6). The pore space connectivity was determined by using a ‘‘burning’’ algorithm [65–66]. In order to overcome these difficulties we developed an empirical approach that combined XCMT and SEM image analysis with experimental data to infer the initial permeability field in the samples. The approach involves several steps as follows: 2.2.6.2.1. Step 1: qualitative assessment of pore structures and their permeability. We categorized the representative pore structure patterns observed in the 2-D segmented SEM images for XCMT gray-scale matrix types into high, medium, and low permeability

regimes, each assumed with a constant permeability value, as shown in Fig. 7a. In what follows we use Khigh, Kmed and Klow to denote the assumed high, medium, and low permeability constants. This approach allows the local permeability field within each numerical grid block to be defined and provides the basis to scale the permeability from the pore to continuum scale. 2.2.6.2.2. Step 2: grid-based permeability calculation. We used different methods to estimate their equivalent permeability values on a grid-block basis for the Vuggy and the Marly cores because their pore space connectivity and heterogeneity differ significantly (Fig. 6). The Vuggy samples contained many medium- and large-size vugs that had a limited impact on flow and transport behavior because most of them were isolated and disconnected from each other (Fig. 6(e)–(h)). However, some distinct macropore clusters or fractures spanned distances of several continuum grid cells. These pore structures may have served as potential preferential paths linking large pores or vugs, and thereby transferring most of the flow and controlling the permeability. These observations suggest that permeability distribution might be more influenced by pore connectivity at a broad range of length scales. For this reason we evaluated the connectivity of different types of pore structures identified in step 1, especially distinguishing between high, medium, and low permeability regions throughout the core. We treated both high and medium permeability regions that formed

Fig. 6. Macro-pore connectivity in (a) M-0.5; (b) M-1; (c) M-2; (d) M-3 dolostones, and (e) V-0.5; (f) V-1; (g) V-2; (h) V-3 limestones.

Y. Hao et al. / Advances in Water Resources 62 (2013) 388–408

globally-connected flow pathways as one continuum, and low permeable (or impermeable) rock matrix as the other continuum to create a dual-continuum model. The low permeability value, Klow, was directly assigned to less permeable matrix continuum. However, it was a challenge to determine grid-based effective permeability tensor distribution for more permeable continuum because high and medium permeability zones were closely clustered or overlapping with one another, leading to complex heterogeneity and strong permeability anisotropy. In order to preserve geometric features of connected pore regions and their effects on permeability distribution we mapped clusters of vugs, fractures, and highly permeable regions observed in tomography data onto a three-dimensional regularly spaced grid (Fig. 8(a)) in a manner analogous to mapping discrete fractures onto a continuum system, as was proposed by Botros et al. [67] and Reeves et al. [68]. Then we evaluated the effective permeability tensors on a grid block basis, following an approach similar to that used in fault or fracture transmissibility calculation for reservoir simulation models [67– 69], which is summarized below and also depicted in Fig. 8(b). As described in [67–69] the volumetric Darcy flux Qi,j between adjacent grid blocks i and j in three dimensional space is written as

Q ij ¼

K ij

l

Aij

Dpij ; L

ð15Þ

where Dpij, Kij, Aij, and L are the pressure difference, equivalent inter-grid permeability, cross-sectional area and distance between

395

two grid blocks, respectively. The volumetric Darcy flux Qi,j can be further expressed as the sum of flow contributions from high and medium permeable regions, Q hij and Q m ij ,

Q ij ¼ Q hij þ Q m ij ¼

 Dpij  K high Ahij þ K med Am ij ; lL

ð16Þ

where Ahij and Am ij represent the total cross-sectional areas of high and medium permeability zones, respectively, which pass through the interface between grid blocks i and j. We estimated them geometrically by mapping inter-connected permeable pathways onto a regularly spaced grid system. Finally combining Eq. (15) and (16) yields equivalent permeability for i  j grid connection,

K ij ¼

 1  K high Ahij þ K med Am ij : Aij

ð17Þ

A more uniform distribution of macro-pore connections was observed in Marly samples, highlighting the homogeneity of these samples. Therefore we employed a simplified renormalization technique [70] to transform local gray-scale defined permeability distribution (Fig. 7(b)) into three-dimensional permeability tensors on a grid-block basis. 2.2.6.2.3. Step 3: experimental calibration. The grid-based effective permeability values were further calibrated by using a global permeability multiplier to match numerical to experimental initial pressure differences across cores. The above steps were repeated by adjusting initial gray-scale specific permeability values (e.g.,

Fig. 7. Estimation of permeability due to microporosity for each XCMT gray scale from 2D segmented SEM images. (a) Three permeability regimes were identified based on qualitative assessment of pore structures and connectivity from SEM data, each defined over a specific XCMT gray-scale range. Khigh, Kmed, and Klow denote constant permeability values assumed for high, medium and low permeable regions, respectively. Note that in the segmented SEM images green is pore space, red is dolomite, and black is calcite. (b) Estimated permeability distribution across a sample cross-section (left image), and inside a continuum grid (right image).(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

396

Y. Hao et al. / Advances in Water Resources 62 (2013) 388–408

Khigh, Kmed and Klow) in step 1 until the corresponding simulated carbonate dissolution fronts, transient responses of pressure and chemistry solutions fit both the experimental measurements and characterization data. 3. Results 3.1. Marly dolostone: development of stable dissolution fronts Figs. 9–12 compare the measured and simulated final dissolution fronts, time-dependent output solution chemistry, and mass transfer rates for Marly dolostones upon reaction with pCO2 = 0.5, 1, 2, and 3 MPa brines. The mass transfer rates for dolomite and calcite are defined in [27]. The solution chemistry provided a key constraint to the model because it was a measure of the mass transfer from solid to solution that was responsible for new pore space within the rock. Solution pH is a consequence of the input brine pCO2 as well as carbonate mineral dissolution, which partially neutralizes the initially acidic pH levels. Both model and experiment demonstrated the ability of calcite and dolomite dissolution to rapidly buffer the solution composition. It is noted that the simulation under-predicted the final pH for the M-0.5 case by about 0.25 units when compared to the experimental results. This deviation may be attributed to the fact that the measured total inorganic carbon (TIC) data was lower than expected pCO2 level for this experiment [27]. The agreement between modeled and experimental solution chemistry translated to similarly reasonable agreement for estimates of calcite and dolomite mass transfer rates, in particular for experiments conducted at higher pCO2 (Fig. 12). The three-dimensional numerical model results captured the general trend of the carbonate dissolution fronts observed via XCMT imaging. The combined use of a cubic porosity–permeability

relationship (n = 3, Eq. (11)) and effective mapping of the initial permeability field reproduced the alteration patterns observed in three of the four reacted dolostone cores (M3-0.5, 2, 3), as shown in Figs. 9(a), 11(a) and 12(a). These patterns consisted of relatively stable reaction fronts with fairly evenly distributed pore space, calcite and dolomite distributions, especially along the axial direction of each core sample, and small but visible dissolution fingers formed downstream of the stable dissolution fronts. The model reproduced the flow path developed within sample M-1 (Fig. 10(a)), which remained preferentially confined to the more porous dolomitic half of the core (see Fig. 4) with a higher porosity–permeability correlation of n = 8 (Eq. (11)). 3.2. Vuggy limestones: unstable dissolution fronts Figs. 13–16 compare the measured and simulated final dissolution fronts, and time-dependent differential core pressure, output solution chemistry, and mass transfer rates for the relatively impermeable Vuggy limestone cores upon reaction with pCO2 = 0.5, 1, 2, and 3 MPa brines. The three-dimensional numerical model predicted the formation of highly porous and conductive channels spanning inlet to outlet of the core, as well as localized branching channels. The wormholes are largely attributed to high heterogeneity and porosity/permeability contrast of these samples, as will be discussed below. The pressure change across these cores was a useful constraint for the form of the porosity–permeability function used in the model. The Vuggy cores required a higher-order correlation between permeability and porosity (n = 6–8) than most porous Marly dolostone cores (n = 3), as well as accurate mapping of effective permeability at the grid scale, in order to match the differential pressure decreases and the final shape and orientation of the wormholes. Modeled pressure curves matched the observed pressure data

Fig. 8. Visual representation of the empirical method used to scale continuum grid permeability for Vuggy limestones by (a) performing pore connectivity analysis (denoted by arrow 1), and mapping inter-connected macropore clusters onto a regularly spaced continuum grid system (denoted by arrow 2) with different colors representing different connected pore clusters, and (b) estimating inter-grid effective permeability.

Y. Hao et al. / Advances in Water Resources 62 (2013) 388–408

for all the Vuggy cores. The final wormhole shapes and orientation predicted by the simulations are largely consistent with those observed in XCMT data, with the exception of the V-2 case, in which the simulated wormholes developed along different preferential paths to observed paths. A possible explanation of this discrepancy between model results and tomography data is that the wormhole development was triggered by small-scale heterogeneities that were below the resolution of the tomography data, and thus were not properly represented by the initial permeability estimation. This also highlights the strong dependence of model results on the initial core characterization. The model results for solution chemistry were in good agreement with experimental measurements, in particular for times prior to pressure breakthrough, except for the V-0.5 case where the simulation overestimated solution acidity and underestimated calcium concentrations. As discussed by Smith et al. [27] the discrepancy observed for the case of V-0.5 case may be attributed with difficulties controlling very low pCO2 conditions during the experiments. Both experimental and simulation results confirmed that pressure breakthrough (related to the development of a lowresistance, through-going dissolution feature) occurred between 25–50 h after CO2 injection and was accompanied by a drop in pH except the V-0.5 case. The slight differences between simulated and experimental post-breakthrough pH (Fig. 16(c)) suggest more carbonate mineral dissolution where the simulated pH is higher and less carbonate mineral dissolution where the measured pH is lower.

397

4. Discussion In this section we discuss the interplay between heterogeneity, carbonate reactivity, flow velocity, and time on the development and consequences of stable and unstable dissolution fronts. We also discuss the limitations of applying continuum-scale models to the evolution of porosity and permeability in carbonate rocks. 4.1. The role of heterogeneity on the initiation of dissolution front The carbonate rocks used in this study exhibited heterogeneity at multiple length scales and allowed common theories that relate flow, mineral reactivity, and porosity/permeability contrast on the development of stable and unstable dissolution fronts to be tested. Previous work has correlated the shape of dissolution fronts to the heterogeneity, and as well as relative dominance of advection, diffusion, and/or reactivity using Péclet (Pe) and Damköhler (Da) dimensionless numbers (e.g., [1,7,10]). The theoretical work of Szymczak and Ladd [64] bears distinction because it addressed the importance of the porosity contrast across the reactive front in promoting the spontaneous development of an unstable dissolution from the inlet. These authors proposed the mechanism for the initiation and development of unstable dissolution fronts in highly porous and permeable materials (as is the case of the Marly dolostone), as well as within less permeable formations with small initial porosity (e.g., Vuggy limestones). Our experimental and simulation results confirmed the important roles of porosity

Fig. 9. Comparison between the measured and simulated (a) final shape of dissolution fronts, time-dependent output solution (b) pH and (c) species concentration, and (d) mass transfer rates for M-0.5 case. The right plot in (a) shows the predicted porosity change with the iso-surface threshold value as 0.03.

398

Y. Hao et al. / Advances in Water Resources 62 (2013) 388–408

contrast and sample heterogeneity on the development of dissolution front types. In our experiments, the range of mineral reactivity produced by varying pCO2 from 0.5 to 3 MPa did not sufficiently alter the reaction kinetics for a given rock type, and the type of dissolution front was not strongly associated with variations in pCO2. The dissolution front was associated with rock type, sample heterogeneity and porosity/permeability contrast as proposed by Ortoleva et al. [63] and Szymczak and Ladd [64]. In this section we use the simulation results to discuss the effects of heterogeneity, porosity contrast, and pore flow velocities on the development of dissolution fronts. The combined observation of ‘‘stable’’ dissolution fronts with small but visible dissolution fingers in three of the four Marly dolostone samples (M-0.5, M-2, and M-3) is consistent with the notion that instabilities can develop in the presence of seemingly small variations within an otherwise homogeneous matrix. The small dissolution fingers were caused by non-uniform reaction over the core cross-section, and can be understood as the interplay between flow and reaction at macro- and micro-pore levels. Micropores within Marly samples constituted more than 90% of pore space, shared a relatively narrow size distribution between 0.5 and 5 microns, and were uniformly distributed and well-connected, preserving the overall homogeneity of this samples formation. The macropores on average accounted for less than 10% of total pore space, varied in size from several tens of microns up to one millimeter, and were not fully connected throughout the core even at

SEM-scale resolution. The presence of macropores, particularly when clustered together to form more conductive pathways, introduced local-scale heterogeneity in otherwise homogeneous cores, leading to variations in the transport of reactive fluid and the initiation of dissolution ‘‘fingering’’. Such small-magnitude instabilities or fingering tended to grow slowly in more porous Marly formations that have small porosity contrast. A steadily advancing dissolution front was observed during the periods of most Marly experiments (M3-0.5, M-2, and M-3), and prior to the development of substantial instability or fingering in M3-1 core (Figs. 9a–12a, 17(a)–(c)). These observations are in a qualitative agreement with predictions made by linear stability theory [64]. The positive feedback between mineral dissolution and localized permeability can eventually lead to the evolution of a prominent preferential flow path. We extended the simulation beyond the duration of the M-3 experiment to illustrate this phenomenon (Fig. 17(d)). Continued flow of CO2-rich brine led to the formation of a conically shaped wormhole that diverted fluid flow and locally increased pore velocities. The effect of increasing pore velocity on the development of dissolution front is also observed in the M-1 core (Fig. 10(a)), which possessed a bimodal distribution of impermeable calcite and more permeable dolomite as seen in Fig. 4. During reactive brine infiltration, the injected fluid largely bypassed the calcite-rich zone in spite of its high reactivity and flowed instead through the more porous dolomite-rich half of the core. Fluid flow through the more porous half reduced the effective cross-sectional

Fig. 10. Comparison between the measured and simulated (a) final shape of dissolution fronts, time-dependent output solution (b) pH and (c) species concentration, and (d) mass transfer rates for M-1 case. The right plot in (a) shows the predicted porosity change with the iso-surface threshold value as 0.05.

Y. Hao et al. / Advances in Water Resources 62 (2013) 388–408

399

Fig. 11. Comparison between the measured and simulated (a) final shape of dissolution fronts, time-dependent output solution (b) pH and (c) species concentration, and (d) mass transfer rates for M-2 case. The right plot in (a) shows the predicted porosity change with the iso-surface threshold value as 0.03.

area, and increased the local flow velocity and the local Pe number, consistent with the notion that higher flow velocity tends to lead to the development of preferential flow paths. Additionally, heterogeneous distributions of pore space and carbonate minerals in the dolomitic half resulted in uneven mineral dissolution fronts, illustrated by the different lengths of the conically shaped wormholes (Fig. 10(a)). The development of wormholes in the Vuggy samples is consistent with the notion that unstable dissolution fronts initiate instantaneously in carbonate rocks with high porosity contrast associated with the heterogeneous distribution of pore space [64]. Although micropores comprised more than 80% of pore space in the Vuggy cores, they possessed a much wider pore size distribution (between 0.05 and 100 lm) than in the Marly cores, and were not uniformly distributed throughout the full length of the cores. Vuggy heterogeneity was further accentuated by a similarly wide range of macropore sizes. The role of permeability/porosity contrast and heterogeneity on the development of wormholes in Vuggy limestones, as a consequence of reactive transport and carbonate mineral dissolution, is best illustrated in Fig. 17(e)–(h), showing the evolution of the wormhole over time. In this example (V-3), the initial macropore connection map (Fig. 6(h)) shows a few small fractures spanning the first few millimeters into this low-permeability core, forming narrow but highly conductive flow paths from the inlet. At early times, injected CO2-equilibrated brine was channeled more rapidly

into these fractured domains than into low-permeability matrices, leading to uneven flow and mineral dissolution across the core cross-section. As injection and dissolution progressed, the higher porosity/permeability contrast between reacted and unaltered zones focused more reactive fluid into these flow paths, confining the bulk flux through a smaller cross-sectional area. Branching dissolution fingers competed for reactive brine until a single channel broke through to the large vugs present in the downstream half of the core, essentially forming a preferential flow channel through the core. The simulated wormhole was completely formed after 25 h of injection, in agreement with the time necessary to affect a complete decrease in core differential pressure (see Fig. 16(b)). The simulated evolution of this wormhole clearly illustrates that unstable dissolution was initiated in zones with high pore connectivity and grew through the eventual dissolution of very low permeable matrices present between regions of higher permeability. 4.2. Carbonate reactivity Reactive transport models require rate equations that are tied to mineral surface area and mineral equilibrium to describe the evolution of pore space in highly reactive systems such as the carbonate core flood experiments described here. We used the simulations to fit experimental data assuming separate acid and neutral mechanisms. The rate constant values were derived from iteratively fitting the simulation to the solution chemistry, pressure,

400

Y. Hao et al. / Advances in Water Resources 62 (2013) 388–408

and tomography data, using the initial distribution of carbonate minerals and surface areas from the characterization data as well as activation energies reported in the literature as parameter input for forward calculations. The fitted results are compared with literature values at 25 °C in Table 3 and at 60 °C in Table 4. The fitted rate constants for dolomite were the same for both the Vuggy limestone and Marly dolostone experiments and have values of k(25 °C)acid, dol = 104.19 and k(25 °C)neutral, dol = 107.53 mol m2 s1. The fitted rate constants were about 8 times lower than those determined for single mineral systems in CO2-poor fluids [55] and 2 times lower than those derived in CO2-rich fluids at conditions far from equilibrium [71] (Table 4). The fitted rate constants for calcite for the Vuggy limestone experiments have values k(25 °C)acid, cal = 100.43 and k(25 °C)neutral, cal = 105.81 mol m2 s1. These values yield net dissolution rates that are within 6 times those reported in the literature when extrapolated to 60 °C (Table 4). This was not the case for the Marly dolostone experiments, which required acid and neutral constants that were 200 and 30 times smaller than their counter parts for the Vuggy experiments (Table 3). The large mismatch between fitted calcite rate constants for the Marly dolostone experiments and published data highlights some difficulties in capturing pore scale phenomena in the continuum model even when constrained by high resolution characterization and time dependent solution chemistry and pressure data. The mismatch may reflect an overestimation of effective surface areas

available for reaction under transport conditions, when averaging pore-scale heterogeneity from the XCMT characterization dataset onto the larger continuum grid scale, as has been found in other modeling studies (e.g., [36,72]). This effect appears to be amplified at lower mineral abundances. An independent modeling comparison study found that net reaction was dependent on heterogeneous distribution of reactive minerals at a pore scale, and that continuum models failed to capture pore-scale model results at lower abundances of the reactive mineral [33]. This explanation is consistent with our results which show good agreement for calcite rate constants derived from the core flood and single mineral dissolution experiments when calcite makes up 55–63% by volume in the Vuggy limestone and poor agreement when calcite makes up 8–20% calcite by volume in the Marly dolostone. The derived rate constant may also be sensitive to the form of the surface area equation (Eq. (12)), but we were unable to test more complex models (e.g., [59]). The simulated results suggest that acid and neutral rate constants coupled with a reaction affinity term can be used to predict mineral reactivity in CO2-EOR and CO2 storage environments and possibly acid stimulation operations, with some model calibration used to reduce the uncertainty associated with measuring reactive surface area and appropriately scale heterogeneous distribution of less abundant reactive minerals. In these cases, it is possible to fit a single rate constant that averages the effects of detailed mechanisms and still allows the rate to slow as equilibrium is approached

Fig. 12. Comparison between the measured and simulated (a) final shape of dissolution fronts, time-dependent output solution (b) pH and (c) species concentration, and (d) mass transfer rates for M-3 case. The right plot in (a) shows the predicted porosity change with the iso-surface threshold value as 0.03.

Y. Hao et al. / Advances in Water Resources 62 (2013) 388–408

and achieve similar fidelity between simulated and experimental solution chemistry [26]. 4.3. Effects of dissolution on permeability change Darcy-scale continuum models rely on porosity–permeability correlations to quantify the effects of dissolution on permeability change. We employed an empirical power law in the continuum models to relate local permeability change to porosity as described

401

by Eq. (11), where the exponent n was calibrated to fit the experimental data. The value of n appears to be strongly biased by sample anisotropy caused by initial pore connectivity (or lack thereof). All cores that exhibited unstable dissolution fronts were relatively more heterogeneous and higher values of n were required to fit the simulations to experimental data. The best-fit model value was n = 3 for M-0.5, M-2, and M-3; n = 6 for V-1 and V-2; n = 8 for M1, V-0.5, and V-3 experiments. These results indicate that higher values of n are indicative of a transition from a homogeneous to

Fig. 13. Comparison between the measured and simulated (a) final shape of dissolution fronts, (b) differential pressure values, time-dependent output solution (c) pH and (d) species concentration, and (e) mass transfer rates for V-0.5 case. The right plot in (a) shows the predicted porosity change with the iso-surface threshold value as 0.25.

402

Y. Hao et al. / Advances in Water Resources 62 (2013) 388–408

heterogeneous or unstable dissolution pattern as evolving porosity focuses the reaction pathway [19,24]. More heterogeneous rocks require higher n values, because less carbonate dissolution is required to connect the larger vugs or small fractures. The use of higher n values required for porosity–permeability correlations in the continuum models is also consistent with the trend observed experimentally by Smith et al. [27] in which the measured porosity

and permeability data were fitted to the power law function with exponent n ranging from 7 to 21. We also examined the simulated porosity and permeability values for all the cases, and found that at continuum grid blocks porosity evolved between 0 and 1, and permeability evolved between 5  1021 and 2  107 m2. These findings indicate that the simulated porosity or permeability values appeared to be in a reasonable range even for large n values.

Fig. 14. Comparison between the measured and simulated (a) final shape of dissolution fronts, (b) differential pressure values, time-dependent output solution (c) pH and (d) species concentration, and (e) mass transfer rates for V-1 case. The right plot in (a) shows the predicted porosity change with the iso-surface threshold value as 0.2.

Y. Hao et al. / Advances in Water Resources 62 (2013) 388–408

The experimental results of the M-0.5, M-2, and M-3 cases were well fitted with n = 3, implying that the traditional Kozeny–Carman correlation was adequate to describe permeability changes in these relatively homogeneous carbonate rocks. These Marly dolostone experiments were of limited duration because the uniform dissolution tended to compromise the integrity of the heatshrink core jacketing. One question posed by linear stability analysis is if preferential pathways would eventually develop given a longer experimental timescale. Extension of the simulation beyond

403

experiment run times allowed for qualitative assessment of the value of n and its potential impact beyond the empirical calibration we employed in this current work. Extending the simulation for the Marly dolostones yielded a conical shaped wormhole after about 5 days of reaction (Fig. 17(a)–(d)). Use of a higher n value (n = 8) still captured the early stages of dissolution front (Fig. 18(a)–(c); however the dissolution fingers were more pronounced and the eventual wormhole advanced more rapidly along the core walls (Fig. 18(d)).

Fig. 15. Comparison between the measured and simulated (a) final shape of dissolution fronts, (b) differential pressure values, time-dependent output solution (c) pH and (d) species concentration, and (e) mass transfer rates for V-2 case. The right plot in (a) shows the predicted porosity change with the iso-surface threshold value as 0.4.

404

Y. Hao et al. / Advances in Water Resources 62 (2013) 388–408

4.4. Dimensional analysis: dimensionless numbers As noted before the formation and development of unstable dissolution fronts are the result of the interplay between core heterogeneity, advection, diffusion, and carbonate reactivity. The relationship between these physical quantities can be further characterized and quantified by using the non-dimensional variables

and parameters such as Péclet (Pe) and Damköhler (Da), kinetic (Nk) and porosity contrast (D) numbers as was done in Smith et al. [27] prior to calibrating the reactive-transport model to fit the experimental results. The Péclet number, Pe, which measures the relative importance of advection to diffusion, is typically defined as

Pe ¼

uL : D

ð18Þ

Fig. 16. Comparison between the measured and simulated (a) final shape of dissolution fronts, (b) differential pressure values, time-dependent output solution (c) pH and (d) species concentration, and (e) mass transfer rates for V-3 case. The right plot in (a) shows the predicted porosity change with the iso-surface threshold value as 0.4.

405

Y. Hao et al. / Advances in Water Resources 62 (2013) 388–408

Fig. 17. Spatial and temporal development of wormhole and fingering in M-3 dolostones after CO2 flooding of (a) 4, (b) 8, (c) 12, and (d) 120 h; and in V-3 limestones after CO2 flooding of (e) 4, (f) 12, (g) 24, and (h) 40 h. The figures show the predicted porosity change with time.

Here u is an average flow velocity that is related to the Darcy velocity q by u = q//, D is diffusion coefficient, and L is a characteristic length scale, selected as either the core length or a characteristic pore-scale length. In order to describe the magnitude of carbonate reaction versus advective transport we define the Damköhler number by using the following formulation [27,28,73,74],

Da ¼

AkL ¼ /uC eq

     Acal kcal Adol kdol L wþ ð1  wÞ ; /u C eq;cal C eq;dol

ð19Þ

where k is a kinetic rate constant of calcite or dolomite reaction, A is reactive mineral specific surface area, Ceq is the solubility of each mineral, and w is a calcite-abundance weighting factor. The kinetic number also referred to as diffusive Damköhler number [27,28,74], is expressed as

Nk ¼ Pe  Da ¼

     2 Acal kcal Adol kdol L wþ ð1  wÞ C eq;cal C eq;dol /D

ð20Þ

and used to compare the effects of carbonate reactivity and diffusive transport. The porosity contrast number is another important dimensionless parameter that influences the dissolution front formation and instability growth [64],

D ¼ ð/f =/i Þ  1;

dimensionless numbers estimated, respectively. The Pe numbers are 1 for all the experiments, indicating that the advective transport is dominant over diffusion at a core scale. Da and Nk calculated with calibrated rate constants are about 3 to 60 times higher than those estimated in our companion paper with comparable Pe values [27]. The fact that Da  1 and Nk  1 also suggests that the characteristic times of reaction kinetics are much faster than those of advective and diffusive transport over a core length scale. Vuggy limestones have very large Da and Nk numbers that are about two or three orders of magnitude larger than those of Marly dolostones, because of high calcite reactivity (combining kinetic rate constants and specific surface areas). According to linear stability theory [64], large Damköhler (Da) numbers should yield sharp dissolution fronts that cause perturbations with short wavelengths to produce

Table 4 Comparison of fitted kinetic rate constants for calcite and dolomite at 60 °C, pCO2 = 3 MPa, and pH = 5 with those reported by Palandri and Kharaka [55] and Pokrovsky et al. [71]. Rater constants

a

Marly cores a Vuggy cores a Palandri and Kharaka [55] b Pokrovsky et al. [71]

ð21Þ

where /i and /f denote porosities ahead of and behind the dissolution front, respectively. Tables 5 and 6 summarize the major parameters used in the dimensionless-number calculations, and values of the

a b

Combined acid plus neutral mechanisms. pCO2 Mechanism.

kcalcite

kdolomite

106.40 104.27 105.08 103.52

105.92 105.92 105.02 105.54

406

Y. Hao et al. / Advances in Water Resources 62 (2013) 388–408

Fig. 18. Spatial and temporal development of wormhole and fingering in M-3 dolostones, using the porosity–permeability relationship with n = 8, after CO2 flooding of (a) 4, (b) 8, (c) 12, and (d) 120 h. The figures show the predicted porosity change with time.

Table 5 Major parameters used in dimensionless-number calculations [27]. Parameter

Value

Core length Darcy velocity Injected brine pH

0.03 m 4.716  106 m/s 4.3; 4.5; 4.8; 4.9 (for 3; 2; 1; 0.5 MPa) 105.86; 106.03; 106.27; 106.34 mol/m2 s (for Marly cores 3; 2; 1; 0.5 MPa) 103.60; 103.79; 104.08; 104.18 mol/m2 s (for Vuggy cores 3; 2; 1; 0.5 MPa) 105.62; 105.71; 105.84; 105.88 mol/m2 s (for Marly cores 3; 2; 1; 0.5 MPa) 105.62; 105.71; 105.84; 105.88 mol/m2 s (for Vuggy cores 3; 2; 1; 0.5 MPa) 0.053; 0.049; 0.047; 0.044 mol/L (for 3; 2; 1; 0.5 MPa) 0.019; 0.018; 0.017; 0.0164 mol/L (for 3; 2; 1; 0.5 MPa) 2  109 m2/s

a

kcalcite

a

kdolomite

Ceq,calcite Ceq,dolomite Diffusion coefficient

Note: Reactive mineral specific surface areas were obtained from [27]. a Calibrated carbonate reaction rates combining both acid and neutral mechanisms (Table 3, this work).

unstable fronts when advection is dominant. Our results are consistent with linear stability analysis. Unstable reactions fronts developed in Vuggy samples (Da  104) are much sharper and narrower than those in Marly samples (Da  102). Reaction yielded a longer-wavelength disturbance superimposed on steadily advanced dissolution fronts in the Marly cores, which eventually led to the development of a major preferential flow path (Fig. 17). In contrast to Marly cores, instabilities in the Vuggy cores

Table 6 Dimensionless number values. Core ID

M-0.5 M-1 M-2 M-3 V-0.5 V-1 V-2 V-3

Dimensionless number values Pe

Da

Nk

D

2.0  102 2.3  102 2.3  102 2.0  102 4.7  102 5.4  102 4.4  102 4.7  102

2.1  102 2.9  102 3.8  102 2.4  102 1.6  104 2.0  104 3.2  104 6.1  104

4.4  104 6.6  104 8.7  104 4.8  104 7.5  106 1.1  107 1.4  107 2.9  107

1.9 2.2 2.2 1.8 5.7 6.7 5.3 5.7

Note: The core length was selected as the characteristic length scale L in the calculations.

developed instantaneously and aggressively. These different dissolution behaviors and patterns are not only due to larger porosity contrast and higher Damköhler numbers as explained by linear stability analysis [64], but also may be attributed to strong heterogeneity present in Vuggy limestones as compared to the Marly dolostones. The flow velocity varies greatly due to strong spatial variations in porosity and permeability in the samples, thereby leading to significant variability in local Péclet and Damköhler numbers. Therefore a simple non-dimensional number may not adequately characterize the ultimate unstable dissolution patterns for a more complex and heterogeneous rock sample like Vuggy limestone [28].

4.5. Limitations of the continuum description of reactive transport and wormhole development A fundamental assumption for continuum-scale reactive transport models is that chemical species are well mixed at the porescale, leading to a uniform distribution of species concentrations and reaction progression over the larger REV that forms the building blocks of the model. Theoretical modeling studies require diffusion to be the dominant transport mechanism at the pore scale to satisfy this well-mixed assumption (e.g., [75]). The dimensional analysis performed by Smith et al. [27] showed that at most pore scales, diffusion dominates over reaction with kinetic number Nk 1 in which Nk measures the relative importance of reaction over diffusion at a pore scale, leading to uniform species concentrations for most of the pore dimensions found in our dolostone and limestone samples. Low Péclet (Pe) values [27], indicating limited effects of advection compared to diffusion, also imply a complete pore-scale mixing in the cores we considered [34]. However, the transition from porous flow to channel flow (particularly after wormhole development in the Vuggy limestones) resulted in a larger pore space or effective flow cross-section (the wormhole diameter as opposed to pore throat dimensions). At this point, the assumption of a well-mixed pore fluid may no longer be satisfied and fluid concentration gradients from the center of the channel to the walls may develop. Therefore the use of a Darcy flow approximation may not be adequate to describe flow and transport within a channel/wormhole where reactive fluid no longer has sufficient time and exposure to carbonate surfaces to buffer pH. This may explain the over-prediction of the final pH values by the model once breakthrough was attained (Fig. 16(c)).

Y. Hao et al. / Advances in Water Resources 62 (2013) 388–408

407

5. Conclusions

References

We have compared three-dimensional modeling results with solution chemistry, pressure, and characterization data from core flood experiments to quantify the impacts of fluid flow, carbonate reactivity, and formation heterogeneity on the development of carbonate dissolution front. The good agreements between model and experiment showed that continuum scale models can be used to describe evolution of porosity and permeability caused by the dissolution of carbonate minerals in dolostone and limestone cores of varying heterogeneity. We have examined the role of heterogeneity on the initiation and development of dissolution front. Estimation of the effective permeability within and between representative element volumes was necessary in order to accurately model the evolution of porosity and permeability in response to carbonate mineral dissolution. The model could reproduce stable dissolution fronts with small dissolution fingers that eventually evolved into conical shaped wormholes when the simulation was extended beyond the experimental run time in largely homogenous Marly cores. As heterogeneity increased and the effective permeability was more variable in Vuggy limestone samples, the model was capable of producing conical and/or well-defined wormholes only when the initial anisotropic permeability was explicitly described and a higher order dependence on porosity (n = 6–8) was used. Although the model parameter calibration is specific to the Marly dolostone and Vuggy limestone flow units found at the IEAGHG Weyburn-Midale CO2 Monitoring and Storage Project (Saskatchewan, Canada), several aspects may be transferable to other sites with similar mineralogy and heterogeneity. The numerical work presented in this study provides a calibrated benchmark model, not only allowing for an improved understanding of the chemical and physical behaviors arising from carbonate dissolution, particularly in lower-permeability carbonate formations, but also developing a useful basis for upscaling flow, transport and reaction properties from the core to intermediate or reservoir scale. The various dissolution features could be adequately described by (1) volume averaging pore space, mineralogy, and accessible surface area from the micron- to the grid scale; (2) assessing effective permeability within and between grids; (3) using reaction affinity-based kinetic equations; and (4) using power-law porosity–permeability and porosity-surface area/mineral abundance correlations. Application of this model to other carbonate sites may require additional calibration to constrain the reaction rates and powerlaw correlations used in reactive-transport simulations to account for changes to porosity (and permeability) as a consequence of calcite and dolomite dissolution. Although the simulation results shown here are for small cores, differences in the two rock types suggest that it may be possible to develop heterogeneity indices that could be used to constrain parameter input at larger scales.

[1] Hoefner ML, Fogler HS. Pore evolution and channel formation during flow and reaction in porous media. AlChE J 1988;34:45–54. http://dx.doi.org/10.1002/ aic.690340107. [2] Hung KM, Hill AD, Sepehrnoori K. A mechanistic model of wormhole growth in carbonate matrix acidizing and acid fracturing. J Petrol Technol 1989;41:59–66. http://dx.doi.org/10.2118/16886-PA. [3] Daccord G, Lenormand R, Liétard O. Chemical dissolution of a porous medium by a reactive fluid – I. Model for the ‘‘wormholing’’ phenomenon. Chem Eng Sci 1993;48:169–78. http://dx.doi.org/10.1016/0009-2509(93)80293-Y. [4] Daccord G, Liétard O, Lenormand R. Chemical dissolution of a porous medium by a reactive fluid – II. Convection vs reaction, behavior diagram. Chem Eng Sci 1993;48:179–86. http://dx.doi.org/10.1016/0009-2509(93)80294-Z. [5] Huang T, Hill AD, Schechter RS. Reaction rate and fluid loss: the keys to wormhole initiation and propagation in carbonate acidizing. In: SPE 37312, Int Symp on Oilfield Chemistry. Houston, TX; Feb 18–21, 1997. http://dx.doi.org/ 10.2118/37312-MS. [6] Fredd CN, Fogler HS. Influence of transport and reaction on wormhole formation in porous media. Fluid Mech Transp Phenom 1998;44:1933–49. http://dx.doi.org/10.1002/aic.690440902. [7] Golfier F, Zarcone C, Bazin B, Lenormand R, Lasseux D, Quintard M. On the ability of a Darcy-scale model to capture wormhole formation during the dissolution of a porous medium. J Fluid Mech 2002;457:213–54. http:// dx.doi.org/10.1017/S0022112002007735. [8] Golfier F, Bazin B, Lenormand R, Quintard M. Core-scale description of porous media dissolution during acid injection – Part I: Theoretical development. Comput Appl Math 2004;23(2–3):173–94. http://dx.doi.org/10.1590/S010182052004000200005. [9] Golfier F, Quintard M, Bazin B, Lenormand R. Core-scale description of porous media dissolution during Acid Injection – Part II: Calculation of the effective properties. Comput Appl Math 2006;25:55–78. [10] Panga MKR, Ziauddin M, Balakotaiah V. Two-scale continuum model for simulation of wormholes in carbonate acidization. AlChE J 2005;51:3231–48. http://dx.doi.org/10.1002/aic.10574. [11] Maheshwari P, Ratnakar RR, Kalia N, Balakotaiah V. 3-D simulation and analysis of reactive dissolution and wormhole formation in carbonate rocks. Chem Eng Sci 2013;90:258–74. http://dx.doi.org/10.1016/j.ces.2012.12.032. [12] Singurindy O, Berkowitz B. The role of fractures on coupled dissolution and precipitation patterns in carbonate rocks. Adv Water Resour 2005;28:507–21. http://dx.doi.org/10.1016/j.advwatres.2005.01.002. [13] Kalia N, Balakotaiah V. Modeling and analysis of wormhole formation in reactive dissolution of carbonate rocks. Chem Eng Sci 2007;62:919–28. http:// dx.doi.org/10.1016/j.ces.2006.10.021. [14] Kalia N, Balakotaiah V. Effect of medium heterogeneities on reactive dissolution of carbonates. Chem Eng Sci 2009;64:376–90. http://dx.doi.org/ 10.1016/j.ces.2008.10.026. [15] Cohen CE, Ding D, Quintard M, Bazin B. From pore scale to wellbore scale: impact of geometry on wormhole growth in carbonate acidization. Chem Eng Sci 2008;63:3088–99. http://dx.doi.org/10.1016/j.ces.2008.03.021. [16] Izgec O, Zhu D, Hill, AD. Models and methods for understanding of early acid breakthrough observed in acid core-floods of vuggy carbonates. In: SPE 122357, 2009 SPE European Formation Damage Conference; Scheveningen, The Netherlands; May 27–29, 2009. http://dx.doi.org/10.2118/122357-MS. [17] McDuff D, Jackson S, Shuchart C, Postl D. Understanding wormholes in carbonates: unprecedented experimental scale and 3D visualization. J Petrol Technol 2010;62:78–81. http://dx.doi.org/10.2118/129329-MS. [18] Izgec O, Demiral B, Bertin H, Akin S. CO2 injection into saline carbonate aquifer formations. I: Laboratory investigation. Transp Porous Med 2008;72:1–24. http://dx.doi.org/10.1007/s11242-007-9132-5. [19] Izgec O, Demiral B, Bertin H, Akin S. CO2 injection into saline carbonate aquifer formations. II: Comparison of numerical simulations to experiments. Transp Porous Med 2008;73:57–74. http://dx.doi.org/10.1007/s11242-007-9160-1. [20] Grigg RB, McPherson BJ., Svec RK. Laboratory and model tests at reservoir conditions for CO2-brine-carbonate rock systems interactions. In: The Second Annual DOE Carbon Sequestration Conference, Washington, DC; May 5–8, 2003. [21] Luquot L, Gouze P. Experimental determination of porosity and permeability changes induced by injection of CO2 into carbonate rocks. Chem Geol 2009;265:148–59. http://dx.doi.org/10.1016/j.chemgeo.2009.03.028. [22] Bacci G, Korre A, Durucan S. An experimental and numerical investigation into the impact of dissolution/precipitation mechanisms on CO2 injectivity in the wellbore and far field regions. Int J Greenhouse Gas Control 2011;5:579–88. http://dx.doi.org/10.1016/j.ijggc.2010.05.007. [23] Ellis B, Peters C, Fitts J, Bromhal G, McIntyre D, Warzinski R, Rosenbaum E. Deterioration of a fractured carbonate caprock exposed to CO2-acidified brine flow. Green Gases Sci Technol 2011;1:248–60. http://dx.doi.org/10.1002/ ghg.25. [24] Gouze P, Luquot L. X-ray microtomography characterization of porosity, permeability and reactive surface changes during dissolution. J Contam Hydrol 2011;120–121:45–55. http://dx.doi.org/10.1016/j.jconhyd.2010.07.004. [25] Smith MM, Sholokhova Y, Hao Y, Carroll SA. Evaporite caprock integrity: an experimental study of reactive mineralogy and pore-scale heterogeneity during brine-CO2 exposure. Environ Sci Technol 2013;47(1):262–8. http:// dx.doi.org/10.1021/es3012723.

Acknowledgments This work was performed under the auspices of the US Department of Energy by Lawrence Livermore National Laboratory under Contract DE-AC52-07NA27344. Funding for this work was provided by the Petroleum Technology Research Centre (PTRC), IEA GHG Weyburn CO2 Monitoring and Storage Project and the US DOE, Office of Fossil Energy, Carbon Sequestration Program. The authors would like to thank the editor and anonymous reviewers for their constructive comments and suggestions.

408

Y. Hao et al. / Advances in Water Resources 62 (2013) 388–408

[26] Carroll S, Hao Y, Smith M, Sholokhova Y. Development of scaling parameters to describe CO2 -rock interactions within Weyburn-Midale carbonate flow units. Int J Greenhouse Gas Control 2013. http://dx.doi.org/10.1016/ j.ijggc.2012.12.026. [27] Smith M, Sholokhova Y, Hao Y, Carroll S. CO2-induced dissolution of low permeability carbonates. Part I: Characterization and Experiments. Adv Water Res. 2013;62:370–87. http://dx.doi.org/10.1016/j.advwatres.2013.09.008. [28] Steefel CI. Geochemical kinetics and transport. In: Brantley SL, Kubicki JD, White AF, editors. Kinetics of water–rock interaction. New York: Springer Science; 2008. p. 545–89. http://dx.doi.org/10.1007/978-0-387-73563-4_11. [29] Steefel CI, DePaolo DJ, Lichtner PC. Reactive transport modeling: an essential tool and a new research approach for the earth sciences. Earth Planet Sci Lett 2005;240:539–58. http://dx.doi.org/10.1016/j.epsl.2005.09.017. [30] Szecsody JE, Zachara JM, Chilakapati A, Jardine PM, Ferrency AS. Importance of flow and particle-scale heterogeneity on CoII/III EDTA reactive transport. J Hydrol 1998;209:112–36. http://dx.doi.org/10.1016/S0022-1694(98)00114-0. [31] White AF, Brantley SL. The effect of time on the weathering of silicate minerals: why do weathering rates differ in the laboratory and field? Chem Geol 2003;190:69–89. http://dx.doi.org/10.1016/j.chemgeo.2003.03.001. [32] Li L, Peters CA, Celia MA. Upscaling geochemical reaction rates using porescale network modeling. Adv Water Resour 2006;29:1351–70. http:// dx.doi.org/10.1016/j.advwatres.2005.10.011. [33] Li L, Peters CA, Celia MA. Effects of mineral spatial distributions on reaction rates in porous media. Water Resour Res 2007;43:W01419. http://dx.doi.org/ 10.1029/2005WR004848. [34] Li L, Steefel CI, Yang L. Scale dependence of mineral dissolution rates within single pores and fractures. Geochem Cosmochem Acta 2008;72:360–77. http://dx.doi.org/10.1016/j.gca.2007.10.027. [35] Li L, Gawande N, Kowalsky MB, Steefel CI, Hubbard SS. Physicochemical heterogeneity controls of uranium, bioreduction rates at the field scale. Environ Sci Technol 2011;45(23):9959–66. http://dx.doi.org/10.1021/ es201111y. [36] Molins S, Trebotich D, Steefel CI, Shen C. An investigation of the effect of pore scale flow on average geochemical reaction rates using direct numerical simulation. Water Resour Res 2012;48:W03527. http://dx.doi.org/10.1029/ 2011WR011404. [37] Cirpka OA, Valocchi AJ. Two-dimensional concentration distribution for mixing-controlled bioreactive transport in steady state. Adv Water Resour 2007;30:1668–79. http://dx.doi.org/10.1016/j.advwatres.2006.05.022. [38] Navarre-Sitchler A, Steefel CI, Yang L, Tomutsa L, Brantley SL. Evolution of porosity and diffusivity associated with chemical weathering of a basalt clast. J Geophys Res Earth Surf 2009;114:F02016. http://dx.doi.org/10.1029/ 2008JF001060. [39] Werth CJ, Cirpka OA, Grathwohl P. Enhanced mixing and reaction through flow focusing in heterogeneous porous media. Water Resour Res 2006;42:WR004511. http://dx.doi.org/10.1029/2005WR004511. [40] Knutson C, Valocchi A, Werth C. Comparison of continuum and pore-scale models of nutrient biodegradation under transverse mixing conditions. Adv Water Resour 2007;30:1421–31. http://dx.doi.org/10.1016/ j.advwatres.2006.05.012. [41] Wolery TJ. EQ3/6, a software package for geochemical modeling of aqueous systems: package overview and installation guide (version 7.0). Livermore, CA: Lawrence Livermore National Laboratory; UCRL-MA-110662-PT-1; 1992. [42] Johnson JW, Oelkers EH, Helgeson HC. SUPCRT92: A software package for calculating the standard molal thermodynamic properties of minerals, gases, aqueous species, and reaction from 1 to 5000 bar and 0° to 1000 °C. Comput Geosci 1992;98:131–9. http://dx.doi.org/10.1016/0098-3004(92)90029-Q. [43] Helgeson HC. Thermodynamics of hydrothermal systems at elevated temperatures and pressures. Am J Sci 1969;267:729–804. http://dx.doi.org/ 10.2475/ajs.267.7.729. [44] Drummond SE. Boiling and mixing of hydrothermal fluids: chemical effects on mineral precipitation. PhD dissertation, Pennsylvania State University; 1981. [45] Nitao JJ. Reference manual for the NUFT flow and transport code, version 2.0. Livermore, CA: Lawrence Livermore National Laboratory; UCRL-MA-130651; 1998. [46] Hao Y, Sun Y, Nitao JJ. Overview of NUFT: a versatile numerical model for simulating flow and reactive transport in porous media. In: Zhang F, Yeh G, Parker JC, editors. Groundwater reactive transport models. Bentham Science Publishers; 2012. p. 219–32. http://dx.doi.org/10.2174/ 97816080530631120101. [47] Buscheck TA, Glascoe LG, Lee KH, Gansemer J, Sun Y, Mansoor K. Validation of the multiscale thermohydrologic model used for analysis of a proposed repository at Yucca Mountain. J Contam Hydrol 2003;62–3:421–40. http:// dx.doi.org/10.1016/S0169-7722(02)00157-2. [48] Glassley WE, Nitao JJ, Grant CW. Three-dimensional spatial variability of chemical properties around a monitored waste emplacement tunnel. J Contam Hydrol 2003;62–3:495–507. http://dx.doi.org/10.1016/S01697722(02)00153-5. [49] Johnson JW, Nitao JJ, Knauss KG. Reactive transport modeling of CO2 storage in saline aquifers to elucidate fundamental processes, trapping mechanisms and sequestration partitioning. In: Baines SJ, Worden RH, editors. Geological storage of carbon dioxide, vol. 223. Geological Society, London: Special Publications; 2004. p. 107–28. [50] Carroll SA, Hao Y, Aines R. Geochemical detection of carbon dioxide in dilute aquifers. Geochem Trans 2009;10:4. http://dx.doi.org/10.1186/1467-4866-104.

[51] Morris JP, Hao Y, Foxall W, McNab W. A study of injection-induced mechanical deformation at the in Salah CO2 storage project. Int J Greenhouse Gas Control 2011;5:270–80. http://dx.doi.org/10.1016/j.ijggc.2010.10.004. [52] Morris JP, Detwiler RL, Friedmann SJ, Vorobiev OY, Hao Y. The large-scale geomechanical and hydrogeological effects of multiple CO2 injection sites on formation stability. Int J Greenhouse Gas Control 2011;5:69–74. http:// dx.doi.org/10.1016/j.ijggc.2010.07.006. [53] Lu C, Sun Y, Buscheck TA, Hao Y, White JA, Chiaramonte L. Uncertainty quantifi cation of CO2 leakage through a fault with multiphase and nonisothermal effects. Green Gases Sci Technol 2012;2:445–59. http://dx.doi.org/10.1002/ ghg.1309. [54] Sun Y, Tong C, Trainor-Guitton WJ, Lu C, Mansoor K, Carroll SA. Global sampling for integrating physics-specific subsystems and quantifying uncertainties of CO2 geologic sequestration. Int J Greenhouse Gas Control 2013;12:108–23. http://dx.doi.org/10.1016/j.ijggc.2012.10.004. [55] Palandri JL, Kharaka YK. A compilation of rate parameters of water-mineral interaction kinetics for application to geochemical modeling. Menlo Park, CA: US Geological Survey; Open File Report 2004-1068; 2004. [56] Drever JI. The geochemistry of natural waters. Englewood Cliffs, NJ: Printice Hall; 1988. [57] Land LS. Failure to precipitate dolomite at 25 °C from dilute solution despite 1000-fold oversaturation after 32 years. Geochemistry of low temperature processes. Aquat Geochem 1998;4(3-4):361–8. http://dx.doi.org/10.1023/ A:1009688315854. [58] Civan F. Scale effect on porosity and permeability: kinetics, model, and correlation. AlChE J 2001;47:271–87. http://dx.doi.org/10.1002/ aic.690470206. [59] Noiriel C, Luquot L, Madé B, Raimbault L, Gouze P, van der Lee J. Changes in reactive surface area during limestone dissolution: an experimental and modelling study. Chem Geol 2009;265:160–70. http://dx.doi.org/10.1016/ j.chemgeo.2009.01.032. [60] Emmanuel S, Berkowitz B. Mixing-induced precipitation and porosity evolution in porous media. Adv Water Resour 2005;28:337–44. http:// dx.doi.org/10.1016/j.advwatres.2004.11.010. [61] Lichtner PC. The quasi-staionary state approximation to coupled mass transport and fluid-rock interaction in a porous medium. Geochim Cosmochim Acta 1988;52:143–65. http://dx.doi.org/10.1016/00167037(88)90063-4. [62] Ortoleva P, Merino E, Moore C, Chadam J. Geochemical self-organization. I: Reaction-transport feedbacks and modeling approach. Am J Sci 1987;287:979–1007. http://dx.doi.org/10.2475/ajs.287.10.979. [63] Ortoleva P, Chadam J, Merino E, Sen A. Geochemical self-organization. II: The reactive-infiltrative instability. Am J Sci 1987;287:1008–40. http://dx.doi.org/ 10.2475/ajs.287.10.1008. [64] Szymczak P, Ladd AJC. Instabilities in the dissolution of a porous matrix. Geophys Res Lett 2011;38:L07403. http://dx.doi.org/10.1029/2011GL046720. [65] Stauffer D. Introduction to percolation theory. London, Philadelphia: Taylor and Francis; 1985. [66] Bentz DP, Garboczi EJ. Percolation of phases in a three-dimensional cement paste microstructural model. Cem Concr Res 1991;21(2):325–44. http:// dx.doi.org/10.1016/0008-8846(91)90014-9. [67] Botros FE, Hassan AE, Reeves DM, Pohll G. On mapping fracture networks onto continuum. Water Resour Res 2008;44:W08435. http://dx.doi.org/10.1029/ 2007WR006092. [68] Reeves DM, Benson DA, Meerschaert MM. Transport of conservative solutes in simulated fracture networks. 1: Synthetic data generation. Water Resour Res 2008;44:W05404. http://dx.doi.org/10.1029/2007WR006069. [69] Manzocchi T, Walsh JJ, Nell P, Yielding G. Fault transmissibility multipliers for flow simulation models. Petrol Geosci 1999;5:53–63. http://dx.doi.org/ 10.1144/petgeo.5.1.53. [70] Renard P, LeLoc’h G, Ledoux E, de Marsily G, Mackay R. A fast algorithm for the estimation of the equivalent hydraulic conductivity of heterogeneous media. Water Resour Res 2000;36(12):3567–80. http://dx.doi.org/10.1029/ 2000WR900203. [71] Pokrovsky OS, Golubev SV, Schott J, Castillo A. Calcite, dolomite and magnesite dissolution kinetics in aqueous solutions at acid to circumneutral pH, 25 to 150 °C and 1 to 55 atm pCO2: New constraints on CO2 sequestration in sedimentary basins. Chem Geol 2009;265:20–32. http://dx.doi.org/10.1016/ j.chemgeo.2009.01.013. [72] Lichtner PC, Kang Q. Upscaling pore-scale reactive transport equations using a multi-scale continuum formulation. Water Resour Res 2007;43:W12S15. http://dx.doi.org/10.1029/2006WR005664. [73] Steefel CI, Lasaga AC. Evolution of dissolution patterns: Permeability change due to coupled flow and reaction. In: Melchior D, Bassett RL, editors. Chemical Modeling of Aqueous Systems II, ACS Symposium Series No. 416. Washington, DC: American Chemical Society; 1990. p. 212–25. [74] Steefel CI, Maher K. Fluid–rock interaction: a reactive transport approach. Rev Mineral Geochem 2009;70:485–532. http://dx.doi.org/10.2138/ rmg.2009.70.11. [75] Dentz M, Le Borgne T, Englert A, Bijeljic B. Mixing, Spreading, and reaction in heterogeneous media: a brief review. J Contam Hydrol 2011;120–1:1–17. http://dx.doi.org/10.1016/j.jconhyd.2010.05.002.