Numerical simulation of CO 2 -enhanced oil recovery in fractured shale reservoirs using discontinuous and continuous Galerkin finite element methods

We employed a discrete fracture-matrix model to simulate CO 2 ﬂ ooding in fractured shale reservoirs, utilizing both discontinuous Galerkin (DG) and continuous Galerkin (CG) ﬁ nite element methods. The DG-CG FEM ’ s accuracy was validated against the McWhorter problem, ensuring the reliability of the simulation results. Our model also considered various factors, including reservoir heterogeneity, fracture permeability, CO 2 injection volume


Introduction
The effective exploitation of unconventional oil and gas is of great strategic significance to alleviate the contradiction between oil and gas supply and demand, promoting the low-carbon transformation of energy structure (JIN et al., 2021a).In China, shale oil resources are abundant and exhibit extensive geographical distribution.Nonetheless, notable challenges persist, notably concerning the low recovery and limited production experienced by individual wells (JIN et al., 2021b).The shale oil development goal of China is to achieve a production output of 6.5 million metric tons by 2025, with increasing shale oil recovery rate serving as a crucial support for realizing the objective (Yang and Huang, 2019).Since supercritical CO 2 can rapidly penetrate into the microporosity of reservoir rocks, CO 2 injection for enhanced recovery has been applied industrially in the field of oil and gas development (LU et al., 2021;MA et al., 2017).Not only does CO 2 injection enhance the oil and gas recovery rate, but it also facilitates CO 2 geological storage, mitigates the global greenhouse effect, and aids China in reaching its "carbon peak, carbon neutral" goals.Therefore, the development of CO 2 injection shale oil technology has a broad application prospect and strategic significance.
At present, CO 2 enhanced oil recovery from shale mainly focuses on two aspects: shale core flooding test and numerical simulation of gas injection enhanced recovery (MEI et al., 2018;JIA et al., 2019;FAN et al., 2022;Li et al., 2022;Wan et al., 2022;Huang et al., 2023).In terms of experimental research, core-scale CO 2 shale oil replacement tests have been carried out (ALHARTHY et al., 2018;ELWEGAA et al., 2019;FAKHER and IMQAM, 2020;LANG et al., 2021) and the results show that CO 2 injection can improve shale oil recovery.Moreover, the efficiency of this replacement process is influenced by factors such as the duration of CO 2 exposure, shale permeability, porosity, and the maturity of organic matter.Chen et al. (Cheng et al., 2014), Zhu et al. (ZHU et al., 2018) and Yu et al. (YU et al., 2021) adopted numerical simulations to show that volumetric fracturing penetrates the internal fractures of the reservoir, and then injecting CO 2 effectively improve the recovery rate.If CO 2 is directly injected into the reservoir with low permeability, the recovery rate of shale oil is reduced.The existing studies focus on the impact of secondary fractures generated by fracturing on gas injection and oil recovery, while ignoring the impact of primary fractures widely present inside the reservoir.
Fractures in reservoirs profoundly affect fluid flow paths and mass transport.The models describing fluid flow in fractured rock are generally categorized into equivalent continuous models, discrete fracture network models, and discrete fracture-matrix models (GLÄSER et al., 2017;BERRE et al., 2019).Among them, the equivalent continuum model is subdivided into single-pore medium model, dual-pore medium model, and multi-pore medium model (BERRE et al., 2019).The assumptions of the single-pore medium model are based on the effective medium theory, which ignores the properties of the fracture network and calculates the effective permeability of the fracture network based on the shape, size, pore size and orientation distribution of the fracture as well as the matrix permeability.The dual-pores medium model and the multi-pore medium model mathematically describe the pores and fracture, and the fluid-mass balance equations are established within the fracture and matrix systems respectively.At any point in space, hydraulic parameters (permeability, porosity, etc.) have single or multiple values.Fracture-matrix interactions are represented by fluid exchange terms that incorporate microscale effects at the macroscale and do not geometrically characterize the fracture.Therefore, it is possible to carry out relevant studies directly using the method of simulating the flow in porous media (WANG et al., 2000;LIU and ZHANG, 2008;Yang et al., 2008).Currently, the common oil and gas numerical simulation software ECLIPSE, TOUGH2 and CMG are using equivalent continuous models (PRUESS et al., 1999;GUIDE, 2002;LAW et al., 2002;LAW et al., 2003;GeoQuest, 2010).Although this method is computationally efficient, it is more difficult to geometrically realize discrete fractures in porous media and their effects on local fluid paths.In the discrete fracture network and discrete fracture matrix models, the matrix and the fracture are retained as separate geometrical objects, and the fracture region is explicitly created within the overall framework.The discrete-fracture model ignores the permeability of the matrix and considers only the fluid behavior in the fracture network.In the discrete fracture matrix model, diffusion, desorption, or two-phase seepage of the fluid within the matrix is ab considered.Meanwhile, the model downscales to deal with fractures and explicitly constructs low-dimensional fractures.The discrete fracture matrix model allows the effect of fractures on the flow topology to be explicitly modeled compared to the dual-porosity and dual-permeability model in an equivalent continuous medium (KHOEI et al., 2016a;CHEN et al., 2017;ZHANG et al., 2017;MENG et al., 2018;CUSINI et al., 2021).
Based on the discrete fracture matrix model, this paper establishes a two-phase fluid flow model for CO 2 enhanced shale oil recovery.The model incorporates the two-phase seepage governing equations for CO 2 and shale oil, accounting for their behavior in both the porous regions and low-dimensional discrete fracture within the reservoir, including the fluid-mass transfer between the fractures and the pores.We apply this established model to a randomly fractured shale reservoir for numerical simulation purposes.Through this application, we analyze the effects of reservoir inhomogeneity, fracture permeability, CO 2 injection volume, and the injection scheme on the field, as well as on the efficiency and volume of shale oil extraction.

