Numerical Approach of a Coupled Pressure‑Saturation Model Describing Oil‑Water Flow in Porous Media

Two-phase flow in porous media is a very active field of research, due to its important applications in groundwater pollution, CO 2 sequestration, or oil and gas production from petroleum reservoirs, just to name a few of them. Fractional flow equations, which make use of Darcy’s law, for describing the movement of two immiscible fluids in a porous medium, are among the most relevant mathematical models in reservoir simulation. This work aims to solve a fractional flow model formed by an elliptic equation, representing the spatial distribution of the pressure, and a hyperbolic equation describing the space-time evolution of water saturation. The numerical solution of the elliptic part is obtained using a finite-element (FE) scheme, while the hyperbolic equation is solved by means of two different numerical approaches, both in the finite-volume (FV) framework. One is based on a monotonic upstream-centered scheme for conservation laws (MUSCL)-Hancock scheme, whereas the other makes use of a weighted essentially non-oscillatory (ENO) reconstruction. In both cases, a first-order centered (FORCE)- 𝛼 numerical scheme is applied for intercell flux reconstruction, which constitutes a new contribution in the field of fractional flow models describing oil-water movement. A relevant feature of this work is the study of the effect of the parameter 𝛼 on the numerical solution of the models considered. We also show that, in the FORCE- 𝛼 method, when the parameter 𝛼 increases, the errors diminish and the order of accuracy is more properly attained, as verified using a manufactured solution technique.


Introduction
Reservoir simulation was a very active field of research since the beginning of petroleum engineering in the 1930s. Mathematical modeling and numerical simulation play a fundamental role in the analysis of oil reservoir behavior and the study of multiphase flow in porous media. First computer codes for reservoir simulation were developed in the decades of the 1940s and 1950s, when carrying out relevant research in the numerical solution of flow equations. Some details on the historical evolution of reservoir simulation can be found in [12].
Research on multiphase flow in porous media attracts currently great attention, as it can be applied to a number of situations in the field of fluid dynamics, such as groundwater pollution, reservoir simulation, water-oil-gas systems, oil and gas recovery, or CO 2 sequestration in geological structures. Interesting situations commonly encountered are oil-water systems, which constitute the aim of the present work. A great number of works devoted to multiphase flow in porous media can be found in the literature such as the classical references [1,5,19,38,39], or more recent ones as [7,16,26,41].
In this work, we focus our attention on the mathematical modeling and numerical simulation of the underground movement of oil and water. This movement is favored by connected pores and also fractures that are usually present in the geological medium. The typical mathematical models devoted to this kind of applications are systems of parabolic equations with saturation and pressure of the different phases as unknowns. It is worth mentioning here the frequently used black-oil model, see for instance [37] for details, and the so-called compositional models. A comparison of the results obtained with both models can be found in [3], where numerical simulation was performed by means of a discontinuous Galerkin method for mass transport, coupled with an implicit mixed hybrid finiteelement (FE) scheme for computing pressure and velocity fields. A flexible general purpose model for multiphase isothermal or thermal compositional flow simulation in porous media was introduced in [40], where numerical simulation results were validated using classical benchmark problems based on black-oil and compositional models. Most of the numerical schemes commonly used in this context are based on finite-difference and FE numerical methods. We refer the reader to the seminal work [27] for a strong background in reservoir simulation, comprising physical principles, mathematical models, and numerical resolution by finite-difference formulations. There is also some literature on the application of finitevolume numerical methods in the field of reservoir simulation, see for instance [6,11,14,18,20,23].
A wise combination of statistical techniques and streamlines-based methods in oil reservoir applications was described in [25], where a multilevel Montecarlo method combined with a stream-based solver was utilized to perform two-phase Buckley-Leverett transport simulations of an oil reservoir with uncertain heterogeneous permeability. Going ahead with streamline methods for flow simulation in porous media we refer the interested reader to [29] where, in addition, other geological applications are put forward.
In this work, we consider a hyperbolic-elliptic fractional flow model, which is a system with an elliptic equation for the pressure and a hyperbolic equation for the saturation, as described in [8,13,22]. In [8], analytical and numerical study of this type of system was conducted. A local existence and uniqueness study, based on the Arzela-Ascoli theorem applied to oil recovery processes was introduced in [28].
In this work, we solve the elliptic equation using an FE scheme and the hyperbolic equation by means of FV methods, due to its better performance in the presence of 1 3 discontinuities. The pressure gradient is computed at cell interfaces, from the nodal values of the pressure, using piecewise linear interpolation, selecting the stencil that minimizes the absolute value of the slope, as performed in essentially non-oscillatory (ENO) reconstruction. Flow velocities are then obtained by making use of Darcy's law, which relates the volumetric flow rate to the pressure gradient.
Prominent studies and detailed analyses of the processes of flow and transport through porous media, with emphasis on applications of Darcy's law, could be found in the classical reference [5], and also in the more recent one [4].
A relevant feature of this work is the study of the effect of the parameter appearing in the first order centered (FORCE)-scheme for numerical flux reconstruction, in the context of problems arising in oil-water movement in porous media. These FORCE-techniques were introduced in [35] for conservative problems, and after that extended to nonconservative systems in [9]. A deep study of this technique, with a special focus on the parameter and applications to Baer-Nunziato equations for compressible two-phase flow was conducted in [36].
As a first approach, we use in this work the MUSCL-Hancock method (MHM) which is of particular interest to solve the hyperbolic part of this model, that is the saturation equation, since it is a simple, but still effective, procedure and second order accurate in space and time. In order to increase the order of accuracy, we resort to fifth-order WENO reconstruction with the RK3-TVD approach for time integration, which can give better results, especially for long-time simulations.
In Table 1, we summarize the main parameters and variables used in the mathematical model, together with their physical meaning.
The hyperbolic-elliptic model is not commonly used in reservoir modeling, but has several advantages derived from the fact that hyperbolic equations represent adequately the displacement of water and oil and, from the point of view of computational cost, the numerical resolution of this hyperbolic-type equation is less expensive than the usual parabolic models, considering stability restrictions.
The main contribution of this work concerns the use of the FORCE-numerical schemes applied to reservoir simulation and oil-water flow in porous media. This technique allows obtaining a family of numerical methods, depending on the parameter , with good convergence properties. The rest of this document is structured as follows: in Sect. 2, the two-phase flow mathematical model based on an elliptic-hyperbolic coupled formulation is described. Section 3 is devoted to a brief overview of the numerical schemes including the assessment of the numerical accuracy, using a manufactured solution. Section 4 is focused on the numerical resolution of the coupled model, including the extension to a 2D situation. Finally, conclusions and future research are given.

