Analytical solution of Buckley-Leverett equation for gas flooding including the effect of miscibility with constant-pressure boundary

Abstract

This paper presents an analytical solution of Buckley-Leverett equation for gas flooding with constant-pressure boundary including the effect of miscibility on the viscosity and relative permeability. First, a relative permeability model and a viscosity model with consideration of miscibility are used to describe the variations of relative permeability and viscosity of oil and gas. Then, based on the fractional-flow theory, the Buckley-Leverett equation for gas flooding with constant-pressure boundary including the effect of miscibility is constructed and solved analytically. From the analytical solution, the saturation and pressure profiles, the total volumetric flux and the breakthrough time are determined. To verify the theory, the analytical solution is compared with the numerical solution. The comparison shows that the analytical solution is in reasonable agreement with numerical solution. Through the study on the influential factors, it can be concluded that total volumetric flux is increasing with the increases of permeability and pressure and decrease of gas viscosity. The increase of total volumetric flux accelerates breakthrough of the injected gas. Furthermore, with the pressure increase, there are remarkable reduction in residual oil saturation and improvement of relative permeability, resulting in higher gas saturation and oil displacement efficiency. The analytical solution presented in this paper provides guidance on analyzing the distribution of saturation and pressure profiles, predicting the gas production and oil recovery efficiency of oil well.

Introduction

For its remarkable achievement in oil recovery and storage capacity for the greenhouse gases, the gas flooding technology has been attracting more and more concern and discussion among the EOR methods (Kamali and Cinar, 2014; Zhao et al., 2016). The gas injection is divided into two categories, consisting of hydrocarbon gas injection, such as natural gas and produced gas, and non-hydrocarbon gas injection, such as CO₂, N₂ and flue gas. The injected gas not only dissolves into oil leading oil viscosity reduction, volumetric expansion, and interfacial tension (IFT) reduction but also plays a role in dissolved gas drive (Green and Willhite, 1998; Holm, 1986). The contribution of these mechanisms is different for oil recovery under different reservoir and fluid condition (Klins, 1984). The primary mechanism of gas injection is through mass transfer to reduce IFT and increase the miscibility of oil and gas, especially under high pressure (Sheng, 2013), which allows for the residual oil in water flooding to be displaced and recovered.

The theory of fractional flow reveals the basic law of non-piston displacement for two-phase flow, based on some indispensable and reasonable assumptions. The theory is generalized, starting with the Buckley-Leverett theory for water flooding and has been applied to various EOR methods, such as gas flooding and polymer flooding (Buckley and Leverett, 1942; Ehlig-Economides and Economides, 2010; Ghanbarnezhad and Lake, 2012; Hirasaki, 1981; Johansen et al., 2016; Pope, 1980). The theory involves one-dimensional flow of two incompressible fluids in porous media. To easily obtain the analytical solution, necessary and reasonable assumptions are made. Typically, the major assumptions are following: (1) the porous media is homogeneous, isotropic, and isothermal, (2) dispersion of fluids is negligible, (3) gravity and capillarity are negligible, (4) Darcy’s law is applied, (5) and a continuous injection of constant composition is injected starting at time zero (Pope, 1980; Rui et al., 2018). The fractional flow for the displacing fluid is the ratio of the volume flux of the displacing fluid to the total volume flux; by substituting the Darcy’s Law into the volume flux, the function of the fractional flow is determined. By substituting the function of the fractional flow into the conservation equation, the new equation is known as Buckley–Leverett equation. Afterwards, the front advance velocity, the breakthrough time and the pressure and saturation profiles, indicators of sweep efficiency and oil recovery are deduced with boundary conditions. Yortsos et al. presented a mathematical model of Buckley-Leverett equation for water injection and obtained its analytical and numerical solution (Yortsos et al., 1999). The log-log plot of water/oil ratio and time was used to predict the well pattern and the relative permeability characteristics. Despite the simplification of the mathematical model, the theory retains the key feature of the two-phase flow and provides a highly valuable fundamental understanding of the displacement mechanisms. The theory has the best reservoir engineering principles and practices and is efficient in precisely investigating the two-phase flow problems.

Many researchers have noticed the advantage of the theory of fractional flow in gas flooding. Welge proposed a method to predict oil recovery for displacement of oil by enriched gas in a three-component, two-phase system (Welge et al., 1961). The method considered the mechanisms of wet gas displacing oil, including volumetric expansion, viscosity reduction and condensation of intermediate components in the gas. Wu et al. constructed a one-dimensional non-Darcy model of immiscible fluids, in which non-Darcy flow was described with the Barree and Conway model (Wu et al., 2011). They derived its analytical solution and compared the Barree and Conway model with the Forchheimer model. The results showed that the non-Darcy flow was dominated by the non-Darcy relations and the parameters of the injection and production. Laforce and Johns (2005) established and solved a three-phase four-component flow model in a one-dimensional porous medium by the method of characteristics. It concluded that the presence of the three-phase region added the complexity of the types of composition routes and profoundly affected on oil recovery, which decreased with displacing-fluid enrichment. Johansen et al. (2005) presented a comprehensive analysis of mathematical model of four-component gas/oil displacement systems in which components partition between two phases is governed by the geometry of a set of equilibrium tie lines. The paper pointed out that when a system contained global triangular structure, the wave analysis indicated that the shock and rarefaction surfaces coincided. LaForce extended the parameterization and presented an analytical solution for flow of four components in three phases with a realistic phase behavior model. The study of displacements at pressures both below and above the minimum miscibility pressure (MMP) revealed that solution structure for these displacements is qualitatively similar to that of the constant K-value system (LaForce and Orr, 2009; Rezaei et al., 2013). Yang (2014) constructed a mathematical model for carbonated water injection under constant pressure boundaries and developed its numerical solution using the IMPES method. Through the model, they predicted the oil recovery, the saturation distribution and the pressure distribution from injector to producer.

