Effect of Reservoir Heterogeneity on CO₂ Flooding in Tight Oil Reservoirs

Abstract

Carbon dioxide (CO₂)-enhanced oil recovery (EOR) has great potential and opportunity for further development, and it is one of the vital carbon capture, utilization, and storage (CCUS) technologies. However, strong heterogeneity is one of the several challenges in developing reservoirs, especially for China’s continental tight oil reserves. This study investigates the effects of heterogeneous porosity and permeability on CO₂ flooding evolution in low-permeable tight formation. We simulated CO₂-EOR using a numerical model developed on the platform of TOUGH2MP-TMVOC to evaluate the effect of different levels of heterogeneity on oil production, gas storage, and flow behaviors in a tight reservoir, controlled by standard deviation and correlation length. A comparison of nine cases reveals that porosity heterogeneity commonly intensifies flow channeling, and there is an oil production decline with higher standard deviation and longer correlation length of porosity field. In addition, the porosity correlation length has a negligible effect on reservoir performance when the standard deviation is relatively low. Furthermore, strong heterogeneity also has a negative impact on the storage capacity of CO₂ and oil production. Notably, as the standard deviation was raised to 0.1, a small sweep region arose with the early CO₂ breakthrough, which led to a worse flooding effect. Finally, this study exemplifies that a higher injection/production rate and CO₂ alternating N₂ injection strategies can improve oil recovery in highly heterogeneous reservoirs.

1. Introduction

As an important carbon emission reduction method, CCUS technologies can help with achieving the global carbon neutrality target [1,2,3,4,5]. In particular, CCUS is a competitive approach to offset the costs by enhancing the oil recovery (CO₂-EOR), coalbed methane (CO₂-ECBM), natural gas (CO₂-EGR), and geothermal energy [2,3,6,7,8,9]. Since the 1970s, CO₂ has been utilized for tertiary recovery as commercial-scale oil production in the U.S. and China [6]. Until the end of 2021, more than 100 pilot-scale engineering CCUS projects have been established globally, and most of these are CO₂-EOR or pure CO₂ capture projects. China’s tight gas/oil development has entered a new stage over the last decade, especially in the Songliao, Sichuan, and Ordos basins. However, unlike the marine sedimentary complexes in North America, complex petroleum geology exists at various basins in China, especially for basins with continental deposits [10,11,12,13,14,15,16]. The properties of continental tight oil reservoirs are greatly influenced by deposition, resulting in poor reservoir and fluid characteristics, including severe heterogeneity, high viscosity, and low fluidity. The above factors promote significant challenges in developing continental tight oil in China [17,18,19,20,21,22]. In this regard, CO₂ has been increasingly used in gas injection and enhanced oil recovery in China due to the capability to decrease crude oil viscosity, extraction of light components, and competitive adsorption. Several CO₂-EOR pilot projects have been constructed in Songliao, Bohai Bay, and Ordos basins by the major Chinese petroleum companies [3,6,7,23].

Many parameters have been studied during the CO₂-EOR process, including the safety of CO₂ transport, CO₂ properties, and CO₂ injection strategies [24,25,26,27,28,29,30,31]; however, limited attention has been given to analyzing heterogeneity’s impact on the performance in CO₂ flooding processes. For the first time, Cheng [32] investigated heterogeneity’s effect on improved hydrocarbon recovery in the CO₂ huff-n-puff method through the UT-COMP reservoir simulator. Nevertheless, this study did not imply a sensitive relation between recovery and correlation length. Yu [33] studied the effects of permeability heterogeneity on horizontal well production and suggested that the heterogeneity is favorable for the CO₂ huff-n-puff. It is because a more heterogeneous reservoir has a higher residual oil saturation; more oil can be activated through carbon dioxide diffusion. Bao et al. [34] built 3D reservoir models with 12 hydraulic fractures to study the effects of correlation length of the natural fracture in shale reservoirs during the CO₂ huff-n-puff process. They concluded that heterogeneity is unfavorable during primary production, while a long correlation length is beneficial to oil production. Moreover, a small correlation length favors CO₂ injectivity in the short term and vice versa. Wang et al. [35] researched reservoir heterogeneity and other influencing factors using numerical simulations of CO₂ flooding in glutenite reservoirs and concluded that stronger reservoir heterogeneity causes the uneven spread of CO₂ flooding and ultimately leads to lower reservoir efficiency. In addition, the large pores and preferential seepage pathways could raise the possibility of rapid gas breakthrough. The effect of core-scale heterogeneity on the CO₂ flooding was estimated by Al-Bayati [36] via a series of man-made core samples. The experimental results from different displacement situations showed that the level of heterogeneity in layered samples influences the recovery considerably. Generally, a high permeability ratio gives rise to a lower ultimate recovery factor. Ding [37] investigated the effects of reservoir heterogeneity through indoor CO₂ flood experiments. He observed that oil production in homogeneous cores is generally more prominent than in heterogeneous cores. In addition, more CO₂ consumption occurs in heterogeneous cores to achieve the same oil production. Their work demonstrated that oil recovery is sensitive to heterogeneity for both immiscible and miscible flooding. Wu et al. [38] evaluated the effect of permeability heterogeneity on the CO₂ and N₂ injection process in tight oil reservoirs using the numerical simulation method. They concluded that reservoir performance would become poor with the strong heterogeneity. Furthermore, nitrogen injection is more favorable in maintaining formation pressure than pure carbon dioxide injection.

