# CO2 Permeability Analysis of Caprock Containing a Single Fracture Subject to Coupled Thermal-Hydromechanical Effects

Authors: Qian Yin; Hongwen Jing; Haijian Su; Huidong Wang
Journal: Mathematical Problems in Engineering
 (2017)

---

## Abstract

Coupled THM (thermal-hydromechanical) processes have become increasingly important in studying the issues affecting subsurface flow systems. CO2 permeability of the fracture in caprock is a key factor that affects sealing efficiency of caprock. A new model associated with coupled THM processes that shows a good reliability was derived. Then, based on the COMSOL multiphysics software, a series of numerical calculations were performed on caprock models with a single fracture subject to coupled THM effects. Transmissivity of the fracture as a function of fracture angle, overburden pressure, fluid pressure difference, injected CO2 temperature, and the initial fracture aperture was elucidated, respectively. Average transmissivity of the fracture undergoes an increase by 1.74 times with the fracture angle (45°–90°), 2-3 orders of magnitude with the fluid pressure difference (5–30 MPa), and 4-5 orders of magnitude with the initial fracture aperture (0.05–0.5 mm), while it decreases by 3-4 orders of magnitude as overburden pressure increases from 30 to 80 MPa. Injected CO2 temperature has a small impact on the fracture permeability. This work provides an alternative tool to enrich the numerical modeling for the assessment of CO2 caprock sealing efficiency.


---

## Body

## 1. Introduction

To address the increasing concerns regarding carbon dioxide emission and its impact on climate change, CO2 geological sequestration has become a promising approach [1–5]. To reduce the atmosphere emissions, a large amount of CO2 is injected into the geological storage reservoirs, as shown in Figure 1, which may be gradually accumulated at the bottom of caprocks and lead to stress field changes in caprock. However, if the reservoir pressure is high enough to cause mechanical failure in caprock and connected pathways are created through fractures, a potential hazard of CO2 leakage will occur [6–8].Figure 1
Schematic of CO2 injection in the presence of fractures within the caprock layer.Research on subsurface CO2 flow systems involves thermal (T), hydrodynamic (H), and mechanical (M) processes. In fact, these processes are interrelated and affect each other and are referred to as “coupled THM effects” [9], which has a significant influence on sealing efficiency of caprock [10]. Permeability of the fracture in caprock is the key safety issue for CO2 geological sequestration in storage reservoirs.Numerical simulations have been widely used to evaluate CO2 multiphase flow, diffusion, geomechanics, and chemical reactions during CO2 injection and storage. In the multiphase flow research field, TOUGH2 codes, which consider the cross-coupling of TH and THC processes for multiphase flow, were developed [11]. Two existing well-established codes, TOUGH2 and FLAC3D, have been adopted as a pragmatic approach for modeling coupled THM processes [12]. Besides, TOUGHREACT and FLAC3D have been linked together to simulate coupled THMC processes [13]. A novel fully coupled flow and geomechanics model TOUGH2-EGS in enhanced geothermal reservoirs based on average Navier equation was developed [14]. FEHM finite element codes were also applied to simulate coupled THM processes in subsurface fractured media [15]. A mechanical simulation module TOUGH2Biot, which was based on the extended Biot consolidation model and finite element method, was developed and applied to CO2 sequestration simulation [9]. COMSOL multiphysics software, which can be used to simulate ground water flow subject to coupled THM effects if a suitable template and relationship between the coupled processes are constructed, has been widely employed recently [16]. In previous studies, even though substantial efforts have been devoted to estimation and prediction of the CO2 sequestration performance, the theory for fracture permeability in caprock and the corresponding sealing efficiency of caprock have so far not been fully developed due to dual complexity of THM coupled processes and geological conditions.This paper is organized as follows. A new coupled thermal-hydromechanical model for CO2 flow through a single fracture in caprock was first derived. The governing equations were linked with COMSOL multiphysics software, and the reliability of the model was verified using a sample problem. Finally, several numerical calculations on caprock models with a single fracture subject to coupled THM effects were performed, and CO2 permeability of the fracture with respect to different fracture angle, overburden pressure, fluid pressure difference, injected CO2 temperature, and initial aperture was, respectively, evaluated. In this study, these models were calculated under simplified conditions of single-phase flow and heat conduction alone in thermal field for brevity.


## 2. Governing Equations

In the following, a set of field equations are defined which govern the deformation of caprock matrix, the fluid flow through the fracture, and the heat conduction process. These derivations are obtained based on the following assumptions. (i) Caprock matrix is a kind of homogeneous, isotropic, and elastic continuum. (ii) Strains are much smaller than the length scale. (iii) No crack propagation happens to the caprock matrix and no dislocation occurs between the matrix blocks. (iv) The matrix is impermeable, and CO2 flows through the fracture alone. The fluid flow behavior can be described using Darcy’s law. (v) Heat effect induced by fluid flow through the fracture is negligible, and heat conduction within the matrix follows Fourier’s law. (vi) Density and viscosity of the supercritical CO2 vary with the temperature and pressure.


### 2.1. Governing Equation for Caprock Matrix Deformation