Mathematical Model
We consider in this work a coupled model pressure-saturation, describing the immiscible movement of oil in direction of the production well. This movement is favored by the injection of water into the reservoir, usually done either to keep the pressure in the reservoir or to drive the oil towards the production well. Saturation stands for the fraction of the pore volume which is occupied by a particular fluid (oil, gas, or water). The oil can flow if a critical saturation is exceeded, otherwise, the oil remains in the pores and is not able to flow. In practical applications, usually, the oil is displaced from the reservoir by injecting water or gas. Nonetheless, some amount of oil will remain without being displaced, which gives rise to residual oil saturation.
The addition of saturations of wetting phase, s w , and the non-wetting phase, s n , must verify s w + s n = 1 . We note that the non-wetting phase usually refers to oil, while the wetting phase commonly refers to water. The usual way to identify wetting and non-wetting phases is by measuring the contact angle between the fluid and the solid surface. When > 90 • phases are non-wetting, whereas phases with < 90 • are wetting. See a sketch in Fig. 1.
The flow velocity of a particular phase is obtained through Darcy's law, which in the case of multiphase flow adopts the form where the relative permeability of phase i is given by k ri . K represents the total permeability, corresponding to rock permeability which, in this model, has a value of 300 mD (millidarcies). These data come from previous works comprising the study of the medium and rock-fluid characterization performed by geologists, geophysicists, and petrophysicists.
We consider in this work a hyperbolic-elliptic model describing the oil-water movement in the porous medium. Despite there are some references on this type of mathematical formulation, such as [8,22] in oil reservoir simulation or [24] in surfacesubsurface flow, it is not very frequent in the literature. The model we are dealing with is • pressure equation (elliptic) • saturation equation (hyperbolic) In (2), (s w ) stands for the total mobility of both phases combined, that is, while Q t is a specific volumetric injection/production rate term given by Q t = with s cw being the connate water saturation, which represents the amount of water adsorbed on the surface of the grains of the rock divided by the pore volume, and s or is the residual oil saturation. The value of the exponent W depends upon the model considered, see [2] for different choices. In Sect. 4 of this work, details are given on the values of the parameters involved. Fig. 1 Water as wetting phase, oil as non-wetting phase. When > 90 • , the phase is non-wetting, while if < 90 • , the phase is wetting We have also introduced in (3) v w representing the flow velocity given by Darcy's law. Equations (2) and (3) are coupled through the mobility, (s w ) , which is a function of the saturation, as well as the pressure, P, that depends on the flow velocity.