The above studies show that heterogeneity can lead to different effects on oil recovery and CO₂ breakthroughs; however, its effects in a low permeable reservoir are still unexplored. In addition, more oilfields are currently conducting CO₂ enhanced oil recovery pilot projects and storage in low-permeability reservoirs. Thus, the heterogeneity effect in the tight reservoir is evaluated in this work by adopting a numerical modeling approach and random porosity generation method. The simulations of nine cases with different porosity distributions have been carried out to investigate the heterogeneity influence on fluid flow and characteristics, CO₂ flooding effect, and CO₂ storage capacity in a low-permeable reservoir. In addition, this study also explored how to restrain the influence of strong heterogeneity in the process of gas injection.

2. Simulation Concept and Method

2.1. Simulation Approach

Some simulators such as TOUGH2-ECO2N, TOUGH2-ECO2M, TOUGH2-TMGAS, and MUFTE can be used to simulate the CO₂ flow process in subsurface porous media [39,40,41,42]. However, the main limitation of these frameworks is the use of two-phase conditions for modeling. Particularly, CO₂-enhanced oil recovery involves a series of physical and chemical changes of gas and oil components; thus, the mechanisms differ from the pure CO₂ sequestration in aquifer or coal beds. This work aims to evaluate the effect of heterogeneity during the CO₂-EOR process. The numerical simulator, TOUGH2MP-TMVOC, which belongs to the family of the TOUGH2 program, is used to carry out this study. TMVOC is an advanced numerical simulator that can simulate the flow behavior of water, non-condensable gases (NCGs), and non-aqueous phase liquids (NAPLs). Regarding eight gases, 18 NAPLs and water are included; thus, multi-phase multi-component non-isothermal flow of hydrocarbon mixtures in saturated/unsaturated heterogeneous media can be modeled with reasonable accuracy. In addition, the multi-phase seepage parameters in reservoirs, e.g., relative permeability, capillary pressure, and diffusion coefficient can be established and adjusted. The updating of fluid property parameters such as saturation, viscosity, and density for each phase can be executed in every iteration step.

2.2. Mathematical Model and Governing Equations

The physical processes that involve CO₂, oil, and water flow in porous media are governed by mass, momentum, and energy conservation laws. The mathematical model is based on the following assumptions and considerations: (a) transport and flow mechanisms include multi-phase Darcy flow and multi-component diffusion; (b) heat transfer in the porous medium; (c) the involvement of several components such as CO₂, water, and oil with some pseudo-oil components; (d) interphase mass transfer including dissolution of CO₂ in oil; (e) phase behavior of multi-phase and multi-component; (f) adsorption and precipitation of asphaltene. These physical laws can be represented mathematically by partial differential or integral equations at the macroscopic level.

For the multi-phase multi-component fluid flow in the porous media, the mass flux can be calculated as the sum of injection and flux from the boundaries, which can be written as:

$$\frac{\partial M^k}{\partial t}=-\overrightarrow{\nabla }\left({{\displaystyle \sum }}_{\beta }{\overrightarrow{F}}_{\beta }{x}_{\beta }^k\right)+q^k$$

(1)

where k is the index for the mass components; M^k is the mass accumulation term of component k; F is the mass or heat flux; β represents the liquid/gas/NAPL phase; $x_{\beta }^k$ is the component k’s fraction in phase β; and q is the source or sink.

The mass of a component k is the sum of its presence in all phases. Mass accumulation term M^k can be realized as the following equation:

$$M^k=\varphi {{\displaystyle \sum }}_{\beta }{S}_{\beta }{\rho }_{\beta }{x}_{\beta }^k$$

(2)

where $\varphi$ is the porosity [-]; $S_{\beta }$ is the saturation of phase β [-]; and ${\rho }_{\beta }$ is the density of phase β [mol/m³].

As for NAPL, adsorption on rock face is also considered as the following equation:

$$M^k=\varphi {{\displaystyle \sum }}_{\beta }{S}_{\beta }{\rho }_{\beta }{x}_{\beta }^k+\left(1-\phi \right){\sigma }_r{\sigma }_w{x}_w^k{k}_d$$

(3)

where ${\sigma }_r$, ${\sigma }_w$ represent the density of the rock and water, respectively; $x_w^k$ is the mole fraction of volatile oil components in aqueous phase, and $k_d$ is distribution coefficient for the water phase.

The total advective mass flux can be calculated as the sum of fluid flow of all the phases and written as:

$$F^k={{\displaystyle \sum }}_{\beta }{x}_{\beta }^k{F}_{\beta }$$

(4)

Here, the fluid mass flow can be obtained based on Darcy’s law:

$${\overrightarrow{F}}_{\beta }=-k\frac{k_{r\beta }{\rho }_{\beta }}{{\mu }_{\beta }}\left(\overrightarrow{\nabla }{P}_{\beta }-{\rho }_{\beta }\overrightarrow{g}\right)$$

(5)

where k represents the absolute permeability (m²); $k_{r\beta }$ is the phase β relative permeability (-); ${\mu }_{\beta }$is the phase β viscosity (Pa.s); $P_{\beta }$ is the phase β fluid pressure (Pa); g is the gravitational acceleration vector (m/s²);

The heat conservation includes conductive and convection, which can be expressed as:

$$F^{NK+1}=-K\nabla T+{{\displaystyle \sum }}_{\beta }{h}_{\beta }{F}_{\beta }$$

(6)

where K is thermal conductivity (K/m); $\nabla T$ denotes the temperature gradient (W/(m · K)); $h_{\beta }$ is the specific enthalpy in phase β (J/kg).

The determination of relative permeability and capillary pressure is important for the simulation of fluid flow in the porous media. Several empirical models can be used to calculate the relative permeability. Regarding three-phase fluid flow, the modified version of Stone’s three-phase relative permeability is used in this work and is described in the following equations.

$$k_{rg}={\left[\frac{S_g-S_{gr}}{1-S_{wr}}\right]}^n$$

(7)

$$k_{rw}={\left[\frac{S_w-S_{wr}}{1-S_{wr}}\right]}^n$$

(8)

$$k_{ro}=\left[\frac{1-S_g-S_w-S_{or}}{1-S_g-S_{wr}-S_{or}}\right]\left[\frac{1-S_{wr}-S_{or}}{1-S_w-S_{or}}\right]{\left[\frac{\left(1-S_g-S_{wr}-S_{or}\right)\left(1-S_w\right)}{\left(1-S_{wr}\right)}\right]}^n$$

(9)

where $k_{rg}$ is the gas phase relative permeability; $k_{rw}$ is the water phase relative permeability; $k_{ro}$ is the NAPL phase relative permeability; $S_g$ is the gas phase saturation; $S_w$ is the water phase saturation; $S_o$ is the NAPL saturation ($S_o=1-S_g-S_w$); $S_{gr}$ is the irreducible gas saturation; $S_{wr}$ is the irreducible water saturation; and $S_{or}$ is the irreducible NAPL saturation.

Similarly, the three-phase capillary pressure function of Parker et al. [43], which is also a built-in function in TMVOC, is utilized to calculate the capillary pressure between phases. It is given in the following equation:

$$P_{cgo}=-\frac{{\rho }_{w}g}{{\alpha }_{go}}{\left[{\overline{S}}_l{}^{-1∕n}-1\right]}^{1∕m}$$

(10)

$$P_{cgw}=-\frac{{\rho }_{w}g}{{\alpha }_{ow}}{\left[{\left({\overline{S}}_W\right)}^{-1/m}-1\right]}^{1/n}-\frac{{\rho }_{w}g}{{\alpha }_{go}}{\left[{\left({\overline{S}}_l\right)}^{-1∕m}-1\right]}^{1∕n}$$

(11)

$$P_{cow}=P_{cgw}-P_{cgo}$$

(12)

where ${\overline{S}}_w=\frac{\left(S_w-S_m\right)}{\left(1-S_m\right)},$ ${\overline{S}}_l=\frac{\left(S_w+S_n-S_m\right)}{\left(1-S_m\right)}$, $m=1-\frac{1}{n}$, $P_{cga}$ is the gas–NAPL capillary pressure, and$P_{cgw}$ is the gas–water capillary pressure.

Molecular diffusion also plays a vital role in mass transport, especially for the small advective velocity. Typically, diffusion flux is related to the concentration gradient of the diffusing component and can be expressed according to Fick’s law as follows:

$$f=-d\nabla C$$

(13)

where d denotes the effective diffusivity. This parameter is dependent on the porous media, the pore fluid, and the diffusing component normally. The concentration variable C can be chosen by mass and mole fraction.

The concentration can be expressed by the various models, such as De Marsily’s (1986) function [44]. Consider the multi-components diffusion in

$$f_{\beta }^k=-\varphi {\tau }_0{\tau }_{\beta }{\rho }_{\beta }{d}_{\beta }^k\nabla X_{\beta }^k$$

(14)

where $\varphi$ denotes the porosity [-]; ${\tau }_0$ is a factor dependent on the porous media [-]; ${\tau }_{\beta }$ represents the coefficient that depends on the saturation of phase β [-]; ${\tau }_0{\tau }_{\beta }$ denotes the tortuosity [-]; ${\rho }_{\beta }$ is the density [kg/m³]; and $d_{\beta }^k$ represents the molecular diffusion coefficient for component k in phase β [kg/m³].

Therefore, a single effective multi-phase diffusion coefficient that combines all material constants and tortuosity factors can be expressed as:

$${{\displaystyle \sum }}_{\beta }^k=\varphi {\tau }_0{\tau }_{\beta }{\rho }_{\beta }{d}_{\beta }^k$$

(15)

For the two-phase conditions, the diffusive flux can be written as:

$$f^k=-{{\displaystyle \sum }}_1^k\nabla X_1^k-{{\displaystyle \sum }}_g^k\nabla X_g^k$$

(16)

It is worth mentioning that in TMVOC, NCG dissolution in the aqueous phase and NAPL phase is described by Henry’s law. However, the exact solubility of CO₂ in oil depends on oil components, pressure, and temperature for the unconventional hydrocarbon reservoirs. In addition, due to the difficulty of determining the equilibrium solubility by correlating at any temperature and pressure accurately, the experimental data of CO₂ solubility in a specific crude oil can be used in the numerical model. In this work, the experimental solubility is adopted from the work of Ju [45] to realize the gas flooding process under the condition of immiscible flooding.

In TOUGH2MP-TMVOC, the integral finite difference method has been used to discretize the continuum equations for solving multi-flow in the porous medium. Therefore, the integral of Equation (1) can be written in the following equation.

$$\frac{d}{dt}{{\displaystyle \int }}_{V_n}{M}^{k}d{V}_n=-{{\displaystyle \int }}_{{\mathsf{\Gamma }}_n}{F}^k·n d{\mathsf{\Gamma }}_n+{{\displaystyle \int }}_{V_n}{q}^{k}d{V}_n$$

(17)

2.3. Numerical Model Validation

The accuracy of the simulation method and model used in this work have been verified in the previous research [9,46,47,48], which can guarantee the reliability of numerical simulations. However, to ensure the applicability of the simulation method for CO₂-EOR, model validation is also carried out. One 1D model was generated for comparison with the slim-tube experimental data conducted by Jilin Oilfeld, CNPC [49].

Table 1. Main parameters data of 1D model.

Parameter	Value	Parameter	Value
Grid block length (m)	10 × 1 × 1	Water saturation	0.2
Grid numbers	10	Oil saturation	0.8
Porosity	0.3	CO2 injection rate (kg/s)	0.0001
Initial temperature (°C)	90	Initial pressure (MPa)	30

summarizes the main parameters utilized in the simulation that are the same as the experiments. The mesh model with injection and production locations is demonstrated in Figure 1.

Figure 2a presents the relation between the oil recovery and CO₂ injection volume. It can be observed that when the injected CO₂ volume reaches 0.4 times the pore volume (0.4 PV), the recovery obtained from the simulation started to be lower than the experimental value until the injected volume reaches about 0.82 PV. Furthermore, due to the gas breakthrough, it also can be found that the enhancement in recovery after 1.0 PV injection of CO₂ slowed down. In addition, the simulated recovery was lower than the experimental results after 1.0 PV injection, which is probably since the miscible behavior occurred in the experiment. Except for a few differences between the simulated and experimental results, the recovery matches well, which validates the reliability of our numerical model.

