A novel numerical simulation of CO₂ immiscible flooding coupled with viscosity and starting pressure gradient modeling in ultra-low permeability reservoir

Abstract

CO₂ immiscible flooding is an environmentally-friendly and effective method to enhance oil recovery in ultra-low permeability reservoirs. A mathematical model of CO₂ immiscible flooding was developed, considering the variation in crude oil viscosity and starting pressure gradient in ultra-low permeability reservoirs based on the non-Darcy percolation theory. The mathematical model and numerical simulator were developed in the C++ language to simulate the effects of fluid viscosity, starting pressure gradient, and other physical parameters on the distribution of the oil pressure field, oil saturation field, gas saturation field, oil viscosity field, and oil production. The results showed that the formation pressure and pressure propagation velocity in CO₂ immiscible flooding were lower than the findings without considering the starting pressure gradient. The formation oil content saturation and the crude oil formation viscosity were higher after the consideration of the starting pressure gradient. The viscosity of crude oil considering the initiation pressure gradient during the formation was higher than that without this gradient, but the yield was lower than that condition. Our novel mathematical models helped the characterization of seepage resistance, revealed the influence of fluid property changes on seepage, improved the mathematical model of oil seepage in immiscible flooding processes, and guided the improvement of crude oil recovery in immiscible flooding processes.

Graphical abstract

Keywords

viscosity / starting pressure gradient / flow simulation / CO₂ immiscible flooding / ultra-low permeability reservoir

1 Introduction

In recent years, ultra-low permeability and tight reservoirs have attracted much attention. CO₂ injection is a clean and effective method considered to improve the recovery of ultra-low permeability reservoirs or depleted and water-flooded reservoirs (Ju et al., 2020; Cui et al., 2022). Previous studies showed that there were two basic characteristics in the production process of CO₂ in ultra-low permeability reservoirs (Javadpour, 2009; Civan, 2010; Pertsin and Grunze, 2004). The first property was the change in oil viscosity during CO₂ flooding. The second characteristic was the variation in the oil starting pressure gradient. These two properties made the fluid flow a typical low-speed non-Darcy flow in the CO₂ development process of ultra-low permeability and tight reservoir. Numerous studies and numerical experiments have been performed on the effects of low-velocity non-Darcy flows (Wang, 2021; Allen and Sun, 2012; Battiatoa and Tartakovskya, 2011).

CO₂ immiscible flooding is a process of replacing crude oil by dissolving some CO₂ in the crude oil. It causes a decrease in the viscosity of crude oil and the interfacial tension while the volume of crude oil expands. The application of CO₂ immiscible flooding mainly includes 1) restoring the pressure of a depleted reservoir by CO₂, 2) gravity stabilization of CO₂ immiscible flooding, and 3) exploitation of thick oil reservoirs (IPCC, 2005; Firoozabadi and Myint, 2010; Michael et al., 2011; Wang, 2011).

After the CO₂ dissolution in crude oil, the viscosity of crude oil decreases significantly depending on the pressure, temperature, and initial viscosity of crude oil. Generally, the higher the initial viscosity of crude oil, the higher the percentage of viscosity reduction after CO₂ dissolution. Further enhancement in pressure after the saturation of oil with CO₂ increases the oil viscosity due to compression.

The starting pressure gradient of crude oil in CO₂ immiscible flooding is closely related to the viscosity of crude oil. The starting pressure gradient (Gao et al., 2021; Huang, 1998; Zhang et al., 1998) is the pressure gradient when the pressure drawdown begins. It acts on fluid porous media and has attained a certain level to help the fluid overcome the viscous force and start to flow. If the pressure gradient is small, the flow velocity increases slowly and obeys the nonlinear law. However, the flow velocity rises quickly and obeys the linear law if the pressure gradient exceeds the starting pressure gradient. In the last decade, Chinese researchers have conducted a large number of mathematical simulations and empirical studies on starting pressure gradient and low-velocity non-Darcy flow.

Many accomplishments have been obtained through the numerical model of low-velocity non-Darcy flow for ultra-low permeability reservoirs (An et al., 2017; Li et al., 2016). (Wang and Sheng 2017a) deduced a non-Darcy flow model according to the results of Huang (1998). The single author citation format was written incorrectly and has been re-cited for the corresponding reference In this regard, the performance of a vertical well and a multi-fractured horizontal well was studied by adopting this model. The terminal oil recovery from non-Darcy flow was almost 48% of that from Darcy flow for a vertical well, and 80% of a multi-fractured horizontal well in their simulation. Lu et al. (2017) pointed out that fluid flow in a low permeability reservoir was non-Darcy flow with the starting pressure gradient. The new mathematical models for boundary-dominated flow under two different conditions were established. The pressure derivative was a concave curve during transient flow, which went upwards. The bigger the starting pressure gradient, the higher the flow resistance. Davarpanah (2020) focused on the driving mathematical equation of recovery and established a mathematical model of polymer-assisted nanoparticle migration.

The fine mesh and precise prediction fluid properties play an important role in underground conditions due to the complexity of underground systems and involved physical processes. In the literature review, oil viscosity and oil starting pressure gradient in CO₂ immiscible flooding were not taken into account. The previous mathematical model was insufficient for changes in the physical and seepage properties of the fluid during CO₂ immiscible flooding in the low-permeability reservoir.

In this study, the influence of oil viscosity and fluid starting pressure gradient was mainly considered for CO₂ immiscible flooding. The independent programming implementation of a numerical simulator for CO₂ immiscible flooding was presented on a microcomputer by coupling with the corrected models of viscosity and starting pressure gradient for a simulated reservoir. The main purpose of this work was to study the mechanisms of immiscible gas flooding to more accurately reveal the characteristics of fluid flow and oil displacement in immiscible gas flooding. It was investigated by considering the change in fluid physical properties in the seepage system. In this regard, the medium or high permeability reservoir was regarded as a special case of the ultra-low permeability reservoir.

2 Model description

2.1 Viscosity and starting pressure gradient corrected models

The relative permeability of oil and gas (k_ro and k_rg) is a function of gas saturation (s_g). These parameters can be obtained by fitting the phase permeability curve to field-measured cores:

(1)k_{\mathrm{r}\mathrm{g}}=0.06\cdot {\left(\frac{s_{\mathrm{g}}}{1-s_{\mathrm{o}\mathrm{r}}}\right)}^{2.39} ,

(2)k_{\mathrm{r}\mathrm{o}}=1.232\cdot {\left(1-\frac{s_{\mathrm{g}}}{1-s_{\mathrm{o}\mathrm{r}}}\right)}^{3.353} .

The oil phase starting pressure gradient (G_o) is a function of oil mobility \left({\displaystyle \frac{k_{\mathrm{o}}}{{\mu }_{\mathrm{o}}}}\right) (Zhang et al., 1998; Lv et al., 2002). It can be determined by the regression of the core starting pressure gradient experimental data (Chi, 2017):

(3)G_{\mathrm{o}}=0.586\times {\left(\frac{k_{\mathrm{o}}}{{\mu }_{\mathrm{o}}}\right)}^{-0.857}.

The crude oil viscosity (µ_o) is a function of pressure (P) and temperature (T):

(4){\mu }_{\mathrm{o}}={\mu }_{\mathrm{o}}(P,T ).

In the process of displacement, the physical properties of the fluid were constantly changing during the formation due to the constant dissolution of CO₂ in crude oil. Therefore, the oil phase viscosity and oil phase starting pressure gradient needed correction as in (Chung et al., 1988):

1) The CO₂ solubility in the crude oil formation

The solubility of CO₂ (R_s) in crude oil formation is defined as the standard volume unit of dissolved CO₂ in surface crude oil under reservoir conditions. The solubility (R_s) is a function of formation pressure (P), temperature (T), and crude oil density (γ), in which (P) and (T) have the highest effects, but (γ) has a low influence:

(5)R_s=\{0.178[a_1{\gamma }^{a_2}{T}^{a_7}+a_3{T}^{a_4}\exp (-a_{5}p-a_6/p)]{\}}^{-1},

among them

(6)T=T_1\times \frac{\text{9}}{\text{5}} +\text{32},

(7)p=p_1\times 145,

(8)a_1=0.4936\times {10}^{-2},

(9)a_2=4.0928,

(10)a_3=0.571\times {10}^{-6},

(11)a_4=1.6428,

(12)a_5=0.6763\times {10}^{-3},

(13)a_6=-781.1334 ,

(14)a_7=-\mathrm{0.2499.}

2) The corrected viscosity of crude oil (µ_om) can be determined using the temperature and pressure during formation by Eq. (15):

