Effects of geochemical reaction on double diffusive natural convection of CO2 in brine saturated geothermal reservoir

International Journal of Heat and Mass Transfer 77 (2014) 519–528

Contents lists available at ScienceDirect

International Journal of Heat and Mass Transfer journal homepage: www.elsevier.com/locate/ijhmt

Effects of geochemical reaction on double diffusive natural convection of CO2 in brine saturated geothermal reservoir Akand Islam ⇑, Aboulghasem Kazemi Nia Korrani, Kamy Sepehrnoori, Tad Patzek Department of Petroleum & Geosystems Engineering, The University of Texas at Austin, TX 78712, USA

a r t i c l e

i n f o

Article history: Received 21 December 2013 Received in revised form 19 May 2014 Accepted 20 May 2014

Keywords: Geochemical reaction Double diffusion Natural convection Geothermal effect Heterogeneity CO2 sequestration

a b s t r a c t In order for mitigation of greenhouse gas emissions, as we know, the storage of CO2 in deep geological formations is considered as a promising option. Besides convective mixing an understanding of geochemical reaction that affects the long term storage process in deep saline aquifers is of great importance. In this study we numerically investigate double diffusive natural convection of CO2 in brine saturated geothermal reservoir which is anisotropic in permeability variations, impervious from the sides, and is open to CO2 at the top. Our primary objective is to understand effects of geochemical reactions of different rates and orders on density driven natural convection of CO2 due to concentration and geothermal gradients. We present propagation of CO2 plumes over long period of time by analyzing different combination of problem parameters: Damkohler number (0.01 6 Da 6 105) with variation of 1st, 2nd, and 3rd order reactions; solutal Rayleigh number (500 6 Ras 6 2000); the buoyancy ratio (2 6 N 6 100); Dykstra– Parsons coefficient (0.35 6 Vdp 6 0.85); and fixed Lewis and Prandtl numbers. Reaction order is set from relevant stoichiometric ratio. In each case results are quantified in term of percentage of CO2 deposition. It is found that mineral interactions make traceable difference of depositions as reaction rate increases, especially when Da > 103. In one hand, compared to very minimum geochemical interaction (Da  0.01) the strong reaction effects (Da  104) can make difference of more than 5% in the period of 500 years of trapping. On the other hand, at a fixed equivalent Da reaction order also makes substantial distinction as deposition time passes. Heterogeneity plays a vital role. Geothermal gradient has very negligible effect, however only after very long time (>500 years). Ó 2014 Elsevier Ltd. All rights reserved.

1. Introduction

CO2 ðgÞ ¼ CO2 ðaqÞ

Deep saline aquifers have the highest storage potential to bridge the gap between climate change mitigation and CO2 emissions [30,15,1,18,22,2,20]. Physical, residual, solubility, and mineralogic trapping are the proposed mechanisms by which CO2 is hold in the saline aquifers and does not escape to the surface [16,29,2,24]. Among these mechanisms, mineralogic trapping is the only technique that sequesters the CO2 while others are referred as storage mechanisms [10]. Hence, modeling of geochemical processes is a necessary predictive tool for long-term CO2 disposal in deep aquifers. When CO2 is injected into a saline aquifer, it can react with other chemical species in aqueous phase as well as constituents on the solid phase. The geochemical changes caused by CO2 injection into aquifers include acidification and carbonation of the native brine and potential mineral dissolution and precipitation reactions driven by the aqueous composition changes [25]:

H2 O þ CO2 ðaqÞ ¼ H2 CO3

⇑ Corresponding author. Tel.: +1 5124365343. E-mail address: [email protected] (A. Islam). http://dx.doi.org/10.1016/j.ijheatmasstransfer.2014.05.040 0017-9310/Ó 2014 Elsevier Ltd. All rights reserved.

H2 CO3 ¼ Hþ þ HCO3 CaCO3 ðSÞ þ Hþ ¼ Caþ2 þ HCO3 High HCO 3 concentration in the aqueous phase is a driving force for CO2 to be dissolved more in the aqueous phase [21]. Such reactions can change the hydrological and chemical properties of subsurface systems, alter flow pathways, and eventually affect contaminant transport and behavior significantly [33]. Mineral dissolution and precipitation reactions are important for evaluating the potential CO2 storage capacity and can also influence the performance of the injection well. Many researchers investigated the geochemistry in the CO2 sequestration, mostly related to kinetics problem [8]. Noh et al. [25] proposed a mathematical formulation by which geochemistry is integrated with the multi-phase flow. Local equilibrium was assumed in the formulation. The theory led to graphical

520

A. Islam et al. / International Journal of Heat and Mass Transfer 77 (2014) 519–528

Nomenclature c Cp D Da g h H k L n p k Le Pe Pr Ra S T t u V x z

concentration [mol/m3] heat capacity at constant pressure [J kg1 K1] Diffusion coefficient [m2/s] Damkohler number gravitational acceleration [m/s2] enthalpy [J kg1] porous medium height [m] permeability [m2] porous medium length [m] number of nodes pressure [Pa] wave number [–] Lewis number [–] Peclet number [–] Prandtl number [–] Rayleigh number [–] average deposition [%] Temperature [°C] time [s] velocity [m/s] Dykstra–Parsons coefficient distance along x-axis distance along z-axis