Mathematical model
The discrete fracture matrix model utilizes explicit lowdimensional interfaces to equivalently replace the highdimensional regions occupied by fractures.The matrix and lowdimensional fractures are retained as separate geometric objects, and the discrete fracture network is explicitly modeled within the overall framework, as shown in Figure 1.The two-dimensional computational domain Ω is equivalently described as a twodimensional matrix region Ω m and a one-dimensional fracture region Ω f .Therefore, the method is also known as mixeddimensional or hybrid-dimensional discrete fracture model.The fracture computational domain contains any set of fully or partially interconnected fractures, and thus Ω f is expressed as Ω f ∪ i 1 Ω fi , i being the total number of fractures.where ∂Ω is the boundary of the computational domain, and τ and n are the tangent and normal directions of the fractures.

Governing equations for two-phase flow
The basic assumptions of the mathematical model are as follows: 1) the fluid process is isothermal; 2) CO 2 and shale oil are immiscible fluids; 3) the velocities of the free-state fluids in the matrix and the fracture satisfy Darcy's law; 4) the mass exchange of CO 2 and shale oil in the matrix and fracture satisfies the linear mass-transfer equations; and 5) the differential capillary pressure effect between the matrix and the fracture is neglected.Based on the above assumptions, the governing equations for the two-phase fluids of CO 2 and shale oil in the reservoir matrix Ω m are expressed as (ZHANG et al., 2017;MA et al., 2021a): The two-phase flow in the reservoir fracture Ω f is expressed as (LAW et al., 2003;MENG et al., 2018): where superscript α ∈ m, f denotes matrix and fractures.Subscript β ∈ o, g denotes shale oil and CO 2 .ϕ α is the porosity, S α β is the fluid saturation, is the fluid density, c β is the fluid compressibility coefficient, p α β is the fluid pressure, and p α c is the capillary pressure.Absolute permeability tensor k α k α I, k α is the absolute permeability, I is the unit matrix, k α rβ is the relative permeability, μ β is the fluid viscosity; d f is the fracture aperture.The source-sink term ] describes the jump in fluid normal flux in the matrix grids adjacent to the fracture, and the normal velocity of the fluid at the matrix grid is expressed as (Brenner et al., 2018): The weak forms of the discontinuous Galerkin method for Equations 1, 2 are shown below (MA et al., 2021a;MA et al., 2021b): The weak form of the continuous Galerkin method in Equations 3, 4 is shown below (MA et al., 2021b): where p α g and S α o are the trial functions, , The penalization factor δ o is set to 0.01, and h is the grid size.The proposed model can consider the effect of low-permeability barriers on fluids, whereas we focus on high-permeability fractures.In other words, the fractures described herein have a Frontiers in Energy Research frontiersin.orghigher permeability compared to the matrix, resulting in a relatively small gradient of pore pressure in the direction normal to the fracture.Thus, the pore pressures of matrix adjacent to the fractures is approximated equivalently as continuous when crossing the fracture.