Unfortunately, most of these studies have focused on the problems involved component and phase, from two phases to three phases and from three components to four components. However, the changed physical properties of fluids during the gas injection have been overlooked in the development of analytical solution of Buckley-Leverett equation. Actually, as the recovery mechanism stated previously, since the gas is injected into the porous media, the complex interaction of oil and gas will lead to changes of physical properties such as viscosity and relative permeability, the key characteristics of multiphase seepage. In addition, most previous studies were based on the assumption of constant viscosity and relative permeability. Meanwhile, it is noticed that the classical fractional-flow theory is based on the assumption of constant volumetric flux. However, it is also common that the injection pressure and the production pressure are both kept constant over time (Cui et al., 2018b; Johansen et al., 2016). Under the constant-pressure boundary conditions, the volumetric flux is not constant and the volumetric flux at different time is unknown. Although the one-dimensional two-phase flow model with constant-pressure boundary and dynamic properties can be solved by the numerical method, the numerical simulation is time-consuming, and there are numerical errors in the numerical discretization which may result in physically unreal oscillation and wrong evaluation of oil reservoirs (Aziz, 1979; He et al., 2017, 2018). Thus, in order to ensure the accuracy and improve the computational efficiency, an analytical solution is indispensable (He et al., 2018).

The purpose of this paper is to develop an analytical solution of the Buckley-Leverett equation for gas flooding with constant-pressure boundary conditions considering the changed viscosity and the changed relative permeability. To describe the variation of relative permeability, we improve the analytical model proposed by Coats based on the mechanism of the gas flooding. The relative permeability under different miscible condition is determined by the interpolation method. Moreover, we introduce the Todd-Longstaff model to calculate the effective viscosities of the oil and gas. The analytical solution presented in this paper provides the saturation and pressure profile, the total volumetric flux and the breakthrough time. Furthermore, the influence of permeability, gas viscosity and pressure on process of gas flooding is fully discussed through the analytical results.

Mathematical model

Before the construction of the mathematical model of Buckley-Leverett equation, firstly we study the varying physical properties of fluids. It has been pointed out that the relative permeability and viscosity will change in the process of gas injection. Thus, it is necessary to describe this variation with a proper model. There are increasingly sophisticated or efficient methods to model the relative permeability and viscosity of fluids in the gas flooding. To obtain the analytical solution of Buckley-Leverett equation, simple and practical models with proper parameters are needed.

Relative permeability of gas/oil system

Relative permeability is one of the fundamental parameters and key characteristics of multiphase seepage. In the process of CO₂ flooding, the physical and chemical properties of the fluids change drastically, especially at the interface (Cui et al., 2018a; Orr, 2007). Under such condition, the influence of the thermodynamic condition on the relative permeability is by no means negligible. IFT, generally through the capillary force, plays an important role in determining the flow behavior of hydrocarbon fluids in porous rocks.

Some scholars have carried out researches from a mathematical point of view to predict relative permeability at different pressure. All of these models adopt an interpolation function to interpolate between base curves (at low pressure where the injected gas and the crude oil are immiscible) and straight lines representing the miscible condition at high pressure (generally higher than MMP). Coats first proposed the model of this type for individual-phase relative permeability with the interpolation which is the function of IFT (Coats, 1980). The model employs linear interpolation to modify the relative permeability curve at different values of IFT, which is shown as follows

$$K_{\mathit{ro}}=F_k\cdot K_{\mathit{ro}}^{\mathit{imm}}+\left(1-F_k\right)\cdot K_{\mathit{ro}}^{\mathit{mis}}$$

(1)

$$K_{\mathit{rg}}=F_k\cdot K_{\mathit{rg}}^{\mathit{imm}}+\left(1-F_k\right)\cdot K_{\mathit{rg}}^{\mathit{mis}}$$

(2)

$$S_{\mathit{or}}=F_k\cdot S_{\mathit{or}}^{\mathit{imm}}$$

(3)

$$S_{\mathit{gi}}=F_k\cdot S_{\mathit{gi}}^{\mathit{imm}}$$

(4)

F_k is the relative permeability interpolation parameter, given by the following formula

$$F_K=\min \left[1,{\left(\frac{\sigma }{{\sigma }_0}\right)}^N\right]$$

(5)

where σ is the IFT, and σ₀ is the “initial” tension corresponding to the read-in capillary pressure curve. $K_{\mathit{ro}}^{\mathit{imm}}$ and $K_{\mathit{rg}}^{\mathit{imm}}$ are immiscible relative permeabilities. $K_{\mathit{ro}}^{\mathit{mis}}$ and $K_{\mathit{rg}}^{\mathit{mis}}$ are miscible relative permeabilities generally taken to be a linear function of saturation.

The model further revised by many researchers and widely used in commercial simulators, such as Eclipse and CMG (Al-Nuaimi et al., 2018; Coats, 1980; GeoQuest, 2014; Hustad and Browning, 2010). This is a simple and practical approach to evaluate the effect of miscibility on the relative permeability. Whitson et al., Pope et al. and Blom et al. proposed mathematical models in which the relative permeability varied with capillary number rather than IFT (Blom and Hagoort, 1998; Pope et al., 2000; Whitson and Fevang, 1997). Laroche developed a 1D, horizontal, immiscible two-phase flow model taking the capillary end effect into consideration. Through the model, the exponent of the relative permeability near its residual value was determined from a segment of the curve of the outlet fraction flow (Laroche et al., 2001). Abbas (2016) improved JBN method for calculating relative permeability using the cubic spline numerical modeling technique. Compared with the traditional way, the technique showed good results and the overall numerical error was reduced. The interpolation parameter is the function of IFT in these models (Al-Nuaimi et al., 2018; Alzayer et al., 2018; Hustad and Browning, 2010; Mott et al., 2000). Mott et al. and Al-Nuaimi et al. established empirical formulas to determine the interpolation parameter by fitting the experimental data (Al-Nuaimi et al., 2018; Mott et al., 2000). The limitation of the method is that it needs to perform experiments which are time-consuming and costly. Alzayer et al. (2018) adopted the Gibbs free energy to determine the interpolation parameter, which required a complex calculation process. Coats pointed out that the treatment given here was not based on any theory or experimental evidence. In conclusion, on one hand, it is necessary to prove the correctness of the interpolation method; on the other hand, an effective and simple method to determine the interpolation parameter is needed.

