Flow path development in different CO2 storage reservoir scenarios

Engineering Geology 127 (2012) 54–64

Contents lists available at SciVerse ScienceDirect

Engineering Geology journal homepage: www.elsevier.com/locate/enggeo

Flow path development in different CO2 storage reservoir scenarios A critical state approach J. Alonso a,⁎, V. Navarro a, B. Calvo b a b

Geoenvironmental Engineering Group, Civil Engineering School, University of Castilla-La Mancha, Spain Aragón Institute of Engineering Research (I3A), University of Zaragoza, Spain

a r t i c l e

i n f o

Article history: Received 29 July 2011 Received in revised form 20 December 2011 Accepted 5 January 2012 Available online 12 January 2012 Keywords: CO2 storage Geomechanics, coupled fluid flow and geomechanical modeling Elastic–plastic behaviour Critical state model

a b s t r a c t This study deals with the pressure-dependence of porous rock deformation, the permeability changes associated with this deformation, and the consequent effects of fluid flow during the CO2 injection phase in different configurations of deep saline aquifers. The study is approached from the standpoint of critical state soil mechanics, analyzing the effect exerted by different aspects that influence the geomechanical behaviour of porous rocks during CO2 storage. Through the study of the hydromechanical response of the environment, this paper shows how to assess the advisability of the use of CO2 storage, also identifying those cases and environments that would be less likely to lead to the creation or reopening of potential leakage paths. © 2012 Elsevier B.V. All rights reserved.

1. Introduction The geological storage of CO2 is deemed by the Intergovernmental Panel on Climate Change (IPCC, 2005) to be a viable strategy to mitigate climate change caused by anthropogenic factors. Pilot experiments are currently underway and some analogies have been made with other activities related to gas injection into geological formations, such as Enhanced Oil Recovery (Emberley et al., 2005; Jessen et al., 2005; Blunt et al., 2010). However, in view of the complex interaction between the fluid phases involved (Juanes et al., 2006; Pruess and Spycher, 2007), the large volumes of fluid considered in injection operations (Holloway, 2005), as well as the limitations in terms of the possible dissipation of excess water pressures (Zhou et al., 2008), it is necessary to conduct a specific study on the physical and chemical processes that take place in the context of deep geological CO2 storage. From a geomechanical point of view, CO2 injection will cause an impact due to changes in the stress field that result from changes to the pore pressure, buoyant pressure and volume of the rock (Orlic, 2009). The decreasing isotropic effective stresses that develop during CO2 injection entail an increase in volume, and therefore, in the porosity of the porous rocks that make up the storage reservoir. The increased porosity causes, in turn, a rise in the intrinsic permeability of the porous medium (Davies and Davies, 2001), favouring the

⁎ Corresponding author. E-mail address: [email protected] (J. Alonso). 0013-7952/$ – see front matter © 2012 Elsevier B.V. All rights reserved. doi:10.1016/j.enggeo.2012.01.001

migration of the fluid phases. Moreover, the reduction of the isotropic effective stresses during CO2 injection is usually accompanied by a variation in the deviatoric stresses, which may even generate plastic strains if the yield surface is reached. If this were to occur, then it might lead to the development of narrow tabular zones of localized strain, referred to as deformation bands (Aydin et al., 2006). The expected behaviour of the deformation band depends on the location of the stress point on the yield surface at the moment of localization. If it takes place on the dilatant side of the yield surface, it would be expected to produce a dilatant shear band (Fig. 1), where there would be a more pronounced increase in porosity and permeability than in the porous medium that is unaffected by the deformation band. Therefore, the dilatant shear bands could act as flow pathways, whereas the ones that were generated in a compactive regime, where a decrease in porosity and permeability would occur, would act as flow barriers. It is not common to use elastic–plastic models when simulating the geomechanical behaviour of CO2 reservoirs. However, they should be used to analyze the potential development of the flow barriers or pathways. Of all the different plastic models in existence, critical state models would perhaps be the most advantageous. These models allow for the reproduction of the most important aspects of the porous rocks (Cuss et al., 2003; Sheldon et al., 2006) that make up the deep saline aquifers considered as CO2 storage systems in this paper. The importance of identifying preferential flow paths is even greater when geomechanical coupling is taken into account. The development of conduit type structures in which CO2 leakage can be located requires the occurrence of dissolution reactions, which are

J. Alonso et al. / Engineering Geology 127 (2012) 54–64

Dilatancy

Dilatancy Fig. 1. Schematic dilatant shear band in microscopic scale. Solid arrows show the flow vectors (without scale).

only quantitatively important if a flow regime is associated (Gaus, 2010). Hence, the identification of flow pathways may help locate the areas in which leakage paths are most likely to develop. This would explain why it is of interest to characterize the different heterogeneities that may be introduced into the flow by the deformations caused by injection pressure. The purpose of this paper is to analyze the effect exerted by the flow boundary conditions, the tectonic regime, the depositional and diagenetic history of the system, as well as the geometry of the storage in the development of these heterogeneities, using a critical state model to characterize the mechanical behaviour.

