A fluid pressure and deformation analysis for geological sequestration of carbon dioxide

Computers & Geosciences 46 (2012) 31–37

Contents lists available at SciVerse ScienceDirect

Computers & Geosciences journal homepage: www.elsevier.com/locate/cageo

A fluid pressure and deformation analysis for geological sequestration of carbon dioxide Zhijie Xu a,n, Yilin Fang b, Timothy. D. Scheibe b, Alain Bonneville c a

Computational Mathematics Group, Fundamental and Computational Sciences Directorate, Pacific Northwest National Laboratory, Richland, WA 99352, USA Hydrology Technical Group, Energy and Environment Directorate, Pacific Northwest National Laboratory, Richland, WA 99352, USA c Earth System Science Division, Energy and Environment Directorate, Pacific Northwest National Laboratory, Richland, WA 99352, USA b

a r t i c l e i n f o

abstract

Article history: Received 23 January 2012 Received in revised form 21 April 2012 Accepted 21 April 2012 Available online 30 April 2012

We present a fluid pressure and deformation analysis for geological sequestration of carbon dioxide based on a simplified hydro-mechanical model. This model includes the geomechanical part that relies on the theory of linear elasticity, while the fluid flow is based on the Darcy’s law. Two parts are coupled together using the linear poroelasticity theory. For a typical geological sequestration in a semi-infinite geometry with diminishing pressure and deformation fields at infinity, the Helmholtz decomposition can be applied to the displacement vector. Hence, the flow equation can be decoupled from the equation of linear elasticity. Solutions for fluid pressure were obtained for this typical scenario and solutions for ground deformation were obtained using the method of Green’s function. Finally, solutions were compared against numerical results using a finite element method for aquifers with two different thicknesses. General agreement can be obtained between analytical and numerical solutions. The model is useful in estimating the temporal and spatial variation of fluid pressure and the mechanical deformation during the entire injection period. & 2012 Elsevier Ltd. All rights reserved.

Keywords: Geological sequestration Geomechanics Poro-elasticity Hydro-mechanical

1. Introduction The purpose of this paper is to present an analysis of the geomechanical response subjected to the injection of substantial carbon dioxide into a deep geological formation. It is well known that carbon dioxide sequestration in deep saline aquifers could be an effective mitigation strategy for the reduction of CO2 emissions into the environment. In order to achieve a substantial reduction of CO2 emissions and significantly alleviate the atmospheric carbon problem, massive scale carbon sequestration is required involving millions of tons of carbon dioxide injection over several decades. Carbon dioxide in the supercritical state enables efficient geological storage due to its high density relative to the gas state. However, supercritical CO2 is less dense than brine formation waters, so the buoyancy driven vertical flow of supercritical CO2 requires that the aquifer be capped by a layer of low permeability rock (‘‘caprock’’) to prevent the vertical migration of CO2 to the ground surface (shown in Fig. 1). Injection of CO2 into the deep saline aquifer changes the entire stress state and leads to the pressure buildup under the caprock. The changes in the effective stress and the induced deformation could potentially damage the caprock and create new fractures or reactivate existing sealed

n

Corresponding author. Tel.: þ1 509372 4885; fax: þ1 509 372 4720. E-mail address: [email protected] (Z. Xu).

0098-3004/$ - see front matter & 2012 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.cageo.2012.04.020

fractures. This fracturing process could open new flow paths for supercritical CO2 and substantially reduce the sequestration effectiveness. Therefore, a coupled hydro-mechanical model is essential for the caprock integrity analysis and the risk analysis associated with the potential CO2 leakage. Rigorous analysis of CO2 geological sequestration almost always relies on computational models that provide mathematical descriptions for the most relevant physical and chemical processes (Celia and Nordbotten, 2010). For example, the combined fluid flow and transport process and geochemical reactions have important roles for carbon sequestration. CO2 can react with minerals to form carbonates as part of a ‘‘weathering’’ process that occurs naturally, is thermodynamically favorable, and can be mathematically described by the ‘‘solute precipitation and/or dissolution’’ process (Xu and Meakin, 2008, 2011; Xu et al., 2012) and the ‘‘corrosion process’’ (Xu et al., 2011) that have been studied previously. Some of the most relevant processes for carbon sequestration include multi-phase flow, multi-component transport, geochemical reaction, geomechanical response, and heat transport for non-isothermal problems. The multiple time and space scales associated with each relevant process significantly complicate the problem. A fully coupled three-dimensional multi-phase, multi-physics, multicomponent computational model, though provides a relatively more complete description of the problem, is sometimes fragile due to the difficulties in model parameter identification, process

