A new mathematical model and experimental validation on foamy-oil flow in developing heavy oil reservoirs

To model foamy-oil flow in the development of heavy oil reservoirs, three depletion experiments were conducted with foamy oil treated as a pseudo-single-phase flow. In this pseudo single phase, dispersed bubbles are viewed as a part of the oil, and the redefined effective permeability varies with the changes of pressure depletion rate, oil viscosity, and gas saturation. A mathematical expression for the effective permeability was developed based on experiments, where the viscosity of foamy oil is assumed to be approximately equal to the saturated oil under equivalent conditions. The compressibility coefficient of foamy oil is treated as a volume-weighted compressibility coefficient of that of oil and gas phases. A new mathematical model for foamy-oil flow was proposed with consideration of foamy-oil supersaturation. To validate the mathematical model, the oil recovery and the production gas-oil ratio (GOR) calculated by the new model, conventional black oil model, supersaturation model and pseudo-bubble-point (PBP) model were all compared with those of the experimental data. The new model provided a substantially improved fit to the experimental data compared with the rest three models, which verifies the suitability of the mathematical model presented for simulating foamy-oil flow in the development of heavy oil reservoirs.


Experiments
This experimental part has two functions, the first one is to help establish a mathematical expression of foamy-oil relative permeability, in which the effective permeability was measured through experiments, and the second function is to help validate the accuracy of the model we proposed in this paper by comparing the recovery and production gas-oil ratio (GOR) of experiments with the new model and the other three models which have been proposed by other researchers. Experimental apparatus. Figure 1 illustrated the equipment used in the laboratory depletion experiments.
An intermediate container is used to transfer the saturated oil sample from the sample container to the long sand pack. The long sand pack is constructed from a special rough-walled steel tube of 5.0 cm diameter and 80.0 cm 1-pump working fluid; 2-ISCO pump; 3-valve; 4-sample container; 5-thermostatic system; 6-oil volume gauge; 7-gas volume gauge; 8-back pressure pump; 9-long sand pack; 10-six-port valve; 11-intermediate container.
length with a maximum burden pressure of 32.0 MPa. A back pressure pump with a maximum operating pressure of 70.0 MPa is used to control the pressure of the outlet of the sand pack. An International Standard Classification of Occupations (ISCO) pump is used for driving the fluid into the sand pack and for pressurizing the system at a constant pressure or constant flow rate. The maximum flow rate is 25.0 ml/min and the maximum pressure is 70.0 MPa. An electronically controlled thermostat is used to ensure a constant temperature in the range of room temperature to 200.0 °C. A metering system including gas and oil volume gauges is used to record the amount of gas and oil production. A vacuum pump is used to initialize the sand pack to the designated pressure at the beginning of the experiments.
Material. The design of sand packs used in the experiments were according to the properties of layers in MPE-3 Block in Venezuela, the sand packs were comprised of clean quartz sand, and the compressibility coefficient was 6.5 × 10 −4 MPa −1 . The other parameters are listed in Table 1.
Two heavy oil samples were prepared for the experiments. The first oil sample (No. 1) was crude oil derived from the MPE-3 Block in Venezuela, and the second oil sample (No. 2) was a mass fraction of 90% Venezuelan crude oil mixed with 10% kerosene. Methane was used as the dissolved gas in these experiments. Under the designated conditions of temperature at 50.5 °C and pressure at 12.2 MPa, the specific parameters of the oil samples are shown in Table 2.
According to the formation water properties, the water used in the experiments was a NaHCO 3 type in which the total salinity was 12,500 ppm, the concentrations of HCO − 3 and Cl − were 2,450 ppm and 10,350 ppm, respectively, and the PH value was 7. 35 The experimental procedures were conducted as follows: (1) Before starting the experiments, the sand pack was initialized from vacuum to the designated pressure and then saturated with brine water. The saturated water volume was measured for determining the porosity. The absolute permeability was measured by injecting water at different rates through the sand pack. (2) Live oil (the oil contains dissolved gas) was prepared in the sample container by mixing the methane gas and dead oil to the desired solution gas-oil ratio (GOR) at the designated pressure. Prior to mixing, 5-10 mm diameter glass beads were placed in the sample container to accelerate the dissolution of gas in the oil. The inlet was closed after injection and the mixture was stirred at a fixed rate at constant temperature for about 10 to 30 days to ensure homogeneous dissolution of the gas in the crude oil. (3) The intermediate container was initialized from vacuum to the designated pressure and the saturated oil was subsequently transferred thereto at high pressure. The oil sample was injected from intermediate container into the sand pack from the inlet end at a rate of 1.0 ml/min until the sand-pack pressure reached the desired value. Then, the oil and irreducible water saturation were calculated. The sand pack was allowed to equilibrate for at least 24 hours before the experiment was initiated. (4) During each depletion experiment, the valve at the injection end was kept closed while the pressure at the production end was decreased at a constant rate under the control of the back pressure pump. The cumulative oil and gas productions were measured at regular time intervals. Each experiment was halted when the outlet did not produced oil any more.  (5) The oil recovery was the result of the cumulative oil production divided by the total oil within the sand pack, whereas the production GOR was the result of the cumulative gas production divided by the cumulative oil production.

