# Coupled Geomechanical-Flow Assessment of CO2 Leakage through Heterogeneous Caprock during CCS

Authors: Guan Woo Kim; Tae Hong Kim; Jiho Lee; Kun Sang Lee
Journal: Advances in Civil Engineering
 (2018)

---

## Abstract

The viability of carbon capture sequestration (CCS) is dependent on the secure storage of CO2 in subsurface geologic formations. Geomechanical failure of caprock is one of the main reasons of CO2 leakage from the storage formations. Through comprehensive assessment on the petrophysical and geomechanical heterogeneities of caprock, it is possible to predict the risk of unexpected caprock failure. To describe the fracture reactivation, the modified Barton–Bandis model is applied. In order to generate hydro-geomechanically heterogeneous fields, the negative correlation between porosity and Young’s modulus/Poisson’s ratio is applied. In comparison with the homogeneous model, effects of heterogeneity are examined in terms of vertical deformation and the amount of leaked CO2. To compare the effects of heterogeneity, heterogeneous models for both geomechanical and petrophysical properties in coupled simulation are designed. After 10-year injection with petrophysically heterogeneous and geomechanically homogeneous caprock, CO2 leakage is larger than that of the homogeneous model. In contrast, heterogeneity of geomechanical properties is shown to mitigate additional escape of CO2. Vertical displacement of every heterogeneous model is larger than homogeneous model. The model with compressive tectonic stress shows much more stable trapping with heterogeneous caprock, but there is possibility of rapid leakage after homogeneous caprock failure.


---

## Body

## 1. Introduction

CO2 sequestration in aquifer is an effective and a verified method to reduce the atmospheric CO2 airborne fraction. The world’s first commercial approach of CO2 storage project is attempted at the Sleipner gas field in the North Sea, and CO2 sequestration in a saline aquifer has been regarded as the feasible technology [1]. For the geological storage of CO2, there should be an impermeable caprock which acts as a physical barrier to prevent the leakage and migration of the buoyant CO2 [2]. One of the considerable risks is the leak of CO2 incurred by a mechanical failure of caprock. Injected CO2 can leak out and flow throughout the opened fractures and faults of caprock; this is an important point in the caprock’s sealing mechanism. Thus, it is important to monitor the geomechanical change of storage formation to evaluate the stability of CO2 injection. In the aftermath of large volume of CO2 injection, fluid pressure in the storage formation will rise, and increase of mechanical stresses will also occur. If injected CO2 induce too high fluid pressure, the induced stresses will exceed the summation of least principal effective stress and tensile strength and affect irreversible mechanical displacements of storage rock, reactivating natural fractures or creating additional fractures. Geological displacements of CO2 in storage formation have been observed at Weyburn CO2 monitoring and storage project [3], and surface uplift caused by increased subsurface pressure has been detected at the In Salah project [4]. These changes could generate new high-permeable zone as leakage path, thereby substantially decreasing the efficiency of CO2 storage.To ensure the feasibility of potential storage site, intensive estimation of formation properties should be preceded [5]. Hurter et al. [6] emphasized the importance of numerical simulation as preinjection studies. Several numerical simulation studies have been conducted to assess the stability of CO2 injection with geomechanical module. Lucier et al. [7] applied a geomechanical analysis of a potential injection zone in the Ohio Valley region to evaluate the potential risks for injection-induced stresses. Rutqvist et al. [8] conducted coupled reservoir-geomechanical simulations to evaluate the shear and tensile failure in order to accurately estimate the potential of fracture reactivation for a geomechanical CO2 storage formation, and they suggested that long-lasting observation is needed to geomechanical changes resulting from the injected fluid pressure. Park et al. [9] evaluated the effects of geomechanical properties on the fracture reactivation by using two-dimensional simulation model, and they also examined the effects of the heterogeneity of geomechanical properties during CO2 storage.For a much more accurate evaluation of geomechanical change, the coupled mechanism of geomechanical and petrophysical properties was considered. Cappa and Rutqvist [10] tried to reveal the interaction between fluid flow and mechanical deformation during CCS and demonstrated how to analyze coupled hydromechanical processes. Khajeh et al. [11] investigated the impact of coupled geomechanical and petrophysical properties including heterogeneity and range of uncertainties. Previous researches have considered this coupling, but they did not focus on the correlation among these properties that could affect the accuracy of numerical modeling.In this study, we constructed a CCS model with heterogeneous caprock and compared the flow and geomechanical responses to CO2 injection with particular focus on the risks to storage security posed by geomechanical deformation. The correlations between geomechanical and petrophysical properties were regarded as a function of porosity, and several numbers of numerical simulations were conducted to evaluate the effects of the heterogeneous caprock during CCS. The amount of stored CO2 and vertical deformation of caprock within a two-dimensional formation models were compared for several rock properties including parameters such as Poisson’s ratio, Young’s modulus, porosity, and permeability. These properties were distributed in randomly generated heterogeneity conditions. The effects of heterogeneity on caprock were revealed by comparing the average values of heterogeneous model simulation results with homogeneous caprock model. Considering correlation between geomechanical and petrophysical properties can help assessing formation stability and mitigating geological risk. Additionally, insitu tectonic stress is known to play an important role in fracture propagation [12, 13]. Horizontal tectonic stress also influences the porosity and permeability of formation [14]. However, there appears to be a lack of information regarding the influence of in situ tectonic stress on long-time CO2 storage performance. In order to address this issue, an investigation of heterogeneous permeability fields including additional horizontal stress has been conducted.


## 2. Theoretical Background

### 2.1. Fracture Reactivation