The study of discretization and grid convergence is performed to demonstrate the numerical method for the CO₂-EOR. In this section, six different grid block sizes (a different number of cells) are generated, with 5, 10, 20, 40, 80, and 160 cells, respectively. All the parameters for the six simulations are the same as those of the model used for the comparison with experimental data above. It can be observed that the recovery degree gradually improves with the increase in the number of cells (Figure 2b). The maximum final recovery value is 96.68% with 160 cells, which is close to the recovery of 96.31% obtained with 80 cells. Similarly, the results obtained by 40 cells, 20 cells, and 10 cells are 95.56%, 94.62%, and 93.68%, respectively. It is worth mentioning that the value of recovery does not deviate from each other significantly in the five cases. However, the final oil recovery is maintained at 89.04% when the total cells are only 5; thus, increasing the grid block size (more than 1 × 1 × 1 m) will decrease the final recovery significantly. This study demonstrated that when the grid size is less than 1 × 1 × 1 m, the numerical model has well grid convergence and consistency. It also proved the capability of the method in handling CO₂-EOR simulation. The results present that there is not much difference when the cells increase from 10 to 160 and regarding computational time, a grid size of 1 × 1 m in the x, y-axis is employed in the next simulation work.

3. Model Generation

3.1. Heterogeneous Porosity and Permeability Fields

Typically, the formation porosity follows a normal or log-normal distribution, which can be represented by a spatially autocorrelated random field. The variogram model and the correlation length can characterize the spatial autocorrelation features. The variance model describes how the semivariance (the difference between the statistical variance of porosity and the covariance) changes with the distance between two locations. The exponential and spherical variogram models have been broadly used for geostatistics. The spherical variance function model is used in this study to realize the two-dimensional log-normal porosity distributions. In order to obtain the porosity fields that follow the specified probability distribution and spatial autocorrelation for the simulations, the “gstat” package in R programming language [50,51] was employed for the heterogeneous modeling.

In this study, the degree of heterogeneity is evaluated by the standard deviation (SD) and the correlation length (CL) of porosity values, since the spatial variation in porosity determines flow pattern and reservoir permeability. After obtaining the random porosity distribution, the distribution of permeability is simulated according to the KC equation. Permeability is related to porosity, tortuosity, and other parameters via this semi-empirical equation. In addition, the corresponding modified KC equation is used for the low permeability reservoirs, which allows obtaining more accurate permeability. Generally, it has the following form [52,53]:

$${\mathrm{K}}_{\mathrm{K}\mathrm{C}}=\frac{8}{35\left(3+2\sqrt{2}\right)}\frac{{\varphi }^3}{k{\left(1-\varphi \right)}^2}{\lambda }^2$$

(18)

where d is the grain diameter, and k = cτ² is the K-C constant, which is dependent on the material.

Table 2. Simulation plan of the porosity field.

	Mean	Standard Deviation (-)	Correlation Length (m)
Set 1:	0.07	0.025	12.5
	0.07	0.025	25
	0.07	0.025	50
	0.07	0.05	12.5
Set 2:	0.07	0.05	25
	0.07	0.05	50
	0.07	0.1	12.5
Set 3:	0.07	0.1	25
	0.07	0.1	50

summarizes the plan for this investigation. The simulation consists of three sets with different values of SD and CL, while the porosity with a mean value of 0.07 in each case.

3.2. Numerical Model

A brick model of 160 × 100 × 10 m in the x, y, and z directions is used for simulation, forming part of a common tight reservoir well pattern with a confined boundary modified from Wang et al. [54]. The computation domain is illustrated in the red rectangle with the injection well on the upper left corner and the production well on the lower right (Figure 3). The half-length of one fracture is 80 m with the porosity value of 0.5, while the average value of the matrix porosity is 0.07. The initial distribution of components and properties of the oil phase are shown in

Table 3. Oil phase composition and properties.

Component	Mole Fraction, %	Critical Temperature, K	Critical Pressure, MPa	Boiling Point, K	Molecular Weight, g/mol	Density (289 K), kg/m3
COMP1	0.36	630.3	3.73	417.6	106.2	880.0
COMP2	0.64	591.8	4.1	383.8	92.1	867.1

[55], and the input parameters applied in the model for simulations are listed in

Table 4. Basic parameters data of simulation model.

Property	Value	Unit
Grid block length	100 × 160 × 10	m
Mean matrix porosity	0.07	dimensionless
Fracture length	80	m
Fracture porosity	0.5	dimensionless
Initial temperature	90	°C
Injection temperature	40	°C
Initial pressure	20	Mpa
Injection rate	0.01	kg/s
Production pressure	15	Mpa
Water saturation	0.4	dimensionless
Oil saturation	0.6	dimensionless
Production time	10	year

[56].

The porosity and permeability distribution fields in the simulated model are generated by the interpolation method of Section 3.1 (Figure 4). It can be observed that if the SD and CL values are small, the porosity distribution is more uniform, and the larger variation appears in a small range. As the SD value increases, the homogeneity of porosity distribution declines. Moreover, as the CL value continues to increase, the porosity distribution in the same set is no longer uniform, and the areas with larger values appear as a cluster. Along with the CL enhancement, the area of the clustered aggregation increases (Figure 4i). The permeability distribution almost corresponds to the porosity and maintains the same distribution pattern, as shown in Figure 5.

