A Picard-Newton approach for simulating two-phase flow in petroleum reservoirs

In this work, we perform a comparative study of some of the most well-known approaches for solving the system of algebraic equations, obtained by discretizing the governing equations using the Finite Volume Method, for a three-dimensional two-phase (water-oil and water-gas) flow in an oil reservoir. We consider that the flow is isothermal, the fluids immiscible, and we take into account the compressibility of the fluid and the porous matrix. We also use a model of well-reservoir coupling for specified flow rates of injection and production. The solution strategies considered are the Fully Implicit Method, the IMPES Method, the Sequential Method, and a Picard-Newton Method, which represents the main contribution of this work. To illustrate the accuracy of the methods, we considered a two-phase flow in slab geometry, two-phase flow in a five-spot arrangement well, and gas production in a reservoir. For the cases simulated here, the Picard-Newton Method was able to correctly reproduce the flow physics with accuracy comparable to the other three methods.


I. INTRODUCTION
The problems involving transport phenomena in porous media find applications in different areas of scientific and technological knowledge. To name a few examples, we can cite the chemical industry, environmental sciences (transport of contaminants in soil), geology, mechanical engineering, and medicine [9]. This work focuses on a complex problem, which is the flow in petroleum reservoirs [23].
The mathematical models, for fluid flow in porous media, consist of partial differential equations (PDEs). Together with the appropriate boundary and initial conditions (and the constitutive equations), these equations describe the flow through reservoirs [16]. As the governing equations are nonlinear, and the dimensions of the reservoir are on the order of kilometers, we have a complex phenomenon. Also, as the properties generally vary within the reservoir, analytical methods of solution can only be used if physical and mathematical simplifications are assumed. Nev-ertheless, the validity of these simplifications is limited, especially when we consider applications such as real three-dimensional flows [43]. Therefore, for more realistic cases, it is necessary to use numerical methods to solve the governing partial differential equations. In this context, since the 1960s, numerical simulation of oil reservoirs has grown in importance in the oil and gas industry [32].
In the 1950s the first simulators of oil reservoirs appeared, and the petroleum industries started to use computers and numerical analysis as design tools. Nowadays, we use numerical simulation as a tool for management planning, production forecasting, and decision making [24]. The term "numerical reservoir simulation" became commonly used in the 1960s, and we have since applied the numerical simulation as a tool of reservoir engineering to predict the performance of reservoirs under different production scenarios. This branch of engineering deals with oil reserves and the improvement of oil recovery in the oil and gas industry, thus facilitating the decision-making process in financial management [31]. Concerning the numerical simulation, its main objective is to predict reservoir and well pressures, saturation (for multiphase flows), and composition (for a compositional flow) when recovering hydrocarbons.
The nature of the fluid contained in oil reservoirs strongly depends on the stage of the recovery process. At the early stage, called primary recovery, the reservoir generally contains a single-phase, gas or oil, and the production is made by decompression of the fluid and the rock without energy supplementation, by an external source (although in practice this already happens).
This stage ends when the reservoir pressure is no longer sufficient to bring oil to the production wells. Primary recovery usually leaves 70% to 85% of hydrocarbons in the reservoir. A fluid (usually water) is injected into the reservoir (through injection wells) to displace oil toward the production wells aiming to recover some amount of this oil. This procedure serves to maintain the pressure level and the flow in the reservoir. In this process, called secondary recovery, if the reservoir pressure is above the bubble point pressure of the oil phase, there will be two immiscible phases (water and oil). Otherwise, a gaseous phase will also appear and will be in thermodynamic equilibrium with the oil phase. In this case, we often assume that there is no mass transfer between the water phase and the other two phases (oil and gas). Nevertheless, there will be mass transfer between the oil and gas phases [16].
However, unfortunately, the secondary recovery process has been successful in recovering up to 50% of the initial oil amount [16]. Even if there is a large amount of oil in the reservoir, at some point, the production has to be interrupted since secondary enhanced recovery ceases to be economically viable [29]. This situation is aggravated in the case of heavy and viscous oil reservoirs since water is highly mobile [16]. As is well known, in these cases, water tends to flow through the paths of less resistance, forming viscous fingers that can seriously reduce sweep efficiency and compromise the recovery due to the earlier breakthrough [3].
To try to recover the oil that remains in the reservoir after secondary recovery, have been developed several methods. We classify these methods as mi-crobiological, miscible, thermal, and chemical groups. They are known as a tertiary recovery process, and the flow is typically compositional [16].
Currently, even before the decline in production, it is usual to inject water or gas at the beginning of production, aiming to maintain high reservoir pressure [33]. It is noteworthy that, presently, secondary and tertiary terminologies are employed depending on the type of process used. Thus, secondary recovery means the injection of water or gas in the immiscible case, and tertiary recovery corresponds to the other processes [33].
The main objective of this work is the comparative study of some different numerical methods applied to multiphase flow in porous media and a proposal of a new approach. There are three well known traditional approaches to solving nonlinear system of equations for two-phase flows in oil reservoirs [51]: 1. The Fully Implicit (FI) Method (or Newton's Method).
The Fully Implicit method [7,22,36,67,77,56,53,4] is the most widely used in commercial simulators in the oil industry [43]. In this method, transmissibilities, capillary pressures, and source terms (wells) are linearized and all unknowns are obtained simultaneously at each time step [14]. This numerical method is unconditionally stable and allows the use of large time steps, except for precision constraints. However, this method leads us to solve large system of linear algebraic equations resulting in a large computational effort, particularly for highly refined meshes. The computational implementation of the FI method is not a trivial task compared to those of IMPES and Sequential methods.
When capillary effects are not preponderant (weakly coupled equations) or when we use small time steps, the decoupling of equations can lead to numerical methods with better convergence properties [43]. Sheldon et al. [70], and later Fagin and Stewart [34], developed the IMPES method [14,56,45,53], which consists of implicitly solving the pressure equation and explicitly the saturation equation [21]. Capillary pressures and transmissibilities are evaluated explicitly in time n or time n + 1 in a known iteration (n + 1, v). It has the lowest computational cost per iteration due to the decoupling of pressure and saturation equations and the explicit saturation calculation. Have been developed extensions of the IMPES method [45,20,81,76,15,39] and others, based on it, have been proposed [21,10,12] for the Black Oil Model [16,7,44,78,35,42] and the Compositional Model [16,22,41]. However, all these methods suffer from numerical instability due to the explicit calculation of the saturation and are subject to a time step restriction of Courant-Friedrichs-Lewy (CFL) type [18,19,38]. We also mention here the Adaptive Implicit Method (AIM) that employs different levels of implicitness in a different part of the domain. We may choose certain variables to be considered explicitly or implicitly, independently of the choices in other regions, and this procedure may vary from one time-step to the next [22,72,80,37]. Redondo et al. [68] recently discussed the efficiency of the IMPES method for two-phase problems considering high-performance computing applications.
The Sequential method [43,51,2,64,46,50] was first applied to the Black Oil model by Watts [76]. The pressure is calculated in the same way as in the IM-PES method, followed by an implicit saturation calculation to improve the stability properties of the numerical method. Recently, Correa and Murad [25] applied a sequential method on the numerical simulation of three-phase immiscible flow, including geomechanical effects.
As the main contribution of this work, we apply a Picard-Newton approach (named Hybrid Method) to solve two-phase flow in a reservoir. In this method, the pressure is determined implicitly as in the IMPES method, but now saturation is obtained using a fully implicit linearization [43,32,7], a Newton's method approximation. Differently from the recent contribution due to Wong et al. [79], here Newton's method is only applied for the water saturation equation, and the pressure equation is solved using the Picard's method. Moreover, contrary to Wong et al. [79], no particular hypothesis is considered concerning fluid properties along with the iterative solution process.

II. PHYSICAL MODEL
Two-phase flow through a heterogeneous porous medium is modeled by equations that represent the conservation of mass, momentum, and energy. For an isothermal flow, we use only the continuity equation and the modified Darcy's law [7]. We also need to consider that fluid and rock properties are a function of pressure and saturation [33]. Now, we obtain the transport equation for each phase. Thus, by using the modified Darcy's law in the continuity equation we obtain where α stands for the wetting (β) and non-wetting (γ) phases, v is the seepage velocity, k rα is the relative permeability, k is the permeability tensor, µ is the viscosity, p is the pressure, ρ is the density, g is the gravity magnitude, D is the vertical coordinate, φ is the porosity, S is the phase saturation, andq αm is the source term. We can also rewrite Equation (3) as (4) where B α = ρ αsc /ρ α is the formation volume factor (FVF), ρ αsc is the density at standard conditions, q αsc =q αm /ρ αsc , and α = ρ α g. We assume that the medium is saturated (S β + S γ = 1), and pc = p β − p γ where pc is the capillary pressure. We usually write these equations in terms of nonwetting phase pressure and wetting phase saturation [16] and Time derivative terms in Eqs. (5) and (6) can be rewritten as and ∂ ∂t where φ = dφ dp γ , and we consider that

Initial and boundary conditions
The mathematical model represented by Eqs. (5) and (6) has no unique solution unless appropriate boundary and initial conditions are specified.
The initial condition is given for an arbitrary time instant t 0 and, for example, it could be [16] where p 0 and S 0 represent the pressure and saturation values in the entire domain Ω at the initial time.
For numerical reservoir simulations, we can prescribe initial pressures at a given reference depth, and we use hydrostatic gradients and capillary effects to determine the respective initial values at the other depths [32].
We know that boundary conditions can be of three types: Dirichlet, Neumann or Robin. In the Dirichlet condition, we specify pressure or saturation at the boundary ∂Ω [16] p = p spe (x, y, z, t), S = S spe (x, y, z, t).
On the other hand, in the Neumann condition, when we impose the flow across the boundary [16] ρv · n = f (x, y, z, t) (11) where n indicates the normal outward unit vector. For an impermeable boundary f (x, y, z, t) = 0. Lastly, the Robin boundary condition is a combination of the first two [16] f p p + f v ρv · n = f (x, y, z, t) where f p , f v , and f are known functions.

Properties of fluid and rock
Important rock properties include porosity and permeability. Among the fluid properties of interest, we can mention viscosity and formation volume factor. It is also necessary to consider the variables arising from fluid-rock and fluid-fluid interaction, such as capillary pressure and relative permeability.

Fluid isothermal compressibility
In this work, water compressibility c w (psi −1 ) is estimated by [58] c w = 1 7.033p w − 537T + 541.5s a + 403.3 (13) where T is the temperature ( o F), p w is the water pressure (psi), and s a is the water salinity (g/L). Oil compressibility c o (psi −1 ) is calculated by [65] c o = 1.705 × 10 −7 R 0.69357 where R sob is the solubility rate at the bubble point (scf/STB), γ g = ρ gsc /ρ asc is the gas phase density in solution, ρ gsc and ρ asc are the density at standard conditions of the gas and the air, respectively, p o the oil pressure and is the API gravity, and γ o = ρ osc /ρ wsc is the oil phase density, with ρ osc and ρ wsc standing for oil and water densities at standard conditions.

Formation volume factor
For slightly compressible fluids, we use where the subscript b indicates properties at the bubble point. Still, for gas, we use the real gas law [7], such that where Z is the compressibility factor and we consider that Z sc ≈ 1. This factor can be determined using [30] Z = 1+c 1 (T r )ρ r +c 2 (T r )ρ 2 r −c 3 (T r )ρ 5 r +c 4 (T r , ρ r ), (18)  where the subscript r represents the reduced values of density and temperature and where A i are constants [30] to be provided, and we use a Newton-Raphson method to calculate Z in Eq. (18).

Viscosity
For undersaturated case (p o > p ob ), oil viscosity µ o (cp) can be determined by [66] where c µo is the coefficient of variation of viscosity. For the viscosity of the water µ w (cp), we can consider the following correlation [52]: where µ wb is the water viscosity at 1 atm and reservoir temperature, and c µw is the coefficient of variation of viscosity.
In this work, we calculate the viscosity of the natural gas µ g (cp) using the correlation of Lee et al. [48]: where K = (9.379 + 0, 01607M ) 1.5 209.2 + 19.26M + T , and for ρ g in g/cm 3 , T in o R and M in lbm/mol.

Porosity
We can model the variation of porosity as a function of pressure from a first-order approximation expression [16] where φ 0 is the porosity at reference pressure and c φ the rock compressibility.
2.2.5. Relative permeabilities and capillary pressure Relative permeabilities are determined using a Corey modified model [33] for the wetting phase and for the non-wetting phase. In these equations, S β and S iβ are the saturation and irreducible saturation of the wetting phase, respectively, S γrβ is the residual saturation of the non-wetting phase and k rβmax = Finally, the capillary pressure is estimated as a function of the saturation of the wetting phase by a power-law relation [33] pc(S β ) = pc max where pc max = pc(S iβ ).

III. FINITE VOLUME METHOD
The numerical solution of partial differential equations involves three main parts. Firstly, we discretize the physical domain into small cells (finite volumes). Then, we need to convert partial differential equations and auxiliary conditions (initial and boundary conditions) into algebraic equations. Finally, we must solve a system of algebraic equations.
The governing equations can be discretized using conventional numerical methods, such as the Finite Difference Method (FDM) [59], the Finite Volume Method (FVM) [60,74] and the Finite Element Method (FEM) [28], among others.
In this work, we choose to use the finite volume method. This method is conservative since we integrate the balance equations over each finite volume and, therefore, local and global conversations are assured regardless of the number of finite volumes [49].

Space and time discretizations
Here, we emphasize that the governing equations of the slightly compressible two-phase flow are integrated over each finite volume to obtain the discretized form of these equations.
For the sake of simplicity, we employ a compact notation introduced by Patankar [60]. We consider a three-dimensional finite volume, as shown in Fig. 1. The node (or center of the cell) P has coordinates (i, j, k) and its "neighbours" have coordinates (i − 1, j, k) = W (West), (i + 1, j, k) = E (East), (i, j − 1, k) = S (South), (i, j + 1, k) = N (North), (i, j, k − 1) = B (Bottom) and (i, j, k + 1) = T (Top). We identify the faces of the finite volumes by lower case letters: Next, we perform a time integration, from t n to t n+1 = t n + ∆t, and over the finite volume (depicted in Fig. 1) of governing Eq. (5) and we finally obtain the discretized form of the transport equation (34), for p = p γ and S = S β , where we determine pressure and saturation at the P center of finite volumes. The coefficients a γF are given by the transmissibilities T γf at control volume faces where F stands for W, S, B, E, N and T , d stands for x, y and z and f stands for w, s, b, e, n and t depending on the node, direction and interface. For example, if we only consider the x-direction and we use a similar procedure for the other directions. The other coefficients are given by where δx, δy and δz are the distances between the nodes in their respective directions x, y and z, and A x , A y and A z are the cross-sectional area of control volume faces, and V P = ∆x∆y∆z. On the other hand, we also integrate Eq. (6): and we obtain, for pc = pc(S β ), p = p γ and S = S β ,

IV. NUMERICAL AND COMPUTATIONAL ASPECTS
As is generally known [7,13,16], non-linear partial differential equations describe two-phase flows in porous media. Thus, the discretization of these equations results in a coupled system of non-linear algebraic equations. This system of non-linear equations must be linearized if we want to use numerical methods appropriate for solving systems of linear equations [43].

Linearization
The non-linear terms, present in the governing equations, include transmissibilities, capillary pressures, coefficients of the transient terms and terms of production/injection. We classify, in general, these nonlinearities as weak or strong [32].
Formation volume factor, viscosity, and porosity are associated with weak nonlinearities for liquids, while relative permeability and capillary pressure are related to strong nonlinearities. Therefore, transmissibilities and productivity indexes incorporate both nonlinearities. For two-phase flows containing a gas phase, strong nonlinearities also include the coefficients that depend on the pressure [32].
Here, we define the transmissibilities as Also, we can use the general form: where G f is a geometric term (does not depend on the properties of the fluid), F n+1 pf depends on the pressure and F n+1 Sf depends on the saturation. To illustrate the Picard-type linearization method [57] and the Fully Implicit method [73], we consider F n+1 pe and F n+1 Se terms. We begin with the Picard iteration, for the time n + 1 and the previous iteration v, and this method is conditionally stable [7]. For the Fully Implicit method we have where the properties were considered to depend only on the values in the adjacent cells (P and E). If we use a first-order upwind method, one of the derivatives will be null. The Fully Implicit method is unconditionally stable [7]. The same procedure can be employed to determine the approximations at the faces w, n, s, t, and b for both methods.
We employ the Picard method to evaluate the nonlinear terms that appear in the transient terms. We should linearize the source terms by a method compatible with the linearization adopted for the nonlinearities present in transmissibilities [32].

Approximation of properties at the finite volume
faces An interpolation technique is adopted since we determine the transmissibilities on cell faces and the pressure and saturation on cell centers. An arithmetic mean, which is a second-order approximation [7], is used to evaluate F pf . For example, also known as the central difference scheme [74]. For saturation-dependent terms, F Sf , central difference scheme leads to the appearance of instability, and upwind methods can be applied instead [7]. For a first-order upwind technique, where v xe is the flow velocity at the direction x at face e. Higher-order upwind or TVD (Total Variation Diminishing) non-oscillatory methods can also be thought of as alternatives [49].
A harmonic mean is used for G f to determine the absolute permeability [7] and, at the face e, where A x is the area perpendicular to the flow and We can apply similar expressions to interpolate the absolute permeability at the other faces.

Two-phase wellbore-reservoir coupling
One of the difficulties associated with wellbore modeling in oil reservoir simulations is that the region in the vicinity of the well has the highest pressure gradient, and it is much smaller than the usual size of a finite volume in a computational mesh [17]. To overcome this difficulty, we employ a well-reservoir coupling technique to relate reservoir grid pressure to well pressure by introducing, for example, source or sink terms and the productivity index. The most widely used model for oil reservoir simulations is the equivalent radius technique [62].
In this work, we assume that the well passes through some cells (finite volumes) ψ well in the vertical direction and that each cell c represents a part of the well (c ∈ ψ well ). To describe the well pressure distribution as a function of time we need to know the cell mean pressure (p γc ), the flow pressure at the face of the formation (p wellc ), and the production/injection rate for each phase: q γscc and q βscc .
We consider here production/injection wells with a specified flow rate or pressure. Thus, in the cells containing the well, there are additional unknowns: 1. The flow pressure (p wellc ) and the flow rate of each of the phases (q γscc and q βscc ) if the total production flow (q scsp ) is specified, 2. The production flow for each of the phases (q γscc and q βscc ) and the total flow if we impose the pressure, 3. The flow pressure (p wellc ) if we specify the injection flow (q βsc sp ), 4. The injection rate (q βscc ) if the flow pressure is imposed on the injection well, and we need an additional equation that relates the unknowns in the cell and the well. For each cell c, c ∈ ψ well , the flow rate of nonwetting and wetting phases can be written, respectively, as [5] and where J α representing the productivity index For vertical wells (parallel to the z-direction), the geometric factor G wellc is given by [63] G wellc = 2π∆z c k xc k yc ln req c r well (65) and r eq c = 0.28 where r eq c is the equivalent radius and r well is the well radius.
The total well flow rate, for each phase, is obtained from the sum of all layers [32], that is, and The total flow for the production wells corresponds to the sum of all layers for both phases [32]: For injection wells, we assume that the mobility of the injected fluid (water) is equal to the mobility of the fluid in the reservoir [32]. Thus: where J i is the productivity index for the injection Once again, we obtain the total flow from the sum of all the layers: For wells passing through multiple layers, we adopt a reference level for the computation of the well pressure and relate that pressure to the pressures of the adjacent cells, considering the gravitational effects and neglecting losses due to friction and inertia [32] where p wellc is the well pressure at D c level, p well ref is the well reference pressure at reference level D ref and¯ αc is the fluid specific weight determined with the mean pressure. Substituting Eq. (73) into Eqs. (62), (63) and (69), we obtain, respectively, and Using now Eq. (73) in Eqs. (70) and (72) we obtain, respectively, and

Conservative approaches to time derivatives
Nonconservative schemes can result in numerical instability for strongly non-linear PDE's (such as those for two-phase flows) and will lead to a violation of the mass balance [32]. Therefore, the difference and we can write for a conservative expansion as We use Eq. (80) in the discretization of transient terms in Eqs. (5) and (6), leading to discretized Eqs. (35) and (45).

Time step estimation
The use of variable time steps (the highest possible) would be a wise strategy to reduce the number of iterations to achieve convergence.
Let it n−1 be the number of iterations required for the convergence with the time step ∆t n−1 , [55] uses to calculate the next time step for the IMPES method. A similar but more flexible criterion is where n decr and n incr are the rates of decrease and increase of the time step, it decr is the previous minimum number of iterations required for the convergence, and it incr is the maximum number of iterations (it incr < it decr ). The time step is limited to a maximum value, ∆t n ≤ ∆t max . If the simulation does not converge, the time step is decreased (for example, halved), and the simulation continues.
4.6. Method of solving linear systems An efficient method to solving a linear system whose coefficient matrix is not symmetric is the Preconditioned Biconjugate Gradient Stabilized Method [75]. Algorithm 1 shows the sequence of operations associated with the implementation of this method [8]. The solution vector x and the residue vector r are iteratively updated. The iteration scheme continues until the highest absolute value of the residue is less than tol = 1.10 −8 .
For the Fully Implicit method, an ILU Factor(1) [69] factorization method allowing fill-in is employed, which is a more effective preconditioner although more computationally expensive than the ILU(0).

V. FULLY IMPLICIT METHOD
The most widely studied and known method for solving systems of nonlinear equations is the Newton's (or Newton-Raphson's) method [11]. We use a fully implicit linearization within Newton's method to solve Eqs. (35) and (45), leading to the following system of equations [43] where δX n+1,v+1 = X n+1,v+1 − X n+1,v , J represents the Jacobian matrix, X = (X 1 , X 2 , . . . , X N ) T is the unknown vector, where X i = (p γi , S βi ) T , and Once the linear system of Eqs. (83) has been solved, we calculated If we interpolate the terms evaluated at the cell faces using only the known values in adjacent cells, for a natural ordering [1], the Jacobian matrix has a block structure where each block is a submatrix [7]: where n = P , and m = {A, N, W, P, E, S, B}. Next, we put Eqs. (35) and (45) in a residual form for each cell (inner cells): and From Eqs. (86) and (87), the system of Eqs. (83) can be written, for each interior cell, in the form The linear system that must be solved iteratively within Newton's method has 2N unknowns and equations for a domain composed of N cells. The number of equations can be changed if we consider the wellreservoir coupling.