32

Z. Xu et al. / Computers & Geosciences 46 (2012) 31–37

z o

Caprock Aquifer

r

Z2 = -1400m Z1 = -1500m

Fig. 1. A schematic plot of a typical geological sequestration of carbon dioxide.

coupling, severe heterogeneity and nonlinearity, and many other computational challenges. Effective models based on homogenization and upscaling can be obtained for simple physics (Xu, 2012). However, simplifications to the fully coupled three-dimensional models for carbon sequestration are almost always necessary for a practically tractable problem (Celia and Nordbotten, 2009). For example, the similarity solutions developed by Nordbotten and Celia (Nordbotten and Celia, 2006) focus on the multi-phase flow process and prediction of the CO2 migration during the entire injection period. This paper focuses on the development of solutions for the geomechanical response to large scale CO2 injection. Accordingly, some simplifications are necessary in order to obtain analytical solutions. In the current model, the heat transport equation is eliminated by assuming that the injected CO2 is at a temperature close to that of the formation. Geochemical responses are assumed to take much longer time scales in contrast to the geomechanical responses (Celia and Nordbotten, 2010). Geochemical reactions between CO2 and caprock mineralogy could happen from several years to hundred years. The porosity change due to the mineral dissolution and/or precipitation can affect the fluid flow and hence the geomechanical response (Gaus et al., 2005). This effect is not considered in the current model. Moreover, the two-way coupling between fluid flow and geomechanics was simplified by assuming the injection induced stress and deformation do not significantly change the fluid flow properties (Celia and Nordbotten, 2010). We do not distinguish two fluid phases, namely the brine water phase and the supercritical CO2 phase, in order to study the geomechanical response (stress and deformation) that largely depends on the fluid pressure distribution. The original multiphase flow problem was treated as a single phase flow with effective fluid properties that obeys the Darcy’s law in general. Therefore, a single pressure field p was used in both CO2 and brine water phases by neglecting the capillary pressure pc between two phases, which could be small compared to the total injection pressure. In principle, the simplified problem can be potentially applied to many other problems in various areas, for example, the geomechanical effect due to the fluid migration and heat transport in oil/gas field and geothermal reservoirs that has been recently studied by Chen (Chen, 2011), and subsidence due to groundwater withdraw that can be best understood with the consolidation theory (Gambolati and Freeze, 1973; Gambolati et al., 1974). Another similar problem is the deformation due to the eruption of volcanoes that was modeled by the change of hydrodynamic pressure in a small sphere (magma reservoir) in the semi-infinite elastic solid (Mogi, 1958). The purpose of this paper is to establish pressure and deformation solutions of typical geometry and injection parameters with a focus on the geological sequestration of carbon dioxide.

As the first step in the paper, we will focus on calculating the deformation with reasonable injection parameters and typical geometry. A comparison with full multiphase flow model could be in place in order to show the difference. The paper is organized as follows. We first present a hydro-mechanical model in Section 2, followed by the solutions for fluid pressure and displacement in Section 3, the numerical comparison in Section 4, and the conclusion in Section 5.