(15)\ln {\mu }_{\mathrm{o}\mathrm{m}}=X_{\mathrm{o}}\ln {\mu }_{\mathrm{o}}+X_{\mathrm{s}}\ln {\mu }_{\mathrm{g}},

among them

(16)X_{\mathrm{s}}=V_{\mathrm{s}} /(\alpha V_{\mathrm{o}}+V_{\mathrm{s}}),

(17)X_{\mathrm{o}}=1-X_{\mathrm{s}},

(18)\alpha =0.25r^{-4.16}{T}_{\mathrm{r}}^{1.85}\frac{\exp (7.36)-\exp (7.36-7.36p_{\mathrm{r}})}{\exp (7.36)-1},

(19)T_{\mathrm{r}}=(1.8T +32)/547.57,p_{\mathrm{r}}=0.1354p,

where V_o is the volume fraction of the oil phase, V_s is the gas volume fraction, µ_g and µ_o are the viscosity of the gas phase and oil phase, respectively.

X_s can be calculated from Eq. (20):

(20)X_{\mathrm{s}}=\frac{1}{\alpha F_{{\mathrm{C}\mathrm{O}}_2}/(5.618F_{\mathrm{o}}{R}_{\mathrm{s}})+1} =\frac{F_{\mathrm{o}}{F}_{\mathrm{s}}}{\alpha +F_{\mathrm{o}}{F}_{\mathrm{s}}-1}.

where F_{{\mathrm{C}\mathrm{O}}_2} is the ratio of CO₂ volume under standard conditions to that under reservoir temperature and pressure, and F_o is the ratio of the volume of crude oil at reservoir temperature and standard pressure of 0.1 MPa to the volume at reservoir temperature and pressure.

The corrected crude oil viscosity (µ_om) is obtained:

(21){\mu }_{\mathrm{o}\mathrm{m}}={\mu }_{\mathrm{o}\mathrm{m}}(P,T).

Therefore, the oil phase starting pressure gradient (G_o) can also be corrected as (G_om)

(22)G_{\mathrm{o}\mathrm{m}}=G_{\mathrm{o}\mathrm{m}}\left(\frac{k_{\mathrm{o}}}{{\mu }_{\mathrm{o}\mathrm{m}}}\right).

Equations (1)–(22) constitute the corrected mathematical models of viscosity and starting pressure gradient.

2.2 The CO₂ immiscible flooding model coupled with the viscosity and start-up pressure gradient model

The subsurface-fluid properties are dependent on the reservoir conditions such as pressure, temperature, and fluid composition. The reservoir conditions are also affected by the fluid properties. Therefore, the prediction of the fluid properties under reservoir conditions plays a key role in the reservoir simulation. In this work, the viscosity and start-up pressure gradient models provided the correction of corresponding parameters coupled with the numerical model of CO₂ immiscible flooding to evaluate the fluid properties under certain reservoir conditions. It increased the accuracy of seepage simulation. The coupler is shown in Fig.1.

The viscosity model and the start-up pressure gradient model described the change in fluid properties. They were coupled with the numerical simulation model of carbon dioxide immiscible flooding, which could better describe the seepage process. Additionally, the modules like the stress-sensitive module of porous media were further associated with this coupling framework.

2.2.1 Assumption conditions

The primary assumptions at the beginning of this study were as follows:

1) the fluid in the reservoir was isothermal seepage;

2) the rock was a micro-compressible non-homogeneous pore medium;

3) the fluid was a compressible fluid;

4) the change in crude oil viscosity was considered;

5) the variation in fluid initiation pressure gradient was evaluated;

6) gravity and capillary forces were not taken into account.

2.2.2 The seepage equation and auxiliary equations

The equation \mathrm{\nabla }={\displaystyle \frac{\mathrm{\partial }}{\mathrm{\partial }x}}+{\displaystyle \frac{\mathrm{\partial }}{\mathrm{\partial }y}} is used for the convenience of writing. The two-dimensional CO₂ immiscible flooding seepage mathematical model is as follows.

Oil phase seepage equation

(23)\mathrm{\nabla }\cdot \left[\frac{k{k}_{\mathrm{r}\mathrm{o}}}{{\mu }_{\mathrm{o}}}{\rho }_{\mathrm{o}}(\mathrm{\nabla }{p}_{\mathrm{o}}-G_{\mathrm{o}})\right]+ q_{\mathrm{o}}=\frac{\mathrm{\partial }}{\mathrm{\partial }t}(\varphi {\rho }_{\mathrm{o}} s_{\mathrm{o}}),

Gas phase seepage equation

(24)\begin{array}{rl}& \mathrm{\nabla }\cdot \left[\frac{R_{\mathrm{s}\mathrm{o}}k{k}_{\mathrm{r}\mathrm{o}}}{{\mu }_{\mathrm{o}}}{\rho }_{\mathrm{o}}(\mathrm{\nabla }{p}_{\mathrm{o}}-G_{\mathrm{o}})\right]\text{+}\mathrm{\nabla }\cdot \left[\frac{k{k}_{\mathrm{r}\mathrm{g}}}{{\mu }_{\mathrm{g}}}{\rho }_{\mathrm{g}}\mathrm{\nabla }{p}_{\mathrm{g}}\right]+q_{\mathrm{g}}\\ & =\frac{\mathrm{\partial }}{\mathrm{\partial }t}\left[\varphi ({\rho }_{\mathrm{o}}{s}_{\mathrm{o}}{R}_{\mathrm{s}\mathrm{o}}+{\rho }_{\mathrm{g}}{s}_{\mathrm{g}})\right],\end{array}

where s_o is the oil saturation, s_g is the gas content saturation, ρ_o is the oil phase density, ρ_g is the gas phase density, R_so is the ratio of dissolved gas to oil, B_o is the oil phase volume factor, B_g is the gas phase volume factor, k_ro is the relative permeability of the oil phase, k_rg is the relative permeability of the gas phase, µ_o is the oil phase viscosity, µ_g is the gas phase viscosity, p_o is the oil phase pressure, p_g is the gas phase pressure, q_o is the volume of oil output per time and per volume, and q_g is the volume of injected gas per time and per volume.

Auxiliary equations:

(25)s_{\mathrm{o}}+s_{\mathrm{g}}=1.

Capillary force equation:

(26)p_{\mathrm{o}}=p_{\mathrm{g}}-p_{\mathrm{c}\mathrm{g}\mathrm{o}}.

The equation for relative permeability in the oil phase:

(27)k_{\mathrm{r}\mathrm{o}}=f( s_{\mathrm{g}}).

The equation for relative permeability in the gas phase:

(28)k_{\text{rg}}=f(s_{\mathrm{g}} ).

The equation for CO₂ solubility in crude oil:

(29)R_{\mathrm{s}}=\{0.178[a_1{\gamma }^{a_2}{T}^{a_7}+a_{\text{3}}{T}^{a_4}\exp (-a_{5}p-a_6/p)]{\}}^{-1}.

The calculation of the CO₂ volume fraction integration factor and the correction of crude oil viscosity:

(30)X_{\text{s}} =\frac{1}{\alpha F_{\text{c}{\text{o}}_{\text{2}}}/(5.618F_{\text{o}}{R}_{\text{s}})+1}=\frac{F_{\text{o}}{F}_{\text{s}}}{\alpha +F_{\text{o}}{F}_{\text{s}}-1},

(31)\ln {\mu }_{\text{om}}=X_{\text{o}}\ln {\mu }_{\text{o}} +X_{\text{s}}\ln {\mu }_{\text{g}}.

The correction of the fluid starting pressure gradient:

(32)G_{\text{i}} =a{\left(\frac{K_{\text{i}}}{{\mu }_{\text{i}}}\right)}^{\text{b}}.

2.2.3 Fixed solution conditions

The solution conditions were divided into initial conditions and boundary conditions. The boundary conditions were classified into outer boundary conditions and inner boundary conditions. The outer boundary conditions indicated the pressure at the reservoir boundary and whether it was closed or not. However, the inner boundary conditions represented the state of the injection and extraction wells.

The initial conditions are

(33)\{\begin{array}{c}{p}_o(x,y,0)=p_{\mathrm{o}\mathrm{i}}\\ s_g(x,y,0)=s_{\mathrm{g}\mathrm{i}}\end{array}\begin{array}{cc}& \end{array}(0⩽x⩽L_x,0⩽y⩽L_y).

where p_oi is the initial oil pressure, and s_gi is the initial gas saturation.

There are two types of external boundary conditions.

1) Closed outer boundary