The overburden stress, known as the total stress, induced by external loading, is sustained by the rock matrix and the formation fluid pressure. To evaluate the fluid pressure required to fracture reactivation in caprock, Terzaghi [15] suggested the law of effective stress as follows:(1)σ=σ′+p0,where σ is the total stress (kPa), σ′ is the effective stress (kPa), and p0 is the fluid pressure (kPa).The tensile failure potential is affected by applying the assumption that tensile fracture propagation can occur when the fluid pressure exceeds the summation of least principal effective stress and tensile strength [16]:(2)pfc≥σ3+σt,where pfc is the critical fluid pressure (kPa), σ3 is the least principal effective stress (kPa), and σt is the tensile strength (kPa). Most of sedimentary rocks have a relatively low tensile strength, typically few MPa or less. Fjaer et al. [17] indicated that it is a standard approximation for numerous applications that the tensile strength is zero. In the case that the pressure build-up by injected fluid exceeds the least principal effective stress, the fracture propagates in the direction being perpendicular to least principal effective stress. This process is similar to the tensile reactivation of natural fractures.


### 2.2. Modification of Fracture Permeability

A hydromechanical fracture permeability model was based on the modified Barton–Bandis model, as shown in Figure1 [18]. Initially, effective stress (σ′) is much larger than fluid pressure (point A), and there are no activated fractures in the rock matrix. The effective stress is decreased along the CO2 injection (path AB). On this path, the range of fracture permeability is small, and the reactivation is reversible. When the effective stress is below point B, the fracture opens abruptly, and the fracture permeability increases drastically from the initial value to maximum value (khf) on path BC. The threshold value (point B) is named the fracture opening stress (frs). Although the effective stress decreases consistently, the fracture permeability keeps its maximum value (path DE). Effective stress increment and the fracture aperture reduction due to decline of fluid pressure lead to the decrease of fracture permeability. However, fracture permeability change does not occur on the previous path BC. Instead, it follows a new hyperbolic curve (path FG) because the fracture cannot be closed completely. Fracture closure permeability is always bigger than the residual value of fracture closure permeability (krcf). If the increment of fluid pressure induces another reduction of the effective stress, fracture permeability will increase along the path GFED to the maximum fracture permeability. Fracture permeability (kf) can be calculated as follows [18, 19]:(3)kf=kccfee04,e=e0−Vj,Vj=σ′kni+σ′/Vm,Vm=e01−krcfkccf1/4,where kccf is the fracture closure permeability (md), e0 is the initial fracture aperture (m), Vj is the joint closure under a normal fracture effective stress, kni is the initial normal fracture stiffness (kPa/m), and Vm is the maximum fracture closure (m).Figure 1
Modified Barton–Bandis model for the change in fracture permeability.


### 2.3. Permeability-Porosity Relationship

Walsh and Brace [20] indicated that permeability is regarded as proportional to simple integer powers of the related pore geometry parameters such as hydraulic radius, porosity, and tortuosity. These parameters are normally regarded to be related to each other through the power-law relationship (k∝ϕα), leading to power-law dependence of permeability on porosity [21]. Experimental evidence indicates that power-law exponent does not always match porosity changes. However, for elastic or plastic compaction, chemical alteration, or microcracking, it is usually regarded that there is a power-law relationship between these parameters [22].


### 2.4. Porosity-Geomechanical Properties Relationship

To analyse the petrophysical and geomechanical properties, sedimentary rocks were studied using samples cored from several wells from the Krishna-Godavari basin on the east coast of India [23]. Figure 2 shows best fit linear curves for uniaxial compressive strength, effective porosity, and Young’s modulus (E) of core samples. As could be seen that the linear relationship between Young’s modulus strength, and uniaxial compressive strength has reasonable correlation coefficients (R2 = 0.93). Similarly, the correlation coefficients of the best fit linear curves between uniaxial compressive and effective porosity are high (R2 = 0.89). Using linear equations, it is possible to predict the geomechanical properties of sedimentary rocks from petrophysical properties. The results of the regression analysis indicate that the petrophysical and geomechanical properties have a positive correlation to each other. Figure 3 shows the relationship between the geomechanical properties of E and Poisson’s ratio (ν) [24]. Rock physics templates show the mapping of lines of mineralogical mixtures and constant porosity on the crossplots of geomechanical properties [25]. Solid lines represent the clay content, and dashed lines indicate constant values of porosity. The intersection of these lines is overlaid with points, representing chosen values of porosity, from 0 to 0.20. Each color of dots represents constant clay content, explained by the legend. From each constant porosity line, it can be seen that the increase in porosity will decrease the values of ν and E. This may imply that porosity has negative correlation with ν and E.Figure 2
Linear relationship and correlation coefficient between petrophysical and geomechanical properties of cores samples from wells in the Krishna-Godavari basin. (a) Young’s modulus and uniaxial compressive strength. (b) Uniaxial compressive strength and effective porosity.
(a)
(b)Figure 3
Crossplots of Young's modulus and Poisson’s ratio, showing the degree of clay. Mineralogical substitution occurs between (a) quartz and clay and (b) dolomite and clay.
(a)
(b)


## 2.1. Fracture Reactivation

The overburden stress, known as the total stress, induced by external loading, is sustained by the rock matrix and the formation fluid pressure. To evaluate the fluid pressure required to fracture reactivation in caprock, Terzaghi [15] suggested the law of effective stress as follows:(1)σ=σ′+p0,where σ is the total stress (kPa), σ′ is the effective stress (kPa), and p0 is the fluid pressure (kPa).The tensile failure potential is affected by applying the assumption that tensile fracture propagation can occur when the fluid pressure exceeds the summation of least principal effective stress and tensile strength [16]:(2)pfc≥σ3+σt,where pfc is the critical fluid pressure (kPa), σ3 is the least principal effective stress (kPa), and σt is the tensile strength (kPa). Most of sedimentary rocks have a relatively low tensile strength, typically few MPa or less. Fjaer et al. [17] indicated that it is a standard approximation for numerous applications that the tensile strength is zero. In the case that the pressure build-up by injected fluid exceeds the least principal effective stress, the fracture propagates in the direction being perpendicular to least principal effective stress. This process is similar to the tensile reactivation of natural fractures.


