Coupled Numerical Simulation of CO₂ - EOR Flooding: Integrating Multiphysics Interactions

Abstract

Numerical simulation of CO₂-enhanced oil recovery (CO₂-EOR) critically advances hydrocarbon field development by quantifying multiphase saturation dynamics during CO₂-driven displacement. This study establishes a cylindrically configured 1D triple-phase mathematical model integrating Darcy’s law and mass conservation principles, employing an implicit pressure-explicit saturation (IMPES) finite-difference scheme to resolve spatiotemporal evolution of aqueous, oleic, and gaseous phase saturations. Systematic incorporation of chemical reaction kinetics, viscosity-pressure coupling, and dynamic relative permeability effects yields a novel computational framework for immiscible displacement analysis. Simulations reveal two governing mechanisms: 1) CO₂-saturated fluid/rock interactions induce pore-throat dilation, amplifying effective flooding radii; 2) Wellbore-formation pressure differentials (ΔP) dictate CO₂ plume propagation, where elevated ΔP expands repulsion radii and oil displacement annulus thickness. Paradoxically, increased reservoir porosity reduces annular confinement while diminishing displacement efficiency despite enhanced volumetric throughput, with simulations confirming persistent oil-rich annuli characteristic of non-miscible regimes. These findings provide actionable guidelines for optimizing CO₂-EOR injectivity parameters and ensuring long-term carbon sequestration integrity in heterogeneous formations, bridging theoretical modeling with field-scale implementation strategies.

Key words
CCUS; CO₂-EOR; CCO₂ flooding; numerical simulation; geophysical-chemical reactions; viscosity

Introduction

Carbon Capture, Utilization and Storage (CCUS) technology is one of the key technologies to reduce CO₂ emissions and mitigate the greenhouse effect (Zhang et al. 2024). While conventional CO₂ mitigation strategies primarily focus on emission reduction through energy conservation, CCUS technology demonstrates substantially greater decarbonization efficacy and scalability potential when strategically optimized geological reservoir selection and storage protocols are implemented (Liao et al. 2022). So far, there are over 200 CCUS projects worldwide (Zuloaga-Molero et al. 2016) that leverage CO₂ flooding technology to enhance crude oil properties and reservoir porosity. This is achieved by injecting CO₂ into the reservoir as an innovative replacement agent through injection wells. This approach not only substantially elevates the oilfield recovery rate but also concurrently accomplishes CO₂ emission reduction and sequestration, thereby fostering a mutually beneficial outcome (Wu et al. 2022). The numerical simulation of CO₂ flooding represents a pivotal technology in the realm of oil and gas field development, enabling enhanced understanding and optimization of extraction processes (Ramadhan et al. 2023). To gain a deeper comprehension of the mechanisms and influencing variables of CO₂ flooding, optimizing the flooding strategy, augmenting recovery rates, minimizing costs, and extending the productive lifespan of oilfields, it is exceedingly significant to delve into the mathematical models underpinning the numerical simulation of CO₂ flooding. This endeavor holds immense value and potential for advancing the field.

Figure 1
Schematic diagram of CCUS technology and main types.

Since its inaugural adoption in the Sacroc field of the United States in 1972, CO₂ flooding has garnered extensive research attention from petroleum engineers worldwide. This has led to the formulation of numerous diverse mathematical models aimed at simulating and optimizing the process of CO₂ (Bou-Mikael 1996). Comparative experiments and mathematical simulations conducted in the North Sea oilfield of the United Kingdom, analyzing CO₂, CH₄ and dissolved gas, revealed that the application of CO₂ as a solvent outperformed both CH₄ and dissolved gas solutions by a significant margin (Cardenas et al. 1984). The fluid flow model devised by Gerritsen & Durlofsky (2005) for intricate multiphase and multicomponent systems employs detailed calculations of reservoir heterogeneity to describe the gas injection process within a miscible phase. This model demonstrates that miscible-phase injection can enhance both recovery rates and CO₂ sequestration efficiency (Gerritsen & Durlofsky 2005). Toi et al. (2010) established a mathematical model for the diffusion coefficient of CO₂ in hydrocarbons under ambient temperature and pressure conditions, based on experimental studies involving CO₂ diffusion within hydrocarbons and subsequent measurements of the diffusion coefficients under those standard conditions (Toi et al. 2010). Domestic scholars have conducted extensive research in the theoretical realm of oil repulsion. Cheng et al. (1998) formulated a pioneering 3D multiphase multicomponent reservoir simulator, explicitly incorporating CO₂ aqueous-phase solubility dynamics to resolve oil-gas-water phase interactions in compositional reservoir systems (Cheng et al. 1998). Hou Jian developed a practical mathematical model for CO₂ mixed-phase flooding building upon the refined black oil model. This model employs the adjustment to the permeability and effective viscosity of oil and gas to accurately simulate the mixed-phase process (Hou 2004). Xu Geyuan et al. established a comprehensive multi-phase and multi-component zonal seepage model incorporating the critical factor of diffusion into the analysis (Xu 2011). Gao Ran et al. conducted a simulation encompassing the distribution of CO₂ across oil, gas and water phases, the impact of dissolved CO₂ on reservoir properties, as well as the entire process of CO₂-induced oil displacement and sequestration. This simulation integrated the two-phase flashing of oil and gas, the concurrent dissolution of CO₂ in water, and the reciprocal effects of these phenomena (Gao et al. 2021).

Currently, the majority of reservoir numerical simulation software available in the market can be used to simulate the CO₂ displacement process (Cui et al. 2016), they often fall short in accurately accounting for the multifaceted effects of CO₂. Specifically, these software tools may overlook the influence of CO₂ as a hydrocarbon component in oil and gas phases on the overall crude oil mobility (Wang et al. 2023a, b). Additionally, they may not comprehensively consider the effects of CO₂ dissolution in water, the chemical reactions between carbonic acid and formation minerals, and the subsequent impact of the mineralization process on reservoir permeability (Zhao & Liu 2023). To better simulate the CO₂ flooding process, it is imperative to integrate the enhanced crude oil mobility resulting from CO₂ dissolution, the geophysical and chemical reactions occurring within the formation under specific conditions, and the dynamics of the CO₂ flooding process itself (Wang et al. 2024). This holistic approach ensures a more accurate portrayal of the intricate interactions and effects involved in CO₂-enhanced oil recovery.

In this study, a one-dimensional three-phase mathematical model is formulated under the assumption of a columnar formation. This model serves as the foundational basis for subsequent solution and analysis, enabling the simultaneous simulation of CO₂ distribution within oil and gas phases, the dissolution of CO₂ in water, as well as the reactions that occur under specific formation conditions during the displacement process. Furthermore, the study employs advanced programming techniques to optimize parameters based on the numerical simulation results.

STRATIGRAPHIC PHYSICAL MODEL AND BOUNDARY CONDITIONS

The physical modeling and assumed conditions of the strata in this study are as follows:

1. (1)

   A borehole of radius $r_0$ in the center of an isotropic cylindrical formation is divided into cells of radius $r_e$using an isobaric progression grid.

2. (2)

   There are three-phase isothermal seepage of water, oil, and CO₂ gas in the formation. Oil and water are insoluble in each other, and CO₂ gas can be dissolved in water and oil, and all three phases conform to Darcy’s law;

3. (3)

   The formation rock and fluid are compressible, the influence of capillary force is ignored, and gravity is ignored.

The initial and boundary conditions in this study are as follows:

Distribution of pressure, water, oil, and CO₂ gas three-phase saturation at various locations in the formation at the beginning of the CO₂ flooding.

Pressure:

$P_r{|}_{t=0}=P_0$ , the pressure of the units in the formation, except for the borehole, is equal to the initial pressure of the formation $P_0$;

$P_r{|}_{r=r_0,t=t_0}=P_d$, the pressure of the unit in the borehole is the bottomhole pressure $P_d$.

Saturation:

$S_w\left(r\right){|}_{t=t_0}=S_{w0}$, the water saturation of all units in the formation is equal to the initial water saturation of the formation $S_{w0}$;

$S_o\left(r\right){|}_{t=t_0}=S_{o0}$, the oil saturation of all units in the formation is equal to the initial oil saturation of the formation $S_{o0}$ ;

$S_m\left(r\right){|}_{t=t_0}=S_{m0}=0$, CO₂ gas saturation of all units in the formation is equal to the initial CO₂ gas saturation in the formation $S_{m0}=0$;

Figure 2
Schematic diagram of formation modeling.

MATHEMATICAL MODELING

Derivation of Continuity Equation

The continuity equation is established under the condition of continuity, the so-called continuity condition means that when the fluid flows in the space where it is located, there is no void formation inside, and the establishment of the continuity equation should also follow a basic law of nature - the law of conservation of mass.

Take a cylindrical unit, let it be unit $j$, its height is $H$,inner and outer diameter are $r_{j-1}$, $r_j$, the direction of fluid flow is: from inside to outside. Then along the reservoir flow direction, in the standard environment, in 4 time from the inner face of the water, oil, CO₂ gas three-phase flow into the unit of mass are:

$${\left({\rho }_w{v}_w\right)|}_{r_{j-1}}2\pi r_{j-1}Hdt$$

$${\left({\rho }_o{v}_o\right)|}_{r_{j-1}}2\pi r_{j-1}Hdt$$

$${({\rho }_g{v}_g+R_{so}{\rho }_o{v}_o+R_{sw}{\rho }_w{v}_w)|}_{r_{j-1}}2\pi r_{j-1}Hdt$$

The mass of each of the three phases flowing out of the unit cell from the outer end face is:

$${\left({\rho }_w{v}_w\right)|}_{r_j}2\pi r_{j}Hdt$$

$${\left({\rho }_o{v}_o\right)|}_{r_j}2\pi r_{j}Hdt$$

$${({\rho }_g{v}_g+R_{so}{\rho }_o{v}_o+R_{sw}{\rho }_w{v}_w)|}_{r_j}2\pi r_{j}Hdt$$

Then it can be obtained that after $dt$ time the mass change of the three phases of water, oil and CO₂ in the cell and out of the cell are:

$${\left({\rho }_w{v}_w\right)|}_{r_{j-1}}2\pi r_{j-1}Hdt-{\left({\rho }_w{v}_w\right)|}_{r_j}2\pi r_{j}Hdt$$

$${\left({\rho }_o{v}_o\right)|}_{r_{j-1}}2\pi r_{j-1}Hdt-{\left({\rho }_o{v}_o\right)|}_{r_j}2\pi r_{j}Hdt$$

$${({\rho }_g{v}_g+R_{so}{\rho }_o{v}_o+R_{sw}{\rho }_w{v}_w)|}_{r_{j-1}}2\pi r_{j-1}Hdt-{({\rho }_g{v}_g+R_{so}{\rho }_o{v}_o+R_{sw}{\rho }_w{v}_w)|}_{r_j}2\pi r_{j}Hdt$$

And the overall change in the three-phase fluid mass of water, oil, and CO₂ gas in the cell in $dt$ time, respectively:

$$[{\left(⌀S_w{\rho }_w\right)|}_{t_j}-{\left(⌀S_w{\rho }_w\right)|}_{t_{j-1}}](r_j^2-r_{j-1}^2)\pi H$$

$$[{\left(⌀S_o{\rho }_o\right)|}_{t_j}-{\left(⌀S_o{\rho }_o\right)|}_{t_{j-1}}](r_j^2-r_{j-1}^2)\pi H$$

$$[{(⌀S_g{\rho }_g+R_{so}⌀S_o{\rho }_o+R_{sw}⌀S_w{\rho }_w)|}_{t_j}-{(⌀S_g{\rho }_g+R_{so}⌀S_o{\rho }_o+R_{sw}⌀S_w{\rho }_w)|}_{t_{j-1}}](r_j^2-r_{j-1}^2)\pi H$$

The mass conservation continuity equation is derived as follows, using the CO₂ gas phase as an example:

From the law of conservation of mass: total change in mass in the cell = mass flowing into the cell - mass flowing out of the cell. Then it can be obtained:

$${({\rho }_g{v}_g+R_{so}{\rho }_o{v}_o+R_{sw}{\rho }_w{v}_w)|}_{r_{j-1}}2\pi r_{j-1}Hdt-{({\rho }_g{v}_g+R_{so}{\rho }_o{v}_o+R_{sw}{\rho }_w{v}_w)|}_{r_j}2\pi r_{j}Hdt=$$

$$[{(⌀S_g{\rho }_g+R_{so}⌀S_o{\rho }_o+R_{sw}⌀S_w{\rho }_w)|}_{t_j}-{(⌀S_g{\rho }_g+R_{so}⌀S_o{\rho }_o+R_{sw}⌀S_w{\rho }_w)|}_{t_{j-1}}](r_j^2-r_{j-1}^2)\pi H$$

The left side of the equation $\approx -\frac{\partial }{\partial r}\left[({\rho }_g{v}_g+R_{so}{\rho }_o{v}_o+R_{sw}{\rho }_w{v}_w)r\right]\Delta r2\pi Hdt$

The right side of the equation $\approx \frac{\partial }{\partial t}(⌀S_g{\rho }_g+R_{so}⌀S_o{\rho }_o+R_{sw}⌀S_w{\rho }_w)\pi \Delta r2r_{j}Hdt$

The equation is simplified to give:

$$-\frac{1}{r}\frac{\partial }{\partial r}\left[({\rho }_g{v}_g+R_{so}{\rho }_o{v}_o+R_{sw}{\rho }_w{v}_w)r\right]=\frac{\partial }{\partial t}(⌀S_g{\rho }_g+R_{so}⌀S_o{\rho }_o+R_{sw}⌀S_w{\rho }_w)$$(1)

Similarly, the mass conservation continuity equations for the water and oil phases can be obtained as:

$$-\frac{1}{r}\frac{\partial }{\partial r}\left[{\rho }_w{v}_{w}r\right]=\frac{\partial }{\partial t}\left(⌀S_w{\rho }_w\right)$$(2)

$$-\frac{1}{r}\frac{\partial }{\partial r}\left[{\rho }_o{v}_{o}r\right]=\frac{\partial }{\partial t}\left(⌀S_o{\rho }_o\right)$$(3)

In the process of CO₂-EOR flooding, Oil, gas and water three-phase fluids flow in the pore medium, and their flows obey Darcy’s law. The equation of Darcy’s law in Cartesian coordinate system takes the form:

$$v=-\frac{K{K}_r}{\mu }\frac{\partial p}{\partial r}$$

Substituting into the mass conservation continuity equation:

$$\frac{1}{r}\frac{\partial }{\partial r}\left[{\rho }_w\frac{K{K}_{rw}}{{\mu }_w}r\frac{\partial p_w}{\partial r}\right]=\frac{\partial }{\partial t}\left(⌀S_w{\rho }_w\right)$$(4)