Well-reservoir coupling
For specified flow rates, the fully implicit treatment of the source terms requires that we evaluate wellbore reference pressure at the current time and iteration, p n+1,v+1 well ref , and we should consider it as an additional unknown. The well-reservoir coupling equations are linearized using a fully implicit method [32] and are given by where, for production wells, we have and, for injection wells, we write Once the system of linear algebraic equations represented by Eq. (83) is solved, including the terms from the well-coupling, p n+1,v+1 well ref can be calculated by being this formulation the same as that found in [32].
Therefore, we will increase the original system (Eq. (83)) by n well (total number of production and injection wells) equations and unknowns, and we solve it iteratively until we achieve a prescribed stopping criterion for each time step.

VI. IMPLICIT PRESSURE EXPLICIT SATURATION
The Implicit Pressure Explicit Saturation (IMPES) method is one of the most employed in simulating twophase flow in oil reservoirs. As well known, the pressure and saturation equations are solved separately using implicit and explicit schemes, respectively [14]. Capillary pressures and transmissibilities are evaluated explicitly in time n or time n + 1 in a known iteration (n + 1, v).
An equation for the pressure of the non-wetting phase, in each cell, can be obtained by combining Eqs. (35) and (45) where q Isc ≡ B γP q γsc + B βP q βsc .
The main idea behind this formulation is that the total transmissibility remains approximately constant for some time, although the individual transmissibilities of each phase may vary strongly as a function of the saturation of the wetting phase (water). Therefore, the pressure field remains stable, although the saturation of the wetting phase may vary strongly [64].