To elucidate the mechanical response of caprock containing fractures under coupled THM effects and to evaluate the permeability of fractures within the caprock, a typical mechanical model is built and shown in Figure2(a). The lower boundary of the model is displacement constraint. Due to continuous CO2 injection, pressure is built up at the lower boundary of the caprock. Buried depth of caprock provides a vertical pressure of q at the upper boundary, and both lateral sides are subjected to an equal pressure of ql. Then, stress analysis of a microunit chosen from the caprock matrix (the dashed region in Figure 2(a)) is conducted, as shown in Figure 2(b). T0 and TC, respectively, denote the temperature of caprock matrix and injected CO2. For a homogeneous, isotropic, and elastic medium, the strain-displacement relation of the matrix can be expressed as (1)εijT+εijq=12ui,jT+ui,jq+uj,iT+uj,iq,where εT and εq are component of the total strain tensor produced by temperature and applied load, respectively. uT and uq are the component of displacement. The equilibrium equation is defined as(2)σij,iT+σij,iq+fj=0,where σT and σq denote the component of the total stress tensor produced by the temperature and applied load. fj denotes the component of the body force.Figure 2
Mechanical model of caprock containing fractures.
(a)
(b)Then, the constitutive relation for the deformed caprock matrix becomes(3)σij=2Gεij+λεkkδij-αEΔT1-2νδij,where λ and G are Lame constants, ν is Poisson’s ratio, E is Young’s modulus, ΔT is the temperature variation, α is the thermal expansion coefficient of rock, and δij is the Kronecker delta. Combination of (1)–(3) yields the equilibrium differential equation, written as(4)G∇2ui+λ+Guj,jj-αE1-2νΔT,i+fi=0.Equation (4) is the governing equation for caprock deformation, from which stress distribution of caprock can be calculated.Natural fractures are often subjected to field stresses or mechanical displacements, which have a direct influence on the fracture aperture and hence the permeability of the fractured rock. The fracture aperture may increase due to shear dilation [20–23] or decrease in response to the normal loads [24–27]. At present, relations between fracture permeability and the applied normal loads are clear and have been widely accepted by scholars in the field. However, no unified understanding has been obtained for effect of shear stress on the fracture permeability. Therefore, in this study, the role of shear stress on CO2 flows through the fracture is not taken into account.Figure3 shows a typical mechanical model of a fracture. When CO2 flows in the fracture, a fluid pressure of Pf is applied on the fracture surface. σfn and kn, respectively, denote the normal contact stress and normal stiffness between two blocks.Figure 3
Mechanical model of the fracture.According to the definition of effective stress in porous media, the effective stress in the fracture can be written as(5)σfne=σfn-Pf.By using the distinct element code of UDEC, a simple description of the relation betweenσfne and the mechanical aperture bm was given [28], as indicated by the solid lines in Figure 4. In the range of residual aperture bres and maximum aperture bmax, σfne is linearly proportional to bm, written as(6)bm=bm0-σfnekn,where bm0 is the initial fracture aperture with no applied load.Figure 4
Fracture aperture as a function of the effective normal stress.From Figure4, when σfne>0, fracture closure happens, and with the increase of σfne, bm continues to decrease linearly until bres. When bm=bres, the fracture aperture keeps a constant value which no longer depends on σfne. When σfne<0, bm increases with the increase of -σfne until bmax. When bm=bmax, bm remains unchanged.In this study, forσfne<0, the linear relation between σfne and bm in Figure 4 is still adopted to describe the variation of bm. However, for σfne>0, the Barton-Bandis equation is utilized to evaluate the fracture closure behavior, as shown by the dashed line in Figure 4. The equation is a kind of hyperbolic model put forward by Bandis et al. [29] through experiments to describe the fracture deformation with the effective normal stress, written as(7)ΔVjc=σfnekn0+σfne/bm0,where ΔVjc is the fracture closure and kn0 is the initial normal stiffness with no applied load.According to (7), the normal stiffness of the fracture can be determined as(8)kn=∂σfne∂Vjc=kn01-σfnekn0bm0+σfne-2and the fracture aperture can be described as(9)bm=bm0-ΔVjc.


### 2.2. Governing Equation for Fluid Flow

For the two dimensional model described in Figure2, fluid flow through the fracture can be solved as a one-dimensional problem. Figure 5 presents the hydraulic calculation model for a microunit in the fracture i.Figure 5
Hydraulic model of fluid flow through a fracture.In terms of the local coordinate system of a fracture, for continuously saturated fluid flow, the mass conservation equation regardless of the source sink can be written as(10)-∂ρCvibh∂xidxi=∂ρCbhdxi∂t,where ρC is the density of CO2, vi denotes the flow velocity, t is the time, and bh is the equivalent hydraulic aperture.According to Darcy’s law,vi is given by(11)vi=-Ki∂hi∂xi,where Ki is the permeability of the fracture.By neglecting the velocity head, the relationship between total headhi and the osmotic pressure Pi can be written as(12)Pi=hi-ziρCg,where hi is the head distribution of the fracture and zi is the position head of the fracture corresponding to the global coordinate system.Taking a derivative of (12) yields the following equation:(13)∂hi∂xi=1ρCg∂Pi∂xi+∂zi∂xi-PiρC2g∂ρC∂xi.Substituting (13) into (11) and then into (10), we obtain(14)∂∂xiTfig∂Pi∂xi+ρCg∂zi∂xi-PiρC∂ρC∂xidxi=∂ρCbhdxi∂t,where Tfi is the transmissivity of fracture and xi is the local coordinate of fracture.During the fracture deformation process (closure/opening), the length ofdxi is unchanged, while the fluid density and fracture aperture vary. Equation (14) can be rewritten as(15)∂ρCbhdxi∂t=∂ρC∂tbhdxi+∂bh∂tρCdxi.Assuming that the variation of fluid concentration is negligible, compression coefficient and the expansion coefficient of the fluid can be, respectively, written as [30](16)βP=1ρC∂ρCP∂P,Tis constantβT=-1ρC∂ρCT∂T,Pis constant,where ρCP and ρCT, respectively, denote the variation of fluid density induced by external pressure P and temperature T. Therefore, the total variation of density can be described as(17)∂ρC∂t=∂ρCP∂t+∂ρCT∂t=ρCβP∂P∂t-βT∂T∂t.The compression coefficient of the fracture can be written as(18)α=-1bhdxi∂bhdxi∂σfne.Since the fracture length does not vary withσfne, (18) can be rewritten as(19)α=-1bh∂bh∂σfne=δnbh,where δn is the normal flexibility coefficient of the fracture, which can be written as(20)δn=1kn=∂ΔVjc∂σfne.When the external load is ensured, the total stress of the fracture will keep constant. The relationship between the effective normal stress and fluid pressure can be written as(21)dσfne=-dPi.Substituting (21) into (19), we obtain(22)∂bh∂t=δn∂Pi∂t.Combination of (14), (15), (17), and (22) yields the following governing equation for fluid flow, expressed as(23)ρCgbhβP+δn∂P∂t+∂∂xiTf∂Pi∂xi+ρCg∂zi∂xi-PiρC∂ρC∂xi=ρCgbhβT∂T∂t.Assuming thatbh is equal to bm and fluid flow through the fracture can be described using Darcy’s law, permeability of the fracture can be written as (24)Kfi=ρCgbm02kn0+σfnebm0-1/bm0kn0+σfne212μC.The transmissivity can be written as(25)Tfi=Kfibhi=ρCgbm02kn0+σfnebm0-1/bm0kn0+σfne312μC,where μC is the dynamic viscosity of CO2 and g presents the gravitational acceleration.


### 2.3. Governing Equation for Heat Conduction

The injection of CO2 could result in variation of the temperature in caprock, which would then produce thermal stress [31]. Besides, the temperature alterations also lead to change of the physical properties of CO2. Therefore, the temperature has a significant influence on permeability of the fracture. In this study, thermal effects produced by CO2 flow through the fracture are neglected. Instead, the temperature of the injected CO2 is regarded as a known condition which is imposed on the lower boundary of caprock. When the temperature of CO2 is different from that of caprock, heat conduction phenomenon occurs. The governing equation can be written as(26)∂T∂t=λsρscs∂2T∂x2+∂2T∂y2,where T denotes temperature variables and λs, cs, and ρs, respectively, denote the heat conduction coefficient, heat capacity, and density of rock matrix.


### 2.4. Property Parameters of Supercritical CO2

Property parameters of the supercritical CO2 injected in the storage layer vary with the pressure and temperature. In this study, the density and dynamic viscosity of the supercritical CO2 are, respectively, determined according to the empirical models put forward by Span and Wagner [17] and Vesovic et al. [18], as shown in Figure 6.Figure 6
Variations of (a) density [17] and (b) dynamic viscosity [18] of supercritical CO2 with the temperature and pressure.
(a)
ρ
C (T: 300–400 K, P: 5–50 MPa)
(b)
μ
C (T: 260–400 K, P: 5–50 MPa)


## 2.1. Governing Equation for Caprock Matrix Deformation