2. Hydromechanical behaviour model In this section we define how the hydromechanical behaviour of the rocks is going to be modeled. Storage in deep saline formations entails the introduction of a new fluid phase rich in CO2 into a medium originally made up by a solid phase (formed by an aggregate of solid particles with varying amounts of mineral cement that acts as a binder) and an aqueous wetting phase that is rich in water (brine), containing varying amounts of dissolved salts. During the injection process, the CO2 rich, non-wetting phase, owing to its buoyancy, takes on the shape of a plume that partially displaces the brine in the aquifer. However, given the relatively low compressibility of the brine, overpressure can develop (especially during the injection phase but also thereafter) inducing large pressure gradients and reactivating natural pressure as well as groundwater flow regimes (Gaus, 2010), which would change the preexisting effective stresses in the geological media over a very extensive domain. The hydromechanical response of the medium is studied by means of a finite element model developed in Comsol Multiphysics software. The model makes use of three types of conservation equations, including mass balance equations for each fluid species and a momentum equation describing rock deformation forces. The following paragraphs describe briefly the most important aspects of the hydromechanical model.

55

2.1. Flow model Fluid phase flow is controlled by mass conservation equations. The mass balance of the species in their different phases is defined in Eqs. (1) and (2) in Table 1. In Eq. (1), D/Dt is the material derivative, vS is the velocity of the solid skeleton, ρϕ and qϕ are the averaged density and flow from phase ϕ (W for the liquid phase or “wetting”, and NW for the “non-wetting” phase rich in CO2), and m K is the mass of each species K (L, liquid, or CO2). Symbol “∇” defines the gradient operator. In Eq. (2), n is the porosity, SW is the degree of saturation in the wetting phase, and SNW corresponds to the non-wetting phase, where ρW and ρNW are their densities. The solubility is defined through the mass fractions of each substance, XϕK (XϕK ≡ m K/mϕ). Moreover, in Table 1, Eqs. (3) and (4) define the flow. In Eq. (4), the first term on the right describes the advective transport of the species in the phase, found by using Darcy's Law. Here, because isotropic permeability is assumed, Kϕ = Kδij, with Kϕ being the intrinsic permeability matrix, K the scalar value of permeability and δij the delta kronecker. μϕ, Pϕ and ρϕ are, respectively, the average dynamic visK cosity, the pressure and the average density of the phase. krϕ is the relative permeability of species K in phase ϕ. ∇z expresses the gradient of the vertical coordinate z. The second term in Eq. (4) describes the dispersive mass flow of the species with respect to the average phase flow, calculated by means of Fick's law. DV is an effective coefficient of molecular diffusion in a porous medium. It depends on the temperature, gas pressure, tortuosity of the medium and the saturation. ∇XϕK defines the gradient of the mass fraction of species K in phase ϕ. In keeping with Pruess and Spycher (2007), here the density and viscosity values of the CO2 rich phase are approximated by the corresponding properties of pure CO2, without water present. The density of the CO2 was calculated by means of the equation formulated by Redlich and Kwong (1949), according to the method proposed by Spycher et al. (2003), whereas the dynamic viscosity of the CO2 was obtained by means of the explicit expression reported in McPherson et al. (2008). Also in accordance with Pruess and Spycher (2007), the solubilities between CO2 and water in the wetting and non-wetting phases, respectively, were obtained on the basis of the correlations of Spycher et al. (2003). 2.2. Geomechanical model Geomechanical modeling of the geological storage of CO2 commonly uses linear elastic constitutive models (Rutqvist et al., 2002; Rutqvist et al., 2008). Some papers consider the generation or reactivation of fractures through a failure envelope defined by means of a Mohr–Coulomb model (Rutqvist et al., 2007; Chiaramonte et al., 2008; Orlic, 2009; Shi and Durucan, 2009; Vidal-Gilbert et al., 2009). However, soft porous rocks such as many sandstones (see Sheldon et al., 2006 and Coop and Willson, 2003, for example) or chalks (see Collin et al., 2002) show a mechanical behaviour that differs substantially from linear elastic behaviour, which could be characterized by the following aspects: i) non-linear bulk modulus, ii) the presence of a yield surface that defines a change in trend between the elastic regime and a regime of non-recoverable deformations and iii) mechanisms of plastic deformation associated with the degradation of the

Table 1 Balance equations used for the fluid phases. Description

Equation

Mass balance of species K

DS mK Dt K

Mass of species K Flow mass of species K Flow mass of species K in phase ϕ

#


þ mK ∇⋅vS þ ∇ ρK qK ¼ 0

(1)

K K m = nSWρWXW + nSNWρNWXNW

(2)

K K K K qW + ρNW qNW ρKqK = ρW  K  K K K K φ krφ ∇P φ þ ρφ g∇z −ρφ Dv ∇X Kφ ρφ qφ ¼ −ρφ μ

(3)

φ

(4)

56

J. Alonso et al. / Engineering Geology 127 (2012) 54–64

Table 2 Equations used for the definition of the elastic–plastic constitutive model.

Table 3 Suction bonding definition in the elastic–plastic model.

Description

Equation

#

Equation

#

Elastic behaviour

dεev ¼ p″κ⋅υ ⋅dp″

(5)

dq dε eS ¼ 3G