2. A hydro-mechanical model for CO2 geological sequestration The hydro-mechanical model includes both fluid flow and the elasticity equations that describe the geomechanical response to the fluid injection and flow. The system of equations for the hydro-mechanical model first includes the flow equations of Darcy’s law in terms of the total pore fluid pressure ptot:   @p k c yb tot ¼ rU ðDptot rgÞ þ , ð1Þ @t m r where ptot(x,t)¼pini þ Dp is the pressure field at position x and time t and r is the Laplace operator. Eq. (1) describes the Darcy flow for single phase considering that capillary pressure is small compared to the injection pressure. pini is the initial (reference) pressure field and Dp represents the pressure change at position x and time t before and after the injection. y is the porosity of the aquifer. b, m, and r are all fluid and flow properties representing the fluid compressibility, viscosity, and density. k is the aquifer permeability and c is the external source term with a unit of kg/ (m3 s). Since By substituting expression pini ¼ rgz (the hydrostatic pressure, where z is the coordinate in depth direction) into Eq. (1), we have the equation for Dp (the pressure change, we use p¼ Dp through the end for a simpler representation),   @p k c yb ¼ rU rp þ : ð2Þ @t m r When it comes to the coupled hydro-mechanical model, the earliest theory was developed in 1923 by Terzaghi (Terzaghi, 1923) to account for the effect of pore fluid on the soil consolidation, followed by the linear theory of poroelasticity developed by Biot in 1930’s (Biot, 1935, 1941), and further improved in 1950’s (Biot, 1955, 1956, 1962), @p @ðrUuÞ k cM þ aB M ¼ Mr2 p þ , @t @t m r

ð3Þ

where aB ¼1  (K/Ks) is the Biot’s parameter and M ¼(BKu/aB) is the Biot modulus. K is the drained bulk modulus, Ks is the bulk modulus of solid constituent (Matrix), and Ku is the undrained bulk modulus. B is the Skempton pore pressure coefficient. In the limiting case where the bulk modulus K is negligible compared to that of solid constituent Ks (namely, K 5Ks, aB ¼ 1 and M¼1/(yb) in Eq. (3)), the coupled hydro-mechanical model reads @p 1 @ðrUuÞ k c þ ¼ r2 p þ , @t yb @t ymb yrb

ð4Þ

and ðl þ GÞrðrUuÞ þ Gr2 u ¼ rp:

ð5Þ

Eq. (4) is a Darcy’s flow equation in terms of the pressure field p and Eq. (5) is the displacement formulation (the Navier–Lame´ equation) of elasticity, where u is the displacement vector. The two-way coupling between fluid flow and geomechanics was implemented through the divergence term rUu (volumetric strain) in Eq. (3) and the gradient term rp acting as a source in Eq. (5). G is the shear modulus and l is Lame’s constant in the theory of elasticity, both of which can be related to the Young’s

Z. Xu et al. / Computers & Geosciences 46 (2012) 31–37

modulus and the Poisson’s ratio by Gð3l þ2GÞ l and n ¼ : E¼ l þG 2ðl þ GÞ

ð6Þ

Eqs. (4) and (5) are valid for arbitrary geometry and boundary conditions. For practical CO2 geological sequestration, a typical geometry considered here is shown in Fig. 1, where CO2 is injected into a confined aquifer formation with a position from z1 to z2 in depth at a constant injection rate c. The pressure dissipation into the caprock and bottom in z-direction would be negligible due to very low permeability. The stress boundary conditions for this particular geometry are

szz ðz ¼ 0Þ ¼ sxz ðz ¼ 0Þ ¼ syz ðz ¼ 0Þ ¼ 0,

ð7Þ

where sij is the ij component of stress tensor. For a homogeneous, isotropic, linear elastic solid, sij is related to the strain component eij through

sij ¼ lekk dij þ 2Geij ,

ð8Þ

where dij is Kronecker delta and ekk follows the standard Einstein summation. The strain component eij is further related to the displacement vector u through the compatibility equation for small deformation, where

eij ¼ ð@i uj þ@j ui Þ=2:

ð9Þ

Eqs. (4) and (5) together with the boundary condition (7)–(9) provide a complete description of the hydro-mechanical model for the geometry in Fig. 1 that is solvable by any standard numerical method. In order to obtain insights on the fluid pressure and deformation fields during the entire injection period, it is important and possible to find the analytical solutions for this typical geometry.

3. Solutions for a typical CO2 geological sequestration To find the analytical solution, we first make use the identity

r2 u ¼ rðrUuÞr  ðr  uÞ,

ð10Þ