To elucidate the mechanical response of caprock containing fractures under coupled THM effects and to evaluate the permeability of fractures within the caprock, a typical mechanical model is built and shown in Figure2(a). The lower boundary of the model is displacement constraint. Due to continuous CO2 injection, pressure is built up at the lower boundary of the caprock. Buried depth of caprock provides a vertical pressure of q at the upper boundary, and both lateral sides are subjected to an equal pressure of ql. Then, stress analysis of a microunit chosen from the caprock matrix (the dashed region in Figure 2(a)) is conducted, as shown in Figure 2(b). T0 and TC, respectively, denote the temperature of caprock matrix and injected CO2. For a homogeneous, isotropic, and elastic medium, the strain-displacement relation of the matrix can be expressed as (1)εijT+εijq=12ui,jT+ui,jq+uj,iT+uj,iq,where εT and εq are component of the total strain tensor produced by temperature and applied load, respectively. uT and uq are the component of displacement. The equilibrium equation is defined as(2)σij,iT+σij,iq+fj=0,where σT and σq denote the component of the total stress tensor produced by the temperature and applied load. fj denotes the component of the body force.Figure 2
Mechanical model of caprock containing fractures.
(a)
(b)Then, the constitutive relation for the deformed caprock matrix becomes(3)σij=2Gεij+λεkkδij-αEΔT1-2νδij,where λ and G are Lame constants, ν is Poisson’s ratio, E is Young’s modulus, ΔT is the temperature variation, α is the thermal expansion coefficient of rock, and δij is the Kronecker delta. Combination of (1)–(3) yields the equilibrium differential equation, written as(4)G∇2ui+λ+Guj,jj-αE1-2νΔT,i+fi=0.Equation (4) is the governing equation for caprock deformation, from which stress distribution of caprock can be calculated.Natural fractures are often subjected to field stresses or mechanical displacements, which have a direct influence on the fracture aperture and hence the permeability of the fractured rock. The fracture aperture may increase due to shear dilation [20–23] or decrease in response to the normal loads [24–27]. At present, relations between fracture permeability and the applied normal loads are clear and have been widely accepted by scholars in the field. However, no unified understanding has been obtained for effect of shear stress on the fracture permeability. Therefore, in this study, the role of shear stress on CO2 flows through the fracture is not taken into account.Figure3 shows a typical mechanical model of a fracture. When CO2 flows in the fracture, a fluid pressure of Pf is applied on the fracture surface. σfn and kn, respectively, denote the normal contact stress and normal stiffness between two blocks.Figure 3
Mechanical model of the fracture.According to the definition of effective stress in porous media, the effective stress in the fracture can be written as(5)σfne=σfn-Pf.By using the distinct element code of UDEC, a simple description of the relation betweenσfne and the mechanical aperture bm was given [28], as indicated by the solid lines in Figure 4. In the range of residual aperture bres and maximum aperture bmax, σfne is linearly proportional to bm, written as(6)bm=bm0-σfnekn,where bm0 is the initial fracture aperture with no applied load.Figure 4
Fracture aperture as a function of the effective normal stress.From Figure4, when σfne>0, fracture closure happens, and with the increase of σfne, bm continues to decrease linearly until bres. When bm=bres, the fracture aperture keeps a constant value which no longer depends on σfne. When σfne<0, bm increases with the increase of -σfne until bmax. When bm=bmax, bm remains unchanged.In this study, forσfne<0, the linear relation between σfne and bm in Figure 4 is still adopted to describe the variation of bm. However, for σfne>0, the Barton-Bandis equation is utilized to evaluate the fracture closure behavior, as shown by the dashed line in Figure 4. The equation is a kind of hyperbolic model put forward by Bandis et al. [29] through experiments to describe the fracture deformation with the effective normal stress, written as(7)ΔVjc=σfnekn0+σfne/bm0,where ΔVjc is the fracture closure and kn0 is the initial normal stiffness with no applied load.According to (7), the normal stiffness of the fracture can be determined as(8)kn=∂σfne∂Vjc=kn01-σfnekn0bm0+σfne-2and the fracture aperture can be described as(9)bm=bm0-ΔVjc.


## 2.2. Governing Equation for Fluid Flow

For the two dimensional model described in Figure2, fluid flow through the fracture can be solved as a one-dimensional problem. Figure 5 presents the hydraulic calculation model for a microunit in the fracture i.Figure 5
Hydraulic model of fluid flow through a fracture.In terms of the local coordinate system of a fracture, for continuously saturated fluid flow, the mass conservation equation regardless of the source sink can be written as(10)-∂ρCvibh∂xidxi=∂ρCbhdxi∂t,where ρC is the density of CO2, vi denotes the flow velocity, t is the time, and bh is the equivalent hydraulic aperture.According to Darcy’s law,vi is given by(11)vi=-Ki∂hi∂xi,where Ki is the permeability of the fracture.By neglecting the velocity head, the relationship between total headhi and the osmotic pressure Pi can be written as(12)Pi=hi-ziρCg,where hi is the head distribution of the fracture and zi is the position head of the fracture corresponding to the global coordinate system.Taking a derivative of (12) yields the following equation:(13)∂hi∂xi=1ρCg∂Pi∂xi+∂zi∂xi-PiρC2g∂ρC∂xi.Substituting (13) into (11) and then into (10), we obtain(14)∂∂xiTfig∂Pi∂xi+ρCg∂zi∂xi-PiρC∂ρC∂xidxi=∂ρCbhdxi∂t,where Tfi is the transmissivity of fracture and xi is the local coordinate of fracture.During the fracture deformation process (closure/opening), the length ofdxi is unchanged, while the fluid density and fracture aperture vary. Equation (14) can be rewritten as(15)∂ρCbhdxi∂t=∂ρC∂tbhdxi+∂bh∂tρCdxi.Assuming that the variation of fluid concentration is negligible, compression coefficient and the expansion coefficient of the fluid can be, respectively, written as [30](16)βP=1ρC∂ρCP∂P,Tis constantβT=-1ρC∂ρCT∂T,Pis constant,where ρCP and ρCT, respectively, denote the variation of fluid density induced by external pressure P and temperature T. Therefore, the total variation of density can be described as(17)∂ρC∂t=∂ρCP∂t+∂ρCT∂t=ρCβP∂P∂t-βT∂T∂t.The compression coefficient of the fracture can be written as(18)α=-1bhdxi∂bhdxi∂σfne.Since the fracture length does not vary withσfne, (18) can be rewritten as(19)α=-1bh∂bh∂σfne=δnbh,where δn is the normal flexibility coefficient of the fracture, which can be written as(20)δn=1kn=∂ΔVjc∂σfne.When the external load is ensured, the total stress of the fracture will keep constant. The relationship between the effective normal stress and fluid pressure can be written as(21)dσfne=-dPi.Substituting (21) into (19), we obtain(22)∂bh∂t=δn∂Pi∂t.Combination of (14), (15), (17), and (22) yields the following governing equation for fluid flow, expressed as(23)ρCgbhβP+δn∂P∂t+∂∂xiTf∂Pi∂xi+ρCg∂zi∂xi-PiρC∂ρC∂xi=ρCgbhβT∂T∂t.Assuming thatbh is equal to bm and fluid flow through the fracture can be described using Darcy’s law, permeability of the fracture can be written as (24)Kfi=ρCgbm02kn0+σfnebm0-1/bm0kn0+σfne212μC.The transmissivity can be written as(25)Tfi=Kfibhi=ρCgbm02kn0+σfnebm0-1/bm0kn0+σfne312μC,where μC is the dynamic viscosity of CO2 and g presents the gravitational acceleration.