$$\frac{1}{r}\frac{\partial }{\partial r}\left[{\rho }_o\frac{K{K}_{ro}}{{\mu }_o}r\frac{\partial p_o}{\partial r}\right]=\frac{\partial }{\partial t}\left(⌀S_o{\rho }_o\right)$$(5)

$$\frac{1}{r}\frac{\partial }{\partial r}\left[({\rho }_g\frac{K{K}_{rg}}{{\mu }_g}+R_{so}{\rho }_o\frac{K{K}_{ro}}{{\mu }_o}+R_{sw}{\rho }_w\frac{K{K}_{rw}}{{\mu }_w})r\frac{\partial p_g}{\partial r}\right]=\frac{\partial }{\partial t}(⌀S_g{\rho }_g+R_{so}⌀S_o{\rho }_o+R_{sw}⌀S_w{\rho }_w)$$(6)

Pressure Solving

Neglect capillary forces:$p_o=p_w=p_g=P$, this is obtained by substituting ${\lambda }_o=\frac{K{K}_{ro}}{{\mu }_o},{\lambda }_w=\frac{K{K}_{rw}}{{\mu }_w},{\lambda }_o=\frac{K{K}_{rg}}{{\mu }_g}$ into the mass conservation continuity equation:

$$\frac{1}{r}\frac{\partial }{\partial r}\left({\rho }_o{\lambda }_{o}r\frac{\partial P}{\partial r}\right)=\frac{\partial }{\partial t}\left(⌀S_o{\rho }_o\right)$$(7)

$$\frac{1}{r}\frac{\partial }{\partial r}\left({\rho }_w{\lambda }_{w}r\frac{\partial P}{\partial r}\right)=\frac{\partial }{\partial t}\left(⌀S_w{\rho }_w\right)$$(8)

$$\frac{1}{r}\frac{\partial }{\partial r}\left[({\rho }_g{\lambda }_g+R_{so}{\rho }_o{\lambda }_o+R_{sw}{\rho }_w{\lambda }_w)r\frac{\partial P}{\partial r}\right]=\frac{\partial }{\partial t}(⌀S_g{\rho }_g+R_{so}⌀S_o{\rho }_o+R_{sw}⌀S_w{\rho }_w)$$(9)

$$令L_O=\frac{\partial }{\partial t}\left(⌀S_o{\rho }_o\right),L_w=\frac{\partial }{\partial t}\left(⌀S_w{\rho }_w\right),L_g=\frac{\partial }{\partial t}(⌀S_g{\rho }_g+R_{so}⌀S_o{\rho }_o+R_{sw}⌀S_w{\rho }_w)$$

In the mass conservation continuity equation density$({\rho }_o,{\rho }_w,{\rho }_g)$, gas solubility$(R_{sw},R_{so})$, porosity$\left(⌀\right)$ are all functions of pressure$\left(P\right)$, which can be obtained by expanding the derivative phase in $(L_o,L_w,L_g)$ by time:

$$L_O=\frac{\partial }{\partial t}\left(⌀S_o{\rho }_o\right)=⌀{\rho }_o\frac{\partial S_o}{\partial t}+(S_o{\rho }_o\frac{\partial ⌀}{\partial P}+⌀S_o\frac{\partial {\rho }_o}{\partial P})\frac{\partial P}{\partial t}$$(10)

$$L_w=\frac{\partial }{\partial t}\left(⌀S_w{\rho }_w\right)=⌀{\rho }_w\frac{\partial S_w}{\partial t}+(S_w{\rho }_w\frac{\partial ⌀}{\partial P}+⌀S_w\frac{\partial {\rho }_w}{\partial P})\frac{\partial P}{\partial t}$$(11)

$$\begin{array}{c}\hfill L_g=\frac{\partial }{\partial t}(⌀S_g{\rho }_g+R_{so}⌀S_o{\rho }_o+R_{sw}⌀S_w{\rho }_w)\\ \hfill =⌀{\rho }_g\frac{\partial S_g}{\partial t}+(S_g{\rho }_g\frac{\partial ⌀}{\partial P}+⌀S_g\frac{\partial {\rho }_g}{\partial P})\frac{\partial P}{\partial t}+\\ \hfill R_{so}⌀{\rho }_o\frac{\partial S_o}{\partial t}+(R_{so}{S}_o{\rho }_o\frac{\partial ⌀}{\partial P}+R_{so}⌀S_o\frac{\partial {\rho }_o}{\partial P}+⌀S_o{\rho }_o\frac{\partial R_{so}}{\partial P})\frac{\partial P}{\partial t}+\\ \hfill R_{sw}⌀{\rho }_w\frac{\partial S_w}{\partial t}+(R_{sw}{S}_w{\rho }_w\frac{\partial ⌀}{\partial P}+R_{sw}⌀S_w\frac{\partial {\rho }_w}{\partial P}+⌀S_w{\rho }_w\frac{\partial R_{sw}}{\partial P})\frac{\partial P}{\partial t}\end{array}$$(12)

There is $S_o+S_w+S_g=1$, this is obtained by substituting $\frac{\partial S_g}{\partial t}=-\frac{\partial S_o}{\partial t}-\frac{\partial S_w}{\partial t}$ into Eq. (12):

$$L_g=(R_{so}⌀{\rho }_o-⌀{\rho }_g)\frac{\partial S_o}{\partial t}+(R_{sw}⌀{\rho }_w-⌀{\rho }_g)\frac{\partial S_w}{\partial t}+\left(\begin{array}{cc}\hfill & S_g{\rho }_g\frac{\partial ⌀}{\partial P}+⌀S_g\frac{\partial {\rho }_g}{\partial P}+R_{so}{S}_o{\rho }_o\frac{\partial ⌀}{\partial P}+R_{so}⌀S_o\frac{\partial {\rho }_o}{\partial P}+\hfill \\ \hfill & ⌀S_o{\rho }_o\frac{\partial R_{so}}{\partial P}+R_{sw}{S}_w{\rho }_w\frac{\partial ⌀}{\partial P}+R_{sw}⌀S_w\frac{\partial {\rho }_w}{\partial P}+⌀S_w{\rho }_w\frac{\partial R_{sw}}{\partial P}\hfill \end{array}\right)\frac{\partial P}{\partial t}$$

In this way, $\frac{\partial S_o}{\partial t},\frac{\partial S_w}{\partial t}$can be eliminated through$(\frac{1}{{\rho }_o}-\frac{R_{so}}{{\rho }_g})L_o+(\frac{1}{{\rho }_w}-\frac{R_{sw}}{{\rho }_g})L_w+\frac{1}{{\rho }_g}{L}_g$So there’s only one unknown $P$ in the new formula that needs to be solved:

$$(\frac{1}{{\rho }_o}-\frac{R_{so}}{{\rho }_g})L_O+(\frac{1}{{\rho }_w}-\frac{R_{sw}}{{\rho }_g})L_w+\frac{1}{{\rho }_g}{L}_g=$$

$$[\frac{\partial ⌀}{\partial P}+\frac{⌀S_g}{{\rho }_g}\frac{\partial {\rho }_g}{\partial P}+\frac{⌀S_o}{{\rho }_o}\frac{\partial {\rho }_o}{\partial P}+\frac{⌀S_w}{{\rho }_w}\frac{\partial {\rho }_w}{\partial P}+⌀S_o{\rho }_o{\rho }_g\frac{\partial R_{so}}{\partial P}+⌀S_w{\rho }_w{\rho }_g\frac{\partial R_{sw}}{\partial P}]\frac{\partial P}{\partial t}=⌀C_t\frac{\partial P}{\partial t}$$