## 2.2. Modification of Fracture Permeability

A hydromechanical fracture permeability model was based on the modified Barton–Bandis model, as shown in Figure1 [18]. Initially, effective stress (σ′) is much larger than fluid pressure (point A), and there are no activated fractures in the rock matrix. The effective stress is decreased along the CO2 injection (path AB). On this path, the range of fracture permeability is small, and the reactivation is reversible. When the effective stress is below point B, the fracture opens abruptly, and the fracture permeability increases drastically from the initial value to maximum value (khf) on path BC. The threshold value (point B) is named the fracture opening stress (frs). Although the effective stress decreases consistently, the fracture permeability keeps its maximum value (path DE). Effective stress increment and the fracture aperture reduction due to decline of fluid pressure lead to the decrease of fracture permeability. However, fracture permeability change does not occur on the previous path BC. Instead, it follows a new hyperbolic curve (path FG) because the fracture cannot be closed completely. Fracture closure permeability is always bigger than the residual value of fracture closure permeability (krcf). If the increment of fluid pressure induces another reduction of the effective stress, fracture permeability will increase along the path GFED to the maximum fracture permeability. Fracture permeability (kf) can be calculated as follows [18, 19]:(3)kf=kccfee04,e=e0−Vj,Vj=σ′kni+σ′/Vm,Vm=e01−krcfkccf1/4,where kccf is the fracture closure permeability (md), e0 is the initial fracture aperture (m), Vj is the joint closure under a normal fracture effective stress, kni is the initial normal fracture stiffness (kPa/m), and Vm is the maximum fracture closure (m).Figure 1
Modified Barton–Bandis model for the change in fracture permeability.


## 2.3. Permeability-Porosity Relationship

Walsh and Brace [20] indicated that permeability is regarded as proportional to simple integer powers of the related pore geometry parameters such as hydraulic radius, porosity, and tortuosity. These parameters are normally regarded to be related to each other through the power-law relationship (k∝ϕα), leading to power-law dependence of permeability on porosity [21]. Experimental evidence indicates that power-law exponent does not always match porosity changes. However, for elastic or plastic compaction, chemical alteration, or microcracking, it is usually regarded that there is a power-law relationship between these parameters [22].


## 2.4. Porosity-Geomechanical Properties Relationship

To analyse the petrophysical and geomechanical properties, sedimentary rocks were studied using samples cored from several wells from the Krishna-Godavari basin on the east coast of India [23]. Figure 2 shows best fit linear curves for uniaxial compressive strength, effective porosity, and Young’s modulus (E) of core samples. As could be seen that the linear relationship between Young’s modulus strength, and uniaxial compressive strength has reasonable correlation coefficients (R2 = 0.93). Similarly, the correlation coefficients of the best fit linear curves between uniaxial compressive and effective porosity are high (R2 = 0.89). Using linear equations, it is possible to predict the geomechanical properties of sedimentary rocks from petrophysical properties. The results of the regression analysis indicate that the petrophysical and geomechanical properties have a positive correlation to each other. Figure 3 shows the relationship between the geomechanical properties of E and Poisson’s ratio (ν) [24]. Rock physics templates show the mapping of lines of mineralogical mixtures and constant porosity on the crossplots of geomechanical properties [25]. Solid lines represent the clay content, and dashed lines indicate constant values of porosity. The intersection of these lines is overlaid with points, representing chosen values of porosity, from 0 to 0.20. Each color of dots represents constant clay content, explained by the legend. From each constant porosity line, it can be seen that the increase in porosity will decrease the values of ν and E. This may imply that porosity has negative correlation with ν and E.Figure 2
Linear relationship and correlation coefficient between petrophysical and geomechanical properties of cores samples from wells in the Krishna-Godavari basin. (a) Young’s modulus and uniaxial compressive strength. (b) Uniaxial compressive strength and effective porosity.
(a)
(b)Figure 3
Crossplots of Young's modulus and Poisson’s ratio, showing the degree of clay. Mineralogical substitution occurs between (a) quartz and clay and (b) dolomite and clay.
(a)
(b)


## 3. Numerical Simulation

GEM was used to assess the fracture propagation during the CO2 injection, which is a compositional simulator developed by the Computer Modeling Group [18]. The geomechanics module in GEM was applied to compute the deformation in the rock and the stress change. All the rocks in the matrix were assumed to be elastic during the injection process. To describe the fracture permeability change, dual permeability system and Barton–Bandis model were used. For the in situ stress, the weight of the overburden rock causes the vertical stress. The initial fluid pressure was induced by the hydrostatic gradient of 9.81 MPa/km with an atmospheric pressure of 0.1 MPa from the surface. Constant pressure condition is supposed at the horizontal outer boundary to minimize the boundary effects. All of the simulations were organized in a constant temperature. The static geomechanical and hydraulic properties of the base case are shown in Table 1, with the typical rock types of the permeable caprock and the formations representing sandstone and shale [26].Table 1
Hydraulic and geomechanical properties for homogeneous model.
Properties
Overburden rock
Caprock
Aquifer
Permeability (md)
25
0.0001
100
Porosity
0.21
0.01
0.25
Rock density (kg/m3)
2,260
2,260
2,260
Young’s modulus,E (GPa)
5
5
5
Poisson’s ratio,ν
0.25
0.285
0.25
Maximum fracture permeability,kf  max (md)
—
35
—
Minimum fracture permeability,kf  min (md)
—
5
—