4. Simulation Results and Discussion

The homogeneous fields have also been simulated for comparison. The porosity value of the homogeneous field is 0.07. The gas, oil saturation, and reservoir pressure during the CO₂ injection process have been measured. Figure 6 shows the evolution of the reservoir phase saturation and pressure in the half, second, sixth, and tenth years. The CO₂ displacement front advances evenly during the injection process, and the gas and oil saturations alter correspondingly. However, the pressure field changes significantly in the early injection stage, and the pressure near the injection well increases. The pressure near the injection well declines with continuous production, and the overall distribution remains uniform. Moreover, the final cumulative oil recovery of the homogeneous model is 22.23%.

4.1. Effect on Oil/Gas Saturation and Pressure

Figure 7 describes the variations in oil/gas saturation and reservoir pressure using Set 1 during the CO₂ flooding process. The results of three cases for short, medium, and long correlation lengths (i.e., 12.5 m, 25 m, and 50 m) are presented. It can be seen that with the enhancement in the heterogeneity of the reservoir, uneven distribution appears in the CO₂ displacement front, which is clearly different from the homogeneous case. When SD is equal to 0.025, the final sweep area of each case is roughly the same. Nevertheless, when CL is 50 m, the sweep region is slightly reduced, mainly in the lower-left and upper-right corners. With the continuous injection of CO₂, the gas saturation near the injection well increases and reaches 0.5 by the end of 10-year production. The reservoir pressure near the injection well increases remarkably and reaches about 30 MPa in a 6-month period. After two years of production, the overall reservoir pressure has stabilized, and the pressure near the production well has dropped markedly.

Figure 8 demonstrates the changes of oil and gas saturation and reservoir pressure in the case with an SD value of 0.05. The heterogeneity in porosity becomes large, making the difference in the range of saturation variation more obvious. In other words, the direction of CO₂ flow and morphology are significantly influenced by the heterogeneity in porosity and permeability, and the distinction between sweep regions is large. When CL is 12.5 m, the porosity and permeability near the injection well are relatively low, and the CO₂ front proceeds uniformly. In addition, the displacement effect in the right side of the reservoir is better, leading to the largest sweep area. In the cases where CL equals 25 m and 50 m, the CO₂ advances mainly to the downward and right sides, and the formation pressure increases accordingly. The maximum formation pressure is close to 35 MPa in these two cases during the entire process. While the case with CL equals 12.5 m, the maximum formation pressure reaches 42 MPa due to the poor seepage conditions near the injection well. Afterward, the formation pressure decreases uniformly as the CO₂ enters the region with better seepage conditions having smaller correlation lengths and better homogeneity.

The oil/gas saturation and the pressure change of Set 3 are shown in Figure 9. In Set 3, the standard deviation of porosity reaches 0.1, which means the heterogeneity is the strongest. The saturation distribution maps show that the swept area of CO₂ is tremendously reduced, and the flow pattern of carbon dioxide is also quite different from the first two sets at different times. In other words, the overall displacement efficiency is significantly reduced. It can be seen that when the CO₂ flooding is carried out for six months, the change in gas saturation range is significantly smaller compared to the first two sets. For the cases where CL equals 12.5 m and 25 m, CO₂ mainly flows along the left side. However, in the case with CL equal to 50 m, CO₂ mainly spreads horizontally along the upper area. In the case of CL equal to 50 m, CO₂ has the smallest swept area by the second year, and CO₂ has already reached the location of the production well. However, it reaches the production well in approximately four years in the other two cases, indicating an earlier CO₂ breakthrough in the strongest heterogeneity fields.

4.2. Effects on CO₂ Storage and Oil Production

This section compares net CO₂ storage and oil production capacity at different heterogeneity levels for all cases. Figure 10 shows the results obtained from these ten heterogeneous cases. The total amount of net CO₂ stored and oil produced varies with standard deviation and correlation length. For the homogenous field, the storage capacity of CO₂ is large and has the same value when SD is 0.025. In the cases where SD equals 0.025, the storage capacity decreases when CL increases from 12.5 to 50 m. It is observed that there is no significant influence of CL at the SD of 0.025 on the storage capacity. When SD increases to 0.05 and 0.1, a noticeable difference in CO₂ storage capacity appears. In the case where SD equals 0.05 and CL equals 12.5 m, the lowest CO₂ storage amount is found. It is because the porosity and permeability near the injection well are quite low, and the seepage condition is poor, which leads to complexity in CO₂ injection. Instead, it has the maximum CO₂ storage amount when the CL is 25 m in Set 2. Moreover, in the case where SD is 0.1, the storage capacity gradually decreases as the CL increases from 12.5 m to 50 m. In the case of oil production potential, it is also observed that the homogeneous field contributes to the largest production, and there is a slight decline in oil production as SD increases. It is evident from Figure 10 that in Set 1 and Set 2, the standard deviations were 0.025 and 0.05, respectively, and the increase in CL resulted in a decrease in oil production. Comparatively, in Set 3 with a standard deviation of 0.1, enhancing CL did not decline production. It suggests that when the standard deviation increases to a certain level, changing the correlation length has little effect on oil production; however, it still has a large impact on net CO₂ storage.