and rewrite the Navier–Lame´ Eq. (5) as ðl þ2GÞrðrUuÞGr  ðr  uÞ ¼ rp,

ð11Þ

where r  is the curl operation on the vector field. Next, the displacement u is a smooth and decaying vector field that is vanishing at infinity. It can be decomposed into u ¼ rf þ r  u

ð12Þ

using the Helmholtz theorem, where f is a scalar field and u is a divergence free vector field satisfying rUu ¼0. By this way, the displacement field u is split into an irrotational field and a rotational field. Substituting of Eq. (12) into Eqs. (3) and (11), and making use the identities rU(r  u)¼0, r  (ru)¼0, r  (r  u)¼ r(rUu)  r2u, and r  r(rUu)¼0 lead to the following equations 2

@p @ðr fÞ k c þ aB M ¼ M r2 p þ M, @t @t m r

ð13Þ

and 2

2

ðl þ2GÞrðr fÞ þ Gr  r u ¼ rp:

ð14Þ

Applying the divergence operator (rU) on both sides of Eq. (14) and using the identity rU(r  r2u)¼0 yield the relationship 2

2

2

ðl þ2GÞr ðr fÞ ¼ r p:

ð15Þ

The integration of Eq. (15) leads to

r2 f ¼

1

l þ 2G

p þ A,

ð16Þ

33

where A is an integration constant. Since the displacement vector u and pressure p vanishes at infinity (u9N ¼9pN ¼0) for this particular geometry, we have the constant A¼ r2f9N ¼(rUu)9N ¼0. Substitution of Eq. (16) back into the flow Eq. (13) leads to the decoupled flow equation in terms of the fluid pressure p, @p k=m c=r ¼ : r2 p þ @t 1=M þ aB =ðl þ 2GÞ 1=M þ aB =ðl þ2GÞ

ð17Þ

For the geometry shown in Fig. 1 with a constant line injection rate at the centerline (r ¼0) across the entire aquifer, Eq. (17) can be written in the polar coordinate system (there is no pressure variation along the z-direction due to the constant line injection), where r is the radial coordinate, ! @p @2 p 1 @p ¼ Dr 2 p ¼ D þ : ð18Þ @t r @r @r 2 An equivalent diffusion coefficient D was defined to describe the pressure diffusion, D¼

k=m : 1=M þ aB =ðl þ2GÞ

ð19Þ

It is clear that the effects of fluid and flow properties (permeability, viscosity) and solid mechanical properties (modulus) on pressure diffusion can be lumped into a single parameter D. Large modulus and permeability or small viscosity leads to faster pressure diffusion. A similar expression has been derived for D oedometric depletion, where D ¼ ððk=mÞ=ð1=MÞ þ ða2B =ðl þ 2GÞÞÞ (Charlez, 1997). More specifically, two results are exactly the same for aB ¼1. Here the pressure equation is decoupled from the coupled model (Eqs. (3) or (4) and (5)) by use of Helmholtz decomposition of the displacement field. In this particular injection scenario, a line injection is assumed along the centerline (injection well radius is zero), and the boundary condition for pressure Eq. (18) for is written as   @p  c0 =r 2prD and p9r ¼ 1 ¼ 0, ¼ ð20Þ  @r r ¼ 0 yb þ 1=ðl þ 2GÞ where c0 is a line injection rate with a unit of kg/(m s). The solution of p due to the combined Eqs. (18) and (20) can be obtained as ~ ZÞ, p ¼ pFð

ð21Þ

where the dimensionless number Z is defined as Z ¼(r2/(4Dt)) and F(Z) is an exponential integral function defined as, Z 1 0 FðlÞ ¼ et =t 0 dt 0 : ð22Þ Z

The scaling factor (with a unit of pressure) p~ is defined as, 0

p~ ¼

1 cm : 4p rk

ð23Þ