Relative permeability and capillary pressure
The same relative permeability and capillary pressure models are selected for the matrix and fractures in shale reservoirs with the following expressions (MA et al., 2021a): In the above equation: k α ro and k α rg are the relative permeability of shale oil and CO 2 ; S α e is the effective saturation; p α c is the capillary pressure; ω α is the van Genuchten coefficient, and λ α is the pore size distribution index.

Initial and boundary conditions
The initial conditions in the matrix and fractures of the shale reservoir are set as follows: where p α gi and S α oi denote the initial CO 2 pressure and oil phase saturation of the reservoir.
The production well is set up as Dirichlet boundary conditions with constant pore pressures and saturation, which are expressed as follows: The injection well is set as a constant flow boundary: where p b g and S b o are the CO 2 pressure and oil-phase saturation at the production well, Q g is the CO 2 injection volume.
The difference in initial values of the variables between the production well and the shale reservoir increases the nonlinear characteristics of the two-phase flow equations.To further improve the convergence and stability of the model, the method of adding penalty function is introduced in this study to set Dirichlet boundary conditions (MA et al., 2021a), then equations Eq. 15 and Eq.16 are rewritten as: 3 Model implementation We use ADFNE software to generate the coordinate information of random fractures.Subsequently, we utilize COMSOL Livelink in conjunction with MATLAB to import this random fracture information into the COMSOL Multiphysics software, allowing us to generate the corresponding geometric model.Following this, we incorporate the two-phase flow equations ( 6) to (9) for both pores and fractures into the built-in PDE weak form and lowdimensional PDE weak form modules of the COMSOL Multiphysics finite element software, respectively.We also set the model boundary and initial conditions.In the time domain, we choose to discretize the equations using the q-order backward difference method.For solving the nonlinear algebraic system, we selecte the MUMPS direct solver based on LU decomposition, applying the damped Newton method.To achieve higher computational accuracy, we set the relative tolerance to 0.001.The algorithm employs an adaptive approach for time discretization, in contrast to the pre-established time step method which selects time steps randomly.By default, this method initiates with the first step being 0.1% of the total end time.

Model validation
In this section, the model and its numerical results are validated through the McWhorter problem, in which the capillary effects of immiscible and incompressible two-phase flow in porous media are considered.The geometry, boundary and initial conditions of the model are shown in Figure 2. The length of the model is 2.6 m.At the initial moment, the simulated region is fully saturated by the non-wetting phase fluid and the reservoir pressure is 2 × 10 5 Pa.The left boundary has the water saturation of 1 and the pressure is 2 × 10 5 Pa, while the other boundaries are no-flow boundaries.Brooks- Corey functions for capillary pressure and relative permeability are implemented in the model.The relevant parameter settings in the simulation are shown in Table 1.The comparisons between the semi-analytical solution of the McWhorter problem and the numerically calculated water saturation curves are given in Figure 3.The results show that the water is transported from the left boundary by about 1.5 m after 10,000 s.Meanwhile, the numerical and theoretical results match well, proving the applicability and accuracy of the numerical method in two-phase flow.