### 3.1. Model Description

In order to assess the effects of the geomechanical rock properties of the fracture reactivation, several simulations were conducted on two-dimensional cross-sectional model. The model considers permeability and porosity of rock matrix, Poisson’s ratio, and Young’s modulus, which are regarded to be the most important properties for the investigation of the geomechanical state. The formation model was constructed, as shown in Figure4, including an aquifer, caprock, and overburden rock. The formation model locates vertically from 1,000 to 1,330 m depth. Model of 2,020 × 20 × 330 m3 was divided into 101 × 1 × 33 grid blocks. Pure CO2 was injected through an injection well at 1,300 m depth and a constant rate of 120,000 m3/day for 10 years. Table 2 shows component properties of CO2 used in the simulation [27]. CO2 existed as a supercritical state at the given temperature and pressure (T = 50°C and p=12,853kPa). At first, the injected CO2 was sequestered under the caprock and then leaked through the fracture in the caprock after the fracture reactivation occurred. During the simulation, the vertical deformation and effective stress were measured.Figure 4
Model description: aquifer, caprock, and overburden. The injector is perforated in bottom three layers of the aquifer.Table 2
CO2 component properties.
Property
Value
p
c (atm)
72.8
T
c (K)
304.2
Molecular weight
44.01Porosity-permeability relationship is calculated using a power law. To generate geostatistical hydrogeomechanically heterogeneous fields, the negative correlation between porosity and Young’s modulus is applied. Poisson’s ratio is negatively related to porosity, too. Table3 shows physical data of geological caprocks from saline CO2 storage in the Paradox Basin [28–30]. Geomechanical properties of simulation values are modeled based on this formation data.Table 3
Physical data of geological seals considered for saline CO2 sequestration in the Paradox Basin.
Parameter
Value
Depth (m)
1,707–1,768
Thickness (m)
9–12
Porosity (%)
1.6–4.1
Vertical permeability,−lnkv (m2)
39.16–36.86
Horizontal permeability,−lnkh (m2)
44.32–44.12(44.21 ± 0.09)
Compressive strength,σ (MPa)
166–311(246.58 ± 54.74)
Poisson’s ratio,ν
0.24–0.33(0.27 ± 0.04)
Young’s modulus,E (MPa)
38,562–56,613(47,315 ± 6,447)
Basin areal extent (km2)
28,490


### 3.2. Heterogeneity Description

In this study, a stochastic simulation method is used for heterogeneous reservoir model. To describe formation heterogeneity, Dykstra–Parsons coefficient (VDP), which is the most common measure of heterogeneity in petroleum industry, is used [31, 32]. It is defined as(4)VDP=k50−k84.1k50,where k50 is the median permeability value and k84.1 is the permeability at 84.1% probability (one standard deviation). The model with same VDP also has the same correlation length (Table 4). Averages of all variables are same as homogeneous model (Table 1). Each model has same mean and variance of permeability, but all of the realizations have their own permeability distribution. The result graphs are drawn from the average of 25 realizations. The number of realizations was chosen after testing the sensitivity of results. Though variations were sometimes observed, 25 were chosen for the final number of realizations for analysis due to computational constraints [11].Table 4
Numerical attributes considered for caprock permeability.
V
DP
Mean (md)
Variance
0.3
0.0001
0.25
0.5
0.0001
1
0.7
0.0001
3To evaluate how heterogeneity considerations for petrophysical and geomechanical properties affect output results, four cases are designed as below:Case A: heterogeneous petrophysical and homogeneous geomechanical propertiesCase B: homogeneous petrophysical and heterogeneous geomechanical propertiesCase C: heterogeneous petrophysical and geomechanical propertiesCase D: homogeneous petrophysical and geomechanical propertiesThe effects of heterogeneity are revealed in comparison to homogeneous model. The results are compared in terms of vertical deformation and the amount of leaked CO2.


## 3.1. Model Description

In order to assess the effects of the geomechanical rock properties of the fracture reactivation, several simulations were conducted on two-dimensional cross-sectional model. The model considers permeability and porosity of rock matrix, Poisson’s ratio, and Young’s modulus, which are regarded to be the most important properties for the investigation of the geomechanical state. The formation model was constructed, as shown in Figure4, including an aquifer, caprock, and overburden rock. The formation model locates vertically from 1,000 to 1,330 m depth. Model of 2,020 × 20 × 330 m3 was divided into 101 × 1 × 33 grid blocks. Pure CO2 was injected through an injection well at 1,300 m depth and a constant rate of 120,000 m3/day for 10 years. Table 2 shows component properties of CO2 used in the simulation [27]. CO2 existed as a supercritical state at the given temperature and pressure (T = 50°C and p=12,853kPa). At first, the injected CO2 was sequestered under the caprock and then leaked through the fracture in the caprock after the fracture reactivation occurred. During the simulation, the vertical deformation and effective stress were measured.Figure 4
Model description: aquifer, caprock, and overburden. The injector is perforated in bottom three layers of the aquifer.Table 2
CO2 component properties.
Property
Value
p
c (atm)
72.8
T
c (K)
304.2
Molecular weight
44.01Porosity-permeability relationship is calculated using a power law. To generate geostatistical hydrogeomechanically heterogeneous fields, the negative correlation between porosity and Young’s modulus is applied. Poisson’s ratio is negatively related to porosity, too. Table3 shows physical data of geological caprocks from saline CO2 storage in the Paradox Basin [28–30]. Geomechanical properties of simulation values are modeled based on this formation data.Table 3
Physical data of geological seals considered for saline CO2 sequestration in the Paradox Basin.
Parameter
Value
Depth (m)
1,707–1,768
Thickness (m)
9–12
Porosity (%)
1.6–4.1
Vertical permeability,−lnkv (m2)
39.16–36.86
Horizontal permeability,−lnkh (m2)
44.32–44.12(44.21 ± 0.09)
Compressive strength,σ (MPa)
166–311(246.58 ± 54.74)
Poisson’s ratio,ν
0.24–0.33(0.27 ± 0.04)
Young’s modulus,E (MPa)
38,562–56,613(47,315 ± 6,447)
Basin areal extent (km2)
28,490