bT /

l j q w

coefficient of thermal expansion [K1] porosity [–] viscosity [kg m1 s1] thermal conductivity density [kg/m3] stream function [m3 m1 s1]

Superscript ⁄ dimensionless quantity Subscripts 0 initial value a stoichiometric coefficient of CO2 brine phase b stoichiometric coefficient of minerals dp Dykstra–Parsons i node in x – direction j node in z – direction m minerals s solutal x x-coordinate z z-coordinate

Greek letters a thermal diffusivity bc coefficient of density increase by concentration [m3/mol]

solution from which it is easy to see when and under what conditions mineralization occurs during the injection. Okamoto et al. [28] used GEM-GHG simulator and performed a sensitivity study of CO2 mineralization to investigate effective parameters on CO2 behavior over a long time. They used a real reservoir model of Nagaoka pilot test site in Japan. Reactive surface area and component of minerals were reported as essential parameters. Fischer et al. [6] conducted an experiment to study the significance of the CO2brine-rock interactions during the geological storage of the CO2. In this work, CO2-treated samples were compared against untreated samples to investigate the CO2 effects on the mineralogy of a sandstone rock. Their experimental results showed the stabilization or precipitation of albite and dissolution of calcium-rich plagioclase, K-feldspar and anhydrite. Wolf et al. [39] investigated the in situ reactions involved in CO2 sequestration through a micro-reaction system. In their study synchrotron X-ray Diffraction and Raman Spectroscopy were used to characterize the carbonation reaction. O’Connor et al. [27] also conducted an experiment to study mineral sequestration of CO2. Several suggestions were documented to accelerate the slow geological processes including carbon dioxide activity increase in solution or the reactive surface area and imperfections into the crystal lattice through high energy attrition grinding, activation of minerals through thermally, etc. Wellman et al. [38] investigated the CO2-brine-rock interactions by some corefloods, and then used TRANS-TOUGH simulator (combination of TOUGH-EOSCO2 with the chemistry code TRANS) to model their experimental work and to perform some sensitivity analysis. Lithology type and the magnitude of fluid flushed through the media were documented to be the two major controlling factors in this process. Rosenbauer et al. [31] showed that more CO2 might be sequestered in a deep saline aquifer not only through the mineral trapping but also by the formation of bicarbonate ion (HCO 3 ) in the aqueous phase. Druckenmiller et al. [3] reported the main

parameters in terms of the geochemistry for a CO2 sequestration process as, pressure; temperature; pH; brine concentration; and rock composition. It has been also included that the temperature effect is more pronounced compared to the pressure effect. Xu et al. [41] developed a reactive geochemical transport model in TOUGHREACT for a sandstone-shale system under high CO2 pressures. The model was then used to interpret the mass transfer of aqueous chemical components, the alteration pattern of minerals, sequestration of CO2, and changes of petrophysical properties (e.g., porosity) in a Gulf Coast aquifer. Schumacher and Anne [34] used Geochemist’s Workbench software (GWB) to analyze the CO2–water–rock interactions occurring during a carbon sequestration pilot test into a Mississippian oil reservoir. Several core and cutting samples as well as electric logs were used to model the static and dynamic of the pilot during the CO2 injection. Simulation results show net dissolution during the injection stages while net precipitation for post-CO2 injection. In recent studies [13,14] double diffusive natural convection of CO2 in brine saturated homogeneous and heterogeneous permeable geothermal reservoir were investigated. Here we have coupled the geochemistry and convective mixing in anisotropic homogeneous medium and have presented CO2 settlement results to understand the role of these processes on the long term fate of the sequestration. 2. The reservoir model and governing equations The hypothetical geothermal reservoir is considered to be a two-dimensional rectangular cavity filled with brine saturated porous medium with height H, and length L as shown in Fig. 1 Permeability varies spatially, i.e., k = k(x, z). There is no flow condition across all boundaries and no solute flux across lateral and bottom boundaries are assumed. There is no heat flux across the lateral

521

A. Islam et al. / International Journal of Heat and Mass Transfer 77 (2014) 519–528

  kz @p uz ¼   qg : l @z

ð3Þ

(c) Concentration transport equation:

@c @c @c @2c @2c / þ ux þ uz ¼ /D þ @t @x @z @x2 @z2

!

0

 k ca cbm :

ð4Þ

(d) Reaction kinetics equation:

@cr 0 / ¼ k ca cbm : @t

ð5Þ

where reaction is presented by:

aðCO2  brineÞ þ bðmÞ ! productðanother solid materialÞ

ð6Þ

(e) Energy equation:

!   @T @T @T @2T @2T þ Q; ¼j q0 C p þ ux þ uz þ @t @x @z @x2 @z2

Fig. 1. Hypothetical geothermal reservoir model.