(34)\{\begin{array}{c}{\left({\displaystyle \frac{\mathrm{\partial }{p}_o}{\mathrm{\partial }x}}\right)}_{x=0}=0\\ {\left({\displaystyle \frac{\mathrm{\partial }{p}_o}{\mathrm{\partial }x}}\right)}_{x=L_x}=0\\ {\left({\displaystyle \frac{\mathrm{\partial }{p}_o}{\mathrm{\partial }y}}\right)}_{y=0}=0\\ {\left({\displaystyle \frac{\mathrm{\partial }{p}_o}{\mathrm{\partial }y}}\right)}_{y=L_y}=0\end{array}\begin{array}{cc}& \end{array}(t>0).

2) Outer boundary with constant pressure

(35)\{\begin{array}{rl}& p(0,y,t)=p_{\mathrm{e}}\\ & p(L_x,y,t)=p_{\mathrm{e}}\\ & p(x,0,t)=p_{\mathrm{e}}\\ & p(x,L_y,t)=p_{\mathrm{e}}\end{array}\begin{array}{cc}& \end{array}(t>0).

The external boundary of constant pressure is shown in Fig.2.

The inner boundary conditions were also divided into two types.

1) The fixed yield

(36)Q_{\mathrm{v}\mathrm{l}}=\mathrm{C}\mathrm{o}\mathrm{n}\mathrm{s}\mathrm{t}\mathrm{a}\mathrm{n}\mathrm{t}\mathrm{s},\mathrm{l}=\mathrm{o},\mathrm{g}.

2) The constant pressure at the hole’s bottom flow

The p_gf is known for producing wells, and p_igf is used for gas injection wells. The generation of producing wells can be indicated as follows:

(37)Q_{\mathrm{v}\mathrm{o}\mathrm{i},\mathrm{j}}=P{I}_{\mathrm{o}}(p_{\mathrm{o}\mathrm{i},\mathrm{j}}- p_{\mathrm{w}\mathrm{f}})

The injection of gas can be expressed as

(38)Q_{\mathrm{v}\mathrm{g}\mathrm{i},\mathrm{j}}=G{I}_{\mathrm{g}}(p_{\mathrm{i}\mathrm{w}\mathrm{f}-}{p}_{\mathrm{g}\mathrm{i},\mathrm{j}}),

where PI_o is the oil phase production index and GI_g is the gas phase injection index.

The oil phase seepage equation, the gas phase seepage equation, auxiliary equations, initial conditions, and boundary conditions constitute together a mathematical model of CO₂ immiscible flooding seepage considering the variation in crude oil viscosity and the change in fluid initiation pressure gradient.

Equations (1)–(38) constitute the CO₂ immiscible flooding model coupled with the models of viscosity and starting pressure gradient.

3 Solution methodology

The developed mathematical model was solved by numerical methods. The solution was carried out using the implicit pressure explicit satiation (IMPES) method, which was a sequential solution method to solve multiphase seepage flows. The continuous partial differential equations were first discretized using the finite difference method. The discretized equations are linearized and solved by solving the linear system of equations.

The capillary force Eq. (26) is substituted into the gas phase percolation Eq. (24). The oil phase and gas phase percolation equations are rewritten as

(39)\mathrm{\nabla }\cdot \left[\frac{k{k}_{\mathrm{r}\mathrm{o}}}{{\mu }_{\mathrm{o}}}{\rho }_{\mathrm{o}}(\mathrm{\nabla }{p}_{\mathrm{o}}-G_{\mathrm{o}})\right]+ q_{\mathrm{o}}=\frac{\mathrm{\partial }}{\mathrm{\partial }t}(\varphi {\rho }_{\mathrm{o}} s_{\mathrm{o}}),

(40)\begin{array}{rl}& \mathrm{\nabla }\cdot \left[\frac{k{k}_{\mathrm{r}\mathrm{o}}}{{\mu }_{\mathrm{o}}}{\rho }_{\mathrm{o}}{R}_{\mathrm{s}\mathrm{o}}(\mathrm{\nabla }{p}_{\mathrm{o}}-G_{\mathrm{o}})\right]\text{+}\mathrm{\nabla }\cdot \left[\frac{k{k}_{\mathrm{r}\mathrm{g}}}{{\mu }_{\mathrm{g}}}{\rho }_{\mathrm{g}}(\mathrm{\nabla }{p}_{\mathrm{o}}+\mathrm{\nabla }{p}_{\mathrm{c}\mathrm{g}\mathrm{o}})\right]+q_{\mathrm{g}}\\ & =\frac{\mathrm{\partial }}{\mathrm{\partial }t}\left[\varphi ({\rho }_{\mathrm{o}}{s}_{\mathrm{o}}{R}_{\mathrm{s}\mathrm{o}}+{\rho }_{\mathrm{g}}{s}_{\mathrm{g}})\right].\end{array}

The flow coefficients are defined:

(41){\lambda }_{\mathrm{o}}=\frac{k{k}_{\mathrm{r}\mathrm{o}}{\rho }_{\mathrm{o}}}{{\mu }_{\mathrm{o}}},

(42){\lambda }_{\mathrm{g}}=\frac{k{k}_{\mathrm{r}\mathrm{g}}{\rho }_{\mathrm{g}}}{{\mu }_{\mathrm{g}}}.

Therefore, the oil phase and gas phase percolation equations can be simplified as

(43)\mathrm{\nabla }\cdot [{\lambda }_{\mathrm{o}}(\mathrm{\nabla }{p}_{\mathrm{o}}-G_{\mathrm{o}})]+ q_{\mathrm{o}}=\frac{\mathrm{\partial }}{\mathrm{\partial }t}(\varphi {\rho }_{\mathrm{o}} s_{\mathrm{o}}),

(44)\begin{array}{rl}& \mathrm{\nabla }\cdot [{\lambda }_{\mathrm{o}}{R}_{\mathrm{s}\mathrm{o}}(\mathrm{\nabla }{p}_{\mathrm{o}}-G_{\mathrm{o}})\text{+}{\lambda }_{\mathrm{g}}(\mathrm{\nabla }{p}_{\mathrm{o}}+\mathrm{\nabla }{p}_{\mathrm{c}\mathrm{g}\mathrm{o}})]+q_{\mathrm{g}}\\ & =\frac{\mathrm{\partial }}{\mathrm{\partial }t}\left[\phi ({\rho }_{\mathrm{o}}{s}_{\mathrm{o}}{R}_{\mathrm{s}\mathrm{o}}+{\rho }_{\mathrm{g}}{s}_{\mathrm{g}})\right].\end{array}

The oil phase percolation equations are first discretized in time and space using a finite difference unequal distance mesh:

(45)\begin{array}{rl}& \frac{1}{\mathrm{\Delta }{x}_i}[{\lambda }_{\mathrm{o}xi+\frac{1}{2},j}^n\left(\frac{p_{\mathrm{o}i+1,j}^{n+1}-p_{\mathrm{o}i,j}^{n+1}}{\mathrm{\Delta }{x}_{i+\frac{1}{2},j}}-G_{\mathrm{w}xi+\frac{1}{2},j}^n\right)\\ & -{\lambda }_{\mathrm{o}xi-\frac{1}{2},j}^n\left(\frac{p_{\mathrm{o}i,j}^{n+1}-p_{\mathrm{o}i-1,j}^{n+1}}{\mathrm{\Delta }{x}_{i-\frac{1}{2},j}}-G_{\mathrm{w}xi-\frac{1}{2},j}^n\right)]\\ & +\frac{1}{\mathrm{\Delta }{y}_i}[{\lambda }_{\mathrm{o}yi+\frac{1}{2},j}^n\left(\frac{p_{\mathrm{o}i+1,j}^{n+1}-p_{\mathrm{o}i,j}^{n+1}}{\mathrm{\Delta }{y}_{i+\frac{1}{2},j}}-G_{\mathrm{w}yi+\frac{1}{2},j}^n\right)\\ & -{\lambda }_{\mathrm{o}yi-\frac{1}{2},j}^n\left(\frac{p_{\mathrm{o}i,j}^{n+1}-p_{\mathrm{o}i-1,j}^{n+1}}{\mathrm{\Delta }{y}_{i-\frac{1}{2},j}}-G_{\mathrm{w}yi-\frac{1}{2},j}^n\right)]\\ & +q_{\mathrm{o}}^{n+1}=\frac{1}{\mathrm{\Delta }t}[(\varphi {\rho }_{\mathrm{o}}{s}_{\mathrm{o}}{)}^{n+1}-(\varphi {\rho }_{\mathrm{o}}{s}_{\mathrm{o}}{)}^n].\end{array}

Both sides of the equation are multiplied by the same mesh volume. The volume of mesh (i, j) is

(46)V_{ij}=x_{\mathrm{i}}.y_{\mathrm{i}}.h.

Equation (45) can be written as