Model setup
Given the challenge of directly obtaining the distribution characteristics and attribute information of fractures within the reservoir, we employ the ADFNE program to generate two sets of two-dimensional random fractures (Alghalandis, 2017)

Variable
Value Unit The left panel illustrates the spatial distribution of water saturation at various time steps, whereas the right panel presents the line distribution of water saturation from numerical (lines) and semi-analytical (dots) solution.Geometric model, boundary conditions and initial conditions of a fractured shale reservoir.
the gas injection wells increase, while those near the extraction wells decrease, resulting in an expanding range of pressure fluctuations.The higher permeability of the fractures, compared to the surrounding matrix, facilitates a more substantial flow of CO 2 into the interior of the reservoir along the fracture regions.Pore pressure distribution at t = 50th, 300th and 600th day for CO 2 injection into fractured reservoirs.Saturation distribution at t = 50th, 300th and 600th day for CO 2 injection into fractured reservoirs.

Effects of heterogeneity in shale reservoirs
In shale oil reservoirs, a considerable number of weak structural surfaces, such as laminae and natural fractures, are developed (JIN et al., 2021a;JIN et al., 2021b).Although fractures account for a small portion of the volume of the underlying shale oil reservoir, the fractures contribute to the heterogeneity of reservoir, profoundly affecting fluid transport and shale oil extraction efficiency (LEI et al., 2021).To accurately understand the impact of shale oil reservoirs heterogeneity properties on internal fluid transport, we investigate the effects of CO 2 injection into fractured and homogeneous reservoirs.Figures 7, 8 shows the distribution of reservoir pore pressure and CO 2 saturation along the monitoring line at different time steps (t = 50th, 300th and 600th day) during CO 2 injection into both the fractured and homogeneous reservoirs.The fractures enhance the overall permeability of the shale oil reservoir, increasing the mobility of both CO 2 and shale oil.As a result, the CO 2 saturation within the fractured reservoir at the same location is significantly higher than that in the homogeneous reservoir.Due to the lower permeability of matrix in the homogeneous reservoir, CO 2 tends to accumulate near the well, leading to slightly higher CO 2 saturation in its vicinity than at the wellhead of the fractured reservoir.The pore pressure near the injection well in the homogeneous reservoir is about 14.5 MPa at 600 days into the injection; however, in the fractured reservoir, the pore pressure near the injection well is about 20.1 MPa, due to the injection of a larger volume of CO 2 .
Figure 9 illustrates the shale oil extraction rate and cumulative production resulting from CO2 injection into both fractured and homogeneous reservoirs.At the outset of gas injection and production, the production rate of shale oil surges rapidly due to the significant pressure difference between the gas extraction well and the reservoir.During the process of production, the production rate of shale oil follows a general trend of increasing initially and then decreasing.On the 50th day of gas injection, the production rate of shale oil in the fractured reservoir experiences a rebound, attributable to the injection of a larger quantity of CO 2 .By the 600th day of gas injection, the production rates of shale oil in the homogeneous and fractured reservoirs are approximately 0.09 kg/d and 0.17 kg/d, respectively, while the cumulative production amounts to approximately 60.8 kg and 90.4 kg, respectively.Throughout the entirety of the simulation process, the presence of natural fractures contributes to a 48.9% increase in the cumulative production of shale oil.