## 2.3. Governing Equation for Heat Conduction

The injection of CO2 could result in variation of the temperature in caprock, which would then produce thermal stress [31]. Besides, the temperature alterations also lead to change of the physical properties of CO2. Therefore, the temperature has a significant influence on permeability of the fracture. In this study, thermal effects produced by CO2 flow through the fracture are neglected. Instead, the temperature of the injected CO2 is regarded as a known condition which is imposed on the lower boundary of caprock. When the temperature of CO2 is different from that of caprock, heat conduction phenomenon occurs. The governing equation can be written as(26)∂T∂t=λsρscs∂2T∂x2+∂2T∂y2,where T denotes temperature variables and λs, cs, and ρs, respectively, denote the heat conduction coefficient, heat capacity, and density of rock matrix.


## 2.4. Property Parameters of Supercritical CO2

Property parameters of the supercritical CO2 injected in the storage layer vary with the pressure and temperature. In this study, the density and dynamic viscosity of the supercritical CO2 are, respectively, determined according to the empirical models put forward by Span and Wagner [17] and Vesovic et al. [18], as shown in Figure 6.Figure 6
Variations of (a) density [17] and (b) dynamic viscosity [18] of supercritical CO2 with the temperature and pressure.
(a)
ρ
C (T: 300–400 K, P: 5–50 MPa)
(b)
μ
C (T: 260–400 K, P: 5–50 MPa)


## 3. Coupled Model Validation

A 2D single fracture model, which has been studied before by Bower and Zyvoloski [19], is built and calculated to verify the effectiveness of the coupled model discussed above. For the coupled THM model in this study, the thermal field just plays a role in thermal stress within the rock matrix and property parameters of the fluid, with a clear physical process. Thus, we just make a comparison of the hydromechanical calculation results with those reported by Bower and Zyvoloski [19]. Numerical model is shown in Figure 7, with the input parameters listed in Table 1.Table 1
Input parameters of the validation model [19].
Parameter
Value
Note
Length, m
25
Model
Width, m
1
Density, kg/m3
2716
Matrix
Young’s modulus, MPa
1000
Poisson’s ratio
0.0
Normal stiffness, MPa/m
1 × 106
Fracture
Porosity
1.0
Initial aperture, m
1 × 10−5
Residual aperture, m
1 × 10−30
Maximum aperture, m
0.002
Density, kg/m3
1000
Fluid
Dynamic viscosity, Pa⋅s
0.001
Compression coefficient, 1/Pa
0.0Figure 7
Validation model [19].An initial stress field ofPi = 21.0 MPa is imposed in both the matrix and fracture, and the initial fracture aperture is set to be 1 × 10−5 m. The fluid is continuously injected in the fracture from the left side with Pl = 21.9 MPa, and Pr at right side is kept constant of 21.0 MPa. The left, upper and lower model boundaries are all displacement constraint. Normal stiffness of the fracture is kept unchanged of 1 × 106 MPa/m.For this case, numerical and analytic solutions of aperture variation along the fracture length att = 500 and 2000 days have been given by Bower and Zyvoloski [19]. The numerical results were calculated using the FEHM codes, and the analytic solutions were obtained using the method put forward by Wijesinghe [32]. Then, the fracture aperture was recalculated by solving the hydromechanical coupled model derived in this study with the COMSOL multiphysics. The comparison results are shown in Figure 8. To evaluate how these curves fit well with each other, an evaluation coefficient μ is put forward [33]:(27)μ=∑i=1nΔbhc-Δbh22n,where n is the total number of measuring points, Δbhc is the fracture aperture change calculated in our study, and Δbh2 refers to fracture aperture change calculated using other methods.Figure 8
Comparison of analytic, numerical (FEHM) aperture [19] with the calculated results using COMSOL multiphysics in our study.By using (27), μ between analytic solution and the calculated results in our study is only 3.36 × 10−7 and 1.63 × 10−7 m, respectively, at t = 500 and 2000 days. Besides, μ between the FEHM results and our results is only 3.79 × 10−7 and 3.18 × 10−7 m at t = 500 and 2000 days, respectively. Obviously, the calculation results agree well with the numerical and analytic results of Bower and Zyvoloski [19], which indicates a good reliability of the coupled model in our study.


## 4. CO2 Permeability Analysis of Single Fracture in Caprock

### 4.1. Model Setup

To quantitatively analyze permeability of the single fracture in caprock, a conceptual model is set up, as shown in Figure9. The model is square with the size of 10 m × 10 m. The fracture connects the upper and lower model boundary and passes through the center position of the model. The right, left, and lower model boundary are displacement constraint. A vertical load of q is applied on the upper model boundary. Fluid pressure of Pb and Pt is, respectively, applied at the lower and upper side of the fracture. Temperature of Tb and Tt is imposed on the lower and upper model boundary. Input parameters are listed in Table 2.Table 2
Parameters used in 2D THM coupled model.
Parameter
Value
Note
Density, kg/m3
2500
Rock matrix
Young’s modulus, GPa
100
Poisson’s ratio
0.3
Heat conductivity, W/(m × K)
2.57
Specific heat capacity, J/(m × K)
710
Thermal expansion, 1/K
6.0 × 10−7
Initial normal stiffness of fracture closure, GPa/m
60
Fracture
Normal stiffness of fracture opening, GPa/m
100
Initial aperture, m
b
m
0
Residual aperture, m
1 × 10−6
Maximum aperture, m
0.002
Fluid pressure at the lower side, MPa
P
b
Fluid pressure at the upper side, MPa
P
t
Density, kg/m3
Figure6(a)
Supercritical CO2
Dynamic viscosity, Pa⋅s
Figure6(b)
Compression coefficient, 1/Pa
Calculated with Figure6(a)
Thermal expansion, 1/K
Calculated with Figure6(b)Figure 9
2D conceptual model of a single fracture in caprock subjected to coupled THM effects, in whichβ denotes included angle between the fracture and horizontal direction, T0 is the initial temperature of rock matrix, and Tb is equal to the temperature of injected CO2.


### 4.2. Results and Discussion

First, variation ofσfne, Pf, and Tf of the fracture, as well as ρC and μC of the supercritical CO2 under THM coupled effects at different times for a vertical fracture (β = 90°) were analyzed. Then, permeability of the fracture in caprock with respect to different q, Pb, bm0, β, and Tb was, respectively, studied.


#### 4.2.1. THM Coupled Effects on a Vertical Fracture