(6)

p″S(y) = (p″S(0))A(y) ⋅ exp B(y) AðyÞ ¼ λðλ−κ yÞ−κ

(13) (14)

(7)

BðyÞ ¼ NλððyyÞ−N Þ−κ

(15)

λ(y) = λ ⋅ ϕ (y) N(y) = N ⋅ ϕ (y) ϕ(y) = 1-a ⋅ (1-exp(b ⋅ y)) y = (1 − SW) ⋅ Y(s)

(16) (17) (18) (19)

s=P ATM Y ðsÞ ¼ 1 þ 10:7þ2:4⋅s=P ATM

(20)

″

Þ⋅p ⋅υ G ¼ 3ð1−2⋅ν 2ð1þν Þ⋅κ



f ≡g≡ p″ þpp ″ þp″ T S M



−η M⋅ð1−α Þ

1−α M⋅α ⋅η

− exp þ1 Yield loci and plastic ⋅ where potential ″ ″ p¼p þp T q η¼ p ∂f Flow rule dεP ¼ dΛ⋅ ∂σ ″  ‘
 ″ ″ Hardening rule dp″ S ¼ ∂p∂p″ ðS0Þ ⋅ λ−κ ⋅ dεpV þ ξ⋅dεpS þ ∂p∂yS ⋅dy




α ð1−α Þ2

S

(8)

(9) (10)

dp″M = − p″M ⋅ ρM ⋅ |dεvp|

(11)

p″T = ρTp″M

(12)

cement and the generation of water meniscuses between the solid particles. All of these aspects have been taken into account in the mechanical model used in this article by means of the critical state model with associated plasticity outlined in the following paragraphs. A more detailed description is given in Navarro et al. (2010). Elastic behaviour is defined by a hypoelastic model based on what is experimentally observed in triaxial compression tests. Elastic volume and deviation deformations are presented in Eqs. (5) and (6) of Table 2, where κ is the slope for isotropic unloading or reloading paths, G is the shear modulus (whose expression is found in Eq. (7) of Table 2), and ν, Poisson´s ratio. Like the plastic potential, the yield surface is defined by Eq. (8) of Table 2, where p″S, p″T and p″M are the plastic variables used in the model, and α and M would be the parameters controlling the shape of the yield surface. In turn, p″ is the effective mean stress, and q is the deviatoric stress. In the equations, the effective stresses are defined as σ″ = σ − Pδij, where P is the total fluid pressure (P = PWSW + PNWSNW). Of the three plastic variables considered in the model, p″S introduces the effect of the granular structure and the impact of the effects exerted by the water meniscuses on the solid skeleton. On the other hand, p″M and p″T are related to the presence of cement between the solid particles, as proposed by Nova et al. (2003). Through the variables p″M and p″T it is possible to extend the yield surface initially characterized by p″S as can be seen in Fig. 2. The cancelling out of variables p″T and p″M would lead to the complete destruction of the cement, which means that the yield surface would be determined by the effect of the suction and the granular skeleton. The plastic strains are calculated in consideration of the classic theory of plasticity, by means of the flow rule described in Eq. (9) of Table 2, where dΛ

would be the plastic multiplier. The hardening rules are presented in Eqs. (10) to (12) of Table 2, and they consider the evolution of the plastic variables. As can be observed in Eq. (12), p″T and p″M would be inter-related and would define a first type of “bonding” by intergranular cementation. The second type of “bonding”, in keeping with the idea proposed by Gallipoli et al. (2003), and adopted by Borja (2004), is associated with the presence of water in the form of meniscus lenses and is characterized by a bonding variable (Gallipoli et al., 2003), y. This variable causes an increase in p″S, the internal variable which characterizes macroscopic fabric behaviour. In accordance with Gallipoli et al. (2003), this effect may be characterized by Eqs. (13), (14) and (15) of Table 3. In these equations, parameters λ and κ are the slopes obtained by plotting the void ratio against the natural logarithm of the isotropic effective stress for the saturated condition of the normal compression line and the unloading–reloading lines, respectively. N is the intercept of the saturated normal compression line at p″ = 1 MPa. The compliance λ(y) and the intercept N(y) at “suction bonding” y are captured by means of Eqs. (16), (17) and (18) of Table 2 (Gallipoli et al., 2003), where a and b are two material parameters of the model. In Eq. (13) p″S (0) is the hardening parameter at a saturated state. It is assumed to include the dependence of the evolution of p″S with regard to the plastic strains (Eq. 10, Table 2). The interparticle bonding y is defined as the product of two factors: the degree of saturation of the non-wetting fluid, 1 − SW, and a function Y(s) which accounts for the increment with the increase in the stabilizing interparticle force exerted by a single meniscus (Gallipoli et al., 2003). This results in the multiplicative structure introduced in Eq. (19) of Table 3. In this table, for each suction (s = PW − PNW), function Y(s) is proportional to the squared average particle diameter (Cornelis et al., 2004). If there is an average diameter of 1 μm, then Eq. (20) in Table 3 can be used, in keeping with the approach proposed by Borja (2004).

q DILATION

COMPACTION

p p

T

p

S

p

M