ISSN: 2349-6495(P) | 2456-1908(O)
After the calculation of the pressure from Eq. (93), this pressure is used to determine the saturation explicitly from the linearized form of Eq. (45) Equations (93) and (102) are solved until we reach convergence.
The IMPES method to be stable requires often prohibitively small time steps due to the explicit calculation of saturation, especially for long simulation times, refined meshes, and certain classes of problems [16]. To circumvent such a situation, we can use different time steps for pressure and saturation calculations since the non-wetting pressure varies more slowly than the saturation of the wetting phase [39]. In fact, the implicit calculation of Eq. (93) demands a computational effort which is higher than that required for the solution of Eq. (102) explicitly.

Well-reservoir coupling
To avoid the numerical instability in the computation of p well ref (for specified flow rate), we adopt an implicit formulation where we couple its solution with the non-wetting pressure, Eq. (93). For the set of cells c ∈ ψ well that contains the production wells, from Eqs. (74) and (75), we have The substitution of Eq. (103) into Eq. (93) introduces an additional unknown, p n+1,v+1 well ref , which can be calculated by incorporating Eq. (76) into the system. For injection wells (q γscc = 0) and from Eq. (77): for production wells, or from for injection wells.

VII. SEQUENTIAL METHOD
As pointed out earlier, the IMPES method is conditionally stable [61,54]. Therefore, the Sequential method [51] aims to overcome this issue and introduces an implicit saturation procedure, but the pressure and saturation equations are solved separately. Thus, we obtain the pressure as in the IMPES method, and subsequently, we calculate the saturation implicitly [71].
Using the chain rule, we can show [64] that we can express the discretized equation (45) for the wetting phase in terms of saturation. Then, an implicit calculation and a Picard linearization result in (for S = S β ) where a SF ≡ a βF dpc dS f , (108) Therefore, Eq. (93) allows for an implicit determination of the pressure. Once the pressure is known, saturation is also calculated implicitly from Eq. (107). Until we attain the convergence criterion, we solve these equations iteratively.

