Open AccessArticle

Fracture Evolution during CO₂ Fracturing in Unconventional Formations: A Simulation Study Using the Phase Field Method

Bing Yang

^1,2,3,*,

Qianqian Ren

^1,3,

Hai Huang

^3,4,*,

Haizhu Wang

⁵,

Yong Zheng

^3,4,

Liangbin Dou

^1,3,

Yanlong He

^3,4,

Wentong Zhang

⁶,

Haoyu Chen

^3,4 and

Ruihong Qiao

The Key Laboratory of Well Stability and Fluid & Rock Mechanics in Oil and Gas Reservoir of Shaanxi Province, Xi’an Shiyou University, Xi’an 710065, China

Shaanxi Key Laboratory of Carbon Dioxide Sequestration and Enhanced Oil Recovery, Xi’an 710065, China

College of Petroleum Engineering, Xi’an Shiyou University, Xi’an 710065, China

⁴

Shaanxi Key Laboratory of Advanced Stimulation Technology for Oil & Gas Reservoirs, Xi’an Shiyou University, Xi’an 710065, China

⁵

State Key Laboratory of Petroleum Resources and Prospecting, China University of Petroleum, Beijing 102249, China

⁶

College of New Energy, Xi’an Shiyou University, Xi’an 710065, China

Authors to whom correspondence should be addressed.

Processes 2024, 12(8), 1682; https://doi.org/10.3390/pr12081682

Submission received: 23 July 2024 / Revised: 6 August 2024 / Accepted: 8 August 2024 / Published: 12 August 2024

(This article belongs to the Special Issue Innovations in Hydraulic Fracturing Technology for Unconventional Reservoirs)

Download

Browse Figures

Versions Notes

Abstract

With the introduction of China’s “dual carbon” goals, CO₂ is increasingly valued as a resource and is being utilized in unconventional oil and gas development. Its application in fracturing operations shows promising prospects, enabling efficient extraction of oil and gas while facilitating carbon sequestration. The process of reservoir stimulation using CO₂ fracturing is complex, involving coupled phenomena such as temperature variations, fluid behavior, and rock mechanics. Currently, numerous scholars have conducted fracturing experiments to explore the mechanisms of supercritical CO₂ (SC-CO₂)-induced fractures in relatively deep formations. However, there is relatively limited numerical simulation research on the coupling processes involved in CO₂ fracturing. Some simulation studies have simplified reservoir and operational parameters, indicating a need for further exploration into the multi-field coupling mechanisms of CO₂ fracturing. In this study, a coupled thermo-hydro-mechanical fracturing model considering the CO₂ properties and heat transfer characteristics was developed using the phase field method. The multi-field coupling characteristics of hydraulic fracturing with water and SC-CO₂ are compared, and the effects of different geological parameters (such as in situ stress) and engineering parameters (such as the injection rate) on fracturing performance in tight reservoirs were investigated. The simulation results validate the conclusion that CO₂, especially in its supercritical state, effectively reduces reservoir breakdown pressures and induces relatively complex fractures compared with water fracturing. During CO₂ injection, heat transfer between the fluid and rock creates a thermal transition zone near the wellbore, beyond which the reservoir temperature remains relatively unchanged. Larger temperature differentials between the injected CO₂ fluid and the formation result in more complicated fracture patterns due to thermal stress effects. With a CO₂ injection, the displacement field of the formation deviated asymmetrically and changed abruptly when the fracture formed. As the in situ stress difference increased, the morphology of the SC-CO₂-induced fractures tended to become simpler, and conversely, the fracture presented a complicated distribution. Furthermore, with an increasing injection rate of CO₂, the fractures exhibited a greater width and extended over longer distances, which are more conducive to reservoir volumetric enhancement. The findings of this study validate the authenticity of previous experimental results, and it analyzed fracture evolution through the multi-field coupling process of CO₂ fracturing, thereby enhancing theoretical understanding and laying a foundational basis for the application of this technology.

Keywords:

unconventional resource; CO₂ fracturing; multi-field coupling; phase field

1. Introduction

China boasts abundant tight oil and gas resources with enormous development potential. In the process of developing unconventional resources, hydraulic fracturing technology occupies an irreplaceable position [1,2,3]. While hydraulic fracturing offers significant advantages for enhancing oil and gas production, it still faces numerous unavoidable challenges in reservoir management and environmental protection. Primarily, large-scale volumetric fracturing consumes substantial water resources, which poses significant constraints in the arid and water-deficient regions of northwestern China [4,5,6]. Furthermore, the addition of a large quantity of chemicals to fracturing fluids can cause contamination to the upper aquifers and surface environments. Moreover, significant investment is required for purification of returned fracturing fluids. Additionally, the presence of water can induce clay mineral expansion, increasing the risk of reservoir contamination and affecting the effectiveness of fracturing for production enhancement [7,8]. These issues directly constrain the development of unconventional oil and gas resources.

In recent years, with the worsening global climate and the gradual introduction of “dual carbon” goals, the application of CO₂ in oil and gas development is considered a promising approach [9,10,11]. Given the formation depth and geological conditions of unconventional oil and gas, liquid CO₂ injected into reservoirs would convert to a supercritical state. Consequently, researchers have shifted their focus to utilizing supercritical CO₂ (SC-CO₂) for the development of unconventional resources [12].

Utilizing SC-CO₂ as a fracturing fluid for the stimulation of unconventional reservoirs offers numerous advantages [13]. CO₂, being a water-free working fluid, can mitigate clay expansion, thus preventing reservoir contamination [10]. SC-CO₂ readily penetrates the micropores in a reservoir, enhancing the complexity of artificial fractures. Additionally, SC-CO₂ exhibits stronger adsorption on rock surfaces compared with methane, which assists in displacing natural gas [14]. Post CO₂ fracturing, it is easier to flow back or even unnecessary to do so, allowing for direct in situ carbon storage. This method would contribute to achieving the “dual carbon” goals [15,16,17]. Numerous researchers have conducted numerous hydraulic fracturing experiments, particularly focusing on the initiation and propagation of fractures using SC-CO₂. While some scholars have engaged in numerical simulation studies, there are challenges such as simplified model treatments and the lack of consideration for multi-field coupling.

For compressible fluids like CO₂, their properties are significantly influenced by changes in temperature and pressure. When conducting simulation studies, it is essential to consider these variations in fluid properties, which pose new demands on related research efforts. In 2014, Fang et al. [18] conducted a comparative study using the Universal Distinct Element Code (UDEC) on the fracturing processes of supercritical water, slickwater, and CO₂ in shale. The results indicated that the viscosity of fracturing fluids affects the formation of fractures. Zhang and Liu et al. [19,20,21], based on finite element algorithms and incorporating partial secondary development, established a numerical model for CO₂ fracturing. They conducted a comparative study on the effects of heavy oil, water, nitrogen, and CO₂ on fracturing. The study found that SC-CO₂ fracturing exhibits lower initiation pressure and induces more complex fractures. Subsequently, Zhang et al. [22] developed a mechanical damage model based on finite elements, combining the maximum tensile stress criterion and Mohr–Coulomb criterion, to simulate the fracture propagation behavior in heterogeneous rocks. They investigated the effects of rock heterogeneity, the Biot coefficient, confining pressure, interfacial tension, and the fracturing fluid type on the fracturing process. In 2019, Yan et al. [23] employed the extended finite element method (XFEM) to establish a fluid-structure coupling model for simulating the process of SC-CO₂ fracturing in coal seams. Using this model, they delineated the fracturing stages into SC-CO₂ injection and CO₂ phase change-induced fracture propagation stages. Also in 2019, Mollaali et al. [24] utilized phase field methods to simulate CO₂ fracturing processes in porous media. The model incorporated modified Darcy’s law to simulate fluid flow and validated several numerical examples. Then, He et al. [25] integrated the discontinuous deformation method (DDM) and finite volume method (FVM) to construct a three-dimensional model for fracture propagation and reservoir rock deformation. They analyzed the influences of geological and structural factors on fracture propagation during SC-CO₂ fracturing.

From the literature review, it is evident that there are relatively few numerical studies on SC-CO₂ fracturing. Numerical research methodologies for CO₂ fracturing generally fall into four categories: (1) finite element-based fracturing simulations; (2) the extended finite element method (XFEM); (3) phase field methods; and (4) the discontinuous deformation method (DDM). The simulation processes often overlook heat transfer between CO₂ and the formations, as well as changes in the fluid properties of CO₂. Research into the fracture propagation characteristics of SC-CO₂ fracturing based on thermo-hydro-mechanical coupling remains an area requiring further investigation.