## 3.2. Heterogeneity Description

In this study, a stochastic simulation method is used for heterogeneous reservoir model. To describe formation heterogeneity, Dykstra–Parsons coefficient (VDP), which is the most common measure of heterogeneity in petroleum industry, is used [31, 32]. It is defined as(4)VDP=k50−k84.1k50,where k50 is the median permeability value and k84.1 is the permeability at 84.1% probability (one standard deviation). The model with same VDP also has the same correlation length (Table 4). Averages of all variables are same as homogeneous model (Table 1). Each model has same mean and variance of permeability, but all of the realizations have their own permeability distribution. The result graphs are drawn from the average of 25 realizations. The number of realizations was chosen after testing the sensitivity of results. Though variations were sometimes observed, 25 were chosen for the final number of realizations for analysis due to computational constraints [11].Table 4
Numerical attributes considered for caprock permeability.
V
DP
Mean (md)
Variance
0.3
0.0001
0.25
0.5
0.0001
1
0.7
0.0001
3To evaluate how heterogeneity considerations for petrophysical and geomechanical properties affect output results, four cases are designed as below:Case A: heterogeneous petrophysical and homogeneous geomechanical propertiesCase B: homogeneous petrophysical and heterogeneous geomechanical propertiesCase C: heterogeneous petrophysical and geomechanical propertiesCase D: homogeneous petrophysical and geomechanical propertiesThe effects of heterogeneity are revealed in comparison to homogeneous model. The results are compared in terms of vertical deformation and the amount of leaked CO2.


## 4. Results and Discussion

### 4.1. CO2 Storage

Figure5 shows gas saturation of single realization and CO2 leakage regime. The figure is about Case C which has petrophysically heterogeneous properties and geomechanically heterogeneous properties and Dykstra–Parsons coefficient of 0.7. Increased formation pressure causes fractures in caprock, and CO2 escapes from the aquifer to overburden rock. Figure 6 depicts the cumulative amount of remaining CO2 in aquifer during 10-year injection. These graphs show the total amount of gas sequestered in the aquifer. Each curve is generated by averaging results from 25 heterogeneous models. The rest of gas has escaped to overburden formation. There could be further escape after injection, until the system reaches equilibrium. Figure 6(a) indicates simulation model with petrophysically heterogeneous and geomechanically homogeneous properties. By investigating the results of Case A, it can be seen that heterogeneous model shows less CO2 moles in the aquifer than homogeneous model. Also, the amount of leaked CO2 depends on heterogeneity. In other words, more petrophysically heterogeneous model presents more CO2 leakage. In heterogeneous caprock, there could be relatively weak shale which has high porosity and permeability. Injected fluids tend to flow through high permeable zone and partially increase formation pressure. Concentration of pressure incurs fracture reactivation and earlier leakage of CO2. Figure 6(b) shows simulation model with petrophysically homogeneous and geomechanically heterogeneous properties. In early stage of injection, heterogeneous models show more leakage than homogeneous model. However, as time goes on, the leakage is mitigated compared to homogeneous model. This difference is caused by imbalance of geomechanical properties. Being different from the heterogeneous caprock, the failure of homogeneous caprock occurs abruptly at the center. In heterogeneous caprock, the leakage does not occur above the injector. Initially, CO2 gas leaks through the location region with high Young’s modulus and Poisson’s ratio because this spot could not elastically disperse the stress. As injection continues, further escape of CO2 is mitigated by the heterogeneity of caprock. Compared to homogeneous model, heterogeneous model shows unstable increase of cumulative CO2 moles. Also, more heterogeneous model has larger difference compared with homogeneous model. Figure 6(c) shows the cumulative CO2 moles in petrophysically heterogeneous and geomechanically heterogeneous caprock. This model also shows more leakage in more heterogeneous model at early stage of injection, but after 4 years of injection, more heterogeneous model shows slightly more CO2 storage under the influence of geomechanical heterogeneity. After that, each model reveals similar leakage along continuous injection due to integrated effect of heterogeneous properties. Heterogeneous models show less difference compared to above two cases due to combined effect of petrophysical and geomechanical properties.Figure 5
Gas saturation of single realization forVDP = 0.7 and Case C.Figure 6
Cumulative in-place amounts of CO2 in aquifer at different times: Cases A, B, and C.
(a)
(b)
(c)In order to assess the range on uncertainty, a box-whisker plot format is used to illustrate minimum, 25th, 50th, 75th percentiles, and maximum values for sequestered CO2 after 10-year injection. As could be seen from Figure 7(a), the range of uncertainty is getting larger along the Dykstra–Parsons coefficient in Case A, and the amount of sequestered CO2 is getting smaller compared with homogeneous model. On the contrary to Figures 7(a) and 7(b) shows that geomechanical heterogeneity has favorable effect in CO2 trapping. Similar to Case A, the range of uncertainty is also increasing with Dykstra–Parsons coefficient getting bigger. Likewise, Case C shows the increase of uncertainty range as shown in Figure 7(c).Figure 7
Box-whisker plots of sequestered CO2 in aquifer after 10-year injection: Cases A, B, and C.
(a)
(b)
(c)Figure8, which shows the box-whisker plot comparing four cases, presents the effect of heterogeneity of each property after 10-year injection for most heterogeneous models (VDP = 0.7). Case A indicates 3.9% lower sequestered CO2. Heterogeneities of permeability and porosity have a negative influence on CO2 trapping because CO2 leakage occurs through the highly permeable shale first in petrophysically heterogeneous formation. However, in geomechanically heterogeneous model, leakage also occurs through the weak point in first 4 years. Then, geomechanical heterogeneity mitigates further escape of CO2 because geomechanically heterogeneous caprock includes elastic formation which has low Young’s modulus and Poisson’s ratio, and this section helps caprock endure the stress elastically. After 10 years of injection, Case B shows 4.6% higher amount of stored CO2. Due to combined effects of these two properties, Case C reveals only 0.3% difference compared with homogeneous case because composite effects of geomechanical and petrophysical properties are offset.Figure 8
Box-whisker plots of sequestered CO2 in aquifer after 10-year injection for all considered cases (VDP = 0.7).