Fig. 2. Yield surface considered in the geomechanical model.

J. Alonso et al. / Engineering Geology 127 (2012) 54–64

stress and of the deviation that takes place in the elastic–plastic stress–strain paths.

2.3. Hydromechanical coupling As David et al. (2001) points out, it is of fundamental importance in geophysics to understand the interplay between deformation, microstructural evolution and fluid flow in porous rocks. In the context of the geological storage of CO2, the correct definition of this relation influences the storage capacity itself and even the possible generation of leakage paths or flow barriers. In this article the intrinsic permeability K is correlated to porosity according to the function proposed by Rutqvist et al. (2002), and modified from Davies and Davies (2001): K ¼ K 0 exp½cðn=n0 −1Þ

57

ð21Þ

where n is the porosity, n0 and K0 are the porosity and the permeability at zero stress, and c is a fitting parameter that should be experimentally determined. The evolution of the porosity is related to the volumetric strains, which, in the mechanical constitutive model adopted, may develop as a consequence of variations in the mean

3. Study for different storage system

boundary

conditions

in

the

CO2

The effect of the flow contour conditions on the geomechanical behaviour of the storage reservoirs is studied by means of a modification of the three synthetic models presented in Zhou et al. (2008). The modification carried out in this article considers a system comprising two lithologies (Fig. 3), a formation of sandstones with high permeability and 100 m thickness, and a seal formation of shales of 50 m thickness. Plane strain and plane flow conditions is assumed. The injection takes place at point P located in the centre of the sandstone formation. An initial stress distribution is assumed in which the pressures of the wetting phase and the vertical stress are defined by the weight of a sediment column of 1350 m over the injection point. The initial stresses are related by σ″1 / σ″3 = 1.25 and σ″1 / σ″2 = 1.10, i.e. a compressive stress regime. The liquid degree of saturation is defined by means of the retention curve of Van Genuchten

a

b

c

Fig. 3. Schematic of model geometry with configurations (a) “open”, (b) “closed” and (c) “semi-closed”, modified from Zhou et al. (2008). Reference points P and Q are shown.

58

J. Alonso et al. / Engineering Geology 127 (2012) 54–64

Table 4 Hydraulic and mechanical parameters used in the analysis. Material

Aquifer

Cap

M α λ κ ξ ν a b N ρM ρT p″S (MPa) p″M (MPa) p″T (MPa) n0 P0 (kPa) λ SrW SrNW c K0 (m2)

1.00 0.5 0.16 0.01 0.1 0.3 0.2 0.1 1.65 12.0 0.095 p″(t = 0) 15 0.5 0.3 19.6 0.457 0.3 0.05 22.2 1 × 10− 13

1.00 0.5 0.16 0.01 0.1 0.3 0.2 0.1 1.65 12.0 0.095 p″(t = 0) 1 0.05 0.01 3100 0.457 0.3 0.05 22.2 1 × 10− 17

(1980), while the relative permeabilities of the wetting and nonwetting phase are calculated using the expressions proposed by Corey (1954). The hydraulic and mechanical parameters used in the

models are presented in Table 4. In this table, parameters P0 and λ define the retention curve, whereas the irreducible saturations of the wetting (SrW) and non-wetting (SrNW) phase are used to calculate the relative permeabilities. In all cases an injection of 0.05 kg/s per linear meter is assumed, which would be equal to 0.025 kg/s per linear meter, taking the symmetry of the problem into account. In this work, we have modelled three different boundary conditions. The first consideration is an “open” (Fig. 3a) configuration, broadly extended in a horizontal direction, which allows the fluid phases to flow in this direction, but not vertically. This configuration refers to the existence of a low permeability geological formation over the seal formation where the injected CO2 plume displaces horizontally to the water-rich phase. The velocity at which the two fluid phases are displaced is mainly determined by the relationship between the injection volume and the permeability of the medium. The second consideration is a “closed” system (Fig. 3b), where the storage takes place in a medium confined, both vertically and horizontally, by low permeability materials. The last case under consideration is a “semi-closed” medium (Fig. 3c), which permits the vertical circulation of the fluid phases, but at a velocity that depends on the permeability and thickness of the rock above the shales. This effect has been idealized by assuming that the pressure of the liquid corresponds to the weight of the water column. The imposition of this Dirichlett type condition implies that the medium located over the storage system has an infinite permeability. This simplification has

Fig. 4. Total fluid pressure P evolution at points P (grey) and Q (black) in configurations (a) “open”, (b) “closed”, (c) “semi-closed”, (d) “extensional”, (e) “overpressured” and (f) “thick”.

J. Alonso et al. / Engineering Geology 127 (2012) 54–64