Many previous studies have addressed issues of gas/oil relative permeability as a function of IFT through laboratory experiments. Bardon and Longeron measured gas and oil curves of methane/heptane system under different IFT using unsteady method (Longeron, 1980). IFT was controlled by varying the equilibrium pressure of the mixture. The experimental results showed that the shape of the curves varied considerably depending on the physical properties of the phases and a linear relation between the relative permeability and IFT, seen in Figure 1(a). Further analysis revealed that the major change in relative permeability seemed to occur in the vicinity of σ = 0.04 mN/m and the relative permeability curves tended towards straight lines as IFT approached 0.001 mN/m. Asar and Handy (1988) did a similar study using a highly volatile methane/propane system to measure steady-state relative permeability as the function of IFT. Similar effect was observed with a critical IFT greater than 0.18 mN/m for a gas-condensate system. Munkerud and Torsaeter used a six-component system and several reference systems to examine the effect of IFT on the relative permeability. A typical figure of relative permeability for different IFT is shown in Figure 1(b), which indicates that the relative permeabilities are almost straight lines when σ = 0.03 and curve more with the increase of IFT. Residual oil saturation and relative permeability determined from the displacement tests are affected strongly by IFT (Blom et al., 2000; Cinar and Orr, 2005; Jahanbakhsh et al., 2016). These experimental studies demonstrate that oil/gas relative permeabilities are affected significantly by IFT, and a linear relationship between the relative permeability and IFT is observed. Especially, when IFT is less than a threshold, the impact of IFT on the relative permeability is more pronounced.

Figure 1.

(a) Relation between IFT and liquid relative permeability for different vapor saturation (Longeron, 1980); and (b) relative permeability for different IFT in a six-component gas condensate system, Berea core (Munkerud and Torsaeter, 1995).

Gas injection is a widely practiced mean for improving oil recovery because of its high displacement efficiency and relatively low cost. More favorable miscibility can be developed by higher pressure, lower temperature, lighter oil and richer gas. CO₂ injection process is generally operated at pressure above the MMP, which is the lowest pressure where the injected gas becomes miscible with crude oil. There are two methods to determine MMP from the viewpoint of miscibility mechanisms: based on the high oil recovery, more specifically, small increment of oil recovery per pressure, such as slim tube test (SL) (Yellig and Metcalfe, 1980), or based on the zeros IFT, such as vanishing-IFT technique (VIT) (Ayirala and Rao, 2011; Rao, 1997) and rising bubble apparatus (Christiansen and Haines, 1987).

However, an inconsistency may exist between them. Figure 2(a) and (b) shows the MMP from SL and VIT separately. It can be seen that the MMP from VIT is higher than that from SL, and data from literatures also reveal similar results, as shown in Figure 3 (Al Riyami and Rao, 2015; Jaubert et al., 2002; Rao and Lee, 2000; Tathed et al., 2010). The important thing to conclude from this analysis is the miscibility criteria from VIT are more stringent. It is generally believed that miscibility may repeatedly develop and break down in a reservoir due to dispersion arising from viscous fingering and reservoir heterogeneity. Thus, many “miscible” processes are in fact a mixture of miscible and low LFT, near-miscible processes (Shyeh-Yung, 1991). The effluent from the SL test consisted of dark in-place oil followed by a red fluid, then an orange fluid, then a yellow fluid, which also demonstrate that the displacement process in the sand packed coil is dynamic. When the reservoir pressure is higher than MMP from VIT test, by contrast, the IFT during miscible displacement process equals to zeros. Consequently, MMP from SL and its corresponding IFT represent the turning point of immiscible and miscible.

Figure 2.

(a) Comparison between MMP from SL and VIT for Rainbow stock-tank oil with propane (Rao, 1997); and (b) comparison between MMP from SL and VIT for Rainbow stock-tank oil with ethane (Rao, 1997).

Figure 3.

MMP from VIT vs. MMP from SL.

According to the above analyses, a methodology to determine the relative permeability interpolation under different miscibility is proposed, considering the complicated interaction between oil and gas. The IFT at the pressure of MMP is defined as critical IFT, σ₀, below which relative permeability change dramatically. As a result, interpolation parameter, F_k, as a function of IFT, is determined to obtain the relative permeability. The Ramey correlation is recommended in the absence of IFT data (Ramey, 1973), seen below. The parameters in the correlation are shown in Appendix 1. The form of Ramey correlation is

$${\sigma }_{\mathit{go}}^{\frac{1}{4}}=P_{\overline{o}}\left[x_{\overline{o}}\left(\frac{{\rho }_o}{M_o}\right)-y_{\overline{o}}\left(\frac{{\rho }_g}{M_g}\right)\right]+P_{\overline{g}}\left[x_{\overline{g}}\left(\frac{{\rho }_o}{M_o}\right)-y_{\overline{g}}\left(\frac{{\rho }_g}{M_g}\right)\right]$$

(6)

The relative permeability of under immiscible condition where $\sigma >{\sigma }_0$ can be obtained by Corey-Brooks model (Brooks and Corey, 1964), shown as followed

$$K_{\mathit{ro}}^{\mathit{imm}}=\mathit{KRO}\cdot {\left(1-S_{\mathit{gn}}\right)}^{n_o}$$

(7a)

$$K_{\mathit{go}}^{\mathit{imm}}=\mathit{KRG}\cdot S_{\mathit{gn}}^{n_g}$$

(7b)