Mathematical Models
In this section, a mathematical model will be established to simulate the development of foamy oil reservoirs. To understand to whole process easily, a workflow was made to illustrate how we develop this model in Fig. 2.

Assumptions.
(1) Only oil and gas phases are present, and foamy-oil flow follows Darcy's law under reservoir conditions. Fluid flow is isothermal. (2) All released solution gas remains entrained in the oil phase until the reservoir pressure drops to the pseudo-bubble-point pressure. Below the pseudo-bubble point pressure, free gas begins to form. (3) The effects of the capillary force acting between liquid and gas phases, and gravity are both ignored. (4) The liquid phase is slightly compressible while the gas phase is strongly compressible.
Determination of the key parameters of the model. The effective permeability in porous media. Because foamy oil is regarded as a pseudo-single phase, and the dispersed gas bubbles are viewed as a part of the oil phase. The bubbles are assumed to flow with the oil 23 . Devereux (1974) assumed that the permeability of a porous media decreases as the fraction of droplets in the emulsion increases, and implied that the permeability reduced with decreasing flow rate and with an increasing drop-size to pore-size ratio 37 . Soo and Radke (1986) also proposed a similar theory 38 . These theories are suitable for foamy oil as well, where the high pressure depletion rate causes that the large bubbles in the oil are divided into smaller bubbles; thus, the higher the pressure depletion rate, the less the permeability reduces 39 . A higher crude oil viscosity ensures that the bubble diameter can be stabilized at a smaller scale, corresponding to a less influence on the permeability reduction in porous media 23,40,41 .
Based on the above analysis, the volume fraction of the dispersed bubbles (or gas saturation), viscosity of crude oil, and the pressure depletion rate have a significant influence on the permeability in porous media. The relationship between the effective permeability of sand pack and above influence parameters will be built in the subsequent process [42][43][44][45][46] .
During the depletion processes of experiments A, B, and C, the redefined effective permeability of the sand pack was measured at different outlet pressure when the pressure was less than the pseudo-bubble-point pressure. When the outlet pressure equals pseudo-bubble-point pressure, the gas saturation is called critical gas saturation (CGS), which is defined as the minimum gas saturation at which gas-phase flow can occur.
Here, redefined relative permeability K r is defined as Where, K r is dimensionless permeability, decimal; K is absolute permeability, μm 2 ; K eff is the effective permeability of foamy-oil-flow in the porous media, μm 2 . Figure 3 shows the relationship between redefined relative permeability and gas saturation based upon the results measured of the above three depletion experiments (experiments A, B, and C).
The redefined effective permeability was calculated according to Darcy's law. The gas saturation was calculated by the Henry's law and the real gas state equation. Henry's law means the dissolved gas oil ratio increases linearly with the increase of pressure, especially in heavy oils, this regulation is much more accurate [47][48][49] . Therefore, the gas saturation near the outlet end was obtained according the outlet pressure, and the redefined effective permeability was computed by the pressure gradient, fluid viscosity, and the flow velocity at the outlet end. Notably, the redefined effective permeability and gas saturation in Fig. 3 represented the corresponding values at outlet ends of sand packs in three experiments. The results of these three experiments also coincide with the analysis above. The corresponding relative permeability increases with the increase of viscosity of oil sample and depletion rate when the gas saturation is smaller than the critical gas saturation (CGS).
From Fig. 3, if K r is a function of the pressure depletion rate, oil viscosity, and gas saturation when the gas saturation is less than CGS, the oil viscosity and pressure depletion rate are in proportion to K r , and the gas saturation is in inverse proportion to K r . To clarify, we define a dimensionless parameter β: Where: β is a dimensionless parameter, dimensionless; α is unit conversion factor, 6 × 10 10 /MPa 2 ; p d is pressure depletion rate, 10 −1 MPa/min; μ o is the oil viscosity, mPa · s; S g is the gas saturation, decimal; S gc is the gas critical saturation, decimal. As shown in Fig. 3, K r increasingly approaches a value of 1 as the value of β increases. As such, K r can be expressed as Where, K r is relative permeability, dimensionless; χ is a constant which need experimental fitting, dimensionless.
It can be proven that, when χ·β tends to infinity, the value of K r approaches 1, indicating that the effective permeability of foamy-oil flow in the sand pack approaches that of the absolute permeability. The value of χ can be obtained from the experimental results through trial and error. Thus, the CGS and χ obtained in experiments A, B, and C are shown in Table 3.
The viscosity of foamy oil. The apparent viscosity of the foamy oil is an important parameter in modeling the behavior of foamy-oil flow. Application of emulsion theories suggests that the foamy dispersion viscosity should be higher than the continuous phase viscosity 38,49 . However, such theories, when concerned with foamy-oil flow applied directly to porous media, remain questionable 22 . Researchers have evaluated foamy oil viscosity using different methods, where the results are higher or lower than saturated oil viscosity under equivalent conditions 2, [50][51][52][53][54][55] . Alshmakhy and Maini used three different viscosity measurement techniques to measure the apparent viscosity of foamy oil, saturated oil, and dead oil 55 . The researchers found the viscosity of foamy oil to be similar to that of saturated oil over a large range of gas saturation. In the present paper, it is assumed that the viscosity of foamy oil is equal to that of saturated oil under equivalent conditions, and this assumption also fitted well with  Table 3. The critical gas saturation (CGS) and χ obtained in experiments A, B, and C.
the experimental measure. So the saturated oil viscosity is a function of the pressure p, the foamy oil viscosity is therefore obtained by the saturated oil viscosity at different pressures according to the following relationship: Where, μ fo is foamy-oil viscosity, mPa · s.
The compressibility coefficient of foamy oil. The compressibility coefficient of foamy oil is composed of two partsthe oil and the gas compressibility coefficients 5 . It complies with the linear volume weighted average compressibility coefficients of oil and gas when the gas concentration is much higher or lower than that of oil 56 . According to the data in above three experiments, the maximum gas saturation was less than 12% when the foamy-oil phenomenon occurred. Therefore, the value of the foamy-oil compressibility coefficient c fo , can be expressed as where, c g at a temperature T is given as The total compressibility coefficient c t is the sum of the rock compressibility coefficient c f and c fo : where, c fo is the compression coefficient of foamy oil, MPa −1 ; c o is the compression coefficient of oil phase, MPa −1 ; c g is the compression coefficient of gas phase, MPa −1 ; c t is the total compression coefficient, MPa −1 ; c f is the compression coefficient of rock, MPa −1 ; S o and S g are the saturation of oil and gas, respectively, decimal; Z is the compressibility factor, dimensionless; T is temperature, K; p is the pressure, 10 −1 MPa.