been chosen so as to focus the study on the analysis of the effect of geomechanical behaviour. This study examines the evolution of the fluid pressure and stress– strain paths in the sandstones (point P) and shales (point Q). During the initial injection phase there is an increase in the average fluid pressure P in points P and Q in all the configurations used (Fig. 4). In configurations “open” and “closed” (Fig. 4a and 4b, respectively) the increase in P takes place continuously throughout the injection process. In contrast, the “semi-closed” system exhibits a more complex behaviour. As the injection process progresses, an ascending plume of CO2 is generated, which after a certain period of time, reaches the lower contour of the shales. At this time, part of the CO2 remains confined in the upper area of the sandstones, since initially the air entry pressure of the shales is not exceeded. As a consequence of the high saturation level in the CO2 rich phase of the sandstones in the vicinity of the injection area, the relative permeability of the liquid is reduced, for which reason the flow of the aqueous phase towards the shales is interrupted and the increase in P is reduced (Fig. 4c). Moreover, the non-wetting fluid pressure PNW imposed in the upper contour of the shales dissipates the possible increases in fluid pressure that may occur at point P. This phenomenon is illustrated in Fig. 5, where the area recording the highest saturation in the non-wetting phase occurring 2000 days from the start of the injection shows the greatest decrease of P in the shales, which increases gradually towards the zone to the right of the CO2 plume, where the saturation in the non-wetting phase is lower. In contrast, the total pressure P, does not exhibit such a strong change in trend at point P located in the storage formation (Fig. 4c). This is due to the fact that, although the flow is also conditioned by the dissipation of the excess fluid pressures as a result of the condition imposed in the

upper contour of the system, it is confined by the contact with the shales in the region where the vertical flow is interrupted. Fig. 6 shows the stress paths and the initial yield function for each of the configurations under consideration, while Fig. 7 exhibits the evolution of the value of the yield surface f (Eq. 8) at point Q versus the time elapsed from the start of the injection, as a criterion to define the beginning of the plastic straining. The development of the plastic strains in the shales begins earlier in model “closed” (Fig. 7). This result is predictable since the dissipation of the increase in fluid pressure is not allowed. Therefore, the effective stresses are reduced more quickly and the stress path reaches the yield surface here before it does in the rest of the configurations. In system “open” (Fig. 6a), on the other hand, the excess fluid pressure is dissipated better since the migration of the fluid phases is not as restricted. The specific evolution of the fluid pressure in case “semi-closed” (Fig. 6c) causes a maximum value of f to be reached in the shales after roughly 1000 days. Until that point in time, given the initial conditions included in the problem, the whole stress path occurs in an elastic regime. Once the total fluid pressure begins to decrease, there is an increase in the effective stresses, which brings the stress path closer to the yield surface on the compactive side. Fig. 4 shows that in the sandstones the increase of the total pressure is quicker than in the shales for all the cases under study. Moreover, in the model proposed, the degree of cementation, characterized by the plastic variable p″M, is much higher in sandstones than in shales. For this reason, despite presenting more unfavourable stress paths, the yield surface is not reached during the entire computation time under consideration (Fig. 7). A decrease of p ″M in the sandstones would imply that the development of plastic deformations could occur previously in the aquifer than in the shales, which may alter the flow regime in the storage reservoir.

SNW

5 km

0.025

59

Q

Cap

P

Aquifer

P (MPa)

t = 2000 days

kg s

1000 m

Fig. 5. Non-wetting phase degree of saturation, fluid pressure contours (MPa) and flow vectors of the wetting phase (brine) in model “semi-closed”, after 2000 days of injection.

60

J. Alonso et al. / Engineering Geology 127 (2012) 54–64

Fig. 6. Stress paths (continuous line) and initial yield surface (dotted line) at points P (grey) and Q (black) in configurations (a) “open”, (b) “closed”, (c) “semi-closed”, (d) “extensional”, (e) “overpressured” and (f) “thick”.

Fig. 8 depicts the evolution of intrinsic permeability, in terms of porosity (Eq. 21). It is possible to see greater increments in permeability in cases “closed” and “open” than in case “semi-closed”, due to the fact that the decrease in isotropic stress in both cases,

which is a consequence of the variation in fluid pressure, is more important. Moreover, in case “closed”, where the yield surface is reached earlier, the plastic strains generate higher volumetric strains and therefore greater variations of intrinsic permeability.

Fig. 7. Yield surface value evolution f at point Q (caprock) for the different analysis performed in the paper. Plastic strains develop when f = 0.

J. Alonso et al. / Engineering Geology 127 (2012) 54–64

1.0E-08 Permeability at pointQ (caprock)

1.0E-09 1.0E-10

K(m2)

1.0E-11 1.0E-12 1.0E-13

Open O Closed C

1.0E-14

Semi-closed SC

1.0E-15 1.0E-16 1.0E-17 0

2000

4000

6000

8000

t(days) Fig. 8. Intrinsic permeability evolution (m2) at Q, in the caprock, for models “open”, “closed” and “semi-closed”.

Although huge variations in permeability are observed in cases “closed” and “open”, the results are comparable to those obtained by considering the parameters of the power law reported by Dong et al. (2010) in shales.

61