Well-reservoir coupling
For a specified flow rate for the set of cells c ∈ ψ well , we treat q n+1,v Iscc (Eq. (93)) in the same way as in the IMPES method. Therefore, we replace Eq. (103) into Eq. (93) for production wells or Eq. (104) for injection wells, and we thus introduce the unknown p n+1,v+1 well ref . Once we have calculated it (see [27]), we can determine q n+1,v βscc in Eq. (107) from Eq. (105) for production wells or from Eq. (106) for injection wells.

VIII. A PICARD-NEWTON APPROACH
We dedicate this section to discussing an approach to solving the governing equations of the twophase flow, Eqs. (35) and (45). As in the Sequential method, the main idea is to improve numerical stability. However, we now calculate saturation using a fully implicit linearization resulting in an approximation of Newton's method. To obtain the numerical solution, we choose a strategy based on an operator splitting technique [73].
As seen in Section VI., the pressure formulation given by Eq. (93) is stable even though the saturation of the wetting phase may vary strongly. Concerning producing and injecting wells, we know that this formulation remains stable since production decreases the pressure of the reservoir while injection increases it and that q Isc contemplates source terms associated with production and injection wells.
Our approach, named the Hybrid method achieves improvements concerning even the IMPES and Sequential methods. In these methods we linearize, by a Picard iteration, the transmissibilities, the capillary pressures, and the productivity indexes. In this case, these properties can vary strongly with the saturation and render the formulation unstable. On the other hand, we can minimize this problem by choosing a fully implicit formulation for linearization.
We begin by calculating pressure at time n + 1 and iteration v + 1 from Eq. (93) and we consider that it is no longer unknown for the determination of the saturation. The fully implicit linearization of Eq. (45) leads to a system of equations, which is an approximation of Newton's method. The term "approximation" is used since the pressure is not known, but only an estimate of it at (n + 1, v + 1). Therefore, we can express this system of equations in the following matrix form: where δS n+1,v+1 ..,S w N ) T is the unknown vector, and R w = (R w1 , R w2 , R w3 , ...,R w N ) T is the residual vector with n = 1, 2, 3, ...N for a reservoir containing N cells.
After we have solved the system of linear algebraic equations (110), saturation is obtained by for each cell of the computational domain. We can write the system (110) for each inner cell as follows: where Accordingly, Eq. (93) is solved implicitly for the pressure followed by a further implicit solution of Eq. (110) for the saturation. We solve these equations iteratively until we reach a prescribed stopping criterion.