Derivation of the mathematical model. Modeling the fluid flow equation.
Emulsion has been typically treated as a single phase in previous studies 37,38 . Similarly, for foamy-oil flow, foamy oil is herein treated as a pseudo-single phase. The mass conservation equation can be expressed as follows: where, ρ fo is the density of foamy oil, g/cm 3 ; v fo is the velocity of foamy oil, cm/s; φ is the porosity of porous media, decimal; ∇p is the pressure gradient, 10 −1 MPa/cm; t is time, s. Eq. (10) is the general formula of foam oil flow in heavy oil reservoirs. Actually, the conventional black oil model is its specific expression.
(1) When the reservoir pressure is higher than bubble-point pressure, K eff becomes the absolute permeability K, μ fo becomes μ o , and the oil compressibility coefficient is very small, the second term on the left side of Eq. (10) can be ignored. The model can be simplified as conventional black oil model. (2) When the reservoir pressure is between bubble-point pressure and pseudo-bubble-point pressure, the order of magnitude of the foamy-oil compressibility coefficient is very large and closed to that of gas, and, therefore, the second term on the left side of Eq. (10) cannot be ignored. It represents foamy oil flow in heavy oil reservoirs. (3) When the reservoir pressure is lower than pseudo-bubble-point pressure, it is the oil-gas two-phase flow; the model becomes dissolved gas-drive process in black oil model 57,58 .
In Eq. (10), K eff is given by Eq. (1) and Eq. (3), μ fo is given by Eq. (4), and c fo and c t are given by Eq. (5) and Eq. (7), respectively. In Eq. (10), the two variables p and S g (from Eq. (5)) require an auxiliary equation to establish their relationship.