where $C_t=C_p+S_g{C}_g+S_o(C_o+{\rho }_o{\rho }_g{R}_{so}{C}_{R_{so}})+S_w(C_w+{\rho }_w{\rho }_g{R}_{sw}{C}_{R_{sw}})$

$C_p=\frac{1}{⌀}\frac{\partial ⌀}{\partial P},C_o=\frac{1}{{\rho }_o}\frac{\partial {\rho }_o}{\partial P},C_w=\frac{1}{{\rho }_w}\frac{\partial {\rho }_w}{\partial P},C_g=\frac{1}{{\rho }_g}\frac{\partial {\rho }_g}{\partial P},C_{R_{so}}=\frac{1}{R_{so}}\frac{\partial R_{so}}{\partial P},C_{R_{sw}}=\frac{1}{R_{sw}}\frac{\partial R_{sw}}{\partial P}$ Then the left end of the continuity equation for a three-phase fluid can also be multiplied by the corresponding coefficients and summed:

$$\begin{array}{c}\hfill (\frac{1}{{\rho }_o}-\frac{R_{so}}{{\rho }_g})\left[\frac{1}{r}\frac{\partial }{\partial r}\left({\rho }_o{\lambda }_{o}r\frac{\partial P}{\partial r}\right)\right]+(\frac{1}{{\rho }_w}-\frac{R_{sw}}{{\rho }_g})\left[\frac{1}{r}\frac{\partial }{\partial r}\left({\rho }_w{\lambda }_{w}r\frac{\partial P}{\partial r}\right)\right]\\ \hfill +\frac{1}{{\rho }_g}\frac{1}{r}\frac{\partial }{\partial r}\left[({\rho }_g{\lambda }_g+R_{so}{\rho }_o{\lambda }_o+R_{sw}{\rho }_w{\lambda }_w)r\frac{\partial P}{\partial r}\right]=⌀C_t\frac{\partial P}{\partial r}\end{array}$$(13)

The derivation of cell-specific pressure-solving equations is carried out using the gas phase of cell No. 1 as an example:

The mass of the gas passing through the rear face of unit $j$ is:

$$M_g^j={[-({\rho }_g{\lambda }_g+R_{so}{\rho }_o{\lambda }_o+R_{sw}{\rho }_w{\lambda }_w)r_{j}2\pi H\frac{\partial P}{\partial r}]|}_{r_j}$$

The mass of the gas passing through the front face of unit $j$ is:

$$M_g^{j-1}={[-({\rho }_g{\lambda }_g+R_{so}{\rho }_o{\lambda }_o+R_{sw}{\rho }_w{\lambda }_w)r_{j-1}2\pi H\frac{\partial P}{\partial r}]|}_{r_{j-1}}$$

The gas mass increment for unit $j$ is:$\Delta M_{gj}=M_g^{j-1}-M_g^j$

The volume of unit $j$ is:$V_j=\pi r_j^{2}H-\pi r_{j-1}^{2}H$

Then there are:

$${\frac{1}{r}\frac{\partial }{\partial r}\left[({\rho }_g{\lambda }_g+R_{so}{\rho }_o{\lambda }_o+R_{sw}{\rho }_w{\lambda }_w)r\frac{\partial P}{\partial r}\right]|}_j=\frac{\Delta M_{gj}}{V_j}$$(14)

The proof is as follows:

$$\frac{\Delta M_{gj}}{V_j}=\frac{M_g^{j-1}-M_g^j}{V_j}=\frac{{{\left[({\rho }_g{\lambda }_g+R_{so}{\rho }_o{\lambda }_o+R_{sw}{\rho }_w{\lambda }_w)r_{j}2\pi H\frac{\partial P}{\partial r}\right]|}_{r_j}-\left[({\rho }_g{\lambda }_g+R_{so}{\rho }_o{\lambda }_o+R_{sw}{\rho }_w{\lambda }_w)r_{j-1}2\pi H\frac{\partial P}{\partial r}\right]|}_{r_{j-1}}}{\pi r_j^{2}H-\pi r_{j-1}^{2}H}$$

$$=\frac{{{\left[({\rho }_g{\lambda }_g+R_{so}{\rho }_o{\lambda }_o+R_{sw}{\rho }_w{\lambda }_w)r_j\frac{\partial P}{\partial r}\right]|}_{r_j}-\left[({\rho }_g{\lambda }_g+R_{so}{\rho }_o{\lambda }_o+R_{sw}{\rho }_w{\lambda }_w)r_{j-1}\frac{\partial P}{\partial r}\right]|}_{r_{j-1}}}{\left(\frac{r_j+r_{j-1}}{2}\right)(r_j-r_{j-1})}$$

$$=\frac{{\frac{\partial }{\partial r}\left[({\rho }_g{\lambda }_g+R_{so}{\rho }_o{\lambda }_o+R_{sw}{\rho }_w{\lambda }_w)r\frac{\partial P}{\partial r}\right]|}_j(r_j-r_{j-1})}{\left(\frac{r_j+r_{j-1}}{2}\right)(r_j-r_{j-1})}\approx {\frac{1}{r}\frac{\partial }{\partial r}\left[({\rho }_g{\lambda }_g+R_{so}{\rho }_o{\lambda }_o+R_{sw}{\rho }_w{\lambda }_w)r\frac{\partial P}{\partial r}\right]|}_j$$

Ditto:

$${\frac{1}{r}\frac{\partial }{\partial r}\left({\rho }_o{\lambda }_{o}r\frac{\partial P}{\partial r}\right)|}_j=\frac{\Delta M_{oj}}{V_j}$$(15)

$${\frac{1}{r}\frac{\partial }{\partial r}\left({\rho }_w{\lambda }_{w}r\frac{\partial P}{\partial r}\right)|}_j=\frac{\Delta M_{wj}}{V_j}$$(16)

In this way, Eq. (13) can discretize the time($t$) from $t_n$ to $t_{n+1}$ and the radius($r$) from $r_{j-1}$ to $r_j$. And the formula can be simplified:

$$\begin{array}{c}\hfill \frac{{\left[\left({\rho }_o{\lambda }_o\right)r_j\frac{\partial P}{\partial r}\right]|}_{r_j}}{{\rho }_{oj}}-\frac{{\left[\left({\rho }_o{\lambda }_o\right)r_{j-1}\frac{\partial P}{\partial r}\right]|}_{r_{j-1}}}{{\rho }_{oj}}+\frac{{(R_{soj}-R_{soj-1})\left[\left({\rho }_o{\lambda }_o\right)r_{j-1}\frac{\partial P}{\partial r}\right]|}_{r_{j-1}}}{{\rho }_{gj}}+\\ \hfill \frac{{\left[\left({\rho }_w{\lambda }_w\right)r_j\frac{\partial P}{\partial r}\right]|}_{r_j}}{{\rho }_{wj}}-\frac{{\left[\left({\rho }_w{\lambda }_w\right)r_{j-1}\frac{\partial P}{\partial r}\right]|}_{r_{j-1}}}{{\rho }_{wj}}+\frac{{(R_{swj}-R_{swj-1})\left[\left({\rho }_w{\lambda }_w\right)r_{j-1}\frac{\partial P}{\partial r}\right]|}_{r_{j-1}}}{{\rho }_{gj}}+\\ \hfill \frac{{\left[\left({\rho }_g{\lambda }_g\right)r_j\frac{\partial P}{\partial r}\right]|}_{r_j}}{{\rho }_{gj}}-\frac{{\left[\left({\rho }_g{\lambda }_g\right)r_{j-1}\frac{\partial P}{\partial r}\right]|}_{r_{j-1}}}{{\rho }_{gj}}=\frac{V_j⌀C_t}{2\pi H}\frac{\partial P}{\partial r}\end{array}$$(17)

