This work uses Computational Fluid Dynamics (CFD) to simulate the two-phase flow (oil and water) through a reservoir represented by a sandbox model. We investigated the influence in the flows of water having higher and lower mobilities than oil. To accomplish this, we also developed a dedicated solver, with the appropriated equations and representative models implemented in the open-source CFD OpenFOAM platform. In this solver, the black-oil model represented the oil. The results show that the Buckley–Leverett water-flood equation is a good approach for the three-dimensional flow. We observe that the water wall front is mixed to some extent with the oil and evolves obeying an exponential law. Water with mobility lower than oil is not common. However, in this case, the oil recovery is improved and the amount of injected water is reduced. The results comparing different mobilities show that a careful economic assessment should be performed before the field development. We have shown that the low water mobility can increase, as in this studied example, the water front saturation from 0.57 to 0.73, giving a substantial improvement in the oil recovery. The reservoir simulation can provide all process information needed to perform an economical assessment in an oil field exploration.

## I. INTRODUCTION

Because of the high investment required in the discovery of new reserves and the costs involved in the exploration, prediction tools became very important for the oil industry. The increasing advance in high-performance computing permitted Computational Fluid Dynamics (CFD) to obtain numerical solutions of the equations that model the different systems in an acceptable time interval, highlighting computer simulations as one of these tools. Several methodologies^{1} related to reservoir management applications have been developed. Most times, however, it is difficult to apply them to real reservoirs possessing many unknowns and particular features. In these cases, the real reservoir is replaced by a reference model with known properties and down-scaled to a laboratory level.

Simplified models, which are a reasonable approach to the complex reality of an oil reservoir, permit the exploration of the fundamentals of the oil flow inside the rock matrix. It also contributes to a better understanding of the influence of the control parameters on the oil recovering process. As proposed by Bear,^{2} we chose the reference model called the sandbox because we can easily recognize the concepts of geometric similarity, kinematic similarity, and dynamic similarity.

Oil reservoir simulation has been cited in many works^{3} as a tool for reservoir evaluation and management to minimize costs and to increase process efficiency. At present, we can find several commercial softwares dedicated to oil reservoir simulations such as STARS of CMG, ECLIPSE-100 of Schlumberger, REVEAL of Petroleum Experts, and UTCHEM from Texas University at Austin. We used the CFD software OpenFOAM^{4} based on the finite volume method of discretization. The main variables here are the pressure and the saturation of the wetting fluid, calculated by the implicit pressure and explicit saturation (IMPES) method,^{5} where water pressure is the implicit variable and water saturation is the explicit variable.

The water-flood equation from Buckley and Leverett^{2,5–9} was made for a one-dimensional flow, but it is a good approximation for a three-dimensional domain. In addition, in a water flood with the water mobility less than that of the oil, we have a substantial reduction of the amount of the injected water, considering the same quantity of recovered oil.

The models used for representing permeabilities and capillary pressure are of fundamental importance for investigating the mobility’s influence in the immiscible displacement of the oil. In this work, we selected for the permeabilities the model from Stone,^{10} as described in Sec. II. This model also applies to the three-phase flow. For the capillary pressure, we used the empirical correlations from Ref. 5.

Enhanced oil recovery aiming to remove the oil still in place by changing the mobility ratio between water and oil has been pursued as a response to economic demands. Many of these efforts were dedicated to how to improve the oil secondary recovery in what is called Chemical Enhanced Oil Recovery (CEOR).^{11} Since a few decades ago until recent years, we can find many works dedicated to developing surfactants and polymers as agents for CEOR.^{12–14} Some others are proposing Microbial Enhanced Oil Recovery (MEOR),^{15} and in more recent developments, the usage of nanoparticles has been the subject of extensive works,^{11,16–20} most of them based on experimental results.

In Sec. II, we give a brief review of the porous media for a two-phase model, and we present the parameters used in the sandbox model. To end this section, we present the equations implemented by the solver in the OpenFOAM software. In Sec. III, we present the results of the numerical simulations. Initially, we consider one- and two-dimensional models and calculate the time series of the saturation. A plot is shown as well as snapshots for times before and after the water-flooding. The step-like appearance of the saturation function of time is discussed in further detail later. Several results compare the flow behavior for water having mobilities lower and higher than oil. Discussing the consequences of that behavior as well as a deeper view of the waterfront pushing oil are presented in Sec. IV.

