Laboratory-to-Field Scale Numerical Investigation of Enhanced Oil Recovery Mechanism for Supercritical CO₂-Energized Fracturing

Abstract

This study systematically performs multi-scale numerical investigation of supercritical CO₂-energized fracturing, widely employed for enhanced oil recovery (EOR) in tight oil and gas reservoirs. Two distinct models, spanning from core scale to field scale, are designed to explore the diffusion patterns of CO₂ into the matrix and its impact on crude oil production at varying scales. The core-scale model employs discrete grid regions to simulate the interaction between fractures and the core, facilitating a comprehensive understanding of CO₂ diffusion and its interaction with crude oil. Based on the core-scale numerical model, the wellbore treatment process is simulated, investigating CO₂ distribution within the core and its influence on crude oil during the well treatment phase. The field-scale model employs a series of grids to simulate fractures, the matrix, and the treatment zone. Additionally, a dilation model is employed to simulate fracture initiation and closure during CO₂ fracturing and production processes. The model explores CO₂ diffusion and its interaction with crude oil at different shut-in times and various injection rates, analyzing their impact on cumulative oil production within a year. The study concludes that during shut-in, CO₂ continues to diffuse deeper into the matrix until CO₂ concentration reaches an equilibrium within a certain range. At the core scale, CO₂ penetrates approximately 4 cm into the core after a 15-day shut-in, effectively reducing the viscosity within a range of about 3.5 cm. At the field scale, CO₂ diffusion extends up to approximately 4 m, with an effective viscosity reduction zone of about 3 m. Results suggest that, theoretically, higher injection rates and longer shut-in times yield better EOR results. However, considering economic factors, a 20-day shut-in period is preferred. Different injection rates indicate varying fracture conduction capabilities upon gas injection completion.

1. Introduction

In unconventional reservoirs, due to low permeability, hydraulic fracturing is required to achieve large-scale production [1]. Horizontal drilling and hydraulic fracturing technologies have significantly increased the oil production from unconventional resources [2,3,4]. The slickwater fracturing fluid system, known for its efficiency and cost effectiveness, has found widespread application in the hydraulic fracturing process [5]. However, it has also introduced a series of issues, such as the potential for clay swelling, resulting in reduced reservoir permeability and wastage of water resources [6].

As gas is a highly compressible fluid, an energized fluid typically refers to a fluid that contains some gas content. When such energized fluids are used in hydraulic fracturing, it is referred to as energized fluid fracturing. The gases employed can include N₂, CO₂, or natural gas. When the fracturing fluid flows back, the gas will expand to provide drive energy to produce reservoir fluids—water, oil, and gas [7]. When using CO₂ as the fracturing fluid, it does not pose the same risk of formation damage as water-based fracturing fluids and exhibits high performance in enhancing hydrocarbon production [8], Meanwhile, along with an urgent demand for mitigating the greenhouse effect [9,10], it becomes increasingly attractive because of its environmentally friendly way to store the CO₂ in a depleted oil/gas reservoir [11]. The concept of Sc-CO2 sand fracturing was first proposed in 1981 and successfully applied to the Glauconite sandstone reservoir in Canada for enhanced oil and gas production through fracturing, with no damage to the reservoir and a significant increase in oil and gas production [12]. Since 2005, the domestic Changqing Oilfield, Zhongyuan Oilfield, Jilin Oilfield, and Yanchang Oilfield have successively carried out liquid CO₂ fracturing to address the problems of low reservoir pressure, severe water sensitivity, large water consumption, and slow fracturing fluid backflow, achieving good production increase results [13,14,15].

CO₂ exists in gaseous, liquid, or supercritical states under varying temperature and pressure conditions, and its phase behavior is primarily influenced by pressure and temperature. Liquid CO₂ exhibits ultra-low viscosity (0.02~0.16 mPa·s), and when the pressure exceeds 7.38 MPa and the temperature surpasses 31.1 °C, CO₂ exists in a supercritical state. In this state, CO₂ possesses high mobility, similar to a gas, while maintaining a high density close to that of a liquid [16]. In theory, the unique properties of supercritical CO₂ allow it to penetrate through any pore size larger than the size of CO₂ molecules [17]. As a result, supercritical CO₂ can infiltrate micro-fractures, including those that conventional fracturing fluids, including liquid CO₂, cannot access. During the process of fracture propagation, supercritical CO₂ experiences minimal flow resistance within the fracture, and its energy is primarily consumed by rock fracturing. This means that the fracture propagation tends to be dominated by the rock’s toughness [18]. However, the expansion of fractures in high-viscosity fluid fracturing typically occurs in a viscosity-dominated phase, where energy is primarily consumed by fluid flow. Experiments conducted on granite rock for water, liquid, and supercritical CO₂ fracturing have shown that, when injecting the lower viscosity supercritical CO₂, the rock fracturing pressure is lower, and acoustic emission monitoring signals are more likely to form a three-dimensional spatial distribution, resulting in more complex fracture geometries [19,20].