Numerical Schemes
The mathematical model described in the previous section is solved in the present work using a finite-element method (FEM) for the pressure equation, and two different FV formulations for the saturation part, which are based on a MUSCL-Hancock approach and a WENO5-RK3TVD scheme. A FORCE-numerical scheme for intercell flux reconstruction is used in both FV formulations. Most of the details on these numerical techniques are skipped here, only some description of the FORCE-scheme is given, since this is a new contribution in reservoir simulation.
As it is well known, Godunov's theorem (see [15]) states that linear numerical schemes for solving partial differential equations, having the property of not generating new extrema (monotone scheme), can be at most first-order accurate. According to this theorem, in order to obtain monotone schemes of order higher than 1, nonlinear numerical schemes must be implemented. One of the numerical techniques used in this work is the MHM, which is nonlinear and second-order accurate. In this case, we use the second-order MHM technique, although we remark that, according to [33], the MUSCL approach allows the construction of very high-order schemes. Further, the WENO methodology improves the order of accuracy of the numerical schemes developed.
In order to perform the numerical resolution of the coupled model, the spatial domain is discretized into subintervals S j = [x j , x j+1 ] , which are finite elements for the case of the pressure equation and control volumes when solving, by the FV method, the saturation equation. Hence, the nodes used for the pressure are intercell boundaries for saturation. The process followed, at each time step, starts by solving the elliptic part, that is the pressure equation, by means of an FE scheme to yield nodal values of the pressure P j , at the nodes with abscissae x j . In this FE approach, linear two-node P 1 elements are used.
Next, at the same time step, an FV scheme is applied to solve the saturation equation, and the nodes, x j , are taken now as cell interfaces of the control volumes. This is useful since the values of the velocity can be straightforward computed at cell interfaces from the nodal values of the pressure according to Darcy's law. These values of the velocity are needed to solve the saturation equation. Since we just have point values of the pressure at nodes locations, it is necessary to perform a reconstruction process to approximate the derivative | . Thus, the velocity of water at the cell interface j can be obtained according to the formula By integrating (3) over the control volume S j and dividing by the length of the control volume, Δx j , as well as assuming is constant in the control volume, we have where s w = s w (x, t) and v w = v w (x, t); Kk r m j .
x j where s w,j = s w,j (t) and v w,j = v w,j (t) with s w,j the cell average of the saturation within cell S j which is given by Equation (9) can be solved using an efficient ODE solver. In this work, we use the RK3-TVD scheme, although different solvers could be applied. We remark that numerical schemes attaining arbitrary high order of accuracy, such as the ADER approach, are very efficient techniques to solve this type of problems. We refer the reader to the relevant publications [30,34] for the classical ADER methodology applied to the numerical solution of hyperbolic problems and also to [10,17] for the application of an interesting variant of the ADER technique, named local space-time DG (LSTDG), specially focused to problems with stiff source terms. As aforementioned, the FORCE-numerical technique is used in this work to carry out flux reconstruction at cell interfaces. Some details on this technique are given in the next subsection.