where $S_{\mathit{gn}}\text{=}\frac{S_g-S_{\mathit{gi}}}{1-S_{\mathit{gi}}-S_{\mathit{or}}}$ . n_o and n_g can be obtained by fitting the measured non-miscible relative permeability curve.

Most of the previous studies have revealed that the relative permeability curves tended towards straight lines when miscible condition is developed. That is to say

$$K_{\mathit{ro}}^{\mathit{mis}}=\left(1-S_{\mathit{gn}}\right)$$

(8a)

$$K_{\mathit{go}}^{\mathit{mis}}=S_{\mathit{gn}}$$

(8b)

However, the derivative of relative permeability with respect to gas saturation under such condition is not equal to zero. Without lack of generality, we can assume

$$K_{\mathit{ro}}^{\mathit{mis}}={\left(1-S_{\mathit{gn}}\right)}^{n_m}$$

(9a)

$$K_{\mathit{go}}^{\mathit{mis}}=S_{\mathit{gn}}^{n_m}$$

(9b)

Generally, it is recommended that n_m = 1.1.

Viscosity of oil and gas

CO₂ solubility results in oil viscosity reduction due to gas injection. Zhou and Rui pointed out that with the mix of the heavy oil and CO₂, the viscosity of the heavy oil can be reduced remarkably (up to 98%) (Zhou et al., 2019). Longeron concluded from the displacement experiments that the ratio of the viscosities varied much less, though the capillary number showed very large variation essentially due to the variation in IFT (Longeron, 1980). In this paper, the effective oil and gas viscosities follow the Todd-Longstaff model (Todd and Longstaff, 1972). They are given by

$${\mu }_{\mathit{oeff}}={\mu }_o^{1-\omega }\cdot {\mu }_m^{\omega }$$

(10a)

$${\mu }_{\mathit{geff}}={\mu }_g^{1-\omega }\cdot {\mu }_m^{\omega }$$

(10b)

where ω is the mixing factor of viscosity. Todd believes that the appropriate value of ω to be used for miscible flood simulation is about 1/3 considering the reservoir condition.

The mixture viscosity is given by a 1/4-power fluid mixing rule, shown as follows

$${\left(\frac{1}{{\mu }_m}\right)}^{\frac{1}{4}}=\frac{S_g^{\prime }}{S_n^{\prime }}{\left(\frac{1}{{\mu }_g}\right)}^{\frac{1}{4}}+\frac{S_o^{\prime }}{S_n^{\prime }}{\left(\frac{1}{{\mu }_o}\right)}^{\frac{1}{4}}$$

(11)

where

$$S_o^{\prime }=S_o-S_{\mathit{or}}$$

$$S_g^{\prime }=S_g-S_{\mathit{gi}}$$

$$S_n^{\prime }=S_o^{\prime }+S_g^{\prime }$$

Buckley-Leverett equation

The fundamental principle that underlies the description of flow in a porous medium is conservation of mass. For one-dimensional seepage of two phase in homogeneous isotropic horizontal porous medium, the following assumptions are made to obtain the analytical solution.

The two-phase flow is linear and horizontal.

The two fluids, crude oil and injected gas, are both incompressible.

The porous media is incompressible.

The effect of capillary and gravity is negligible.

The gas is injected at the inlet, x = 0.

The injection pressure and production pressure are constant.

Subject to the assumptions stated, the conservations of mass for the oil (o) and gas (g) phase become

$$\varphi \frac{\partial S_g}{\partial t}+\frac{\partial v_g}{\partial x}=0$$

(12)

$$\varphi \frac{\partial S_o}{\partial t}+\frac{\partial v_o}{\partial x}=0$$

(13)

Considering $S_o+S_g\text{=}1$ , the addition of equations (12) and (13) gives

$$\frac{\partial v_t}{\partial x}=\frac{\partial \left(v_g+v_o\right)}{\partial x}=0$$

(14)

where $v_t\equiv v_g+v_o$ , which indicates that the total volumetric flux in one-dimensional percolation keeps constant with spatial variation.

The classical fractional fractional-flow theory is based on the assumption of constant volumetric flux, $\frac{\partial v_t}{\partial t}=0$ . In this paper, it is assumed that the injection pressure and production pressure are constant over time. Figure 4 is a schematic diagram of the seepage model. The initial and boundary conditions of saturation and pressure are listed as follows.

Figure 4.

Linear displacement in a porous media.

Initial and boundary conditions of pressure

$$\begin{array}{c}p\left(0,t\right)=p_{\mathit{in}}\\ p\left(L,t\right)=p_{\mathit{out}}\\ p\left(x,0\right)=p_{\mathit{out}}\end{array}$$

(15)

Initial and boundary conditions of saturation

$$S_g\left(0,t\right)=1-S_{\mathit{or}},S_g\left(x,0\right)=S_{\mathit{gi}}$$

(16)

When co-seepaging in the porous media, the two phases are subject to Darcy’s Law. For the oil and gas phase, the volume fluxes are given separately as follows

$$v_o=-\frac{K{K}_{\mathit{ro}}}{{\mu }_{\mathit{oeff}}}\frac{\partial p}{\partial x}$$

(17)

$$v_g=-\frac{K{K}_{\mathit{rg}}}{{\mu }_{\mathit{geff}}}\frac{\partial p}{\partial x}$$

(18)

Adding the two equations above, the total volumetric flux can be expressed as

$$v_t=v_g+v_o=-\left({\lambda }_g+{\lambda }_o\right)\frac{\partial p}{\partial x}\text{=}-{\lambda }_t\frac{\partial p}{\partial x}$$

(19)

where ${\lambda }_t\left(S_g,\Delta p\right)\text{=}\frac{K{K}_{\mathit{ro}}\left(S_g,\Delta p\right)}{{\mu }_{\mathit{oeff}}\left(S_g\right)}$ , representing the total mobility of the fluids at given saturation and pressure.