The models analyzed show that it is in the systems confined by very low permeability covers where the highest increases of P, and therefore the most important reductions in effective stress, occur. The greater the decrease in effective stresses, the greater the likelihood that it will intersect with the yield surface and plastic strains will be generated. Hence, it may be concluded that the first two cases (“closed” and “open” in the models analyzed) favour the development of plastic strains in the shales, with the possibility that preferential flow paths may develop as a result of the generation of location phenomena in a regime exhibiting plastic behaviour. As an example Fig. 9 shows the preferential leak paths of model “closed”, associated with the development of deformation bands in the formation of shales. In contrast, with more permeable geologic formations over shales, the flow regime generates an increase of the value of the suction at the point of contact between the formations, as can be seen in Fig. 10a. This would be accompanied by an increase in the degree of saturation in the non-wetting phase (Fig. 10b) and may cause the air entry pressure to be exceeded. These conditions (system “semi-closed” in the model analyzed), could lead to a diffuse loss of CO2 towards the overlying formations as a consequence of the overcoming of the air entry pressure in the shales. 4. Tectonic regime The influence of the tectonic regime is illustrated by comparing the response of the closed system (Fig. 3a) to the response that this system would have given the initial existence of an extensional tectonic regime, with an initial relation of stresses σ″1/σ″3 = 0.75 and

Fig. 9. Volumetric plastic strain and associated deformation bands for the model “closed”.

62

J. Alonso et al. / Engineering Geology 127 (2012) 54–64

a

150

Shales

125

Height (m)

100

75

Sandstones 50

OOpen CClosed

25

Semi-closed SC

0 0

1000

2000

3000

4000

5000

6000

7000

8000

9000

s (kPa)

b

150

Shales

125

Height (m)

100

75

Sandstones 50

OOpen CClosed

25

Semi-closed SC

0 0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

SNW Fig. 10. (a) Suction and (b) degree of saturation for the non-wetting phase at the left boundary for models “open”, “closed” and “semi-closed” after 2000 days of injection.

σ″1/σ″2 = 0.85. Although the influence exerted by the stress state on geomechanical behaviour has been studied previously by other authors (see Rutqvist et al., 2007, for instance), this paper examines the question using a critical state model approach. The trajectory of stresses obtained in the extensional model is presented in Fig. 6d. While in the compressive regime (“closed”) the stress paths are characterized by a deviatoric component that increases with the decrease in isotropic stress (Fig. 6b), in model “extensional model” the deviatoric stress is reduced initially as the effective pressures drop. Moreover, the evolution of the injection pressure in the case of “extensional” model (Fig. 4d) is slower than in “closed”, and the yield surface is reached later (Fig. 7). As can be observed in Figs. 6b and 6d, in both cases, the stress paths reach the yield surface in the brittle area. However, although in the compressive regime the

intersection of the stress path with the yield surface occurs in a region that is very close to the critical state line, in the extensional tectonic regime the intersecting point is located in a zone towards the left within the yield surface. Therefore, since associated plasticity is adopted, in an extensional tectonic regime there would be a higher degree of volumetric plastic strain, so that the variations in porosity and permeability would be greater than in a compressive regime, where the volumetric deformations did not have such a high level of development. Similarly, as a result of the location of the stresses that could occur in the brittle regime, the deformation bands developed will be wider and will exhibit greater permeability as compared with the compressive regime. This would encourage the migration of the fluids inside, thus reducing the effectiveness of the storage system.

J. Alonso et al. / Engineering Geology 127 (2012) 54–64

5. Burial history/overpressured systems During deposition, shales that enclose sand lenses acquire higher fluid pressures during burial than pure shales (Ozkaya, 1986). Therefore, high excess fluid pressures may develop in shales at depth, where effective stresses decrease. To assess this aspect in terms of the potential geological storage of CO2, we studied, again, a “closed” system, analogous to the one presented in Fig. 3b, but considering here a fluid overpressure of 3 MPa (Model “overpressured”) as compared to the original model. The results are shown in Figs. 4, 6, and 7. In the “overpressured” model (Fig. 6e) the intersection of the stress path with the yield surface barely changes versus model “closed” (see Fig. 6b). Hence the plastic deformations begin at practically the same time as in the case of a hydrostatic distribution of the fluid pressure. Additionally, the evolution of the total fluid pressure at the intersecting point is very similar to what was recorded in the reference case (Fig. 4e). It may be concluded that, under the conditions considered here, the existence of an overpressure of 3 MPa has no major effect on CO2 storage capacity. 6. Aquifer geometry and thickness Synthetic exercises have been carried out in previous sections. These exercises were based in simple geometries, characterized by the horizontal deposition of the geological units. In reality this is not usually the case. Among the geometric aspects that influence the

63

stress–strain paths of storage reservoirs during the geological storage of CO2, a fact that stands out is the disposition of lithological units conditioned by the presence of folds or faults. As an example, Orlic (2009) studied the effect of the geological structure on the evolution of the stress paths during gas extraction operations by means of a simplified hydromechanical model in which a single fluid phase was considered. The study showed the concentration of stresses in regions associated with geometric singularities of the geological structure, where the development of plastic deformations or the generation of breakage mechanisms would be more likely. However, the study of the effect of geometry is complex and should be studied on a case by case basis. Thus, this article has been restricted to the study of the effect of another geometric factor, the aquifer thickness, which also influences the geomechanical behaviour of the geological storage of CO2. To this end, model “closed” has been modified, considering a sandstone thickness of 150 m (Model “thick”) instead of 100 m as considered in the original model (Fig. 3b). Fig. 6f shows the stress path and Fig. 4f the evolution of the fluid pressure. As can be seen, a greater thickness of sandstones allows for the higher mobility of the fluid phases, so that an increase in fluid pressure takes place more slowly in “thick” than in model “closed”, which is thinner. Therefore, the intersection with the yield surface, and, consequently the beginning of plastic straining, occurs later than in the reference model (Fig. 7). This is corroborated in Fig. 11, which shows the scope of the CO2 plume generated after 1620 days for models “closed” and “thick” (from 100 to 150 m thick, respectively).