An example calculation can be made with the values of variables given in Table 1 (in the limiting case aB ¼1, M¼(1/yb)) where an equivalent pressure diffusivity is computed as D ¼0.808 m2/s and the scaling factor p~ ¼ 4:642 MPa. The injection rate was selected to equal the CO2 production rate by a standard 1000-MW coal-fired power plant (Hitchon, 1996). Other parameters are selected according to CO2 properties. With these numbers, Fig. 2 presents the variation of pressure field p across a 40 km region during the entire injection period of two decades. The effect of aquifer permeability on the injection pressure is shown in Fig. 3 at a location of 5 km away from the injection well. ~ but A decrease in permeability k increases the scaling factor p, decreases the diffusivity D. The combined effects lead to an increase in the pressure p. However, a larger modulus means a

34

Z. Xu et al. / Computers & Geosciences 46 (2012) 31–37

20

Table 1 Values of variables used in the example calculations.

18 Property

Unit

Value

Lame’s constant, l Shear Modulus, G Injection rate, c0 Compressibility, b CO2 Density, r Porosity, y Fluid viscosity, m permeability, k Biot’s parameter, aB

GPa GPa kg/(m s) 1/Pa kg/m3 Dimensionless Pa s m2 Dimensionless

14 14 3.5 10  9 600 0.1 10 3 10 13 1

Pressure Pr = 5km (MPa)

16 14 12 10 8 λ = G = 14GPa

6

λ = G = 2GPa

4

λ = G = 50GPa

35

2 t = 5years

30

0

t = 10years

0

25 Pressure (MPa)

t = 20years

4

6

20 15

5 0 0

5

10

15

20 r (km)

25

30

35

40

Fig. 2. Spatial variation of fluid pressure after 5 years, 10 years, 15 years and 20 years since injection.

30 k = 1x10-13m2

k=

14

16

18

20

With the solution for the pressure field p (Eq. (21)), now let us return to seek displacement solutions for Eq. (5) combined with the boundary conditions (7). The more fundamental problem behind Eq. (5) is to find the elastic field of a center of dilation (due to the pressure p at that center) in a semi-infinite solid. The method of Green’s function can then be applied to find the elastic solution for any arbitrary pressure field p(x,t) via superimposition of the elemental solution for each center of dilation. The solution to the elastic field of a center of dilation was first obtained by Mindlin and Cheng (Mindlin and Cheng, 1950a, b) using the method of image and later by Sen (Sen, 1951) for isotropic elastic semi-infinite solid. Mindlin and Cheng expressed the solution in terms of the potential functions of center of dilation and its image in the infinite solid. More specifically (as shown in Fig. 5), the displacement field u at point A (positioned by vector R(x,y,z)) due to a center of dilation with a position vector R0 (x0 ,y0 ,z0 ) can be written as 2 ~ ~ þ2z @ c ! i uðRÞ ¼ rc þð34nÞrc @x@z " # ~ ! ~ ~ ! @2 c @2 c @c þ 2z j þ 2z 2 2ð34nÞ k, @y@z @z @z

k = 2x10-13m2

25

8 10 12 Time (Years)

Fig. 4. Time variation of fluid pressure at 5 km away from the injection point for aquifers with different elastic modulus.

10

Pressure Pr = 5km (MPa)

2

t = 15years

0.5x10-13m2

20

ð24Þ

! ! ! where i , j , k are unit vectors in Cartesian coordinate system. c and c~ are two harmonica functions satisfying the Laplace equation. c is the potential function for the center of dilation,

15

10

cðRÞ ¼ 

1 pðR0 ,tÞdR03 , 4p ðl þ2GÞr0

ð25Þ

~ is the potential function for the image of the center of and c dilation,

5

0 0

2

4

6

8 10 12 Time (Years)

14

16

18

20

Fig. 3. Time variation of fluid pressure at 5 km away from the injection point for aquifers with different permeability.

~ and eventually leads to larger diffusivity D, but has no effect on p, a decrease in the pressure p at given location and time, as shown in Fig. 4. The solution (21) is the transient Theis solution (Theis, 1935) incorporating the mechanical effect.

c~ ðRÞ ¼ 

1 pðR00 ,tÞdR003 : 4p ðl þ 2GÞr00

ð26Þ

The position of center and its image is symmetric about the x–y plane, so that R0 ðx0 ,y0 ,z0 Þ ¼ R00 ðx0 ,y0 ,z0 Þ:

ð27Þ

In cylindrical coordinate system (r,g,z), the distance between point A and the center of dilation is qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ð28Þ r0 ðr,zÞ ¼ 9RR0 9 ¼ ðzz0 Þ2 þ r2 þ r20 2rr 0 cos g,

Z. Xu et al. / Computers & Geosciences 46 (2012) 31–37

Image of center of dilation

R’’

35

0.08 t = 1 year t = 5 years

z

uz (z = 0) (m)

0.06

ρ'’

x

0.04

0.02

A R

r

0.00 0

ρ' R’

r0

t = 10 years

10

20 r (km)

30

40

Fig. 6. Spatial variation of ground uplift after 1 year, 5 years, and 10 years since injection.

Center of dilation

Fig. 5. Time variation of fluid pressure at 5 km away from the injection point for aquifers with different elastic modulus.

0.08

ð29Þ

where r and r0 are the radial coordinate for point A and the dilation center, respectively. g is the difference of azimuthal angle between point A and the dilation center. Mathematical analysis of subsidence behavior above oil reservoir has been implemented by Geertsma (Geertsma, 1973) based on the concept of nucleus-of-strain introduce by Mindlin and Chen (Mindlin and Cheng, 1950a). Hereby we are following the similar procedure. After substitution of expressions (25) and (26) into the displacement solution (24) and applying the method of Green’s function, the vertical displacement solution uz is: Z 1 Z p Z z2 1 uz ¼ 2pðl þ 2GÞ 0 0 z1 " # 0 zz ð14nÞz þ ð34nÞz0 6zðz þ z0 Þ2 pðr 0 ,tÞ   ð30Þ r 0 dz0 dg dr 0 , 005 03 003

r

0.07

0.06

0.05

0.04 0

2

4

6 time (years)

8

10

Fig. 7. Time variation of ground uplift at injection point since injection.

r

r

and the radial displacement ur is: Z 1 Z p Z z2 1 ur ¼ 2pðl þ 2GÞ 0 0 z1   1 ð34nÞ 6zðz þz0 Þ pðr 0 ,tÞðrr 0 cos gÞ 03 þ  r 0 dz0 dg dr 0 : 005 003

r

uz (z = 0, r = 0) (m)

and the distance between point A and the image is qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi r00 ðr,zÞ ¼ 9RR00 9 ¼ ðz þ z0 Þ2 þ r2 þ r20 2rr 0 cos g,

r

r

More specifically, the free surface uplift is: Z 1 Z p Z z2 1 pðr 0 ,tÞz0 uz 9z ¼ 0 ¼ r 0 dz0 dg dr 0 : pðl þGÞ 0 0 z1 r03

the radial displacement at the free surface reads, Z 1 Z p Z z2 1 pðr 0 ,tÞðrr 0 cos gÞ ur 9z ¼ 0 ¼ r 0 dz0 dg dr 0 : pðl þ GÞ 0 0 z1 r03 ð31Þ

ð34Þ

Fig. 8 plots the variation of the free surface radial displacement at a time of 1 year, 5 years, and 10 years since injection from Eq. (34), where the radial displacement is zero at the injection point, and gradually increases to a maximum value, followed by a decrease to zero at r ¼N.

ð32Þ

By inserting the pressure solution (Eq. (21)) into Eq. (32), Fig. 6 shows the spatial variation of free surface uplift at time of 1 year, 5 years and 10 years after injection from Eq. (32). At the injection point (point o in Fig. 1), the ground uplift reads, 0 1 Z 1 1 1 1 B C pðr 0 ,tÞr 0 @qffiffiffiffiffiffiffiffiffiffiffiffiffiffi  qffiffiffiffiffiffiffiffiffiffiffiffiffiffiAdr 0 : ð33Þ uz 9z ¼ 0,r ¼ 0 ¼ lþG 0 2 2 2 2 z2 þ r 0 z1 þ r 0 Fig. 7 presents the variation of the ground uplift at the injection point o during the first 10 years since injection based on Eq. (33), with an annual increase of 2.5 mm in the last four years. Similarly,