The fractional-flow function for the gas phase is defined as follows

$$f_g=\frac{v_g}{v_g+v_o}$$

(20)

By substituting the Darcy equation into the definition formula, the fractional-flow function can be written as a function of gas saturation and pressure

$$f_g(S_g,p)=\frac{1}{1+\frac{K_{\mathit{ro}}(S_g,p)}{K_{r\text{g}}(S_g,p)}\cdot \frac{{\mu }_{\mathit{geff}}(S_g)}{{\mu }_{\mathit{oeff}}(S_g)}}$$

(21)

This equation can be used to determine the fractional-flow function of gas at any position of the reservoir after the gas saturation and pressure are given. In order to obtain analytical solution, the average pressure is used to determine the fractional-flow function. For the displacement process stated above where the pressure at both ends are given, the fractional-flow function is a function of saturation, as follows

$$f_g\left(S_g,p\right)={f_g\left(S_g\right)|}_{p=\overline{p}}$$

(22)

where $\overline{p}\text{=}\frac{p_{\mathit{in}}-p_{\mathit{out}}}{2}$ , the average pressure during the displacement.

By substituting the definition of fractional-flow function into equation (12), it can be rewritten as follows

$$\varphi \frac{\partial S_g}{\partial t}+v_t\frac{\partial f_g}{\partial t}=0$$

(23)

For the gas saturation is the function of time and distance, according to derivative method for compound function, equation (23) can be written as

$$\varphi \frac{\partial S_g}{\partial t}+v_t\frac{\partial f_g}{\partial S_g}\frac{\partial S_g}{\partial x}=0$$

(24)

The equation is referred to as the Buckley-Leverett equation, which constitutes the mathematical model in this paper coupled with the initial and boundary conditions of pressure and saturation.

Analytical solution

Frontal advance equation for gas flooding

The gas saturation on the section is constant, that is $d{S}_g=0$ . It is expressed as

$$d{S}_g=\frac{\partial S_g}{\partial t}dt+\frac{\partial S_g}{\partial x}dx=0$$

(25)

Combining the above equation with equation (24), the front advance velocity of a certain saturation is given as

$${\frac{\mathit{dx}}{\mathit{dt}}|}_{S_g}=\frac{v_t}{\varphi }{f^{\hspace{0.17em}\prime }}_g\left(S_g\right)$$

(26)

Integrating equation (19) between inlet and outlet ends, the following relation can be obtained

$$\Delta p\text{=}{p}_{\mathit{in}}-p_{\mathit{out}}=v_t{\displaystyle {\int }_0^L\frac{dx}{{\lambda }_t}}$$

(27)

The location of the front at a given time is x_f corresponding the frontal saturation S_gf. Dividing the integral into two parts, behind the gas front and ahead of the front, equation (27) can be rewritten as

$$\Delta p=v_t\left({\displaystyle {\int }_0^{x_f}\frac{dx}{{\lambda }_t}}+{\displaystyle {\int }_{x_f}^L\frac{dx}{{\lambda }_t}}\right)\text{=}{v}_t\left({\displaystyle {\int }_0^{x_f}\frac{dx}{{\lambda }_t}}+\frac{L-x_f}{{\lambda }_t\left(S_{\mathit{gi}},\Delta p\right)}\right)$$

(28)

Applying integration by parts twice and rearranging it give

$${\displaystyle {\int }_0^{x_f}\frac{\mathit{dx}}{{\lambda }_t}}=\frac{Q\left(t\right)}{\varphi }{\displaystyle {\int }_{1-S_{\mathit{or}}}^{S_{\mathit{gf}}}\frac{f^{\hspace{0.17em}″}\left(S_g\right)}{{\lambda }_t}\text{d}{S}_g}$$

(29)

where $Q\left(t\right)\text{=}{\displaystyle {\int }_0^t{v}_t\left(t\right)}\text{d}t$ .

The integral on the right side of the above equation is a constant. Let it be

$$\zeta \text{=}{\displaystyle {\int }_{1-S_{\mathit{or}}}^{S_{\mathit{gf}}}\frac{f^{″}\left(S_g\right)}{{\lambda }_t}\text{d}{S}_g}$$

(30)

Therefore, the total volumetric flux is shown as

$$v_t=\frac{\Delta p}{\frac{Q\left(t\right)}{\varphi }\zeta +\frac{L-x_f}{{\lambda }_t\left(S_{\mathit{gi}},\Delta p\right)}}$$

(31)

Integrating the front advance velocity, equation (26), the location of a given gas saturation at time t can be obtained

$$x\left(S_g,t\right)=\frac{Q\left(t\right)}{\varphi }{f^{\prime }}_g\left(S_g\right)$$

(32)

For the saturation of the shock front, S_gf, determined by the Welege’s (1952) method seen in Figure 5, the corresponding location of the shock front is

$$x_f=\frac{Q\left(t\right)}{\varphi }{f^{\prime }}_g\left(S_{\mathit{gf}}\right)$$

(33)

Figure 5.

Graphical determination of front saturation with Welge’s method.

Thus, the total volumetric flux can be rewritten as

$$v_t=\frac{\Delta p}{\left[\frac{\zeta }{{f^{\prime }}_g\left(S_{\mathit{gf}}\right)}-\frac{1}{{\lambda }_t\left(S_{\mathit{gi}},\Delta p\right)}\right]x_f+\frac{L}{{\lambda }_t\left(S_{\mathit{gi}},\Delta p\right)}}$$

(34)

Substituting the above equation into the front advance velocity, the following can be obtained

$$\frac{d{x}_f}{\mathit{dt}}=\frac{1}{\varphi }\frac{\Delta p{f^{\prime }}_g\left(S_{\mathit{gf}}\right)}{\left[\frac{\zeta }{{f^{\prime }}_g\left(S_{\mathit{gf}}\right)}-\frac{1}{{\lambda }_t\left(S_{\mathit{gi}},\Delta p\right)}\right]x_f+\frac{L}{{\lambda }_t\left(S_{\mathit{gi}},\Delta p\right)}}$$