In this study, we establish a multi-field coupling model for CO₂ fracturing based on the phase field differential method. By comparing the simulations of water and CO₂ fracturing processes, the coupling characteristics are analyzed. We further investigate the influence mechanisms of different factors such as the in situ stress, injection temperature, and displacement on CO₂-induced fracturing. The results deepen our understanding of the multi-field coupling mechanisms during fracturing processes, laying a theoretical and parameter-based foundation for the application of this technology.

2. Simulation Methodology of Hydraulic Fracturing

During CO₂ fracturing, fracture propagation is a coupled phenomenon involving fluid flow in porous media and fracture areas, rock deformation, and fracture propagation within the rock mass. This study will derive governing equations for porous media deformation and fracture propagation based on phase field theory, aiming to establish a thermo-hydro-mechanical coupling model for fracturing which accounts for CO₂ property changes and heat transfer.

2.1. Numerical Theory Based on Phase Field Method

(1): Brittle Fracture Theory

The phase-field damage model is closely related to strain-based damage models, but its foundation differs. This model is a regularization expression based on Griffith’s classical brittle fracture theory variational formulation. Assuming an arbitrary body

Ω \subset ℝ^{d} (d \in \{1, 2, 3\}

with an external boundary

\partial Ω

and an internal discontinuous boundary Γ, the displacement of body Ω at time (t) is represented by

u (x, t) \subset ℝ^{d}

, where (

x

) is the position vector. It is also considered that the volume force

b (x, t) \subset ℝ^{d}

acts on the body Ω, and the traction force

f (x, t)

acts on the boundary

\partial Ω_{t_{i}}

[26].

The energy required to create a fracture surface per unit area is defined as the critical fracture energy density G_c, often referred to as the critical energy release rate. The total potential energy Ψ_pt(u, Γ) can be expressed in terms of the elastic energy ψ_ε(ε), fracture energy, and external work energy [27]:

Ψ_{p t} (u, Γ) = \int_{Ω} ψ_{ε} (ε) d Ω + \int_{Γ} G_{c} d S - \int_{Ω} b \cdot u d Ω - \int_{\partial Ω_{t_{i}}} f \cdot u d S

(1)

The linear strain tensor

ε = ε (u)

is given by the following equation:

ε_{i j} = \frac{1}{2} (\frac{\partial u_{i}}{\partial x_{j}} + \frac{\partial u_{j}}{\partial x_{i}})

(2)

Assuming the linear elastic material is isotropic, the elastic energy density ψ_ε(ε) can be represented as follows:

ψ_{ε} (ε) = \frac{1}{2} λ_{L} ε_{i i} ε_{j j} + μ_{L} ε_{i j} ε_{i j}

(3)

where,

λ_{L}

and

μ_{L}

are Lamé coefficients.

(2): Phase Field Approximation of Fracture Energy

We define a phase field scalar ϕ(x, t) ∈ [0, 1] to approximate fracture presence, where ϕ = 1 represents a fracture and ϕ = 0 denotes an undamaged rock matrix. For two-dimensional and three-dimensional problems, the fracture surface density per unit volume is given by [28]

γ (ϕ, \nabla ϕ) = \frac{ϕ^{2}}{2 l_{0}} + \frac{l_{0}}{2} \frac{\partial ϕ}{\partial x_{i}} \frac{\partial ϕ}{\partial x_{i}}

(4)

where

l_{0}

is a parameter controlling the transition of the phase field from 0 to 1, also known as the length scale parameter. This parameter reflects the shape of the fracture; larger values indicate wider fractures, and vice versa.

Based on Equation (4), the fracture energy can be expressed as follows:

\int_{Γ} G_{c} d S = \int_{Ω} G_{c} [\frac{ϕ^{2}}{2 l_{0}} + \frac{l_{0}}{2} \frac{\partial ϕ}{\partial x_{i}} \frac{\partial ϕ}{\partial x_{i}}] d Ω

(5)

The fracture surface energy is derived from elastic energy, indicating that the evolution of the phase field is driven by elastic energy. To ensure fractures are driven solely by tensile loads, the elastic energy is decomposed into tensile and compressive components. The strain tensor ε is decomposed as follows:

\{\begin{cases} ε_{+} = \sum_{a = 1}^{d} {〈ε_{a}〉}_{+} n_{a} \otimes n_{a} \\ ε_{-} = \sum_{a = 1}^{d} {〈ε_{a}〉}_{-} n_{a} \otimes n_{a} \end{cases}

(6)

where ɛ₊ and ɛ₋ are the tensile and compressive strain tensors, respectively, and ε is the principal strain, with n_a as the direction vector. By introducing the decomposed strain tensor, the elastic energy density is expressed as follows:

\{\begin{cases} ψ_{ε}^{+} (ε) = \frac{λ}{2} {〈t r (ε)〉}_{+}^{2} + μ t r (ε_{+}^{2}) \\ ψ_{ε}^{-} (ε) = \frac{λ}{2} {〈t r (ε)〉}_{-}^{2} + μ t r (ε_{-}^{2}) \end{cases}

(7)

By assuming only the tensile part of the elastic energy density is influenced by the phase field, the following equation is then used to simulate stiffness degradation [29]:

ψ_{ε} (ε) = [(1 - k) {(1 - ϕ)}^{2} + k] ψ_{ε}^{+} (ε) + ψ_{ε}^{-} (ε)

(8)

where k is a model parameter which prevents singularities in the elastic energy density as ϕ approaches 1. Additionally, it is required that k > 0 and k ≪ 1.

(3): Governing Equations of Rock Mechanics

Aside from the total potential energy, the kinetic energy

ψ_{k}

of the body

Ω

also needs to be considered:

ψ_{k} (\dot{u}) = \frac{1}{2} \int_{Ω} ρ_{s} {\dot{u}}_{i} {\dot{u}}_{i} d Ω

(9)

where

\dot{u} = \frac{\partial u}{\partial t}

and

ρ_{s}

is the density of the rock in kg/m³.

The total Lagrangian energy functional

ℑ

can be expressed as the sum of the fracture energy (Equation (5)), elastic energy (Equation (8)), kinetic energy (Equation (9)), and the phase field approximation of the external load potential energy:

\begin{array}{l} ℑ = \frac{1}{2} {\int_{Ω} ρ}_{s} {\dot{u}}_{i} {\dot{u}}_{i} d Ω - \int_{Ω} \{[(1 - k) {(1 - ϕ)}^{2} + k] ψ_{ε}^{+} (ε) + ψ_{ε}^{-} (ε)\} d Ω \\ - \int_{Ω} G_{c} [\frac{ϕ^{2}}{2 l_{0}} + \frac{l_{0}}{2} \frac{\partial ϕ}{\partial x_{i}} \frac{\partial ϕ}{\partial x_{i}}] d Ω + \int_{Ω} b_{i} u_{i} d Ω + \int_{\partial Ω_{t_{i}}} f_{i} u_{i} d S \end{array}

(10)

Furthermore, by computing the first variation of the functional

ℑ

and setting it to zero, we obtain

\{\begin{array}{r} \frac{\partial σ_{i j}}{\partial x_{i}} + b_{i} = ρ_{s} {\ddot{u}}_{i} \\ [\frac{2 l_{0} (1 - k) ψ_{ε}^{+}}{G_{c}} + 1] ϕ - l_{0}^{2} \frac{\partial^{2} ϕ}{\partial x_{i}^{2}} = \frac{2 l_{0} (1 - k) ψ_{ε}^{+}}{G_{c}} \end{array}

(11)

where

{\ddot{u}}_{i} = \frac{\partial^{2} u}{\partial t^{2}}

, σ_ij represents the Cauchy stress tensor components, given by

σ_{i j} = [(1 - k) {(1 - ϕ)}^{2} + k] \frac{\partial ψ_{ε}^{+}}{\partial ε_{i j}} + \frac{\partial ψ_{ε}^{-}}{\partial ε_{i j}}

(12)

Equation (12) can be further expressed as follows:

σ = [(1 - k) {(1 - ϕ)}^{2} + k] [λ {〈t r (ε)〉}_{+} I + 2 μ ε_{+}] + λ {〈t r (ε)〉}_{-} I + 2 μ ε_{-}

(13)

where I is the unit tensor

\in ℝ^{d \times d}

To ensure a monotonic increase in the phase field, irreversibility conditions are introduced during the compression or unloading processes. Thus, a strain history field H is introduced:

H (x, t) = \max_{s \in [0, t]} ψ_{ε}^{+} (ε (x, s))

(14)

By replacing

ψ_{ε}^{+}

in Equation (11) with

H (x, t)

, the strong form is derived:

\{\begin{array}{r} \frac{\partial σ_{i j}}{\partial x_{i}} + b_{i} = ρ_{s} {\ddot{u}}_{i} \\ [\frac{2 l_{0} (1 - k) H}{G_{c}} + 1] ϕ - l_{0}^{2} \frac{\partial^{2} ϕ}{\partial x_{i}^{2}} = \frac{2 l_{0} (1 - k) H}{G_{c}} \end{array}

(15)

(4): CO₂ Fluid and Heat Transfer Control Equation

The governing equations for fluid in porous media are derived from mass conservation, momentum conservation, and state equations. The mass conservation equation is given by

\partial t (φ ρ_{f}) + \nabla \cdot (ρ_{f} v) = 0 in Ω

(16)

- ρ_{f} v \cdot n = Q g on Γ

(17)

where φ denotes the porosity of the porous medium (the volume fraction occupied by fluid),

ρ_{f}

represents the fluid density,

v

is the Darcy velocity vector, and Q_g denotes the fluid sources, representing the volumetric flow rate per unit volume. Additionally, assuming the rock is saturated with fluid, the fluid content per unit volume within the rock is denoted as

φ ρ_{f}

The momentum balance is formulated in the form of Darcy’s law, implying a linear relationship between the fluid velocity and hydraulic gradient:

v = - \frac{k_{p}}{μ} \nabla p - ρ_{f} g \nabla z

(18)

where k_p = k(d) represents the rock permeability, μ is the dynamic viscosity of the fluid, and g and z denote the gravitational acceleration and depth, respectively. In this study, we assume that porosity is independent of stress conditions, and hence the

- ρ_{f} g \nabla z

term can be neglected.

The first term on the left side of Equation (16), representing the rate of change of the fluid content, can be expressed as follows:

\partial t (φ ρ_{f}) = ρ_{f} \partial t φ + φ \partial t ρ_{f} = ρ_{f} \partial t ε v + φ \partial t ρ_{f}

(19)

Here, it is assumed that the rate of change of the pore volume equals the rate of volumetric strain, as strain given by

ε v = \nabla \cdot u

Under isothermal conditions, the gas density varies significantly with the pressure and can be determined using an equation of state (EOS). Currently, the widely used EOS for determining the properties of CO₂ is the Span–Wagner (S-W) equation [30], which is based on the Helmholtz free energy. The density and pressure of CO₂ are related by the following formula:

p = (1 + δ ϕ_{δ}^{r}) ρ_{f} R T

(20)

where R is the universal gas constant and

ϕ_{δ}^{r}

is the derivative of the remaining part of the total expression of the Helmholtz energy

ϕ^{r}

with respect to the reduced density

δ

. Here, we have

\begin{array}{l} ϕ_{δ}^{r} = \sum_{i = 1}^{7} n_{i} d_{i} δ^{d_{i} - 1} τ^{t_{i}} + \sum_{i = 8}^{34} n_{i} e^{- δ^{c_{i}}} [δ^{d_{i} - 1} τ^{t_{i}} (d_{i} - c_{i} δ^{c_{i}})] + \\ \sum_{i = 35}^{39} n_{i} δ^{d_{i}} τ^{t_{i}} e^{- α_{i} {(δ - ε_{i})}^{2} - β_{i} {(τ - γ_{i})}^{2}} [\frac{d_{i}}{δ} - 2 α_{i} (δ - ε_{i})] + \sum_{i = 40}^{42} n_{i} [Δ^{b_{i}} (ψ + δ \frac{\partial ψ}{\partial δ}) + \frac{\partial Δ^{b_{i}}}{\partial δ} δ ψ] \end{array}

(21)

where

Δ = {(1 - τ) + Ai {[{(δ - 1)}^{2}]}^{\frac{1}{2 β i}}}^{2} + B i {[{(δ - 1)}^{2}]}^{a i}

In Equation (21),

ρ_{c}

and

T_{c}

denote the critical density and temperature, respectively, while

δ = ρ_{f} / ρ_{c}

and

τ = T / T_{c}

are simplification parameters. For the definitions and numerical values of other constant parameters

n_{i}, d_{i}, t_{i}, c_{i}, α_{i}, β_{i}, γ_{i}, b_{i}, A_{i}, B_{i}, a_{i}

, please refer to Reland Span’s work [30].

By substituting Equations (18)–(20) into Equations (16) and (17), the flow control equation for CO₂ is obtained:

φ \partial t ρ_{f} + ρ_{f} \partial t ε v - \nabla \cdot (ρ_{f} \frac{k (d)}{μ} \nabla p) = 0, in Ω ρ_{f} v \cdot n = - Q g, on Γ p = p D, on Γ p

(22)

where

Γ p = \partial Ω \ Γ

, which represents the entire body domain

Ω

except for the internal injection boundary

Γ

, usually called a pore space domain

Γ p

Assuming the fluid flow and heat transfer process in porous media adhere to local thermal equilibrium conditions, neglecting the effects of fluid work and viscous dissipation on heat transfer, the heat transfer equation between the fluid and the rock in porous media is given by [31]

C_{p, e f f} \frac{\partial T}{\partial t} + ρ_{f} C_{p, f} u \cdot \nabla T + \nabla q = Q

(23)

q = - k_{e f f} \nabla T

(24)

C_{p, e f f} = (1 - φ) ρ_{s} C_{p, s} + φ ρ_{f} C_{p, f}

(25)

where

C_{p, e f f}

is the effective volumetric specific heat capacity at a constant pressure in J/(m³·K);

C_{p, f}

is the constant pressure specific heat capacity of the fluid in J/(m³·K);

C_{p, s}

is the constant-pressure heat capacity of the rock in J/(m³·K);

q

is the heat flux in W/m²; and Q is the heat source in W/m³. In Equation (24),

k_{e f f}

can be expressed as follows:

k_{e f f} = (1 - φ) k_{s} + φ k_{f}

(26)

where

k_{s}

and

k_{f}

represent the thermal conductivity of the rock and fluid, respectively.

(5): Determination of Parameter $l_{0}$

In phase field models, the selection of the length scale parameter remains an open question. Borden provided an analytical solution for the critical tensile stress which a one-dimensional rod can withstand under tension [29]

σ_{c r} = \frac{9}{16} \sqrt{\frac{E G_{c}}{3 l_{0}}}

(27)

where E is the Young’s modulus and G_c is the critical energy release rate. When l₀ approaches zero, this indicates the appearance of a crack tip. There is an obvious singularity which leads to a physically meaningless infinite tensile strength.

When all other parameters except for

l_{0}

are known, it can be concluded that

l_{0} = \frac{27 E G_{c}}{256 σ_{c r}^{2}}

(28)

The critical energy release rate G_c, Young’s modulus E, and critical stress

σ_{c r}^{2}

can be estimated through some conventional experiments.

2.2. Solution Method

The simulation method was implemented in COMSOL Multiphysics 6.2 while incorporating six coupled solving modules: Solid Mechanics, Darcy Flow, Heat Transfer, Phase Field, History of Strain, and Storage. All modules were developed based on spatial domain standard finite element discretization and time domain finite difference discretization, formulated in strong form, and integrated for coupled solving. The Solid Mechanics module utilized a linear elastic material model for displacement calculations. Variations in the CO₂ fluid properties were implemented through embedded S-W control equations. The Darcy Flow module employed transient formulations of Darcy’s law for fluid pressure calculations. The Heat Transfer module utilized a porous media heat transfer model to compute temperature distributions. By modifying the Helmholtz equation in COMSOL, a phase field model was established, and the parameters of

ϕ

and H were solved based on the constructed phase field model and historical strain model. The coupled solution process of the model is shown in Figure 1.

2.3. Verification of Model and Analysis of Grid Dependence

This section is based on the classic two-dimensional flat plate tensile model [32], and in it, we calculate the propagation process of the central crack to verify the feasibility of the model. The schematic diagram of the model is shown in Figure 2. The material parameters were as follows: Young’s modulus E = 210 GPa, Poisson’s ratio

ν

= 0.3, critical energy release rate G_c = 2700 J/m², and length of edge L = 1 mm. To prevent singular stiffness matrices, the chosen model parameter was k = 1 × 10⁻⁹. The length scale parameters were

l_{0}

= 7.5 × 10⁻³ mm and 1.5 × 10⁻² mm, with a grid size h of approximately 3.5 × 10⁻³ mm. A uniform increasing displacement load u_y = 1 × 10⁻⁵ mm/s was applied to the top of the model, while the left and right boundaries had zero displacement. The stretching process was achieved through the phase field method, and the distribution of central cracks at different displacements was obtained under two different length scale parameter conditions. The results are shown in Figure 3. The results are consistent with our expectations; smaller length scale parameters corresponded to narrower crack propagation widths, whereas larger parameters resulted in wider cracks. The obtained crack patterns aligned with those reported in [28,33], where the cracks extended horizontally and reached the boundaries. As depicted in Figure 4, the displacement–load curve during the stretching process closely matched the results from the literature, confirming the model’s accuracy. However, minor errors were observed during model solving, attributable to discrepancies between the algorithms and slight differences in the grid cell sizes compared with the original model.

Furthermore, this hydraulic fracturing model operates on calculations based on grid nodes, and hence the size of the grid cells may influence the simulation results. Typically, smaller grid cell sizes and tolerance factors facilitate better convergence of the model. However, to optimize computational costs, it is often necessary to reasonably increase the grid cell sizes and tolerance factors. During testing of the fracturing model, various grid cell sizes (h = 5 × 10⁻⁶, 1 × 10⁻⁵, 2 × 10⁻⁵ m) and tolerance factors (1 × 10⁻³, 1 × 10⁻⁵, 1 × 10⁻⁷) were employed. The results indicate that these values did not affect the pattern of fracture propagation, demonstrating good consistency in the outcomes. The geometric dimensions of the hydraulic fracturing model were fixed. Grid sizes of h = 2 × 10⁻⁵ m and a tolerance factor of 1×10⁻³ were selected for subsequent simulations. Therefore, the irrelevance of the grid during the simulation process was validated.

3. Simulation of Water and SC-CO₂ Fracturing

3.1. Geometrical Model

According to the stress conditions of the tight reservoirs (the Chang 6 group of the Yanchang oilfield), the vertical principal stress is higher than the horizontal principal stresses. The previous experimental results also showed that the fractures induced by SC-CO₂ fracturing predominantly propagate along the direction perpendicular to the horizontal principal stress and exhibit mostly vertical fractures [34]. Therefore, selecting the fracture distribution of any horizontal cross-section of the rock sample (perpendicular to the vertical principal stress direction) can roughly represent the morphology of the fractures. In addition, to save computational costs, this study will establish a two-dimensional model of CO₂ fracturing which can be regarded as a selected cross-section perpendicular to the vertical principal stress direction. This model has a theoretical basis and scientific significance.

As shown in Figure 5, the outer boundary edge length of the model is L, the center position is the wellbore with a radius r, and the area between the wellbore and the boundary is the rock matrix. The left and lower boundaries of the model are slip boundaries, with the maximum geostress

σ_{m a x}

applied to the top boundary and the minimum geostress

σ_{m i n}

applied to the right boundary. The pore pressure of the formation is

P_{r}

, the bottomhole pressure is

P_{w}

, and in the initial state,

P_{r} = P_{w}

. Fluid was injected at a constant volume flow rate. Assuming the outer boundary of the model is adiabatic, and no seepage occurs, the initial temperature of the formation is

T_{0}

, and the injection temperature of the fluid is

T_{i n}

. The green line is a reference line for studying the evolution of variables such as the temperature in the coupling process and the distance from the wellhole.

3.2. Parameters of Model

The reservoir parameters utilized in this model were derived from field data and various rock property testing experiments. The simulation replicated conditions akin to laboratory fracturing experiments, and the initial reservoir and fluid parameters set for the model are shown in Table 1.

Since water is a micro-compressed fluid; its physical properties are less affected by changes in temperature and pressure. Under the simulation conditions, it was assumed that the properties of the water did not change. However, for CO₂ fluids, its properties are significantly affected by changes in temperature and pressure, and thus the property variations during the fracturing process cannot be ignored. According to the S-W equation [30], the primary properties of CO₂ (density, viscosity, thermal conductivity, and specific heat capacity) exhibited significant variations with changes in temperature and pressure, as illustrated in Figure 6. It can be observed that the CO₂ density, viscosity, and thermal conductivity were highly sensitive to temperature and pressure changes. Within the temperature range from 273 to 473 K and the pressure range from 0 to 100 MPa, the CO₂ density ranged from nearly 0 to close to 1200 kg/m³. The viscosity and thermal conductivity also underwent significant changes on multiple scales.

3.3. Fracture Evolution during Water and SC-CO₂ Fracturing

Figure 7 and Figure 8 illustrate the primary processes of water and SC-CO₂ fracturing in tight reservoirs. It can be observed that both fracturing processes exhibited similarities. Initially (t = 0), due to the influence of the maximum in situ stress in the y direction, a stress concentration appeared near the wellbore with an approximately symmetric distribution along the x direction. As the fluid continued to be injected (t = 0.5t_f, where t_f denotes the total fracturing duration), the stress concentration in the x direction diminished and gradually appeared symmetrically along the y direction. When the fluid pressure reached the rock’s initiation pressure, fractures initiated near the wellbore along the direction of the maximum in situ stress, followed by the propagation of primary fractures along the y direction. This outcome aligns with classical fracture mechanics theory [35]. Additionally, it is noted that the time required from fracture initiation to complete propagation was notably short, constituting less than one tenth of the entire fracturing process duration, which is consistent with the experimental findings in previous research [34,36].

In addition, significant differences were observed in the hydraulic fracturing with water and SC-CO₂. In the water fracturing process (Figure 7), fractures initiated earlier compared with those induced by SC-CO₂, which was due to the fact that water exhibits lower compressibility than SC-CO₂. There were also differences in the induced fractures; water induced a single fracture along the distribution of the maximum in situ stress which did not fully penetrate the entire rock matrix. In contrast, SC-CO₂ induced a fracture which extended along the y direction (as depicted in Figure 8), fully penetrating the rock matrix to the boundary, and it additionally generated a short microfracture along the x direction near the wellbore. These simulation results are consistent with earlier experimental findings, demonstrating that SC-CO₂ exhibits energy storage effects during fracturing, inducing more complex fractures. The induced fractures were not only distributed along the direction perpendicular to the minimum principal stress (vertical fractures) but also along the direction of the horizontal minimum principal stress (horizontal fracture), and the extent of vertical fracture propagation was broader.

3.4. Multi-Field Coupling Characteristic in Water and CO₂ Fracturing

(1): Breakdown Pressure

During the fracturing simulation, the variation in fluid pressure at the top of the wellbore was recorded to analyze the pressure change with dimensionless time (t/t_f), as depicted in Figure 9. It was observed that during water fracturing, the pressure initially increased linearly, followed by a deceleration in the growth rate until reaching the rock’s breakdown pressure. In contrast, during SC-CO₂ fracturing, the pressure exhibited a slow initial rise, followed by a rapid linear increase. As the pressure approached the rock’s fracturing threshold, its growth rate decelerated again until fracturing occurred. This process took longer than that water fracturing, consistent with trends observed in earlier fracturing experiments [36].

Furthermore, the theoretical breakdown pressures of the rock were determined using the previously proposed H-W formula [37] and H-F formula [35,38] and compared with the breakdown pressures obtained from the fracturing model. The results are shown in Table 2. The results indicate that the rock breakdown pressure for water fracturing closely aligned with the values calculated by the H-W formula, whereas for SC-CO₂ fracturing, it was lower than that of water fracturing, being reduced by approximately 27%, and closer to the values predicted by the H-F formula. This reaffirms that SC-CO₂ fracturing helps reduce a rock’s breakdown pressure, with the H-F formula being more suitable for modeling the SC-CO₂ fracturing process. These numerical simulation results also corroborate the experimental findings [34].

(2): Temperature Field of Formation

In this model, due to the initial temperature of the formation being 323.15 K, which is higher than the original temperatures of the water and SC-CO₂ injections, heat transfer occurred between the fluid and the rock during the fracturing process. This study recorded the temperature distributions along the green line in Figure 1 at different times to characterize the reservoir heat transfer. The distribution of reservoir temperatures at distances from the wellbore to the outer boundary is illustrated in Figure 10.

It was observed that during the injection of both fluids, there was a significant reduction in the formation temperature near the wellbore, where the reservoir temperatures closely approached the fluid temperatures. There was a temperature transition zone between the wellbore and the far field formation. In the initial stage (t = 1/10t_f), the length of the transition zone was relatively short, indicating a short distance for heat transfer. As fluid injection continued, the length of this temperature transition zone gradually increased.

In Figure 10a, it can be seen that during water fracturing, this temperature transition zone expanded from 0.01 m to 0.04 m. At t_f/10, the length of the temperature transition zone formed by water fracturing and CO₂ fracturing was similar, being 0.01 m. When the fluid injection time was t_f/4, the temperature transition zone of the water fracturing process was shorter than that of the CO₂ injection process. In the subsequent injection process, there was little difference in the length of the temperature transition zone between the two fracturing processes. The temperature transition zone of water at the completion of fracturing was close to 0.04 m, while the temperature transition zone of CO₂ fracturing was still slightly lower than the formation temperature at 0.04 m. This difference arose because both fluids had injection temperatures lower than the formation temperature, leading to convective and conductive heat transfer within the reservoir. Under the same injection rate and reservoir conditions, water, which has a higher thermal conductivity than SC-CO₂, typically results in a shorter temperature transition zone. However, in this model, water was injected at a temperature of 293.15 K, which was lower than the temperature of SC-CO₂ (308.15 K). With the injection of fluid, this rapid heat transfer effect would offset the influence of the temperature difference, and thus the transition zone would gradually become closer to SC-CO₂ fracturing. Finally, the effect caused by the temperature difference reached stability, and the temperature transition zone was slightly shorter than that of SC-CO₂ when water fracturing was completed.

(3): Formation Displacement Field

Figure 11 illustrates the contour map of subsurface displacement during the SC-CO₂ fracturing process. Initially (t = 0), the displacement of the formation increased gradually along the x direction from the left boundary and reached its maximum at the right boundary. This was due to the initial conditions where the displacements at the left and bottom boundaries were zero, and the rock bore the minimal and maximal stresses only along the other two boundaries. The compressive stress load in the y direction caused the reservoir to experience uniform displacement along the direction of minimum principal stress. As time progressed (t = 0.5t_f, 0.9t_f), the distribution of displacement in the formation gradually shifted. The lower left corner region exhibited lower displacement compared with the upper right corner due to boundary constraints, yet it still showed a progressive increase along the x- direction. However, due to differences in stress loading in two directions, the displacement field did not present a diagonally symmetrical distribution. This uneven displacement distribution is more likely to induce tensile and shear failure in rock.

To compare the displacement characteristics during water and SC-CO₂ fracturing, the wellbore top (intersection of the green line segment and red wellbore boundary in Figure 1) was selected as a reference point. The displacement over time at this location was obtained and plotted, as shown in Figure 12. For most of the duration before the rock was completely fractured, the displacement at the wellbore top gradually increased with continuous fluid injection, culminating in a sudden jump when the rock fractured completely, reaching its maximum value. Furthermore, real-time displacements during water fracturing were consistently higher than those during SC-CO₂ fracturing. Ultimately, the maximum displacement induced at the wellbore top was 0.63 mm for water fracturing and 0.45 mm for SC-CO₂ fracturing. This indicates that water fracturing induced wider fractures compared with SC-CO₂ fracturing, aligning with the experimental results on induced fracture widths.

(4): Porosity of Formation

The porosity (

φ_{p}

) is an important indicator for assessing the effectiveness of hydraulic fracturing, and the porosity of the fractured formation could be obtained through the following formula:

φ_{p} = \frac{\int_{Ω} φ_{i} d Ω}{\int_{Ω} 1 d Ω}

(29)

Here,

φ_{i}

represents the porosity values at any locations within the computational domain

Ω

. Initially, the porosity distribution within the model was integrated and subsequently divided by the area of the computational domain to determine the reservoir porosity at a specific time.

Figure 13 depicts the variation in reservoir porosity at different times during water and SC-CO₂ fracturing. It was observed that both fluid-induced reservoir porosities increased continuously with time, and the porosity increase rate caused by SC-CO₂ fracturing in the initial stage was faster. At equivalent dimensionless times, the reservoir porosity induced by SC-CO₂ was greater than that induced by water fracturing. A sudden increase in porosity occurred at the moment when the rock was completely fractured. Upon completion of fracturing, the reservoir porosities induced by water and SC-CO₂ were 0.142 and 0.1918, which increased by 1.75 times and 2.72 times compared with the original formation (0.0516), respectively. This further demonstrates that under similar conditions, SC-CO₂ fracturing achieves superior reservoir stimulation effects.

(5): Permeability of Formation

The reservoir permeability is also a crucial parameter for evaluating reservoir stimulation effects. In this study, the wellbore top was also selected as a reference point to compare the variation characteristics of reservoir permeability. Figure 14 illustrates the variation in permeability at the wellbore top with dimensionless time. It can be observed that in the initial stage of water and SC-CO₂ injection, the changes in induced reservoir permeability were relatively minor. However, until 0.9t_f for the fracturing time, the reservoir permeability underwent significant growth, culminating in a sudden increase in its maximum value upon the rock completely fracturing. The final reservoir permeabilities induced by water and SC-CO₂ fracturing were 6.34535 × 10⁻¹³ m² and 9.47968 × 10⁻¹³ m², respectively, exhibiting an exponential increase compared with the initial reservoir permeability (6.7 × 10⁻¹⁷ m²). This also indicates that under similar conditions, SC-CO₂ achieves better reservoir modification effects than water.

4. Effect of Various Factors on Coupling Characteristics during SC-CO₂ Fracturing

Based on the thermo-hydro-mechanical coupling model established in Section 2, this section will continue to investigate the influence of the in situ stress, fluid temperature, and injection rate on SC-CO₂ fracturing. The focus was on characterizing the fracture propagation characteristics under different conditions, as well as post-fracturing parameters such as the porosity and permeability of the formation. The aim was to validate the earlier fracturing experiment results and elucidate the effects of various factors on the thermo-hydro-mechanical coupling characteristics of SC-CO₂ fracturing. The model parameters were set as shown in Table 3, while other parameters not mentioned here remained consistent with those in Table 1.

4.1. Impact of In Situ Stresses on Coupling Characteristics of SC-CO₂ Fracturing

In situ stress, as one of the fundamental geological parameters, is a crucial factor influencing the effectiveness of hydraulic fracturing. This section continues the study of the fracture propagation characteristics and post-fracturing reservoir conductivity based on the thermo-hydro-mechanical coupling fracturing model. It aims to validate previous experimental results and quantitatively evaluate the stimulation effects of SC-CO₂ fracturing. Figure 15 illustrates the distribution characteristics of fractures induced by SC-CO₂ under different stress differentials. Under relatively large stress differentials, such as 4 MPa (Figure 15a), SC-CO₂ induced a primary fracture which extended through the wellbore along the direction of maximum principal stress, consistent with classical fracture mechanics theory [36]. As the stress differential decreased to 2 MPa (Figure 15b), in addition to inducing a primary fracture along the direction of maximum principal stress, secondary fractures were also initiated near the wellbore along the direction of minimum principal stress. With a further reduction in the stress differential (Figure 15c,d), the morphology of the fractures became more complex, and their propagation directions became more diverse. This included not only the primary fractures extending along the direction of maximum principal stress but also the generation of two randomly extending fracture branches, resulting in a spatial “Y”-shaped distribution of fractures. The morphological results of the fractures indicate that the in situ stress conditions are the significant factors controlling the fracture propagation patterns during SC-CO₂ fracturing. Larger stress differentials tend to induce fractures oriented more vertically to the direction of minimum principal stress in rock formations, whereas smaller differentials favor the formation of a complex network of fractures which extend in multiple directions.

Furthermore, an investigation was conducted into the porosity and permeability of the rocks post fracturing. Equation (29) was used to obtain the changes in formation porosity induced by SC-CO₂ fracturing under different stress conditions, as depicted in Figure 16. The results indicate a significant enhancement in rock porosity post SC-CO₂ fracturing compared with the initial formation porosity (0.1216). Moreover, as the stress differential decreased, the post-fracturing rock porosity gradually increased. Specifically, at a stress differential of 0, the porosity relative to the original formation increases by 1.18 times, reaching 0.265.

Additionally, using the wellbore top as a reference point, the permeability values at different times were extracted to analyze the impact of varying stress differentials on the permeability of a post-fracturing reservoir, as shown in Figure 17. It can be observed that the pattern of changes in reservoir permeability post SC-CO₂ fracturing aligned with that of the porosity. As the stress differential decreased, rock permeability gradually increased, ultimately reaching 2.14 × 10⁻¹² m², which represents an increase of 10⁴ orders of magnitude compared with the original reservoir. This further underscores that lower stress differentials enhance reservoir conductivity, thereby improving reservoir stimulation effects.

4.2. Impact of Fluid Temperature on Coupling Characteristics of SC-CO₂ Fracturing

In this section, the effects of the fluid temperature on the SC-CO₂ fracturing characteristics are explored. As shown in Table 3, the initial reservoir temperature (353.15 K) was maintained while varying the CO₂ injection temperature (313.15~353.15 K) to investigate the propagation characteristics of the fractures, porosity, and permeability parameters. Since the injection temperatures of the fluids were all above the critical temperature of CO₂, the simulation ensured that the CO₂ fluid entered a supercritical state. The fracture distribution induced by different temperatures of SC-CO₂ is illustrated in Figure 18. It was revealed that with an increasing fluid injection temperature, induced fractures transitioned gradually from multiple branching to single branching. Additionally, influenced by the maximum and minimum in situ stresses (18 and 16 MPa, respectively), fractures predominantly propagated along the direction of maximum stress. At CO₂ injection temperatures of 313.15 K and 323.15 K, two main fractures penetrating the rock and one intersecting secondary micro-fracture were induced, indicating a relatively complex fracture morphology. With an increase in fluid temperature, however, the fracture distribution pattern gradually became more singular, tending to generate a longitudinally extended fracture intersecting the wellbore, and the lower half of the fracture branch gradually shifted to the right.

These findings indicate that the injection temperature of SC-CO₂ affects the fracture propagation pattern. Under constant reservoir temperature conditions, a higher temperature differential favors inducing complex fracture patterns. This is because a greater temperature differential induces thermal stresses in the formation, which not only promotes rock breakdown but also enhances the complexity of fractures. According to the results from Section 2, during fluid injection, a temperature transition zone near the wellbore would induce thermal stresses on the rock, thereby generating more complex fractures near the wellbore. However, the extent of this transition zone is limited, and it is difficult to extend complex fracture into the far formation while solely relying on thermal stresses. As the temperature increases, although fluid diffusivity enhances, aiding in reducing the breakdown pressure of rock, the reduced temperature differential between the fluid and the formation diminishes the thermal stress effects, resulting in a relatively simpler fracture pattern.

Additionally, the porosity and permeability of rocks after SC-CO₂ fracturing at different temperatures was analyzed, as shown in Figure 19 and Figure 20. It can be observed that changes in the reservoir porosity and permeability correlated significantly with the complexity of the fracture patterns. At a temperature of 313.15 K, the porosity of the fractured rock increased nearly 1.8 times, reaching 0.14, while the permeability reached 1.156 × 10⁻¹² m², which is an exponential increase of 10⁴ times relative to the original formation. This further underscores the significant influence of temperature differentials between the formation and the fluid on reservoir stimulation. It is suggested that adjusting the fluid injection temperature can enhance fracturing effectiveness.

4.3. Impact of Injection Rate on Coupling Characteristics of SC-CO₂ Fracturing

The injection rate of fluid, as one of the fundamental engineering parameters, is also a crucial factor influencing the effectiveness of hydraulic fracturing. In this section, based on a thermo-hydro-mechanical coupling model of hydraulic fracturing, we investigated the characteristics of SC-CO₂ fracturing under different injection rates. The induced fracture distribution patterns under various SC-CO₂ injection rates are illustrated in Figure 21. The results indicate that the propagation pattern of the fractures was not obviously influenced by the injection rate. With the change in injection rate, the fracture morphology was almost the same, and it propagated along the direction of maximum principal stress. Increasing the injection rate contributed to enhancing both the width and extension distance of the fractures. These findings are consistent with field hydraulic fracturing test results.

Furthermore, an investigation was conducted into the porosity and permeability of the reservoir after hydraulic fracturing, with the results presented in Figure 22 and Figure 23, respectively. Relative to the original reservoir, significant increases in the porosity and permeability of the fractured rock were observed under different fluid injection rates. At an injection rate of 20 mL/s, the porosity of the fractured formation increased by 1.2 times, and the permeability increased by nearly 10⁴ times. With increasing injection rates of SC-CO₂, there was a trend of increasing porosity and permeability in the fractured reservoir, but the degree of growth was not significant. Compared with an injection rate of 20 mL/s, increasing the rate to 60 mL/s resulted in porosity and permeability increases of 11.5% and 7.3%, respectively, indicating that under these simulation conditions, increasing the fluid injection rate had little effect on the fracturing effect. This can be attributed to the significant influence of in situ stress on the propagation pattern of fractures in rock, which leads to a predominantly unidirectional extension of fractures. In field trials, the extension distance and width of fractures generally increased by increasing the pump rates to achieve better reservoir stimulation effects.

5. Conclusions and Suggestions

This study developed a coupled thermo-hydro-mechanical fracturing model considering the properties of CO₂ and heat transfer, based on theories of brittle fracture, phase-field, governing equations for porous media flow, and solid mechanics. It analyzed the thermo-hydro-mechanical coupling characteristics during hydraulic fracturing processes using water and SC-CO₂, and it further investigated the influence of various geological parameters (such as in situ stress) and engineering parameters (i.e., fluid injection temperature and rate) on SC-CO₂ fracturing. The previous experimental results were verified and compared, and a comparative analysis was carried out. The key findings and investigation suggestions are as follows:

(1): The simulation results of water and SC-CO₂ fracturing reaffirm the effectiveness of SC-CO₂ in reducing reservoir fracturing pressure. Moreover, SC-CO₂ fracturing induced more complex fracture patterns compared with water fracturing, enhancing the conductivity of the post-fractured reservoir. During fluid injection, heat transfer between the fluid and rock created a temperature transition zone near the wellbore. Beyond this transition zone length, the reservoir temperature remained unchanged.
(2): With the injection of SC-CO₂, the displacement field within the reservoir exhibited asymmetric shifts, and abrupt changes in displacement occurred during fracture formation, resulting in a maximum reservoir displacement smaller than that of water fracturing.
(3): Influenced by differential stress, fractures tended to propagate along directions perpendicular to the minimum horizontal stress. Greater differential stress values led to simpler crack patterns, and as these values decreased, the induced fractures from SC-CO₂ became more complex, exhibiting multiple branching and multidirectional expansion forms. Post fracturing, the porosity and permeability of the reservoir significantly increased compared with the original state. With increasing differential stress values, however, the porosity and permeability of the fractured formation showed a declining trend.
(4): The greater the temperature difference between the fluid and the formation, the more complex the induced fracture patterns induced by SC-CO₂ fracturing. However, these complex fracture branches often extended only to the near-well region, which correlated significantly with the temperature transition zone formed near the wellbore. The post-fracturing trends in reservoir porosity and permeability changes aligned with the expansion patterns of the fractures; greater temperature differences enhanced the conductivity of the reservoir after fracturing.
(5): With an increase in the injection rate of SC-CO₂, the induced fractures exhibited greater widths and extended over longer distances. Influenced by the in situ stress, the fractures primarily propagated along directions perpendicular to the minimum horizontal stress, and increasing the injection rate had minimal impact on the distribution of fracture morphology. Increasing fluid injection rates enhanced the conductivity of the reservoir, gradually increasing both the porosity and permeability post fracturing.
(6): Despite enhancing upon previous studies, this research still has limitations in both model development and the research process. For instance, it did not account for the three-dimensional propagation of fractures, the geological properties were set as fixed parameters, and it overlooked the impact of geological heterogeneity. In future studies, our research team aims to address these issues to achieve more realistic simulation results. We also plan to integrate these findings with field applications to optimize the reservoir fracturing parameters.

Author Contributions

Conceptualization, B.Y. and H.H.; methodology, H.W., Q.R. and H.C.; software, B.Y. and R.Q.; validation, L.D., Y.H., W.Z. and Y.Z.; writing—original draft preparation, B.Y. and Q.R.; writing—review and editing, Q.R.; supervision, B.Y. and H.H.; project administration, B.Y. and H.H. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by the National Natural Scientific Foundation of China (No. 52304008), the Shaanxi Province of China Science Foundation for Distinguished Young Scholars (No. 2022JC-37), the Natural Science Basic research Program of Shaanxi Province (No. 2023-JC-QN-0403), the China Postdoctoral Science Foundation (No. 2023MD734223), the Open Foundation of Shaanxi Key Laboratory of Carbon Dioxide Sequestration and Enhanced Oil Recovery, the Key Laboratory of Well Stability and Fluid & Rock Mechanics in Oil and Gas Reservoirs of Shaanxi Province (No. 23JS047), and a Xi’an Science and Technology Association young talent lifting project (959202413078).

Data Availability Statement

The detailed research parameters and results are given in this paper.

Acknowledgments

This paper was funded by the National Natural Scientific Foundation of China (No. 52304008), the Shaanxi Province of China Science Foundation for Distinguished Young Scholars (No. 2022JC-37), the Natural Science Basic research Program of Shaanxi Province (No. 2023-JC-QN-0403), the China Postdoctoral Science Foundation (No. 2023MD734223), the Open Foundation of Shaanxi Key Laboratory of Carbon Dioxide Sequestration and Enhanced Oil Recovery, the Key Laboratory of Well Stability and Fluid & Rock Mechanics in Oil and Gas Reservoirs of Shaanxi Province (No. 23JS047), and a Xi’an Science and Technology Association young talent lifting project (959202413078). We also appreciate the Youth Innovation Team of Shaanxi Universities for their assistance.

Conflicts of Interest

The authors declare no conflict of interest.

References

Soliman, M.Y.; Daal, J.; East, L. Impact of fracturing and fracturing fechniques on productivity of unconventional formations. In Proceedings of the SPE/EAGE European Unconventional Resources Conference and Exhibition, Vienna, Austria, 20–22 March 2012; p. 150949. [Google Scholar]
Longde, S.U.N.; Caineng, Z.O.U.; Ailin, J.I.A.; Yunsheng, W.; Rukai, Z.; Songtao, W.; Zhi, G. Development characteristics and orientation of tight oil and gas in China. Pet. Explor. Dev. 2019, 46, 1073–1087. [Google Scholar]
Kang, Y.; Tian, J.; Luo, P.; You, L.; Liu, X. Technical bottlenecks and development strategies of enhancing recovery for tight oil reservoirs. Acta Pet. Sin. 2020, 41, 467. [Google Scholar]
Small, M.J.; Stern, P.C.; Bomberg, E.; Christopherson, S.M.; Goldstein, B.D.; Israel, A.L.; Jackson, R.B.; Krupnick, A.; Mauter, M.S.; Zielinska, B.; et al. Risks and risk governance in unconventional Shale gas development. Environ. Sci. Technol. 2014, 48, 8289–8297. [Google Scholar] [CrossRef] [PubMed]
Slutz, J.; Anderson, J.; Broderick, R.; Horner, P. Key Shale gas water management strategies: An economic assessment tool. In Proceedings of the International Conference on Health, Safety and Environment in Oil and Gas Exploration and Production, Perth, Australia, 11–13 September 2012; p. 157532. [Google Scholar]
Wang, J.; Liu, M.; Bentley, Y.; Feng, L.; Zhang, C. Water use for shale gas extraction in the Sichuan Basin, China. J. Environ. Manag. 2018, 226, 13–21. [Google Scholar] [CrossRef] [PubMed]
Anderson, R.L.; Ratcliffe, I.; Greenwell, H.C.; Williams, P.A.; Cliffe, S.; Coveney, P.V. Clay swelling—A challenge in the oilfield. Earth-Sci. Rev. 2010, 98, 201–216. [Google Scholar] [CrossRef]
Kargbo, D.M.; Wilhelm, R.G.; Campbell, D.J. Natural gas plays in the Marcellus Shale: Challenges and potential opportunities. Environ. Sci. Technol. 2010, 44, 5679–5684. [Google Scholar] [CrossRef] [PubMed]
Yu, W.; Lashgari, H.R.; Wu, K.; Sepehrnoori, K. CO₂ injection for enhanced oil recovery in Bakken tight oil reservoirs. Fuel 2015, 159, 354–363. [Google Scholar] [CrossRef]
Abedini, A.; Torabi, F. On the CO₂ storage potential of cyclic CO₂ injection process for enhanced oil recovery. Fuel 2014, 124, 14–27. [Google Scholar] [CrossRef]
Li, X.; Elsworth, D. Geomechanics of CO₂ enhanced shale gas recovery. J. Nat. Gas Sci. Eng. 2015, 26, 1607–1619. [Google Scholar] [CrossRef]
Wang, H.; Shen, Z.; Li, G. The development and prospect of supercritical carbon dioxide drilling. Pet. Sci. Technol. 2012, 30, 1670–1676. [Google Scholar] [CrossRef]
Sinal, M.L.; Lancaster, G. Liquid CO₂ fracturing: Advantages and limitations. J. Can. Pet. Technol. 1987, 26, PETSOC-87-05-01. [Google Scholar] [CrossRef]
Duan, S.; Gu, M.; Du, X.; Xian, X. Adsorption equilibrium of CO₂ and CH₄ and their mixture on Sichuan Basin shale. Energy Fuels 2016, 30, 2248–2256. [Google Scholar] [CrossRef]
Hamza, A.; Hussein, I.A.; Al-Marri, M.J.; Mahmoud, M.; Shawabkeh, R.; Aparicio, S. CO₂ enhanced gas recovery and sequestration in depleted gas reservoirs: A review. J. Pet. Sci. Eng. 2021, 196, 107685. [Google Scholar] [CrossRef]
Li, Y.B. The challenge and opportunity of CCUS in the development of unconventional resource. Front. Energy Res. 2023, 11, 1153929. [Google Scholar]
Liu, Z.-X.; Gao, M.; Zhang, X.-M.; Liang, Y.; Guo, Y.-J.; Liu, W.-L.; Bao, J.-W. CCUS and CO₂ injection field application in abroad and China: Status and progress. Geoenergy Sci. Eng. 2023, 229, 212011. [Google Scholar] [CrossRef]
Fang, C.; Chen, W.; Amro, M. Simulation study of hydraulic fracturing using Supercritical CO₂ in Shale. In Proceedings of the Abu Dhabi International Petroleum Exhibition and Conference, Abu Dhabi, United Arab Emirates, 31 October–3 November 2014; pp. D41S–D63S. [Google Scholar]
Zhang, X.; Wang, J.; Gao, F.; Ju, Y. Impact of water, nitrogen and CO₂ fracturing fluids on fracturing initiation pressure and flow pattern in anisotropic shale reservoirs. J. Nat. Gas Sci. Eng. 2017, 45, 291–306. [Google Scholar] [CrossRef]
Liu, L.; Zhu, W.; Wei, C.; Elsworth, D.; Wang, J. Microcrack-based geomechanical modeling of rock-gas interaction during supercritical CO₂ fracturing. J. Pet. Sci. Eng. 2018, 164, 91–102. [Google Scholar] [CrossRef]
Wang, J.; Elsworth, D.; Wu, Y.; Liu, J.; Zhu, W.; Liu, Y. The Influence of fracturing fluids on fracturing processes: A comparison between water, oil and SC-CO₂. Rock Mech. Rock Eng. 2018, 51, 299–313. [Google Scholar] [CrossRef]
Zhang, Q.; Ma, D.; Liu, J.; Wang, J.; Li, X.; Zhou, Z. Numerical simulations of fracture propagation in jointed shale reservoirs under CO₂ fracturing. Geofluids 2019, 2019, 2624716. [Google Scholar] [CrossRef]
Yan, H.; Zhang, J.; Zhou, N.; Li, M. Staged numerical simulations of supercritical CO₂ fracturing of coal seams based on the extended finite element method. J. Nat. Gas Sci. Eng. 2019, 65, 275–283. [Google Scholar] [CrossRef]
Mollaali, M.; Ziaei-Rad, V.; Shen, Y. Numerical modeling of CO₂ fracturing by the phase field approach. J. Nat. Gas Sci. Eng. 2019, 70, 102905. [Google Scholar] [CrossRef]
He, Y.; Yang, Z.; Jiang, Y.; Li, X.; Zhang, Y.; Song, R. A full three-dimensional fracture propagation model for supercritical carbon dioxide fracturing. Energy Sci. Eng. 2020, 8, 2894–2906. [Google Scholar] [CrossRef]
Zhou, S.; Rabczuk, T.; Zhuang, X. Phase field modeling of quasi-static and dynamic crack propagation: COMSOL implementation and case studies. Adv. Eng. Softw. 2018, 122, 31–49. [Google Scholar] [CrossRef]
Li, H.; Xie, L.; Ren, L.; He, B.; Liu, Y.; Liu, J. Influence of CO₂-water-rock interactions on the fracture properties of sandstone from the Triassic Xujiahe Formation, Sichuan Basin. Acta Geophys. 2021, 69, 135–147. [Google Scholar] [CrossRef]
Miehe, C.; Hofacker, M.; Welschinger, F. A phase field model for rate-independent crack propagation: Robust algorithmic implementation based on operator splits. Comput. Methods Appl. Mech. Eng. 2010, 199, 2765–2778. [Google Scholar] [CrossRef]
Borden, M.J.; Verhoosel, C.V.; Scott, M.A.; Hughes, T.J.R.; Landis, C.M. A phase-field description of dynamic brittle fracture. Comput. Methods Appl. Mech. Eng. 2012, 217, 77–95. [Google Scholar] [CrossRef]
Span, R.; Wagner, W. A New Equation of State for Carbon Dioxide Covering the Fluid Region from the Triple-Point Temperature to 1100 K at Pressures up to 800 MPa. J. Phys. Chem. Ref. Data 1996, 25, 1509–1596. [Google Scholar] [CrossRef]
Li, X. Study on Fluid Heat Transfer and Rock Damage Characteristics of Carbon Dioxide Fracturing under Multi-Field Coupling. Ph.D. Thesis, China University of Petroleum, Beijing, China, 2018. [Google Scholar]
Miehe, C.; Welschinger, F.; Hofacker, M. Thermodynamically consistent phase-field models of fracture: Variational principles and multi-field FE implementations. Int. J. Numer. Methods Eng. 2010, 83, 1273–1311. [Google Scholar] [CrossRef]
Hesch, C.; Weinberg, K. Thermodynamically consistent algorithms for a finite-deformation phase-field approach to fracture. Int. J. Numer. Methods Eng. 2014, 99, 906–924. [Google Scholar] [CrossRef]
Yang, B.; Wang, H.; Wang, B.; Shen, Z.; Zheng, Y.; Jia, Z.; Yan, W. Digital quantification of fracture in full-scale rock using micro-CT images: A fracturing experiment with N₂ and CO₂. J. Pet. Sci. Eng. 2021, 196, 107682. [Google Scholar] [CrossRef]
Haimson, B.; Fairhurst, C. Initiation and extension of hydraulic fractures in rocks. Soc. Pet. Eng. J. 1967, 7, 310–318. [Google Scholar] [CrossRef]
Yang, B.; Wang, H.; Shen, Z.; Olorode, O.; Wang, B.; Zheng, Y.; Yan, W.; Jia, Z. Full-sample X-ray microcomputed tomography analysis of supercritical CO₂ fracturing in tight sandstone: Effect of stress on fracture dynamics. Energy Fuels 2021, 35, 1308–1321. [Google Scholar] [CrossRef]
Hubbert, M.K.; Willis, D.G. Mechanics of Hydraulic Fracturing. Trans. AIME 1957, 210, 153–168. [Google Scholar] [CrossRef]
Detournay, E.; Cheng, A. Influence of pressurization rate on the magnitude of the breakdown pressure. In Proceedings of the 33rd U.S. Symposium on Rock Mechanics (USRMS), Santa Fe, NM, USA, 3–5 June 1992; pp. 92–325. [Google Scholar]

Figure 1. COMSOL implementation of phase field modeling for fracturing problems.

Figure 2. Schematic of a cracked square plate under a single-edge notched tension test, where L = 1 mm.

Figure 3. Validation of 2D single-edge notched square subjected to tension. (a) When

l_{0} = 7.5 \times 10^{- 3} mm

, the crack extended in a pattern under displacements of

u = 5.3 \times 10^{- 3} mm

u = 5.6 \times 10^{- 3} mm

and

u = 5.9 \times 10^{- 3} mm

. (b) When

l_{0} = 1.5 \times 10^{- 2} mm

, the crack extended in a pattern under displacements of

u = 5.0 \times 10^{- 3} mm

u = 5.2 \times 10^{- 3} mm

, and

u = 5.4 \times 10^{- 3} mm

Figure 3. Validation of 2D single-edge notched square subjected to tension. (a) When

l_{0} = 7.5 \times 10^{- 3} mm

, the crack extended in a pattern under displacements of

u = 5.3 \times 10^{- 3} mm

u = 5.6 \times 10^{- 3} mm

and

u = 5.9 \times 10^{- 3} mm

. (b) When

l_{0} = 1.5 \times 10^{- 2} mm

, the crack extended in a pattern under displacements of

u = 5.0 \times 10^{- 3} mm

u = 5.2 \times 10^{- 3} mm

, and

u = 5.4 \times 10^{- 3} mm

Figure 4. Load–displacement curves of the 2D single-edge notched tension test.

Figure 5. Geometric model of CO₂ fracturing.

Figure 6. The main properties of CO₂: (a) density, (b) viscosity, (c) heat conductivity, and (d) specific heat capacity changing with temperature and pressure. (a) The density varied with the temperature and pressure. (b) The viscosity distribution of CO₂. (c) The thermal conductivity distribution of CO₂. (d) The specific heat capacity of CO₂.

Figure 7. Main process of water fracturing on tight reservoir.

Figure 8. Main process of CO₂ fracturing on tight reservoir.

Figure 9. The pressure at the top of the wellbore during water and SC-CO₂ fracturing.

Figure 10. Temperature distribution of formation at different times. (a) Formation temperature distribution during water fracturing. (b) Formation temperature distribution during SC-CO₂ fracturing.

Figure 11. Contour of displacement in SC-CO₂ fracturing process.

Figure 12. Displacement changes with time at top of wellbore.

Figure 13. Variation characteristics of formation porosity during water and SC-CO₂ fracturing.

Figure 14. Permeability changes with time at top of wellbore.

Figure 15. Fracture mode induced by SC-CO₂ under different stress conditions ((a) 12/8 MPa; (b) 12/10 MPa; (c) 12/11 MPa; and (d) 12/12 MPa).

Figure 16. Porosity of formation fractured by SC-CO₂ under different in situ stress values.

Figure 17. Permeability of formation fractured by SC-CO₂ under different in situ stress values.

Figure 18. Fracture mode induced by SC-CO₂ under different temperature conditions.

Figure 19. Porosity of formation fractured by SC-CO₂ under different temperatures.

Figure 20. Permeability of formation fractured by SC-CO₂ under different temperatures.

Figure 21. Fracture mode induced by SC-CO₂ under different injection rates.

Figure 22. Porosity of formation fractured by SC-CO₂ under different injection rates.

Figure 23. Permeability of formation fractured by SC-CO₂ under different injection rates.

Table 1. The parameters set in the thermo-hydro-mechanical fracturing model.

Category	Parameters	Value	Unit
Reservoir parameters	Side length (L)	0.3	m
	Initial porosity	0.0516	-
	Rock density	2460	kg/m³
	Initial permeability	6.7 × 10⁻¹⁷	m²
	Specific heat at constant pressure	1097	J/(kg·K)
	Heat conductivity coefficient	2.6	W/(m·K)
	Coefficient of thermal expansion	1.139 × 10⁻⁵	1/K
	Initial temperature	323.15	K
	Young modulus	17.8	GPa
	Poisson’s ratio	0.24	-
	Biot coefficient	0.6	-
	Maximum stress	8	MPa
	Minimum stress	5	MPa
Wellbore parameters	Radius	0.0175	m
Wellbore parameters	Initial pressure	0.1	MPa
Fluid parameters	Injection rate of water	10	mL/s
	Temperature of water	293.15	K
	Injection rate of CO₂	10	mL/s
	Injection temperature of CO₂	308.15	K
Phase field parameters	Length scale parameter ( $l_{0}$ )	0.1	mm
Phase field parameters	Critical energy release rate (G_c)	175	J/m²

Table 2. Breakdown pressure of rock obtained using different calculation methods.

Calculation Method	H-W Formula	H-F Formula	Water Fracturing Simulation Result	SC-CO₂ Fracturing Simulation Result
Breakdown pressure (MPa)	13.0	11.0	14.1	10.3

Table 3. Fracturing parameters under the influence of different factors.

Influence Factor	Parameters	Value	Unit
In situ stress	Maximum stress	12	MPa
	Minimum stress	12/11/10/8	MPa
	Initial porosity	0.1216	-
	Injection rate of fluid	40	mL/s
Fluid temperature	Fluid temperature	313.15/323.15/333.15/343.15/353.15	K
	Initial reservoir temperature	353.15	K
	Maximum stress	18	MPa
	Minimum stress	16	MPa
Injection rate	Injection rate of fluid	20/30/40/50/60	mL/s
	Maximum stress	15	MPa
	Minimum stress	10	MPa

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

© 2024 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

Share and Cite

MDPI and ACS Style

Yang, B.; Ren, Q.; Huang, H.; Wang, H.; Zheng, Y.; Dou, L.; He, Y.; Zhang, W.; Chen, H.; Qiao, R. Fracture Evolution during CO₂ Fracturing in Unconventional Formations: A Simulation Study Using the Phase Field Method. Processes 2024, 12, 1682. https://doi.org/10.3390/pr12081682

AMA Style

Yang B, Ren Q, Huang H, Wang H, Zheng Y, Dou L, He Y, Zhang W, Chen H, Qiao R. Fracture Evolution during CO₂ Fracturing in Unconventional Formations: A Simulation Study Using the Phase Field Method. Processes. 2024; 12(8):1682. https://doi.org/10.3390/pr12081682

Chicago/Turabian Style

Yang, Bing, Qianqian Ren, Hai Huang, Haizhu Wang, Yong Zheng, Liangbin Dou, Yanlong He, Wentong Zhang, Haoyu Chen, and Ruihong Qiao. 2024. "Fracture Evolution during CO₂ Fracturing in Unconventional Formations: A Simulation Study Using the Phase Field Method" Processes 12, no. 8: 1682. https://doi.org/10.3390/pr12081682

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Article metric data becomes available approximately 24 hours after publication online.

Article Menu

Fracture Evolution during CO₂ Fracturing in Unconventional Formations: A Simulation Study Using the Phase Field Method

Abstract

1. Introduction