Boundary and initial conditions are as follows:β = 90°, bm0 = 0.5 mm, q = 60 MPa, Pb = 30 MPa, Pt = 10 MPa, Tb = 333.15 K, Tt = 303.15 K, and T0 = 303.15 K.At initial time, CO2 was continuously injected in the fracture through the lower fracture side with a constant Pb = 30 MPa. At different times, temperature distribution along the whole fracture length is shown in Figure 10, in which, x-axis (L0) represents the distance from the measure point to the lower fracture tip along the fracture direction (in the range of 0–10 m). During the fluid flow test, heat conduction occurs in rock due to higher temperature of injected CO2 (Tb) than that of the rock matrix (T0). The heat transfers from high-temperature to the low-temperature position and finally approaches a stable state. From Figure 10, with the increase of time, heat at the lower fracture segment gradually transfers to the upper position until a stable state is achieved at t = 300 days.Figure 10
Temperature distribution along the fracture at different times.During the heat conduction process, thermal expansion happens to the rock. Since both lateral boundaries (left and right) of the model are displacement constraint, the thermal expansion of rock leads to variation of normal stress in the fracture. Figures11(a)–11(c), respectively, show σfn, Pf, and σfne along the fracture length at different times. By comparing Figures 10 and 11(a), generally, a high temperature in fracture corresponds to a high total normal stress, and a steady stress field is achieved when the thermal field is steady. Variation of σfn in the fracture with time results mainly from the thermal field.Figure 11
Distribution of (a) total normal stress, (b) fluid pressure, and (c) effective normal stress along the fracture at different times.
(a)
σ
f
n
(b)
P
f
(c)
σ
f
n
eFluid pressurePf in the fracture as a function of L0 at different times is shown in Figure 11(b). With the increase of time, CO2 flows along the direction of pressure drops, and, for this case, Pf reaches steady values at t = 30 days. By using (5), variation of σfne of the fracture at different times can be obtained (Figure 11(c)). In a certain small range of L0, Pf⁡ is larger than σfn, which results in negative values of σfne and a larger bh than bm0, as shown by the dashed line in Figures 11(c) and 12. Besides, σfne increases gradually as L0 increases, which leads to corresponding decrease of bh. It can also be seen from Figure 12 that, in the range of L0 from 4 to 10 m, bh first increases and then decreases, finally reaching stable values, as indicated by the red arrows. The main reason for this phenomenon is due to earlier stable state of hydrodynamic field (t = 30 days) than that of the thermal field (t = 300 days).Figure 12
Variation ofbh along the fracture at different times.Distribution ofρC and μC of the supercritical CO2 along the fracture at different times is displayed in Figure 13. For L0 = 3–10 m, both property parameters show an ascending-descending variation before attaining constant values, which are marked by the red arrows.Figure 13
Variation of (a) density and (b) viscosity of supercritical CO2 along the fracture at different times.
(a)
Density
(b)
ViscosityUnder the coupled THM effects and taking variation of property parameters of CO2 with different temperature and pressure into account, the transmissivity Tf of the fracture can be evaluated (Figure 14). Once hydrodynamic field of the fracture is stable, Tf will reach constant values.Figure 14
Transmissivity distribution along the fracture at different times.


#### 4.2.2. Effect of Included Angleβ

Due to complicated geological conditions of the caprock, included angle between fracture and the horizontal direction is various. Thus, it is of great significance to study the effect ofβ on permeability of the fracture under coupled THM effects. Six models with various β (45°, 51°, 59°, 68°, 79°, and 90°) were, respectively, set up. Boundary and initial conditions are as follows: bm0 = 0.5 mm, q = 60 MPa, Pb = 20 MPa, Pt = 10 MPa, Tb = 333.15 K, Tt = 303.15 K, and T0 = 303.15 K. Other input parameters were the same as those listed in Table 2. Figure 15(a) presents variation of average Tf of the fracture with different β. With the increase of time, the average Tf increases gradually and then attains a constant value. Under coupled THM effects, the average Tf experiences an increasing trend with the increase of β. The average Tf of the fracture with included angle β of 45°, 51°, 59°, 68°, 79°, and 90° is 1.66 × 10−5, 2.01 × 10−5, 2.55 × 10−5, 3.32 × 10−5, 4.15 × 10−5, and 4.55 × 10−5 m2/s. The average Tf for β = 90° increases by 1.74 times over that for β = 45°, which indicates a weaker sequestration performance of caprock with a larger β.Figure 15
Effect of (a) included angleβ, (b) overburden pressure q, (c) fluid pressure difference ΔP, (d) temperature Tb, and (e) initial aperture bm0 on average Tf of the fracture subject to coupled THM effects.
(a)
Included angleβ
(b)
Overburden pressureq
(c)
Fluid pressure differenceΔP
(d)
TemperatureTb
(e)
Initial aperturebm0


#### 4.2.3. Effect of Overburden Pressureq

When the storage layer of CO2 is in different buried depths, the overburden pressure q of caprock is various, which would then largely impact the permeability of fracture. Six models with different q (30, 40, 50, 60, 70, and 80 MPa) were built. Boundary and initial conditions are as follows: bm0 = 0.5 mm, β = 90°, Pb = 30 MPa, Pt = 10 MPa, Tb = 333.15 K, Tt = 303.15 K, and T0 = 303.15 K. Average Tf of fracture for caprock with different q is displayed in Figure 15(b). With the increase of q from 30 to 80 MPa, Tf shows 3-4 orders of magnitude reduction, which corresponds to a gradually better sequestration performance of caprock.


#### 4.2.4. Effect of Fluid Pressure DifferenceΔP

For CO2 sequestration in the storage layer, it is easier for CO2 to diffuse with a larger inlet fluid pressure. However, a large inlet pressure would produce a large Pf to the fracture in caprock, which would then result in decrease of σfne and increase of fracture permeability. Six models with different fluid pressure difference (ΔP = 5, 10, 15, 20, 25, and 30 MPa) were, respectively, set up. ΔP denotes the difference between Pb and Pt. Boundary and initial conditions are as follows: bm0 = 0.5 mm, q = 60 MPa, β = 90°, Pt = 10 MPa, Tb = 333.15 K, Tt = 303.15 K, and T0 = 303.15 K. Variation of average Tf with different ΔP is shown in Figure 15(c). In the range of ΔP from 5 to 30 MPa, Tf in stable state shows an increase of 2-3 orders of magnitude.


#### 4.2.5. Effect of TemperatureTC

TemperatureTC of the injected CO2 also plays a role in permeability of fracture in caprock. Three models with different Tb (308.15, 328.15, and 353.15 K) were, respectively, built to analyze the effect of TC. Both T0 and Tt were set to be 308.15 K, with bm0 = 0.5 mm, q = 60 MPa, β = 90°, Pb = 30 MPa, and Pt = 10 MPa. Calculation results are shown in Figure 15(d). It can be seen that the effect of Tb on average Tf of fracture is not as significant as that of β, q, and ΔP discussed above. Generally, when Tb>T0, both σfn and bh of fracture show an increasing trend. However, the results in Figure 15(d) show a relationship of Tf (328.15 K) > Tf (308.15 K) > Tf (353.15 K). This is because permeability of fracture is not only related to bh but also related to property parameters of fluid under coupled THM effects.


