Relevance of analytical Buckley-Leverett solution for immiscible oil displacement by various gases

In order to generate the valid numerical simulation model, the sufficient amount of gathered data from the oil field is required. However, it is not always possible to acquire such data at the initial stage of project development. Buckley and Leverett (1942) developed the analytical solution allowing to assess the oil displacement efficiency. One of the main assumptions of this model is incompressibility of oil and injected fluid. For slightly compressible water and oil such assumption is rational. However, that is not always the case when the gas is injected. This research aims to identify the conditions at which the usage of the incompressible gas model is appropriate. Likewise, the cases when the model of compressible gas is required are also evaluated. To accomplish the goals of this research, the comparative analysis between the injection of compressible and incompressible gases was undertaken using the numerical solution of the correspondent reservoir engineering problem. The validation of the numerical model was undertaken showing that it matches the analytical Buckley-Leverett solution. The findings of this research indicate that the relative and absolute density change with the pressure of the injected gas has the profound impact on the convergence between two models under consideration. With the increase in the injection pressure, the discrepancy between the models of compressible and incompressible gas raises for all the considered injection fluids (CO2, CH4 and N2). Due to a steep slope of 'density-pressure' curve for CO2 at low initial reservoir pressure, the incompressible model cannot accurately predict the oil displacement efficiency by this gas at any reasonable injection pressure. All 1D results are also representative for 2D simulation. However, the mismatch between two models increases considerably for 2D simulation scenarios.


Introduction
The extraction of hydrocarbons from oil and gas reservoirs is a complex process which depends on fluid and rock properties as well as on reservoir driving mechanisms. Primary oil recovery, in which natural reservoir energy is used to produce oil, is followed by the secondary stage where the reservoir pressure is typically maintained by water or gas flooding. Albeit, the injection of such fluids might commence from the beginning of the field life span due to technical and economic reasons (?).
According to Craig and Bray (1971), the gas injection projects can be traced to 1917, whereas Lake (2007) reported the beginning of immiscible lean hydrocarbon gas flooding in the US from the 1930th. A number of gases are currently used for pressure maintenance. The most widespread among them are nitrogen (N 2 ), carbon dioxide (CO 2 ) and lean hydrocarbon gas (mainly CH 4 ) (?). The gas is typically injected into the overlying oil interval gas cap or into the oil column. In some cases CO 2 might be injected into deep saline aquifers in order to reduce greenhouse gas emissions (??).
The fundamental manuscript Mechanism of Fluid Displacement in Sands was presented by Buckley and Leverett (1942). This work represents the analytical solution of two-phase immiscible incompressible fluid transport through the porous medium allowing to estimate the dynamics of reservoir fluid displacement by the injectant. The proposed model was based on several assumptions. Among them are the incompressibility of all reservoir fluids, and the absence of capillary and compositional effects. The initial pressure and saturation in a reservoir have to be constant, and the injection flow rate has to be steady.
The further development of Buckley-Leverett model was offered by Welge (1952) who found the analytical solution for a saturation shock front. These studies assume that the oil displacement by injected gases can be estimated with sufficient accuracy, assuming gas incompressibility. This concept is also supported by several other research papers, such as Kern (1952), Shreve and Welch (1956), Pirson (1977), etc.
When the complexity of the oil reservoir does not correspond to the analytical Buckley-Leverett solution, it is more practical to use the numerical reservoir simulation techniques, allowing to implement complicated development conditions. However, in some cases express analyses are required, for instance, as an alternative for comparison with numerical simulation results. In the circumstances when the sufficient production data is not yet available, some methods allowing to estimate the oil displacement efficiency are vital.
Importantly, there is an approach allowing to estimate the oil displacement by compressible gases analytically (Bedrikovetsky 2013). One chapter of the work conducted by Professor Bedrikovetsky scrutinises the effects of compressibility on two-phase displacement providing valuable findings. However, the implemented analytical methods are significantly more sophisticated and, consequently, less appropriate for the first simple estimation of the oil recovery than the aforementioned Buckley-Leverett solution. Moreover, using such analytical models of the compressible gas injection, it is impossible to specify the classical boundary conditions, such as the constant injection rate and fixed initial reservoir pressure at the production well.
Thus, a comparative analysis of numerical compressible and incompressible gas injection models is essential. The aim of the paper is to evaluate at which circum-stances the Buckley-Leverett solution, or incompressible gas numerical model generally, might be applicable for accurate estimation of the oil displacement efficiency. In order to scrutinize the research questions, the injection scenarios of most widely used gases will be simulated at different reservoir and injection pressures.
Clearly, the numerical diffusion can to some extent affect the results of numerical calculations. If this is a case, it can be inappropriate to compare such numerical solutions with analytical models. Hence, the two numerical models accounting for gas compressibility and incompressibility are compared under the scope of this research, where the incompressible model represents the numerical analogue of the analytical Buckley-Leverett solution.
When the incompressible gas is considered, its density is held constant. If the average or maximum reservoir pressure is used for density estimation, the saturation profile depends on the distance between the production and injection wells. As this distance increments, the maximum and average reservoir pressures raise accordingly. Thus, as the first gas reached the bottom hole of the injection well, the density value for the incompressible case was fixed at the initial reservoir pressure. Such approach is also suitable as the reservoir pressure is constant around the production well in accordance with the Buckley-Leverett analytical model. The rest of the paper is structured as follows. Section 2 describes the physical and mathematical representation of research objectives in the differential and numerical forms. Likewise, the geometrical characteristics of 1D and 2D models are also presented. Section 3 lists the results of the numerical analysis, while section 4 scrutinises and discusses the main findings of the research. Subsequently, the main conclusions are summarised in section 5.