(35)

Let

$$\begin{array}{l}2a=\frac{\zeta }{{f^{\prime }}_g\left(S_{\mathit{gf}}\right)}-\frac{1}{{\lambda }_t\left(S_{\mathit{gi}},\Delta p\right)}\\ b=\frac{L}{{\lambda }_t\left(S_{\mathit{gi}},\Delta p\right)}\\ c=\frac{\Delta p{f^{\prime }}_g\left(S_{\mathit{gf}}\right)}{\varphi }\end{array}$$

(36)

Therefore, the front advance velocity can be rewritten as

$$\frac{d{x}_f}{\mathit{dt}}=\frac{c}{a{x}_f+b}$$

(37)

Integrating the above equation from the inlet to the front, the location of the front is

$$x_f\text{=}\frac{-b+\sqrt{b^2+4ac\cdot t}}{2a}$$

(38)

Meanwhile, the total volumetric flux is determined as

$$v_t=\frac{\Delta p}{\sqrt{b^2+4ac\cdot t}}$$

(39)

When the front breaks through at the outlet, x_f in equation (38) equals the length of the porous media, which means x_f = L. Thus, the breakthrough time can be obtained

$$t_{\mathit{BT}}=\frac{a{L}^2+\mathit{bL}}{c}$$

(40)

The profiles of saturation and pressure

By integrating the front advance velocity of the corresponding saturation, the location of a certain saturation can be determined

$$x(S_g,t)=\frac{\Delta p{f}^{\hspace{0.17em}\hspace{0.17em}\prime }(S_g)}{\varphi }{\displaystyle {\int }_0^t\frac{1}{\sqrt{b^2+4ac\cdot t}}\text{d}t}$$

(41)

The pressure profile can be determined by integrating equation (19). The front is the dividing point of fluid distribution and the properties of fluid in the two parts are different leading to different pressure profile behind and ahead of the front.

$${\displaystyle {\int }_{p_{\mathit{in}}}^p\text{d}p=-{\displaystyle {\int }_0^x{\lambda }_t\text{d}x}}$$

(42)

Behind the gas front, the gas has swept the zone where the oil and gas coexist. Therefore, the pressure profile of this part can be determined by integrating from the inlet to the gas front. We can get

$$p\left(x,t\right)=p_{\mathit{in}}-v_t\frac{Q\left(t\right)}{\varphi }{\displaystyle {\int }_{1-S_{\mathit{or}}}^{S_g}\frac{f^{\hspace{0.17em}\hspace{0.17em}″}\left(S\right)}{{\lambda }_t\left(S\right)}\text{d}S}$$

(43)

Since the properties of the fluid vary along the path, the pressure profile is nonlinear which can be observed from equation (43).

Ahead of the gas front, for the single phase of the oil phase flowing in the porous media, the pressure profile which is linear, can be determined by integrating from the gas front to the outlet.

$$p\left(x,t\right)=p_{\mathit{out}}+\frac{v_t}{{\lambda }_t\left(S_{\mathit{gi}}\right)}\left(L-x\right)$$

(44)

Model verification

In this section, the verification is performed by comparing the analytical solution with the numerical solution. First, based on the gas flood model, the relative permeabilities of gas and oil under different condition are determined. Then, the saturation and pressure profiles of gas flooding are calculated from the analytical model stated above and numerical solution with finite difference method (FDM). In the numerical solution, the same models of viscosity and relative permeability are adopted to determine the varying properties of the two fluids.

The physical flow model is a one-dimensional linear porous medium, which is at first saturated uniformly with oil (S_oi = 0.95) and gas (S_gi = 0.05). A constant inlet and outlet pressure is imposed, starting from t = 0. The properties of porous media and fluids are shown in Table 1. The immiscible relative permeability curve is measured experimentally and the analytical expression is obtained by fitting the relative permeability index in the Corey-Brook model, which is convenient for the determination of first and the second derivative of fractional-flow fraction. The relative permeabilities under different miscible condition are determined through the model in the first part of the mathematical model and the fractional-flow is obtained by equation (21). Thus, the saturation and pressure profiles are calculated from equations (41), (43) and (44) separately.

Table 1.

The parameters used in the calculation.

Properties of the reservoir		Properties of the fluids
φ	0.2	μo (mPa·s)	1.81
K (mD)	20	μg (mPa·s)	0.035
L (m)	200	ρo (g/cm3)	0.7248
pi (MPa)	15.78	ρg (g/cm3)	0.4333
T (°C)	84	pb (MPa)	10.48
pin (MPa)	22	KRG	0.6181
pout (MPa)	10	KRO	0.9
		Sor	0.24
		Sgi	0.05

Relative permeability is an exponential function of saturation in the Corey-Brooks model. The relative permeability index can be obtained by fitting the double logarithmic curve of the relative permeability experimental data, the result of which is n_g = 2.1079 and n_o = 2.9852, shown in Figure 6(a) and (b). The relative permeability calculated by the Corey-Brooks model is in good agreement with the experimental data, as shown in Figure 7.

Figure 6.

(a) and (b) Log-log plot of gas saturation and relative permeability for gas and oil separately.

Figure 7.

The relative permeability in experiment and Corey model.

The oil recovery varying with pressure reaches its inflection at the pressure of 16.5 MPa, corresponding to the IFT of 1.97 cP, as shown in Figure 8. Therefore, the critical IFT is 1.97 cP. Effective miscible displacement could not be formed between oil and gas where the IFT is higher than critical IFT. The diffusivity of the gas was found to increase with temperature and pressure (Jamialahmadi et al., 2006). Shu et al. (2017) concluded that the pressure is the most important factor affecting the mass transfer during the gas injection. They pointed that the interface concentration of CO₂ in the oil phase increased with the pressure. The more intense the mass transfer between the fluids, the lower the IFT. When the IFT is lowered than 1.97 cP, the miscible condition of oil and gas will have a drastic change and the relative permeability will vary with pressure. Based on the relative permeability model of gas flooding established above, the relative permeability under different miscible condition is determined. Then, the fractional-flow curve is obtained from equation (21). Thus, the saturation and pressure profiles are calculated from equations (41), (43) and (44) separately.

