Impact of material properties on caprock stability in CO2 geological storage

Geomechanics for Energy and the Environment 11 (2017) 28–41

Contents lists available at ScienceDirect

Geomechanics for Energy and the Environment journal homepage: www.elsevier.com/locate/gete

Impact of material properties on caprock stability in CO2 geological storage Chao Li *, Lyesse Laloui Laboratory of soil mechanics - Chair ‘‘Gaz naturel’’ Petrosvibri, Swiss Federal Institute of Technology, EPFL, Lausanne, Switzerland

highlights • Numerical simulations are performed assess caprock stability under cooling and pressurization during the injection of CO2 . • Material property ratio R is defined to illustrate the failure potential of caprock. • Coupled effects of material properties on caprock stability are examined.

article

info

Article history: Received 29 November 2016 Received in revised form 27 April 2017 Accepted 20 June 2017 Available online 14 July 2017

a b s t r a c t In this study, thermo-hydro-mechanical coupled simulations are performed to investigate caprock stability under cooling and pressurization during the injection of carbon dioxide. It is shown that failure potential of caprock increases with a higher value of the material parameter R which is defined as a ratio of material properties (the thermal expansion coefficient and the elastic properties) of caprock and aquifer. The study also underlines that the consideration of the Biot coefficient of the aquifer is essential for properly estimating aquifer examination. The analysis shows that the cooling is mainly ) ( of caprock transferred by conduction, and the temperature drops with time t at a rate of erfc

√1 2 Dt

, where the rate

increases if caprock has a higher thermal diffusivity D.

© 2017 Elsevier Ltd. All rights reserved.

1. Introduction The storage of CO2 in geologic reservoirs, particularly those in deep aquifers, has become a mitigation method used to reduce the impact of CO2 and the greenhouse effect.1,2 There are currently several on-going large-scale projects for CO2 storage in deep aquifers, such as the In Salah CO2 storage site in Algeria.3 Aquifers are usually overlaid by very low permeable caprocks, which are of primary importance for preventing CO2 leakage. The integrity of the caprock is most likely to be deteriorated during the initial stage of the injection.4 The hydromechanical (HM) coupled phenomena involved in CO2 injection problems have received particular attention, as reviewed by Rutqvist.5 A large volume of CO2 injection results in an accumulation of overpressure. The overpressure rapidly perturbs the stress field upon injection in the aquifer and extends to the caprock. In addition, the injection temperature may be lower than that of the aquifer, as in the CO2 storage project at In Salah, Algeria,6 adding a thermal factor to the HM coupling. The cooling will contract rocks and reduces stress if constrained, whereas the author. * Corresponding E-mail address: [email protected] (C. Li). http://dx.doi.org/10.1016/j.gete.2017.06.003 2352-3808/© 2017 Elsevier Ltd. All rights reserved.

overpressure will result in an expansion and a stress reduction. Both effects may be counterbalanced or superposed depending on the specific configuration of the problem, which requires a fully coupled thermal-hydro-mechanical (THM) investigation. This has motivated numerous computational studies, including this study, which has two main objectives. The first objective is to assess whether the failure potential is attained when the caprock is subjected to a lower temperature, and the second objective is to highlight the importance of key factors that are involved in the THM processes. Regarding the former, low temperature injection may cause a reduction in compressive stress or even induce tensile stress in the caprock,7,8 whereas an improvement of the caprock stability has been found by Vilarrasa et al.9,10 Vilarrasa et al.9 stated that these studies are contradictory; however, they are actually complementary, as will be demonstrated in our study. Despite similar theoretical frameworks, the problem set-ups were different between these studies, and factors such as material properties also differed, which significantly influenced the results. This leads to the second objective, which is to determine the effect of material properties on caprock stability. Current analyses on material sensitivity mainly focus on the effects of each intrinsic rock property on reservoir performance, such as the monotonic influence on overpressure increase by an aquifer’s permeability and the lack of a

C. Li, L. Laloui / Geomechanics for Energy and the Environment 11 (2017) 28–41

significant influence on caprock stability when the Young’s modulus is varied, as found by Bao et al.11 and Vilarrasa et al.12 However, two key THM coupled properties of rocks are not investigated thoroughly in these studies: (i) the HM coupled property, namely, the Biot coefficient, and (ii) a thermal–mechanical property, namely, the thermal expansion coefficient. In addition, regarding the specific boundary configuration of CO2 injection problems (caprock overlying aquifer), the combined material properties of caprock and aquifer are believed to affect the response regimes and may also provide an explanation for the aforementioned contradictory conclusions. Furthermore, the thermal diffusion properties are also crucial for cooling transfer inside the caprock, the effect of which on caprock stability is not yet understood. The main goal of this paper is to address these gaps such that a more rigorous conclusion can be achieved when assessing caprock stability under low temperature CO2 injection during the initial stage. We first introduce the essential physical processes occurring during the CO2 injection stage. The THM fully coupled theoretical formulation is then presented in detail with the incorporation of CO2 nonlinear state properties, especially with an emphasis on coupled phenomena and thermo-dynamics and the necessary constitutive equations to describe them. Next, results are presented with a reference model, with a particular reference to the THM response in the caprock. The coupled effect of the caprock thermal properties, the Biot coefficients of the caprock and aquifer, and the caprock thermal expansion coefficient are investigated, which has to our knowledge not yet been addressed. These properties are not studied solely from both geomechanical and geometrical perspectives. The geomechanical approach is to couple the investigation with other parameters such as Young’s modulus to achieve more systematic insights. The objective from the geometrical point of view is to present the caprock stability results under the influence of the behaviour of the aquifer. Although numerical efforts place a burden on the feasibility of computing the material sensitivity of the properties, study cases were carefully chosen for their realism and to allow the extraction of novel quantitative information on the coupled effect of material properties on caprock stability. 2. Coupled processes in CO2 geological storage 2.1. Physical and geometrical description Fig. 1 represents the geometry of the model under investigation. The model consists of three geological strata involved in the CO2 storage problem, all of which are porous media: the sealing caprock, storage aquifer and underburden rock. Two phase fluid flow processes in such porous media can be described in the context of the theory of mixtures.13 This theory postulates the mass balance for each phase and each individual species, providing the explicit mass quantity during the phase change, which cancels out in the balance equations of the species under the compositional approach.14,15 To understand the coupling of the overpressure, temperature effect and thermal diffusion to the deformable porous media mechanics, it is necessary to define the appropriate conditions that are encountered in a CO2 storage project. One million tonne of CO2 is injected over one year, which is modelled as a prescribed CO2 mass flow through a vertical well along the thickness of the aquifer in a normal faulting system. The initial temperature of the reservoir is set to 330K. An injection temperature of 30 ◦ C is imposed, which is relevant within the range for onshore industrial projects such as In Salah CO2 storage project in Algeria6 and can also be confirmed by the findings of Lu and Connell.16 The aquifer is considered to be water saturated prior to the injection and acts as a host medium for CO2 . Low permeability and high air-entry value of the caprock prevent threats from potential leakage of

29

CO2 . A constantly distributed stress of 13.5 MPa is applied along the top of the caprock, which is equivalent to the dead load of the overburden. The initial hydrostatic fluid conditions are applied with the model being fully saturated with water. Therefore, the top of the aquifer exhibits a pressure of 8 MPa. The displacement on the right-hand side and bottom of the model is constrained in the perpendicular direction. The two main coupled processes are thus CO2 inflow into the storage aquifer from the injection well and the rapid generation of pressurization, inducing stress variation as the main HM coupling. The cooling occurs simultaneously adjacent to the injection well to cope with pressurization-induced expansion and provokes shrinkage of surrounding materials. The convective thermal flow, which is driven by the CO2 inflow, dominates the heat transfer in the aquifer. An upwards movement of cold flow occurs from the bottom of the caprock, leading to the THM coupling process, which is presented in the following section. 3. Thermal-hydro-mechanical formulation 3.1. Mass balance As stated previously, the compositional approach was employed for this study, as implemented for water and perfect gas by Collin15 using the finite element code Lagamine. This approach has the advantage of writing the mass balance equation for two phase fluids in a straightforward manner. We improved the current code with supercritical fluid properties, including terms for storage of both CO2 and water in liquid and gaseous forms, the advective flows of both fluids, and the non-advective flow of the dissolved CO2 in the water4 . The water phase pressure pw , CO2 phase pressure pc , temperature T and solid displacement field u were chosen as primary state variables to describe the state of the material. The mass balance equations and fluid flows were expressed in the moving current configuration using a Lagrangian updated formulation,17 through which the solid mass is automatically conserved for a reference elementary volume V : d ((1 − n) ρs V )

= 0, (1) dt where t is time, and n is the porosity, ( ∂ u ) and a volumetric deformation ≈ div . ρs is the solid density, for which rate is defined as V1 dV dt ∂t the equation of state is adapted from Lewis and Schrefler18 1 dρs

ρs dt

=

[

1

1 d

(Sw pw + (1 − Sw ) pc ) ( )] dT ∂u − 3 (b − n) αs − (1 − b) div , dt ∂t 1−n

(b − n)

Ks dt

(2)

where b is the Biot coefficient, Ks is the bulk modulus of the solid matrix, and αs is the linear thermal expansion coefficient for the solid. The water phase saturation Sw is a function of the capillary pressure (suction) s, which is defined as the difference between the CO2 phase and water phase pressures s = Pc − Pw . Introducing Eq. (2) into Eq. (1), the solid mass balance equation is then expressed in terms of the variation in the porosity by the primary state variables: dn dt

=

1 − n d ρs

ρs

dt

= (b − n)

[

+ (1 − n) div

(

∂u ∂t

)

1 d

(Sw pw + (1 − Sw ) pc ) ( )] dT ∂u − 3αs + div . dt ∂t Ks dt

(3)

This defines the coupling between the geomechanics and thermo-hydraulics. The change in porosity is explicitly shown as

30

C. Li, L. Laloui / Geomechanics for Energy and the Environment 11 (2017) 28–41

Fig. 1. Geometry of the considered storage system, location of result points.

being proportional to three terms: (a) a volumetrically weighted fluid pressure Pf = Sw Pw + (1 − Sw ) Pw , (b) temperature and (c) matrix volumetric deformation. The mass balance equation for the water specie and CO2 specie are written, respectively as

∂ (nSw ρw ) + div (ρw qw ) = 0,    ∂t

and

(4)

liquid water

supercritical/liquid CO2

(5)

dissolved CO2 in the water

Where the dissolution of CO2 in the water is accounted for but no water dissolution in the CO2 is involved, there is one more component in the mass balance of CO2 specie, i.e., Eq. (5), than in that of water, i.e., Eq. (4). It can also be noted that the mass exchange associated with the dissolution is not visible in Eq. (5) as a consequence of employing the compositional approach. Among these terms, the liquid water and supercritical/liquid CO2 flow qw and qc are governed by the generalized Darcy’s law for porous media to determine the advective flow: qw,c = −

( ) ] kkr w,rc [ grad pw,c + ρw,c g , µw,c

ρw = ρw0 [1 + κT (pw − pwr ) − βw (T − Tr )] ; 1 ∂ρw 1 ∂ρw = κT ; − = 3βw . (7) ρw 0 ∂ p w ρw 0 ∂ T ρw0 is the reference density for a given salinity at a reference pressure pwr and reference temperature Tr . κT is the isothermal water compressibility, and βw is the linear thermal expansion coefficient of water. The CO2 involved in storage cannot be considered to be a perfect gas for its high pressure high temperature condition. A compressibility factor Z is often used to describe the deviation from a perfect gas density to a real gas:

∂ (n (1 − Sw ) ρc ) + div (ρc qc ) ∂t    ∂ (nSw ρdc ) + + div (ρdc qw ) + div (idc ) = 0.  ∂t  

density ρw is

(6)

where k is the intrinsic permeability, and kr w and krc are the water and CO2 relative permeabilities, which are geomaterial dependent parameters. µw is the dynamic viscosity of water, which is calculated according to Thomas and King,19 and µc is the dynamic viscosity of CO2 following the relationship of Fenghour et al.,20 which is valid for the range of temperatures and pressures considered here. The liquid water is considered as a compressible and dilatant fluid, for which definition of the linearized relationship for the

1 ∂ρc ; = Z R T ρc ∂ pc ( ) 1 ∂ρc 1 1 dZ =− + . ρc ∂ T T Z dT

ρc =

1 M pc

(

1 pc

−

1 dZ Z dpc

) ; (8)

The factor Z is a pressure- and temperature-dependent parameter which is detailed in the Appendix. The advantage of expressing the density change as a function of Z is that one can express the compressible behaviour of gaseous, liquid and supercritical CO2 with appropriate equations of state (EOSs). In this study, the factor Z was calculated using Peng and Robinson’s EOS21 because of its accuracy and computational efficiency.22 The Henry law that describes the amount of a given perfect gas that dissolves in a volume of water is extended for a real gas with the gas fugacity coefficient 8.23 The dissolution amount is represented by the dissolved mass in a unit volume of water as:

ρdc =

8Z ρw RT co

Keq,2g −l Mw

ρc

(9)

co

where Keq,2g −l is the temperature dependent Henry’s constant determined according to Ref. 24 This relationship is used in the CO2 diffusion law, which is based on Fick’s law in a tortuous medium: idc = −nSw τ Dc ρw grad

(

ρdc ρw

)

,

(10)

C. Li, L. Laloui / Geomechanics for Energy and the Environment 11 (2017) 28–41

where Dc is the diffusion coefficient of the dissolved CO2 in the water phase, and τ is the tortuosity of the porous media.

31

3.3. Momentum balance The balance of linear momentum is written as

3.2. Energy balance

div (σ) + ρ g = 0,

Assuming a thermal equilibrium, a unique temperature is defined for the mixture. The energy balance equation of the mixture has the following form:

where σ is the total Cauchy stress tensor, with tensile stress taken as negative, and ρ is the bulk density, which is defined as

∂H ∂t 

+ div (Γ) + div (fT ) = 0,       conduction

Heat storage



(11)

convection





Heat transfer

where H is the enthalpy of the whole medium, and heat transfer is comprised of Γ, which is the heat conduction, and fT , which is the convection. The mixture enthalpy can then be defined as the sum of the heat of each constituent, neglecting the contribution of the dissolved CO2 in the water: H = (1 − n) ρs cp,s (T − T0 ) + nSr ρw cp,w (T − T0 )

+ n (1 − Sr ) ρc hc ,

(12)

where cp,α is the heat capacity of the component α . The term hc corresponds to the specific enthalpy of CO2 , which is temperature and pressure dependent. It can be determined through the fundamental thermodynamics relationship by employing the Peng and Robinson EOS: hc =

1 MC

( Z+

⎛ ln ⎝

Tc T

)

√

φ0 −

2a

[

4bM

(

c √ ) ⎞

(

√ ) ⎠,

Z + 1+ Z + 1−

κT 1+ √ α TTc

]

2 B

(13)

2 B

where a, b, B, α and κ are the Peng and Robinson EOS parameters, Tc is the critical temperature of CO2 (Tc = 304.25 K), and φ 0 is the Helmholtz energy for CO2 under the ideal gas condition as given by Span and Wagner.25 This formulation explicitly accounts the Joule– Thomson effect, of which the development is detailed in.4 The heat transfer is governed by the heat conduction and convection:

+ cp,w (T − T0 ) ρw qw + hc ρc qc .

ρ = (1 − n) ρs + nSw ρw + n (1 − Sw ) ρc .

(14)

The first term corresponds to the heat conduction, where the mixture conductivity is considered as a function of the volume ratio of solid, liquid water and CO2 phases. It can be observed that thermo-hydraulics and thermos-mechanical coupling are present through the porosity and saturation. As stated previously and presented later in the article, the most important thermal feature in the caprock is cast in the form of heat conduction. If considering a particular configuration such that the vertical one dimensional flux occurs from the bottom of the caprock, the energy balance equation Eq. (11) can be simplified under a water saturated condition:

∂T ∂ 2T λeq nλw + (1 − n) λS = D 2 and D = ( )eq = . (15) ∂t ∂z (1 − n) ρs cp,s + nρw cp,w ρ cp [ ] D m2 /s is the thermal diffusivity, which is a measure of the thermal inertia and the controlling parameter for heat conduction eq process, as defined as the ratio of ( the )eqmatrix conductivity λ and the volumetric heat capacity ρ cp . The cooling is faster if the material has a higher thermal diffusivity, and the mechanical stability will therefore be altered more quickly.

(17)

The behaviour of the solid matrix is assumed to be governed by the generalized effective stress tensor σ ′ through a combination of total stress and fluid pressures:

σ ′ = σ − b (Sw pw + (1 − Sw ) pc ) I,

(18)

where I is the identity matrix, and b is the Biot coefficient that is defined as b = 1 − KKm , where Km is the bulk modulus of the s

porous matrix, and Ks the bulk modulus of the solid matrix. An average fluid pressure is defined as follows and is weighted by the saturation of each phase: pf = Sw pw + (1 − Sw ) pc .

(19)

The effective stress in Eq. (15) thus becomes

σ ′ = σ − b p f I.

(20)

An important ingredient that contributes to the thermomechanical coupling is accounted for the definition of strain due to the phenomenon of thermal expansion. The following description of thermo-elastic strains is used: dε = E−1 dσ ′ − αs IdT ,

(21)

where dε is the total strain tensor increment, E is the linear elastic tensor, and αs is the thermal expansion coefficient. Using Young’s modulus E and Poisson’s ratio ν to characterize the elastic geomaterial, we can write this expression in a more explicit form using indexed notation: dεij = or dσij′ =

Γ + fT = − (nSw λw + n (1 − Sw ) λC + (1 − n) λS ) gradT

(16)

1+ν E E 1+ν

dσij −

{

ν E

dεij +

′ dσkk δij − αs dT δij

(22)

} E αs dT ν dεkk δij + δij . 1 − 2ν 1 − 2ν

From Eqs. (20)–(22), we find that during a temperature decrease at constant volume, the relationship can be described by dσ ′ = dσ − bdpf I and

dσ ′ = Eαs dT

or dσij′ =

E αs dT 1 − 2ν

δij . (23)

When evolving under a decrease of temperature alone and in the absence of fluid pressure variation, the stress will reduce as the product of the elastic properties, temperature variation and thermal expansion coefficient is increased. This highlights the mechanical influence of the choice of the thermal expansion coefficient. A large temperature increase and high magnitude of the elastic modulus will certainly reduce more stress, as illustrated by Preisig and Prévost.7 On the other hand, when a constraint on the temperature is imposed, as in the case of isothermal injection or in the area that the temperature front has not yet reached, the effect of the Biot coefficient on the stress reduction is evident for a given fluid pressure increase. The resolution follows a monolithic approach to formulate the stiffness matrix of the finite element method,17 of which the diagonal of the matrix is the uncoupled terms and the remainder is the coupled terms. The fully coupled stiffness matrix, with regard to the uncoupled stiffness matrix, is

32

C. Li, L. Laloui / Geomechanics for Energy and the Environment 11 (2017) 28–41

completely modified at each iteration by the Newton–Raphson algorithm15 ; all derivatives are computed, and all the coupled effects are therefore updated to converge the solution to an equilibrium state. 4. Model characteristics It is necessary to describe the materials involved in this study. Sealing caprocks usually have very low permeabilities and porosities, whereas aquifers are much more porous, and we therefore consider clayey shale as the caprock and sandstone as the storage rock, as in most on-going CO2 storage projects.3 The fluid thermodynamic properties were taken from Jamieson et al.26 and Wilson,27 and the geomaterial parameters were taken as typical averaged values from Rutqvist and Tsang28 and Vilarrasa et al.10 for the reference model. The mechanical constitutive law is chosen as thermal-elasticity as the aforementioned studies. The modelling of the problem focuses on the thermal-hydro-elastic analysis of the interaction between the storage aquifer, the caprock and the low temperature CO2 . The choice of elasticity is pragmatic and relevant for the range of stresses where the caprock and the aquifer are not expected to behave plastically and without development of fractures.29–31 The rocks are initially saturated and later come into contact with the invading CO2 , and the water and CO2 pressures begin to generate with the development of suction. Therefore, the aquifer will be subjected to desaturation through changes in relative permeability, thereby influencing pressure generation. A van Genuchten function is used to describe the retention behaviour of rocks, and a power law is used for the relative permeability: Sw = 1 + (s/Pr )1/(1−m)

(

)−m

, and

CKW kr w = Sw krc = (1 − Sw )CKC ,

(24)

(25)

where m and Pr are a material parameter and a reference pressure, and CKW and CKC are material parameters. The parameters in the reference simulation are summarized in Table 1. The variations in the coupled parameters and thermal properties are specified in the corresponding analysis. The model consisted of 8723 quadrilateral elements and was run as an axisymmetric model. The initial stress equilibrium was obtained between the application of the body force and the lithostatic stress along the top of the caprock. The horizontal stress was the lateral earth pressure ratio, which was taken as 0.6 as for a normal stress faulting system. 5. Non-isothermal CO2 injection Two simulation scenarios are considered here: (i) isothermal injection scenario where the temperature are kept constant at all points during the simulation and the injection temperature is assumed as the same as reservoir temperature; (ii) non-isothermal injection scenario where the injection temperature is set as 30 ◦ C and the degree of freedom on temperature is free, allowing the variation of temperature. The two simulations described in Fig. 2 correspond to the stress and displacement evolutions of both the isothermal injection and low temperature injection. For the isothermal injection (Fig. 2(a)–(d)), injection-induced overpressure is the only driving mechanism for caprock stability reduction. The fluid pressure increased rapidly at the very beginning of the injection, whereupon both the horizontal and vertical effective stresses dropped quickly in the first dozens of days (Fig. 2(a) and (b)). An expansion in the aquifer was observed and pushed the caprock upwards, and a slight compression regime was thus triggered in the caprock because the vertical displacement was

restricted due to the bending moment,32,33 which was observed as a transfer delay of displacement as shown between −950 m to −1000 m in Fig. 2(d). The occurrence of this compression contributed a small addition of stress, which unfortunately remained at a level negligible to improve caprock stability. This also indicates that the stress state in the caprock depended on the expansion level of the aquifer. The results of the low temperature supercritical CO2 injection are shown in Fig. 2(e)–(h) to illustrate the thermal contribution to the reservoir mechanical performance. It can be noted that the cooling had a major influence on the horizontal components of the stress and displacement (Fig. 2(e) and (g)) and a moderate influence on the vertical components (Fig. 2(f) and (h)). The cooling significantly reduced the vertical displacement at 300 days (Fig. 2(d) and (h)) and altered the horizontal displacement regime from expansion to contraction (Fig. 2(c) and (g)). Although a free contraction upon cooling was not fully allowed due to the boundary condition on the left-hand side, whereupon such a deformation constraint turns into a significant horizontal stress reduction upon cooling (recall Eq. (25)). There was nearly twice as much stress reduction in the aquifer upon cooling as that without cooling in the horizontal direction and even more reduction at the aquifer– caprock interface over both the short and long term, as can be clearly observed when comparing Fig. 2(a) and (e). On the other hand, the vertical displacement was less constrained, although the level of stress reduction was as significant as the horizontal one in the aquifer at 10 days, as illustrated by the results in Fig. 2(b) and (f). This can be explained by the fact that the expansion was sustained by the overpressure, which was being generated inside the aquifer over the short term; the deformation upon cooling was thus mostly restricted by such an expansion and was converted to a stress reduction in return. Over the 300 days of injection, part of the effective stress reduction was recovered in the aquifer due to an overpressure drop. The overall decrease in the vertical effective stress was thereby reduced in the aquifer. In the lower portion of the caprock, the vertical effective stress continued to decrease, which was a direct consequence of the cooling-induced stress reduction under a constrained deformation condition. Fig. 2(a), (b), (e) and (f) clearly show that, in this particular study with the chosen parameters in Table 1, the effect of cooling is greater than that of overpressure on the effective stress reduction, being twice that of the decrease in the effective stress within the aquifer and more than triple that on the interface of the aquifer–caprock, especially in the horizontal direction. Fig. 3 shows more precisely the development of the horizontal effective stress of a section 5 m from the well, with an enlargement in the stress and temperature distribution of the first 25 m of the caprock. After 10 days of low temperature injection, locations above the caprock–aquifer interface were subjected to a 23 ◦ C decrease in temperature. The horizontal effective stress decreased approximately 2 MPa at that interface, whereas it was reduced by 1.2 MPa in the aquifer in which temperature decreased 30 ◦ C. This remarkable sharp stress reduction in the caprock with respect to that of the aquifer was found to depend on the ratio of the thermoEα between the caprock and aquifer, mechanical coupled terms 1− 2ν which will be investigated later on in the Influence of thermal– mechanical coupling parameters section. The upward movement of the cooling continued from the aquifer to the caprock, reducing the stress as shown in the enlargement of Fig. 3. After 100 days, the stress variation remained unchanged at the interface, where the temperature had dropped to the injection temperature. As shown in the enlargement of Fig. 3, the caprock underwent a stress reduction as the temperature decrease spread upward and tended to a steady state. The stress reduction exhibited a negligible delay with respect to the cooling front, and both had a nearly identical trend. This corresponds to the

C. Li, L. Laloui / Geomechanics for Energy and the Environment 11 (2017) 28–41

33

Table 1 Material parameters for the reference model. Thermal parameters

Symbol

Unit

Seal

Aquifer

Underburden

Solid thermal conductivity Water thermal conductivity CO2 thermal conductivity Solid specific heat capacity Water specific heat capacity Solid thermal expansion coefficient Water compressibility Water thermal expansion coefficient

λs λw λa

W/(m K) W/(m K) W/(m K) J/(kg K) J/(kg K) K−1 Pa−1 K−1

1.50 0.67 0.08 950 4183 1.0 × 10−5 3.3 × 10−10 4.5 × 10−5

2.50 0.67 0.08 850 4183 1.0 × 10−5 3.3 × 10−10 4.5 × 10−5

1.50 0.67 0.08 950 4183 1.5 × 10−5 3.3 × 10−10 4.5 × 10−5

kint krc kr w m Pr n0

τ

m2 – – – MPa – –

1 × 10−18 Sc6 6 Sw 0.80 0.60 0.01 0.50

1 × 10−13 Sc3 3 Sw 0.50 0.02 0.10 0.50

1 × 10−18 Sc6 6 Sw 0.80 0.60 0.01 0.50

ρs ρw

kg/m3 kg/m3

2700 1000

2400 1000

2700 1000

E

GPa – –

5.0 0.30 0.60

2.5 0.30 0.60

5.0 0.30 0.60

cp , s cp,w

αs κT βw

Flow parameters Intrinsic permeability CO2 relative permeability Water relative permeability Van Genuchten parameter Van Genuchten parameter Initial porosity Tortuosity Other parameters Solid specific mass Water specific mass Mechanical parameters Young modulus Poisson ratio Initial stress factor

ν K0

Fig. 2. Results of a hydromechanical simulation and thermo-hydro-mechanical simulation for (a), (e) Horizontal effective stress variation; (b), (f) vertical effective stress variation; (c), (g) horizontal displacement and (d), (h) vertical displacement at a section 5 m from the injection well. Displacement is positive in upward and radially outward directions.

proportionality between the two terms dσ ′ and dT , as indicated in Eq. (24), revealing that injection-induced pressurization has been compensated for by the cooling effect in the caprock and that the thermal contraction had a major influence on the stress reduction there. In this study, the Joule–Thomson effect has not been observed as the injection pressure raises not fast enough to form the adiabatic condition, i.e., energy has time to dissipate. It is expected that

the effect will be found significant in the injection scenario that the permeability is lower and injection rate is considerably higher. 6. Influence of the thermal transport properties The thermal conduction was dominant in the caprock as observed in the simulation. It was therefore important to investigate the role of the thermal diffusivity in controlling the heat

34

C. Li, L. Laloui / Geomechanics for Energy and the Environment 11 (2017) 28–41

Fig. 3. Horizontal effective stress variations at a section 5 m from the injection well with an enlargement of the stress variation and temperature in the first 25 m in the caprock.

Fig. 4. Effect of the thermal diffusivity on the horizontal effective stress variations and temperature profiles for three diffusivities at t = 300 days at a section 5 m from the injection well.

transfer process. Experimental studies on the thermal conductivity of shales show that the conductivity can vary from 0.7 to 2.1 W/(m K),34,35 whereas the specific capacity is in the range 850–1100 J/(kg K)35–38 and the density is measured from 2400– 2700 kg/m3 , as shown in Ref. 39. The diffusivity is the ratio between the thermal conductivity and the production of the density and specific heat capacity, as shown in Eq. (15). We suggest three thermal diffusivities for the caprock: D = 3.35E −7(λ = 0.95, cp = 1050, ρ = 2700), 5.85E −7(λ = 1.5, cp = 950, ρ = 2700), and

1.03E −6 m2 /s (λ = 2.1, cp = 850, ρ = 2400). Note that the highest diffusivity was calculated from the maximum combination of the experiment data, whereas a lower density was chosen, which changed the initial stress state of the caprock. The three simulations described for Fig. 4 correspond to a higher thermal expansion coefficient of the caprock, αc = 2.0E −5, to better illustrate the effect. As shown in the temperature profile of Fig. 4(d), the cooling front moved faster upwards, with a higher thermal diffusivity. Correspondingly, in the presence of restricted

C. Li, L. Laloui / Geomechanics for Energy and the Environment 11 (2017) 28–41

deformation, the horizontal effective stress was reduced according to the distribution of the temperature decrease. For the simulation with the highest value of thermal diffusivity, D = 1.03E − 6, the cooling had an impact on the first 15 m of the caprock for a period of 300 days, which was twice that of D = 3.35E − 7. The acceleration in the stress reduction upon cooling can be clearly observed with the increase in the diffusivity in Fig. 4(a)–(c), for which the rate was proportional to the value of the thermal diffusivity. Therefore, a prediction of the thermal transfer can be helpful to foresee a cooling breakthrough of the caprock. The heat transfer inside the caprock was governed by the thermal conduction. Based on that observation, we propose to approximate the vertical movement of the cooling front in a section of the caprock by solving a simplified 1D heat transfer equation as in Eq. (15). As seen in Fig. 4, the reduction in the effective stress mostly followed the temperature distribution in the caprock. Therefore, the temperature evolution could be used as an estimation tool to evaluate the time and distance that were affected by the cooling-induced stress reduction. If the aquifer–caprock interface is subjected to an injection temperature Tinj = 300 K at an instant t = 0 and the far side of the vertical section remains at the initial temperature of T0 = 330 K, the analytical solution to Eq. (15) can be used to estimate the vertical temperature distribution: T (z , t ) = Tinj − T0 erfc

(

)

(

z

√ 2 Dt

)

+ T0 .

(26)

The Complementary Error Function erfc was solved using the built-in MATLAB function. The results in Fig. 5 correspond to the three simulations of Fig. 4. The analytical calculation was in excellent agreement with the numerical simulation on the temperature evolution for simulation times of 100 days and 300 days. This demonstrates that the heat transfer in the vertical direction was governed by heat conduction. In addition, we calculated the temperature distribution for another 4 years after the injection. It can be observed that the cooling transfer was initially ( rapid ) but slowed quickly later on, following the nature of erfc

√1 2 Dt

.

The cooling front remained at the vertical position in the caprock over the long term. The calculation also confirmed the numerical observation of Vilarrasa et al.9 for their long term simulation. Note that this solution is useful for the zone close to the injection well, where the temperature of the caprock base attains the injection temperature quickly after the beginning of injection. For the far zone, in which only a minor drop of temperature was found, Tinj refers to the temperature at the aquifer–caprock interface. 7. Influence of the thermo-mechanical coupling parameters The thermal expansion coefficient governs the level of shrinkage of rocks upon cooling; a small literature has addressed this aspect for CO2 problems, and a value of unity was usually considered in those studies. Only Vilarrasa et al.9 found that the caprock stability is reduced in a normal faulting system from one of their simulation scenarios where the thermal expansion of solid phase is neglected. Cooling-induced stress reduction is sensitive to the thermal expansion coefficient, as shown in Eq. (23). Experimental investigations have shown that the coefficient differs from unity and is in the range of 1.0E − 5 to 2.0E − 5.40,41 The quantitative influence of the coefficient should be investigated. Fig. 6 shows graphs of the horizontal stress reduction and displacement for the three thermal expansion coefficients for the caprock, αc = αa , 1.5αa and 2αa , with αa = 1.0E − 5 in the three cases. For the sake of clarity, because its overall behaviour varied slightly in the later period of injection, as shown in Fig. 2(e), the horizontal stress variation is only shown at 100 days, whereas the displacement is shown at 4 times in the first year. Fig. 6(d)– (f) show that the horizontal displacement increased very slightly

35

Fig. 5. Analytical predictions (black lines) of the temperature profiles with a comparison to the numerical results (red lines) for t = 100 days and t = 300 days at a section 5 m from the injection well. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

in the negative direction with an increase in the thermal expansion coefficient. The displacement was still very constrained, as mentioned previously. Consequently, the increase in the thermal expansion coefficient can only trigger a severe reduction in stress under such deformation constrained conditions. The magnitude of this effective stress reduction is proportional to the value of the thermal expansion coefficient for a given stiffness of rock. Note that for the case where the thermal expansion coefficients of the caprock and aquifer are equal, αc = αa , the horizontal stress drops more in the caprock than in the aquifer because the caprock has a stiffness twice that of the aquifer. This observation is confirmed by Eq. (22); the cooling-induced stress reduction is indeed controlled Eα by the combination of three hydro-mechanical parameters as 1− 2ν for a same temperature decrease. At the aquifer–caprock interface, c αc the ratio R = 1E− / Ea αa could indicate the stress reduction 2ν 1−2ν difference between the point in the caprock and that in the aquifer, if both are subjected to the same temperature decrease. The stiffness of the caprock was thus investigated to discrimEα of the caprock and aquifer. As shown in inate the effect of 1− 2ν the Fig. 7, three caprock stiffnesses, Ec = 1, 2 and 5 GPa,were considered for a specified thermal expansion coefficient αc = 1.5αa , Ea = 2.5 GPa for the aquifer, and ν = 0.3 for both caprock and aquifer. The results are shown for 100 days and 300 days after the beginning of injection, where the temperature at the interface had already reached the injection temperature. For all three cases, the vertical effective stress experienced a similar magnitude in reduction, which was not significantly influenced by variations in Young’s modulus, whereas Young’s modulus had a major impact on the horizontal effective stress reduction. For the lowest value of c αc Ec = 1 GPa, the ratio R = 1E− / Ea αa between the caprock and 2ν 1−2ν aquifer was small: R = 0.6. The decreased horizontal effective stress in the caprock upon cooling was therefore lower than that in the aquifer. The decrease in the vertical effective stress was higher than the horizontal one with this combination inside the caprock. In the case that the stiffness of the caprock was equal to Ec = 2 GPa, as shown in Fig. 7(b), which is, however, smaller than Eα that of the aquifer, 1− of the caprock became slightly higher due 2ν to a higher thermal expansion coefficient, leading to R = 1.2. The

36

C. Li, L. Laloui / Geomechanics for Energy and the Environment 11 (2017) 28–41

Fig. 6. Effect of the thermal expansion coefficient of the caprock on the horizontal stress variations and the horizontal displacement at 5 m from the injection well. Displacement is positive in upward and radially outward directions.

Fig. 7. Effect of the Young modulus of the caprock on the vertical and horizontal effective stress variations at 5 m from the injection well. Displacement is positive in upward and radially outward directions.

decrease in horizontal effective stress was thereby lightly greater in the caprock than that in the aquifer, which was less than the effective stress reduction of the vertical one at the caprock–aquifer interface. As long as the ratio was high, R = 3 as illustrated in Fig. 7(c), the horizontal stress experienced a major reduction and was greater than the decrease in the vertical effective stress, and the sharp difference in the stress reduction was observed.

8. Influence of the hydro-mechanical coupling parameter The Biot coefficient b governs the percentage of stress reduction which is caused by the overpressure. For a soft and unconsolidated rock, b = 1, and for a stiff rock, 0.5 < b < 1. The Biot coefficient is assumed to be 1 in current numerical studies,7,12,42 but the coefficient usually differs from 1 for sandstones from experimental

C. Li, L. Laloui / Geomechanics for Energy and the Environment 11 (2017) 28–41

studies.43–46 Note that a Biot coefficient of 0.5 can reduce the overpressure effect by half. A recent study has shown that a very low porous shale has a Biot coefficient close to 0.9.47 Fig. 8 illustrates the effect of the Biot coefficient on the vertical effective stress variation 1σv′ and the vertical displacement. Aquifers are usually represented as a sandstone-like material, for which the solid compressibility is comparable to the matrix compressibility, leading to a lower value of the Biot coefficient. A comparison of the graphs in Fig. 8(a) and (c) reveals that the Biot coefficient of the aquifer had a moderate influence on the vertical displacement and a minor influence on the stress reduction. Recalling Eq. (18), the effect of the overpressure on the effective stress—the hydromechanical coupling effect is reduced with a lower value of the Biot coefficient, whereupon the aquifer undergoes a smaller order of magnitude of expansion. The displacement reduces 30% at 10 days when the Biot coefficient of the aquifer changes from 1.0 to 0.6. Accordingly, the caprock undergoes the same reduction because of displacement continuity. Because the initial expansion driven by overpressure was reduced with a lower Biot coefficient, whereas the subsequent cooling continued to have a dominant impact on deformation, a negative vertical displacement was eventually observed at 300 days. This suggests that the caprock may induce settlement for long term cooling when the hydromechanical coupling term becomes weak. The results in Fig. 8(b)–(d) demonstrate that the consideration of the caprock’s Biot coefficient had a negligible effect on the caprock deformation. The level of displacement in the lower portion of the caprock was essentially controlled by the expansion of the aquifer, which implies that it was primarily influenced by the consideration of the Biot coefficient of the aquifer. The subsequent reduction in the vertical stress was solely affected by cooling, not by the hydromechanical coupling. 9. Discussion of the failure threshold due to the cooling injection of CO2 Despite numerical studies of CO2 injection that covered a series of material properties, as reviewed by Bao et al.,11 attempts to investigate the links between thermal-hydro-mechanical coupled parameters and thermal properties and the failure threshold are missing. Fig. 9 displays the Mohr circle movement of simulations with different thermal–mechanical coupled coefficient settings of caprock. Five Mohr circles after 100 days of injection for ratios c αc R = 1E− / Ea αa ranging from 0.6 to 4 are synthesized in Fig. 9. The 2ν 1−2ν decreases in the vertical effective stresses in all cases were similar in magnitude, although the stress decreases in the horizontal direction were quite different from one case to another. This led to a movement of all the Mohr circles towards the failure line with the shrinkage or expansion of the circles. It was observed that the circles were smaller than the initial one for R = 0.6 and R = 1.2. We may thus recover the conclusion of Vilarrasa et al.,10 in which a softer caprock was involved, that the failure potential of caprock is reduced and slightly reduced, respectively for R = 0.6 and R = 1.2. However, the Mohr circle expanded as the ratio R increased, and a pure left-shifting of the circle for R = 2 and a proportional increase in the radius of the Mohr circles with the increase in the ratio R were observed, which raised the failure potential. In the scenario considered in this study, R = 4, the Mohr circle expanded significantly and was very close to the failure line. This also confirms the observation of Preisig and Prévost,7 in which the stiffness of the caprock was 3 times higher than that of the aquifer, and both had equal thermal expansion coefficients and Poisson’s ratios. Fig. 9 reveals that, when subjected to low temperature CO2 injection, the caprock stability can be compromised with a stiffer caprock and a higher thermal expansion coefficient. In addition,

37

the failure potential is sensitive to the ratio R of the thermal– mechanical parameter combination of the Young’s modulus, Poisson’s ratio and thermal expansion coefficient. Given the involvement of a highly stiff and thermal-mechanically sensitive caprock i.e., a higher caprock Young’s Modulus and thermal coefficient, the caprock stability is subjected to a high failure potential and may eventually reach the failure threshold if confronted by a locally weak fractured zone in the caprock. The Mohr circles in Fig. 10 mark the effect of the Biot coefficient on the failure threshold. The four cases described for Fig. 8 are summarized here. All the simulations had ratios of R = 3, and the Mohr circles thereby moved towards the failure lines with expansion at 100 days after injection. Despite the fact that the initial Mohr circle sizes were different because of the presence of the Biot coefficient in the definition of the effective stress in Eq. (18), the changes in size were quite similar for all cases. This is consistent with aforementioned analysis that the effective stress reduction in the caprock is primarily induced upon thermal–mechanical coupling instead of overpressure, whereas the Biot coefficient is responsible for the overpressure-induced stress variations. Nevertheless, one can observe that the failure potential for the lowest Biot coefficient set of the caprock and the aquifer (red circles in Fig. 10) was slightly larger than the other cases, as represented by the slightly higher slope of the red tangent line compared to that of the blue line. This explains that the thermal effect affected the caprock stability more strongly when the overpressure effect was attenuated. Fig. 11 plots the evolution of the Mohr circle with two thermal diffusivities for the case of R = 4. The results presented in the remainder of the figure were obtained at 5 m above (at 995 m depth) the caprock–aquifer interface to clearly show the heat transfer. The superposition of the blue and green dashed circles shows that the stress states after 10 days of injection were equal for the cases of D = 3.35E − 7 and D = 5.85E − 7. As analysed previously, the overpressure was dominant at the beginning of injection, whereupon both circles were identical with almost no influence from the thermal diffusivity. The difference with the comparison to the case of D = 1.03E − 6 on the caprock’s density, which was lower than in the other cases, remained, the Mohr circle thereby was on the left initially, and the subsequent shift was applied. After 300 days of injection, the cooling front propagated at different rates according to the diffusivities, and the temperature in the case of the lowest D = 3.35E − 7 remained approximately 4 ◦ C higher than that in the case of D = 5.85E − 7 and approximately 7 ◦ C higher than for D = 1.03E − 6 (recall Fig. 4). Consequently, the caprock exhibited a major failure potential because the red circle had the greatest expansion, even with its initially small Mohr circle. During significant cooling, the horizontal stress decreased rapidly, and the size of the Mohr circle increased rapidly. The Mohr circle for D = 5.85E − 7 was also an expansion but with a smaller magnitude, whereas the decrease in size of the Mohr circle for D = 3.35E − 7 eventually became an increase as time increased. Based on this observation, the caprock with a lower thermal diffusivity experienced less cooling, and the mechanical stability was less vulnerable for this particular moment. Nevertheless the failure potential remained the same, and a delay existed only between the higher and lower diffusivities. 10. Conclusions We have performed series of fully coupled numerical simulations of CO2 injection into a reservoir. The results are presented for the coupled parameter combination ratios R = 0.6 to R = 4.0, which represent the range of parameters that are suitable for materials associated with CO2 injection projects. Previous studies9,10 showed that low temperature injection may increase caprock stability (for the case R < 1) . On the other hand, when a stiffer

38

C. Li, L. Laloui / Geomechanics for Energy and the Environment 11 (2017) 28–41

Fig. 8. Effect of the Biot coefficient on the vertical effect stress variations and vertical displacement at 5 m from the injection well. Displacement is positive in upward and radially outward directions.

Fig. 9. Influence of the material property combination ratio R on the Mohr circles after 100 days of injection and the Mohr circle of the initial condition at the aquifer–caprock interface and 5 m from the injection well.

caprock with a higher thermal expansion coefficient (for the case R ≥ 4 ), as in the study of Preisig and Prévost,7 a failure is attained with an additional temperature decrease. The present simulations are consistent with the earlier predictions and also emphasize that a caprock with a higher ratio R nevertheless has a high potential of failure, deteriorating the mechanical stability. The current numerical analysis was conducted to assess aquifer performance and surface movement with a Biot coefficient equal to unity and without considering the thermal diffusive properties. With an investigation of experimental evidence, the present results complement earlier studies and extends them in a range of geophysical conditions in which the coupled parameters and heat

transfer parameters are substantially validated in practice. Considering the Biot coefficient of an aquifer is essential for properly estimating the expansion of the aquifer and eventually the surface movement. However, considering the caprock Biot coefficient did not influence the results because HM coupling processes did not occur in the caprock in this study. It is necessary to investigate if a higher permeable caprock is involved in the problem because the overpressure may penetrate, which will lead to the occurrence of HM processes in the caprock. In view of the thermal processes, the heat transfer and thermal-induced stress reduction in the caprock can be evaluated following 1D thermal flow theory; the tempera( ) ture decreases with time according to the nature of erfc

√1

t

. In

C. Li, L. Laloui / Geomechanics for Energy and the Environment 11 (2017) 28–41

39

Fig. 10. Influence of the Biot coefficients on the Mohr circles after 100 days of injection and the Mohr circle of the initial condition at the aquifer–caprock interface and 5 m from the injection well. The blue and red lines are tangent lines to the blue and red Mohr circles in solid, respectively. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 11. Influence of the caprock diffusivities on the Mohr circles at 5 m farther in the caprock and 5 m from the injection well. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

addition, the thermal diffusive coefficient plays a major role in the time delay to prevent the failure of the caprock. The interaction between the low temperature injection of CO2 and the sealing caprock of a deep aquifer has been investigated using parametric studies. The study revealed the relatively high importance of the low temperature injection in relation to the caprock stability, which are influenced by the coupled nature of the thermal, hydraulic, and mechanical phenomena that occur. It has been shown that heat transfer and the thermal–mechanical and hydro-mechanical parameters have a very strong influence on the mechanical stability of caprock. The high degree of coupled complex phenomena observed in this study led to apparently paradoxical results such that a change

in the coupled parameter combination (Young’s modulus and thermal expansion coefficient) may reduce or improve caprock stability, and a change in the heat diffusion coefficient may accelerate or delay shearing failure. For this reason, it is risky to extrapolate effects with a set of parameters to other situations in which the materials may have different properties or geometrical configurations are changed. For the same reason, it is necessary to exercise special caution when averaging material properties over a wide range of values, especially when handling the coupled terms of the problem. Acknowledgements The authors gratefully thank Prof. Charlier and Dr. Collin (University of Liège) for providing the LAGAMINE finite element code.

40

C. Li, L. Laloui / Geomechanics for Energy and the Environment 11 (2017) 28–41

Appendix. Equation of state of peng and robinson The properties of gases and liquids strongly depend on the environmental pressure and temperature conditions. The equation of state that was developed by Peng and Robinson21 is used to describe the relationship of the volume V, pressure P and temperature T of real fluids as presented in the following: P =

RT V −b

−

a

(A.1)

V (V + b) + b (V − b)

where R a

universal gas constant ac α

α

⌈

ac

0.457235R2 Tc2 /Pc 0.37464 + 1.54226ω − 0.26992ω2 if ω ≤ 0.49 0.379642 + 1.48503ω − 0.164423ω2 + 0.016666ω3 if ω > 0.49 Pitzer acentric factor, ω = 0.239 critical pressure, Pc = 7.39 MPa critical temperature, Tc = 304.25 K

κ κ

ω Pc Tc

1+κ 1−

(

√

T /Tc

)⌉

Eq. (A.1) can be rewritten in a cubic polynomial form for the compressibility Z as Z 3 − (1 − B) Z 2 + A − 3B2 − 2B Z − AB − B2 − B3 = 0

(

)

(

)

(A.2)

where A= B= Z =

aP R2 T 2 bP RT PV RT

(A.3)

.

The solution to the cubic equation Eq. (A.2) yields one or three real roots for the compressibility factor Z , depending on the number of phases in the system. The largest root is for the compressibility Z of the gaseous phase, whereas the smallest positive root corresponds to that of the liquid in a two-phase system. References 1. Bachu S. Sequestration of CO2 in geological media: criteria and approach for site selection in response to climate change. Energy Convers Manage. 2000;41(9):953–970. 2. Bryant E. Climate Process and Change. Cambridge, UK: Cambridge University Press; 1997. 3. Metz B, Davidson O, de Coninck H, Loos M, Meyer L, IPCC Special Report on Carbon Dioxide Capture and Storage; 2005. 4. Li C, Laloui L. Coupled multiphase thermo-hydro-mechanical analysis of supercritical CO2 injection: Benchmark for the In Salah surface uplift problem. Int J Greenhouse Gas Control. 2016;51:394–408. 5. Rutqvist J. The geomechanics of CO2 storage in deep sedimentary formations. Geotech Geol Eng. 2012;30(3):525–551. 6. Bissell RC, Vasco DW, Atbi M, Hamdani M, Okwelegbe M, Goldwater MH. A full field simulation of the in Salah gas production and CO2 storage project using a coupled geo-mechanical and thermal fluid flow simulator. Energy Procedia. 2011;4:3290–3297. 7. Preisig M, Prévost JH. Coupled multi-phase thermo-poromechanical effects. Case study: CO2 injection at In Salah, Algeria. Int J Greenhouse Gas Control. 2011;5(4):1055–1064. 8. Gor GY, Elliot TR, Prévost JH. Effects of thermal stresses on caprock integrity during CO2 storage. Int J Greenhouse Gas Control. 2013;12:300–309. 9. Vilarrasa V, Olivella S, Carrera J, Rutqvist J. Long term impacts of cold CO2 injection on the caprock integrity. Int J Greenhouse Gas Control. 2014;24:1–13. 10. Vilarrasa V, Silva O, Carrera J, Olivella S. Liquid CO2 injection for geological storage in deep saline aquifers. Int J Greenhouse Gas Control. 2013;14:84–96.

11. Bao J, Hou Z, Fang Y, Ren H, Lin G. Uncertainty quantification for evaluating impacts of caprock and reservoir properties on pressure buildup and ground surface displacement during geological CO 2 sequestration. Greenhouse Gases Sci Technol. 2013;3(5):338–358. 12. Vilarrasa V, Bolster D, Olivella S, Carrera J. Coupled hydromechanical modeling of CO2 sequestration in deep saline aquifers. Int J Greenhouse Gas Control. 2010;4(6):910–919. 13. Bowen RM. Compressible porous media models by use of the theory of mixtures. Internat J Engrg Sci. 1982;20(6):697–735. 14. Panday S, Corapcioglu MY. Reservoir transport equations by compositional approach. Transp Porous Media. 1989;4(4):369–393. 15. Collin F. Couplages Thermo-Hydro-Mécaniques Dans Les Sols et Les Roches Tendres Partiellement Saturés. Belgium: University of Liège; 2003. 16. Lu M, Connell L. Non-isothermal flow of carbon dioxide in injection wells during geological storage. Int J Greenhouse Gas Control. 2008;2(2):248–258. 17. Charlier R, Radu J-P, Collin F. Numerical modelling of coupled transient phenomena. Rev Fr Génie Civ. 2011;5(6):719–741. 18. Lewis RW, Schrefler BA. The Finite Element Method in the Static and Dynamic Deformation and Consolidation of Porous Media. 2nd ed., Chichester, UK: Wiley; 1998. 19. Thomas HR, King SD. A non-linear, two-dimensional, potential-based analysis of coupled heat and mass transfer in a porous medium. Internat J Numer Methods Engrg. 1994;37(21):3707–3722. 20. Fenghour A, Wakeham WA, Vesovic V. The viscosity of carbon dioxide. J Phys Chem Ref Data. 1998;27(1):31. 21. Peng D-Y, Robinson DB. A new two-constant equation of state. Ind Eng Chem Fundam. 1976;15(1):59–64. 22. Lin CK. Algorithm for determining optimum sequestration depth of CO 2 trapped by residual gas and solubility trapping mechanisms in a deep saline formation. Geofluids. 2008;8(4):333–343. 23. Pruess K, García J. Multiphase flow dynamics during CO2 disposal into saline aquifers. Environ Geol. 2002;42(2–3):282–295. 24. Crovetto R. Evaluation of Solubility CO2-H2O. J Phys Chem Ref Data. 1991;20(3):575– 589. 25. Span R, Wagner W. A new equation of state for carbon dioxide covering the fluid region from the triple-point temperature to 1100 K at pressures up to 800 MPa. J Phys Chem Ref Data. 1996;25(6):1509. 26. Jamieson DT, Tudhope JS, Morris R, Cartwright G. Physical properties of sea water solutions: heat capacity. Desalination. 1969;7(1):23–30. 27. Wilson JV, Approximations for physical properties of sea salt solutions. Oak Ridge, Tennessee; 1973. 28. Rutqvist J, Tsang C-F. A study of caprock hydromechanical changes associated with CO 2 -injection into a brine formation. Environ Geol. 2002;42(2–3):296– 305. 29. Selvadurai APS. Heave of a surficial rock layer due to pressures generated by injected fluids. Geophys Res Lett. 2009;36(14):L14302. 30. Shi J-Q, Sinayuc C, Durucan S, Korre A. Assessment of carbon dioxide plume behaviour within the storage reservoir and the lower caprock around the KB502 injection well at In Salah. Int J Greenhouse Gas Control. 2012;7:115–126. 31. Jaeger JC, Cook NGW, Zimmerman RW. Fundamentals of Rock Mechanics. Blackwell Pub.; 2007. 32. Selvadurai APS, Mechanics of an embedded caprock layer during pressurization of a CO2 storage reservoir. In: Proceedings of ComGeo 09: Computational Geomechanics I; 2008, pp. 466–475. 33. Li C, Barès P, Laloui L. A hydromechanical approach to assess CO2 injectioninduced surface uplift and caprock deflection. Geomech Energy Environ. 2015;4:51–60. 34. Blackwell D, Steele J. Thermal conductivity of sedimentary rocks: measurement and significance. Naeser N, McCulloh T, eds. in: Thermal History of Sedimentary Basins SE - 2. New York: Springer; 1989:13–36. 35. Eppelbaum L, Kutasov I, Pilchin A. Thermal properties of rocks and density of fluids. in: Applied Geothermics. Berlin Heidelberg: Springer; 2014:99–149. 36. Pasquale V, Gola G, Chiozzi P, Verdoya M. Thermophysical properties of the Po Basin rocks. Geophys J Int. 2011;186(1):69–81. 37. Waples DW, Waples JS. A review and evaluation of specific heat capacities of rocks, minerals, and subsurface fluids. Part 2: Fluids and porous rocks. Nat Resour Res. 2004;13(2):123–130. 38. Robertson EC, Thermal Properties of Rocks, US Dep. Inter. Geol. Surv.; 1988, pp. 88–441. 39. Shalabi FI, Cording EJ, Al-Hattamleh OH. Estimation of rock engineering properties using hardness tests. Eng Geol. 2007;90(3–4):138–147. 40. Jobmann M, Polster M. The response of Opalinus clay due to heating: A combined analysis of in situ measurements, laboratory investigations and numerical calculations. Phys Chem Earth Parts A/B/C . 2007;32(8–14):929–936. 41. Gilliam TM, Morgan TM, Shale: Measurement of Thermal Properties (No. ORNL/TM-10499). Oak Ridge National Lab. TN (USA). Oak Ridge, Tennessee; Jul 1987. 42. Orlic B. Some geomechanical aspects of geological CO2 sequestration. KSCE J Civ Eng. 2009;13(4):225–232.

C. Li, L. Laloui / Geomechanics for Energy and the Environment 11 (2017) 28–41 43. Vidal-Gilbert S, Tenthorey E, Dewhurst D, Ennis-King J, Van Ruth P, Hillis R. Geomechanical analysis of the Naylor Field, Otway Basin, Australia: Implications for CO2 injection and storage. Int J Greenhouse Gas Control. 2010;4(5):827– 839. 44. Nur A, Byerlee JD. An exact effective stress law for elastic deformation of rock with fluids. J Geophys Res. 1971;76(26):6414. 45. Siggins AF. Saturation, pore pressure and effective stress from sandstone acoustic properties. Geophys Res Lett. 2003;30(2):10–13.

41

46. Bouteca MJ, Bary D, Piau JM, Kessler N, Boisson M, Fourmaintraux D. Contribution of poroelasticity to reservoir engineering: Lab experiments, application to core decompression and implication in HP-HT reservoirs depletion. in: Rock Mechanics in Petroleum Engineering. 2013. 47. Favero V, Ferrari A, Laloui L. Thermo-mechanical volume change behaviour of Opalinus Clay. Int J Rock Mech Min Sci. 2016;90:15–25.