4. Comparison with FEM solutions To test the validity and assess the accuracy of the solution, coupled model Eqs. (4) and (5) were numerically solved by a finite element method (FEM) to obtain the numerical solutions of the fluid pressure and the deformation field for the same geometry and injection parameters. The FEM model was set with zmin ¼ 2000 m and r max ¼ 100 km. Fig. 9 shows the model setup and a typical FEM mesh consisting of 24,000 triangular elements with adaptive meth refinement very close to the injection well. Sizes of triangular elements range from less than 1 m close to the injection well up to 1000 m. The aquifer layer of thickness

Z. Xu et al. / Computers & Geosciences 46 (2012) 31–37

H¼z2  z1 was used for the geological sequestration of carbon dioxide. Two thicknesses have been studied, namely H¼100 m (z1 ¼ 1500 m and z2 ¼  1400 m) and H¼200 m (z1 ¼ 1500 m, z2 ¼  1300 m). The comparison of analytical solutions against the FEM solutions for the pressure field and the displacement field after one year is presented in Figs. 10 and 11, respectively. In general, both pressure and displacement are shown in good agreement with numerical solutions. In order to check the accuracy after long period of injection and the effect of aquifer thickness, we compare the maximum ground uplift uz(r¼0,z¼0) against the FEM solutions during the first 10 years of injection. The results are presented in Fig. 12, where solid lines represent the analytical solutions from Eq. (33). For both aquifer thicknesses, the analytical solutions are in agreement with the numerical solutions. After 10 years of injection, the analytical solutions overestimate the deformation by 2.5% for H¼100 m and around 5% for H¼200 m. Analytical solutions provide a better estimation of the ground deformation for deeper aquifers, where the center of dilation used in the derivation is far from the ground surface.

decrease in fluid pressure. The displacement solution can be used to estimate the ground deformation during the entire injection period, and the effects of aquifer permeability and modulus on the ground deformation. Ongoing studies include the comparison of the current model with field-measured ground uplift data.

20

10

5

5. Conclusion

0

Solutions for fluid pressure and displacement fields were developed for a typical geological sequestration example based on a hydro-mechanical model. It was shown that an increase in aquifer permeability or decrease in elastic modulus leads to a

0

5

10

15 r (km)

20

25

30

Fig. 10. Comparison of spatial variation of fluid pressure after 1 year with the numerical solution. Solid line presents solutions from Eq. (21).

0.05

0.025 t = 1 year t = 5 years t = 10 years

Eq .(32) at t = 1 Year FEM solution at t = 1 Year

0.04

uz (z = 0) (m)

0.020

ur (z = 0) (m)

FEM solution at t = 1 Year Eq. (21) at t = 1 Year

15 Pressure (MPa)

36

0.015

0.010

0.03

0.02

0.01 0.005

0.00 0.000 0

10

20 r (km)

30

40

Fig. 8. Spatial variation of ground radial deformation after 1 year, 5 years, and 10 years since injection.

0

10

20 r (km)

40

Fig. 11. Comparison of spatial variation of ground uplift after 1 year with the numerical solution. Solid line presents the solution from Eq. (32).

z o r

Z2 = -1400m or -1300m Z1 = -1500m Zmin = -2000m Eq.

30

100km

Fig. 9. Sketch of a Finite element model with a zoom of the mesh configuration close to the injection well.

Z. Xu et al. / Computers & Geosciences 46 (2012) 31–37

uz (z = 0, r = 0) (m)

0.08

0.07

0.06 Eq. (33) for H = 100m FEM for H = 100m Eq. (33) for H = 200m FEM for H = 200m

0.05

0.04 0

2

4

6 time (years)

8

10

Fig. 12. Comparison of the maximum ground uplift at the origin for the first 10 years of injection with the numerical solutions. Solid lines present the solutions from Eq. (33).

