CO 2 Permeability Analysis of Caprock Containing a Single Fracture Subject to Coupled Thermal-Hydromechanical Effects

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 calculationswere performed on caprockmodels 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–30MPa), and 4-5 orders of magnitude with the initial fracture aperture (0.05–0.5mm), while it decreases by 3-4 orders of magnitude as overburden pressure increases from 30 to 80MPa. 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.


Introduction
To address the increasing concerns regarding carbon dioxide emission and its impact on climate change, CO 2 geological sequestration has become a promising approach [1][2][3][4][5].To reduce the atmosphere emissions, a large amount of CO 2 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 CO 2 leakage will occur [6][7][8].
Research on subsurface CO 2 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 CO 2 geological sequestration in storage reservoirs.
Numerical simulations have been widely used to evaluate CO 2 multiphase flow, diffusion, geomechanics, and chemical reactions during CO 2 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 FLAC 3D , have been adopted as a pragmatic approach for modeling coupled THM processes [12].Besides, TOUGHREACT and FLAC 3D 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 CO 2 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 CO 2 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 CO 2 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 CO 2 permeability of the fracture with respect to different fracture angle, overburden pressure, fluid pressure difference, injected CO 2 temperature, and initial aperture was, respectively, evaluated.In this study, these models were calculated under simplified conditions of singlephase flow and heat conduction alone in thermal field for brevity.

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 CO 2 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 CO 2 vary with the temperature and pressure.

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 Figure 2(a).The lower boundary of the model is displacement constraint.Due to continuous CO 2 injection, pressure is built up at the lower boundary of the caprock.Buried depth of caprock provides a vertical pressure of  at the upper boundary, and both lateral sides are subjected to an equal pressure of   .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). 0 and  C , respectively, denote the temperature of caprock matrix and injected CO 2 .For a homogeneous, isotropic, and elastic medium, the straindisplacement relation of the matrix can be expressed as where   and   are component of the total strain tensor produced by temperature and applied load, respectively.  and   are the component of displacement.The equilibrium equation is defined as where   and   denote the component of the total stress tensor produced by the temperature and applied load.  denotes the component of the body force.Then, the constitutive relation for the deformed caprock matrix becomes where  and  are Lame constants, ] is Poisson's ratio,  is Young's modulus, Δ is the temperature variation,  is the thermal expansion coefficient of rock, and   is the Kronecker delta.Combination of (1)-(3) yields the equilibrium differential equation, written as 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][21][22][23] or decrease in response to the normal loads [24][25][26][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 CO 2 flows through the fracture is not taken into account.
Figure 2: Mechanical model of caprock containing fractures.

Block Flow direction
Block Figure 3 shows a typical mechanical model of a fracture.When CO 2 flows in the fracture, a fluid pressure of  f is applied on the fracture surface. fn and  n , respectively, denote the normal contact stress and normal stiffness between two blocks.
According to the definition of effective stress in porous media, the effective stress in the fracture can be written as By using the distinct element code of UDEC, a simple description of the relation between  fne and the mechanical aperture  m was given [28], as indicated by the solid lines in Figure 4.In the range of residual aperture  res and maximum aperture  max ,  fne is linearly proportional to  m , written as where  m0 is the initial fracture aperture with no applied load.From Figure 4, when  fne > 0, fracture closure happens, and with the increase of  fne ,  m continues to decrease linearly until  res .When  m =  res , the fracture aperture keeps a constant value which no longer depends on  fne .When  fne < 0,  m increases with the increase of − fne until  max .When  m =  max ,  m remains unchanged.
In this study, for  fne < 0, the linear relation between  fne and  m in Figure 4 is still adopted to describe the variation of  m .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 where Δ c  is the fracture closure and  n0 is the initial normal stiffness with no applied load.
According to (7), the normal stiffness of the fracture can be determined as and the fracture aperture can be described as

Governing Equation for Fluid Flow.
For the two dimensional model described in Figure 2, 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 .
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 where  C is the density of CO 2 , V  denotes the flow velocity,  is the time, and  h is the equivalent hydraulic aperture.
According to Darcy's law, V  is given by where   is the permeability of the fracture.By neglecting the velocity head, the relationship between total head ℎ  and the osmotic pressure   can be written as where ℎ  is the head distribution of the fracture and   is the position head of the fracture corresponding to the global coordinate system.Taking a derivative of (12) yields the following equation: Substituting ( 13) into ( 11) and then into (10), we obtain where  f is the transmissivity of fracture and   is the local coordinate of fracture.
During the fracture deformation process (closure/ opening), the length of   is unchanged, while the fluid density and fracture aperture vary.Equation ( 14) can be rewritten as Assuming that the variation of fluid concentration is negligible, compression coefficient and the expansion coefficient of the fluid can be, respectively, written as [30] where  C and  C , respectively, denote the variation of fluid density induced by external pressure  and temperature .Therefore, the total variation of density can be described as The compression coefficient of the fracture can be written as Since the fracture length does not vary with  fne , ( 18) can be rewritten as where  n is the normal flexibility coefficient of the fracture, which can be written as 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 Substituting ( 21) into (19), we obtain Combination of ( 14), ( 15), (17), and ( 22) yields the following governing equation for fluid flow, expressed as Assuming that  h is equal to  m and fluid flow through the fracture can be described using Darcy's law, permeability of the fracture can be written as Mathematical Problems in Engineering Figure 6: Variations of (a) density [17] and (b) dynamic viscosity [18] of supercritical CO 2 with the temperature and pressure.
The transmissivity can be written as where  C is the dynamic viscosity of CO 2 and  presents the gravitational acceleration.

Governing Equation for Heat
Conduction.The injection of CO 2 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 CO 2 .Therefore, the temperature has a significant influence on permeability of the fracture.In this study, thermal effects produced by CO 2 flow through the fracture are neglected.Instead, the temperature of the injected CO 2 is regarded as a known condition which is imposed on the lower boundary of caprock.When the temperature of CO 2 is different from that of caprock, heat conduction phenomenon occurs.The governing equation can be written as where  denotes temperature variables and  s ,  s , and  s , respectively, denote the heat conduction coefficient, heat capacity, and density of rock matrix.

Property Parameters of Supercritical CO 2 .
Property parameters of the supercritical CO 2 injected in the storage layer vary with the pressure and temperature.In this study, the density and dynamic viscosity of the supercritical CO 2 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.

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.
An initial stress field of   = 21.0MPa 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   = 21.9MPa, and   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 × 10 6 MPa/m.For this case, numerical and analytic solutions of aperture variation along the fracture length at  = 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   [19] with the calculated results using COMSOL multiphysics in our study.
these curves fit well with each other, an evaluation coefficient  is put forward [33]: where  is the total number of measuring points, Δ hc is the fracture aperture change calculated in our study, and Δ h2 refers to fracture aperture change calculated using other methods.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  = 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  = 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.

CO 2 Permeability Analysis of
Single Fracture in Caprock   At initial time, CO 2 was continuously injected in the fracture through the lower fracture side with a constant  b = 30 MPa.At different times, temperature distribution along the whole fracture length is shown in Figure 10, in which, -axis ( 0 ) 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 CO 2 ( b ) than that of the rock matrix ( 0 ).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  = 300 days.
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.Figures 11(a)-11(c), respectively, show  fn ,  f , 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.
Fluid pressure  f in the fracture as a function of  0 at different times is shown in Figure 11(b).With the increase of time, CO 2 flows along the direction of pressure drops,  and, for this case,  f reaches steady values at  = 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  0 ,  f is larger than  fn , which results in negative values of  fne and a larger  h than  m0 , as shown by the dashed line in Figures 11(c) and 12. Besides,  fne increases gradually as  0 increases, which leads to corresponding decrease of  h .It can also be seen from Figure 12 that, in the range of  0 from 4 to 10 m,  h 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 ( = 30 days) than that of the thermal field ( = 300 days).Distribution of  C and  C of the supercritical CO 2 along the fracture at different times is displayed in Figure 13.For  0 = 3-10 m, both property parameters show an ascendingdescending variation before attaining constant values, which are marked by the red arrows.
Under the coupled THM effects and taking variation of property parameters of CO 2 with different temperature and pressure into account, the transmissivity  f of the fracture can be evaluated (Figure 14).Once hydrodynamic field of the fracture is stable,  f will reach constant values.15(e).Clearly, the larger  m0 , the larger average  f .In the range of  m0 from 0.05 to 0.5 mm, the average  f at stable state undergoes 4-5 orders of magnitude increase, which would then degrade the CO 2 sequestration performance.Besides, longer time is needed to attain a stable  f with the decrease of  m0 , as indicated by the dashed lines.

Conclusions
In this study, a new coupled CO 2 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 CO 2 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 at  = 300 days, which corresponds to a steady stress state.CO 2 flows along the direction of pressure drops, and  f reaches steady values at  = 30 days.With the increase of  0 ,  fne increases while  h decreases gradually.For  0 = 3-10 m, both  C and  C of the supercritical CO 2 show an ascendingdescending variation before attaining constant values.
(3) For all tested cases, with the increase of time, transmissivity  f of fracture increases before approaches a stable value.With the increase of , average  f of fracture shows an increasing trend. f decreases with the increase of overburden pressure.In the range of fluid pressure difference from 5 to 30 MPa, stable  f shows an increase of 2-3 orders of magnitude. f is less dependent on  C of the injected CO 2 compared with that for initial aperture which shows 4-5 orders of magnitude increase in the range of  m0 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 CO 2 sequestration in caprock.Besides, FE models of fracture networks will be set up and solved by the computational simulation methods.

2 Injection of supercritical C 2 Figure 1 :
Figure 1: Schematic of CO 2 injection in the presence of fractures within the caprock layer.

Figure 3 :
Figure 3: Mechanical model of the fracture.

Figure 4 :
Figure 4: Fracture aperture as a function of the effective normal stress.

Figure 5 :
Figure 5: Hydraulic model of fluid flow through a fracture.

Figure 9 :
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,  0 is the initial temperature of rock matrix, and  b is equal to the temperature of injected CO 2 .

Figure 6 2 Figure 10 :
Figure 10: Temperature distribution along the fracture at different times.

Figure 11 :
Figure 11: Distribution of (a) total normal stress, (b) fluid pressure, and (c) effective normal stress along the fracture at different times.

Figure 12 :
Figure 12: Variation of  h along the fracture at different times.

Figure 13 :Figure 14 :
Figure 13: Variation of (a) density and (b) viscosity of supercritical CO 2 along the fracture at different times.

Figure 15 :
Figure 15: Effect of (a) included angle , (b) overburden pressure , (c) fluid pressure difference Δ, (d) temperature  b , and (e) initial aperture  m0 on average  f of the fracture subject to coupled THM effects.

Table 2 .
Fluid pressure of  b and  t is, respectively, applied at the lower and upper side of the fracture.Temperature of  b and  t is imposed on the lower and upper model boundary.Input parameters are listed in 4.2.Results and Discussion.First, variation of  fne ,  f , and  f of the fracture, as well as  C and  C of the supercritical CO 2 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 ,  b ,  m0 , , and 4.1.Model Setup.To quantitatively analyze permeability of the single fracture in caprock, a conceptual model is set up, as shown in Figure 9.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  is applied on the upper model boundary.b was, respectively, studied.4.2.1.THM Coupled Effects on a Vertical Fracture.Boundary and initial conditions are as follows:  = 90 ∘ ,  m0 = 0.5 mm,  = 60 MPa,  b = 30 MPa,  t = 10 MPa,  b = 333.15K,  t = 303.15K, and  0 = 303.15K.

Table 2 :
Parameters used in 2D THM coupled model.
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:  m0 = 0.5 mm,  = 60 MPa,  b = 20 MPa,  t = 10 MPa,  b = 333.15K,t=303.15K,and0=303.15K.Other input parameters were the same as those listed in Table2.Figure15(a) presents variation of average  f of the fracture with different .With the increase of time, the average  f increases gradually and then attains a constant value.Under coupled THM effects, the average  f experiences an increasing trend with the increase of .The average  f 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 m 2 /s.The average  f for  = 90 ∘ increases by 1.74 times over that for  = 45 ∘ , which indicates a weaker sequestration performance of caprock with a larger .m0=0.5 mm,  = 90 ∘ ,  b = 30 MPa,  t = 10 MPa,  b = 333.15K,t=303.15K,and0=303.15K.Average  f of fracture for caprock with different  is displayed in Figure15(b).With the increase of  from 30 to 80 MPa,  f shows 3-4 orders of magnitude reduction, which corresponds to a gradually better sequestration performance of caprock.4.2.4.Effect of FluidPressure Difference Δ.For CO 2 sequestration in the storage layer, it is easier for CO 2 to diffuse with a larger inlet fluid pressure.However, a large inlet pressure would produce a large  f 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 (Δ = 5, 10, 15, 20, 25, and 30 MPa) were, respectively, set up.Δ denotes the difference between  b and  t .Boundary and initial conditions are as follows:  m0 = 0.5 mm,  = 60 MPa,  = 90 ∘ ,  t = 10 MPa,  b = 333.15K,t=303.15K,and0=303.15K.Variation of average  f with different Δ is shown in Figure15(c).In the range of Δ from 5 to 30 MPa,  f in stable state shows an increase of 2-3 orders of magnitude.4.2.5.Effect of Temperature   .Temperature  C of the injected CO 2 also plays a role in permeability of fracture in caprock.Three models with different  b (308.15,328.15,and353.15K)were,respectively,built to analyze the effect of  C .Both  0 and  t were set to be 308.15K,withm0=0.5 mm,  = 60 MPa,  = 90 ∘ ,  b = 30 MPa, and  t = 10 MPa.Calculation results are shown in Figure15(d).It can be seen that the effect of  b on average  f of fracture is not as significant as that of , , and Δ discussed above.Generally, when  b >  0 , both  fn and  h of fracture show an increasing trend.However, the results in Figure15(d) show a relationship of  f (328.15K)>f (308.15K)>f (353.15K).This is because permeability of fracture is not only related to  h but also related to property parameters of fluid under coupled THM effects.4.2.6.Effect of InitialAperture  0 .From (25),  f of fracture is closely related to  m0 .Four more models with different  m0 (0.05, 0.1, 0.3, and 0.5 mm) were set up to elucidate the effect of  m0 on permeability of fracture under coupled THM effects, with q = 60 MPa,  = 90 ∘ ,  b = 30 MPa,  t = 10 MPa,  b = 333.15K,  t = 303.15K, and  0 = 303.15K.  f as a function of time with different  m0 is shown in Figure 4.2.2.Effect of Included Angle .