Figure 8.

IFT between oil and gas and oil recovery in the experiment of slim tube versus different pressure.

In order to verify the analytical solution, the numerical solution is obtained with FDM using the same reservoir condition and physical properties. The geological model is turned into a discrete system by subdividing the x-axis into N grid blocks, as illustrated in Figure 9. The pressure is solved implicitly and the saturation is solved explicitly, and the calculated results are shown in Figures 10 and 11. To avoid the physically unreal oscillation caused by the central difference algorithm, the upwind style with first-order precision is used in the numerical solution, which unfortunately results in significant numerical dispersion and a smear front (Aziz, 1979; Wu, 2001), shown in Figure 12. It can be observed that the numerical solution deviates only slightly at the shock front. This deviation can be improved by refining the grid spacing. To reduce the effect of discretization on numerical simulation results, very fine, uniform mesh spacing is chosen. As can be seen from the figure, the impact of the discretization is very small when Δx = 0.5 m.

Figure 9.

Schematic of numerical simulation with FDM.

Figure 10.

Numerical simulation result for saturation profile.

Figure 11.

Numerical simulation result for pressure profile.

Figure 12.

Numerical solution of the Buckley-Leverett equation with different grid spacing.

The comparison between the analytical and numerical results is shown in Figures 13 and 14. The solid lines are the analytical results and the dashed lines are the numerical results. The results show that the analytical solution is in reasonable agreement with numerical solution. The shock front of saturation advances toward the outlet at a constant speed, which is consistent with the case of constant volumetric flux (Arabzai and Honma, 2013). The pressure profile is nonlinear behind the front and linear ahead of the front.

Figure 13.

Comparison between numerical and analytical results for saturation profile at different time.

Figure 14.

Comparison between numerical and analytical results for pressure profile at different time.

Results and discussion

For the gas flooding in porous media, the properties of porous media and fluids and the injection and production parameters dramatically impact the flow of the fluids, which determines the efficiency of gas displacement. Sensitivity analyses are conducted to investigate the effect of these parameters, including the permeability of the porous media, the viscosity and the injected pressure.

In the Buckley-Leverett equation, these parameters directly determine the fractional flow. As shown in equation (21), the fractional flow is a function of saturation, pressure, and viscosity. In order to study the effect on the fractional flow, the injection pressure is set to 14, 22 and 30 and the viscosity of the gas is set to 0.035, 0.105 and 0.175. Then, based on the relative permeability model of gas flooding established above, the relative permeability curve and fractional flow curve of different pressure and gas viscosity are determined, seen in Figures 15 and 16. With the increase of pressure, remarkable reduction in residual oil saturation and improvement of relative permeability are observed. Experimental data reveal the same variation rule, which is discussed thoroughly above. The essence of fractional-flow function is the mobility ratio between the displacement fluid and the displaced fluid. The results by Rui et al. indicated that the mobility ratio is a key element for determining the displacement efficiency (Rui et al., 2017). The typical fractional flow curve has a characteristic S-shape. With the increase of gas viscosity, the mobility of the displacing fluid increases and the mobility ratio is improved. Thus, the increase of fractional flow varying with gas saturation becomes slower and the front saturation decreases, similar to the result by Luo et al. (2017). For less viscous gas, once the gas breaks through, the well will produce a large amount of gas and the gas/oil ratio will rise sharply. The effect of pressure on fractional flow is complicated. Under the condition of low gas saturation, the fraction flow increases with the increase of pressure, contrary to the condition that the gas saturation is high.

Figure 15.

The relative permeability curve of oil and gas with different injection pressure.

Figure 16.

The fractional flow curve with different injection pressure and oil viscosity.

The mobility ratio of oil and gas is not affected by permeability, but the mobility does. Thus, the effect of viscosity, pressure, and permeability on the process of gas flooding from the aspects of saturation and pressure profiles and the breakthrough time will be discussed in detail.

Effect of permeability

Since the derivation of the classical Buckley-Leverett theory is based on the assumption of constant volumetric flux, the permeability has no effect on the front advance velocity. However, the situation is quite different in the case of constant pressure displacement at both ends. To study the effect of the permeability, the permeability is set from 20 to 120, while other parameters keep constant.

The total volume flux for different permeabilities is shown in Figure 17. With the increase of the permeability, the total flux increases correspondingly and this increase is more obvious over time. The increase of the permeability results in the improvement of the mobility, which makes the fluids flow more easily. Figure 18 is the saturation profile at 500 h for different permeabilities. As shown in the figure, the advance distance of the displacement front is longer as the permeability increases. The front advance velocity is proportional to the total flux. Therefore, the front advance velocity increases and the advance distance of the front becomes longer varying with permeability. In addition, it can be seen that the front saturation and residual oil saturation are not affected by permeability. Because the average pressure is constant, the miscible state does not change, neither does residual oil saturation. Since the mobility ratio is independent of permeability, the front saturation dominated by the mobility ratio is not affected by permeability. Consequently, the displacement efficiency is same, but the oil production rate is larger for higher permeability. The breakthrough time varying with permeability is shown in Figure 19. With the increase of permeability, the breakthrough time is shortened and the decrease of the breakthrough time becomes slower. Permeability has little effect on the breakthrough time when the permeability increases to a certain extent. In other words, the breakthrough time is more sensitive to the permeability for low or ultra-low permeability reservoirs. From equation (40), it can be observed that the breakthrough time is inversely proportional to permeability. The effect of the permeability has its two sides: in a reservoir with higher permeability, the fluids are easier to flow and the oil production is higher; the other side is that the breakthrough time will be shortened. Once the gas channeling is formed, the oil production rate will decline sharply and the gas–oil ratio will experience a rapid rise.