boundaries while the top and bottom boundaries are fixed at isothermal cold and hot temperatures, respectively. Pressure change due to dissolution is not accounted. We assume that the reservoir contains only a single phase fluid (brine). This consideration is valid until the free-gas layer dissolves at which point there is no further source of solute. Capillary transition between gas and brine phases is ignored and therefore only the liquid phase is modeled. The presence of the gas phase is separated by a perturbed boundary on the top. The perturbed condition is discussed in next section. There is no CO2 dissolved in the modeled aquifer except in the interface. The Boussinesq approximation and Darcy flow model are followed. As a major concentration of this study we assume there is interaction between fluid phase and minerals in term of geochemical reaction. The stoichiometric coefficients of CO2 saturated brine and minerals are a and b, respectively and thus the reaction order is determined as a + b. For instance, equation of 0 2nd order reaction is r(c, cm) = k c1 c1m . Reactions may be more complicated in reality, however. The effect of velocity-induced dispersion is ignored. Since this is a natural convection problem Rayleigh number is low enough to expect only the laminar flow regime. Temperature gradient acts opposite to concentration difference. Density gradient due to simultaneous concentration and temperature differences are the driving force to impose natural convection process. Governing equations of flow- and concentration-field read as (a) Continuity equation:

@ q @ðqux Þ @ðquz Þ þ þ ¼ 0: @x @z @t kx @p ; l @x

where 0

Q ¼ k ca cbm Dh:

ð8Þ

0

k represents reaction rate constant. For geochemical reactions the heat of reaction (Dh) is assumed to be negligible [8]. Therefore the heat rate, Q becomes zero. For small density change due to temperature and concentration differences at constant pressure, the brine density is assumed to be a linear function of temperature and concentration of solute. This reads as

q ¼ q0 ½1 þ bc ðc  c0 Þ  bT ðT  T 0 Þ;

ð9Þ

where

bc ¼

  1 @q q0 @c T

and bT ¼ 

  1 @q : q0 @T c

ð10Þ

Differentiating with respect to x we obtain

  @q @c @T : ¼ q0 bc  bT @x @x @x

ð11Þ

After eliminating the pressure by cross-differentiation of Eqs. (2) and (3) we get

  1 @uz 1 @ux g q0 @c @T  ¼ bc  bT kz @x kx @z @x @x l

ð12Þ

ð1Þ

2.1. Dimensionless form of the equations

ð2Þ

We consider the cavity height, H, as the characteristic length and /D as the characteristic velocity. The dimensionless variables H are defined (shown for 2nd order reaction),

(b) Darcy’s law:

ux ¼ 

ð7Þ

Fig. 2. Logarithmic permeability distributions for Vdp of (a) 0.35, (b) 0.55, and (c) 0.85.

522

A. Islam et al. / International Journal of Heat and Mass Transfer 77 (2014) 519–528

Fig. 3. Concentration profiles for Da = 0.1, Ras = 500 at (a) t = 4, (b) t = 50, (c) t = 100, and (d) t = 500 years.

Table 1 CO2 deposition over time (Ras = 500). t (year)

4 20 100 500

S Da = 1000

Da = 10,000

0.18 0.65 2.37 9.28

0.18 0.75 2.79 10.45

x z Hux  Huz  H2  H2  w x ¼ ; z ¼ ; ux ¼ ;u ¼ ;k ¼ ;k ¼ ;w ¼ ; H H /D /D z /D x kx z kz Dt c  c0  c  c0  T  T 0  @w  @w ; cm ¼ ;T ¼ ;u ¼ ; t  ¼ 2 ; c ¼ 0 ; ux ¼  0 0 c0  c c0  c @z z @x T0  T H 0

Ras ¼

Dqc gkH Dq gkH a k H2 Dc ; RaT ¼ T ; Le ¼ ; Da ¼ : /Dl la /D uD ð13Þ

523

A. Islam et al. / International Journal of Heat and Mass Transfer 77 (2014) 519–528 Table 2 CO2 deposition over time (Ras = 2000). t (year)

S N = 100

N=2

Da

4 20 100 500

0.1

10

1000

10,000

0.1

10

1000

10,000

0.18 1.0 4.94 19.12

0.18 1.0 4.95 19.13

0.18 1.02 5.39 20.99

0.18 1.08 6.11 23.98

0.18 1.0 4.94 19.44

0.18 1.0 4.95 19.46

0.18 1.02 5.39 21.17

0.18 1.08 6.11 24.15

After applying the Boussinesq approximation the final forms of dimensionless equations are (eliminating ⁄)

kz

  @2w @2w @c 1 @T ; þ k ¼ Ra  x s @x2 @z2 @x N @x

ð14Þ

@c @c @c @ 2 c @ 2 c þ ux þ uz ¼ þ  Daccm ; @t @x @z @x2 @z2

ð15Þ

@cr ¼ Daccm ; @t

ð16Þ

! @T @T @T 1 @2T @2T ; þ uux þ uuz ¼ þ @t @x @z Pr @x2 @z2

ð17Þ q c D