(47)\begin{array}{rl}& \mathrm{\Delta }{y}_j[{\lambda }_{\mathrm{o}xi+\frac{1}{2},j}^n\left(\frac{p_{\mathrm{o}i+1,j}^{n+1}-p_{\mathrm{o}i,j}^{n+1}}{\mathrm{\Delta }{x}_{i+\frac{1}{2},j}}-G_{\mathrm{w}xi+\frac{1}{2},j}^n\right)\\ & -{\lambda }_{\mathrm{o}xi-\frac{1}{2},j}^n\left(\frac{p_{\mathrm{o}i,j}^{n+1}-p_{\mathrm{o}i-1,j}^{n+1}}{\mathrm{\Delta }{x}_{i-\frac{1}{2},j}}-G_{\mathrm{w}xi-\frac{1}{2},j}^n\right)]\\ & +\mathrm{\Delta }{x}_i[{\lambda }_{\mathrm{o}yi+\frac{1}{2},j}^n\left(\frac{p_{\mathrm{o}i+1,j}^{n+1}-p_{\mathrm{o}i,j}^{n+1}}{\mathrm{\Delta }{y}_{i+\frac{1}{2},j}}-G_{\mathrm{w}yi+\frac{1}{2},j}^n\right)\\ & -{\lambda }_{\mathrm{o}yi-\frac{1}{2},j}^n\left(\frac{p_{\mathrm{o}i,j}^{n+1}-p_{\mathrm{o}i-1,j}^{n+1}}{\mathrm{\Delta }{y}_{i-\frac{1}{2},j}}-G_{\mathrm{w}yi-\frac{1}{2},j}^n\right)]\\ & +q_{\mathrm{o}}^{n+1}{V}_{ij}=\frac{V_{ij}}{\mathrm{\Delta }t}[(\varphi {\rho }_{\mathrm{o}}{s}_{\mathrm{o}}{)}^{n+1}-(\varphi {\rho }_{\mathrm{o}}{s}_{\mathrm{o}}{)}^n]\end{array}.

The oil phase conduction coefficients in the x and y directions are defined, respectively:

(48)T_{\mathrm{o}xi\pm \frac{1}{2}}^n=\frac{\mathrm{\Delta }{y}_j}{\mathrm{\Delta }{x}_{i\pm \frac{1}{2},j}}{\lambda }_{\mathrm{o}xi\pm \frac{1}{2},j}^n,

(49)T_{\mathrm{o}yj\pm \frac{1}{2}}^n=\frac{\mathrm{\Delta }{x}_j}{\mathrm{\Delta }{y}_{i,j\pm \frac{1}{2}}} {\lambda }_{\mathrm{o}yi,j\pm \frac{1}{2}}^n.

Therefore, the oil phase percolation equation becomes

(50)\begin{array}{rl}& T_{\mathrm{o}xi+\frac{1}{2},j}^n\left[(p_{\mathrm{o}i+1,j}^{n+1}-p_{\mathrm{o}i,j}^{n+1})-G_{\mathrm{o}xi+\frac{1}{2},j}^n\mathrm{\Delta }{x}_{i+\frac{1}{2},j}\right]\\ & +T_{\mathrm{o}xi-\frac{1}{2},j}^n\left[(p_{\mathrm{o}i,j}^{n+1}-p_{\mathrm{o}i-1,j}^{n+1})-G_{\mathrm{o}xi-\frac{1}{2},j}^n\mathrm{\Delta }{x}_{i-\frac{1}{2},j})\right]\\ & +T_{\mathrm{o}yi,j+\frac{1}{2}}^n\left[(p_{\mathrm{o}i,j+1}^{n+1}-p_{\mathrm{o}i,j}^{n+1})-G_{\mathrm{o}yi,j+\frac{1}{2}}^n\mathrm{\Delta }{y}_{i,j+\frac{1}{2}}\right]\\ & +T_{\mathrm{o}yi,j-\frac{1}{2}}^n\left[(p_{\mathrm{o}i,j+1}^{n+1}-p_{\mathrm{o}i,j}^{n+1})-G_{\mathrm{o}yi,j-\frac{1}{2}}^n\mathrm{\Delta }{y}_{i,j-\frac{1}{2}}\right]\\ & +{\text{Q}}_{\mathrm{o}}^{n+1}=\frac{V_{ij}}{\mathrm{\Delta }t}\left[(\varphi {\rho }_{\mathrm{o}}{s}_{\mathrm{o}}{)}^{n+1}-(\varphi {\rho }_{\mathrm{o}}{s}_{\mathrm{o}}{)}^n\right],\end{array}

where the extraction of mesh (i, j) is expressed by

(51)Q_{\mathrm{o}}^{n+1}=q_{\mathrm{o}}^{n+1}{V}_{ij}.

From the above equation, the starting pressure gradient term is combined to obtain:

(52)\begin{array}{rl}& \left[T_{\mathrm{o}xi+\frac{1}{2},j}^n(p_{\mathrm{o}i+1,j}^{n+1}-p_{\mathrm{o}i,j}^{n+1})+T_{\mathrm{o}xi-\frac{1}{2},j}^n(p_{\mathrm{o}i,j}^{n+1}-p_{\mathrm{o}i-1,j}^{n+1})\right]\\ & -\left(T_{\mathrm{o}xi+\frac{1}{2},j}^n{G}_{\mathrm{o}xi+\frac{1}{2},j}^n\mathrm{\Delta }{x}_{i+\frac{1}{2},j}+T_{\mathrm{o}xi-\frac{1}{2},j}^n{G}_{\mathrm{o}xi-\frac{1}{2},j}^n\mathrm{\Delta }{x}_{i-\frac{1}{2},j}\right)\\ & +\left[T_{\mathrm{o}yi,j+\frac{1}{2}}^n\left(p_{\mathrm{o}i,j+1}^{n+1}-p_{\mathrm{o}i,j}^{n+1}\right)+T_{\mathrm{o}yi,j-\frac{1}{2}}^n\left(p_{\mathrm{o}i,j+1}^{n+1}-p_{\mathrm{o}i,j}^{n+1}\right)\right]\\ & -\left(T_{\mathrm{o}yi,j+\frac{1}{2}}^n{G}_{\mathrm{o}yi,j+\frac{1}{2}}^n\mathrm{\Delta }{y}_{i,j+\frac{1}{2}}+T_{\mathrm{o}yi,j-\frac{1}{2}}^n{G}_{\mathrm{o}yi,j-\frac{1}{2}}^n\mathrm{\Delta }{y}_{i,j-\frac{1}{2}}\right)\\ & +{\text{Q}}_{\mathrm{o}}^{n+1}=\frac{V_{ij}}{\mathrm{\Delta }t}\left[(\varphi {\rho }_{\mathrm{o}}{s}_{\mathrm{o}}{)}^{n+1}-(\varphi {\rho }_{\mathrm{o}}{s}_{\mathrm{o}}{)}^n\right].\end{array}

For simplicity, the following difference operator is introduced:

(53){\mathrm{\Delta }}_x\eta ={\eta }_{i+1,j}-{\eta }_{i,j},

(54){\mathrm{\Delta }}_y\eta ={\eta }_{i,j+1}-{\eta }_{i,j},

(55)\mathrm{\Delta }( {\xi }_x{\mathrm{\Delta }}_x \eta )={\xi }_{xi+\frac{1}{2},j}( {\eta }_{i+1,j}- {\eta }_{i,j})+ {\xi }_{xi-\frac{1}{2},j}( {\eta }_{i,j}-{\eta }_{i-1,j}),

(56)\mathrm{\Delta }( {\xi }_y{\mathrm{\Delta }}_y \eta )={\xi }_{yi,j+\frac{1}{2}}({\eta }_{i,j+1}- {\eta }_{i,j})+ {\xi }_{yi,j-\frac{1}{2}}( {\eta }_{i,j}-{\eta }_{i,j-1}).

Therefore, a linear difference operator can be presented:

(57)\mathrm{\Delta }( T_{\mathrm{o}x}{\mathrm{\Delta }}_{x}p)=T_{\mathrm{o}xi+\frac{1}{2}}(p_{i+1}-p_i)+T_{\mathrm{o}xi-\frac{1}{2}}(p_i-p_{i-1}),

(58)\mathrm{\Delta }( T_{\mathrm{o}y}{\mathrm{\Delta }}_{y}p)=T_{\mathrm{o}yj+\frac{1}{2}}(p_{i+1}-p_i)+T_{\mathrm{o}yj-\frac{1}{2}}(p_j-p_{\mathrm{j}-1} ).

The oil phase difference equation can be abbreviated as