Well-reservoir coupling
As with IMPES and Sequential methods for a specified flow rate, we determine p well ref by employing an implicit formulation, and we calculate the well pressure in conjunction with the reservoir pressure obtained from Eq. (93).
Once known p n+1,v+1 well ref , we can determine q n+1,v βsc and their derivatives, which appear in Eq. (110), from In this formulation, we adopt the fully implicit linearization for both saturation of the wetting phase and its source term.

IX. NUMERICAL RESULTS
This section presents a comparative study of the computational efficiency of the different approaches used here to solve the discretized equations that govern a slightly compressible (water-oil) and compressible (water-gas) isothermal two-phase flow of immiscible fluids in a reservoir. Therefore, to be able to verify the accuracy of the Picard-Newton method, we perform a comparative study considering the following methods: 1. Fully Implicit (Newton's method); 2. IMPES; 3. Sequential; 4. Hybrid.
As applications, we address the numerical simulation of: 1. water-oil and water-gas flow in a slab-type geometry; 2. water-oil flow in a five-spot pattern; 3. water-gas flow in a reservoir connected to an aquifer.
Both homogeneous and heterogeneous cases are considered here for the first two cases.

Slab geometry
We now turn to the case of three-dimensional flow in a slab reservoir, with dimensions L x , L y , and L z . We specify a water flow rate and a saturation value at x=0 and a prescribed pressure at x = L x . For the other boundaries, we have no flux conditions. Therefore, the main flow is along the x-axis from left to right (Fig. 2). 9.1.1. Case 1: water-oil flow in a homogeneous medium It is injected water into a homogeneous reservoir containing oil and an initial amount of water (S w0 = 0.2) through the boundary at x=0. The pressure at the top of the reservoir (level of reference) is p=5,800 psi, and determined from hydrostatic in the other cells. The temperature of the reservoir is kept constant at T =150 o F. Tables 1-3 show oil, water, and gas properties. Table 4 contains the properties of the rock and other essential properties for two-phase flow: absolute permeabilities, initial porosity, the compressibility of rock, relative permeabilities, among others. We employ these properties in all simulations performed here unless otherwise stated. We impose a uniform pressure gradient, ∂p/∂x = -0.5 psi/ft in the yz-plane, on the left border of the reservoir, and we stipulate a total simulation of t ∼ = 5,000 days (due to the use of a variable time step).   We performed numerical simulations considering a three-dimensional reservoir with dimensions L x =2,000 ft, L y =1,000 ft, and L z = 100 ft and meshes with n x , n y , and n z finite volumes in the three spatial coordinates directions. For the mesh refinement experiment, we considered three computational meshes with a different number of finite volumes in the xdirection (the direction of the flow), Table 5. Figure 3 shows the water saturation profiles along the centreline in the yz-plane (y = 500 ft and z = 50 ft) for t=1,500 days, t=3,000 days, and t=4,500 days, for the three meshes, a constant time step ∆t = 0.1 days, and the Hybrid method. We also see the corresponding oil pressure profiles in Fig. 4. We observe that results do converge when we refine the meshes. We obtain the same results with the Fully Implicit, IMPES, and Sequential methods, and we do not represent them for the sake of brevity. Table 4: Fluid-rock properties.

Parameter
Value Unit 4.0 psi (water-gas) pc max 5.8 psi (water-oil) Mesh n x n y n z 1 125 50 5  2  250 50 5  3  500 50 5 As expected, we see a rarefaction followed by a smoothed discontinuity curve for saturation values (Fig. 3) as a result of compressibility and capillary pressure effects [47]. Although the first-order upwind introduces significant numerical diffusion for unrefined meshes, it does not introduce spurious oscillations [49]. For a compressible flow, the saturation front is delayed concerning the saturation front for an incompressible flow (Buckley-Leverett's problem) since part of the fluid energy is used to compress it. The pressure profiles (Fig. 4) also behave as forecast [47] changing its profile as the front moves forward.
For the accuracy considered, we did not find any significant differences in the values obtained with the other methods: Fully Implicit, IMPES, and Sequential. It is worth mentioning that as we use variable time steps, we did not expect to obtain the same solution since higher time steps lead to higher numerical truncation errors, mainly for a first-order scheme [59]. Figure 5 shows the water saturation field obtained with the Hybrid method and Mesh 3. As can be seen, the front moves maintaining the same uniform and flat profile (without numerical instabilities) along the yzplane, according to the expected physical behavior since the medium is homogeneous, and the injection rate is uniform.  9.1.2. Case 2: water-oil flow in a heterogeneous medium In the following test problem, water is injected in a reservoir maintaining the same physical conditions of Case 1, but now the reservoir is heterogeneous with a permeability field given by (adapted from [26]) where with G(x, y) = sin 6πx L x cos 12πy L y (117) and the max function provides the maximum value between two real numbers. Figure 6 shows an illustration of the permeability field that we have created (Eq. (115)). For this second problem test, we use the same three-dimensional slab geometry. About the domain dimensions, and the number of finite volumes, the values can be found in Table 6. In this case, we construct the computational mesh such as to obtain the spatial increments ∆x and ∆y equal to the space increment ∆x of Mesh 2 utilized in the previous simulation. The maximum simulation time is t ∼ = 6,000 days. Dimensions (ft) Mesh L x = 1,500 n x = 188 L y = 1,000 n y = 125 L z = 100 n z = 5 As in the previous case, we did not find significant differences between the numerical solutions determined using the different methods discussed in this work. For the simulations using the Hybrid method, we can see the displacement of the water saturation field in Fig. 7 for different times. The saturation front moves along the yz-plane with a characteristic profile as a consequence of the particular heterogeneous permeability field. Due to the difference in the permeability values, the flow initially bypasses the regions of low permeability (see Fig. 6). However, as time goes by, these regions will also be filled with the injected fluid.

Buckley-Leverett's problem
To verify the efficiency of our numerical code, we have solved the Buckley-Leverett problem [6]. For this specific problem, we employ the same reservoir dimensions and physical properties of Case 1. However, as well know, Buckley-Leverett's problem considers a incompressible two-phase flow without taking into account the effects of capillarity.
The water saturation profile, along the centerline in the yz-plane, can be viewed in Fig. 8 for Mesh 2 (see Table 5) and t=800 days, t=1,600 days, and t=2,400 days. We obtained the results for the water saturation profile using only the Hybrid, IMPES, and Sequential methods keeping constant the time step: ∆t=0.1 days. We did not apply the Fully Implicit method because according to Kardale [43], difficulties may arise, concerning convergence, if we use this specific method to solve Buckley-Leverett's problem.
When compared to analytical results, from the results shown in Fig. 8 we can verify that the Hybrid, IMPES, and Sequential methods present accurate results. As we can check, we obtained the same saturation profile with the three methods. Considering now the accuracy of these results, when compared to the analytical solution, we can notice the appearance of numerical diffusion [49], although the results are very close to the theoretical ones. We can explain the emergence of numerical diffusion by the use of a first-order upwind and Euler schemes and the size of the finite volumes and time steps. Nevertheless, there are no spurious oscillations and numerical dispersion, and we reproduce the shock (a sharp discontinuity) numerically by only two points. Therefore, we think that our results are in good agreement with the analytical solution [49].

Case 3: water-gas flow in a homogeneous medium
We now turn to the problem of injecting water into a homogeneous three-dimensional reservoir containing gas. As initial conditions, we take for the water saturation S w =0.2 and P ref equal to 4,000 psi at the top of the reservoir. The reservoir temperature is constant (T ref = 150 o F) and the parameters used for calculating the physical properties are shown in Tables 3 and 4. A uniform pressure gradient (∂p/∂x = -0.5 psi/ft) was set in the yz-plane at the left border and t ≈ 7,000 days is the maximum simulation time.
A mesh refinement study was performed for the meshes presented in Table 5 and considering that L x =2,000 ft, L y =1,000 ft and L z = 100 ft.
For the Hybrid method, we used an initial time step equal to 0.1 days and Fig. 9 shows the water saturation profile along the centreline in the yz-plane while Fig. 10 shows the corresponding oil pressure profiles, both for t = 1,500 days, t = 3,000 days and t = 4,500 days. The results are very close to those obtained with Mesh 3, although slight differences can be observed for Mesh 1. Differently from water-oil results, Fig. 3 and Fig. 4, they are now closer to each other. Similar profiles were obtained with the Fully Implicit, IMPES and Sequential methods. However, for the sake of brevity, they are not presented.
Numerical results behave as physically expected ( Fig. 9) and instead of a shock, which would appear for an incompressible flow, we have a jump in saturation after the rarefaction due to compressibility and capillary effects and also to numerical diffusion coming from the first order schemes. By comparing the profiles in Figs. 3 and 9, the effects of the compressibility on forwarding front velocity can be verified. As gas has compressibility higher than that of oil, there is a delay for water saturation profiles, for water-gas flow, when compared to those of water-oil flow for the same flow conditions. We also observe that water saturation values at the end of rarefaction are higher than those for water-oil flow. Regarding pressure (Fig. 10), its variation is following the expected behavior with a sudden change at the end of rarefaction. Figure 11 presents water saturation values for different instants of time, Mesh 3, and the Hybrid method. Although not shown, we obtained similar results for the other methods studied here. The flow occurs in a homogeneous medium at a constant in-jection rate so that it consists of a uniform water front moving with a constant velocity in the yz-plane, and no numerical instability appears until the end of the simulation.   9.1.5. Case 4: water-gas flow in a heterogeneous medium The next case, still adopting the slab geometry, addresses a water-gas flow in a heterogeneous reservoir, under the same conditions as the previous problem, but having now a permeability field described by the function where F (x, y) is given by Eq. (116). We can see the distribution of the permeability values of the generated heterogeneous field in Fig. 12. We can find the dimensions of the reservoir and the mesh used in Table 6. The simulation was performed considering a final time t ∼ = 8,000 days. Figure 13 shows, for six different instants of time, water flowing through the heterogeneous porous medium. We calculated the saturation field with the Hybrid method, and we did not find significant differences between these values and those obtained with the other methods. So, we conclude that the generated flow is realistic and compatible with the permeability field employed in the simulation (Fig. 12).

Five-spot injection and production arrangement
We widely use injection in a five-spot pattern to verify two-phase flow numerical simulations. It consists of four injection wells and a single production well. For a square domain, injection wells are positioned at the vertices with production well at the center. Due to symmetry of the problem only a quarter of the five-spot pattern is considered in the simulations, see Fig. 14. 9.2.1. Case 5: water-oil flow in a homogeneous medium In the following, we focus on the water injection process in a homogeneous reservoir containing oil and an initial amount of water. We use here the same initial conditions and properties of Case 1. Table 7 display the geometric data and computational meshes that we use in this case. Injection and production wells are positioned perpendicular to xy-plan at (x, y)=(0,0) ft and (x, y)=(1,000,1,000) ft, respectively, and both occupy the full extent of the z-axis.
We impose a production flow rate of -200 STB/D and an injection flow rate of 200 STB/D. For the variable time step we use: n decr = 0.75; n incr = 1.1; it decr = 25; it incr = 15; ∆t max = 10 days and ∆t 0 = 0.01 days. The simulations end when time reaches approximately 7,000 days. Dimensions (ft) Mesh L x = 1,000 n x = 125 L y = 1,000 n y = 125 L z = 100 n z = 5 We can see the variation of the water saturation over time in Fig. 15. Water moves according to Darcy's law and shows the characteristic profile well known in the literature [53], from the injection well towards the production well, and may be more or less pronounced depending on the properties of the reservoir and fluids [33]. Even though we have determined these values with the Hybrid method, we have not found differences between these saturation values and those calculated with the other methods.  9.2.2. Case 6: water-oil flow in a heterogeneous medium Our attention turns now to a five-spot pattern with the same two-phase flow, parameters, and physical properties of the previous problem, but now the reservoir is no longer homogeneous. There are two regions within the domain with different permeability values (one hundred times smaller) from that of the rest of the reservoir, Fig. 16 [26]. By setting x 0 , y 0 , and z 0 the origins of these regions and l x , l y , and l z the lengths along the coordinate axes of these regions, Table 8 shows their respective values. Region (x 0 , y 0 , z 0 ) (l x , l y , l z ) 1 (100,100,0) (400,400,100) 2 (600,600,0) (300,300,100) Fig. 16: The quarter five-spot configuration of injection and production wells with two block inclusions. Figure 17 contains the graphical representation of the water infiltration process for a five-spot pattern with the reservoir containing the two regions with different permeability values. We calculated the saturation values using the Hybrid method. As well known, it is natural that water flows initially through the regions that offer less resistance to flow, that is, those with the highest values of absolute permeability. Thus, the main flow tends to deviate from the two regions. We remind the reader that although there is a deviation of the flow, bypassing the regions of low permeability, these regions are not impervious to fluid flow. Therefore, both phases will be present in all reservoir regions, and over time the main flow will be directed to the producing well. Due to the nature of the flow, water saturation values change very slowly, and we can observe at the end of simulations that both regions contain water with a saturation very close to its initial value (0.2).

Production of gas by combined effects
In the last example, we study the production of gas for a reservoir in contact with an aquifer. Therefore, by a combination of effects we can produce gas, that is, gas and rock decompression and the pressure exerted by the water contained in the aquifer (constant pressure condition at the lower boundary), Fig. 18.
The production of gas leads to a pressure drop in the reservoir. Nevertheless, the expansion of the water of the aquifer can counterbalance it. As a consequence, there is an inflow of water into the reservoir. This fact contributes to the preservation of the reservoir pressure [33]. 9.3.1. Case 7: water-gas flow in a homogeneous medium As we had mentioned earlier, we are dealing with gas production in a homogeneous reservoir (containing gas and water initially) with a water inflow due to the existence of an aquifer. The properties used, as well as the initial conditions, are the same as those described for Case 3. We can find the computational meshes employed and the domain dimensions in Table 9. A vertical production well is positioned perpendicular to xy-plan located at (x, y) = (500,500) ft, with 75 ft length, and we produce gas at a constant rate: -200,000 SCF/D. We also use a variable time step: n decr = 0.75; n incr = 1.1; it decr = 20; it incr = 10; ∆t max = 15 days and ∆t 0 = 0.1 days and we specify a maximum simulation time of approximately 5,000 days. Here, we apply the same well-reservoir techniques introduced in Subsection 4.3.  Dimensions (ft) Mesh L x = 1,000 n x = 125 L y = 1,000 n y = 125 L z = 150 n z = 30 We can visualize the water saturation values for simulations performed with the Hybrid method, for six different instants of time, in Fig. 19. From the observation of these figures, we can see the water inflow from the aquifer towards the porous medium as we recover the gas. Since the reservoir and aquifer are connected, water inflow towards the reservoir tends to limit reservoir pressure drops during gas production.

X. CONCLUSION
This work aimed to perform a comparative study of some of the main techniques that we usually use in the numerical solution of two-phase flows in the context of reservoir simulation. With this purpose, we proposed a Picard-Newton approach, referred to as the Hybrid method. We have accomplished this task using a three-dimensional numerical simulator, programmed in C language, based on the finite volume method.
The motivation for proposing a Picard-Newton approach for the calculation of the pressure and saturation has assumed that pressure is typically a global variable, whereas saturation is a local variable. Therefore, when capillary effects are small, the decoupling of the system of governing equations leads to better convergence properties. We obtain the equation for the pressure from the combination of continuity and Darcy's law equations. The goal is to introduce the sum of the transmissibilities of the wetting and non-wetting phases, eliminating the saturation in this equation. Although the individual transmissibilities of each phase vary strongly with the saturation of the wetting phase, their sum remains approximately unchanged over a time interval. In the calculation of the saturation of the wetting phase, a Fully implicit linearization is adopted to improve the numerical stability of the Hybrid method concerning the IMPES method, which explicitly solves saturation, and the Sequential method, which implicitly solves saturation using a Picard linearization.
To carry out this comparative study we have considered: water-oil and water-gas flows in a slab geometry for a homogeneous and heterogeneous medium; water injection in a five-spot pattern for a homogeneous medium and a homogeneous medium containing two regions with low permeability values; and gas production in a homogeneous reservoir connected to an aquifer.
We have considered the Buckley-Leverett problem, whose analytical solution is known, to verify our numerical code. We have obtained numerical values very close to each other with the IMPES, Sequential, and Hybrid methods. Besides, the results have presented a good agreement with the analytical solution. Although we have perceived some numerical diffusion, introduced by the use of the first-order upwind and forward time schemes, it has not compromised the sharpness of the profile, and we have satisfactorily captured the shock after the rarefaction.
All results obtained in this work are physically correct, and in none of the analysed cases, we have not found significant differences concerning the precision of the results for the Fully Implicit, IMPES, Sequential, and Hybrid methods. It is important to note that since we have used variable and distinct time steps, we have not expected to obtain the same solution for all methods, as a consequence of the use of higher time steps. As we have employed first-order time schemes, this also leads to higher truncation errors.