Establishment of the auxiliary equation.
During the process of foamy-oil flow stage, the production gas-oil ratio (GOR) is approximately equal to the solution GOR. Therefore, the velocities of the dispersed bubbles and the oil are equivalent, and the amount of total gas in the control volume (solution gas plus free gas) is a constant. Based on this fact, an equation of state has been derived 28 . It is supposed that, under a given pressure and temperature, a sealed control volume contains an oil volume V o and a free gas volume V g , given as According to Henry's Law, the amount of solution gas in the control volume, given as the solution GOR R si , is proportional to the pressure: Eq. (12) can be expressed as another form where, the subscript 'sc' refers to the volume at the standard pressure (0.1 MPa) and temperature (20 °C). V gd is the volume of gas that the system drops from a certain pressure and temperature to the standard state; V o is the oil phase volume at the standard state, cm 3 . The real gas law is applicable; the number of moles of the gas can be expressed as where, N is the number of moles of the gas, mol; R is gas constant, J/ (mol · K), given as 8.314. Under any pressure and temperature conditions, the sum of the number of moles of gas bubbles and solution gas equals the total moles of gas. Applying the standard pressure and temperature, and replacing V g by V gd , Eq. (14) can be used to calculate the number of solution gas moles at the bubble point pressure p b , and using Eqs. (12) and (13)  where, is a solubility factor, dimensionless; which is given by Eq. (17) gives the relationship between S g and p, and, when p equals the pseudo-bubble point pressure, it can be used to calculate the CGS. In fact, when p decreases, the foamy-oil supersaturation phenomenon always occurs, which means the solution gas in the oil phase is more than the value calculated by equilibrium thermodynamics. However, Eq. (17) does not take into account the supersaturation phenomenon 36 . With the generation of bubbles, the system gradually is restored to a balanced state 12 . The dimensionless pressure and dimensionless gas saturation can be defined as follows.
When the supersaturation phenomenon was taken into account, the two dimensionless parameters were subtracted to define the following supersaturation equation 36 .
Where, p* is dimensionless pressure, dimensionless; S g * is dimensionless gas saturation, dimensionless. We consider the case where the system pressure is suddenly reduced from p b to a lower value. Because the formation of bubbles requires a finite length of time, S g = 0 at the beginning, and, hence, S g * remains initially zero while p* > 0 of Eq. (20). Therefore, if f (p*, S g *) is positive, the system is supersaturated. The supersaturation pressure decreases as time increases, whereas the rate of decrease is proportional to the supersaturation and inversely proportional to the time; thus, the process can be represented by Scientific REPORts | 7: 8534 | DOI:10.1038/s41598-017-08882-2 where, τ is the time required to transform from a non-thermodynamic equilibrium state to a thermodynamic equilibrium state, s. It varies with the properties of the crude oil and the relevant operating conditions, which can be determined by experiments.
When the system is in a supersaturated state, the right side of Eq. (21) is a positive value, where the gas saturation begins to rise, and Eq. (21) describes the variation of the foamy-oil flow process with respect to time. Eq. (21) and Eq. (10) together with the initial and boundary conditions can constitute solvable partial differential equations, which represent the new mathematical model for foamy-oil flow. The corresponding initial conditions are, i and the boundary condition is where, L is location of the sand pack, cm; p i is initial pressure in sand pack, 10 −1 MPa.