FORCE-˛ Numerical Flux
The FORCE numerical scheme was introduced by Toro [31,32] in the context of hyperbolic equations. It was proposed as a deterministic version of Glimm's or random choice method. Since it is a centered scheme, it does not make use of the eigenvalues of the Jacobian matrix, as it happens in other classical numerical fluxes, such as upwind or Rusanov's ones. The extension to multi-dimensional problems was proposed in [9,35,36] where the parameter was introduced to take into consideration several spatial dimensions. The FORCE-numerical schemes emerged then as a family of numerical methods, -dependent, that were successfully applied to different hyperbolic problems. The FORCE-intercell numerical flux, in the context of 1D fractional flow models reads for the left interface, where , as indicated above, is a real number, usually related to the number of space dimensions, that is = 1 for 1D, = 2 for 2D, and = 3 for 3D. However, any value of can be used. In the reference [36], there is a complete study on stability, monotonicity, and numerical viscosity of FORCE-schemes taking values of up to 200. Numerical results given in that reference indicate that the CFL number must decrease as increases, in order to ensure stability. This conclusion is also attained in the present work.
The expression (11) is an arithmetic mean of the version of the Lax-Friedrichs numerical flux which reads and the version of the Lax-Wendroff numerical flux, which is obtained by means of where

Assessment of the Order of Accuracy of the Schemes
We focus our attention on the assessment of the order of accuracy of the numerical schemes. We consider the following problem for the evolution of water saturation: whose solution is s w (x, t) = cos(x(3 − x)t) and F s is a forcing term. The velocity is set here to v w = 0.5 and the fractional flow, f w (s w ) , is obtained according to (4) which takes into account the mobilities of both water and oil.
The problem is solved using the numerical schemes MUSCL-Hancock and the RK3-TVD WENO5, in order to assess their convergence and order of accuracy. Flux reconstruction is fulfilled via the FORCE-technique and different values of the parameter are considered. Results of the application of MUSCL-Hancock approach are displayed in Fig. 2 where both the L 2 -error norm and the numerical orders of accuracy are shown for the values = 1, 2, 3, 14 . It can be observed that the theoretical order of accuracy is successfully attained. Also, we notice that, as the value of increases, the error norm tends to be smaller.
As for the WENO5-RK3TVD approach, convergence results are depicted in Fig. 3, showing that, as the value of the parameter increases, the third order of accuracy is more clearly attained (right frame) and also the error norms become smaller and the evolution of the norm L 2 of the error tends to improve for higher values of (left frame).

Numerical Examples
We recall that the problem under consideration is the fractional flow oil-water hyperbolicelliptic model given by (2) and (3). In these first examples, we omit the source term to deal with the homogeneous system. We use the values s cw = 0.25 , s or = 0.2 and the exponent W = 4 in the case of water and W = 2 when dealing with oil and = 1 . See [21] for details (12) f LF, on the Brooks and Corey model related to water and oil permeabilities, connate saturation, and effective saturation. Details on the physical principles of this type of models can be found in [2]. In Fig. 4, water saturation and oil saturation profiles are shown for different values of the parameter . The velocity has been assumed constant and its value is v w = 4.0 . The show that the solution improves when the value of increases, as the rarefaction (transition zone) and the shock wave (representing the flood front) are better resolved. However, the price to pay is the size of the time step, since it must be reduced, as increases, which agrees with the conclusions reported in [36]. The size of the time step taken here, in order to ensure the stability for all cases considered is Δt = 10 −3 . In Fig. 5, x-t contour plots of water (left frames) and oil (right frames) are depicted for the case of constant velocity v w = 4.0 and for two different values of the parameter namely = 1 (top row) and = 14 (bottom row). The results show that the shock wave representing the flood front appears more clearly when the value of is higher. We now consider the situation where the velocity field is computed from the pressure drop in the reservoir. Dirichlet boundary conditions are assumed for pressure computation, namely P(0)=390, P(50)=186. Results for water and oil saturation are displayed in Fig. 6. It can be verified that the profiles are better defined as the value of the parameter increases.
In Fig. 7, x-t contour plots of water (left frames) and oil (right frames) are depicted for the case of velocity obtained by Darcy's law from the pressure gradient and, as before, for = 1 (top row) and = 14 (bottom row). In Fig. 8, fractional flow function (left frame), water and oil relative permeabilities (right frame) are depicted; numerically obtained for variable velocity, = 14 and time t = 3.