Eq. (17) is logarithmically transformed ($r=r_w{e}^x$) and discretized in the direction of radius:

$$\begin{array}{c}\hfill ({\lambda }_{oj}+{\lambda }_{wj}+{\lambda }_{gj})P_{j+1}^{n+1}+\left(\begin{array}{cc}\hfill & -{\lambda }_{oj}-\frac{{\rho }_{oj-1}{\lambda }_{oj-1}}{{\rho }_{oj}}+\frac{(R_{soj}-R_{soj-1}){\rho }_{oj-1}{\lambda }_{oj-1}}{{\rho }_{gj}}-{\lambda }_{wj}-\frac{{\rho }_{wj-1}{\lambda }_{wj-1}}{{\rho }_{wj}}\hfill \\ \hfill & +\frac{(R_{swj}-R_{swj-1}){\rho }_{wj-1}{\lambda }_{wj-1}}{{\rho }_{gj}}-{\lambda }_{gj}-\frac{{\rho }_{gj-1}{\lambda }_{gj-1}}{{\rho }_{gj}}-\frac{\Delta x{V}_j⌀C_t}{2\pi H\Delta t}\hfill \end{array}\right)P_j^{n+1}\\ \hfill +(\frac{{\rho }_{oj-1}{\lambda }_{oj-1}}{{\rho }_{oj}}+\frac{{\rho }_{wj-1}{\lambda }_{wj-1}}{{\rho }_{wj}}+\frac{{\rho }_{gj-1}{\lambda }_{gj-1}}{{\rho }_{gj}})P_{j-1}^{n+1}=-\frac{\Delta x{V}_j⌀C_t}{2\pi H\Delta t}{P}_j^n\end{array}$$(18)

Saturation Solving

The difference equation for saturation can be obtained by substituting Eq. (14)(15)(16) into Eq. (9)(7)(8). The newly obtained difference formula for solving saturation can discretize the time ($t$) from $t_n$ to $t_{n+1}$. And the formula can be simplified:

$$\frac{{{\left[\left({\rho }_o{\lambda }_o\right)2\pi H\frac{\partial P}{\partial x}\right]|}_{r_j}-\left[\left({\rho }_o{\lambda }_o\right)2\pi H\frac{\partial P}{\partial x}\right]|}_{r_{j-1}}}{\pi r_j^{2}H-\pi r_{j-1}^{2}H}=\frac{[{\left({\rho }_o⌀S_o\right)}_j^{n+1}-{\left({\rho }_o⌀S_o\right)}_j^n]}{\Delta t}$$(19)

$$\frac{{{\left[\left({\rho }_w{\lambda }_w\right)2\pi H\frac{\partial P}{\partial x}\right]|}_{r_j}-\left[\left({\rho }_w{\lambda }_w\right)2\pi H\frac{\partial P}{\partial x}\right]|}_{r_{j-1}}}{\pi r_j^{2}H-\pi r_{j-1}^{2}H}=\frac{[{\left({\rho }_w⌀S_w\right)}_j^{n+1}-{\left({\rho }_w⌀S_w\right)}_j^n]}{\Delta t}$$(20)

$$\begin{array}{c}\hfill \frac{{\left[({\rho }_g{\lambda }_g+R_{so}{\rho }_o{\lambda }_o+R_{sw}{\rho }_w{\lambda }_w)2\pi H\frac{\partial P}{\partial x}\right]|}_{r_j}-{\left[({\rho }_g{\lambda }_g+R_{so}{\rho }_o{\lambda }_o+R_{sw}{\rho }_w{\lambda }_w)2\pi H\frac{\partial P}{\partial x}\right]|}_{r_{j-1}}}{\pi r_j^{2}H-\pi r_{j-1}^{2}H}\\ \hfill =\frac{[{(⌀{\rho }_g{S}_g+R_{so}⌀S_o{\rho }_o+R_{sw}⌀S_w{\rho }_w)}_j^{n+1}-{(⌀S_g{\rho }_g+R_{so}⌀S_o{\rho }_o+R_{sw}⌀S_w{\rho }_w)}_j^n]}{\Delta t}\end{array}$$(21)

Sorting out can be obtained water, oil, CO₂ gas three-phase saturation solution formula:

$$S_{oj}^{n+1}=\frac{{\left({\rho }_o⌀S_o\right)}_j^n}{{\left({\rho }_o⌀\right)}_j^{n+1}}+\frac{2\Delta t[{\left({\rho }_o{\lambda }_o\right)}_j^{n+1}\frac{P_{j+1}^{n+1}-P_j^{n+1}}{\Delta x}-{\left({\rho }_o{\lambda }_o\right)}_{j-1}^{n+1}\frac{P_j^{n+1}-P_{j-1}^{n+1}}{\Delta x}]}{{\left({\rho }_o⌀\right)}_j^{n+1}(r_j^2-r_{j-1}^2)}$$(22)

$$S_{wj}^{n+1}=\frac{{\left({\rho }_w⌀S_w\right)}_j^n}{{\left({\rho }_w⌀\right)}_j^{n+1}}+\frac{2\Delta t[{\left({\rho }_w{\lambda }_w\right)}_j^{n+1}\frac{P_{j+1}^{n+1}-P_j^{n+1}}{\Delta x}-{\left({\rho }_w{\lambda }_w\right)}_{j-1}^{n+1}\frac{P_j^{n+1}-P_{j-1}^{n+1}}{\Delta x}]}{{\left({\rho }_w⌀\right)}_j^{n+1}(r_j^2-r_{j-1}^2)}$$(23)

$$\begin{array}{c}\hfill S_g^{n+1}=\frac{{({\rho }_g⌀S_g+R_{so}{S}_o{\rho }_o⌀+R_{sw}{S}_w{\rho }_w⌀)}_j^n-{(R_{so}{S}_o{\rho }_o⌀+R_{sw}{S}_w{\rho }_w⌀)}_j^{n+1}}{{\left({\rho }_g⌀\right)}_j^{n+1}}+\\ \hfill \frac{2\Delta t[{({\rho }_g{\lambda }_g+R_{so}{\rho }_o{\lambda }_o+R_{sw}{\rho }_w{\lambda }_w)}_j^{n+1}\frac{P_{j+1}^{n+1}-P_j^{n+1}}{\Delta x}-{({\rho }_g{\lambda }_g+R_{so}{\rho }_o{\lambda }_o+R_{sw}{\rho }_w{\lambda }_w)}_{j-1}^{n+1}\frac{P_j^{n+1}-P_{j-1}^{n+1}}{\Delta x}]}{{\left({\rho }_g⌀\right)}_j^{n+1}(r_j^2-r_{j-1}^2)}\end{array}$$(24)

Geophysical-chemical Reactions

In the process of expulsion, the special reactions between carbonic acid generated by the dissolution of CO₂ in water and formation minerals, as well as the effect on formation porosity during mineralization (Dong 2022), all of which are altered by geochemical reactions. The change rule of porosity of the formation medium with effective volume strain is:

$$⌀=\frac{1}{1+\epsilon }[(1+{\epsilon }_0){⌀}_0+{\alpha }_{\epsilon }(\epsilon -{\epsilon }_0)]$$(25)

where $\epsilon ={\epsilon }_V+\frac{P}{K_S}-{\epsilon }_{s1}-{\epsilon }_{s2}$ is the effective volumetric strain, which is the combined result of geochemical volumetric strain (${\epsilon }_{s2}$), adsorption strain (${\epsilon }_{s1}$), pore pressure-induced strain ($\frac{P}{K_S}$), and volumetric strain (${\epsilon }_V$), and ${⌀}_0$ and $⌀$ are the initial and current porosities, respectively. ${\epsilon }_0=\frac{P_0}{K_S}-{\epsilon }_{s1,0}$ is the initial effective volumetric strain. The volumetric strain and chemical volumetric strain in the initial state are both zero. $P_0$ is the initial pore pressure.

For $\epsilon <<1,{\epsilon }_0<<1$ , Eq. (25) is simplified into:

$$⌀={⌀}_0+{\alpha }_{\epsilon }({\epsilon }_V+\frac{P-P_0}{K_s}-{\epsilon }_{s1}-{\epsilon }_{s2})$$(26)

It is clear from Eq. (26) that formation porosity is influenced by the volumetric strain of the formation medium and is related to pore pressure, geochemical volumetric strain, and adsorption-induced volumetric strain (medium expansion).

Correction of Viscosity with Saturation

CO₂ is extremely easy to dissolve in crude oil, so that the viscosity of crude oil is significantly reduced (Shikai 2023). Moreover, the larger the initial viscosity of crude oil is, the more obvious the range of decrease is. The main function of CO₂ to reduce the viscosity of crude oil is to improve the fluidity of crude oil, so that the oil driving efficiency can be achieved with a small amount of driving agent to achieve a certain effect; or with a quantitative amount of driving agent can realize a higher driving effect.

The effect of formation pressure on crude oil viscosity can be broken down into two components: the effect of pressure changes on viscosity in the absence of CO₂ dissolution, and the effect of pressure changes on crude oil viscosity when CO₂ is dissolved and precipitated in the crude oil. These two components do not respond to pressure changes in the same way, which leads to saturation pressure. When the formation pressure is higher than the saturation pressure, the volume change of the crude oil has the main part of the effect on the viscosity of the crude oil, and the volume of the crude oil decreases with the increase of the pressure, which leads to the increase of the viscosity. And when the formation pressure is lower than the saturation pressure, the effect of CO₂ dissolution in crude oil on the viscosity of crude oil accounts for the main part, and as the pressure decreases, the dissolved CO₂ gas precipitates out of the crude oil, leading to the coalescence of heavy hydrocarbon fractions, which manifests itself as a sharp increase in the viscosity of crude oil. In this way, an empirical formula for the variation of viscosity with saturation can be obtained:

$${\mu }_o={\mu }_o^0[\frac{C_{\mu }^a}{P+C_{\mu }^b}+C_{\mu }^{c}P+C_{\mu }^d]$$(27)

Where $C_{\mu }^a$, $C_{\mu }^b$, $C_{\mu }^c$, $C_{\mu }^d$ are the coefficients of variation of crude oil viscosity with the solubility of the CO₂ gas in the crude oil.

Iteration of Penetration

Relative permeability is an intrinsic property used to characterize the ease of fluid flow in porous media. As an important physical parameter, the calculation of relative permeability plays an important role in the numerical simulation of CO₂-EOR flooding. In the case of multiphase fluids seeping through a formation gap (Dong 2020a), the saturation of each phase of the fluid jointly determines the distribution of the flow path (Dong 2020b). This allows the relative permeability to be corrected by the saturation degree and iteratively calculated. The relationship between relative permeability and saturation is as follows:

$$K_{rw}=K_{rworgr}{\left(\frac{S_w-S_{wc}}{1-S_{wc}-S_{or}-S_{gr}}\right)}^{n_w}$$(28)

Where $K_{rworgr}$ is the relative permeability of the water phase of residual oil saturation and retained gas saturation, $S_{or}$ is the residual oil saturation, $S_{gr}$ is the retained gas saturation, $S_{wc}$ is the bound water saturation. The denominator $(1-S_{wc}-S_{or}-S_{gr})$ is the maximum range of water saturation that can be changed, and the numerator $(S_w-S_{wc})$ is the mobile water saturation. $n_w$ is the calculation coefficient of the relative permeability of the water phase.

$$K_{ro}=K_{rocwgr}{\left(\frac{S_w-S_{or}}{1-S_{wc}-S_{or}-S_{gr}}\right)}^{n_o}$$(29)

Where $K_{rocwgr}$ is the relative permeability of oil phase under bound water saturation and retained gas saturation. The denominator $(1-S_{wc}-S_{or}-S_{gr})$ is the maximum range of oil saturation that can be changed, and the numerator $(1-S_w-S_{or})$ is the mobile oil saturation. $n_o$ is the calculation coefficient of the relative permeability of the oil phase.

$$K_{rg}=K_{rgcwor}{\left(\frac{S_g-S_{gr}}{1-S_{wc}-S_{or}-S_{gr}}\right)}^{n_g}$$(30)

Where, $K_{rgcwro}$ is the residual oil saturation and the relative permeability of the gas phase under mature confinement. The denominator $(1-S_{wc}-S_{or}-S_{gr})$ is the maximum range of gas saturation that can be changed, and the numerator $(S_g-S_{gr})$ is the mobile gas saturation. $n_g$ is the calculation coefficient of the relative permeability of gas phase.

Flowchart for Model Solving.

NUMERICAL SIMULATION RESULTS AND ANALYSIS

Initial and boundary conditions

The boundary conditions of this simulation are constant-pressure inner boundary, constant-pressure outer boundary, and the outer boundary pressure 500m away from the center of the borehole is the initial formation pressure. This simulation analyzes the saturation degree of the three of phases water, oil and CO₂ gas in the CO₂-EOR flooding process and the parameters that have influence on the CO₂-EOR flooding efficiency regularly. The parameters and parameter conditions of this simulation are shown in the following Table I:

Table I
Physical parameters of CO₂–EOR flooding simulation

Parameters	Values	Parameters	Values
Initial formation pressure $P_0$	5 MPa	Borehole radius $r_0$	0.1 m
Borehole pressure $P_w$	15 MPa	Viscosity of oil ${\mu }_o$	3.03 mPa·s
Absolute formation permeability $K$	50 md	Viscosity of water ${\mu }_w$	1 mPa·s
Formation porosity $\phi$	0.15	Viscosity of CO₂ ${\mu }_g$	0.2 mPa·s
Initial water saturation $S_{\mathrm{w0}}$	0.35	Bound water saturation $S_{\mathrm{wc}}$	0.15
Initial oil saturation $S_{\mathrm{o0}}$	0.65	Residual oil saturation $S_{\mathrm{or}}$	0.1
Initial CO₂ saturation $S_{\mathrm{g0}}$	0	Stranded CO₂ saturation $S_{\mathrm{gr}}$	0
Formation thickness $H$	5 m	Pore compression coefficient $C_P$	1×10⁻⁴ 1/Pa
Water relative permeability index $n_w$	1.3	Water compressibility $C_w$	1×10⁻⁴ 1/Pa
Oil relative permeability index $n_o$	1.3	Oil compressibility $C_o$	1×10⁻⁴ 1/Pa
CO₂ relative permeability index $n_m$	1.3	CO₂ compressibility $C_g$	1×10⁻⁴ 1/Pa
Strain parameter ${\alpha }_{\epsilon }$	0.5	CO₂ viscosity compression coefficient $C_m$	1×10⁻⁴ 1/Pa
Volume strain ${\epsilon }_V$	0.066	Crude oil viscosity coefficient $C^a{}_{}$	1.06466
Volume compression modulus of particles $K_S$	160.25	Crude oil viscosity coefficient $C^b$	-1.90388
Adsorption strain${\epsilon }_{\mathrm{s1}}$	0.027	Crude oil viscosity coefficient $C^c$	0.00856
Geochemical volumetric strain ${\epsilon }_{\mathrm{s2}}$	0.031	Crude oil viscosity coefficient $C^d$	0.45314