### 4.2. Vertical Displacement

In Figure9, profiles of vertical displacements of caprock for three cases are compared. Continuous injection of CO2 increases fluid pressure in formation. Once the fluid pressure exceeds the least principal effective stress, the mechanical failure of caprock occurs, and CO2 leaks out of the storage formation. The vertical deformation of the caprock was calculated at the top of the caprock after 2 years of injection, just before the caprock cracks. Again, each curve indicates the average of 25 models. During CO2 injection into the formation, increased stress makes the surrounding formation rocks expand and push the caprock upwards. The caprock is bended more in the region just above the injector. As seen in Figure 9(a), for the case with heterogeneous petrophysical properties, all of the curves show higher vertical displacements than homogeneous case. Especially, values of vertical displacements are in order of VDP. This result suggests that top of the caprock of homogeneous model collapses earlier than other models because of pressure concentration. Similarly, Figure 9(b) shows higher vertical displacement of geomechanically heterogeneous model than homogeneous model regardless of VDP. Using petrophysically heterogeneous and geomechanically heterogeneous properties, Figure 9(c) depicts the combined effect of two properties.Figure 9
Vertical displacement of top of caprock: Cases A, B, and C.
(a)
(b)
(c)Figure10 reveals that the distinction of vertical displacement depends on VDP at the center of the caprock. Heterogeneous models at Figure 10(a) show 4.4%, 5.1%, and 6.3% higher vertical displacement compared with homogeneous case, respectively. Also, Case B shows 4.3%, 3.9%, and 4.1% higher displacement, and Case C shows 4.5%, 5.5%, and 6.2% higher displacement as could be seen in Figures 10(b) and 10(c). In comparison among the effect of heterogeneity for most heterogeneous models (VDP = 0.7), Cases A and C show similar range of uncertainty, and Case B shows less deformation but it is still larger than the homogeneous model. As could be seen in Figure 10(a), the vertical displacement of caprock increases along VDP, and the range of uncertainty is also increasing. Figure 10(b) also indicates that geomechanically heterogeneous caprock shows higher value of vertical displacement than homogeneous model, but uncertainty range and vertical displacement of heterogeneous caprock are relatively constant compared with Case A. It is well known that geomechanical properties have influence on rock deformation. Khajeh et al. [11] indicated that heterogeneous geomechanical properties result in significant differences in range of uncertainty in geomechanical response. In this study, geomechanical properties are correlated with petrophysical properties as a function. In this limited situation, heterogeneity in rock mechanical properties has less effect on the geomechanical change of formation in comparison to heterogeneity consideration of petrophysical properties.Figure 10
Box-whisker plots of vertical displacement: Cases A, B, and C.
(a)
(b)
(c)The amounts of sequestered CO2 and range of uncertainty in vertical displacement are larger for models with higher Dykstra–Parsons coefficient. It indicates that heterogeneity for petrophysical properties result in unfavorable effects in flow output variables. In contrast, the heterogeneity of rock mechanical properties can reduce additional escape of stored CO2. In terms of geomechanical change, heterogeneities of both petrophysical and geomechanical properties can affect the uplift of caprock. Highly heterogeneous caprock will be transformed more. Despite the heterogeneity of geomechanics affects the displacement of caprock, Case B shows relatively small range of uncertainty. It is obvious that geomechanical properties affect mechanical change of formation. However, there is limitation of displacement when reservoir properties are correlated. Since both geomechanical and flow responses are of importance to predict and optimize CO2 injection performance, ignoring heterogeneity effects for storage site properties may result in inaccurate analyses.


### 4.3. Effects of Horizontal Tectonic Stress

In order to assess the effect of tectonic stresses, horizontal stress is added to highly heterogeneous model (c) (VDP = 0.7) and homogeneous model along the x-axis. The simulation model has initial stress of 15,503.6 kPa with hydrostatic gradient. To add a 5,000 kPa at x-axis, horizontal compression is loaded to whole grids. Figure 11 indicates the sequestered CO2. Compared to cases considering normal stress only, the model with tectonic stress shows 24% more CO2 storage capacity just before the caprock failure during 4-5 years of injection. Due to the effect of additional horizontal compressive stress, the propagation of fractures is mitigated and a larger amount of CO2 is trapped. After the leakage occurs, the storage efficiency rapidly decreases for a homogeneous model with tectonic stress. Comparing two models including tectonic stress, the difference is up to 8% after 10 years of injection, and it is possible to get bigger. Vertical deformation at the center of the caprock is also increased by 0.6% at heterogeneous model and 1.1% at homogeneous model, respectively (Figure 12). Thus, sealing capacity could be better in more heterogeneous rocks under the high horizontal tectonic stress condition, but unpredicted high stress could make early fracture reactivation of caprock and bring unfavorable effect on long-term CO2 storage relatively homogeneous formation.Figure 11
Cumulative in-place amounts of CO2 in aquifer at different times under the additional tectonic stress (VDP = 0.7).Figure 12
Vertical displacement of top of caprock under the additional tectonic stress (VDP = 0.7).