## II. THEORETICAL REVIEW

### A. The porous media model

The mass continuity equation for a fluid of density *ρ* through an isotropic porous medium is

where *ϕ* stands for the porosity, *q* is the specific discharge, meaning the rate of the flow per unit area (*Q*/*A*), and **u** is the flowing fluid velocity.

Darcy’s equation^{2} is a linear response transport relation in porous media,

which, in general, is valid as long as the Reynolds number based on an average grain diameter does not exceed some value between 1 and 10. In this equation, *K* is the permeability of the medium, ∇*p* is the pressure gradient, *g* is the gravity, and ∇*z* is the gradient of the *z* coordinate.^{2} Substituting this equation into the continuity equation, we obtain

Equation (3) governs the single-phase fluid flow. In the two-phase flow, the two fluids entirely fill the porous volume. The saturation of the wetting fluid and the non-wetting fluid corresponds to the volume percentage of each fluid in the total porous volume. Thus, if $Sw$ is the percentage of the wetting fluid (water) and *S*_{o} is the volume percentage of the non-wetting fluid (oil), we have

From now on, we use subscripts *w* for water and *o* for oil. Both fluids have a critical saturation, called irreducible, or a residual saturation that is the minimum saturation value that could be achieved for each one. So, we name subscript *r* as the endpoint of the corresponding variable. Then, $Srw$ and *S*_{ro} are endpoint saturations of water and oil, respectively.

Relative permeabilities^{10} are defined according to the following equations:

In these equations, $krowo$, $krogo$, $krwo$, $krgo$, and the exponents $row$, *rog*, $rw$, and *rg* are empirical parameters. They correspond to the endpoint of the two-phase oil relative permeability flowing with water, endpoint of the two-phase oil relative permeability flowing with gas, two-phase oil exponent when the other phase is water, two-phase oil exponent when the other phase is gas, endpoint relative permeability of phase water, and endpoint relative permeability of phase gas. The residual saturations are $Srw$ for water, $Sorw$ for oil, and *S*_{Lrg} for the total liquid. The calculated parameters are the two-phase water relative permeability, the two-phase oil relative permeability, the water relative permeability, the gas relative permeability, and the oil relative permeability. Despite these equations are applied for three-phase permeabilities, in our case, as there is no gas flowing, *S*_{g} = *S*_{rg} = 0. In Ref. 10, the results from these equations are compared with the experimental data from Corey, with the following assumptions:

The water permeability during the oil/water flow is the same as the oil permeability during the oil/gas flow.

The oil permeability during the oil/water flow is the same as the gas permeability during the oil/gas flow.

The two-phase relative permeability parameters such as residual saturations, endpoint relative permeabilities, and exponents are based on the two-phase data of Corey. References 21 and 22 propose correlations for capillary pressure and permeabilities based on experimental data. In this work, we used the correlation from Ref. 5 The results were more consistent when compared to those in Ref. 23. Below, we show the equation used to obtain the capillary pressure correlation. In this model, the equation depends only on water saturation,

In Eq. (10), $pcow$ is the capillary pressure between the oil and water, and in Eq. (11), $pcowmax$ and $pcowmin$ are the maximum and minimum limits for $pcow$. By replacing the appropriate variables in Eq. (1), we may write the continuity equation for the wetting and non-wetting fluids,

After an algebraic manipulation, the time derivative of *S*_{o} can be eliminated. Then, the continuity equations divided by $\rho w$ for the wetting phase depend on variables $pw$ and *S*$w$,

where the mobilities are defined as the ratio between the relative permeability and the viscosity,

### B. Model scale

The dimensions of a real oil reservoir are very large, ranging in a hundred meters. To simulate the flow within the reservoir, one would expect to consider real dimensions, rock geophysical parameters, and fluid properties. Considering this scenario, the simulation runtime can be too long. To overcome this problem, models can be built to reduce the runtime to a few hours. In a scaled model, reservoir dimensions, fluid properties, and rock properties are scaled for the laboratory model so that the ratios of various forces in the reservoir and the physical model are the same.^{24} Physical models can be classified into two categories: (a) scaled model and (b) elemental model. In a scaled model, reservoir dimensions, fluid properties, and rock properties are scaled for the laboratory model so that the ratio of various forces in the reservoir and the physical model are the same.^{7} Therefore, it is essential to manipulate the order of magnitude of the variables involved to reduce the time required for the simulations. A crucial point is to seek a non-dimensional number shared by the model and the real reservoir. Similar to the approach of Brooks and Corey,^{23} we arrived at two dimensionless parameters that should be equal for both the model and the real reservoir, provided that the variables *T*_{0}, Δ*p*_{0}, *K*_{0}, *ρ*_{0}, *μ*_{0}, and *L*_{0} are given in proper scales for the time, pressure, permeability, density, viscosity, and length,

Table I shows the values of these parameters for a typical reservoir.

Parameter . | Value . |
---|---|

Water density | 1000 kg/m^{3} |

Oil density | 800 kg/m^{3} |

Water viscosity | 0.0025 kg/ms |

Oil viscdosity | 0.004 kg/ms |

Gravity | 9.81 m/s^{2} |

Absolute permeabillity | 1000 md |

Porosity | 0.2 |

Initial water saturation, $Swi$ | 0.2 |

Residual water saturation, $Srw$ | 0.15 |

Residual oil saturation, S_{ro} | 0.81 |

Parameter . | Value . |
---|---|

Water density | 1000 kg/m^{3} |

Oil density | 800 kg/m^{3} |

Water viscosity | 0.0025 kg/ms |

Oil viscdosity | 0.004 kg/ms |

Gravity | 9.81 m/s^{2} |

Absolute permeabillity | 1000 md |

Porosity | 0.2 |

Initial water saturation, $Swi$ | 0.2 |

Residual water saturation, $Srw$ | 0.15 |

Residual oil saturation, S_{ro} | 0.81 |

The present simulations were performed in a sandbox model in the shape of a cube with dimensions of 30 × 30 × 30 m^{3}. The injection and production wells were located at the same level at 17 m from the bottom. After some test runs, we concluded that the entire domain should be split in 80 × 80 × 50 cells, in a total of 320 000. This is necessary to have a consistent material balance and a very small change in results with the cell number increasing. The time step was 0.001 s. Using the parameters in Table I, *L*_{0} = 30 m, and considering the result of the time for the breakthrough of *T*_{0} = 60 s for Δ*p* = 300psi, *N*_{r1} = 0, 34 and *N*_{r2} = 4.0 × 10^{−6}. At *t* = 0, $Sw$ = 0.2 and $pw$ = 2.3 × 10^{7} Pa. The pressure bottom hole (pbh) for the production well was set to 2.1 × 10^{7} Pa. For Δ*p* = 400, 350, and 300 psi, the pbh’s for the injection wells were set to 2.376 × 10^{7}, 2.342 × 10^{7}, and 2.307 × 10^{7} Pa, respectively.

### C. Solver for the biphasic immiscible flow

We developed a new solver for the system of coupled differential Eqs. (14) and (15), and we used the CFD software OpenFOAM^{4} to program the solver code. OpenFOAM CFD uses classes in the *C*^{++} language that are developed specifically aimed at this type of application, where a system of coupled differential equations needs to be solved numerically using the finite volume method.^{25} The mathematical solution is calculated by the IMPES method [5]. In Eq. (14), the water pressure is the implicit variable and the water saturation is the explicit variable. We implemented a black-oil model for representing the oil behavior. The calculation was deployed in such a way that, at each time step, all pressure-dependent properties must be updated for the entire domain.

We can define the water fractional flow in the following manner:

Provided that the interfacial tensions between the fluids are small and can be neglected and, in addition, neglecting the gravity effects in Eq. (2), we obtain

The displacement of oil from a porous medium by water plays an important role in producing petroleum.^{7} In both natural water-driven and secondary water-flooding, this process is fundamental. In the case of the one-dimensional flow of incompressible immiscible fluids, we can formulate a simple mathematical description of this process. As the water flows, it creates a front water saturation that is higher than the initial reservoir saturation, and it establishes a discontinuity or an abrupt change in the water saturation.

The fractional flow derivative obeys the following relationship:

*S*_{i} is the initial water saturation of the oil field. The fractional flow derivative related to the water saturation will have its maximum at $Sw$ < $Swf$. Considering a differential element of a unidimensional porous media and applying the continuity equation, we obtain

In this equation, *q*_{t} is the total water injected, *A* is the cross section, and *x* is the distance. This equation tells us how fast a given saturation $Sw$ travels inside the domain. In spite, we have introduced the front displacement concept, depending upon the fluid mobilities, where the flow interface may not be sharp. As shown by Dake,^{8} the most favorable condition is when we have a piston-like flow. For this to be true, we need the flow ahead of the interface to be oil flowing at the initial water saturation or connate water, $Swi$. At this point, we have the highest oil mobility, and we can define

In this case of low water mobility, behind the interface, water flows at the oil saturation *S*_{o} = *S*_{ro}. Thus,

For the flow to be piston-like,

where *M*_{ep} is called the endpoint mobility ratio. If *M*_{ep} ≤ 1, it means that, under an imposed pressure unbalance, the oil can travel in a velocity equal to, or greater than, that of the water. Since the water is pushing the oil, there is no tendency for the oil to be by-passed, which results in the sharp interface between the fluids. However, *M*_{ep} ≥ 1 is more common. Here, water can travel faster than the oil, and then, the latter will be by-passed, leading to an unfavorable water saturation profile.

## III. RESULTS

We have performed simple procedures to validate the model before proceeding to a 3-D simulation. First, we built a domain with only one row and 30 columns representing, as close as possible, a 1-D domain. We assumed the water to have a lower mobility than the oil.

Figure 1 shows the time series of the water saturation at the production well, with an applied Δ*P* of 400 psi. We can see an abrupt change in the saturation after 62-time units, and then, from this moment on, only water is produced. Dake^{8} shows that we can calculate the movable oil volume (MOV) according to the equation

The displacement efficiency *E*_{d} is defined by Eq. (26),

All the contained oil in the reservoir was recovered in that 1D-model by the injection of 1 MOV, considering the displacement efficiency equal to one. In the simulated case, the MOV was 7.9 m^{3}. The total injected water was 10.3 m^{3} with a displacement efficiency of 0.73 and, so, in good agreement with the theoretical MOV. Figure 2 shows domains for the one-dimensional and two-dimensional cases and snapshots of the flow simulations before and after the water-flooding.

The plots of $fw$ vs $Sw$ for one-, two-, and three-dimensional cases are coincident showing that they only depend on the saturation. Figure 3 shows the well-known *S* shape fractional flow, the waterfront saturation having the same behavior in all dimensions. The trivial solution of the Buckley–Leverett equation gives us two possible saturations for each different location. Naturally, this is not possible. To overcome this problem, Pinder^{9} states that such an anomaly arises because of the dependency of *dx*/*dt* on saturation that allows a fast-moving saturation to overtake a slow-moving saturation. To eliminate this contradiction, we require that there must be a shock front at which there is a discontinuity in the $Sw$ function, this way, creating the concept of a shock-front wave.

We performed the three-dimensional simulations with three values of the pressure unbalance Δ*P* (300, 350, and 400 psi) and for water mobility higher and lower than that of the oil. In every iteration, we calculated the variables $pw$, $Sw$, $fw$, $dfw$/$dSw$, and *E*_{d}. $S\xafw$ is the reservoir average saturation and $Swi$ is the initial water saturation. Figure 4 shows plots of the time series *E*_{d} for each Δ*P*, in the cases of high and low mobilities. For high mobility, the plot was performed until finishing the possible recoverable oil. Here, for the same time unit, less oil is recovered than in the case with low mobility, and the simulation lasts much longer than the former case. We show the accumulated and instantaneous oil production in Fig. 5.

We can calculate the waterfront saturation from Eq. (20). Figure 6 shows the comparison of the fractional flow for high and low mobilities. The results are $Swf$ = 0.73 and 0.57 for lower and higher water mobilities, respectively. The former shows a better efficiency displacement because in this case, the water is pushing the oil and not by-passing it, as shown in Fig. 6.

The *S* shape curves for the fractional water flow $fw$ as a function of the water saturation $Sw$ do not depend on the pressure unbalance. However, as Fig. 6(a) shows, they are different for water mobilities lower and higher than those of the oil. The same happens with the derivative $dfw$/$dSw$. Figures 7(a) and 7(b) demonstrate the advantage of the oil recovery on low water mobility. In (a), we see the plot of the oil and water permeabilities. The permeability curves are coincident for low and high water mobilities because they only depend on the saturation.^{10} However, mobilities also depend on the viscosities, as can be shown in Figs. 7(a) and 7(b). If we reduce the water mobility, we have the oil flowing ahead of the waterfront saturation with mobility higher than that of the water, which benefits the recovering oil efficiency.

According to the one-dimensional flow model developed by Buckley and Leverett, a water wall front is created at the injection well, pushing the oil ahead. We have shown that the model developed by Buckley and Leverett gives good results even in a three-dimensional domain. A question arises about how important are capillary and gravity effects. Including these terms requires a great effort calculating the fractional flow derivative. We make a considerable simplification if we consider these contributions negligible,

## IV. DISCUSSIONS

In most cases, water mobility is higher than that of the oil in the reservoir. This brings undesirable effects on the flow because as the water moves faster, it should by-pass the oil, reaching more easily the producing well. However, in the opposite case, if we consider the one-dimensional flow, for instance, we should expect that the water pushing the oil ahead would create a sharp waterfront with a water saturation of $Swf$ = 1 − *S*_{ro}. This was confirmed, at first glance, by Fig. 1. However, expanding the time scale to better observe how the waterfront evolves, what we obtain is shown in Fig. 8. By fitting this plot into an exponential curve, we obtain

These results suggest that there is a region where we should have oil and water mixing to some extent. This is possible if tiny droplets of water and oil coexist in the wall front.

Our simulations show that the water creates a spherical front that is deformed, as shown by the snapshots of Figs. 9(a) and 9(d). As the front approaches the producing well, they form something like a water-spout. When the water-spout reaches the producing well, it establishes the breakthrough and reduces the oil recovery rate. In the very beginning, the water spherical front is deformed by gravity. As it moves on, a spout is formed in the gradient pressure direction. With water having higher mobility, the waterfront travels faster and creates an interface with higher pressure, as shown in Fig. 9(a). In the case of low water mobility, the waterfront travels slower and also at a slower pressure allowing for a more efficient sweeping effect like in Fig. 9(d). In both figures, the contour filter from *Paraview* in the post-processing was applied to the *S*_{wf} level. We can see in Fig. 9(d) that the oil volume displaced by the water is greater than that in Fig. 9(a). Both figures show the moment where the waterfront saturation $Swf$ touches the producing well.

In Fig. 9, we demonstrate how much efficient is the water-flooding at low water mobility. In comparison, Fig. 9(e) exhibits a higher water saturation and a greater swept volume than Fig. 9(b), both at the time of the water breakthrough.

Beyond the breakthrough, we see that, in the case of high water mobility, the simulation shows that after the water creating a preferential path, oil no longer flows into the producing well, and it is almost useless keeping the water injection, as seen in Figs. 5(a) and 5(c). Conversely, it would be expected that the oil flowing lasts much longer in the case of low water mobility, by noting that in Figs. 5(b) and 5(d), oil was still flowing through the producing well at the time that the simulation was interrupted. We made another profile comparison at the endpoints of Figs. 5(a)–5(d). Figures 9(b)–9(f) reinforce what was noted in Figs. 5(a)–5(d).

In a water-flooding with high water mobility, the ratio *oil*–*produced*/*water*–*injected* decreases quickly after the water breakthrough, approaching to zero after a certain time. This ratio decrease is faster as higher is the pressure difference applied. In the water-flooding with low water mobility, the recovered oil is greater because, after the water breakthrough, oil still keeps flowing for a long time. In Figs. 10(a) and 10(b), the ratio *oil*–*produced*/*water*–*injected* is plotted for instantaneous and accumulated oil production.

## V. CONCLUSIONS

Even though greater oil recovery in the water-flooding with low water mobility was achieved, we should perform a very detailed economic assessment before the field development. Depending on the commercial commitment with oil production, we should consider the lowest possible pressure differential. A cost-benefit analysis is anticipated as a tool for providing a suitable configuration for the field. The reservoir simulation can provide all process information needed to perform this analysis. Most of the simulators are commercial packages. In this work, we are using an open-source simulator. We introduced a tool in such a way that we can add any model for the black-oil relative permeability and capillary pressure. This way anyone can freely configure and customize his/her model for fluid and rock properties. OpenFOAM CFD comprises a very large users’ community that is always committed to how to develop and deploy solutions for specific problems. So, by delivering this work, we are contributing with a tool that may be further improved by aggregating new functionalities. As an open-source tool, the proposed solver can be adjusted and parameters can be customized, as well as the capillary pressure and relative permeability models can be fine-tuned, providing a better representation of the studied phenomena.

## ACKNOWLEDGMENTS

A. F. Britto acknowledges INCT-GP for a fellowship DTI-C from CNPq-Brazil.