Ras where N ¼ bbTc DDTc ¼ LeRa is the buoyancy ratio and Pr ¼ 0jp . The T boundary conditions are shown in Fig. 1 Initial conditions within reservoir geometry include c = 0 and cm = 1. Eqs. (14)–(17) are to be solved for obtaining w, c, cm, and T.

3. Numerical solution methods As we know numerical simulation of density-driven transport problem is sensitive to discretization errors. Hence, the following CFL condition was fulfilled throughout the simulation:

uDt P E 6 /Dx 2

ð18Þ

Dx where PE ¼ uD/ : The 81  81 grid cells were used, and the time step (Dt) was taken as 106 (equivalent to 29 days based on H = 100 m). We used a sequential implicit solution procedure to solve the elliptic Poisson equation (14), the concentration equation. (15), the reaction Eq. (16), and the energy equation (17). Eq. (14) was solved by Point Gauss–Seidel method with convergence criteria  iþ1 i  P of abs w ww 6 106 . The superscript i denotes iteration number. i

Eqs. (15)–(17) were solved by Alternating Direction Implicit (ADI) method where the convection terms were discretized using upwind differencing and the diffusion terms were discretized using central differencing. The developed code checked against literature results [36,12,13,14] was in satisfactory agreement. Grid independency test was also performed. In order to observe downward finger-like propagation of the CO2 fronts the interface concentration was required to be perturbed. Therefore a sinusoidal perturbation was used as the form

cðx; z ¼ 0; t ¼ 0Þ ¼ 1 þ 0:01 sinð2pkxÞ:

ð19Þ

The wave number, k, 24 was used. This is just for initiation of calculation. The long term behavior is independent of initial perturbation [5], however. In reality, pore-level perturbations and thermo-mechanical fluctuations cause the perturbation to start the finger-like plumes [19,11]. For imposing double diffusion the buoyancy ratio (N) of 2–100 were used. This represents very low to high geothermal gradients [13]. Because temperature at the bot-

tom of the reservoir is higher than at the top the lighter fluid moves upward. On the other hand, the flow of heavier fluid due to concentration gradient is downward, makes mixing favorable. Therefore geothermal reservoir, in general, would appear as a relatively preferable option to sequester CO2. The solutal Rayleigh number, Ras, was varied from 500 to 2000. This is the typical range of geologic formation [8]. The heterogeneous permeability fields were generated by our in-house simulator with an average permeability 100 mD for different Dykstra–Parson coefficients, Vdp = 0.35–0.85. This is to cover usual heterogeneities [4]. In general Damkohler number (Da) varies from 102 to 105 for CO2 sequestration process based on the type of aquifer, temperature, pressure and salinity of the formation brine [8]. Hence to obtain complete picture of geochemical reaction effects we performed simulations from low Da = 101 to high end Da = 104. As we have already assumed the heat of reaction is too low to be considered, temperature change due to reaction is ignored. The reactions would occur in hundreds to thousands years (see comprehensive review by Gaus et al. [7]) periods and thereby this is much unlikely there would be any visible change in short time span. It is noteworthy that, for simplicity, series of complex geochemical reactions can be regarded in a single form because there is always a controlling reaction during dissolution and precipitation of minerals. The reaction kinetics also depends on specific surface of minerals which we do not consider to be changed throughout the convection periods. This is due to the fact that the change in the concentration of minerals is not comparable with CO2 dissolution [7,8]. The porosity (u) of a typical reservoir, 0.30 [26] and diffusion coefficient (D), 4  109 m2/s [23] were applied. Prandtl number, Pr = 0.0062 was used for the energy equation [35]. The thermal diffusivity (a) was considered constant, 3.7  107 m2/s. The fixed Lewis number, Le = 310 was used [17]. Again, there is no change of permeability and porosity thanks to reaction [40,42]. 4. Results and discussions We quantify CO2 deposition (or saturation of the reservoir) owing to dissolution and geochemical reactions simultaneously by the formula

Pni 2 Pnj 2 i

j

ci;j

ðni 2Þðnj 2Þ

S ¼

þ

b a

1 1 þ ba

Pni 2 Pnj 2 i

j

cmi;j

!

ðni 2Þðnj 2Þ

 100:0

ð20Þ

where a and b represent stoichiometric coefficients of reaction between CO2 saturated brine and rock formation, respectively. The equation fulfills the initial (0%) and complete saturation (100%) conditions of the aquifer. In numerator the first term accounts for deposition due to convection (solubility trapping) and the second term within parenthesis contributes for chemical reactions (mineral trapping). All deposition values in subsequent discussions are reported in percentage. For first, second, and third order reaction values of (a, b) are (1, 0), (1, 1), and (1, 2), respectively. In the case of first order reaction quantitatively minerals are assumed to be too high that the term Dacbm would remain constant for entire periods of convection. The geothermal effects are accounted through varying N. It is apparent from Eq. (14) that N is high for less thermal effect, and vice versa. Different levels of heterogeneities are applied. Fig. 2 shows permeability distributions. Vdp = 0.35 represents low variations differ by little more one-order where Vdp = 0.85 displays huge heterogeneities varying by more than five-orders. CO2 transport in heterogeneous media is different than in homogeneous because permeability variations make velocity fluctuations time dependent and cause flow regimes to be distinguishable as fingering, dispersion, and channeling [37].