## 4.1. CO2 Storage

Figure5 shows gas saturation of single realization and CO2 leakage regime. The figure is about Case C which has petrophysically heterogeneous properties and geomechanically heterogeneous properties and Dykstra–Parsons coefficient of 0.7. Increased formation pressure causes fractures in caprock, and CO2 escapes from the aquifer to overburden rock. Figure 6 depicts the cumulative amount of remaining CO2 in aquifer during 10-year injection. These graphs show the total amount of gas sequestered in the aquifer. Each curve is generated by averaging results from 25 heterogeneous models. The rest of gas has escaped to overburden formation. There could be further escape after injection, until the system reaches equilibrium. Figure 6(a) indicates simulation model with petrophysically heterogeneous and geomechanically homogeneous properties. By investigating the results of Case A, it can be seen that heterogeneous model shows less CO2 moles in the aquifer than homogeneous model. Also, the amount of leaked CO2 depends on heterogeneity. In other words, more petrophysically heterogeneous model presents more CO2 leakage. In heterogeneous caprock, there could be relatively weak shale which has high porosity and permeability. Injected fluids tend to flow through high permeable zone and partially increase formation pressure. Concentration of pressure incurs fracture reactivation and earlier leakage of CO2. Figure 6(b) shows simulation model with petrophysically homogeneous and geomechanically heterogeneous properties. In early stage of injection, heterogeneous models show more leakage than homogeneous model. However, as time goes on, the leakage is mitigated compared to homogeneous model. This difference is caused by imbalance of geomechanical properties. Being different from the heterogeneous caprock, the failure of homogeneous caprock occurs abruptly at the center. In heterogeneous caprock, the leakage does not occur above the injector. Initially, CO2 gas leaks through the location region with high Young’s modulus and Poisson’s ratio because this spot could not elastically disperse the stress. As injection continues, further escape of CO2 is mitigated by the heterogeneity of caprock. Compared to homogeneous model, heterogeneous model shows unstable increase of cumulative CO2 moles. Also, more heterogeneous model has larger difference compared with homogeneous model. Figure 6(c) shows the cumulative CO2 moles in petrophysically heterogeneous and geomechanically heterogeneous caprock. This model also shows more leakage in more heterogeneous model at early stage of injection, but after 4 years of injection, more heterogeneous model shows slightly more CO2 storage under the influence of geomechanical heterogeneity. After that, each model reveals similar leakage along continuous injection due to integrated effect of heterogeneous properties. Heterogeneous models show less difference compared to above two cases due to combined effect of petrophysical and geomechanical properties.Figure 5
Gas saturation of single realization forVDP = 0.7 and Case C.Figure 6
Cumulative in-place amounts of CO2 in aquifer at different times: Cases A, B, and C.
(a)
(b)
(c)In order to assess the range on uncertainty, a box-whisker plot format is used to illustrate minimum, 25th, 50th, 75th percentiles, and maximum values for sequestered CO2 after 10-year injection. As could be seen from Figure 7(a), the range of uncertainty is getting larger along the Dykstra–Parsons coefficient in Case A, and the amount of sequestered CO2 is getting smaller compared with homogeneous model. On the contrary to Figures 7(a) and 7(b) shows that geomechanical heterogeneity has favorable effect in CO2 trapping. Similar to Case A, the range of uncertainty is also increasing with Dykstra–Parsons coefficient getting bigger. Likewise, Case C shows the increase of uncertainty range as shown in Figure 7(c).Figure 7
Box-whisker plots of sequestered CO2 in aquifer after 10-year injection: Cases A, B, and C.
(a)
(b)
(c)Figure8, which shows the box-whisker plot comparing four cases, presents the effect of heterogeneity of each property after 10-year injection for most heterogeneous models (VDP = 0.7). Case A indicates 3.9% lower sequestered CO2. Heterogeneities of permeability and porosity have a negative influence on CO2 trapping because CO2 leakage occurs through the highly permeable shale first in petrophysically heterogeneous formation. However, in geomechanically heterogeneous model, leakage also occurs through the weak point in first 4 years. Then, geomechanical heterogeneity mitigates further escape of CO2 because geomechanically heterogeneous caprock includes elastic formation which has low Young’s modulus and Poisson’s ratio, and this section helps caprock endure the stress elastically. After 10 years of injection, Case B shows 4.6% higher amount of stored CO2. Due to combined effects of these two properties, Case C reveals only 0.3% difference compared with homogeneous case because composite effects of geomechanical and petrophysical properties are offset.Figure 8
Box-whisker plots of sequestered CO2 in aquifer after 10-year injection for all considered cases (VDP = 0.7).


## 4.2. Vertical Displacement

