Introduction {#Sec1}
============

In solution gas drive reservoirs, gas is released from solution when the reservoir pressure is below the bubble-point pressure. Gas initially exists in the form of tiny bubbles which is immobile due to the existence of capillary forces. Gas flows only after the bubbles become sufficiently large to coalesce into a continuous gas phase. A result of this process is that the production gas-oil ratio (GOR) rapidly increases after the critical gas saturation has been exceeded^[@CR1],[@CR2]^.

Field observations in some heavy-oil reservoirs in Venezuela, Canada, and China, however, are not consistent with the solution gas drive description because the production GOR remains relatively low^[@CR3],[@CR4]^. In addition, the oil production rates and oil recovery were reported to be significantly higher than the oil production rates and oil recovery that are expected from solution gas drive theory^[@CR5]--[@CR7]^. The most plausible explanation of the unusual production behavior is "foamy-oil" flow^[@CR8]--[@CR12]^, which is a different form of two-phase flow, in which the released gas flows in the form of dispersed bubbles instead of flowing as a continuous phase.

Since the conventional oil resources are dwindling, efficient exploitation of this type of heavy oil reservoir is imperative. Numerous efforts have been made to study foamy solution-gas drive processes^[@CR13]--[@CR22]^. Understanding the foamy oil mechanism and accurately modeling the foamy solution-gas drive processes in heavy oil reservoirs would improve the planning of the recovery of heavy oil.

In their experiments, Kumar *et al*.^[@CR13]^ concluded that high oil recovery can be achieved with high depletion rates. Zhang *et al*.^[@CR14]^ experimentally investigated foamy solution-gas drive processes with various temperatures and discovered that an intermediate temperature is beneficial for generating high oil recovery. Tang *et al*.^[@CR15]^ investigated the effect of overburden pressure by comparing foamy solution-gas drive in consolidated rock and foamy solution-gas drive in unconsolidated rock. They discovered that overburden pressure influences the features of flowing gas bubbles, critical gas saturation and oil recovery. Alshmakhy and Maini^[@CR16]^ reported an experimental study of the effect of gravity on oil recovery. They determined that gravity has a beneficial effect on recovery in some conditions because it increases the mobilizing force that is exerted on gas bubbles and encourages foamy oil flow.

In addition to these factors, oil, gas and rock properties also influence foamy solution-gas drive. High permeabilities and long sandpacks cause high supersaturation, high critical gas saturation, and high oil recovery^[@CR17]^. An increase in initial water saturation causes a decreases in gas-bubble density, which produces low oil mobility and high gas mobility^[@CR18]^. The oil viscosity and asphaltenes help in maintaining foamy oil flow by delaying bubble coalescence. However, they have a minor influences on oil recoveries with high depletion rates^[@CR19],[@CR20]^. In addition, foamy solution-gas drive processes depend on the saturated gas types. The oil recoveries are achieved with CH~4~ saturated oil are higher than the oil recoveries achieved with CO~2~ saturated oil^[@CR21],[@CR22]^.

Numerical simulation offers another effective pathway for the study of foamy solution-gas drive processes. In the early stages of research, conventional solution-gas drive models with empirical adjustments were employed to simulate foamy oil phenomena^[@CR23]^. These models can be used to obtain a history match of reservoir performance. However, unrealistic parameters are often obtained due to unreasonable adjustments that account for the contributions of foamy-oil flow to oil recovery. Therefore, the predictions from these models are questionable.

Then, three types of equilibrium models were proposed based on the assumption of complete equilibrium between the gas phase and the oil phases: pseudo bubble point models^[@CR24]--[@CR26]^, modified fractional-flow models^[@CR27]--[@CR29]^ and reduced viscosity models^[@CR30]^. Foamy oil is a non-thermodynamically stable system^[@CR31]^. If time is sufficient and an environment is helpful, the dispersed gas will separate from the oil phase^[@CR32],[@CR33]^. Thus, the equilibrium models are inherently incapable of reflecting this thermodynamically unstable characteristic of foamy oil and cannot be employed to simulate foamy oil flow.

Some kinetic models have been employed in foamy oil simulation studies to capture the time-dependent characteristic of a foamy oil system. For example, Coombe and Maini^[@CR34]^ proposed a kinetic model that accounts for the changes of foamy oil morphology. Sheng *et al*.^[@CR29]^ described a foamy oil model that includes two kinetic equations that control the transfers among three gas components: solution gas, evolved gas and free gas. Shahvali and Pooladi-Darvish^[@CR35]^ proposed a one-dimensional dynamic model that considers the effects of viscous forces on gas mobility. In recent studies, Liu *et al*.^[@CR36]^ proposed a new mathematical model for foamy-oil flow that consider foamy-oil supersaturation. Oskouei *et al*.^[@CR37]^ developed a kinetic model that applies a non-equilibrium mass transfer approach to model gas evolution in heavy oil.

These models are capable of simulating the one-dimensional kinetic behavior of foamy oil flow by the history matching production behavior of pressure depletion tests. However, they are unable to predict foamy solution-gas drive processes on an oilfield scale. In addition, the models are unable to simultaneously simulate multiple characteristics of foamy solution-gas drive. The effect of microscopic behavior, such as gas bubble nucleation and growth, on foamy solution-gas drive processes is not well understood. Therefore, a new three-dimensional (3D) foamy oil model that captures many important characteristics of foamy solution-gas drive (non-equilibrium behavior, gas evolution kinetics, and the effect of viscous forces on gas mobility) based on pressure depletion tests was developed in this study. Then, the effects of some important parameters on foamy solution-gas drive were systematically investigated, and a model application was conducted in a typical foamy oil reservoir.

Experimental and Modeling {#Sec2}
=========================

Experimental materials {#Sec3}
----------------------

A crude oil from the Orinoco heavy oil belt, which is located east of Venezuela was employed in the pressure depletion tests. The live oil was obtained by mixing crude oil and gas (11.1% carbon dioxide +88.9% methane). The methane gas and carbon dioxide both have purities of +99%. The physical properties of Venezuelan oil are listed in Table [1](#Tab1){ref-type="table"}. As previously stated, the heavy oil in this area has foamy oil tendencies.Table 1Summary of the Venezuelan oil used in the experiments.Fluid propertyGOR (m^3^/m^3^)15Oil density at 20 °C (g/cm^3^)1.013Oil viscosity at 50 °C (Pa•s)24.7Oil FVF at 54.2 °C and 8.65 MPa (m^3^/m^3^)1.173Density at 54.2 °C and 8.65 MPa (g/cm^3^)0.957Acid number (mg/g)4.95Bubble point pressure (MPa)4.95 MPaPseudo bubble point pressures (MPa)3.44 when the resting time is 24 hours2.74 when the resting time is 12 hours1.89 when the resting time is 1 hourMonophasic fluid (mol%)Nitrogen0.13Carbon Dioxide2.8Methane22.43Propane0.08I-Butane0.04N- Butane0.04I-Pentane0.12N- Pentane0.49C~7+~73.87MW418.76SARA analysis (wt.%)Asphaltenes7.78Saturates17.19Aromatics39.17Resins35.86

Experimental setup {#Sec4}
------------------

A schematic of the pressure depletion experimental setup is illustrated in Fig. [1](#Fig1){ref-type="fig"}. The length and diameter of the sandpack was 50 cm and 3.8 cm, respectively; the sandpack was horizontally located in a thermostatic convection oven. The accuracy of the temperature control was ±0.1 °C. Three pressure taps monitored the pressure distribution of the sandpack. A computer was used to record the pressures from the transducers during each test.Figure 1Schematic of the pressure depletion experimental setup.

Brine (0.5% NaCl) was contained in a transfer vessel, and live oil was obtained by saturating the crude oil with the gas in a mixer connected to the inlet end of the sandpack. The brine and live oil were injected into the sandpack by a displacement pump to saturate it prior to each test.

A displacement pump that connects a backpressure regulator (BPR) through a buffer piston vessel with water was used to control the pressure at the outlet end of the sandpack. The fluid was produced and flowed into an oil-gas separator that was located on an electronic balance to measure the produced oil. The separated gas was monitored by a gas meter that was connected to the oil-gas separator.

Experimental procedures {#Sec5}
-----------------------

According to the solution GOR (15 m^3^/m^3^) of the heavy oil, live oil at reservoir conditions (54.2 °C and 8.65 MPa) was prepared with crude oil and gas. The gas Z factor for the preparation of live oil was 0.915; it was measured by several separator tests.

After packing, the sandpack was evacuated by a vacuum pump. Then, the porosity, permeability and oil saturation of the sandpack were obtained by filling it with brine and live oil. The pressure depletion test was started after the sandpack was saturated with the live oil. The pressure was continuously decreased by discharging water from the buffer piston vessel that was connected to the displacement pump. The water discharged at a constant volume flow rate which provided a fixed depletion rate. Due to the decrease in the sandpack pressure, the fluids produced from the sandpack flowed into the oil-gas separator, and the weight of oil and the volume of gas were periodically determined by the electronic balance and the gas meter. The test was finished when the sandpack pressure declined to a relatively low value (near atmospheric pressure).

During this experimental process, the pressures at different locations were monitored with the pressure transducers. The oven temperature was maintained at a reservoir temperature of 54.2 °C. To investigate the effect of depletion rate on foamy solution-gas drive processes, three series of tests with various depletion rates were performed in our study. The properties of the sandpacks and experimental conditions are listed in Table [2](#Tab2){ref-type="table"}.Table 2Properties of the sandpacks and the experimental results at different depletion rates.Test No.Porosity (%)Oil saturation (%)Permeability (D)Depletion rate (MPa/h)Cumulative oil produced (g)Maximum oil production rate (g)Maximum cumulative GOR (m^3^/m^3^)Oil recovery (%)142967.160.437.45.0790.517.2242967.110.844.313.0770.219.1341987.061.66317.8745.624.6

Modeling {#Sec6}
========

Model assumptions {#Sec7}
-----------------

The major assumptions in deriving the 3D foamy oil model are as follows:The foamy oil system is considered to have two phases (oil and gas) and three components (solution gas, free gas, and oil). After gas bubble nucleation, liberation of solution gas in the oil phase is transferred to free gas in the gas phase. This process is described in the continuity equations by the sink/source term *q*~sf~, which is calculated using an equation that considers the kinetics of gas bubble nucleation and growth.The density of the solution gas is equivalent to the density of the free gas, and the velocity in the solution gas in oil phase is equivalent to the velocity in the velocity in the oil phase.The oil and gas relative permeabilities obey Corey-type functions, which incorporate the effects of the velocity and the oil viscosity on gas mobility.The gravity effect is considered in the model, whereas the capillary effect is neglected.The reservoir rock is slightly compressible.The reservoir temperature is isothermal.

Mathematical model {#Sec8}
------------------

According to the mass balance for each component, the continuity equations are derived as follows:$$\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$${\rm{Oil}}\,{\rm{component}}:\nabla \cdot [\frac{k{k}_{ro}}{{\mu }_{o}{B}_{o}}\nabla (p-{\rho }_{o}\,gH)]+{q}_{vo}=\frac{\partial }{\partial t}(\frac{\varphi {S}_{o}}{{B}_{o}})$$\end{document}$$$$\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$${\rm{Solution}}\,{\rm{gas}}\,{\rm{component}}:\nabla \cdot [{R}_{s}\frac{k{k}_{ro}}{{\mu }_{o}{B}_{o}}\nabla (p-{\rho }_{o}\,gH)]-{q}_{sf}+{R}_{s}{q}_{vo}=\frac{\partial }{\partial t}(\frac{\varphi {R}_{s}{S}_{o}}{{B}_{o}})$$\end{document}$$$$\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$${\rm{Free}}\,{\rm{gas}}\,{\rm{component}}:\nabla \cdot [\frac{k{k}_{rg}}{{\mu }_{g}{B}_{g}}\nabla (p-{\rho }_{g}gH)]+{q}_{sf}+{q}_{vg}=\frac{\partial }{\partial t}(\frac{\varphi {S}_{g}}{{B}_{g}})$$\end{document}$$where *k* is the absolute permeability, D; *k*~rg~ and *k*~ro~ are the gas relative permeability and oil relative permeability, respectively; *p* is the reservoir pressure, MPa; *ρ*~g~ and *ρ*~o~ are the gas density and oil density, respectively, kg/m^3^; g is the acceleration of gravity, m/s^2^; *H* is the standard height, m; *μ*~g~ and *μ*~o~ are the gas viscosity and oil viscosity, mPa·s; *S*~g~ and *S*~o~ are the gas saturation and oil saturation, respectively; *ϕ* is the porosity; *B*~g~ and *B*~o~ are the gas formation volume factor and the oil formation volume factor, respectively, m^3^/m^3^; *R*~s~ is the solution GOR, m^3^/m^3^; *q*~vo~ and *q*~vg~ are the oil volumetric rate and the gas volumetric rate, respectively, per unit bulk volume, m^3^/(m^3^·s); *t* is time, s; and *q*~sf~ is the rate of gas liberation from solution, cm^3^/s.

The oil and gas saturations satisfy the following equation:$$\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$${S}_{g}+{S}_{o}=1$$\end{document}$$The necessary initial and boundary conditions are specified as follows:$$\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$${\rm{The}}\,{\rm{initial}}\,{\rm{conditions}}:\,p(t){|}_{t=0}={p}^{0}{S}_{o}(t){|}_{t=0}={S}_{o}^{0}{S}_{g}(t){|}_{t=0}={S}_{g}^{0}$$\end{document}$$$$\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$${\rm{The}}\,{\rm{outer}}\,{\rm{boundary}}\,{\rm{condition}}:\,\frac{\partial p}{\partial n}{|}_{{\rm{\Gamma }}}=0$$\end{document}$$$$\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$${\rm{The}}\,{\rm{inner}}\,{\rm{boundary}}\,{\rm{condition}}:\begin{array}{rcl}p({x}_{w},{y}_{w},{z}_{w},t) & = & {p}_{wf}({x}_{w},{y}_{w},{z}_{w})\\ {\rm{or}}\,\,\,Q({x}_{w},{y}_{w},{z}_{w},t) & = & Q({x}_{w},{y}_{w},{z}_{w})\end{array}$$\end{document}$$where *P*^0^ is the initial reservoir pressure, MPa; *S*~o~^0^ and *S*~g~^0^ are the initial oil saturation and gas saturation, respectively; Γ represents the outer boundary; *n* represents the normal direction of the boundary; *p*~wf~ is the bottom-hole pressure of the production wells, MPa; *Q* is the production rate of the production wells, m^3^/d; and *x*~w~, *y*~w~, and *z*~w~ are the 3D coordinates of the production wells in the *x* direction, *y* direction and *z* direction, respectively. Eq. ([1](#Equ1){ref-type=""}) -- Eq. ([7](#Equ7){ref-type=""}) complete the mathematical description of foamy oil model.

Calculation of *q*~sf~ {#Sec9}
----------------------

Gas bubble nucleation is assumed to occur when the supersaturation Δ*p* exceeds the capillary forces of a crevice^[@CR38]^. Therefore, a crevice of radius *r*~p~ is activated when the inequality defined in the following equation is satisfied:$$\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$${\rm{\Delta }}p={p}_{g}-{p}_{o}={p}_{e}-p\ge \frac{2\gamma }{{r}_{p}}$$\end{document}$$where *γ* is the interfacial tension, dyne/cm; *p* is the actual pressure within an oil-filled pore, MPa; *p*~e~ is the equilibrium pressure, MPa; and *p*~g~ and *p*~o~ are pressures in the gas and the oil, respectively, MPa. In this model, *r*~p~ is assumed to be randomly distributed in the porous media, which is quantitatively described using the Rayleigh distribution function *F*(*r*~p~):$$\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$F({r}_{p})=\frac{{r}_{p}}{{\sigma }^{2}}\exp (-\frac{{{r}_{p}}^{2}}{2{\sigma }^{2}})$$\end{document}$$In this function, the Rayleigh distribution deviation *σ* is an important parameter to determine the distribution of *r*~p~. Assume that the value of supersaturation at a specific time is represented by Δ*p*~1~, at this level of supersaturation, the pores that contain activated sites are pores with a radius larger than *r*~pmin~ (*r*~pmin~ is equal to 2γ/Δ*p*~1~ based on Eq. ([8](#Equ8){ref-type=""})). Gas bubble nucleation occurs in pores with radii between *r*~pmin~ and the maximum radius *r*~pmax~. Therefore, the number of bubbles per unit bulk volume *N* can be obtained by the following equation:$$\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$N=\frac{\varphi {\int }_{{r}_{pmin}}^{{r}_{pmax}}F({r}_{p})d{r}_{p}}{{V}_{s}}$$\end{document}$$where *V*~s~ is the average volume of one pore (assuming the existence of one nucleus bubble in each pore). The growth of the activated gas bubbles is controlled by diffusion of the solution gas in the supersaturated oil. Therefore, the rate of gas bubble growth is^[@CR39]^:$$\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\frac{dr}{dt}=\frac{DRT{\rm{\Delta }}C}{pr}$$\end{document}$$where *r* is the gas bubble radius, m; *D* is the diffusion coefficient, m^2^/s; R is the universal gas constant, J/(mol·K); Δ*C* is the supersaturation expressed as a difference in concentration, mol/m^3^; and *T* is the temperature, K.

If the radius of a gas bubble is increased from *t*^n^ to *t*^n+1^, *r* at time *t*^n+1^ is:$$\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$r({t}^{n+1})=r({t}^{n})+\frac{dr}{dt}{\rm{\Delta }}t$$\end{document}$$where Δ*t* is the increment of time from *t*^n^ to *t*^n+1^, s.

If the gas bubbles are assumed to be spherical, the total volume of the gas bubbles per unit bulk volume is:$$\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$V=N\frac{4}{3}\pi {r}^{2}$$\end{document}$$The volume increase of the gas bubbles from *t*^n^ to *t*^n+1^ is:$$\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\frac{dV}{dt}=\frac{N4\pi {r}^{2}dr}{dt}=\frac{N4\pi DRT{\rm{\Delta }}C(t)r({t}^{n})}{p}$$\end{document}$$For a pre-existing bubbles nucleation model, where new bubbles at different population densities continuously activate, Eq. ([14](#Equ14){ref-type=""}) should be changed to the following form:$$\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\frac{dV}{dt}=\frac{N4\pi {r}^{2}dr}{dt}=\frac{4\pi DRT{\rm{\Delta }}C(t)}{p}\sum _{j}{N}_{j}{r}_{j}$$\end{document}$$where the subscript *j* denotes the *j*^*th*^ time step in which new nuclei at a density of *N*~*j*~ are activated; and *r*~*j*~ is the radius of a gas bubble at the *j*^*th*^ time step.

The volumetric rate of gas liberation from solution per unit bulk volume at standard conditions can be calculated by:$$\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$${q}_{sf}=\frac{dV}{dt}\frac{1}{{B}_{g}}=\frac{4\pi DRT{\rm{\Delta }}C(t)}{p{B}_{g}}\sum _{j}{N}_{j}{r}_{j}$$\end{document}$$The gas concentration can be related to the solution GOR using the following equations:$$\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$${\rm{\Delta }}C(t)=C-{C}_{eq}=\frac{{R}_{s}{\rho }_{gsc}}{{B}_{o}}-\frac{{R}_{seq}{\rho }_{gsc}}{{B}_{o}}$$\end{document}$$where *C* is the average gas concentration in the liquid phase, mol/m^3^; *C*~*eq*~ is the equilibrium gas concentration, mol/m^3^; *R*~seq~ is the equilibrium solution GOR, m^3^/m^3^; and *ρ*~gsc~ is the gas molar density at standard conditions.

Eq. ([18](#Equ18){ref-type=""}) is the final equation applied in the model to account for gas evolution from solution:$$\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$${q}_{sf}=\frac{4\pi DRT(R-{R}_{eq}){\rho }_{gsc}}{p{B}_{o}{B}_{g}}\sum _{j}{N}_{j}{r}_{j}$$\end{document}$$

Oil-gas relative permeabilities {#Sec10}
-------------------------------

During foamy solution-gas drive processes, gas is dispersed in the oil phase to generate foamy oil due to the high viscous forces of heavy oil, which cause a decrease in the gas mobility^[@CR40]--[@CR42]^. Thus, oil and gas relative permeabilities depend not only on gas saturation but also on viscous forces. In this study, the effect of viscous forces on gas mobility are represented by oil and gas relative permeabilities, which are given by the following equations^[@CR38]^:$$\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$${k}_{ro}={k}_{ro}^{0}{(\frac{1-{S}_{g}}{1-{S}_{gc}})}^{{n}_{o}}\,{S}_{gc}={S}_{gci}+\beta v{{\mu }_{o}}^{1/2}$$\end{document}$$$$\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$${k}_{rg}={k}_{rg}^{0}{(\frac{{S}_{g}-{S}_{gc}}{1-{S}_{gc}})}^{{n}_{g}}\,{k}_{rg}^{0}={k}_{rgi}\exp (-\alpha v{{\mu }_{o}}^{\frac{1}{2}})$$\end{document}$$where *S*~*gc*~ is the critical gas saturation; *n*~*o*~ is the Corey exponent for oil; *n*~*g*~ is the Corey exponent for gas; *k*~*rg*~^0^ and *k*~*ro*~^0^ are the end-point gas and oil relative permeability, respectively; *k*~*rgi*~ is the maximum value of gas relative permeability; *S*~*gci*~ is the minimum value of critical gas saturation; *α* and *β* are proportionality constants, s/m; and *v* is the oil phase velocity, m/s.

Other properties {#Sec11}
----------------

Oil and gas production rates:$$\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$${q}_{vo}=\frac{Q}{{V}_{B}{B}_{o}(1+(\frac{{k}_{rg}}{{k}_{ro}})(\frac{{\mu }_{o}}{{\mu }_{g}}))}\,{q}_{vg}=\frac{Q(\frac{{k}_{rg}}{{k}_{ro}})(\frac{{\mu }_{o}}{{\mu }_{g}})}{{V}_{B}{B}_{g}(1+(\frac{{k}_{rg}}{{k}_{ro}})(\frac{{\mu }_{o}}{{\mu }_{g}}))}$$\end{document}$$$$\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$${\rm{Gas}}\,{\rm{FVF}}:\,{B}_{g}=\frac{T{p}_{sc}Z}{p{T}_{sc}}$$\end{document}$$$$\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$${\rm{Oil}}\,{\rm{FVF}}:{B}_{o}={B}_{ob}[1+{C}_{o}(p-{p}_{b})]\,\,p\ge {p}_{b}$$\end{document}$$$$\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$${B}_{o}=1+\frac{{B}_{ob}-1}{{p}_{b}}p\,p < {p}_{b}$$\end{document}$$where *V*~B~ is the volume of a grid block, cm^3^; *Z* is the compressibility factor; *p*~sc~ is the standard pressure, MPa; *p*~*b*~ is the bubble point pressure, MPa; *B*~*ob*~ is the oil formation volume factors at bubble point pressure; *C*~*o*~ is the oil compressibility coefficient, 1/MPa; and *T*~sc~ is the standard temperature, °C.

Solution of the mathematical model {#Sec12}
----------------------------------

The model Eqs ([1](#Equ1){ref-type=""}--[7](#Equ7){ref-type=""}) were solved using the IMPES method. First, the continuity equations for the three components multiplied by appropriate coefficients were merged into a pressure equation. In a 3D model with *Nx*, *Ny*, and *Nz* blocks, 3*N* pressure equations were generated. The system was solved by a conjugate gradient method for the pressure of each block. Once the pressures were obtained, *S*~g~ and *R*~s~ were explicitly obtained, and *S*~o~ was calculated by Eq. ([4](#Equ4){ref-type=""}). To efficiently and accurately solve the mathematical model, we developed a 3D two-phase flow numerical simulator.

Results and Discussion {#Sec13}
======================

Experimental results {#Sec14}
--------------------

Figure [2](#Fig2){ref-type="fig"} and Table [2](#Tab2){ref-type="table"} provide the experimental results for different depletion rates. The pressure in Fig. [2](#Fig2){ref-type="fig"} is the average pressure of the sandpack. Figure [2](#Fig2){ref-type="fig"} shows three typical stages for the foamy solution-gas drive in all tests. In stage 1, when the pressure exceed *P*~b~ (4.95 MPa), the cumulative oil produced, the cumulative GOR and the oil production rate are very low. Because the solution gas is not released from the oil and the fluid in the sandpack flowed as a single oil phase in this stage. The oil is primarily produced by the expansion of the low compressibility oil phase. However, once the pressure is below *P*~b~ (Stage 2), the cumulative GOR remains nearly constant for a period of time, and the cumulative oil produced and oil production rate rapidly increase. The unusual performances are caused by the effect of foamy oil flow. When the sandpack pressure was below *P*~b~, the gas was released from solution and dispersed in the oil as tiny bubbles rather than coalescing to form a continuous gas phase that causes fast gas breakthrough, energy loss and decreased oil production. The oil with dispersed bubbles (foamy oil phenomena) causes oil volume swelling to provide energy and prevent a decrease in pressure and oil production. Since the gas does not become continuous in this stage, its mobility remains low, which causes stable cumulative GOR.Figure 2Cumulative oil produced, cumulative GOR and oil production rate vs. pressure for the three pressure depletion tests: (**a**) 0.4 MPa/h; (**b**) 0.8 MPa/h; (**c**) 1.6 MPa/h.

As pressure depletion continues, at a certain pressure, the critical gas saturation (*S*~gc~) is attained, and the cumulative GOR starts to linearly increase with pressure, which indicates that the gas bubbles in oil begin to coalesce and form a continuous gas phase and foamy oil phenomena begin to disappear. Note that the production of foamy oil continues although free gas starts to flow at this pressure. This finding is consistent with the observations reported by Alshmakhy and Maini^[@CR16]^.

In stage 3, the cumulative oil produced and oil production rate become low, and the cumulative GOR continually increases due to the complete disappearance of foamy oil phenomena in this stage. The gas bubbles trapped in the oil phase completely coalesce and generate free gas. Because the gas mobility is higher than the oil mobility, which causes gas breakthrough, a significant decrease in the cumulative oil produced and oil production rate, and an increase in the cumulative GOR. The three typical stages can be used to determine whether a heavy oil exhibits significant foamy oil tendencies during solution-gas drive processes.

As shown in Table [2](#Tab2){ref-type="table"}, a higher depletion rate produces a higher cumulative oil produced, higher maximum oil production rate, higher oil recovery and lower maximum cumulative GOR. High depletion rates benefit foamy solution-gas drive processes. The mechanism responsible for this pronounced effect of depletion rate is ambiguous and can be related to the depletion rate dependence of *S*~gc~. The highest *S*~gc~ is achieved at the highest depletion rate (Fig. [2](#Fig2){ref-type="fig"}), which decreases the free gas mobility. In addition, high pressure depletion rates generate a large amount of supersaturation. Additional nucleation sites are activated and smaller gas bubbles are formed.

Modeling {#Sec15}
--------

### Verification of the new developed model {#Sec16}

The proposed foamy oil model was validated by comparing the experimental data of the pressure depletion tests in this study with the experimental data of the pressure depletion tests in Kumar's study^[@CR14]^. The rock and fluid properties corresponding to the experimental data are listed in Tables [2](#Tab2){ref-type="table"} and [3](#Tab3){ref-type="table"}. During the history-matching processes, the adjusted parameters in the model were parameters related to the oil and gas relative permeabilities (*S*~*gci*~, *k*~*rgi*~, *α*, and *β*), gas bubble nucleation and growth (*D*, σ, and *r*~pmax~). The results from the developed foamy model and the experiments are compared in Fig. [3](#Fig3){ref-type="fig"}. As depicted by Fig. [3](#Fig3){ref-type="fig"}, fairly reasonable matches are obtained.Table 3Values in the models before and after matching for this study and Kumar's study.ParametersThis studyKumar's studyBefore matchingAfter matchingBefore matchingAfter matchingMaximum *k*~*ro*~^0^1111Maximum *k*~*rgi*~0.10.060.080.05*S* ~*gci*~0.020.0250.0250.02*n* ~*o*~2222*n* ~*g*~2222*r*~p*max*~ (m)1.0 × 10^−7^1.2 × 10^−7^0.8 × 10^−7^1.0 × 10^−7^*σ* (m)2.0 × 10^−8^1.7 × 10^−8^1.0 × 10^−8^1.5 × 10^−8^*α* (s/m)1.0 × 10^6^1.2 × 10^6^1.0 × 10^5^8.4 × 10^6^*β* (s/m)1.0 × 10^4^1.5 × 10^4^1.0 × 10^4^1.2 × 10^4^*D* (m^2^/s)1 × 10^−9^3 × 10^−9^1 × 10^−9^3 × 10^−9^Figure 3Comparison of the experimental data and simulated results: (**a**) this study; (**b**) Kumar's study.

In addition, when pressure exceeds *P*~b~, no gas bubbles were released from the oil phase to form foamy oil. At this stage, only single oil flow exists in the reservoir. Thus, the accuracy of the foamy oil model for the flow process above *P*~b~ can be also validated against the known analytical solution of a one-dimensional single-phase flow model^[@CR6]^ and Eclipse's black oil simulator E100.

The simulation results of pressures above *P*~b~ at different times are compared with the simulation results obtained by the analytical solution, which is presented in Fig. [4](#Fig4){ref-type="fig"}. The analytical solution matches the calculated results of the model. Comparisons of the solutions of the 3D foamy oil model and solutions obtained from E100 indicate that the average relative error in the pressure values in each grid between the 3D developed model and the E100 model is 0.103% after the simulation calculations (t = 30 days). In addition, Fig. [5](#Fig5){ref-type="fig"} implies agreement between the calculated pressure values in a 3D reservoir obtained by the developed model and E100. The simulated results of the developed model are consistent with the results of the analytical solution and E100.Figure 4Comparison of the pressures above Pb vs. distance at different times.Figure 5Comparison of the pressure distributions and specific pressure values.

### Study of affecting factors {#Sec17}

The model after the experimental history-match was employed as the base case; its parameters are listed in Tables [2](#Tab2){ref-type="table"} and [3](#Tab3){ref-type="table"}. To align the range of the factors analyzed with the actual conditions of the Venezuelan reservoir, the model parameters were adjusted with a constant factor of 25% based on the base case (Table [4](#Tab4){ref-type="table"}).Effect of critical gas saturation (*S*~gc~)As shown in Fig. [6](#Fig6){ref-type="fig"}, the higher is the oil recovery, the lower is the cumulative gas produced and the longer is the predicted time for gas produced for a higher *S*~gc~. However, the pressure is not significantly affected by *S*~gc~. Note that the definition of *S*~gc~ is the gas saturation when the gas starts to flow as a continuous phase. Thus, a high *S*~gc~ causes delay in the time for gas produced, which produces high oil recovery and low cumulative gas produced.Figure 6Effect of Sgc on foamy solution-gas drive processes: (**a**) oil recovery; (**b**) cumulative gas produced; (**c**)average pressure.In addition, the calculated *N* (the number of bubbles per unit bulk volume) and *r* (the radius of gas bubbles) in every time step from each run have the same values (the maximum *N* and minimum *r* in each run are 2.04 × 10^10^ and 8.63 × 10^−8^ m, respectively). Thus, *S*~gc~ never affects the gas bubble nucleation and growth behavior (the gas bubble nucleation and growth behavior is primarily related to the pressure), which explains why *S*~gc~ does not significantly influence the pressure profiles (Fig. [6c](#Fig6){ref-type="fig"}).Effect of oil viscosity (*μ*~o~)As revealed by Fig. [7](#Fig7){ref-type="fig"}, increasing *μ*~o~ causes delay of the time in gas produced, a decrease in cumulative gas produced, and an increase in the oil recovery. The results are consistent with the results of an experimental study by Pooladi-darvish and Firoozabadi^[@CR19]^. As shown in Fig. [7c](#Fig7){ref-type="fig"}, increasing the oil viscosity affects the early time stage of the pressure profiles. The pressure rebound is small when the oil viscosity is high because high oil viscosity inhibits the gas diffusion process, which decreases the values of *q*~sf~ (the volumetric rate of gas liberation from solution). For example, when the oil viscosities are 2086.5 mPa·s and 6259.5 mPa·s, the calculated maximum *q*~sf~ are 4.81 × 10^−4^ m^3^/(m^3^⋅s) and 2.68 × 10^−4^ m^3^/(m^3^⋅s), respectively. Thus, increasing the oil viscosity affects the early time stage of foamy solution-gas drive processes. As shown in Fig. [7b](#Fig7){ref-type="fig"}, high oil viscosity decreases the gas mobility. Increasing the oil viscosity from 2086.5 mPa·s to 6259.5 mPa·s increases the critical gas saturation *S*~gc~ from 4.09% to 5.26%, respectively. Therefore, increasing oil viscosity can also affect the late-time stage of foamy solution-gas drive processes due to the decrease in gas mobility.Figure 7Effect of μo on foamy solution-gas drive processes: (**a**) oil recovery; (**b**) cumulative gas produced; (**c**)average pressure.These studies explain how oil viscosity influences foamy solution-gas drive processes. However, their findings do not explain whether a reduction in *q*~sf~ or a reduction in gas mobility enables high viscosity oil to enhance the performance of foamy solution-gas drive processes. Because the calculated value of *q*~*sf*~ continuously changes during the simulation processes, the effect of *q*~*sf*~ on foamy solution-gas drive processes can not be directly investigated by comparing the results of different runs with fixed *q*~sf~ values. According to Eq. ([18](#Equ18){ref-type=""}), the diffusion coefficient *D* is positively proportional to *q*~sf~, and it does not affect gas mobility. Therefore, *D* can be applied to indirectly investigate the effect of *q*~sf~ on foamy solution-gas drive processes.The results indicate that an increase in *D* decreases the non-equilibrium effects and the supersaturation level in the early time of the average pressure profiles because a high *D* facilitates the growth of gas bubbles, and gas liberates from solution at a high rate, which is equivalent to increasing the value of *q*~sf~. However, the cumulative gas produced and oil recovery are insensitive to changes in *D* for the range of the study. Therefore, we infer that oil viscosity affects foamy solution-gas drive by decreasing gas mobility, rather than inhibiting the growth of gas bubbles (decreasing the *q*~sf~).Effect of depletion rateThe simulated results indicate that efficient oil recovery is obtained at higher depletion rates, which indicates that high depletion rates are beneficial to foamy solution-gas drive processes. These predictions are consistent with the experimental observations for the following reasons:First, the higher is the pressure depletion rate is, the larger is the formed supersaturation. Therefore, a larger number of nucleation sites are activated and smaller gas bubbles are formed. For example, when the pressure depletion rates are 1.5 m^3^/h and 4.5 m^3^/h, the calculated maximum *N* (the number of bubbles per unit bulk volume) are 8.13 × 10^9^ and 8.66 × 10^10^, respectively, and the calculated minimum *r* (the radius of gas bubbles) are 9.14 × 10^−8^ m and 8.01 × 10^−8^ m, respectively. A larger number of activated nucleation sites cause a greater expansion energy of heavy oil, and smaller bubbles are easily trapped in the oil phase. Therefore, foamy oil is more effective and stable. In addition, high depletion rates cause an increase in *S*~gc~. Increasing the depletion rate from 1.5 m^3^/h to 4.5 m^3^/h causes an increase in *S*~gc~ from 3.63% to 5.88%. Thus, the free gas mobility decreases, which causes a delay in the time of gas produced and increase the cumulative gas produced. We know that the pressure depletion rate not only influences the early time stage of gas bubble nucleation and growth but also affects the late-time stage of gas flow.Effect of scaleThe effect of scale was investigated by simulating pressure depletion tests in sandpacks of different lengths. Figure [8a](#Fig8){ref-type="fig"} presents the calculated oil recoveries from the sandpacks of different lengths, which shows that the oil recovery increases as the length decreases. The length primarily affects the pressure gradient within the sandpacks. In the short sandpack, a high pressure gradient develops in the entire region of the sandpack. Therefore, the pressure depletion rate for the shorter sandpack is faster than the pressure depletion rate for the longer sandpack (Fig. [8c](#Fig8){ref-type="fig"}), which generates more efficient oil recovery. In addition, the higher pressure depletion rate for the shorter sandpack causes fast gas bubble nucleation, growth and gas produced. This finding explains why the cumulative gas produced in the shorter sandpack is higher than the cumulative gas produced in the longer sandpack (Fig. [8b](#Fig8){ref-type="fig"}).Figure 8Effect of sandpack length on foamy solution-gas drive processes: (**a**) oil recovery; (**b**) cumulative gas produced; (**c**) average pressure.Effect of solution GORTable 4Affecting factors investigated in this study and their effects on foamy solution-gas drive processes.Affecting factorsRange of affecting factorsChange of factors (times)Change of oil recovery (%)Change of cumulative gas produced (%)Depletion rate (cm^3^/h)1.52.2533.754.5345.7953.28Scale (length) (m)0.250.3750.50.6250.75325.89122.33Critical gas saturation0.01250.018750.0250.031250.0375321.4624.35Oil viscosity (mPa•s)2086.53129.7541735216.256259.539.528.02Solution GOR (m^3^/m^3^)7.811.715.619.523.434.825.86Diffusion coefficient (m^2^/s)1 × 10^−11^5 × 10^−11^1 × 10^−10^5 × 10^−10^1 × 10^−9^1000.640.26

The simulated results of depletion tests with oil of varying solution GOR reveal that a high solution GOR causes low oil recovery, high cumulative gas produced and high pressure. The solution GOR exerts a significant influence on the sandpack pressures. However, the effect of solution GOR on oil recovery and cumulative gas produced is not distinct. When the heavy oil has a high solution GOR, the rate of gas bubble growth is high, and more gas released from solution starts to flow as a continuous phase. As a result, the gas saturation in the sandpack increases at a faster rate. When the gas saturation attains *S*~gc~, the free gas begins production. Therefore, the high solution GOR causes low oil recovery and high cumulative gas produced, which indicates that maintaining the foamy oil flow becomes more difficult when the solution GOR of heavy oil is too high. However, when the heavy oil has a high solution GOR, the high rate of gas bubble growth causes significant expansion of the heavy oil, which is beneficial for maintaining the sandpack pressure.

### Sensitive factors {#Sec18}

The previous discussion of the affecting factors indicates that foamy solution-gas drive processes are very complicated, and can be affected by many factors. However, a systematic and quantitative study has determine the sensitive factors for foamy solution-gas drive processes. In this study, all affecting factors were adjusted with a constant factor of 25%. Thus, the variation in the factors investigated in this study remains constant (Table [4](#Tab4){ref-type="table"}). The maximum values are increased by three times compared with their minimum values. In this manner, a comparison of the performance of these cases can provide information about the sensitive factors for foamy solution-gas drive.

Depending on the extent of change shown in Table [4](#Tab4){ref-type="table"}, the simulated results indicate that the corresponding change in oil recovery for the depletion rate, length, and *S*~gc~, *μ*~o~, and GOR are 45.79%, 25.89%, 21.46%, 9.52%, and 4.82% respectively, and the change in cumulative gas produced for the depletion rate, length, and *S*~gc~, *μ*~o~, and GOR are 53.28%, 122.33%, 24.35%, 8.02% and 5.86%, respectively. Therefore, the oil recovery is sensitive to the depletion rate, length, and *S*~gc~, and the cumulative gas produced is sensitive to the length, depletion rate, and *S*~gc~. Oil viscosity and solution GOR are not sensitive factors for foamy solution-gas drive processes. Although the values of the diffusion coefficient increased by 100 times, this coefficient has a minimal effect on foamy solution-gas drive processes (the change in oil recovery and the change in cumulative gas produced are 0.64% and 0.26%, respectively).

### Model application in M reservoir {#Sec19}

The M reservoir is located along the southern margin of the Eastern Venezuela Basin. The developed foamy oil model was used to predict the performance of a production platform in the M reservoir in foamy solution-gas drive processes (Fig. [9](#Fig9){ref-type="fig"}). The controlled region of the platform is 3200 m × 1200 m with a pay zone thickness of 200 m. The values of the porosities and permeabilities vary within the range of 3819 mD--9511 mD and 0.14--0.20, respectively. The region consists of three layers with four wells in each layer. The well spacing and row spacing are 600 m. This region is simulated with 3780 gridblocks. Timesteps are automatically controlled with a maximum allowable timestep size of one day. This limit is chosen due to accuracy considerations.Figure 9Results of history matching and performance prediction for foamy solution-gas drive of a typical production platform.

Because the pressure depletion tests conducted in this study were performed to simulate the cold production of this region, the oil sample was obtained from this region, and the foamy oil model parameters for this simulation are directly applied after an experimental match was obtained (Table [3](#Tab3){ref-type="table"}). The only uncertain parameters in this model are geologic parameters. To ensure the accuracy of these parameters, production history matching was conducted using the developed model, which is shown in Fig. [9](#Fig9){ref-type="fig"}.

As shown in Fig. [9](#Fig9){ref-type="fig"}, the predicted oil recovery rapidly increased, and the producing GOR remained relatively stable. Therefore, the predicted results calculated by the developed model can reflect the production characteristics of foamy solution-gas drive process. Thus, the model is effective for simulating the foamy solution-gas drive process for an actual reservoir.

After the history matching process, the performance prediction of the platform was conducted by the model. The results reveal that the ten-year oil recovery of foamy solution-gas drive production was 9.51%. The increased extent of oil recovery gradually reduces, and the producing GOR sharply increased from 66.42 m^3^/m^3^ to 127.65 m^3^/m^3^ in the next ten years because gas bubbles trapped in oil phase coalesce to form free gas with a decrease of reservoir pressure, which causes the disappearance of foamy oil mechanism. The foamy oil flow will become a two phase flow of oil and gas in the reservoir in the next ten years. Therefore, enhanced oil recovery methods are needed to extend the foamy solution-gas drive process and increase the oil recovery of the platform.

Conclusions {#Sec20}
===========

Three typical stages for foamy solution-gas drive can be employed to determine whether a heavy oil exhibits substantial foamy oil tendencies. Higher depletion rates generate in more activated nucleation sites, smaller gas bubbles, and higher *S*~gc~, which benefits foamy solution-gas drive processes.High oil viscosity is beneficial in foamy solution-gas drive processes. Oil viscosity affects foamy solution-gas drive primarily by decreasing gas mobility rather than inhibiting the growth of gas bubbles.A high pressure gradient develops in the short sandpack, which enable high oil recovery and simultaneously causes fast bubble nucleation, bubble growth and gas produced in the shorter sandpack.The solution GOR exerts a strong influence on sandpack pressures. However, the effect of solution GOR on oil recovery and cumulative gas produced are not distinct.Foamy solution-gas drive processes are sensitive to the depletion rate, length, and critical gas saturation. Oil viscosity, solution GOR and diffusion coefficient are not sensitive factors.Our experimental results are based on ideal conditions. The real practice of developing oil in heavy oil reservoirs requires further analysis. In addition, the new developed model was only verified by pressure depletion tests and a typical foamy oil reservoir. Additional studies are needed to verify the validity of the model in more complex conditions.

**Publisher\'s note:** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

This work was financially supported by National Natural Science Foundation of China (No. 51604293), Shandong Provincial Natural Science Foundation, China (No. ZR2016EEB30), the Fundamental Research Funds for the Central Universities (No. 17CX02009A), Qingdao Applied Basic Research Program (No. 17-1-1-32-jch) and Scientific Research Foundation of China University of Petroleum for Talent Introduction (No. YJ201601093).

X. Sun designed the work, analyzed the results, and wrote the manuscript. Y. Zhang revised the manuscript. Z. Song and P. Li conducted all experiments. S. Wang developed the new model. C. Wang prepared the figures and tables. All authors discussed the results and commented on the manuscript at all stages.

Competing Interests {#FPar1}
===================

The authors declare no competing interests.