524

A. Islam et al. / International Journal of Heat and Mass Transfer 77 (2014) 519–528

Fig. 4. Concentration profiles for Da = 10,000, Ras = 500 at (a) t = 4, (b) t = 50, (c) t = 100, and (d) t = 500 years.

Based on the different level of heterogeneities the mixing zone exhibits different characteristics in the respective regimes. First we show propagation of CO2 fronts over 500 years in Fig. 3 for low Da = 0.1 at Ras = 500. The results are produced for 0.55

heterogeneity and applying buoyancy ratio of 100. It is observed that during early periods, dissolution is the primarily diffusion dominated and therefore CO2 propagates in brine extremely slowly. After about 4 and 20 years the deposition is found to be

A. Islam et al. / International Journal of Heat and Mass Transfer 77 (2014) 519–528

525

Table 3 CO2 deposition over time for different heterogeneities (Ras = 1000, N = 50). t (year)

4 20 100 500

S Vdp = 0.35

Vdp = 0.85

0.18 0.63 2.1 9.47

0.99 4.09 16.27 43.28

0.18 and 0.63, respectively. As time passes the deposition gets little acceleration. Initial perturbations of the wavelength of 24 were imposed and therefore convection cells started to evolve by 23 plumes. The convection cells are surrounded by multiple diffusive shear layers acting as the feeding sites. The cells are attached together through these sites. The merging process starts from sides rather than at the tip of the fingers and thus they spread to the rest of the system. The merging of fingers involves horizontal pressure gradients created by differences in the upward flow, making the fingers unstable to perturbations in the horizontal wavelength [9]. As convection mixing becomes pronounced more cells bend making fluid movement rigorous. Some fingers grow faster than the others. The finger wavelength amplifies due to random but stronger nonlinear interactions. Many of the fingers seem to undergo strong interactions while others in the process of disappearing. Indeed, there is some tenacity of the initial behavior and the convection pattern experienced persists for few moments before the number of fingers starts to decrease reflecting intrinsic properties. In order to understand the dynamics of convective mixing more clearly readers are recommended to watch one of our simulation results from YouTube (http://www.youtube.com/ watch?v=mn9oOWOtSks). Because the top-left corner is relatively of high permeable CO2 plumes get favorable paths to transport and therefore fingers move faster close to the sidewalls than the central regions. By 100 years of mixing process we see two big plumes survive making reservoir 2.09 deposited. After long 100 years of deposition the convection process tend to become mixed diffusive-convective. Because the convection is too slow even after 500 years deposition is far less than one-tenth of complete saturation. The saturation reaches to 7.35. The unstable region is extended to larger area. For one-hundred times higher value of Da = 10, despite the geochemical interaction is expected to be intense numerical results remain almost unchanged just increasing the deposition slightly (7.41 after 500 years). Dynamics of mixing remains same. Small Ras indicates low density difference causing slow convection and therefore there is not much CO2 saturated brine available to react with minerals. Only if the reaction occurs in high enough rates can render some momentum. When Da is increased to 1000 or so the difference is traceable as we see in Table 1. Deposition can be enhanced to 10.45 provided reaction rate increases such that it varies from 102 to 104. The mineral interactions vanish saturated CO2 driving more to convect down resulting in overall trapping boost. Now, if the geothermal gradient is assumed to increase substantially, meaning the value of N is reduced to 2, the simulation results exhibit almost no difference. This turns out that at fixed low Ras, only Da, provided variations by orders-of-magnitude, can play role in differentiating the deposition. The effect of N may be traceable, however negligibly, only after long period of time at high Ras [17,14]. Table 2 shows results of deposition for higher Ras = 2000 from which we can clearly differentiate the impacts compared to the previous case (Ras = 500). The level of geochemical reaction becomes more effective as Ras increases. The concentration maps of Da = 104 is shown in Fig. 4 It shows that the fingering patterns

Fig. 5. Concentration profiles for Da = 100, Ras = 1000 at (a) t = 4, (b) t = 50, (c) t = 100, and (d) t = 500 years.

are dramatically slowing down indicating the change is strong, however makes the overall deposition fast. The consumption of

526

A. Islam et al. / International Journal of Heat and Mass Transfer 77 (2014) 519–528

Fig. 6. Deposition at different reaction orders (Ras = 1000, N = 2, Vdp = 0.55).

dissolved CO2 due to reactions leads to stratification of concentration gradients. At high Da the boundary layer remains about diffusive for long period of times. The flow system would linger the instability as the diffusion time for the large Da increases. Depositions can be more than 5% high when reaction with formation rocks reaches extreme rendering favorable trapping. The difference of 5% over 500 years period is significant sensu lato where the geochemical reactions as well as the CO2 trapping occur in the span of thousands to hundreds of thousands of years [7]. To examine geochemical effects in different heterogeneous permeability levels depositions computed are reported in Table 3. Results are shown for Da = 1000. It seems obvious that the geochemical interactions will change a lot from a reservoir of low to high heterogeneous. Overall trapping will be enhanced much in a reservoir of which local permeability differs by orders-of-magnitude. For heterogeneity increased the convection cells get more dynamicity because of availability of some local high permeable favorable paths and thus CO2 dissolution swells. This brings more CO2 saturated brine in contact with minerals helping settle down fast. For large heterogeneity dispersion phenomenon starts in a short time [4] which also a good reason to obtain more deposition than the case of low heterogeneity. The CO2 plumes are dominated by the permeability distributions pattern. Fig. 5 shows concentration profiles where we can see some plumes spread to more than 80% depth of the reservoir just within 20 years and by 500 years the deposition reaches more than 43. The results of 0.35 heterogeneity in the same time fluctuates by more than 4.5 times less. Interestingly at lower times the differences are more severe because when heterogeneity is low convection occurs by mostly gravity fingering. However, in the case of high heterogeneity dispersive convection tend to appear at short time period and as time passes fluids get motion for all cases. As reactions may occur in different orders, we investigate relative influences of so by calculating depositions. The orders can vary based on the type of reacting minerals, concentration of minerals as well as of dissolved CO2, etc. For instance, the precipitation of calcite is projected to be a first order reaction against the complicated carbonate mineral trapping process in saline aquifers [32]. On the other hand, if the concentration of minerals in the medium is much higher than the dissolved CO2 then reactions may not occur in stoichiometric ratios. This is also the case of first order reaction. It is noteworthy that non-stoichiometric ratios can be related to stoichiometric weight through the definition of Da [8].

Fig. 6 shows deposition results of different reaction orders. Results presented are for Da ffi 100. Because mineral concentrations are considered to be too high in the case of first order reaction the deposition of CO2 is higher than other orders. The heavy mineral interactions will cause to consume much CO2 convected from early times. The second order reaction where the stoichiometric ratio, 1:1, is the slowest process for the long run. Mineral interactions are not much like in the case of first order and concentration changes will be the same always (i.e., exponents of c and cm are 1). We have assumed third order reaction such that with each mole of CO2 two moles of formation rocks react. Because of higher order reaction the rate will increase with time and subsequently the consumption will be more pronounced. Our results show that until 100 years of deposition owing to third order reaction is slower than other two orders. Though scope of interaction is high but there is not much CO2 available in early time to react with. As convection increases the reaction gets more dynamicity to consume more CO2 and thereby starts enhancing the deposition. After 1000 years this order of reaction becomes dominating in settling down the CO2. The relative effects of reaction orders can be more or less based on the concentration difference (Dc) imposed initially (t = 0). We have tested one simulation case (Da = 100) assuming second-order reaction for a large number of time steps to trace dissolution of longer time periods. Fig. 7 shows distributions of CO2-brine concentration, minerals concentration, and the temperatures profile. As seen, after couple of hundreds of years advancement of concentration fingers follow dispersive convection which in turn results in slow movement of the plumes, however wetting large volume of aquifer spaces. This makes overall deposition favorable. Fig. 7c shows minerals distribution which is in consistent with concentration maps of Fig. 7b. The highest activity of CO2-brine at top of the reservoir reduces minerals density to almost zero and where the convection cells have not reached much shows maximum rock formation being unreacted. Because geothermal gradients play very minor role in convection process the temperature distributions remain nearly unchanged. In order to get an idea of after how long time complete deposition can be achieved we extrapolated the numerical simulation values and have found that at t  11,000 years the CO2 saturated brine will be 100% settled down. This results shown in Fig. 8 represents for a typical geologic formation (Ras = 1000) with moderate heterogeneity (Vdp = 0.55).

A. Islam et al. / International Journal of Heat and Mass Transfer 77 (2014) 519–528

527

Fig. 8. Deposition over long periods (Ras = 1000, N = 2, Vdp = 0.55).

5. Conclusions Density driven natural convection in a hypothetical heterogeneous two dimensional geothermal reservoir subject to vertical concentration and temperature gradients, geochemical interactions between CO2 saturated brine and formation rocks, was numerically investigated. We studied propagation of the plumes and their resulting dissolution into brine in term of percentage of deposition over long period of times. For investigating impacts of geochemical reactions on deposition processes we covered the results of wide range of Da and on top of that influence of reaction orders (up to the 3rd order) were analyzed. In practice, reaction formula can be complicated; however we have shown in a simple form (Eq. (6)) representing an average of possible reactions. As a general consequence geochemical reactions boost overall deposition process because mineral trapping leads to stratification of density gradients imposing more CO2 to convect. The intensity of enhancement depends mainly on the order of reaction rate, in another word, on Da based on level of heterogeneity and Ras considered. In summary, following conclusions can be drawn:

Fig. 7. Contours of CO2-brine concentration at (a) 500 and (b) 2500 years; (c) minerals concentration at 2500 years; (d) temperature distributions at 2500 years (Ras = 1000, N = 2, Vdp = 0.55).

At low Da (<1000) the geochemical interactions improve deposition by maximum 2%. When Da > 1000 the deposition can be increased by more than 5% by 500 years’ period. This is significant because dissolved CO2 has trapping time ranging from thousands to hundreds of thousands of years. The influence of mineral trapping is increased with increasing Ras. As the heterogeneity increases geochemical reactions play more role of settling down the CO2. At very high Vdp (=0.85) in moderate Ras (=1000) the simultaneous dissolution and mineral trapping can even reach to 50% by only 500 years. Because mineral concentrations is assumed to be too high in the case of 1st order reaction, so results in higher deposition than other orders beginning from early times. In the case of 3rd order reaction the impact increases with time and after 1000 years it will dominate the depositions. The 2nd order reaction will perform the slowest trapping after 100 years. For a tested case (Da = 100, 2nd order, Ras = 1000, Vdp = 0.55, N = 2) it is found by extrapolating of simulation results that complete saturation can be reached to by 11,000 years. In any case the geothermal gradients play minor effect in depositing the CO2.

528

A. Islam et al. / International Journal of Heat and Mass Transfer 77 (2014) 519–528

Conflict of interest The authors declare they have no conflict of interest with any person or organization. Acknowledgements The work was conducted with the support of the Reservoir Simulation Joint Industry Project, a consortium of operating and service companies for Petroleum and Geosystems Engineering at The University of Texas at Austin. References [1] S. Bachu, CO2 storage in geological media: role, means, status and barriers to deployment, Prog. Energy Combust. Sci. 34 (2008) 254–273. [2] Y. Cinar, A. Riaz, H.A. Tchelepi, Experimental study of CO2 injection into saline formations, SPE J. Dec. (2009) 588–594. [3] M.L. Druckenmiller, M. Maroto-Valer, M. Hill, Investigation of carbon sequestration via induced calcite formation in natural gas well brine, Energy Fuels 20 (2006) 172–179. [4] R. Farajzadeh, P. Ranganathan, P.L.J. Zitha, J. Bruining, The effect of heterogeneity on the character of density-driven natural convection of CO2 overlying a brine aquifer, Adv. Water Res. 34 (2011) 327–339. [5] R. Farajzadeh, H. Salimi, P.L.J. Zitha, J. Bruining, Numerical simulation of density-driven natural convection in porous media with application for CO2 injected projects, Int. J. Heat Mass Transfer 50 (2007) 5054–5064. [6] S. Fischer, A. Liebscher, M. Wandrey, CO2–brine–rock interaction – first results of long-term exposure experiments at in situ P–T conditions of the Ketzin CO2 reservoir, Chem. Erde 70 (2010) 155–164. [7] I. Gaus, P. Audigane, L. Andre, J. Lions, N. Jacquement, P. Durst, Isabel Czernichowski-Laurel, M. Azaroual, Geochemical and solute transport modeling for CO2 storage, what to expect from it?, Int J. Greenhouse Gas Control 2 (2008) 605–625. [8] K. Ghesmat, H. Hassanzadeh, J. Abedi, The impact of geochemistry on convective mixing in a gravitationally unstable diffusive boundary layer in porous media: CO2 storage in saline aquifers, J. Fluid Mech. 673 (2011) 480– 512. [9] L.L. Green, T.D. Foster, Secondary convection in a Hele Shaw cell, J. Fluid Mech. 71 (1975) 675–687. [10] Z. Gu, A geological compositional simulator for modeling CO2 sequestration in geological formations (PhD dissertation), The University of Utah, 2010, 1980. [11] R.D. Gunn, W.B. Krantz, Reverse combustion instabilities in tar sands and coal, Soc. Petrol. Eng. J. 20 (1980) 267–277. [12] H. Hassanzadeh, M. Pooladi-Darvish, D.W. Keith, Modeling of convective mixing in CO2 storage, J. Can. Pet. Technol. 44 (2005) 43–51. [13] A.W. Islam, M.A.R. Sharif, E.S. Carlson, Numerical investigation of double diffusive convection of carbon dioxide in a brine saturated geothermal reservoir, Geothermics 48 (2013) 101–111. [14] A.W. Islam, H.R. Lashgari, Kamy. Sepehrnoori, Double diffusive natural convection of CO2 in brine saturated geothermal reservoir: study of nonmodal growth of perturbations and heterogeneity effects, Geothermics 51 (2014) 325–336. [15] O. Izgec, B. Demiral, H. Bertin, S. Akin, Experimental and numerical investigation of carbon sequestration in saline aquifers, Paper SPE 94697STU, 2005. [16] M. Javaheri, J. Abedi, H. Hassanzadeh, Onset of convection in CO2 sequestration in deep inclined saline aquifers, J. Can. Pet. Technol. 48 (2009) 22–27. [17] M. Javaheri, J. Abedi, H. Hassanzadeh, Linear stability analysis of doublediffusive convection in porous media, with application to geological storage of CO2, Transp. Porous Media 84 (2010) 441–456.

[18] R. Juanes, C.W. MacMinn, Upscaling of capillary trapping under gravity override: application to CO2 sequestration in aquifers, Paper SPE 113496, 2008. [19] L.D. Landau, E.M. Lifshitz, Statistical Physics: Course of Theoretical Physics, Pergamon Press, 1969. [20] Q. Li, Z.S. Wu, X.C. Li, Prediction of CO2 leakage during sequestration into marine sedimentary strata, Energy Convers. Manage. 50 (2009) 503–509. [21] P. Mandalaparty, Reaction chemistry in carbon dioxide sequestration (PhD dissertation), The University of Utah, 2012. [22] D.A.C. Manning, Biological enhancement of soil carbonate precipitation: passive removal of atmospheric CO2, Mineral Mag. 72 (2008) 639–649. [23] D.L. Newell, H.S. Viswanathan, P.C. Litchner, J.W. Carey, J.P. Kaszuba, Experimental determination of supercritical CO2 mass transfer rates into brine, in: American Geophysical Union Fall Meeting, Sanfrancisco, 2011. [24] L. Nghiem, V. Shrivastava, D. Tran, B. Kohse, M. Hassan, C. Yang, Simulation of CO2 storage in saline aquifers, in: SPE 125848, 2009. [25] M. Noh, L.W. Lake, S.L. Bryant, A. Araque-Martinez, Implication of coupling fractional flow and geochemistry for CO2 injection in aquifers, in: SPE 89341, 2004. [26] C.M. Oldenburg, K. Pruess, Layered thermohaline convection in hypersaline geothermal systems, Transp. Porous Media 33 (1998) 29–63. [27] W.K. O’Connor, D.C. Dahlin, D.N. Nielsen, Continuous studies on direct aqueous mineral carbonation for CO2 sequestration, in: 27th international conference on coal utilization & fuel systems, clear water, 2002. [28] I. Okamoto, S. Mito, T. Ohsumi, A sensitivity study of CO2 mineralization using GEM-GHG simulator, Energy Procedia 1 (2009) 3323–3329. [29] L. Ostrowski, B. Ulker, Minimizing risk of gas escape in gas storage by in-situ measurement of gas threshold pressure and optimized completion solutions, in: SPE 113509, 2008. [30] K. Pruess, T. Xu, J. Apps, J. Garcia, Numerical modeling of aquifer disposal of CO2, SPE J. 8 (2003) 49–60. [31] R.J. Rosenbauer, T. Koksalan, J.L. Palandri, Experimental investigation of CO2– brine rock interactions at elevated temperature and pressure: implications for CO2 sequestration in deep-saline aquifers, Fuel Process. Technol. 86 (2005) 1581–1597. [32] Q. Li, Coupled reactive transport model for heat and density driven flow in CO2 storage in saline aquifers, J. Hazard Toxic Radioact. Waste 15 (2011) 251–258. [33] K.P. Saripalli, M.M. Sharma, S.L. Bryant, Modeling the injection well performance during deep-well injection of liquid wastes, J. Hydrol. 227 (2000) 41–55. [34] Schumacher, M. Anne, Modeling of CO2–water–rock interactions in Mississippian sandstone reservoir of Kentucky (Master theses), University of Kentucky, 2013. [35] M.H. Sharqawy, V.J.H. Lienhard, S.M. Zubair, Thermophysical properties of seawater: a review of existing correlations and data, Desalin. Water Treat. 16 (2010) 354–380. [36] M.J. Simpson, T.P. Clement, Theoretical analysis of the worthiness of Henry and Elder problems as benchmarks of density-dependent groundwater flow models, Adv. Water Resour. 26 (2003) 17–31. [37] J.R. Waggoner, J.L. Castilo, L.W. Lake, Simulation of EOR process in stochastically generated permeable media, SPE Form. Eval. 7 (1992) 173–180. [38] T.P. Wellman, R.B. Grigg, B.J. McPherson, R.K. Svec, P.C. Lichtner, Evaluation of CO2–brine–reservoir rock interaction with laboratory flow tests and reactive transport modeling, in: Paper SPE 80228, 2003. [39] G.H. Wolf, A.V.G. Chizmeshya, J. DiefenBacher, M.J. McKelvy, Insitu observation of CO2 sequestration reactions using a novel micro-reaction system, Environ. Sci. Technol. 38 (2004) 932–936. [40] X. Xu, S. Chen, D. Zhang, Convective stability analysis of the long-term storage of carbon dioxide in deep saline aquifers, Adv. Water Resour. 29 (2006) 397– 407. [41] T. Xu, J.A. Apps, K. Pruess, Mineral sequestration of carbon dioxide in a sandstone-shale system, Chem. Geol. 217 (2005) 295–318. [42] B. Zerai, B. Saylor, G. Matisof, Computer simulation of CO2 trapped through mineral precipitation in the Rose-Run sandstone, Ohio, Appl. Geochem. 21 (2006) 223–240.