Solving method.
If the reservoir pressure is greater than the bubble-point pressure and less than the pseudo-bubble-point pressure, the black oil model is employed for these two cases. When the pressure is higher than the bubble-point pressure, It is very easy to obtain numerical solution by the finite-difference method to solve Eq. (10) (the second term on the left side of Eq. (10) can be ignored). When the pressure is lower than the pseudo-bubble-point pressure, the model becomes dissolved gas-drive process in black oil model. To solve this stage of the model, it is necessary to establish two equations of oil and gas, and the relation between these two equations is that gas saturation plus oil saturation equal unit. For gas equation, Perrine-Martin approximation can be used to address the square term of the pressure derivative. Therefore, the model can be solved by coupling gas and oil equations. If the reservoir pressure is between the bubble-point pressure and pseudo-bubble-point pressure, the new model can be employed for this case. Eq.

Results and Discussion
Experimental results. Measurements of the cumulative oil and gas production were employed to calculate the oil recovery and production GOR, and these two production indicators were used to evaluate the extent of agreement between the models and experiments. The results in experiments A, B, and C are shown in Table 4.  Model validation. The depletion experiment of the sand pack can be divided into three stages.
The first stage represents single-phase oil flow. In this stage, both the oil recovery and the production GOR are very low. The sand-pack pressure is greater than the bubble point pressure estimated by the black oil model (the model we used is a commercial numerical simulator named as Computer Modelling Group (CMG), which was developed by a Canadian company) in this stage.
The second stage is foamy-oil flow. This stage represents the major oil-production stage of the depletion experiment. The sand-pack pressure lies between the bubble point pressure and the pseudo-bubble point pressure. The production GOR was slightly higher than the solution GOR. The foamy-oil recovery was found to increase with the increase of pressure depletion rate and oil viscosity. To better verify the simulation accuracy, the new model, the PBP model 26 , the supersaturation model 36 , and the conventional black oil model (CMG) were used for comparison with the data of experiments A, B and C. Based on experimental results, the relaxation time τ = 6 × 10 5 s, 4.7 × 10 5 s, and 4.2 × 10 5 s in the new model and the supersaturation model for the conditions in experiments A, B, and C, respectively. In the PBP model, the phenomenon of supersaturation was not considered, that is to say, Eq (17), not Eq (21), was used to compute in the model. In order to compare the results more obviously, the black oil model was also employed in this stage.
The third stage is the oil-gas two-phase flow. Upon completion of the foamy-oil flow stage, the sand-pack pressure is lower than the pseudo-bubble point pressure, and the gas phase becomes a continuous phase and begins to flow. The oil recovery was low and the production GOR increased sharply, the black oil model (CMG) was used to simulate the production performance in this stage.
Analysis and discussion. The parameters used to solve the models were obtained from the depletion experimental data. From Figs 4(a), 5(a), and 6(a), the oil recoveries simulated by the new model agree well with the experimental data, except in the late stages. This is because in the actual flow process, the possibility of bubble coalescence and the formation of narrowing flow path channels increases with the increase of gas saturation and results in reduced oil mobility, the simulated recoveries are slightly higher than that indicated by the experimental data. However, the PBP model does not take the supersaturation phenomenon into account, so the gas saturation in foamy oil is higher than the real value. It has two opposite effects, on the one hand, high gas saturation increases the system's compression coefficient; on the other hand, it reduces the effective permeability of the core pack. The result of the comprehensive effect is that the oil recovery is higher than experimental data in the early stage of the foam oil flow and lower than the experimental data in the late stage. Though the supersaturation model takes the supersaturation phenomenon into consideration, it addresses the foamy-oil flow as oil-gas two-phase flow, which makes the mobility of gas higher than real situation, thereby producing lower oil recoveries and higher production GOR than the experimental data. The black oil model does not take the foamy-oil properties into account, and addresses the foamy-oil flow process as conventional oil-gas two-phase flow, which makes simulation results quite different from the experimental data.
From Figs 4(b), 5(b), and 6(b), the simulated production GOR by the new model and PBP model agree well with the experimental data. The reason why production GOR remain low levels when the pressure higher than pseudo-bubble-point pressure is that the interfacial tension between heavy oil and gas are too strong, which makes the heavy oils capture and restrain the gas as separate gas bubbles and cannot form continuous free gas channel immediately 59 . Finally, free gas will form when the pressure drops to pseudo-bubble point and production GOR will soar to a big value. However, the results obtained from the supersaturation model and the black oil model indicated that the production GOR increased significantly when the pressure was below the bubble point pressure, especially the black oil model (CMG), which is inconsistent with the experimental results.
Taken together, it is readily apparent that the experimental data is in better agreement with the results simulated by the new model than the other three models, which represents we obtain progress in simulating and evaluating foamy oil's performances and characteristics on foamy-oil flow in porous media for heavy oil reservoirs.