(59)\begin{array}{rl}& {\mathrm{\Delta }}_x(T_{\mathrm{o}x}^n{\mathrm{\Delta }}_x{p}_{\mathrm{o}}^{n+1})-[(T_{\mathrm{o}x}{G}_{\mathrm{o}x}{)}_{i+\frac{1}{2},j}^n\mathrm{\Delta }{x}_{i+\frac{1}{2},j}\\ & +(T_{\mathrm{o}x}{G}_{\mathrm{o}x}{)}_{i-\frac{1}{2},j}^n\mathrm{\Delta }{x}_{i-\frac{1}{2},j}]\\ & +{\mathrm{\Delta }}_y(T_{\mathrm{o}y}^n{\mathrm{\Delta }}_y{p}_{\mathrm{o}}^{n+1})-[(T_{\mathrm{o}y}{G}_{\mathrm{o}y}{)}_{i,j+\frac{1}{2}}^n\mathrm{\Delta }{x}_{i,j+\frac{1}{2}}\\ & +(T_{\mathrm{o}y}{G}_{\mathrm{o}y}{)}_{i,j-\frac{1}{2}}^n\mathrm{\Delta }{y}_{i,j-\frac{1}{2}}]\\ & +{\text{Q}}_{\mathrm{o}}^{n+1}=\frac{V_{ij}}{\mathrm{\Delta }t}[(\varphi {\rho }_{\mathrm{o}}{s}_{\mathrm{o}}{)}^{n+1}-(\varphi {\rho }_{\mathrm{o}}{s}_{\mathrm{o}}{)}^n]\end{array}.

In the two-dimensional plane, the linear difference operator can be expanded as

(60)\begin{array}{rl}& \mathrm{\Delta }(\xi \mathrm{\Delta }\eta {)}_{i,j}={\xi }_{xi+\frac{1}{2},j}({\eta }_{i+1,j}-{\eta }_{i,j})+{\xi }_{xi-\frac{1}{2},j}({\eta }_{i,j}-{\eta }_{i-1,j})\\ & \text{+}{\xi }_{yi,j+\frac{1}{2}}({\eta }_{i,j+1}-{\eta }_{i,j})+{\xi }_{yi,j-\frac{1}{2}}({\eta }_{i,j}-{\eta }_{i,j-1}),\end{array}

therefore

(61)\mathrm{\Delta }T\mathrm{\Delta }p={\mathrm{\Delta }}_x{T}_x{\mathrm{\Delta }}_{x}p+{\mathrm{\Delta }}_y{T}_y{\mathrm{\Delta }}_{y}p.

Equation (59) is used to simplify Eq. (62):

(62)\mathrm{\Delta }( T_{\mathrm{o}}^n\mathrm{\Delta } p_{\mathrm{o}}^{n+1}) \text{+}{R}_{\mathrm{o}}^n{\text{+ Q}}_{\mathrm{o}}^{n+1}=\frac{V_{ij}}{\mathrm{\Delta }t}[(\varphi {\rho }_{\mathrm{o}}{s}_{\mathrm{o}}{)}^{n+1} -(\varphi {\rho }_{\mathrm{o}}{s}_{\mathrm{o}}{)}^n ].

The starting pressure gradient of the oil phase percolation equation is constant. It can be combined as

(63)\begin{array}{rl}{R}_{\mathrm{o}}^n=& -\left[(T_{\mathrm{o}x}{G}_{\mathrm{o}x}{)}_{i+\frac{1}{2},j}^n\mathrm{\Delta }{x}_{i+\frac{1}{2},j}\text{+}(T_{\mathrm{o}x}{G}_{\mathrm{o}x}{)}_{i-\frac{1}{2},j}^n\mathrm{\Delta }{x}_{i-\frac{1}{2},j}\right]\\ & -[(T_{\mathrm{o}y}{G}_{\mathrm{o}y}{)}_{i,j+\frac{1}{2}}^n\mathrm{\Delta }{x}_{i,j+\frac{1}{2}}\text{+}(T_{\mathrm{o}y}{G}_{\mathrm{o}y}{)}_{i,j-\frac{1}{2}}^n\mathrm{\Delta }{y}_{i,j-\frac{1}{2}}].\end{array}

In the same way, as the oil phase seepage equation is derived, the gas phase seepage differential equation is according to the following form:

(64)\begin{array}{rl}& \mathrm{\Delta }[(T_{\mathrm{o}}^n{R}_{\mathrm{s}\mathrm{o}}\mathrm{\Delta }{p}_{\mathrm{o}}^{n+1}\text{+}{T}_{\mathrm{g}}^n(\mathrm{\Delta }{p}_{\mathrm{o}}^{n+1}\text{+}\mathrm{\Delta }{p}_{\mathrm{c}\mathrm{g}\mathrm{o}}^n\text{)] +}{R}_{\mathrm{o}\mathrm{g}}^n{\text{+ Q}}_{\mathrm{g}}^{n+1}\\ & =\frac{V_{ij}}{\mathrm{\Delta }t}[(\varphi {\rho }_{\mathrm{o}}{s}_{\mathrm{o}}{R}_{\mathrm{s}\mathrm{o}}+\varphi {\rho }_{\mathrm{g}}{s}_{\mathrm{g}}{)}^{n+1}-(\varphi {\rho }_{\mathrm{o}}{\mathrm{s}}_{\mathrm{o}}{R}_{\mathrm{s}\mathrm{o}}+\varphi {\rho }_{\mathrm{g}}{s}_{\mathrm{g}}{)}^n].\end{array}

where T_o is the oil phase conduction coefficient and T_g is the gas phase conduction coefficient.

The gas injection of mesh (i, j) is

(65)Q_{\mathrm{g}}^{n+1}=q_{\mathrm{g}}^{n+1}{V}_{ij}.

The starting pressure gradient of the oil-gas mixture in the gas phase percolation equation is varied. It can be combined as

(66)\begin{array}{rl}{R}_{\mathrm{o}\mathrm{g}}=& -\left[(T_{\mathrm{o}x}{G}_{\mathrm{o}x}{)}_{i+\frac{1}{2},j}^n\mathrm{\Delta }{x}_{i+\frac{1}{2},j}\text{+}(T_{\mathrm{o}x}{G}_{\mathrm{o}x}{)}_{i-\frac{1}{2},j}^n\mathrm{\Delta }{x}_{i-\frac{1}{2},j}\right]\\ & -\left[(T_{\mathrm{o}y}{G}_{\mathrm{o}y}{)}_{i,j+\frac{1}{2}}^n\mathrm{\Delta }{x}_{i,j+\frac{1}{2}}\text{+}(T_{\mathrm{o}y}{G}_{\mathrm{o}y}{)}_{i,j-\frac{1}{2}}^n\mathrm{\Delta }{y}_{i,j-\frac{1}{2}}\right].\end{array}

The differential equation of oil phase seepage (62) and gas phase seepage (64) are added. The total pressure Eq. (67) is obtained for the unknown quantity p_o as follows:

(67)\begin{array}{rl}& \mathrm{\Delta }(T_{\mathrm{o}}^n\mathrm{\Delta }{p}_{\mathrm{o}}^{n+1})\text{+}\mathrm{\Delta }[(T_{\mathrm{o}}^n{R}_{\mathrm{s}\mathrm{o}}\mathrm{\Delta }{p}_{\mathrm{o}}^{n+1}\text{+}{T}_{\mathrm{g}}^n(\mathrm{\Delta }{p}_{\mathrm{o}}^{n+1}\text{+}\mathrm{\Delta }{p}_{\mathrm{c}\mathrm{g}\mathrm{o}}^n\text{)]}\\ & +R_{\mathrm{o}}^n\text{+}{R}_{\mathrm{o}\mathrm{g}}^n{\text{+ Q}}_{\mathrm{o}}^{n+1}{\text{+ Q}}_{\mathrm{g}}^{n+1}\\ =& \frac{V_{ij}}{\mathrm{\Delta }t}[(\varphi {\rho }_{\mathrm{o}}{s}_{\mathrm{o}}{)}^{n+1}-(\varphi {\rho }_{\mathrm{o}}{s}_{\mathrm{o}}{)}^n+(\varphi {\rho }_{\mathrm{o}}{s}_{\mathrm{o}}{R}_{\mathrm{s}\mathrm{o}}+\varphi {\rho }_{\mathrm{g}}{s}_{\mathrm{g}}{)}^{n+1}\\ & -(\varphi {\rho }_{\mathrm{o}}{s}_{\mathrm{o}}{R}_{\mathrm{s}\mathrm{o}}+\varphi {\rho }_{\mathrm{g}}{s}_{\mathrm{g}}{)}^n].\end{array}

The equation is organized by moving the flow rate and starting pressure gradient terms to the right side of the equation:

(68)\begin{array}{rl}& \mathrm{\Delta }(T_{\mathrm{o}}^n\mathrm{\Delta }{p}_{\mathrm{o}}^{n+1})\text{+}\mathrm{\Delta }[(T_{\mathrm{o}}^n{R}_{\mathrm{s}\mathrm{o}}\mathrm{\Delta }{p}_{\mathrm{o}}^{n+1}\text{+}{T}_{\mathrm{g}}^n(\mathrm{\Delta }{p}_{\mathrm{o}}^{n+1}\text{+}\mathrm{\Delta }{p}_{\mathrm{c}\mathrm{g}\mathrm{o}}^n\text{)]}\\ =& \frac{V_{ij}}{\mathrm{\Delta }t}[(\varphi {\rho }_{\mathrm{o}}{s}_{\mathrm{o}}{)}^{n+1}-(\varphi {\rho }_{\mathrm{o}}{s}_{\mathrm{o}}{)}^n\\ & +(\varphi {\rho }_{\mathrm{o}}{s}_{\mathrm{o}}{R}_{\mathrm{s}\mathrm{o}}+\varphi {\rho }_{\mathrm{g}}{s}_{\mathrm{g}}{)}^{n+1}-(\varphi {\rho }_{\mathrm{o}}{s}_{\mathrm{o}}{R}_{\mathrm{s}\mathrm{o}}+\varphi {\rho }_{\mathrm{g}}{s}_{\mathrm{g}}{)}^n]\\ & -(R_{\mathrm{o}}^n\text{+}{R}_{\mathrm{o}\mathrm{g}}^n{\text{+ Q}}_{\mathrm{o}}^{n+1}{\text{+ Q}}_{\mathrm{g}}^{n+1}).\end{array}

The pressure equation is linearized. The system of equations can be solved using Newton’s iterative method to find the pressure iteration value (p_{\mathrm{o}}^{n+1}) at the n + 1 moment. From the capillary pressure auxiliary equation p_{\mathrm{o}}=p_{\mathrm{g}}- p_{\mathrm{c}\mathrm{g}\mathrm{o}}, p_{\mathrm{g}}^{n+1} can be found. The resulting p_{\mathrm{g}}^{n+1} is substituted into the differential equation for gas phase percolation, and the s_{\mathrm{g}}^{n+1} explicitly is calculated.

The differential equation for gas saturation is

(69)\begin{array}{rl}& s_{\mathrm{g}}^{n+1}=\frac{\mathrm{\Delta }t}{V_{ij}\varphi {\rho }_{\mathrm{g}}}\{\mathrm{\Delta }[(T_{\mathrm{o}}^n{R}_{\mathrm{s}\mathrm{o}}\mathrm{\Delta }{p}_{\mathrm{o}}^{n+1}\text{+}{T}_{\mathrm{g}}^n(\mathrm{\Delta }{p}_{\mathrm{o}}^{n+1}\text{+}\mathrm{\Delta }{p}_{\mathrm{c}\mathrm{g}\mathrm{o}}^n\text{)]}\\ & +R_{\mathrm{o}\mathrm{g}}^n{\text{+ Q}}_{\mathrm{g}}^{n+1}\}-[(\varphi {\rho }_{\mathrm{o}}{s}_{\mathrm{o}}{R}_{\mathrm{s}\mathrm{o}}{)}^{n+1}-(\varphi {\rho }_{\mathrm{o}}{s}_{\mathrm{o}}{R}_{\mathrm{s}\mathrm{o}}{)}^n]\\ & +(\varphi {\rho }_{\mathrm{g}}{s}_{\mathrm{g}}{)}^n.\end{array}

After finding s_{\mathrm{g}}^{n+1}, s_{\mathrm{o}}^{n+1} can be calculated from the auxiliary equation for saturation s_{\mathrm{o}}=1-s_{\mathrm{g}}.

The IMPES solution method is conditionally stable with the stability condition:

(70)max\left|\mathrm{\Delta }t\left({\lambda }_{i+\frac{1}{2}}\frac{p_{i+1,j}^{n+1}-p_{i,j}^{n+1}}{\mathrm{\Delta }{x}_{i+\frac{1}{2},j}}+{\lambda }_{j+\frac{1}{2}}\frac{p_{i,j+1}^{n+1}-p_{i,j}^{n+1}}{\mathrm{\Delta }{y}_{i,j+\frac{1}{2}}}\right)\right|⩽\epsilon .

4 Parameter and boundary condition processing

4.1 Parameter processing

1) The rate of absolute permeability

Absolute permeability, as a function of spatial coordinates, is calculated as a summed average:

(71)k_{xi\pm \frac{1}{2},j}^n=\frac{\mathrm{\Delta }{x}_{i\pm 1}+\mathrm{\Delta }{x}_i}{{\displaystyle \frac{\mathrm{\Delta }{x}_{i\pm 1}}{k_{xi\pm 1,j}^n}}+{\displaystyle \frac{\mathrm{\Delta }{x}_i}{k_{xi,j}^n}}},

(72)k_{yi\pm \frac{1}{2},j}^n=\frac{\mathrm{\Delta }{y}_{i\pm 1}+\mathrm{\Delta }{y}_i}{{\displaystyle \frac{\mathrm{\Delta }{y}_{i\pm 1}}{k_{yi\pm 1,j}^n}}+{\displaystyle \frac{\mathrm{\Delta }{y}_i}{k_{yi,j}^n}}}.

2) The flow coefficient

In the flow coefficient ({\lambda }_{\mathrm{l}}={\displaystyle \frac{k{k}_{r\mathrm{l}}{\rho }_{\mathrm{l}}}{{\mu }_{\mathrm{l}}}} (l = o, g)), the absolute permeability (k) is taken as the summed average of two adjacent meshes. Then, {\displaystyle \frac{k_{rl}{\rho }_l}{{\mu }_l}} is weighted using a single upstream point of processing:

(73)

(74)

3) Viscosity treatment of oil and gas mixtures

The oil phase viscosity (µ_o) in the oil phase percolation equation was constant. The viscosity of the oil-gas mixture (µ_mix) in the gas phase percolation equation was modified by Eqs. (29)–(31), and the same single point upstream weighting was used.

4) Initiate pressure gradient processing

The change in viscosity of the oil-gas mixture was considered after the dissolution of CO₂ in the oil phase, leading to the variation in the oil-gas mixture starting pressure gradient. The oil-gas mixture starting pressure gradient was calculated according to Eq. (32). A single point upstream weighting was also used to handle the value.

4.2 Boundary condition processing

1) Outer boundary conditions

For the closed outer boundary condition, a virtual mesh was created outside the closed boundary. The pressure of the boundary mesh was equal to the virtual mesh pressure:

(75)\{\begin{array}{c}{p}_{1,j}^{n+1}=p_{2,j}^{n+1}\\ p_{i,1}^{n+1}=p_{i,2}^{n+1}\end{array},

(76)\{\begin{array}{c}{p}_{n_x,j}^{n+1}=p_{n_x-1,j}^{n+1}\\ p_{i,n_y}^{n+1}=p_{i,n_y-1}^{n+1}\end{array}.

2) Internal boundary conditions

The injection and production wells in a reservoir were special mesh cells containing source and sink terms, with negative production from oil producing wells and positive injection from gas injection wells. The mesh (i, j) had a well with a volume flow rate (Q_v) and a fixed pressure of the bottom flow of the production well. Therefore, the volume flow rate (Q_v) needed to be expressed in terms of mesh pressure (p_ij) and bottom flow pressure (p_wf).

The proposed steady-state equation for the planar radial flow of oil production from a producing well was

(77)Q_{\mathrm{v}\mathrm{o}}=-PI D{\lambda }_{\mathrm{o}}[p_{\mathrm{o}i,j}-p_{\mathrm{w}\mathrm{f}}- G_{\mathrm{o}}(r_{\mathrm{e}} -r_{\mathrm{w}} )].

The proposed steady-state equation for the planar radial flow of gas injection volume in gas injection wells was

(78)Q_{\mathrm{v}\mathrm{g}}=-WI D{\lambda }_{\mathrm{g}}[p_{\mathrm{g}i,j}-p_{\mathrm{w}\mathrm{f}}+G_{\mathrm{g}}(r_{\mathrm{e}} -r_{\mathrm{w}} )].

5 Results and discussion

An independent programming implementation of the numerical simulator was presented on a microcomputer for CO₂ immiscible flooding simulation by coupling the reservoir simulation model with corrected models of viscosity and starting pressure gradient. The numerical simulator was compiled in the C++ language. The program framework flow chart is shown in Fig.3.

The example calculation in the numerical simulator using the ultra-low permeability reservoir parameters in F block of S oilfield made the calculation more consistent with the actual situation of the field. It could verify the correctness of the developed seepage mathematical model and numerical solution under the condition of variations in crude oil viscosity and oil phase initiation pressure gradient. Also, it provided a solid foundation and reliable data support for the future development of the mathematical model and numerical solution for CO₂ immiscible flooding.

The parameters of block F were described as follows: The reservoir depth was 2800–3400 m, the reservoir temperature was 142°C, the average formation pressure was 26 MPa, crude oil saturation was 72%, non-reducible water saturation was 28%, residual oil saturation was 10%, the formation oil density was 0.76 g·cm⁻³, oil viscosity was 1.5 mPa·s, CO₂ viscosity was 0.05 mPa·s, CO₂ compressibility was 0.7, reservoir permeability was 5 × 10⁻³μm², and porosity was 12.5%.

The working system represented a constant volume flow rate of injection at the inner boundary and constant bottom flow pressure production. It was closed at the outer boundary. In the numerical simulation of the oil reservoir process, the time step was one day, and the space was 20 × 20 mesh. Also, the space step in the X direction was 20 m, the space step in the Y direction was 10 m, and the total displacement time was t = 1000 d.

5.1 The oil pressure field, gas saturation field, and oil viscosity field with the change in oil viscosity and starting pressure gradient

The viscosity of crude oil and the start-up pressure gradient was considered. Therefore, CO₂ immiscible flooding changed the physical parameters of the formation fluid in real time at different times. In this paper, the coupled model of CO₂ immiscible flooding caused the accurate simulation of real-time changes in the different variables of the seepage field. Fig.4–Fig.5 show the profiles of the oil pressure field, gas saturation field, and oil viscosity field with time nodes of 10, 100, 200, 400, 600, 800, and 1000 days.

From Fig.4, the oil pressure rose near the injection well. The original formation pressure remained unchanged across the rest of the grid. The amount of injected CO₂ and the formation pressure increased during the time, which was transferred from the injection well to the production well. The pressure grid gradually rose, and the color variation indicated that the pressure was gently increasing. The formation pressure changed along with the flow process of CO₂ flooding in time nodes of 10, 100, 200, 400, 600, 800, and 1000 days. The oil pressure decreased step by step from the gas injection well to the production well.

Fig.6 shows the variation in the CO₂ saturation field. Near the gas injection well, CO₂ saturation increased significantly at the 10-day node. After the gradual injection of CO₂, the leading edge of CO₂ immiscible flooding advanced toward the production well, and the value of gas saturation rose considerably. At the 1000-day node, the leading edge of CO₂ immiscible flooding basically reached the production well, and CO₂ appeared in the production well around the same time. The degree of change in carbon dioxide saturation was the largest near the gas injection well, and the related graph rose steeply. The closer to the oil production well, the smaller the degree of variation in carbon dioxide saturation.

Fig.5 represents the change in the viscosity field of crude oil with CO₂ injection. At the beginning, the viscosity of the formation crude oil was a constant value. The viscosity of the crude oil near the gas injection well was significantly reduced for about 10 days due to the dissolution and viscosity reduction in the crude oil by CO₂. Around 200 days, the graph formed a “deep pit” because of the viscosity decrease in crude oil in the affected area by CO₂ propulsion. The whole “viscosity plateau” turned into a “viscosity trough” during the time. The viscosity of crude oil gradually increased from gas injection well to production well. It was due to the gradual reduction in CO₂ gas saturation.

5.2 Oil pressure field, oil saturation field, and oil viscosity field with and without considering the changes in crude oil viscosity and starting pressure gradient

For a better comparison between the amounts of the oil pressure field, oil saturation field, and oil viscosity field with and without considering the change in crude oil viscosity and starting pressure gradient, the calculation process took a time step of one day, a spatial mesh of 10 × 10, a spatial step of 50 m in both X and Y directions, and a total repulsion time of t = 1000 d. Also, a plane contour map was used for display.

Fig.7 shows the distribution of the pressure field. After the gas injection process, the formation pressure rose gradually. The formation pressure decreased gently from the gas injection wells to the production wells, and the gradient of pressure reduction was the maximum in the direction of the injection and production wells. Therefore, the pressure field distribution represented the densest contour in the direction of the injection and production wells. The pressure at each point in the pressure field with considering the change in oil phase viscosity was lower than that without considering the variation in oil phase viscosity. As the viscosity of crude oil declined continuously due to the dissolution of CO₂ in crude oil, the seepage resistance of the formation fluid decreased, and it was easier to establish an effective displacement pressure gradient. From the pressure field distribution at 900 days, the pressure propagated gradually from gas injection wells to production wells. Also, the pressure distribution in the whole pressure field gently increased.

Fig.8 represents the distribution of the oil saturation field. The saturation of crude oil at each point in the oil saturation field with considering the change in crude oil viscosity was lower than that without considering the variation in crude oil viscosity. Therefore, the dissolved CO₂ in crude oil reduced the crude oil viscosity, and the crude oil in the rock pores was easily driven out. It indicated that the crude oil viscosity decreased due to the dissolution of CO₂ in the crude oil, and the crude oil in the rock pore space was easily swept out. The crude oil saturation was lower near the gas injection wells because the crude oil viscosity decreased more and the oil drive efficiency increased. However, the crude oil saturation declined more slowly near the production wells. Regarding the distribution of crude oil saturation at t = 900 d, 1) the saturation of crude oil near the gas injection wells was the lowest, 2) the saturation of crude oil near the production wells was the highest, 3) the area affected by the gas was getting larger and larger, and 4) the gas drive front was gradually advancing from the gas injection wells to the production wells.

Fig.9 shows the distribution of the oil phase viscosity field. The distribution of the viscosity field represented that the viscosity of crude oil in the formation gradually increased from gas injection wells to production wells. Compared with the initial formation crude oil viscosity, the viscosity of crude oil in the formation decreased significantly to 32% of the initial formation crude oil viscosity near the gas injection wells and 25% of the initial formation crude oil viscosity near the production wells. At 300 and 900 days, the viscosity field distribution changed greatly. The gradient of crude oil viscosity decrease was larger near the gas injection wells at t = 300 d. The gradient of crude oil viscosity reduction was higher near the production wells at t = 900 d. However, the crude oil viscosity in most of the area between the gas injection wells and production wells decreased to 0.36 mPa∙s. The viscosity reduction in CO₂ on the formation of crude oil was significant. According to the graph at t = 300 d, the contour line near the injection well was dense and its value was small, but the contour line near the production well was sparse and its value was large. At t = 900 d, the contour line near the injection well began to become sparse, and the value was smaller because the viscosity of crude oil was significantly reduced due to the dissolution of carbon dioxide in crude oil.

Fig.10(a) represents that the overall formation pressure and the pressure propagation velocity with considering the oil phase viscosity and the starting pressure gradient were lower than those without considering the starting pressure gradient.

Fig.10(b) shows that the oil saturation during the formation when the initiation pressure gradient was considered was higher than when the initiation pressure gradient was not considered because it was more difficult to establish the effective replacement pressure when the initiation pressure gradient was considered, and the required replacement pressure was larger, while the effective replacement pressure was smaller, so the crude oil swept out of the rock pore space was less, and the remaining crude oil was more, and the oil saturation after the replacement was higher. The saturation of oil content was also higher.

Fig.10(c) represents that the viscosity of crude oil gradually increased from gas injection wells to production wells; because the high CO₂ concentration near gas injection wells and a large amount of dissolved CO₂ in crude oil caused a significant reduction in crude oil viscosity. However, the low CO₂ concentration near production wells and the small amount of dissolved CO₂ in crude oil led to a trivial decrease in crude oil viscosity. Therefore, the viscosity of crude oil near the production well was higher than that of the injection well. The comparisons obviously showed that the overall viscosity of crude oil in the formation considering the initiation pressure gradient was higher than that without the initiation pressure gradient because the formation pressure was low when the initiation pressure gradient was considered. The effect of CO₂ on viscosity reduction in crude oil was related to the formation pressure, formation temperature, and formation crude oil density. It was most affected by the formation pressure. Therefore, the lower formation pressure led to a weaker reduction in the viscosity of crude oil CO₂, which caused the higher viscosity of formation crude oil compared with that without considering the initiation pressure gradient.

The starting pressure gradient was a function of fluid flow, and the higher the fluid viscosity, the larger the starting pressure gradient. As Fig.10(d) represents, the viscosity of crude oil gradually increased from gas injection wells to production wells. Therefore, the starting pressure gradient gradually increased from gas injection wells to production wells. The closer the gas injection well, the better the effect of CO₂ viscosity reduction on the crude oil. The smaller the starting pressure gradient, the lower the fluid percolation resistance, and the easier the establishment of an effective replacement pressure system.