In Figure9, profiles of vertical displacements of caprock for three cases are compared. Continuous injection of CO2 increases fluid pressure in formation. Once the fluid pressure exceeds the least principal effective stress, the mechanical failure of caprock occurs, and CO2 leaks out of the storage formation. The vertical deformation of the caprock was calculated at the top of the caprock after 2 years of injection, just before the caprock cracks. Again, each curve indicates the average of 25 models. During CO2 injection into the formation, increased stress makes the surrounding formation rocks expand and push the caprock upwards. The caprock is bended more in the region just above the injector. As seen in Figure 9(a), for the case with heterogeneous petrophysical properties, all of the curves show higher vertical displacements than homogeneous case. Especially, values of vertical displacements are in order of VDP. This result suggests that top of the caprock of homogeneous model collapses earlier than other models because of pressure concentration. Similarly, Figure 9(b) shows higher vertical displacement of geomechanically heterogeneous model than homogeneous model regardless of VDP. Using petrophysically heterogeneous and geomechanically heterogeneous properties, Figure 9(c) depicts the combined effect of two properties.Figure 9
Vertical displacement of top of caprock: Cases A, B, and C.
(a)
(b)
(c)Figure10 reveals that the distinction of vertical displacement depends on VDP at the center of the caprock. Heterogeneous models at Figure 10(a) show 4.4%, 5.1%, and 6.3% higher vertical displacement compared with homogeneous case, respectively. Also, Case B shows 4.3%, 3.9%, and 4.1% higher displacement, and Case C shows 4.5%, 5.5%, and 6.2% higher displacement as could be seen in Figures 10(b) and 10(c). In comparison among the effect of heterogeneity for most heterogeneous models (VDP = 0.7), Cases A and C show similar range of uncertainty, and Case B shows less deformation but it is still larger than the homogeneous model. As could be seen in Figure 10(a), the vertical displacement of caprock increases along VDP, and the range of uncertainty is also increasing. Figure 10(b) also indicates that geomechanically heterogeneous caprock shows higher value of vertical displacement than homogeneous model, but uncertainty range and vertical displacement of heterogeneous caprock are relatively constant compared with Case A. It is well known that geomechanical properties have influence on rock deformation. Khajeh et al. [11] indicated that heterogeneous geomechanical properties result in significant differences in range of uncertainty in geomechanical response. In this study, geomechanical properties are correlated with petrophysical properties as a function. In this limited situation, heterogeneity in rock mechanical properties has less effect on the geomechanical change of formation in comparison to heterogeneity consideration of petrophysical properties.Figure 10
Box-whisker plots of vertical displacement: Cases A, B, and C.
(a)
(b)
(c)The amounts of sequestered CO2 and range of uncertainty in vertical displacement are larger for models with higher Dykstra–Parsons coefficient. It indicates that heterogeneity for petrophysical properties result in unfavorable effects in flow output variables. In contrast, the heterogeneity of rock mechanical properties can reduce additional escape of stored CO2. In terms of geomechanical change, heterogeneities of both petrophysical and geomechanical properties can affect the uplift of caprock. Highly heterogeneous caprock will be transformed more. Despite the heterogeneity of geomechanics affects the displacement of caprock, Case B shows relatively small range of uncertainty. It is obvious that geomechanical properties affect mechanical change of formation. However, there is limitation of displacement when reservoir properties are correlated. Since both geomechanical and flow responses are of importance to predict and optimize CO2 injection performance, ignoring heterogeneity effects for storage site properties may result in inaccurate analyses.


## 4.3. Effects of Horizontal Tectonic Stress

In order to assess the effect of tectonic stresses, horizontal stress is added to highly heterogeneous model (c) (VDP = 0.7) and homogeneous model along the x-axis. The simulation model has initial stress of 15,503.6 kPa with hydrostatic gradient. To add a 5,000 kPa at x-axis, horizontal compression is loaded to whole grids. Figure 11 indicates the sequestered CO2. Compared to cases considering normal stress only, the model with tectonic stress shows 24% more CO2 storage capacity just before the caprock failure during 4-5 years of injection. Due to the effect of additional horizontal compressive stress, the propagation of fractures is mitigated and a larger amount of CO2 is trapped. After the leakage occurs, the storage efficiency rapidly decreases for a homogeneous model with tectonic stress. Comparing two models including tectonic stress, the difference is up to 8% after 10 years of injection, and it is possible to get bigger. Vertical deformation at the center of the caprock is also increased by 0.6% at heterogeneous model and 1.1% at homogeneous model, respectively (Figure 12). Thus, sealing capacity could be better in more heterogeneous rocks under the high horizontal tectonic stress condition, but unpredicted high stress could make early fracture reactivation of caprock and bring unfavorable effect on long-term CO2 storage relatively homogeneous formation.Figure 11
Cumulative in-place amounts of CO2 in aquifer at different times under the additional tectonic stress (VDP = 0.7).Figure 12
Vertical displacement of top of caprock under the additional tectonic stress (VDP = 0.7).


## 5. Conclusion

Comparisons among storage sites with heterogeneous caprock showed the significance of coupled simulation of fluid flow and geomechanical deformation with respect to CO2 storage integrity. In particular, they have shown the importance of correlation between hydraulic and geomechanical properties in controlling CO2 leakage through caprock. Four major factors affecting mechanical caprock failure considered in this study are permeability, porosity, Young’s modulus, and Poisson’s ratio. These four properties are correlated as linear relationship based on field data. Heterogeneous petrophysical properties have unfavorable effects on CO2 storage, and geomechanical heterogeneity could mitigate further escape of CO2 after leakage occurs. Heterogeneity effects of these two properties can be offset for CO2 trapping. Heterogeneity of these properties also affects on vertical displacement of caprock. Results of analysis show that heterogeneity of four properties induces further displacement of caprock, but geomechanical properties have limited effects when it is correlated with other properties as a function. Subsurface storage needs to be evaluated in terms of heterogeneity to avoid any risk when operating CO2 sequestration. The hydrogeomechanical correlations are applied for more accurate investigation of CO2 leakage. We considered only four properties because correlations between hydraulic properties and other geomechanical properties have not been thoroughly understood. However, investigation including other formation properties such as rock compressibility and yield stress could provide better understanding of geomechanical study. Another improvement on the analysis of caprock stability can be achieved by the evaluation of 3D model. It could be an effective way to evaluate more realistic CO2 flow regime.



---

*Source: 1474320-2018-02-28.xml*