In order to quantify the reservoir performance and the effects of heterogeneity on oil production, we calculated and compared the results of cumulative recovery for all cases. Figure 11 summarizes the evolutions of the cumulative recovery curve. As a reference for comparison, the final cumulative recovery of the homogeneous case is 22.23%. However, the average cumulative recovery in Set 1, Set 2, and Set 3 is 22.09%, 20.28%, and 16.94%, respectively, indicating that the enhanced heterogeneity due to the increased standard deviation of porosity leads to a decrease in the degree of recovery. It should be noted that the highest cumulative recovery, 22.77%, is achieved when the SD is 0.025, and the CL is 12.5 m. In the case where SD equals 0.1 and the associated length equals 50 m, the recovery degree is 16%, which is the lowest value. However, in this case, the cumulative recovery remains the highest for the first two years, and then, the growth rate declines considerably. Moreover, the cumulative recovery showed a downward trend with the increase in correlation length for each set.

4.3. Effects on CO₂ and Oil Flow

This section investigates the impact of heterogeneity on the flow morphology of CO₂ and oil in the reservoir region by calculating and plotting the flow velocity distribution at different time durations. Three sets of oil and gas flow characteristics with different degrees of heterogeneity are shown in Figure 12, Figure 13 and Figure 14. It can be seen that the flow of oil and gas is greatly affected by the degree of heterogeneity. In three cases of Set 1, the maximum flow velocity of CO₂ and oil are located at the injection/production well (Figure 12). The flow velocity and area for each case differ slightly, and the whole distribution is relatively uniform and regular. Figure 13 presents the oil and gas flow characteristics of Set 2, in which the increase in the porosity standard deviation leads to the appearance of unevenly distributed flow. The flow of CO₂ along a particular streamline starts to appear in the three cases. In addition, the flow of oil changes into dispersion, and the streamlines become sparse in the subsequent stages. This phenomenon becomes more prominent when SD increases to 0.1 (Figure 14). The density of the streamline distribution reduces significantly, and the flow of oil and gas is mainly concentrated in some areas. The main reason is that CO₂ connects the area with better pore seepage conditions between the injection and production wells, forming a predominant pathway. In Set 3, as the CL value increases, the effective flow range of oil and gas gradually decreases, especially for CO₂, as shown in Figure 14c, and the streamline forms a relatively narrow channel. Since most of the oil is not swept effectively, the seepage range is limited to the vicinity of the dominant channel. Moreover, due to the continuous injection of CO₂, the flow is even more channelized along these preferential paths. It can be seen that CO₂ bypassed certain areas and could not reach certain areas with poor porosity at all because of these predominant pathways (Figure 14a,b). In addition, with the increase in the degree of heterogeneity, the flow velocity magnitude is also affected. In Set 3, when SD equals 0.1, the maximum flow rate of CO₂ is between 0.2 and 0.32 kg/s, while the maximum flow rate ranges of the remaining two cases are 0.15–0.20 kg/s and 0.2–0.24 kg/s, respectively.

In order to encounter the issue of CO₂ breakthrough directly, the flow velocity of each phase at the production well is monitored (Figure 15). It can be seen that the flow velocity of both CO₂ and oil increases gradually in the early stage. Once the breakthrough of CO₂ occurs, the flow velocity of CO₂ becomes large in a short period. At the same time, the flow velocity of oil decreases rapidly. Therefore, oil flow velocity in each case has experienced a peak (Figure 15a). The peak appears around the third year in Set 1 but slightly earlier in the second and third set as the heterogeneity increases. It is noteworthy that in Set 2 and Set 3, with stronger heterogeneity, the oil flow velocity experiences a rapid increase for a short period before reaching its peak. It is probably related to the fact that oil production increases for a short period when CO₂ displaces into the previously un-swept regions. This phenomenon is more obvious in the reservoir with stronger heterogeneity since the region with better porosity and permeability becomes more heterogeneous. It results in a significantly higher velocity of CO₂ and oil compared to other cases in the early period as SD equal to 0.1 and CL is 50 m. In this case, the oil flow velocity is minimal in the later period. The rapid increase in CO₂ flow velocity and the peak of the oil flow velocity are both earlier than the other cases, indicating that the CO₂ breakthrough has occurred at the earliest.