Methods
A physical model represents immiscible transport of multiphase (oil and water) flow through porous media. The valuable assumptions include the absence of capillary and gravitational forces, incompressible oil, constant viscosity, as well as homogeneous and isotropic reservoir rock. The aforementioned phenomena can be described by the following system of equations (1 -5) (Muskat et al 1937;Barenblatt et al 1960;Dullien 2012;Bear 2013).
where j = g (gas) , o (oil), φ -porosity, ρ -density, Ssaturation, υ -velocity,q -density of source mass flow, k r -relative permeability, k -absolute permeability, µdynamic viscosity, P -reservoir pressure, t -time, Vvolume, Ω -surface. The formulae (1) and (2) represent the law of conservation of mass and Darcy's equation, respectively. The algebraic expressions (4) and (5) describe the applied initial and boundary conditions. When the numerical solution is applied, the variable splitting (P and S g ) is suitable and it is convenient to convert the equations (1) and (2) into the following forms (Chen et al 2006): The model of a block-centred geometry is used for finite difference representation of the integral equations (6) and (7) (Kazemi et al 1976;Aziz and Settari 1979;Jamal et al 2006;Chen et al 2006;Smith 1985;Ames 2014): where q -mass flow rate. Operator ∆Ω denotes the summation over all surface elements ∆Ω of a grid block. Other parameters which do not contain such operator describe the cell itself. Boundary conditions (5) can be presented in the following form: The applied numerical solution scheme refers to the Sequential solution method (SEQ-method) (Aziz and Settari 1979). The parameters containing unknown variables in the numerical model (8 -10) are given by (11,12): where sign '+' or '−' as a subscript in the pressure symbol (P ± ) means the orientation of the grid cell relative to the particular surface element ∆Ω. If the pressure subscript is '+', the direction of the particular grid cell is positive; otherwise, it becomes negative.
The averaging of any parameter, for instance ρ, indicates arithmetic mean except the relative permeability k r for which the upstream weighting (Aziz and Settari 1979) is used. Averaging is performed between two cells that share a common surface element ∆Ω.
The coefficients in the numerical model (8 -10) are given by equations (13, 14): As the comparative analysis is undertaken between two numerical models (compressible and incompressible gases), it is critical to minimise the discrepancy between the analytical Buckley-Leverett model and its numerical analogous. As analytical Buckley-Leverett solution includes the conservation of initial reservoir pressure around the production well, the numerical model has to obey this condition (15).
Such condition might be obtained when analysing the equation derived by Buckley and Leverett in their model (Buckley and Leverett 1942). This equation does not contain reservoir pressure explicitly. However, using the Darcy and diffusivity equations, the appropriate condition for reservoir pressure might be obtained. Physically, the Buckley-Leverett model implies the immediate production of the fluid reaching the bottom hole of the production well. Clearly, this is only possible if the condition (15) is obeyed.
To the best of authors knowledge, the available and relevant reservoir engineering software does not allow to fix the reservoir pressure in the particular grid block. Therefore, all algorithms of numerical simulation were implemented using C++ programming language.
Production rates can be obtained by the following equations (Aziz and Settari 1979) : where f g = krg µg krg µg + kro µo -gas fractional flow.
In order to conduct the quantitative analysis, the defined 1D and 2D grid geometry and constant rock properties were selected and can be found in Table 1.