From Fig.11, the pressure at the (5, 5) mesh in the formation increased gradually with the gas drive process. However, the pressure value when the initiation pressure gradient was considered was smaller than that without considering the initiation pressure gradient. The pressure at a certain point in the seepage field was lower than that without the starting pressure gradient because more pressure had to be consumed to drive the fluid in the rock pore space when the starting pressure gradient was considered. However, as the replacement process continued, the pressure and replacement energy were continuously transferred from the gas injection wells to the production wells. Accordingly, the pressure at one point in the seepage field was still increasing. Therefore, the pressure on the whole formation was increasing.

From Fig.12, the oil saturation at mesh (5, 5) preserved the original oil saturation of 0.98 around t = 500 d, which meant that the front edge of the drive did not reach this mesh before that. After t > 500 d, the front edge of the gas drive reached this mesh, and the crude oil in the rock pore was continuously swept out, and the oil saturation started to decrease. When the gas drive front continued to advance, the oil saturation at the (5, 5) mesh decreased to about 0.72 at t = 1000 d. The oil saturation considering the starting pressure gradient was slightly higher than that of the conditions without considering the starting pressure gradient.

According to Fig.13, the viscosity of the oil phase at (5, 5) mesh gradually decreased with the gas drive process because the gas saturation gradually increased at this mesh, the oil saturation gently declined, the CO₂ concentration gradually rose, and the amount of dissolved CO₂ in the crude oil gradually increased. Therefore, the crude oil viscosity showed a decreasing trend. Analyzing the start-up pressure gradient, the sweep speed of the driving pressure and the advance speed of the carbon dioxide displacement front were slow. Accordingly, at the same time, the red line in Fig.13 was higher than the black line. It meant that the crude oil viscosity in this grid with the start-up pressure gradient was higher than that without it.

From Fig.14, the production well yield gradually increased with time, and the yield at the initial stage before t = 300 d was similarly higher than that without the starting pressure gradient. However, at the middle and late stages at t > 300 d, the yield with the starting pressure gradient was significantly lower than that without the starting pressure gradient, and the yield gap became larger during the time. The existence of the start-up pressure gradient led to the slow advance speed of the gas drive front and low production, but a longer stable production period. Fig.14 represents that the production curve considering the start-up pressure gradient was gentler than after 300 days.

6 Conclusions

Based on the theoretical study of fluid physical properties in reservoir seepage systems, it was concluded that the fluid physical properties changed during the development process. Also, the variation in fluid physical properties profoundly affected the seepage process. In the process of CO₂ flooding seepage, CO₂ dissolved with crude oil, extraction, mass transfer, and even miscible effect reduced the viscosity of crude oil. Simultaneously, the CO₂ gas saturation change caused a variation in the oil phase relative permeability. Therefore, the oil phase change starting pressure gradient and fluid filtration resistance variation made the flooding seepage calculation more complex than the water flooding. The following conclusions were obtained for CO₂ immiscible flooding.

1) The formation pressure and formation oil saturation were low in the conditions considering the oil phase viscosity in comparison with the situation without taking into account the oil phase viscosity. The oil phase viscosity gradually increased from the gas injection well to the production well.

2) The results where oil phase viscosity was considered and the starting pressure gradient was not considered indicated the low values for formation pressure and the crude oil production, but the high amounts for the formation oil saturation and the formation oil phase viscosity. In these conditions, the starting pressure gradient gradually increased from the gas injection well to the production well.

3) Two ideas were adopted to analyze the pressure, saturation, viscosity, and other parameters of the seepage area at different positions at the same time, and at the same position at different times. It resulted in a clearer understanding of the oil displacement law.

In this study, the preliminary novel findings were presented and analyzed to simulate the CO₂ immiscible flooding for an ultra-low permeability reservoir. More research and improvements remain to be performed for future practical projects. A major direction of the next works is to extend the framework, handling more complicated physical characteristics such as compositional flow. Furthermore, more properties with different fluid compositions will be evaluated. In our future study, the influences of different variables on the seepage process in CO₂ immiscible flooding and miscible flooding will be taken into account. These factors include the oil phase viscosity, oil-gas-water starting pressure gradient, pressure sensitivity, and other variables.

References

[1] AllenR, SunS (2012). Carbon Dioxide Sequestration: Modeling the Diffusive and Convective Transport under a CO₂ Cap. In: SPE Saudi Arabia Section Technical Symposium and Exhibition, AlKhobar, Saudi Arabia, 8–11 April. SPE-160881-MS..

[2] An W Q, Yue X, Feng X, Fu J Y, Fang X, Zou J, Fang W (2017). Non-Klinkenberg slippage phenomenon at high pressure for tight core floods using a novel high pressure gas permeability measurement system.J Petrol Sci Eng, 156: 62–66

[3] Battiatoa, I., Tartakovskya, D. M., Tartakovsky, A. M.(2011). Hybrid models of reactive transport in porous and fractured media. Adv Water Resour. 34 (9): 1140–1150.

[4] Lv G Y, Wang J, Sun Z G (2002). An experimental study on starting pressure gradient of fluids flow in low permeability sandstone porous media.Pet Explor Dev, 29(2): 86–89

[5] Chi J B Lyu G, Zhang X, Wang J (2017). A computational method of critical well spacing of CO₂ miscible and immiscible concurrent flooding.Pet Explor Dev, 44(5): 815–823

[6] ChungF T H, JonesR A, HaiT N (1988) Measurements and correlations of the physical properties of CO₂ heavy crude oil mixtures. SPE 15080

[7] Civan F (2010). Effective correlation of apparent gas permeability in tight porous media.Transp Porous Media, 82(2): 375–384

[8] CuiC Z, YAND W, YAOT Y, Wang J, ZHANGC B, WUZ W (2022). Migration law of CO₂ flooding front and prediction method of gas channeling: a case study of G89–1 block in Shengli Oilfield. Oil Gas Reservoir Evaluation Develop, 12(05): 741–747+763

[9] Davarpanah A (2020). Parametric study of polymer-nanoparticles-assisted injectivity performance for axisymmetric two-phase flow in EOR processes.Nanomaterials (Basel), 10(9): 1818

[10] Firoozabadi A, Myint P C (2010). Prospects for subsurface CO₂ sequestration.AIChE J, 56(6): 1398–1405

[11] HaungY Z (1998). Percolating Flow Mechanism of Low Permeability Reservoir. Beijing: Petroleum Industry Press: 80–86

[12] IPCC (Intergovernmental Panel on Climate Change) (2005). IPCC Special Report on Carbon Dioxide Capture and Storage. New York City: Cambridge University Press

[13] Javadpour F (2009). Nanopores and apparent permeability of gas flow in mudrocks (shales and siltstone).J Can Pet Technol, 48(08): 16–21

[14] JuB S, YuJ B, LvG Z, Cao W D (2020). CO₂ flooding of low permeability reservoir numerical simulation method and application. Geol Oil Recovery Oil gas, 27 (01): 126–133

[15] Li D L, Zha W, Liu S, Wang L, Lu D (2016). Pressure transient analysis of low-permeability reservoir with pseudo threshold pressure gradient.J Petrol Sci Eng, 147: 308–316

[16] Lu J, Dai F, Rahman M M, Escobar F H (2017). Boundary dominated flow in low permeability reservoir with threshold pressure gradient.RPN J Eng Appl Sci, 12(23): 6834–6843

[17] Michael, K., Golab, A., Shulakova, V (2011). Geological storage of CO₂ in saline aquifers—a review of the experience from existing storage operations. Int J Greenhouse Gas Contr, 4(4): 659–667.

[18] Pertsin A, Grunze M (2004). Water-graphite interaction and behavior of water near the graphite surfaced.J Phys Chem B, 108(4): 1357–1364

[19] Gao R , Lv C Y, Lun Zengmin , Wang R (2021). Numerical simulation of carbon dioxide displacement and storage.Special Oil Gas Reserv, 28(02): 102–107

[20] Wang X K, Sheng J J (2017a). Effect of low-velocity non-Darcy flow on well production performance in shale and tight oil reservoirs.Fuel, 190: 41–46

[21] WangY (2021). Experimental study of nonlinear seepage law in tight oil reservoir and its application in digital modeling software. The dissertation for Doctoral Degree. Anhui: University of Science and Technology of China.

[22] ZhangJ G, LeiG L, ZhangY Y (1998). The Oil and Gas Flow Through Porous Media. Dongying: China University of Petroleum Press

RIGHTS & PERMISSIONS

Higher Education Press