4.4. Effect of CO₂ Injection Rate and Production Pressure

Two cases with the strongest and the weakest heterogeneity have been selected, and a comparative study of overcoming the impact of heterogeneity during production is carried out. The production enhancement of these measurements in the development of reservoirs with different degrees of heterogeneity is explored by analyzing the recovery under different CO₂ injection rates and production pressures. The enhancement in gas injection rate results in a significant improvement in recovery (Figure 16a). However, the final cumulative recovery increases about 16.43% in the weak heterogeneity case compared to 8.58% in the stronger heterogeneity case. Nevertheless, reducing the production pressure in the simulation can also increase the recovery effectively, but no noticeable difference is found related to the heterogeneity. The reduction in production pressure from 20 to 7.5 MPa results in the increase in recovery by 4.34% and 5.38% for the two cases, respectively, and with the weaker heterogeneity case improving slightly less (Figure 16b). The above investigation shows that increasing the injection and production rate will increase the oil production and contribute to positive effects.

4.5. Effect of CO₂-Alternating-N₂ (CAN) Injection

Nitrogen gas has good elastic energy during flooding to maintain the reservoir pressure. It has many advantages, such as being abundant in the atmosphere, relatively low-cost, well-established harnessing/capture technology, and non-corrosive property [57,58,59,60]. These advantages make nitrogen gas another significant candidate to improve tight oil reservoirs. Several studies have been carried out to investigate the efficiency of N₂ gas injection on tight oil reservoirs, which suggest that the mixed injection of nitrogen and carbon dioxide can achieve better results than pure carbon dioxide injection [38,60].

In order to reduce the impact of strong heterogeneity on CO₂ flooding, further simulations have been conducted to investigate the mixed injection of carbon dioxide and nitrogen (CO₂-alternating-N₂) and compare its optimization effect in the most heterogeneous reservoirs. Carbon dioxide and nitrogen are injected alternately into the reservoir with an injection rate of 0.01 kg/s. Firstly, carbon dioxide is injected in the first year, which is followed by a continuous injection of nitrogen in the next year. A two-year injection cycle was generated, and the entire simulation process was carried out for five cycles (Figure 17). The cumulative recovery is also used to measure the oil-flooding effect. As shown in Figure 18, in the case of an alternate injection of carbon dioxide and nitrogen, the recovery degree is higher from the early stage than pure carbon dioxide flooding, and the final recovery degree reaches 19.94%. However, the final recovery degree reaches 16.04% in the pure carbon dioxide case. The results indicate that the alternate injection of CO₂ and N₂ can obtain a better oil recovery. It is due to the lower solubility of nitrogen in oil, which can maintain the gas saturation and the reservoir pressure. In addition, a large amount of injected gas can penetrate the low permeability area due to increased mobility, leading to a large sweep area.

5. Conclusions

We used the advanced TOUGH2MP-TMVOC numerical simulator to study the flow process in engineered CO₂ enhanced oil recovery. The effects of spatial porosity heterogeneity on sweep region, flow pattern, oil production, and carbon storage capacity in the tight reservoir are studied. We investigated the two primary variables, the standard deviation and correlation length of the porosity field, and dozens of realizations were performed to analyze the reservoir performance.

The simulation results show that porosity and permeability play a significant role in flow pattern evolution and oil production in the low-permeability reservoir. Compared with homogeneous porosity and permeability fields, heterogeneous fields that ubiquitously exist in nature tend to exacerbate flow channeling and generally undermine reservoir performance. Flow channeling is inevitable regardless of the initial flow pattern. Higher standard deviation and longer correlation length generally lead to worse reservoir performance. Furthermore, the porosity correlation length has a negligible effect on oil production when the standard deviation is relatively low, while it tends to reduce the amount of oil for a large standard deviation. Especially as the standard deviation was raised to 0.1, the smallest sweep region occurred with early CO₂ breakthrough, which led to a worse flooding effect. Moreover, in order to find a way to deal with strong heterogeneity, some sensitivity assessments relating to production strategies are performed. The results show that higher injection and production rates can improve the flooding effect and enhance oil recovery. Furthermore, compared with pure CO₂ injection, nitrogen-assisted injection yields higher oil production when the heterogeneity is high. The present study provides useful guidelines for developing EOR and CCUS, particularly for tight oil reservoirs with strong heterogeneity in China to achieve carbon-neutral target.