#### 4.2.6. Effect of Initial Aperturebm0

From (25), Tf of fracture is closely related to bm0. Four more models with different bm0 (0.05, 0.1, 0.3, and 0.5 mm) were set up to elucidate the effect of bm0 on permeability of fracture under coupled THM effects, withq = 60 MPa, β = 90°, Pb = 30 MPa, Pt = 10 MPa, Tb = 333.15 K, Tt = 303.15 K, and T0 = 303.15 K. Tf as a function of time with different bm0 is shown in Figure 15(e). Clearly, the larger bm0, the larger average Tf. In the range of bm0 from 0.05 to 0.5 mm, the average Tf at stable state undergoes 4-5 orders of magnitude increase, which would then degrade the CO2 sequestration performance. Besides, longer time is needed to attain a stable Tf with the decrease of bm0, as indicated by the dashed lines.


## 4.1. Model Setup

To quantitatively analyze permeability of the single fracture in caprock, a conceptual model is set up, as shown in Figure9. The model is square with the size of 10 m × 10 m. The fracture connects the upper and lower model boundary and passes through the center position of the model. The right, left, and lower model boundary are displacement constraint. A vertical load of q is applied on the upper model boundary. Fluid pressure of Pb and Pt is, respectively, applied at the lower and upper side of the fracture. Temperature of Tb and Tt is imposed on the lower and upper model boundary. Input parameters are listed in Table 2.Table 2
Parameters used in 2D THM coupled model.
Parameter
Value
Note
Density, kg/m3
2500
Rock matrix
Young’s modulus, GPa
100
Poisson’s ratio
0.3
Heat conductivity, W/(m × K)
2.57
Specific heat capacity, J/(m × K)
710
Thermal expansion, 1/K
6.0 × 10−7
Initial normal stiffness of fracture closure, GPa/m
60
Fracture
Normal stiffness of fracture opening, GPa/m
100
Initial aperture, m
b
m
0
Residual aperture, m
1 × 10−6
Maximum aperture, m
0.002
Fluid pressure at the lower side, MPa
P
b
Fluid pressure at the upper side, MPa
P
t
Density, kg/m3
Figure6(a)
Supercritical CO2
Dynamic viscosity, Pa⋅s
Figure6(b)
Compression coefficient, 1/Pa
Calculated with Figure6(a)
Thermal expansion, 1/K
Calculated with Figure6(b)Figure 9
2D conceptual model of a single fracture in caprock subjected to coupled THM effects, in whichβ denotes included angle between the fracture and horizontal direction, T0 is the initial temperature of rock matrix, and Tb is equal to the temperature of injected CO2.


## 4.2. Results and Discussion

First, variation ofσfne, Pf, and Tf of the fracture, as well as ρC and μC of the supercritical CO2 under THM coupled effects at different times for a vertical fracture (β = 90°) were analyzed. Then, permeability of the fracture in caprock with respect to different q, Pb, bm0, β, and Tb was, respectively, studied.


### 4.2.1. THM Coupled Effects on a Vertical Fracture

Boundary and initial conditions are as follows:β = 90°, bm0 = 0.5 mm, q = 60 MPa, Pb = 30 MPa, Pt = 10 MPa, Tb = 333.15 K, Tt = 303.15 K, and T0 = 303.15 K.At initial time, CO2 was continuously injected in the fracture through the lower fracture side with a constant Pb = 30 MPa. At different times, temperature distribution along the whole fracture length is shown in Figure 10, in which, x-axis (L0) represents the distance from the measure point to the lower fracture tip along the fracture direction (in the range of 0–10 m). During the fluid flow test, heat conduction occurs in rock due to higher temperature of injected CO2 (Tb) than that of the rock matrix (T0). The heat transfers from high-temperature to the low-temperature position and finally approaches a stable state. From Figure 10, with the increase of time, heat at the lower fracture segment gradually transfers to the upper position until a stable state is achieved at t = 300 days.Figure 10
Temperature distribution along the fracture at different times.During the heat conduction process, thermal expansion happens to the rock. Since both lateral boundaries (left and right) of the model are displacement constraint, the thermal expansion of rock leads to variation of normal stress in the fracture. Figures11(a)–11(c), respectively, show σfn, Pf, and σfne along the fracture length at different times. By comparing Figures 10 and 11(a), generally, a high temperature in fracture corresponds to a high total normal stress, and a steady stress field is achieved when the thermal field is steady. Variation of σfn in the fracture with time results mainly from the thermal field.Figure 11
Distribution of (a) total normal stress, (b) fluid pressure, and (c) effective normal stress along the fracture at different times.
(a)
σ
f
n
(b)
P
f
(c)
σ
f
n
eFluid pressurePf in the fracture as a function of L0 at different times is shown in Figure 11(b). With the increase of time, CO2 flows along the direction of pressure drops, and, for this case, Pf reaches steady values at t = 30 days. By using (5), variation of σfne of the fracture at different times can be obtained (Figure 11(c)). In a certain small range of L0, Pf⁡ is larger than σfn, which results in negative values of σfne and a larger bh than bm0, as shown by the dashed line in Figures 11(c) and 12. Besides, σfne increases gradually as L0 increases, which leads to corresponding decrease of bh. It can also be seen from Figure 12 that, in the range of L0 from 4 to 10 m, bh first increases and then decreases, finally reaching stable values, as indicated by the red arrows. The main reason for this phenomenon is due to earlier stable state of hydrodynamic field (t = 30 days) than that of the thermal field (t = 300 days).Figure 12
Variation ofbh along the fracture at different times.Distribution ofρC and μC of the supercritical CO2 along the fracture at different times is displayed in Figure 13. For L0 = 3–10 m, both property parameters show an ascending-descending variation before attaining constant values, which are marked by the red arrows.Figure 13
Variation of (a) density and (b) viscosity of supercritical CO2 along the fracture at different times.
(a)
Density
(b)
ViscosityUnder the coupled THM effects and taking variation of property parameters of CO2 with different temperature and pressure into account, the transmissivity Tf of the fracture can be evaluated (Figure 14). Once hydrodynamic field of the fracture is stable, Tf will reach constant values.Figure 14
Transmissivity distribution along the fracture at different times.


### 4.2.2. Effect of Included Angleβ

Due to complicated geological conditions of the caprock, included angle between fracture and the horizontal direction is various. Thus, it is of great significance to study the effect ofβ on permeability of the fracture under coupled THM effects. Six models with various β (45°, 51°, 59°, 68°, 79°, and 90°) were, respectively, set up. Boundary and initial conditions are as follows: bm0 = 0.5 mm, q = 60 MPa, Pb = 20 MPa, Pt = 10 MPa, Tb = 333.15 K, Tt = 303.15 K, and T0 = 303.15 K. Other input parameters were the same as those listed in Table 2. Figure 15(a) presents variation of average Tf of the fracture with different β. With the increase of time, the average Tf increases gradually and then attains a constant value. Under coupled THM effects, the average Tf experiences an increasing trend with the increase of β. The average Tf of the fracture with included angle β of 45°, 51°, 59°, 68°, 79°, and 90° is 1.66 × 10−5, 2.01 × 10−5, 2.55 × 10−5, 3.32 × 10−5, 4.15 × 10−5, and 4.55 × 10−5 m2/s. The average Tf for β = 90° increases by 1.74 times over that for β = 45°, which indicates a weaker sequestration performance of caprock with a larger β.Figure 15
Effect of (a) included angleβ, (b) overburden pressure q, (c) fluid pressure difference ΔP, (d) temperature Tb, and (e) initial aperture bm0 on average Tf of the fracture subject to coupled THM effects.
(a)
Included angleβ
(b)
Overburden pressureq
(c)
Fluid pressure differenceΔP
(d)
TemperatureTb
(e)
Initial aperturebm0