At subsurface temperatures and pressures, CO₂ dissolves in formation water to form carbonic acid. Under reservoir pressure and temperature conditions, this acidic solution has a significant corrosive effect on minerals such as carbonates and feldspar within the rock. It weakens the inter-particle bonds or etches the mineral grain lattice, thereby altering reservoir properties (such as porosity and permeability) as well as rock mechanical properties. Research on the interaction between CO₂, formation water, and rock during CO₂ sequestration and enhanced oil recovery processes (CO₂–water–rock interaction) has shown that, over extended periods, the mineral particles of the rock matrix framework and cement can significantly dissolve. This results in the formation of large pores and even micro-fractures within the matrix, leading to a significant reduction in rock mechanical strength [21]. It is, therefore, inferred that when there is water present in the reservoir (formation water or hydraulic fracturing injection water), there is likely to be CO₂–water–rock interaction during hydraulic fracturing and in the period following well shut-in. However, the extent, patterns, and primary controlling factors of this chemical interaction remain unclear. In certain stratigraphic intervals of marine shale formations in southern China, water saturation can reach 60% to 95% [22], providing the material conditions for CO₂–water–rock interaction. Furthermore, in order to leverage the advantages of CO₂ fracturing fluid and water-based fracturing fluid and enhance the complexity of fracture networks, a hybrid fracturing technique involving alternating injections of CO₂ and slickwater has emerged. This opens up the possibility of CO₂–water–rock interaction. Additionally, to ensure that CO₂ can thoroughly dissolve in the reservoir’s crude oil, reducing its viscosity, or to displace methane adsorbed in the shale, post-fracturing delayed flowback or even flowback is not necessary, providing ample time for CO₂–water–rock interaction. Typically, the influence of acidic solutions on the pattern of rock fracture propagation is manifested in the following ways [23]: On one hand, it results in the weakening of the rock’s mechanical properties, leading to changes in the characteristics of rock or fracture failure. This alteration in mechanical properties causes the transition of rocks from brittle to ductile behavior, where brittle materials are prone to tensile failure, while ductile materials are more susceptible to shear failure. Fracture friction properties change, leading to reduced strength and easier fracture extension, along with changes in the direction of the extension. On the other hand, it enhances the differential dissolution of rocks, i.e., increases heterogeneity, resulting in different fracture propagation patterns or paths.

Extensive experimental studies [24,25,26] and numerical simulation studies [27] have shown that CO₂ fracturing fluid is beneficial to fracture morphology and fracture scale. In addition to the benefits of CO₂ for fracture propagation and fracture morphology, another huge advantage of CO₂ fracturing is the CO₂–oil interactions. In unconventional fractured reservoirs, CO₂ will flow most rapidly through the major and minor fractures, but not significantly through the unfractured rock matrix due to the characteristics of low permeability and low porosity of the rock matrix [28,29,30,31]. In the early socking stage, the pressure gradient between the fracture and the matrix leads to CO₂ penetrating into the limited rock matrix, mainly known as the solution-gas drive. As the shut-in stage continues, the advective mass transfer gradually weakens, and CO₂ penetrates further into the matrix mainly through diffusive mass transfer [32]. During the whole socking stage, the CO₂ transported into the matrix dissolves into oil and causes oil swelling and viscosity reduction. In addition, the pressure slightly increases in the matrix around fractures, and this creates a local gradient where oil is extracted out of the matrix through fractures. There are few studies on the role of CO₂ diffusion in field-scale simulations. The phenomenon that the distance of CO₂ penetrating the matrix from fractures is usually several meters cannot be captured with the grids of several meters to tens of meters used in the field-scale simulation [33,34,35]. Moreover, the dual-porosity model applied generally to simulate fractured reservoirs uses two independent sets of grids to simulate fractures and the matrix, respectively, which results in the transient mass transfer process not being identified in the matrix adjacent to the fractures [36,37,38,39,40].

In this work, we have established models at both the core scale (Model 1) and field scale (Model 2) to investigate the diffusion process of CO₂ during CO₂ fracturing operations and its impact on reducing crude oil viscosity during injection and soaking processes at different scales in Section 2. The Section 3 of this study are dedicated to core scale and oilfield scale simulations, respectively. In the Section 3, we also further optimize the soaking time and injection rate.

2. Numerical Simulation

We jointly developed two numerical simulation models for simulating the mechanisms of CO₂ Enhanced Oil Recovery (EOR) at both the core scale and field scale. Both models share a common set of assumptions and mathematical formulations.

2.1. Assumption Conditions

Using the Winprop™ module in CMG (Computer Modelling Group, Calgary, AB, Canada) to create a component model for simulating CO₂ injection in the reservoir [41], while satisfying the following assumptions:

• Three-Phase Reservoir Composition: The reservoir consists of three phases, including oil, gas, and water. It contains N_c hydrocarbon components. The water phase behavior is independent of the oil and gas phases, and hydrocarbon substances exist in the oil and gas phases.
• Darcy’s Law and Consideration of Gravity Effects: fluid flow within the reservoir follows Darcy’s law, and gravity effects are considered in the flow process.
• Compressible Fluids and Rocks: fluids and rocks are both considered to be compressible in the model.
• Inter-Phase Mass Transfer: the components exhibit inter-phase mass transfer effects.
• Constant Permeability in All Directions: the permeability at a specific point within the reservoir remains constant in all directions.
• Initial Reservoir Stability: prior to CO₂ injection, the reservoir is assumed to be in a state of overall stability and is not disturbed by neighboring well operations.
• Isothermal Reservoir After CO₂ Injection: following the completion of CO₂ injection operations, the reservoir is assumed to remain in an isothermal state.