a

b

Fig. 11. Non-wetting phase degree of saturation after 1620 days of injection in a) 150 m thickness (thick) aquifer b) 100 m aquifer thickness (closed).

64

J. Alonso et al. / Engineering Geology 127 (2012) 54–64

7. Conclusions The study of the stress–strain paths that occur during the geological storage of CO2 by means of critical state models has made it possible to reproduce the evolution of the porosity and permeability of porous rocks. Moreover, with the development of simplified models, it is easier to predict the compactive or dilatant trends in the medium. Therefore they may be used to assess the possible generation of preferential flow paths in storage reservoirs. On the basis of the analyses carried out, it may be concluded that the boundary flow conditions present in a storage reservoir of CO2 have an important effect on the possible development of plastic strains since they condition the dissipation of the increases in fluid pressure caused by the injection. Hence, systems that facilitate the migration of the fluid phases will not show such striking increases in pore pressure, and therefore they will reach the yield surface and the start of plastic strains later. Otherwise, these boundary conditions will lead to the displacement of the brine at great distances from the point of injection, which may raise other environmental issues. In the model that considers an extensional regime, the continual injection of CO2 may create more secure initial stress–strain paths than in a compressive regime, where the yield surface is reached earlier (Fig. 7). However, once the intersection with the yield surface has taken place, it occurs in the brittle area of the region in which the volumetric plastic strain is greater, leading to the generation or reopening of wider preferential flow paths. In this analysis, the presence of an overpressure in the aquifer does not have a significant long-term effect on the behaviour of the storage system in the model presented. Nevertheless, since the overpressure conditions the origin of the stress paths and, along with the plastic parameters, the relative position and distance with regard to the yield surface of the porous rocks, it would indeed play an important role in other situations. Lastly, in the aquifers having a greater thickness, the increase in fluid pressure will be slower. This means that the beginning of the plastic strains will occur later, thus allowing for a greater storage capacity. Acknowledgements This research was financed in part by the Research Grants of the Spanish Ministry of Science and Innovation PSS-120000-2007-27 and PSS-120000-2008-31. References Aydin, A., Borja, R.I., Eichhubl, P., 2006. Geological and mathematical framework for failure modes in granular rock. Journal of Structural Geology 28, 83–98. Blunt, M., Fayers, F.J., Orr Jr., F.M., 2010. Carbon dioxide in enhanced oil recovery. Energy Conversion and Management 34, 1197–1204. Borja, R.I., 2004. Cam–clay plasticity. Part V: a mathematical framework for threephase deformation and strain localization analyses of partially saturated porous media. Computer Methods in Applied Mechanics and Engineering 193, 5301–5338. Chiaramonte, L., Zoback, M.D., Friedmann, J., Stamp, V., 2008. Seal integrity and feasibility of CO2 sequestration in the Teapot Dome EOR pilot: geomechanical site characterization. Environmental Geology 54, 1667–1675. Collin, F., Cui, Y.J., Schroeder, C., Charlier, R., 2002. Mechanical behaviour of Lixhe chalk partly saturated by oil and water: experiment and modelling. International Journal for Numerical and Analytical Methods in Geomechanics 26 (9), 897–924. Coop, M.R., Willson, S.M., 2003. Behaviour of hydrocarbon reservoir sands and sandstones. Journal of Geotechnical and Geoenvironmental Engineering, ASCE 129 (11), 1010–1019. Corey, A.T., 1954. The interrelation between gas and oil relative permeabilities. Producers Monthly 19 (1), 38–41. Cornelis, W.M., Gabriels, D., Hartmann, R., 2004. A conceptual model to predict the deflation threshold shear velocity as affected by near-surface soil water: I. Theory. Soil Science Society of America Journal 68 (4), 1154–1161.