Results
This section reveals the results obtained by using the numerical simulation model discussed above. Two sets of identical calculations were conducted: 1D and 2D immiscible gas displacement of oil. Subsequently, numerical models of compressible and incompressible gas injection scenarios were simulated. For both 1D and 2D cases, varying parameters included injected gas (CO 2 , CH 4 and N 2 ), as well as initial and maximal reservoir pressures (see Table 2). The validation of the applied numerical model was generated to ensure the convergence between the analytical Buckley-Leverett solution (Buckley and Leverett 1942) and its analogous incompressible numerical model (see Fig. 1). The irreducible oil and gas saturations were absent for the validation scenario to exclude the numerical diffusion around the injection well. Such approach is correct as the validation case aims to prove the principal appropriateness of the numerical solution.
Using Peng-Robinson equation of state (Peng and Robinson 1976), 'density vs. pressure' curves were built for all gases. This ensures that gas compressibility, and thus substantial density variations with pressure, are incorporated into the model (see Fig. 2a). Using equations of the straight line, 'density vs. pressure' relationships were further approximated for all the injected fluids. The linear approximation of density change with pressure was chosen as all the trend lines appropriately matched with approximated curves with insignificant errors. Likewise, such approach allowed to determine and clearly visualise two distinctive zones at the CO 2 'density vs. pressure' curve.
The approximations are depicted in the form of dashed lines in Fig. 2a and listed in Table 2. Knowing that the compressibility is a function of density, one can clearly see (Fig. 2a) how density of the particular gas changes with pressure when the compressible model is considered. Additionally, Corey correlation exponents (Corey 1954) were taken from the research conducted by Parvazdavani et al (2013), generalised and used as an input for numerical flow simulator for all gases under consideration (see Fig. 2b).

1D calculations
Figures 3 -6 highlight several most demonstrative 1D results. The calculation scenarios were divided into several categories.
Firstly, the CO 2 injection case was simulated at a low initial reservoir pressure and small injection pressure, due to a steep slope of CO 2 'density vs. pressure' curve in this range. For more details on pressure interval for this scenario see label '1' in Fig. 2a and in Table  2. According to Fig. 3b, the significant gap between compressible and incompressible saturation fronts can be seen for such CO 2 injection case. Such difference between saturation fronts also affects the well flow rates, leading to the increase in accumulated production of oil, which is considerably higher for incompressible model (see Fig. 3a).
Secondly, the injection of three different gases (CO 2 , CH 4 and N 2 ) was simulated at large initial reservoir pressure with small injection pressure (e.g. for CH 4 case see Fig. 4). For more details on pressure range for this calculation, see number '2' in Fig. 2a and in Table 2. The obtained results provide the minor, although visible, gap between the compressible and incompressible saturation fronts for all of the injected gases. Subsequently, both well oil rates and accumulated production values indicated close compliance between compressible and incompressible models of all the injected gases.
Finally, large injection pressure scenarios were simulated for all the gases at high initial reservoir pressure. For more details on pressure diapason for this scenario, see number '3' for CO 2 and '4' for CH 4 in Fig. 2a and in Table 2. Carbon dioxide injection scenario even at large injection pressure showed proximity between the progress of incompressible and compressible immiscible fronts towards the production well. The same trend can be observed by looking at closely matching values of the well oil rate and its accumulated production (see Fig. 5a).
As can be seen from Fig. 6, the incompressible saturation front significantly outstrips the compressible saturation profile for the CH 4 case. Likewise, a substantial difference of well production rate and accumulated production profiles can also be observed, especially before the gas break through. Similar results were obtained for the not depicted N 2 immiscible injection case.  Table 2). (b) Set of relative permeability curves and their range for oil-gas system. nd * -not graphically depicted calculations, but their results are discussed.

2D calculations
As mentioned above, calculation scenarios for 2D geometry repeated 1D scenarios. For the sake of visibility, two out of seven 2D results are presented in this paper. Figure 7 represents the similar calculation scenario as was undertaken for 1D CO 2 (see Fig. 3) case at the low initial reservoir and injection pressures, while Fig. 8  repeats the 1D CH 4 injection at the high initial reservoir and injection pressures. It is clear that substantial difference between the progress of compressible and incompressible gas fronts towards the production well was also preserved for 2D cases. Similarly, well oil rates and its accumulated productions were also considerably low for compressible cases in comparison with incompressible scenarios. Considering the 1D and 2D cases, one can observe the higher divergence between compressible and incompressible cases for 2D reservoir geometry.