Conclusions
Based on the mathematical model proposed and experimental validation in this paper for foamy-oil flow in porous media for heavy oil reservoirs, the following conclusions can be drawn, and the underlying drawbacks of this model are also pointed.
(1) According to the depletion experimental results, the redefined effective permeability of the sand pack varying with the pressure depletion rate, oil viscosity, and gas saturation, and an appropriate mathematical expression for the effective permeability on foamy-oil flow in porous media was developed. (2) It was assumed that the foamy-oil viscosity is equal to the saturated-oil viscosity under equivalent conditions, which was justified by the results of the model. (3) The compressibility coefficient of foamy oil is a volume-weighted of the oil and gas phase compressibility coefficients. The foamy-oil volume dramatically expands between the bubble point pressure and the pseudo-bubble point pressure, which allows significant displacement of the oil. (4) The particular fluid properties of foamy oil are considered the main cause of high oil recovery during the depletion process in heavy oil. The second stage of foamy-oil flow is the main oil-production stage, where high oil recovery is obtained and production GOR is slightly higher than solution GOR. (5) By taking into account foamy-oil supersaturation, a mathematical model that can reflect foamy-oil characteristics was modeled. As verified by three depletion experiments, the simulation results of the new model agree well with the experimental data, while the results obtained from the pseudo-bubble-point (PBP) model, the supersaturation model, and the black oil model (CMG) were not satisfactory. These results suggest that the new model is suitable for application for foamy-oil flow in porous media of heavy oil reservoirs. (6) The determination of some parameters in the model needs experimental assistance, like pseudo bubble point pressure and relaxation time, which requires relatively high accuracy of experiments. High experimental errors may cause wrong results. In addition, the model was only verified by sand pack experiments, therefore, whether it is suitable for reservoir development or not still needs further study.