Figure 17.

Total volumetric flux varying with time for different permeability.

Figure 18.

Saturation profile at 500 h for different permeability.

Figure 19.

Breakthrough time for different permeability.

Effect of viscosity

The major challenges of the gas injection are the early breakthrough and the low sweep efficiency. Unfavorable mobility ratio resulted from the low viscosity of the injected gas in the reservoir condition is one of the core reasons. The gas viscosity is set from 0.035 to 0.21 to study its influence on the gas injection process.

The total volumetric flux varying with time for different viscosities is shown in Figure 20. With the increase of the gas viscosity, the total flux decreases. The less viscous gas has a greater mobility. At the early stage of gas injection, the flow of oil is dominated in the porous media. As a result, the variation of gas viscosity has little impact on the total flux. As the gas is continuously injected into the porous media, the gas begins to flow. With gas viscosity increases, the total volumetric flux decreases gradually. The longer the time, the more pronounced the effect. Figure 21 is the saturation profile at 3500 h for different gas viscosities. The effect of the gas viscosity on the saturation profile is reflected in two aspects. With the increase of the gas viscosity, on the one hand, the advance distance becomes shorter; on the other hand, the front saturation increases. As mentioned above, the decrease of mobility ratio leads to the change of the shape of fractional flow curve and the increase of the front saturation. The increase of the front saturation and the decrease of the advance distance lead to lower microscopic displacement efficiency of gas flooding. Nevertheless, the effect on the breakthrough time is quite different. The breakthrough time is longer with the increase of gas viscosity, as shown in Figure 22. The more viscous gas reduces the severity of the viscous fingering, prolongs the breakthrough time and improves the macroscopic displacement efficiency. Compared with the effect on the microscopic displacement efficiency, the improvement on the sweep efficiency is paid more attention. Hence, some researchers are engaging in increasing the viscosity of the injected gas by chemical additives.

Figure 20.

Total volumetric flux varying with time for different viscosity.

Figure 21.

Saturation profile at 3500 h for different viscosity.

Figure 22.

Breakthrough time for different viscosity.

Effect of pressure

During gas injection, the injection pressure is a very important parameter affecting the production performance. A higher pressure maintenance can change the miscibility of oil and gas, and improve the oil displacement efficiency substantially. The injection pressure is set from 14 to 30 to study its effect.

The total flux for different injection pressure is shown in Figure 23. With the increase of the pressure, the flux increases correspondingly. Other conditions being equal, higher injection pressure means greater pressure difference and higher injection rate, which result in the increase of the total flux. Figure 24 is the saturation profile at 2000 h for different pressure. The effect of the pressure on the saturation profile is profound. On the one hand, the residual oil saturation decreases with the increase of pressure, resulting in higher microscopic displacement efficiency. On the other, the front saturation decreases varying with the pressure. Moreover, the advance distance of the front is longer with higher injection pressure. When the pressure is increased to or near the MMP, miscible or near-miscible status is achieved, where the microscopic displacement efficiency is very high. The breakthrough time varying with pressure is shown in Figure 25. With the increase of the injection pressure, the breakthrough time is shortened and the decrease of the breakthrough time becomes slower. Considering the remarkable oil recovery, a higher pressure is preferable. The upper limit of the injection pressure is the maximum pressure that the equipment can provide or the fractured pressure of the formation.

Figure 23.

Total volumetric flux varying with time for different pressure.

Figure 24.

Saturation profile at 2000 h for different pressure.

Figure 25.

Breakthrough time for different pressure.

Summary and conclusions

This paper presents an analytical solution for gas flooding with constant-pressure boundary including the effect of miscibility by the relative permeability model and the viscosity model. From the analytical solution, saturation and pressure profiles, the total volumetric flux and the breakthrough time are determined. In order to verify the theory, the analytical solution is compared with numerical solution. Through the analytical solution, the effect of the properties of porous media and fluids on the performance of gas flooding is analyzed. Some insights can be obtained:

It is revealed that the analytic solution is in reasonable agreement with numerical solution, apart from the slightly deviation of numerical result at the shock front, which can be improved by refining the grid spacing. The analytical solution is used to obtain some insight into the effect of the properties of porous media and fluids on the performance of gas flooding.

The main effect of permeability is the front advance velocity since the total flux is increasing with the increase of permeability. In addition, the breakthrough time is shortened with the increase of permeability.

Gas viscosity has an effect on the front saturation and front advance velocity because of variations of gas mobility and mobility ratio with gas viscosity. Moreover, breakthrough time is longer with the increase of gas viscosity.

The effect of pressure is complicated. Permeability and gas viscosity have no effect on miscibility, but pressure does. With the pressure increase, there are remarkable reduction in residual oil saturation and improvement of relative permeability, resulting in higher gas saturation and oil displacement efficiency. Furthermore, higher injection pressure accelerates breakthrough of the injected gas.

Through the analyses of the results, the influence of these factors on gas flooding is clarified. This paper introduces the varying properties into the Buckley-Leverett equation and considers the constant-pressure boundaries, which expand the applied scope of the Buckley-Leverett theory. The analytical solution presented in this paper provides guidance on analyzing the distribution of saturation and pressure profiles and predicting the gas production and oil recovery efficiency of oil well. The model is based on the assumption that the fluids is incompressible However, this is not the case in the condition that the average pressure of the reservoir is low. Therefore, the model is suitable for the high pressure condition. Furthermore, it is assumed that the flow is one-dimensional and the vertical flow is neglected. Thus, the effect of gravity override on the breakthrough time is not taken into consideration. Future work will focus on (1) case study to predict the production performance of a real oil reservoir and (2) extension of the work to complex situations like capillary force, gravity, two-dimensional flow, etc.