Discussion
The two numerical models of immiscible oil displacement by gas injection have been compared. The first of them considered injected gas as the incompressible fluid mimicking the Buckley-Leverett model, while the second accounts for gas compressibility. Any outstanding characteristics identified during the numerical simulation process are discussed below.

Gas 'density -pressure' relationship
To consider gas compressibility, Peng-Robinson equation of state was utilised allowing the construction of 'density vs. pressure' relationship for three most widely used for immiscible oil displacement gases, namely CO 2 , N 2 and lean hydrocarbon gas (CH 4 is assumed). The range of considered reservoir pressures varied from 5 to 40 MPa which corresponds to current field practice in the oil and gas industry (Taber et al 1997). Examining the obtained 'density vs. pressure' plot, several conclusions were made allowing to identify calculation scenarios.
According to Fig. 2a, CO 2 has a steep density increase at the pressure range between 5 and 10 MPa. However, as the pressure approaches 10 MPa , CO 2 'density vs. pressure' curve flattens significantly, until the end of the analysed pressure range. The trend is different for CH 4 and N 2 cases, where the slope of 'density vs. pressure' curves is fairly constant along the whole pressure interval.
Thus, the assumption was made that both CH 4 and N 2 'density-pressure' relationships might be approximated by the equation of a straight line. In turn, two such equations are required for the CO 2 injection case. For more details see Table 2 and Fig. 2a. To construct 'density vs. pressure' curves, geothermal gradient (Lake 2007) and hydrostatic water head data was used to calculate temperature and pressure values for different elevation depths. The results were utilised to obtain gas compressibility factor allowing to calculate gas density using an equation of state. However, in the circumstances when primary recovery precedes the immiscible gas injection, reservoir pressure can be significantly reduced. In such case, reservoir temperature is still consistent with geothermal gradient. Thus, for the particular reservoir pressure, temperature values can vary, affecting gas density. The possible density variation at particular pressure was also taken into consideration and does not affect the main qualitative results. Such variation range is depicted as yellow blurred diapason at Fig. 2a.

Relative permeability curves
Relative permeability data was constructed using Corey correlation exponents taken from the research conducted by Parvazdavani et al (2013). The decision to use this particular data was made on the basis of experimental focus of the paper, as well as the broad range of Corey exponents for reservoirs with different fluid and petrophysical characteristics. As indicated by Ghoodjani andBolouri (2011) andAl-Menhali et al (2015), the shape of relative permeability curves does not considerably change for various gases in a particular reservoir, instead the irreducible oil saturation values differ frequently. Such peculiarity allowed to construct generalised gas/oil rel-ative permeability curves which are valid in the range depicted as yellow blurred area at Fig. 2a.

Injection scenario 1
As can be seen from scenario 1 (see Table 2 and Fig. 3), CO 2 injection at a low initial reservoir pressure leads to a significant discrepancy between compressible and incompressible numerical models even at a low injection pressure. Such phenomenon is characterised by substantial absolute and relative change in CO 2 density at pressure range between 5 and 10 MPa (see '1' at Fig. 2a).
Significant lag in the saturation progress of compressible gas front affects the well oil production profile. As is evident from Fig. 3a, the well oil flow rate in the compressible case is lower, comparing to the incompressible numerical solution. We can state that flow rate is directly proportional to the velocity of the saturation front. The well oil flow rate is constant for incompressible case before the gas breaks through. Subsequently, the well oil flow rate drops after breakthrough, as the well produces both oil and gas.
The situation is different for compressible gas where immiscible gas front accelerates during the progress towards the production well. It should also be noted that flow rates almost become equal at the later stage, as several pore volumes were injected into the oil reservoir.

Injection scenario 2
As was described in the results section, the second scenario implied the pressure maintenance by all three gases at the high initial reservoir and low injection pressures. Importantly, the N 2 and CO 2 injection scenarios were not depicted in this article. Former due to its close similarity to CH 4 case, and latter because of its perfect match between compressible and incompressible solutions. Thus, all the qualitative results obtained for CH 4 are also applicable to N 2 . As can be seen from Fig. 2a, both absolute and relative changes in density are small for all the injected gases at the range between 20 and 23 MPa. As a result, the compressible and incompressible solutions for such scenarios differ insignificantly, though some variations can be observed, for instance, in the case of CH 4 injection (see Fig. 4).