These assumptions serve as the basis for constructing the component model in the Winprop™ module of CMG and guide the simulation of CO₂ injection and its impact on the reservoir’s behavior. The model considers complex fluid flow, phase interactions, and compositional variations, providing insights into the response of the reservoir to CO₂ injection.

2.2. Mathematical Model

(1)

Continuity Equation

According to Darcy’s Law:

$${\overrightarrow{V}}_l=-{\displaystyle \frac{k{k}_{rl}}{{\mu }_l}}\left(\nabla P_l-{\rho }_{l}g\nabla D\right) (l=o,g,w)$$

(1)

In Equation (1), ${\overrightarrow{V}}_l$ represents the fluid velocity, m/s; $k$ stands for the reservoir permeability,μm²; $k_{rl}$ denotes the dimensionless relative permeability of the fluid; $P_l$ represents the fluid pressure, GPa; ${\rho }_l$ is the fluid density, kg/m³; g represents the acceleration due to gravity, m/s²; and D is the depth, with positive values denoting downward direction, m. Here, ‘l’ denotes various fluids, with ‘o’ representing crude oil, ‘g’ representing gas, and ‘w’ representing water.

Continuity Equation for Component i:

$$-\nabla .\left(C_{io}{\rho }_o{\overrightarrow{V}}_o+C_{ig}{\rho }_g{\overrightarrow{V}}_g+C_{iw}{\rho }_w{\overrightarrow{V}}_w\right)+q_i={\displaystyle \frac{\partial }{\partial t}}\left[\phi \left(C_{io}{\rho }_o{S}_o+C_{ig}{\rho }_g{S}_g+C_{iw}{\rho }_w{S}_w\right)\right]$$

(2)

In Equation (2), $C_{io}$, $C_{ig}$, and $C_{iw}$ represent the dimensionless mass fractions of component i in oil, gas, and water, respectively; ${\rho }_o$, ${\rho }_g$, and ${\rho }_w$ denote the densities of oil, gas, and water phases, kg/m³, respectively; ${\overrightarrow{V}}_o$, ${\overrightarrow{V}}_g$, and ${\overrightarrow{V}}_w$ stand for the flow velocities of oil, gas, and water, m/s, respectively; $q_i$ represents the mass flux of component i per unit volume of rock, kg/(m³·s); $\phi$ signifies the dimensionless porosity; and $S_o$, $S_g$, $S_w$, respectively represent the saturation of oil, gas, and water, all dimensionless.

Substituting Equation (1) into Equation (2) and taking into account that $C_{wo}$ and $C_{wg}$ are zero for the water component, we can derive the continuity Equations (3) and (4) for the water and hydrocarbon components:

$$\nabla \cdot \left[{\displaystyle \frac{k{k}_{rw}}{{\mu }_w}}{\rho }_w\left(\nabla P_w-{\rho }_{w}g\nabla D\right)\right]+q_w={\displaystyle \frac{\partial }{\partial t}}\left(\phi {\rho }_w{S}_w\right)$$

(3)

$$\begin{array}{l}\nabla \cdot \left[{\displaystyle \frac{k{k}_{ro}}{{\mu }_o}}{C}_{io}{\rho }_o\left(\nabla P_o-{\rho }_{o}g\nabla D\right)+{\displaystyle \frac{k{k}_{rg}}{{\mu }_g}}{C}_{ig}{\rho }_g\left(\nabla P_g-{\rho }_{g}g\nabla D\right)\right]\\ +q_i={\displaystyle \frac{\partial }{\partial t}}\left[\phi \left(C_{io}{\rho }_o{S}_o+C_{ig}{\rho }_g{S}_g\right)\right]\hspace{1em}\left(i=\mathrm{1,2},\cdots ,N_C\right)\end{array}$$

(4)

(2)

Constraint Conditions

Equilibrium State of Oil and Gas Phases:

$$f_{io}=f_{ig}\left(i=\mathrm{1,2},\cdots ,N_c\right)$$

(5)

$$S_o+S_g+S_w=1$$

(6)

$${\sum }_{i=1}^{N_C}{C}_{io}={\sum }_{i=1}^{N_C}{C}_{ig}=1$$

(7)

In Equations (5)–(7), $f_{io}$ represents the fugacity of component $i$ in oil, and $f_{ig}$ represents the fugacity of component $i$ in gas. $C_{io}$, $C_{ig}$, represent the dimensionless mass fractions of component i in oil and gas, respectively. $S_o$ represents oil phase saturation; $S_g$ represents gas phase saturation; and $S_w$ represents water phase saturation. $N_C$ represents the total number of components. The sum of their saturations equals 1, and the sum of mass fractions of each component in oil and gas phases is also equal to 1.

(3)

Auxiliary Equations

Density and viscosity of oil, gas, and water are dependent on reservoir temperature and pressure conditions as well as fluid composition:

$$\left\{\begin{array}{l}{\rho }_o=f\left(P_o,T,C_{io}\right),{\mu }_o=f\left(P_o,T,C_{io}\right)\\ {\rho }_g=f\left(P_g,T,C_{ig}\right),{\mu }_g=f\left(P_g,T,C_{ig}\right)\\ {\rho }_w=f\left(P_w,T\right),{\mu }_w=f\left(P_w,T\right)\end{array}\right.$$

(8)

(4)

Boundary Conditions

Due to the exclusion of external fluid influences on the reservoir, the outer boundaries are considered as closed boundaries:

$${\left.{\displaystyle \frac{\partial P}{\partial n}}\right|}_{\Gamma }=0$$

(9)

The inner boundary is divided into bottomhole pressure-controlled production and constant-rate production. Bottomhole pressure-controlled production complies with Equation (10), while constant-rate production adheres to Equation (11):

$${\left.P\right|}_{r=r_w}= \mathrm{Const}$$

(10)

$${\left.Q\right|}_{r=r_w}= \mathrm{Const}$$

(11)

(5)

Initial Conditions

When the reservoir is at time zero, the pressure and initial water saturation at each point remain constant:

$${\left.P(x,y,z)\right|}_{t=t_0}=P_i\hspace{1em}{\left.S_w\right|}_{t=t_0}=S_{wi}$$

(12)

2.3. Saturation-Phase Curve

Characteristics of two-phase oil and gas flow are illustrated by the phase saturation curve. The fundamental form of the exponential phase saturation curve is:

$$k_{rog}=k_{rogmax}{\left({\displaystyle \frac{S_o-S_{gcon}}{1-S_{gcon}-S_{ocon}}}\right)}^{m_o}$$

(13)

$$k_{rg}=k_{rgmax}{\left({\displaystyle \frac{1-S_o-S_{gcon}}{1-S_{gcon}-S_{ocon}}}\right)}^{m_g}$$

(14)

where $k_{rog}$ represents the oil phase relative permeability; $k_{rg}$ represents the gas phase relative permeability; $S_o$ represents oil phase saturation; $S_{gcon}$ represents gas phase critical saturation; $S_{ocon}$ represents oil phase critical saturation; $m_o$ represents the oil phase relative permeability exponent; $m_g$ represents the gas phase relative permeability exponent; and $k_{rogmax}$, $k_{rgmax}$ represent the maximum value of relative permeability.

The fundamental characteristic parameters of the phase saturation curve are presented in

Table 1. Fundamental characteristic parameters of the phase saturation curve.

Reservoir Region	kromax	krgmax	SgCOn	SoCOn	mo	mg
Region 2 in Model 1, Matrix Zone in Model 2	0.3	1	0.03	0.4	3	3
Region 1 in Model 1, Fracture Zone in Model 2	1	1	0	0	1	1

. The phase saturation curves are illustrated in Figure 1 and Figure 2.

2.4. Pressure-Sensitivity Curve

Using different pressure-sensitivity curves, we simulate stress-sensitive characteristics for supporting the main fractures, the modification zone, and the matrix. The basic form of the pressure-sensitivity curve is as follows:

$$\varphi ={\varphi }_0\cdot e^{C_f\cdot \left(p-p_{base}\right)}$$

(15)

$$k=k_0\cdot {\left({\displaystyle \frac{\varphi }{{\varphi }_0}}\right)}^m\cdot {\left({\displaystyle \frac{1-{\varphi }_0}{1-\varphi }}\right)}^2$$

(16)

where $\varphi$ represents porosity; ${\varphi }_0$ is the initial porosity; $C_f$ denotes rock compressibility, 1/MPa; $p$ stands for reservoir pressure, MPa; $p_{base}$ is the original reservoir pressure, MPa; $k$ represents permeability, md; $k_0$ is the initial permeability; and m is the permeability sensitivity exponent.

2.5. CO₂ EOR Mechanisms in Core Scale

In this section, a two-dimensional core-scale model employing a cylindrical grid configuration was established to simulate and analyze the diffusion and mass transfer of carbon dioxide at the core scale, as well as the pattern of reducing oil viscosity upon contact with carbon dioxide.

In this study, numerical models were constructed using the GEM module in CMG (Computer Modelling Group, Calgary, AB, Canada) [41]. Model 1 is delineated by two regions characterized by distinctive properties. We designate the void space of the diffusion cell, excluding the core sample, as Region ①, while the core sample itself is denoted as Region ② (Figure 3). Multiple grids have been introduced along the axial direction to effectively capture the CO₂ diffusion process. Notably, there is a singular grid designated in the circumferential axial direction, complemented by two grids set in the radial direction. Given the hypothetical assumption of symmetry in the diffusion process across both the circumferential and radial directions. Furthermore, a slender 1 mm layer is established at the interface connecting Regions 1 and ②, serving as a contact zone for oil. Consequently, the total number of grids in Model ① amounts to 2 × 1 × 121, resulting in 242 grids.

In accordance with the findings of Tsau and Barati’s research [42], the entry capillary pressure (Pe) value is established at 0.69 kPa, with the primary objective of retaining the oil phase (wetting phase) within the core cell’s central region. The pressure-sensitive characteristic parameters of the model are listed in

Table 2. Pressure-sensitive characteristic parameters of Model 1.

Region	Porosity %	Permeability mD	Compressibility Factor 10−6/kPa
Region ②	9.26	0.192	3
Region ①	90	1000	0
Thin layer	9.26	0.192	3

.

Given the intricate interplay between CO₂ and crude oil, compositional modeling is required. Based on a series of experimental studies by Qiao et al. [6], appropriate component fitting was conducted for the crude oil components to improve computational efficiency. Given the experimental data for the fluid, WinPro module has been used to get the composition model for the CMG-GEM module. The composition and properties of the simulated crude oil components are presented in

Table 3. Setups for component parameters in this numerical model.

Component	Mole Frac (%)	Critical Pressure (atm)	Critical Temperature (K)	Acentric Factor	Molecular Weight (g/mol)
CO2	-	72.8	304.2	0.225	44.01
C1	19.78552	45.4	190.6	0.008	16.043
C2–C3	5.613936	45.00811	339.36685	0.125	37.0835
C4–C6	6.770647	33.390805	480.55144	0.24959292	77.12715
C7–C12	33.58425	26.947594	596.0906	0.3915838	120.42275
C13–C34+	34.24564	14.419591	776.46775	0.84967219	335.75825

. Models 2 and 1 employ the same component model.

The tasks related to compositional modeling, as mentioned earlier, are executed through the utilization of the equation of state multiphase equilibrium property simulator Winprop™ developed by CMG. The Peng–Robinson state equation is chosen as the fundamental model, and the Jossi–Stiel–Thodos correlation model is used for viscosity characterization [43].

2.6. CO₂ EOR Mechanisms in Field Scale

In this section, we have investigated the pressure distribution, CO₂ distribution, and crude oil viscosity distribution throughout the injection stage and the well closure stage following the injection at the reservoir scale in an oilfield context.

To effectively characterize the fracture network while improving computational efficiency, Model 2 was configured with a grid size of 50 × 62 × 1. In the I direction, the grid width was set to 10 m. In the J direction, grid refinement was applied manually to simulate fractures and the modified zones, with grid widths as follows: 0.1 m for simulating fractures and 0.7932 m, 0.52348 m, and 0.29434 m for the modified zones, respectively. As shown in Figure 4, The model comprises a total of 3100 grids, and a schematic diagram of the model is presented in the figure below:

The process of simulating CO₂ pre-fracturing energy storage and production backflow was achieved by configuring different production models for distinct regions:

• Fracture Region (Region ①): In this region, the focus was on modeling the fracture network. Different amounts of CO₂ were injected to determine the fracture’s conductivity or permeability.
• Stimulated Reservoir Volume (SRV) (Region ②): This area encompassed the vicinity of the fractures. An expansion and compaction model was applied, with changes in permeability and porosity as functions of pressure variations.
• Reservoir Matrix Region (Region ③): In this region, a grid-based porous media model was established to represent the reservoir matrix. This involved considering pore-scale permeability and porosity distributions.
• The goal of this approach was to understand and simulate the processes involved in CO₂ pre-fracturing energy storage and subsequent production flowback by accounting for the different characteristics and behaviors of these three distinct regions.
• In Region 1, fractures were simulated using narrower grid spacing to capture the behavior of the fractures. A permeability multiplier curve in rock was utilized to depict the expansion of fractures occurring during the hydraulic fracturing process, the elastic closure of fractures during the shut-in phase, and the formation of proppant-filled fractures (see Figure 5).
• In Region 2, locally refined grids were utilized to simulate the modified zone formed in the vicinity of the fractures. These refined grids allowed for a detailed representation of the changes in permeability and porosity resulting from the hydraulic fracturing process.

The field-scale model incorporates input parameters associated with reservoir properties (see

Table 4. Reservoir parameters in oilfield-scale modeling.

Reservoir Parameters	Values	Fracture Parameters	Values
Initial Reservoir Pressure, MPa	13.5	Fracture Half-Length, m	150
Matrix Permeability, mD	0.25	Fracture Count	4
Reservoir Porosity	0.02	Horizontal Well Segment Length, m	50
Reservoir Thickness, m	50	Fracture Spacing, m	12.5
Fracture Zone Permeability, mD	2500	Fracture Width, m	0.1
Modified Zone Permeability, mD	100	Fracture Conductivity, D·cm	25
Reservoir temperature, °C	46.4	Initial water saturation	0.45

) and fluid properties (see

Table 5. Fluid parameters in oilfield-scale models.

Parameter	Values
Crude oil density (kg/m3) at 13.5 MPa, 46.4 °C	820.1
Crude oil viscosity (mPa·s) at 13.5 MPa, 46.4 °C	2.79
Solution gas-oil ratio (GOR) (m3/m3)	32.84
Saturation pressure (MPa)	7.39

) from the research of Qiao et al. [6]. The subject reservoir operates under conditions of diminished pressure and temperature, leading to crude oil exhibiting undersaturation, a low gas-oil ratio, and reduced viscosity. Given the notably higher solubility of CO₂ in the oil phase compared to its solubility in the water phase, this simulation disregards the interaction between formation water and CO₂.

To simulate the fracturing process, 83.3 × 10⁴ m³ of supercritical CO₂ (46.4 °C, 13 MPa) is intended to inject the fluid into Model 2 over a period of 2 h, with the aim of creating a SRV region measuring 3.2 m in width and extending within the 50 m height of the formation.

3. Results and Discussion

3.1. Results of Permeability and Viscosity Reduction Patterns of CO₂ in the Core

During the shut-in process following the gas injection, CO₂ continues to permeate into the matrix. In Model 1, we simulate the CO₂ permeation process in the core between the fracture (Region ②) and the matrix (Region ①) under a 5 MPa pressure differential after the gas injection. At the beginning of the simulation phase, the pressure in Region 2 is 14.4 MPa, while the pressure in Region 1 is 19.4 MPa. We conducted simulations and observed the distribution of CO₂ in the core and the viscosity distribution of oil at different positions within the core at various shut-in times. Specifically, we observed the distribution of CO₂ in the core after 0, 1, 7, and 15 days of shut-in (see Figure 6).

Clearly, after 1 day of shut-in simulation, CO₂ has already penetrated the core and contacted the fracture, constituting approximately 50% of the composition within the crude oil. After 7 days of shut-in simulation, CO₂ has infiltrated to approximately one-third of the core, but there is very little CO₂ content, approximately less than 10%, near the region one-third into the core. After 15 days of shut-in simulation, CO₂ has permeated into more than half of the core, and the more distant it is from the fracture-interfacing region, the lower the CO₂ content. In approximately one-third of the core length near the region interfacing with the fracture region, the CO₂ content has reached approximately 40%. To analyze the specific patterns of CO₂ infiltration into the core, a more accurate analysis is needed. Therefore, we analyzed the CO₂ content at various positions within the core at different shut-in times (see Figure 7).

From the above graph, it can be observed that after 1 day of shut-in simulation, the CO₂ content at a depth of 0.5 cm into the core at the region in interfacing with the fracture region has already reached 60%. This implies that at this point, CO₂ has displaced the crude oil deeper into the core. With an increase in the shut-in simulation time, the entire curve shifts to the right, and simultaneously, the peak of the curve gradually decreases. By 15 days of shut-in simulation, the region with CO₂ concentration exceeding 1% has filled the entire core.

Figure 8 illustrates the distribution of oil viscosity at different positions within the core after 0, 1, 7, and 15 days of shut-in simulation. It can be observed that after 1 day of shut-in simulation, the viscosity of the oil at the region in contact with the fracture (Region 2) in the core has dropped to less than 0.5 cp. By 15 days of shut-in simulation, approximately one-third of the core’s oil has a lower viscosity than before the shut-in simulation. However, the oil viscosity deeper within the core has experienced a slight increase compared to before the shut-in simulation. This is due to the pressure increase within the core resulting from the shut-in process, causing an upward migration of the oil. To study the pattern of CO₂-induced viscosity reduction in the oil, we analyzed the oil viscosity at various positions within the core at different shut-in times (see Figure 9).

It can be observed that due to the shut-in simulation, a noticeable rise in the oil’s viscosity occurs at greater depths within the core. However, within a 4 cm range from the fracture (Region 2), the oil viscosity significantly decreases, with the lowest oil viscosity reaching 0.2 cp. As the shut-in simulation progresses, the range of reduced oil viscosity gradually expands. Meanwhile, the oil viscosity near the region in contact with the fracture experiences a slight increase. This is because the permeation of CO₂ reduces the CO₂ content in that region, resulting in a weakening of the viscosity reduction effect. Additionally, shut-in simulation simultaneously causes an increase in pressure deeper within the core, leading to an increase in oil viscosity. Changes in pressure differentials can also result in varying CO₂ content within the core’s depth. To address these issues, we conducted an analysis of CO₂ diffusion and viscosity reduction under different pressure and pressure differential conditions.

Figure 10 depicts the distribution of CO₂ content and oil viscosity at various positions within the core after the execution of Scenario 3 at different shut-in times. It can be observed that after one day of shut-in, the CO₂ content at the region where the fracture contacts the core has already reached 70%. After 15 days of shut-in, the CO₂ content within 1 cm from the fracture in the core is all around 70%. This indicates that at this point, most of the oil in this region has been displaced deeper into the core. The oil’s viscosity within the core, compared to the initial viscosity before shut-in in Scenario 1, is approximately 0.4 cp higher than when there is no shut-in. After one day of shut-in, the oil viscosity at the region where the core contacts the fracture drops to around 0.2 cp, and after 15 days, the oil viscosity within 1 cm of the core’s contact with the fracture is also lowered to around 0.2 cp. This is because the majority of the oil within this area is initially displaced by CO₂ to a deeper part, and then CO₂ displaces the lighter components from the deeper oil into this region. After shut-in, as a result of the pressure elevation, the oil viscosity deeper within the core is approximately 0.3 cp higher than the results in Scenario 1. Higher pressure favors the entry of CO₂ into the core and its dissolution in the oil, thereby expanding the effective viscosity reduction range.

Figure 11 displays the distribution of CO₂ content and oil viscosity at different positions within the core after the execution of Scenario 2 at various shut-in times. Comparing the results with Scenario 1, it can be observed that under a larger pressure differential, more oil near the region where the core contacts the fracture is displaced deeper into the core. As a result, after 1 day of shut-in, the CO₂ content within 1 cm from the core’s contact with the fracture reaches approximately 65%, and after 15 days of shut-in, the deeper regions of the core have higher CO₂ content. However, the higher fracture pressure causes a greater increase in pressure within the core after shut-in, resulting in higher oil viscosity deeper within the core and a reduced effective viscosity reduction range.

The entry of CO₂ into the matrix from the fracture has a relatively limited effective range, and this range can be described using three regions. The first region (region B) is the area influenced by CO₂ diffusion, where CO₂ forms a multiphase mixture with the oil and enhances oil mobility. The second region (region C) is where CO₂ interacts with oil and efficiently diminishes the oil’s viscosity. The third region (region A) represents an area untouched by CO₂.

Drawing insights from the simulation outcomes, we established the relationship between the distance of CO₂ penetration into the matrix and time, as shown in Figure 12. Taking into account the previous analysis, we can clearly define the range of CO₂ interaction with the oil after entering the matrix. CO_2’s ability to penetrate into the core is limited, and the distance it enters into the core to enhance oil mobility (r₁) does not exceed 4.4 cm (after 15 days of shut-in). As shown in Figure 13, the effective viscosity reduction range (r₂) does not exceed 3.75 cm (after 15 days of shut-in). As the shut-in time increases, the forward movement speed of the CO₂ front gradually slows down, and the expansion rate of the CO_2’s effective range decreases.

3.2. Results of CO₂ EOR Mechanisms in Field Scale

3.2.1. During Injection Stage

In the model, the injection was carried out with a gas injection rate of 1×10⁷ m³ over a total injection time of 120 min, as depicted in Figure 14. The initial reservoir pressure was 13.5 MPa. As the injection commenced, the pressure began to rise steadily and continued to do so until the end of the 120 min injection period. Subsequently, the well designated for the injection was subsequently closed, and the production well remained shut. Within the first 60 min following the cessation of the injection, there was a gradual decline in reservoir pressure. This phenomenon can be attributed to the rapid buildup of bottomhole pressure after the cessation of the injection, causing pressure waves to propagate outward from the well, resulting in a localized pressure drop near the wellbore. The depiction of reservoir pressure distribution at the conclusion of the injection process can be observed in Figure 15.

3.2.2. During Shut-In Stage

Figure 16 presents the distribution of CO₂ concentration along the fracture direction at a distance of 1.26 m from the fracture and vertically in the direction of the fracture at a distance of 130 m from the injection point, at different shut-in times. When CO₂ injection is completed, within approximately 10 m along the fracture direction from the injection point, the CO₂ concentration reaches 0.8. At this point, crude oil within the rock matrix near this region is displaced deeper into the matrix. Within 1 m vertically in the direction of the fracture, the CO₂ concentration reaches 0.8, indicating that crude oil in this matrix region is also displaced to deeper depths. As the shut-in period progresses, the curves expand outward from the center. After 70 days of shut-in, the CO₂ concentration exceeds 0.2 at a distance of 90 m along the fracture direction from the fracture edge. In the vertical fracture direction, the leading edge with CO₂ concentrations greater than 3% is located 5.2 m away from the fracture. The distribution of CO₂ concentration at different shut-in times is shown in Figure 17.

Figure 18 illustrates the distribution of crude oil viscosity along the fracture direction at a distance of 1.26 m from the fracture and vertically in the direction of the fracture at a distance of 130 m from the injection point, at varying shut-in durations. Owing to the increased formation pressure caused by hydraulic fracturing, the viscosity exhibited by the crude oil within the rock matrix is typically greater than its initial viscosity.

Nevertheless, it remains observable that within 150 m along the fracture direction from the fracture edge, the viscosity exhibited by the crude oil within the rock matrix is significantly reduced, with the lowest viscosity dropping below 0.3 cp. Beyond 150 m from the injection point, crude oil’s viscosity suddenly increases. This is because the fracture has a half-length of 150 m, resulting in higher CO₂ concentrations within this range. Within 1 m vertically in the direction of the fracture at a distance of 130 m from the injection point, the lowest viscosity also drops to below 0.3 cp.

It is worth noting that the crude oil’s viscosity does not consistently rise with distance. An anomalous increase in viscosity occurs between the location of the lowest viscosity and the fracture, where lighter components of crude oil driven by CO₂ penetrate deeper into the matrix, leaving heavier components behind in that position. In conventional oil reservoirs undergoing CO₂ EOR processes, this phenomenon is referred to as CO₂ stripping. Prolonged shut-in periods can mitigate the increase in viscosity caused by fluid component redistribution dominated by molecular diffusion mechanisms [4]. The distribution of crude oil viscosity at various shut-in durations is shown in Figure 19.

3.2.3. Optimal Shut-In Time

By examining the distribution of CO₂ concentration and crude oil viscosity distribution during the shut-in phase, it is believed that increasing the shut-in time will facilitate the penetration of CO₂ into deeper matrix regions, thus further enhancing its Enhanced Oil Recovery mechanism. To explore how shut-in duration affects production performance, an array of simulations were conducted at various shut-in times, building upon the fundamental production model.

Figure 20 illustrates the cumulative oil production within one year at various shut-in times. The graph visually demonstrates that increasing the shut-in time can significantly enhance production. However, once the shut-in time reaches a certain point, further increases in shut-in time have a diminishing impact on production, although there is still improvement. In theory, extending the shut-in time indefinitely is the optimal choice for maximizing reservoir production. However, considering economic factors, shorter shut-in times can be a favorable option. Therefore, based on the curves presented in Figure 16, we favor a shut-in time of 20 days, as it corresponds to the point where the relative increment in production is the greatest.

3.2.4. Optimal Gas Injection Rate

During the research process, it was observed that different injection rates result in varying reservoir pressures at the completion of injection. This, in turn, affects the conductivity of fractures and leads to differences in the concentration of CO₂ penetration into the matrix. Additionally, different reservoir pressures can lead to variations in crude oil viscosity within the reservoir, which ultimately impacts production rates.

To probe into the ultimate effects of injection rates on fracture conductivity, the influence on the crude oil’s viscosity within the reservoir, and their combined impact on the final oil production, consecutive simulations were executed at various injection rates 1,3,5,7,10,13 and 1610⁶ m³/day, building upon the basic model.

Figure 21 presents the distribution of CO₂ concentration at a distance of 1.26 m from the fracture along the fracture direction and vertically in the direction of the fracture at a distance of 130 m from the injection point, 30 days after the completion of the injection, under different injection rate conditions. The graph visually demonstrates that as the injection rate escalates, the CO₂ concentration in the matrix significantly intensifies. When the injection rate reaches 10 × 10⁶ m³, the CO₂ concentration reaches 0.8 within a range of 200 m along the fracture direction from a point 1.26 m away from the fracture and within a range of 5.2 m vertically in the path of the fracture commencing from the injection location at 130 m.

Figure 22 illustrates the distribution of crude oil viscosity at a distance of 1.26 m from the fracture along the fracture direction and vertically in the direction of the fracture at a distance of 130 m from the injection point, 30 days after the completion of the injection, under different injection rate conditions. As previously analyzed, injecting larger volumes of gas leads to higher reservoir pressures, resulting in higher crude oil viscosity within the reservoir. However, it remains apparent that as the injection rate increases, crude oil viscosity noticeably decreases. When the injection rate reaches 10 × 10⁶ m³/day, the viscosity of crude oil drops below 1 cp within a range of 150 m along the fracture direction from a point 1.26 m away from the fracture and within a range of 2 m vertically along the path of the fracture starting commencing from the injection location at 130 m.

Figure 23 depicts the cumulative oil production within one year and the fracture conductivity at the completion of the injection under different injection rate conditions. Visually, it is evident that as the injection rate increases, both the cumulative oil production within one year and the fracture conductivity at the end of injection rise. Combining this observation with the earlier analysis, it is clear that higher injection rates result in greater reservoir pressure at the conclusion of the injection process, which, due to the use of the Dilation model to simulate fracture initiation, leads to increased fracture conductivity with higher injection rates. Considering economic considerations and the analysis of CO₂ concentration distribution in the matrix and the reduction in crude oil viscosity under different injection rates discussed earlier, we tend to favor an injection rate of 10 × 10⁶ m³/day.

4. Conclusions

The characteristics marked by low porosity and permeability of shale oil reservoirs lead to their elastic recovery and recovery is still very low, and there is a pressing requirement to enhance shale oil recovery. The new technology of CO₂-energized fracturing is one of the most promising and environmental-friendly measures to significantly enhance the recovery of shale oil. Through multi-scale numerical simulation of CO₂-energized fracturing in reservoirs containing shale oil, some key findings of the study can be succinctly summarized as follows:

• During the simulation of hydraulic fracturing operations, CO₂ diffusion reaches a maximum distance of approximately 4 cm into the core, effectively reducing the viscosity within a range of about 3.5 cm.
• Simulating hydraulic fracturing operations on a larger scale, CO₂ diffusion extends to a maximum distance of approximately 4 m, significantly reducing oil viscosity within a distance of about 3.5 m.
• As the injection phase ends and the shut-in phase commences, crude oil within the matrix in the immediate vicinity of the fracture is displaced deeper into the matrix. As the shut-in phase progresses, CO₂ continues to diffuse and dissolve in crude oil, ultimately leading to the re-saturation of this portion of the matrix due to the displacement effect of CO₂.
• Extending the shut-in time is advantageous for improving oil recovery. In pursuit of a balance between economic value and high production, a 20-day shut-in time is recommended.
• Different injection rates result in varying reservoir pressures at the end of the injection. Higher pressures increase oil viscosity, enhancing fracture conductivity. Taking into account economic considerations, an injection rate of 1 × 106 m³/day is considered a suitable choice.

These conclusions highlight the critical aspects of CO₂ fracturing and its impact on enhanced oil recovery, emphasizing the importance of well-designed shut-in periods, injection rates, and the associated effects on reservoir and oil behavior. The findings contribute to the understanding and optimization of CO₂-based EOR strategies in the oil and gas industry.