Two-Dimensional (2D) Extension
As an extension to the multi-dimensional case we consider the following problem, written in the divergence form: with suitable boundary and initial conditions. In (16), Q t and Q w are source terms, related to oil and water injection/production rates. Hereafter, 2D Cartesian domains are considered. The numerical scheme applied is the third-order RK TVD method with the FORCE-  s w (x, y, t) = (1 + t)x 2 (x − 1) 2 (y − 1) 2 y 2 for water saturation. At each particular time step, we solve the pressure equation (17) by the FEM using quadrilateral Q 1 elements (see Fig. 9), hence, nodal values of the pressure are computed. We remark that, although the equation for the pressure is elliptic, it is coupled with the water saturation equation via the mobility, (s w ).
The forcing term F P = 200(1 − x 2 − y 2 )exp(x 2 − y 2 ) appears due to the manufactured solution considered. The nodal values of the pressure allow to compute the velocity at FV cell interfaces, applying Darcy's law. Once the velocities are obtained, their values are introduced in the water saturation equation, (18), which is solved by the FV technique, where v w,x and v w,y are the components of water velocity in x and y directions, respectively. The forcing term appearing in (18), F s , is readily obtained using some software of symbolic computation (Fig. 10). The value of the parameter may have an influence on the behavior of the solution, as depicted in Fig. 11 using 20 × 20 control volumes. The results displayed confirm that higher values of the parameter give rise to a most accurate solution. This effect is even more evident, when the number of cells is reduced. We remark that this result was previously put forward in the 1D case.
Here we give a short description on the FORCE-technique in 2D, due to the relevance of this procedure in the current application. More details on this technique can be found, for instance, in [35,36]. The 2D FORCE-scheme applied to a Cartesian mesh, for a control volume indexed as (I, J) reads    in the space-time domain (x, y) ∈ (0, 1) × (0, 1), t > 0 . As detailed before, the spatial grids are staggered for pressure and water saturation. The rectangle displayed in the solid black line in Fig. 9 represents the spatial domain. Concerning the oil volumetric factor its value is B o = 1.235 6 and the oil flow rate Q o = 54 . The coefficient K refers to the medium permeability, which is taken here as K = 300 mD. In this particular case, we take as the water flow rate Q w = 0 . The Dirac's delta at the RHS of the pressure equation is introduced to locate the depression in the upper right coordinate of the domain. The boundary conditions used for pressure computation are P(x, 0, t) = 420 for 0 ⩽ x ⩽ 1 and P(x, 1, t) = 210 for 0 ⩽ x ⩽ 1 , P(0, y, t) = P(1, y, t) = 210 for 0 ⩽ y ⩽ 1 . The initial condition for water saturation is s w (x, y, 0) = 0.8 for 0 ⩽ y ⩽ 1, 0 ⩽ x ⩽ 0.5 and s w (x, y, 0) = 0.2 elsewhere, while boundary conditions are s w (x, y, t) = 0.8 for y = 0, 0 ⩽ x ⩽ 1 and Neumann homogeneous in the rest of the boundary. Results for pressure distribution are displayed in Fig. 12, whereas water (left plot) and oil (right plot) saturations are shown in Fig. 13. As a final example, we modify problem (24) introducing an injection of water Results are shown in Fig. 14, where an increase in water saturation around the injection location is visible and, consequently, a decrease in oil saturation appears.

Conclusions
In this work, we have used different numerical schemes to solve a fractional two-phase flow model in porous media, describing oil-water movement. The model considered here is a hyperbolic-elliptic system, where the elliptic equation represents the distribution of pressure and the hyperbolic part describes the evolution of the saturation front. Both equations are coupled via the mobility which is a function of saturation. This model differs from the typical one based on two parabolic equations and it is also different from the classical Buckley-Leverett model. In order to obtain the numerical solution of the coupled model, an FE approach is used for the elliptic part (pressure equation) and FV techniques are used to solve the hyperbolic part (saturation equation). The FV schemes implemented are a second-order MUSCL-Hancock approach and a WENO5-RK3TVD technique with a FORCE-scheme for intercell flux reconstruction. The numerical results are assessed by means of a manufactured solution technique. The theoretical order of accuracy is attained in practice and the results have revealed that, in the WENO5-RK3TVD FORCE-approach, for higher values of the order of accuracy is more firmly obtained using higher values of the parameter . Also the error norms are lower for higher values of . Extension of the numerical methods to a bi-dimensional situation is also performed for the coupled model pressure-saturation, including a problem with exact solution, aimed to validate the proposed model.
Acknowledgements This work has been carried out under the project: Development of a numerical multiphase flow tool for applications to petroleum industry, with the Spanish oil company REPSOL.
Funding Open Access funding provided thanks to the CRUE-CSIC agreement with Springer Nature.