Injection scenarios 3 and 4
Scenarios 3 and 4 (see Table 2) simulate the injection of CH 4 and CO 2 at the high initial reservoir (20 MPa) and large differential (15 MPa) pressures (see Fig. 6 and 4, respectively). While CO 2 injection scenario still provides close matching between incompressible and compressible numerical models, this is not the case for CH 4 immiscible flooding. Such discrepancy occurred due to the difference in the slope for 'density vs. pressure' curve for these gases at the considered pressure range. The relative and absolute changes in the density of CO 2 are still not sufficient to cause any significant difference between compressible and incompressible saturation fronts. For displacement of oil by CH 4 injection, saturation front for incompressible gas overrides the compressible model for more than 200 meters before the break through. The results of N 2 injection are highly consistent with the CH 4 case and are not presented graphically in this research.

Injection scenarios 1 and 4 (2D)
Scenarios 1 and 4 were chosen for 2D simulation as the ones showing the highest discrepancy between the compressible and incompressible cases. Importantly, the 2D calculations are qualitatively consistent with their 1D analogues. However, some refinements should be noted. There is a much larger mismatch between the compressible and incompressible simulations in terms of oil flow That clearly leads to the even higher gap between the compressible and incompressible saturation fronts. Thus, there is even larger difference between the accumulated productions of the two models. We can also conclude that mismatch between the compressible and incompressible models will increase even more significantly for 3D simulations.

Conclusion
The concluding section revisits the research questions raised in the introductory part, summarizes the main results and findings of this paper, and provides implications and conclusion based on these findings.
The main aim of this study was to explore the possible application of Buckley-Leverett model for the oil displacement by immiscible gases. The broader focus of this paper was on the comparison between the compressible and incompressible gas injection models in the oil fields. In order to scrutinise the research questions, the numerical model was created allowing to repeat the analytical Buckley-Leverett solution, as well as to simulate the injection of compressible fluids. The 'density vs. pressure' curves were built for the most widely used injection gases (CO 2 , CH 4 , N 2 ). Following the qualitative examination, these curves were used as an input parameter to the reservoir flow simulator.
Conducted flow simulations have allowed to summarise the main research findings as follows.
1. The 'density vs. pressure relationships for CH 4 and N 2 might be approximated using the equation of a straight line. However, two such equations are required to approximate CO 2 density due to its nonlinear shape. 2. The assumption that at high reservoir pressures the results from compressible and incompressible gas behaviour are similar, while there is a significant difference for lower reservoir pressures is frequently made looking expectable and trivial. However, as was revealed in the research, this may be or may not be correct depending on the injected gas or the reservoir pressure variation during the gas flooding process. 3. The injection of CO 2 at a low initial reservoir pressure cannot be correctly predicted by assuming the gas as incompressible fluid, even at low injection pressure. The discrepancy between compressible and incompressible cases increases significantly as injection pressure increments. A similar trend can be observed when oil is displaced by CH 4 and N 2 at a high injection pressure. The initial reservoir pressure does not affect the results in this case. Hence, the reservoir performance cannot be accurately predicted using Buckley-Leverett solution, for the cases mentioned above.
4. The injection of N 2 and CH 4 at a low injection pressure and either high or low initial reservoir pressure leads to sufficient match between compressible and incompressible injection scenarios. Considering high initial reservoir pressure and low CO 2 injection pressure, perfect compliance between the two numerical solutions is observed. Additionally, pressure maintenance using CO 2 at high injection pressure does not affect significantly the convergence between compressible and incompressible models at high initial reservoir pressure. Albeit, the mismatch between the two models is higher in comparison with oil displacement by CO 2 at a low injection pressure. 5. Both absolute and relative density change with pressure affects the results of the simulation. As these values increment, the difference between compressible and incompressible models raises respectively. Thus, these attributes should be taken into consideration when the immiscible gas injection is simulated using Buckley-Leverett model. 6. The velocity of compressible gas saturation front increases in the reservoir as it progresses towards the production well. Such phenomenon is uncommon for the Buckley-Leverett solution where saturation front moves with constant velocity. 7. Two-dimensional simulation allows to qualitatively prove 1D results. However, the difference between compressible and incompressible solutions increases in the 2D case in comparison with 1D calculations. it would not have been possible to complete this work. It should also be noted that the C++ template library for linear algebra 'Eigen' (Guennebaud et al 2010) was very helpful and convenient for the solution of linear algebraic equations. This is a pre-print of an article published in Journal of Petroleum Exploration and Production Technology. The final authenticated version is available online at: https://doi.org/10.1007/s13202-018-0516-6 .