References Biot, M.A., 1935. Le proble me de la consolidation des matie res argileuses sous une charge. Annals of Society Of Science Bruxelles B55, 110–113. Biot, M.A., 1941. General theory of three-dimensional consolidation. Journal of Applied Physics 12, 155–164. Biot, M.A., 1955. Theory of elasticity and consolidation for a porous anisotropic solid. Journal of Applied Physics 26, 182–185. Biot, M.A., 1956. Thermoelasticity and irreversible thermodynamics. Journal of Applied Physics 27, 240–253. Biot, M.A., 1962. Mechanics of deformation and acoustic propagation in porous media. Journal of Applied Physics 33. 1482-&. Celia, M.A., Nordbotten, J.M., 2009. Practical modeling approaches for geological storage of carbon dioxide. Ground Water 47, 627–638. Celia, M.A., Nordbotten, J.M., 2010. How simple can we make models for CO2 injection migration, and leakage? Energy Procedia 4, 3857–3864. Charlez, A.P., 1997. uRock Mechanics: petrolem Applications. Editions Technip, Paris.

37

Chen, Z.R., 2011. Poroelastic model for induced stresses and deformations in hydrocarbon and geothermal reservoirs. Journal of Petroleum Science and Engineering 80, 41–52. Gambolati, G., Freeze, R.A., 1973. Mathematical simulation of the subsidence of Venice: 1. Theory. Water Resources Research 9, 721–733. Gambolati, G., Gatto, P., Freeze, R.A., 1974. Mathematical simulation of the subsidence of Venice. 2. Results. Water Resource 10, 563–577. Gaus, I., Azaroual, M., Czernichowski-Lauriol, I., 2005. Reactive transport modelling of the impact of CO2 injection on the clayey cap rock at Sleipner (North Sea). Chemical Geology 217, 319–337. Geertsma, J., 1973. Land subsidence above compacting oil and gas reservoirs. Journal of Petroleum Technology 25, 734–744. Hitchon, B., 1996. Aquifer Disposal of Carbon Dioxide. Geoscience Publishing Ltd, Sherwood Park, Alberta. Mindlin, R.D., Cheng, D.H., 1950a. Nuclei of strain in the semi-infinite solid. Journal of Applied Physics 21, 926–930. Mindlin, R.D., Cheng, D.H., 1950b. Thermoelastic stress in the semi-infinite solid. Journal of Applied Physics 21, 931–933. Mogi, K., 1958. Relations between the eruptions of various volcanoes and the deformations of the ground surfaces around them. Bulletin of the Earthquake Research Institute, University of Tokyo 36, 99–134. Nordbotten, J.M., Celia, M.A., 2006. Similarity solutions for fluid injection into confined aquifers. Journal of Fluid Mechanics 561, 307–327. Sen, B., 1951. Note on the stresses produced by nuclei of thermo-elastic strain in a semi-infinite elastic solid. Quarterly of Applied Mathematics 8, 365–369. Terzaghi, K., 1923. Die berechnung der durchlassigkeitsziffer des tones aus dem verlauf der hydrodynamischen spannungserscheinungen. Sitzungsberichte. Akadichte der Wissenschaften, Wien Mathematische-Naturwissenschaftliche Klasse, Abteilung IIa 132, 105–124. Theis, C.V., 1935. The relation between the lowering of the piezometric surface and the rate and duration of discharge of a well using ground water storage. Transactions-American Geophysical Union 16, 519–524. Xu, Z., Meakin, P., 2008. Phase-field modeling of solute precipitation and dissolution. Journal of Chemical Physics 129, 014705. Xu, Z., Meakin, P., 2011. Phase-field modeling of two-dimensional solute precipitation/dissolution: solid fingers and diffusion-limited precipitation. Journal of Chemical Physics 134, 044137. Xu, Z., Rosso, K.M., Bruemmer, S.M., 2011. A generalized mathematical framework for thermal oxidation kinetics. Journal of Chemical Physics 135, 024108. Xu, Z.J., 2012. Homogenization and upscaling for diffusion, heat conduction, and wave propagation in heterogeneous materials. Communications in Theoretical Physics 57, 348–354. Xu, Z.J., Huang, H., Li, X.Y., Meakin, P., 2012. Phase field and level set methods for modeling solute precipitation and/or dissolution. Computer Physics Communications 183, 15–19.