Effect of fracture permeability
Fracture permeability is typically one of the key factors influencing the efficiency of CO 2 enhanced shale oil recovery (ALFARGE et al., 2018;FENG et al., 2019).Here, we focus on analyzing the effect of fracture permeability on the behavior of fluid flow, as well as on shale oil recovery efficiency and production volume.Figures 10,11 show the distribution of CO 2 saturation and pore pressure on the 600th day under different fracture permeability conditions (k f 1.00 × 10 −14 , 1.00 × 10 −13 , 1.00 × 10 −12 m 2 ).Figures 12,13 demonstrate the distribution of CO 2 saturation and pore pressure along the monitoring line, as well as the shale oil production rate and total production volume under different fracture permeability conditions.As fracture permeability increases, CO 2 is more readily able to penetrate into the interior of the reservoir, leading to an increase in both CO 2 saturation and reservoir pressure near the gas injection wells.A higher fracture permeability facilitates the flow of shale oil, enhancing the production rate and the total volume of shale oil produced.When the fracture permeability increases from 1.00 × 10 −14 m 2 to 1.00 × 10 −13 m 2 , total shale oil production increases from 80.0 kg to 90.4 kg.However, as fracture permeability further increases from 1.00 × 10 −13 m 2 to 1.00 × 10 −12 m 2 , the total shale oil production increases slightly, from 90.4 kg to 92.6 kg, amounting to an increase of 2.4%.The evolutions of shale oil production rate and cumulative production for CO 2 injection into both fractured and homogeneous reservoirs.Distribution of pore pressure and CO 2 saturation along the monitoring line at t = 600th day with different fracture permeabilities.

FIGURE 13
Shale oil production rate and total production with different fracture permeabilities.

FIGURE 14
Distribution of pore pressure and CO 2 saturation along the monitoring line under different CO 2 injection rate conditions.
Frontiers in Energy Research frontiersin.org09 of CO 2 being injected into the reservoir at any given time.When the injection rate increases from 7.5 × 10 −6 kg/m 2 • s to 3.0 × 10 −5 kg/m 2 • s, the pressure near the injection well rises from 9.88 MPa to 14.1 MPa, and the CO 2 saturation also increases.Furthermore, a higher injection rate enhances the efficiency of shale oil displacement, leading to an increase in both the production rate and total extraction of shale oil.The production efficiencies of the three scenarios are 0.14 kg/d, 0.17 kg/d and 0.22 kg/d, while the recoveries are 81.9 kg, 90.4 kg and 107.7 kg, respectively.The simulation results indicate that a higher gas injection rate yields greater production benefits.However, it is crucial to note that higher pore pressure may increase the likelihood of damage or even rupture in the reservoir and caprock.New fractures in the caprock could provide pathways for CO 2 leakage, escalating the risk of gas escape, diminishing the CO 2 sequestration capacity of the reservoir, and potentially contaminating groundwater, among other issues.Therefore, the mechanical properties of both the reservoir and caprock warrant further investigation.

Effect of different gas injection schemes
To examine the impact of various gas injection strategies, three distinct cases are designed, as depicted in Figure 16.In case 1, the CO 2 injection rate undergoes a stepwise reduction, starting from 1.5 × 10 −5 kg and decreasing to 0.25 × 10 −5 kg.In case 2, there is a stepwise increase in the CO 2 injection rate, ranging from 0.25 × 10 −5 kg to 1.5 × 10 −5 kg.Lastly, Case 3 maintains a constant CO 2 injection rate of 8.75 × 10 −6 kg throughout the simulation.These designs ensure that the total gas injection rate remains the same for all three schemes during the entire simulation period.
The spatial distribution of CO 2 saturation at t = 600 days for different gas injection schemes is illustrated in Figure 17.After 600 days of gas injection and extraction, the CO 2 distribution within the shale reservoir is generally similar across all three injection schemes.Figures 18, 19 highlight variations in pore pressure, extraction rate, and total production near the injection and production wells under different gas injection scenarios.The pore pressures near the injection wellheads for case 1 and 3 quickly peak at 13.5 MPa and 11.5 MPa, respectively, at the onset of gas injection and extraction.Thereafter, the pore pressure decreases under the influence of the production wells.In Case 2, the pore pressure gradually increases to 10.2 MPa.Across all three injection cases, the pore pressure near the production wells follows a pattern of rapid decline followed by a rebound.In the middle stage of extraction, case 1 exhibits the highest pore pressure and extraction rate near the injection wells, while Case 2 has the lowest.The total production in Case 1 surpasses that of the other two cases, indicating that the choice of gas injection influences shale oil production, even when the same amount of CO 2 is injected.Therefore, optimizing the gas injection scheme presents tangible economic benefits for enhancing production capacity.Shale oil extraction rate and total production under different CO 2 injection rates.Distribution of CO 2 saturation at t = 600th day in different gas injection cases.
Evolution of pore pressure near injection and production wells in different gas injection cases.