Influence analysis of parameters

Figure 3 depicts saturation profiles of water, oil, and CO₂ gas prior- and post-geophysical-chemical reaction, with radius on the x-axis and saturation on the y-axis. Notably, the front and back sections of all three phases remain unchanged, while the middle sections shift left: water and oil saturations rise, while CO₂ gas saturation decreases. This is attributed to dissolved CO₂ reacting with the formation, enlarging pores and enabling water and oil to erode the formation more effectively, thereby increasing their saturations and reducing CO₂ saturation. This aligns with findings in Wang et al. (2023a).

Figure 3
Three phase saturation prior- and post-modification for geophysical chemical reaction.

Figure 4 juxtaposes saturation profiles of water, oil, and CO₂ gas prior- and post-correction for viscosity, plotting radius against saturation. Water saturation exhibits a leftward shift upfront and marginally midway, with an unchanged increasing rate and unchanged backend. Conversely, oil saturation exhibits a rightward shift in both frontal and intermediate stages, with peak attenuation attributed to enhanced mobility that disrupts complete immiscible-phase interaction. Concurrently, CO₂ saturation displays similar front-loaded displacement patterns accompanied by accelerated depletion rates in early-to-mid phases, while maintaining stable saturation levels in later stages. This phenomenon stems from CO₂ solubility significantly reducing oil viscosity, thereby enhancing fluidity, CO₂ repulsion efficiency, expanding the repulsion radius, and consequently right-shifting the oil saturation curve. Furthermore, the heightened oil mobility accelerates the non-mixed-phase CO₂ flooding in increased water and CO₂ saturations, which is in harmony with the theory presented in (Huang 2017).

Figure 4
The saturation of three phases prior- and post-modification for viscosity formula.

Analysis of simulation results

Figure 5 illustrates saturation of water oil and CO₂ gas of the initial day. The x-axis is formation displacement radius, and the y-axis is the saturations. As CO₂ gas initially contacts formation fluids, non-mixed-phase repulsion commences. Some CO₂ dissolves in crude oil, expanding its volume and boosting oil saturation, the oil saturation bulges evidently. The remaining CO₂ forms an upper phase, partly dissolving in water, with water saturation increasing more in the first half of contact. Most upper-phase CO₂ migrates forward, extracting hydrocarbons from oil, but insufficiently to reach mixed-phase conditions. Early breakthrough of formation fluids by the upper phase diminishes CO₂ flooding efficiency, aligning with findings in (Kang et al. 2014).

Figure 5
Saturation of three phases on the first day of CO₂ flooding.

Figure 6 displays saturation profiles of water, oil, and CO₂ gas on the 15th day. The x-axis is formation displacement radius, and the y-axis is the saturation. The initial section of each curve is linear, indicating bound water, residual oil, and stagnant CO₂ saturation. As intrusion depth increases, water saturation rises to initial levels, oil saturation peaks from residual to initial, then declines, while CO₂ saturation drops to zero. This signifies non-mixed-phase replacement with CO₂ displacing fluids. The oil saturation bulge reflects CO₂ dissolution, expanding oil volume and saturation. The unaffected formation’s rear section maintains initial three-phase saturation levels.

Figure 6
Saturation of three phases on the 15th day of CO₂ flooding.

Figure 7 shows saturation curves of water, oil, and CO₂ gas vs. flooding time. The x-axis is the radius, and the y-axis is the saturation. The curves depict water saturation over various flooding days. With non-mixed-phase replacement, the expelled radius expands with time, but its growth rate decelerates. Water saturation increases from bound to initial levels over longer periods with diminishing rates. Oil saturation rises from residual to peak, then declines to initial levels over longer durations with slower rates, exhibiting a more pronounced exclusion zone effect. CO₂ saturation declines from residual to initial levels over longer spans with slower rates. CO₂ displaces mobile water and oil, leaving bound water, residual oil, and CO₂ in pores. The horizontal parallel segments at curve ends indicate unaffected strata with initial saturation values.

Figure 7
Three-phase saturation on different days during CO₂ flooding.

Figure 8 illustrates the saturation dynamics of water, oil, and CO₂ gas in three phases versus injection time at varying radii during CO₂ flooding. The horizontal axis is injection time, and the vertical axis is the saturations. The figure reveals that saturation changes commence at different times and rates for units at different locations. Initially, near-wellbore units exhibit rapid water saturation decline to bound water levels within 10 days, while distant units take longer with slower declines. Oil saturation near the wellbore surges to a peak before declining to residual levels, with faster changes compared to distant units. CO₂ saturation increases from zero, driving fluids away from the wellbore, eventually filling all pore space except for residual oil and bound water. Larger radii units experience longer unmixed-phase displacement durations. Oil saturation takes longer to reach residual levels compared to water saturation reaching bound levels at the same location.

Figure 8
Three-phase saturation with time at different radii during CO₂ flooding.

Figure 9 depicts the saturation profiles of water, oil, and CO₂ gas in a three-phase system under varying differential pressures between the wellbore and the formation during a simulation of oil displacement. The horizontal axis represents the radius of formation displacement, and the vertical axis indicates the saturation levels of water, oil, and CO₂ gas in these three phases. The differential pressures between the bottomhole flowing pressure and the formation’s initial pressure are set at 1 MPa, 2.5 MPa, 5 MPa, 7.5 MPa, and 10 MPa, respectively. The numerical simulation results clearly indicate that, at a given simulation time, a larger differential pressure between the bottomhole flowing pressure and the formation’s initial pressure leads to a greater radius of CO₂ displacement and a more pronounced intrusion zone. This is attributable to the fact that the pressure differences the two-phase fluids of water and oil to infiltrate into the formation’s pore spaces. When the formation pressure remains constant, an increase in the bottomhole pressure intensifies the rate of CO₂ displacement, resulting in a larger displacement radius.

Figure 9
Saturation of three phases of different injection pressure during CO₂ flooding.

Figure 10 illustrates the saturation profiles of water, oil, and CO₂ gas in varying formation porosities during an oil displacement simulation. The horizontal axis represents the radial displacement within the formation, while the vertical axis indicates the saturations of water, oil, and CO₂ gas. The formation porosities are set at 0.1, 0.15, 0.2, 0.25, and 0.3, respectively. The saturation curves for these three phases reveal that at a given simulation time, a greater pressure differential between the wellbore and the formation leads to a larger radius of CO₂ displacement and a more pronounced intrusion zone. This phenomenon stems from the fact that the pressure differential between the wellbore and the formation facilitates the infiltration of water and oil, two-phase fluids, into the formation’s pore spaces. When the formation pressure remains constant, an increase in the bottomhole flowing pressure intensifies the rate of CO₂ displacement, leading to a larger displacement radius.