### 4.2.3. Effect of Overburden Pressureq

When the storage layer of CO2 is in different buried depths, the overburden pressure q of caprock is various, which would then largely impact the permeability of fracture. Six models with different q (30, 40, 50, 60, 70, and 80 MPa) were built. Boundary and initial conditions are as follows: bm0 = 0.5 mm, β = 90°, Pb = 30 MPa, Pt = 10 MPa, Tb = 333.15 K, Tt = 303.15 K, and T0 = 303.15 K. Average Tf of fracture for caprock with different q is displayed in Figure 15(b). With the increase of q from 30 to 80 MPa, Tf shows 3-4 orders of magnitude reduction, which corresponds to a gradually better sequestration performance of caprock.


### 4.2.4. Effect of Fluid Pressure DifferenceΔP

For CO2 sequestration in the storage layer, it is easier for CO2 to diffuse with a larger inlet fluid pressure. However, a large inlet pressure would produce a large Pf to the fracture in caprock, which would then result in decrease of σfne and increase of fracture permeability. Six models with different fluid pressure difference (ΔP = 5, 10, 15, 20, 25, and 30 MPa) were, respectively, set up. ΔP denotes the difference between Pb and Pt. Boundary and initial conditions are as follows: bm0 = 0.5 mm, q = 60 MPa, β = 90°, Pt = 10 MPa, Tb = 333.15 K, Tt = 303.15 K, and T0 = 303.15 K. Variation of average Tf with different ΔP is shown in Figure 15(c). In the range of ΔP from 5 to 30 MPa, Tf in stable state shows an increase of 2-3 orders of magnitude.


### 4.2.5. Effect of TemperatureTC

TemperatureTC of the injected CO2 also plays a role in permeability of fracture in caprock. Three models with different Tb (308.15, 328.15, and 353.15 K) were, respectively, built to analyze the effect of TC. Both T0 and Tt were set to be 308.15 K, with bm0 = 0.5 mm, q = 60 MPa, β = 90°, Pb = 30 MPa, and Pt = 10 MPa. Calculation results are shown in Figure 15(d). It can be seen that the effect of Tb on average Tf of fracture is not as significant as that of β, q, and ΔP discussed above. Generally, when Tb>T0, both σfn and bh of fracture show an increasing trend. However, the results in Figure 15(d) show a relationship of Tf (328.15 K) > Tf (308.15 K) > Tf (353.15 K). This is because permeability of fracture is not only related to bh but also related to property parameters of fluid under coupled THM effects.


### 4.2.6. Effect of Initial Aperturebm0

From (25), Tf of fracture is closely related to bm0. Four more models with different bm0 (0.05, 0.1, 0.3, and 0.5 mm) were set up to elucidate the effect of bm0 on permeability of fracture under coupled THM effects, withq = 60 MPa, β = 90°, Pb = 30 MPa, Pt = 10 MPa, Tb = 333.15 K, Tt = 303.15 K, and T0 = 303.15 K. Tf as a function of time with different bm0 is shown in Figure 15(e). Clearly, the larger bm0, the larger average Tf. In the range of bm0 from 0.05 to 0.5 mm, the average Tf at stable state undergoes 4-5 orders of magnitude increase, which would then degrade the CO2 sequestration performance. Besides, longer time is needed to attain a stable Tf with the decrease of bm0, as indicated by the dashed lines.


## 4.2.1. THM Coupled Effects on a Vertical Fracture

Boundary and initial conditions are as follows:β = 90°, bm0 = 0.5 mm, q = 60 MPa, Pb = 30 MPa, Pt = 10 MPa, Tb = 333.15 K, Tt = 303.15 K, and T0 = 303.15 K.At initial time, CO2 was continuously injected in the fracture through the lower fracture side with a constant Pb = 30 MPa. At different times, temperature distribution along the whole fracture length is shown in Figure 10, in which, x-axis (L0) represents the distance from the measure point to the lower fracture tip along the fracture direction (in the range of 0–10 m). During the fluid flow test, heat conduction occurs in rock due to higher temperature of injected CO2 (Tb) than that of the rock matrix (T0). The heat transfers from high-temperature to the low-temperature position and finally approaches a stable state. From Figure 10, with the increase of time, heat at the lower fracture segment gradually transfers to the upper position until a stable state is achieved at t = 300 days.Figure 10
Temperature distribution along the fracture at different times.During the heat conduction process, thermal expansion happens to the rock. Since both lateral boundaries (left and right) of the model are displacement constraint, the thermal expansion of rock leads to variation of normal stress in the fracture. Figures11(a)–11(c), respectively, show σfn, Pf, and σfne along the fracture length at different times. By comparing Figures 10 and 11(a), generally, a high temperature in fracture corresponds to a high total normal stress, and a steady stress field is achieved when the thermal field is steady. Variation of σfn in the fracture with time results mainly from the thermal field.Figure 11
Distribution of (a) total normal stress, (b) fluid pressure, and (c) effective normal stress along the fracture at different times.
(a)
σ
f
n
(b)
P
f
(c)
σ
f
n
eFluid pressurePf in the fracture as a function of L0 at different times is shown in Figure 11(b). With the increase of time, CO2 flows along the direction of pressure drops, and, for this case, Pf reaches steady values at t = 30 days. By using (5), variation of σfne of the fracture at different times can be obtained (Figure 11(c)). In a certain small range of L0, Pf⁡ is larger than σfn, which results in negative values of σfne and a larger bh than bm0, as shown by the dashed line in Figures 11(c) and 12. Besides, σfne increases gradually as L0 increases, which leads to corresponding decrease of bh. It can also be seen from Figure 12 that, in the range of L0 from 4 to 10 m, bh first increases and then decreases, finally reaching stable values, as indicated by the red arrows. The main reason for this phenomenon is due to earlier stable state of hydrodynamic field (t = 30 days) than that of the thermal field (t = 300 days).Figure 12
Variation ofbh along the fracture at different times.Distribution ofρC and μC of the supercritical CO2 along the fracture at different times is displayed in Figure 13. For L0 = 3–10 m, both property parameters show an ascending-descending variation before attaining constant values, which are marked by the red arrows.Figure 13
Variation of (a) density and (b) viscosity of supercritical CO2 along the fracture at different times.
(a)
Density
(b)
ViscosityUnder the coupled THM effects and taking variation of property parameters of CO2 with different temperature and pressure into account, the transmissivity Tf of the fracture can be evaluated (Figure 14). Once hydrodynamic field of the fracture is stable, Tf will reach constant values.Figure 14
Transmissivity distribution along the fracture at different times.


