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.

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 [15]. 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 [68].

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 , have been adopted as a pragmatic approach for modeling coupled THM processes [12]. Besides, TOUGHREACT and 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 Figure 2(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 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). and , 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 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 aswhere 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 becomeswhere 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 [2023] or decrease in response to the normal loads [2427]. 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.

Figure 3 shows a typical mechanical model of a fracture. When CO2 flows in the fracture, a fluid pressure of is applied on the fracture surface. and , 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 and the mechanical aperture was given [28], as indicated by the solid lines in Figure 4. In the range of residual aperture and maximum aperture , is linearly proportional to , written aswhere is the initial fracture aperture with no applied load.

From Figure 4, when , fracture closure happens, and with the increase of , continues to decrease linearly until . When , the fracture aperture keeps a constant value which no longer depends on . When , increases with the increase of until . When , remains unchanged.

In this study, for , the linear relation between and in Figure 4 is still adopted to describe the variation of . However, for , 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 aswhere is the fracture closure and is the initial normal stiffness with no applied load.

According to (7), the normal stiffness of the fracture can be determined asand the fracture aperture can be described as

2.2. 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 aswhere is the density of CO2, denotes the flow velocity, is the time, and is the equivalent hydraulic aperture.

According to Darcy’s law, is given bywhere is the permeability of the fracture.

By neglecting the velocity head, the relationship between total head and the osmotic pressure can be written aswhere 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 obtainwhere 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 and , 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 , (18) can be rewritten aswhere 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 is equal to and fluid flow through the fracture can be described using Darcy’s law, permeability of the fracture can be written as

The transmissivity can be written aswhere is the dynamic viscosity of CO2 and 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 aswhere denotes temperature variables and , , and , 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.

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.

An initial stress field of = 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 = 21.9 MPa, 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 × 106 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 these curves fit well with each other, an evaluation coefficient is put forward [33]:where is the total number of measuring points, is the fracture aperture change calculated in our study, and 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.

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 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. Fluid pressure of and is, respectively, applied at the lower and upper side of the fracture. Temperature of and is imposed on the lower and upper model boundary. Input parameters are listed in Table 2.

4.2. Results and Discussion

First, variation of , , and of the fracture, as well as and 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 , , , β, and was, respectively, studied.

4.2.1. THM Coupled Effects on a Vertical Fracture

Boundary and initial conditions are as follows: = 90°, = 0.5 mm, = 60 MPa, = 30 MPa, = 10 MPa, = 333.15 K, = 303.15 K, and = 303.15 K.

At initial time, CO2 was continuously injected in the fracture through the lower fracture side with a constant = 30 MPa. At different times, temperature distribution along the whole fracture length is shown in Figure 10, in which, -axis () 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 () than that of the rock matrix (). 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 , , and 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 in the fracture with time results mainly from the thermal field.

Fluid pressure in the fracture as a function of 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, reaches steady values at = 30 days. By using (5), variation of of the fracture at different times can be obtained (Figure 11(c)). In a certain small range of , is larger than , which results in negative values of and a larger than , as shown by the dashed line in Figures 11(c) and 12. Besides, increases gradually as increases, which leads to corresponding decrease of . It can also be seen from Figure 12 that, in the range of from 4 to 10 m, 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 and of the supercritical CO2 along the fracture at different times is displayed in Figure 13. For = 3–10 m, both property parameters show an ascending-descending variation before attaining constant values, which are marked by the red arrows.

Under the coupled THM effects and taking variation of property parameters of CO2 with different temperature and pressure into account, the transmissivity of the fracture can be evaluated (Figure 14). Once hydrodynamic field of the fracture is stable, will reach constant values.

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: = 0.5 mm, = 60 MPa, = 20 MPa, = 10 MPa, = 333.15 K, = 303.15 K, and = 303.15 K. Other input parameters were the same as those listed in Table 2. Figure 15(a) presents variation of average of the fracture with different β. With the increase of time, the average increases gradually and then attains a constant value. Under coupled THM effects, the average experiences an increasing trend with the increase of . The average 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 for β = 90° increases by 1.74 times over that for β = 45°, which indicates a weaker sequestration performance of caprock with a larger β.

4.2.3. Effect of Overburden Pressure

When the storage layer of CO2 is in different buried depths, the overburden pressure of caprock is various, which would then largely impact the permeability of fracture. Six models with different (30, 40, 50, 60, 70, and 80 MPa) were built. Boundary and initial conditions are as follows: = 0.5 mm, β = 90°, = 30 MPa, = 10 MPa, = 333.15 K, = 303.15 K, and = 303.15 K. Average of fracture for caprock with different is displayed in Figure 15(b). With the increase of from 30 to 80 MPa, 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

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 to the fracture in caprock, which would then result in decrease of 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 and . Boundary and initial conditions are as follows: = 0.5 mm, = 60 MPa, β = 90°, = 10 MPa, = 333.15 K, = 303.15 K, and = 303.15 K. Variation of average with different is shown in Figure 15(c). In the range of from 5 to 30 MPa, in stable state shows an increase of 2-3 orders of magnitude.

4.2.5. Effect of Temperature

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

4.2.6. Effect of Initial Aperture

From (25), of fracture is closely related to . Four more models with different (0.05, 0.1, 0.3, and 0.5 mm) were set up to elucidate the effect of on permeability of fracture under coupled THM effects, with q = 60 MPa, β = 90°, = 30 MPa, = 10 MPa, = 333.15 K, = 303.15 K, and = 303.15 K. as a function of time with different is shown in Figure 15(e). Clearly, the larger , the larger average . In the range of from 0.05 to 0.5 mm, the average 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 with the decrease of , 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 at = 300 days, which corresponds to a steady stress state. CO2 flows along the direction of pressure drops, and reaches steady values at = 30 days. With the increase of , increases while decreases gradually. For = 3–10 m, both and of the supercritical CO2 show an ascending-descending variation before attaining constant values.(3)For all tested cases, with the increase of time, transmissivity of fracture increases before approaches a stable value. With the increase of , average of fracture shows an increasing trend. decreases with the increase of overburden pressure. In the range of fluid pressure difference from 5 to 30 MPa, stable shows an increase of 2-3 orders of magnitude. is less dependent on of the injected CO2 compared with that for initial aperture which shows 4-5 orders of magnitude increase in the range of 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.

Conflicts of Interest

The authors declare that they have no competing interests regarding the publication of this paper.


The financial supports from the State Key Development Program for Basic Research of China (no. 2013CB036003), the Chinese Natural Science Foundation (no. 51374198), and the College Graduate Research and Innovation Projects of Jiangsu Province () are gratefully acknowledged.