Cuss, R.J., Rutter, E.H., Holloway, R.F., 2003. The application of critical state soil mechanics to the mechanical behaviour of porous sandstones. International Journal of Rock Mechanics and Mining Sciences 40, 847–862. David, C., Menendez, B., Zhu, W., Wong, T.F., 2001. Mechanical compaction, microstructures and permeability evolution in sandstones. Physics and Chemistry of the Earth, Part A: Solid Earth and Geodesy 26, 45–51. Davies, J.P., Davies, D.K., 2001. Stress-dependent permeability: characterization and modeling. SPE Journal 6, 224–235. Dong, J.J., Hsu, J.Y., Wu, W.J., Shimamoto, T., Hung, J.H., Yeh, E.C., Wu, Y.H., Sone, H., 2010. Stress-dependence of the permeability and porosity of sandstone and shale from TCDP Hole-A. International Journal of Rock Mechanics and Mining Sciences 47, 1141–1157. Emberley, S., Hutcheon, I., Shevalier, M., Durocher, K., Mayer, B., Gunter, W.D., Perkins, E.H., 2005. Monitoring of fluid–rock interaction and CO2 storage through produced fluid sampling at the Weyburn CO2-injection enhanced oil recovery site, Saskatchewan, Canada. Applied Geochemistry 20, 1131–1157. Gallipoli, D., Gens, A., Sharma, R., Vaunat, J., 2003. An elasto-plastic model for unsaturated soil incorporating the effects of suction and degree of saturation on mechanical behaviour. Geotechnique 53, 123–135. Gaus, I., 2010. Role and impact of CO2–rock interactions during CO2 storage in sedimentary rocks. International Journal of Greenhouse Gas Control 4, 73–89. Holloway, S., 2005. Underground sequestration of carbon dioxide—a viable greenhouse gas mitigation option. Energy 30, 2318–2333. IPCC, 2005. IPCC special report on carbon dioxide capture and storage. In: Metz, B., Davidson, O., de Coninck, H.C., Loos, M., Meyer, L.A. (Eds.), Working Group III of the Intergovernmental Panel on Climate Change. Jessen, K., Kovscek, A.R., Orr, F.M., 2005. Increasing CO2 storage in oil recovery. Energy Conversion and Management 46, 293–311. Juanes, R., Spiteri, E.J., Orr Jr., F.M., Blunt, M.J., 2006. Impact of relative permeability hysteresis on geological CO2 storage. Water Resources Research 42. McPherson, B., Han, W.S., Cole, B.S., 2008. Two equations of state assembled for basic analysis of multiphase CO2 flow and in deep sedimentary basin conditions. Computers & Geosciences 34, 427–444. Navarro, V., Alonso, J., Calvo, B., Sanchez, J., 2010. A constitutive model for porous rock including effects of bond strength degradation and partial saturation. International Journal of Rock Mechanics and Mining Sciences 47, 1330–1338. Nova, R., Castellanza, R., Tamagnini, C., 2003. A constitutive model for bonded geomaterials subject to mechanical and/or chemical degradation. International Journal for Numerical and Analytical Methods in Geomechanics 27, 705–732. Orlic, B., 2009. Some geomechanical aspects of geological CO2 sequestration. KSCE Journal of Civil Engineering 13, 225–232. Ozkaya, I., 1986. Analysis of factors influencing excess heads in shales during burial. Marine and Petroleum Geology 3, 74–78. Pruess, K., Spycher, N., 2007. ECO2N—a fluid property module for the TOUGH2 code for studies of CO2 storage in saline aquifers. Energy Conversion and Management 48, 1761–1767. Redlich, O., Kwong, J.N.S., 1949. On the thermodynamics of solutions. 5. An equation of state-fugacities of gaseous solutions. Chemical Reviews 44, 233–244. Rutqvist, J., Wu, Y.S., Tsang, C.F., Bodvarsson, G., 2002. A modeling approach for analysis of coupled multiphase fluid flow, heat transfer, and deformation in fractured porous rock. International Journal of Rock Mechanics and Mining Sciences 39, 429–442. Rutqvist, J., Birkholzer, J., Cappa, F., Tsang, C.F., 2007. Estimating maximum sustainable injection pressure during geological sequestration of CO2 using coupled fluid flow and geomechanical fault-slip analysis. Energy Conversion and Management 48, 1798–1807. Rutqvist, J., Birkholzer, J.T., Tsang, C.F., 2008. Coupled reservoir-geomechanical analysis of the potential for tensile and shear failure associated with CO2 injection in multilayered reservoir-caprock systems. International Journal of Rock Mechanics and Mining Sciences 45, 132–143. Sheldon, H.A., Barnicoat, A.C., Ord, A., 2006. Numerical modelling of faulting and fluid flow in porous rocks: an approach based on critical state soil mechanics. Journal of Structural Geology 28, 1468–1482. Shi, J.Q., Durucan, S., 2009. A coupled reservoir-geomechanical simulation study of CO2 storage in a nearly depleted natural gas reservoir. Greenhouse Gas Control Technologies 9 (1), 3039–3046. Spycher, N., Pruess, K., Ennis-King, J., 2003. CO2–H2O mixtures in the geological sequestration of CO2. I. Assessment and calculation of mutual solubilities from 12 to 100 degrees C and up to 600 bar. Geochimica Et Cosmochimica Acta 67, 3015–3031. Van Genuchten, M.T., 1980. Closed-form equation for predicting the hydraulic conductivity of unsaturated soils. Soil Science Society of America Journal 44, 892–898. Vidal-Gilbert, S., Nauroy, J.F., Brosse, E., 2009. 3D geomechanical modelling for CO2 geologic storage in the Dogger carbonates of the Paris Basin. International Journal of Greenhouse Gas Control 3, 288–299. Zhou, Q.L., Birkholzer, J.T., Tsang, C.F., Rutqvist, J., 2008. A method for quick assessment of CO2 storage capacity in closed and semi-closed saline formations. International Journal of Greenhouse Gas Control 2, 626–639.