## 4.2.2. Effect of Included Angleβ

Due to complicated geological conditions of the caprock, included angle between fracture and the horizontal direction is various. Thus, it is of great significance to study the effect ofβ on permeability of the fracture under coupled THM effects. Six models with various β (45°, 51°, 59°, 68°, 79°, and 90°) were, respectively, set up. Boundary and initial conditions are as follows: bm0 = 0.5 mm, q = 60 MPa, Pb = 20 MPa, Pt = 10 MPa, Tb = 333.15 K, Tt = 303.15 K, and T0 = 303.15 K. Other input parameters were the same as those listed in Table 2. Figure 15(a) presents variation of average Tf of the fracture with different β. With the increase of time, the average Tf increases gradually and then attains a constant value. Under coupled THM effects, the average Tf experiences an increasing trend with the increase of β. The average Tf of the fracture with included angle β of 45°, 51°, 59°, 68°, 79°, and 90° is 1.66 × 10−5, 2.01 × 10−5, 2.55 × 10−5, 3.32 × 10−5, 4.15 × 10−5, and 4.55 × 10−5 m2/s. The average Tf for β = 90° increases by 1.74 times over that for β = 45°, which indicates a weaker sequestration performance of caprock with a larger β.Figure 15
Effect of (a) included angleβ, (b) overburden pressure q, (c) fluid pressure difference ΔP, (d) temperature Tb, and (e) initial aperture bm0 on average Tf of the fracture subject to coupled THM effects.
(a)
Included angleβ
(b)
Overburden pressureq
(c)
Fluid pressure differenceΔP
(d)
TemperatureTb
(e)
Initial aperturebm0


## 4.2.3. Effect of Overburden Pressureq

When the storage layer of CO2 is in different buried depths, the overburden pressure q of caprock is various, which would then largely impact the permeability of fracture. Six models with different q (30, 40, 50, 60, 70, and 80 MPa) were built. Boundary and initial conditions are as follows: bm0 = 0.5 mm, β = 90°, Pb = 30 MPa, Pt = 10 MPa, Tb = 333.15 K, Tt = 303.15 K, and T0 = 303.15 K. Average Tf of fracture for caprock with different q is displayed in Figure 15(b). With the increase of q from 30 to 80 MPa, Tf shows 3-4 orders of magnitude reduction, which corresponds to a gradually better sequestration performance of caprock.


## 4.2.4. Effect of Fluid Pressure DifferenceΔP

For CO2 sequestration in the storage layer, it is easier for CO2 to diffuse with a larger inlet fluid pressure. However, a large inlet pressure would produce a large Pf to the fracture in caprock, which would then result in decrease of σfne and increase of fracture permeability. Six models with different fluid pressure difference (ΔP = 5, 10, 15, 20, 25, and 30 MPa) were, respectively, set up. ΔP denotes the difference between Pb and Pt. Boundary and initial conditions are as follows: bm0 = 0.5 mm, q = 60 MPa, β = 90°, Pt = 10 MPa, Tb = 333.15 K, Tt = 303.15 K, and T0 = 303.15 K. Variation of average Tf with different ΔP is shown in Figure 15(c). In the range of ΔP from 5 to 30 MPa, Tf in stable state shows an increase of 2-3 orders of magnitude.


## 4.2.5. Effect of TemperatureTC

TemperatureTC of the injected CO2 also plays a role in permeability of fracture in caprock. Three models with different Tb (308.15, 328.15, and 353.15 K) were, respectively, built to analyze the effect of TC. Both T0 and Tt were set to be 308.15 K, with bm0 = 0.5 mm, q = 60 MPa, β = 90°, Pb = 30 MPa, and Pt = 10 MPa. Calculation results are shown in Figure 15(d). It can be seen that the effect of Tb on average Tf of fracture is not as significant as that of β, q, and ΔP discussed above. Generally, when Tb>T0, both σfn and bh of fracture show an increasing trend. However, the results in Figure 15(d) show a relationship of Tf (328.15 K) > Tf (308.15 K) > Tf (353.15 K). This is because permeability of fracture is not only related to bh but also related to property parameters of fluid under coupled THM effects.


## 4.2.6. Effect of Initial Aperturebm0

From (25), Tf of fracture is closely related to bm0. Four more models with different bm0 (0.05, 0.1, 0.3, and 0.5 mm) were set up to elucidate the effect of bm0 on permeability of fracture under coupled THM effects, withq = 60 MPa, β = 90°, Pb = 30 MPa, Pt = 10 MPa, Tb = 333.15 K, Tt = 303.15 K, and T0 = 303.15 K. Tf as a function of time with different bm0 is shown in Figure 15(e). Clearly, the larger bm0, the larger average Tf. In the range of bm0 from 0.05 to 0.5 mm, the average Tf at stable state undergoes 4-5 orders of magnitude increase, which would then degrade the CO2 sequestration performance. Besides, longer time is needed to attain a stable Tf with the decrease of bm0, as indicated by the dashed lines.


## 5. Conclusions

In this study, a new coupled CO2 flow, caprock deformation, and heat conduction finite element model is developed. A series of numerical calculations using COMSOL multiphysics software on caprock models with a single fracture subject to coupled THM effects were conducted. The main purpose is to elucidate the effect of fracture angle, overburden pressure, fluid pressure difference, injected CO2 temperature, and the initial aperture on single fracture permeability in caprock. The conclusions can be drawn as follows.(1)
A 2D single fracture FE model has been applied to verify the performance of the new model under hydromechanical coupled effects. Variation of fracture aperture along the fracture length shows a good agreement compared with the numerical and analytic results of Bower and Zyvoloski [19], which indicates a good applicability of the new coupled model.(2)
For a vertical fracture under coupled THM effects, with the increase of time, heat in fracture achieves a stable state att = 300 days, which corresponds to a steady stress state. CO2 flows along the direction of pressure drops, and Pf reaches steady values at t = 30 days. With the increase of L0, σfne increases while bh decreases gradually. For L0 = 3–10 m, both ρC and μC of the supercritical CO2 show an ascending-descending variation before attaining constant values.(3)
For all tested cases, with the increase of time, transmissivityTf of fracture increases before approaches a stable value. With the increase of β, average Tf of fracture shows an increasing trend. Tf decreases with the increase of overburden pressure. In the range of fluid pressure difference from 5 to 30 MPa, stable Tf shows an increase of 2-3 orders of magnitude. Tf is less dependent on TC of the injected CO2 compared with that for initial aperture which shows 4-5 orders of magnitude increase in the range of bm0 from 0.05 to 0.5 mm.We have tried in this paper to explain the coupled THM effects on transmissivity of a fracture in caprock. Clearly, more in-depth researches remain to be carried out on this issue. Our future works will focus on multiphase flow in fracture subjected to fully coupled THM processes to simulate the CO2 sequestration in caprock. Besides, FE models of fracture networks will be set up and solved by the computational simulation methods.



---

*Source: 1290748-2017-03-28.xml*