Figure 10
Saturation of three phases at different porosities during CO₂ flooding.

CONCLUSIONS

1. (1)

   Based on Darcy’s law and mass conservation principles, this study develops a 1D triple-phase mathematical framework to characterize CO₂-EOR flooding dynamics, systematically investigating the evolution of formation pressure and phase saturations under variable operational conditions. As immiscible displacement progresses, formation pressure exhibits a monotonic increase accompanied by a diminishing formation-borehole pressure differential. Oil saturation undergoes transient elevation to peak oil saturation before reverting to initial reservoir conditions, while aqueous phase saturation manifests a distinct stepwise distribution pattern across radial distances. Concomitantly, CO₂ gaseous saturation progressively diminishes to baseline levels, demonstrating characteristic radial propagation signatures of CO₂ flooding.

2. (2)

   An equation linking porosity changes to effective volumetric strain is incorporated, elucidating the interplay between porosity, pore pressure, and various strains induced by geophysical-chemical reactions. An empirical viscosity-pressure relation is fitted to experimental data, enhancing the characterization of crude oil viscosity’s sensitivity to volume changes, CO₂ dissolution, and other factors. This contributes to a more physically meaningful numerical simulation of CO₂-EOR flooding. Furthermore, a theoretical relative permeability model for water, oil, and CO₂ phases is presented, facilitating accurate representation of three-phase flow in simulations.

3. (3)

   Observations reveal a stepwise water saturation distribution at varying radial distances, indicative of non-mixed-phase zones. Oil saturation exhibits a distinct replacement ring, peaking at full saturation prior declining to formation levels. CO₂ gas saturation gradually diminishes to zero, underscoring the characteristics of radial displacement.

4. (4)

   Reduction of viscosity highlights enhanced oil mobility, expanding the driving radius but inhibiting full CO₂-oil contact, leading to reduced oil saturation peaks. This accelerates unmixed-phase CO₂ flooding, shrinking the annulus and boosting water and CO₂ saturations. CO₂-infused water and oil phases react with the formation, enlarging pores and fostering increased water and oil saturations, ultimately decreasing CO₂ saturation, flooding radius, and recovery efficiency.

ACKNOWLEDGMENTS

This work is supported by Provincial Frontier Technology Research and Development Plan (grant No.BF2024018).

References

• BOU-MIKAEL SAM. 1996. A New analytical method to Evaluate, Predict, and Improve CO₂ Flood Performance in Sandst one Reservoirs. In: SPE/DOE Improved Oil Recovery Symposium. Tulsa, Oklahoma. doi: 10.2118/35362-MS.
• CARDENAS RL, ALSTON RB, NUTE AJ & KOKOLIS GP. 1984. Laboratory Design of a Gravity-Stable Miscible CO₂ Process. J Petrol Technol 3: 111-118.
• CHENG LS, LI CL & CHEN YM. 1998. Numerical simulation for calculating CO₂ solubility in water phase. J Petrol Univ (Natural Science Edition 2(22): 38-41.
• CUI GD, REN SR & DOU B. 2016. Geothermal exploitation from depleted high temperature gas reservoirs via recycling supercritical CO₂: Heat mining rate and salt precipitation effects. Appl Energ 183: 837-852.
• DONG X, SHEN LW & LIU XF. 2020a. NMR characterization of a tight sand’s pore structures and fluid mobility: an experimental investigation for CO₂ EOR potential. Mar Petrol Geol 118: 104460.
• DONG X, SHEN LW & NASER G. 2020b. How N2 injection improves the hydrocarbon recovery of CO₂ HnP: An NMR study on the fluid displacement mechanisms. Fuel 278: 118286.
• DONG X, SHEN LW & BO L. 2022. Fracture’s Impact on the Recovery of Hydrocarbon from Low-Permeability Rock’s Pores: New Insights from 1H Nuclear Magnetic Resonance Experiment. SPE J 27: 2913-2925.
• GAO R, CY LV & LUN ZM. 2021. Integrated Numerical Simulation of Carbon Dioxide Displacement and Sequestration. Special Oil Gas Reserv 28(02): 102-107.
• GERRITSEN MG & DURLOFSKY LJ. 2005. Modeling fluid flow in oil reservoirs. Annu Rev Fluid Mech 37: 211-238.
• HOU J. 2004. Streamline-Based Mathematical Model for CO₂ Miscible Flooding. Appl Math Mech 25(6): 635-642.
• HUANG BP. 2017. Study on CO₂ immiscible flooding technology in low permeability reservoir of Daqing Oilfield Northeast Petroleum University.
• KANG H, GAO J & SONG XM. 2014. Study on Oil Saturation Distribution of CO₂ Drive in Low Permeability Reservoir. Sci Technol Eng 14(14): 198-201.
• LIAO GZ, HE DB & WANG GF. 2022. Discussion on the limit recovery factor of carbon dioxide flooding in a permanent sequestration scenario. Petrol Explor Develop 49(6): 1463-1470.
• RAMADHAN R, ABDURRAHMAN M & BISSEN R. 2023. Numerical simulation of potential site for CO₂ sequestration in a depleted oil reservoir in northern Thailand. Energy Rep 9(S11): 524-528.
• SHIKAI M. 2023. Numerical simulation study of CO₂ flooding in low permeability reservoirs of Block H. Northeast Petroleum University. doi: 10.26995/d.cnki.gdqsc.2023.000875.
• TOI K, OHORI Y & MAEDA YE. 2010. Effect of molecular weight on the diffusion coefficient of carbon dioxide in polystyrene. J Poly Sci Part A Poly Chem 18(7): 364-375.
• WANG D, LI J & LIAN W. 2024. Impact of CO₂ sequestration on pressure systems in saline aquifer: Simulation of differential pressure evolution in the caprock-reservoir. Geoenerg Sci Eng 236: 212748.
• WANG JG, WANG HM & WANG XL. 2023a. A multiphysical-geochemical coupling model for caprock sealing efficiency in CO₂ geosequestration. Deep Undergr Sci Eng 2(02): 188-203.
• WANG K, MA L & TAYLOR K. 2023b. Microstructure changes as a response to CO₂ storage in sedimentary rocks: Recent developments and future challenges. Fuel 333: 126403.
• WU BZ, YUUAN XY & GUO TZ. 2022. Method for extracting apparent resistivity from transient electromagnetic logging measurements and geosteering potential. J China Univ Petrol, Edition of Natural Science 46: 99-109.
• XU GY. 2011. CO₂-EOR multiphase multicomponent zonal seepage model considering diffusion. Oil Gas Field Surf Eng 2(30): 24-25.
• ZHANG ZC, BER MX & DU SY. 2024. Research progress of CO₂ displacement and storage integration technology. World Petrol Industr 32(01): 100-107.
• ZHAO C & LIU QB. 2023. Investigating effects of structural deformation regimes on mineralization distributions in fluid-saturated rocks: computational simulation approach through generic models. Minerals 13(5): 664.
• ZULOAGA-MOLERO P, YU W & XU YF. 2016. Simulation Study of CO₂-EOR in Tight Oil Reservoirs with Complex Fracture Geometries. Sci Rep 6: 334-345.

Publication Dates

• Publication in this collection
  15 Dec 2025
• Date of issue
  2025

History

• Received
  25 Feb 2025
• Accepted
  01 June 2025