FIGURE 19
Evolution of shale oil recovery rates and total production in different gas injection cases.
Frontiers in Energy Research frontiersin.org11 Ma et al. 10.3389/fenrg.2023.1330290In this paper, we conducted a comprehensive analysis of the fluid dynamics and production characteristics of shale oil in response to CO 2 injection.We presented the discrete fracture matrix model, incorporating both discontinuous and continuous Galerkin finite element methods.The main conclusions are drawn as follows: (1) Fractures in shale oil reservoirs play a crucial role in boosting both the efficiency of shale oil production and the production volumes.When compared to homogeneous reservoirs, fractured reservoirs exhibit an approximate 48.9% increase in cumulative shale oil production.Enhancing the permeability of fractures contributes to improved fluid flow capacity within the reservoir, subsequently fostering an increase in both the shale oil extraction rate and production volume.However, it is important to note that solely increasing the fracture permeability has a limited effect on augmenting the reservoir's shale oil production capacity.(2) A higher CO 2 injection rate enhances both the replacement efficiency and shale oil production.As the CO 2 injection rate increases from 7.5 × 10 −6 kg/m 2 to 3.0 × 10 −5 kg/m 2 • s, the shale oil recovery rate improves from 0.14 kg/d to 0.22 kg/d, and the recovery volume increases from 82 kg to 109 kg.Furthermore, when considering the same total volume of CO 2 injection, the chosen gas injection scheme significantly influences shale oil recovery.Simulation results indicate that gradually decreasing the CO 2 injection rate yields more favorable outcomes for shale oil production compared to a constant flow rate injection.(3) Higher CO 2 injection rates and a stepwise reduction in CO 2 injection result in higher reservoir pore pressures, increasing the possibility of damage or even destruction of the rock near the well.New fractures created could cause CO 2 leakage, which affects the overall gas injection effect and leads to environmental problems.

FIGURE 1
FIGURE 1Geometric diagram of discrete fracture model.

FIGURE 2
FIGURE 2 Geometric model, boundary conditions and initial conditions for the McWhorter problem.The subscripts w and nw refer to the wetting phase and non-wetting phase, respectively.
Figures 5, 6 represent the spatial distribution of pore pressure (p α p α g S α g + p α w S α w ) and CO 2 saturation at different time steps (t = 50th, 300th and 600th day) of CO 2 injection into the fractured reservoir for enhanced shale oil recovery, respectively.During shale oil extraction, pore pressures near FIGURE 3

FIGURE 7
FIGURE 7Distribution of pore pressure along the monitoring line at different time steps.

FIGURE 8
FIGURE 8Distribution of CO 2 saturation at different moments along the monitoring line.

Figures
Figures 14, 15  show the distribution of pore pressure and CO 2 saturation along the monitoring line, as well as the shale oil extraction rate and total amount of production, under the conditions of CO 2 injection rate (Q g 7.5 × 10 −6 , 1.5 × 10 −5 , 3.0 × 10 −5 kg/m 2 • s).Given the same parameters and grid conditions, a larger CO 2 injection rate results in a higher volume

FIGURE 10
FIGURE 10Distribution of CO 2 saturation at t = 600th day with different fracture permeabilities.

FIGURE 11
FIGURE 11Distribution of pore pressure at t = 600th day with different fracture permeabilities.

FIGURE 16
FIGURE 16Evolution of shale oil production in different gas injection cases.

TABLE 1
Parameter settings in McWhorter problems.

TABLE 2
Simulation parameters for CO 2 enhanced shale oil in a fractured reservoir.