Improved IMPES Scheme for the Simulation of Incompressible Three-Phase Flows in Subsurface Porous Media

: In this work, an improved IMplicit Pressure and Explicit Saturation (IMPES) scheme is proposed to solve the coupled partial differential equations to simulate the three-phase ﬂows in subsurface porous media. This scheme is the ﬁrst IMPES algorithm for the three-phase ﬂow problem that is locally mass conservative for all phases. The key technique of this novel scheme relies on a new formulation of the discrete pressure equation. Different from the conventional scheme, the discrete pressure equation in this work is obtained by adding together the discrete conservation equations of all phases, thus ensuring the consistency of the pressure equation with the three saturation equations at the discrete level. This consistency is important, but unfortunately it is not satisﬁed in the conventional IMPES schemes. In this paper, we address and ﬁx an undesired and well-known consequence of this inconsistency in the conventional IMPES in that the computed saturations are conservative only for two phases in three-phase ﬂows, but not for all three phases. Compared with the standard IMPES scheme, the improved IMPES scheme has the following advantages: ﬁrstly, the mass conservation of all the phases is preserved both locally and globally; secondly, it is unbiased toward all phases, i.e., no reference phases need to be chosen; thirdly, the upwind scheme is applied to the saturation of all phases instead of only the referenced phases; fourthly, numerical stability is greatly improved because of phase-wise conservation and unbiased treatment. Numerical experiments are also carried out to demonstrate the strength of the improved IMPES scheme.


Introduction
Numerical simulation for subsurface flows has been extensively applied in industry, such as in the management of groundwater energy and waste pollutants, petroleum engineering, and exploitation of geothermal energy [1][2][3][4][5][6][7][8][9][10][11][12][13]. The main numerical methods in the simulation include the Fully Implicit scheme (FI) [14][15][16][17] and the IMplicit EXplicit scheme (IMEX) [18][19][20]. In the FI scheme, all the unknowns can be derived implicitly. Therefore, it is unconditionally stable. However, it is extremely time-consuming to deal with the nonlinear equations at each time step. Different from the FI scheme, the IMEX scheme processes the linear terms implicitly, while the nonlinear terms explicitly. As a result of this, only conditional stability is guaranteed. More specifically, the IMplicit-Pressure and Explicit-Saturation (IMPES) scheme [15,19,[21][22][23][24][25][26] is regarded as a kind of IMEX scheme, which is extensively leveraged in the three-phase flow simulation in porous media. In the IMPES scheme, the pressure is derived implicitly, then the Darcy velocity and saturation are explicitly derived. The IMPES scheme provides the solution of high accuracy, with less time cost than the FI scheme. Moreover, the IMPES scheme is easy to implement and has, therefore, gained great popularity.
In this work, a modified IMPES scheme is designed to simulate incompressible threephase flows in porous media with the discrete methods of cell-centered finite difference and upwind [27]. An advanced discretized formulation is implied to obtain the total pressure by summing up the discretized conservation equation of each phase, which could guarantee the consistency of the pressure equation. Furthermore, the mass conservation of all the phases can be preserved both locally and globally in this scheme.
This work is organized as follows. Section 2 presents an incompressible and immiscible three-phase flow model. Section 3 reviews the standard IMPES scheme. Section 4 describes the presentation of an improved IMPES scheme. Section 5 shows some applications of our new method.

Mathematical Model
In this work, we consider the flow of three incompressible and immiscible phases in porous media. This three-phase flow system can be used for the simulations in a number of natural and engineering circumstances, in particular, the petroleum reservoirs consisting of water, oil and gas phases if the density of the gas phase does not substantially change in the system. We note that the gas phase needs to be modeled by a compressible phase by using, for example, an equation of state, if its density varies substantially within the domain. The assumptions of incompressibility and immiscibility are for the convenience of presentation and numerical implementation here, but most of our results can be extended to compressible three-phase flows in a straightforward way. We use the subscripts w, o and g to denote the three phases (naming from water, oil and gas phases). The equations of immiscible and incompressible three-phase flow are described by Here, φ is the porosity, k is the absolute permeability, S β , ρ β , p β , k rβ , q β , u β and µ β are the saturation, the density, the pressure, the relative permeability, the injection or production rate, the Darcy velocity and the viscosity of β phase, respectively, while p cij is the capillary pressure, g is the gravity and z is the depth.
The domain Ω ⊂ R 2 is bounded. Γ = ∂Ω is the boundary condition of Ω, and ∂Ω = Γ D Γ N , Γ D Γ N = ∅, where Γ D and Γ N indicate Dirichlet and Neumann boundary conditions. The inlet boundary is presented by Γ in ⊂ Γ, where Γ in = {x ∈ Ω : u t · n < 0}, u t = u w + u o + u g is the total Darcy velocity. The following initial and boundary conditions are added to close the system where n is the unit normal vector.

Standard IMPES Scheme
Following the last section, we can use standard IMPES scheme [24] to solve the PDEs. The key point of the scheme is to decouple pressure and saturation, and the first step is to sum up conservation equation of each phase, leading to where q t = q w + q o + q g is the total injection or production rate and u t = u w + u o + u g is the total Darcy velocity.
In the above relation, we find that if the saturations S m w and S m g are known, the pressure Equation (16) is linear for ψ o . By relating the Darcy velocities of the water and gas phases to the total Darcy velocity, we find By substituting Equations (17) and (18) into (1), the water and gas phases saturation equations can be read as where If the saturations S m g and S m w at time step m are given, we can implicitly solve ψ m+1 o by the pressure Equation (16), and update the saturations S m+1 w and S m+1 g by the saturation Equations (19) and (20). Consequently, the oil phase saturation S m+1 o can be calculated by 1 − S m+1 w − S m+1 g . The saturation-pressure system can be decoupled by the standard IMPES scheme which is summarized in Algorithm 1 for completeness. Algorithm 1 Standard IMPES Scheme 1: Step 1. Given S m w and S m g , calculate u m+1 cow and u m+1 cgo by using Step 2. Given S m w , S m g , u m+1 cow and u m+1 cgo , find u m+1 t,o and ψ m+1 o such that Step 3. Given S m w , S m g , u m+1 t,o and ψ m+1 o , update saturation of each phase by The first and second equations in Step 3 of Algorithm 1 are clearly consistent with Equation (1) when β = w, g. However, the third equation may not satisfy Equation (1). Therefore, the standard IMPES scheme is locally mass conservative for the two reference phases provided that the spatial discretization is also locally conservative, but this local conservation is not preserved for the third phase. In next section, we will propose a fully mass conservative IMPES scheme for this incompressible and immiscible three-phase flow system.

Improved IMPES Scheme
In this section, we will introduce a modified scheme to improve the consistency of the standard IMPES method. The governing equations for the incompressible and immiscible three-phase model have been expressed in Section 2. The detailed discretization procedure is showed below.
We first define w β := −k(∇p β + ρ β g) and u β := λ β w β . We use a rectangular mesh for spatial discretization, and the domain Ω = (0, Lx) × (0, Ly) can be divided into M × N cells, partitioned by δ x × δ y , where We introduce the following standard notation: , , , The cell-centered finite difference method is employed here to discretize pressure gradient ∇p(x, y) at cell centers (x i− 1 2 , y i− 1 2 ), and more details can be found in Rui et al. [28], Weiser and Wheeler [29] and Negaral et al. [21]. The gradient of a general scalar quantity p can be approximated by using a central finite difference in each direction: We can estimate w β at the middle of edges by Here, we only consider the simplified situations, but Starnoni et al. [30] has applied to one-and two-phase flow in porous media for the general circumstances. Subsequently, we will evaluate Darcy velocity by the upwind method, which needs velocity direction information from the last time step. In the very first time step, the direction velocities information is not yet known, and we propose to use of the following formulations: For all later time steps (i.e., m ≥ 1), the Darcy velocities are computed by The saturation S m, * β for all phases can be discretized by the traditional upwind method With the above pressure and velocity formulations, we can obtain the discretized forms of Equations (1)-(5) at each time step m, where, By summing up the above discrete conservation equations of all the three phases, the total conservation equation is presented as The above equation has three unknown pressures p m+1 We need the following two additional equations to implicitly solve p m+1 We consider p m+1 o as the primary unknown for pressures, and we can eliminate the other two pressures p m+1 g and p m+1 w from the system. To do this, we substitute Equations (35) and (36) Finally, the saturation S m β is explicitly updated by using the mass conservation law for each phase. The complete procedure is shown in Algorithm 2. It is clear that the fully conservative IMPES is unbiased toward any of the three phases. We state two methods to update saturation, which is shown in step 3 of Algorithm 2, and the equivalence between these two methods can be proven.   (29).
Proof. Firstly, using Method II, we can obtain the following three equations: With the relation (34), we can obtain Secondly, if Method I is used to this improved IMPES scheme, we have the forms below φ S m+1 Combining Equation (34) with the restriction S k w + S k g + S k o = 1, k = m, m + 1, we can quickly find the following desired result Thus, Method I and Method II are equivalent.
The estimated saturation of three phases meets the discretized mass-conservation Equation (29). For T is a finite partition of Ω, and any cell E ⊂ T , We have the following local conservation of mass for three phases.
Using the divergence theorem, we have, where, n is the outer unit normal vector. When , applying the explicit method for time discretization, the Equation (49) can be rewritten into Since the rectangular cell E i− 1 by the midpoints of cell edges. We also note that the Equation (50) can be consistent with Equation (29)  . Therefore, the improved IMPES scheme is local mass-conservative for three phases.

Numerical Experiments
In this section, four numerical experiments are simulated using our improved IMPES scheme formulated in the previous section. These numerical results demonstrate various features of the proposed scheme. The gravity and capillary pressure are neglected in following three examples for convenience, but the proposed model and numerical scheme in this paper are also applicable to cases with non-zero gravity and capillary pressure. The gravity and capillary pressure effects are considered in example 4, which leads to the counter current flow. The numerical experiment of example 4 is used to verify that the proposed scheme can work well for the situation counter current flow.
Example 1. This example shows a three-phase flow in a homogeneous porous media with a size of 100 × 100 m and a permeability of 100 md. The media are initially saturated with water, oil and gas which have a volume fraction of 10%, 70% and 20% respectively. Water is injected into the porous media at northeast corner of the domain with a constant rate of 6.98×10 −6 m 3 /s. Water-oil-gas mixture is produced at the southwest corner of the domain, where the pressure is fixed at 1 atm. Figure 1 shows the pressure at different time.    Figure 5. Water is injected at the west edge of the domain with a constant rate, and the mixture is produced at the east edge of the domain where pressure is fixed at 1 atm. Figure 6 shows the pressure at different times. Figures 7-9 illustrate the saturation of water, oil and gas at different times. Example 3. This example shows a three-phase flow in a homogeneous porous media with a square domain of [0, 120 m] × [0, 120 m] and a permeability is 120 md. In the simulation of beginning, a block water is located at the center of the domain, there is no oil in the rest of the domain. Water is injected into the porous media at the northeast corner of the domain with a constant rate of 6.98×10 −6 m 3 /s. Water-oil-gas mixture is produced at the southwest corner of the domain, where the pressure is fixed at 1 atm. Figure 10 shows the pressure field at different time. Figures 11-13 illustrate the saturation of water, oil and gas at different time.   Figure 14. The gravity and capillarity is considered in this example. Water is injected in the central part of west boundary at a constant rate of 1.3952 × 10 −4 m −4 /s and the pressure is set as 1 atm on the upper and lower parts of the east boundary where the oil-gas-water mixture is yielded. No flow conditions are enforced on the rest parts of the boundary. The initial distributions of oil, gas, water is shown in Figure 15. Figure  16 shows the Darcy velocity fields of the three phases. Figure 17 illustrates the saturation distributions of water, oil and gas phase after 91,251 days. From the simulation results, we find that the proposed numerical scheme is able to simulate this counter-current flow situation well. Finally, we calculated the residual (mass conservation error) of saturation equation for each phase at every time step in the simulation by l ∞ norm given by formula (51). For the improved IMPES scheme, the calculated mass conservation errors of three phases at all simulation time steps are in the magnitude of 10 −16 , which is considered to be zero without the round-off error. We also ran this example with the same conditions using a conventional IMPES scheme and calculated the residual of mass conservation equation of each phase. The calculated errors of two referred phases (water and gas) at all simulation time steps are in the magnitude of 10 −16 . However, we also calculated the global error for each phase at all simulation time steps and found that the global error of non-referred phase is accumulating from in magnitude of 10 −10 to 10 −4 .

Conclusions
In this work, we have proposed an improved IMPES to simulate three-phase flow systems. Conventional IMPES schemes discretize the pressure equation and two (out of three) saturation equations independently. While this independence could allow one to choose a separate desired numerical algorithm for each of the pressure and saturation equations, the discretized pressure equation is unfortunately inconsistent with the discretized saturation equations. To address this issue, we proposed an improved IMPES scheme, where the discretized pressure equation is fully consistent with the discretized saturation equations. Specifically, the discrete pressure equation in our algorithm is not obtained by discretizing the continuum pressure equation as most people did in the literature; instead, it is obtained by summing up the discretized formulations of the conservation equations for all phases.
Even though the decoupling of the pressure and saturation is majorly based on conventional IMPES schemes, the consistency in our decoupling yields a number of desired numerical features. The most significant improvement in the numerical behavior is the mass conservative property honored in each of the three phases both locally and globally. Conventional IMPES usually solves only two saturation equations numerically (which we refer to as two reference phases), but not the third saturation equation; thus it has a built-in bias, which is eliminated in our improved IMPES scheme. In particular, the conventional IMPES upwinds only the two reference phases, but our scheme is able to apply the upwind scheme to all three phases.
We used the common cell-centered finite difference method for spatial discretization in our numerical examples, but other locally conservative spatial discretizations can also be used in our improved IMPES. Four numerical cases have been carried out to demonstrate the strengths of the improved IMPES scheme. Various flow patterns and a number of different boundary conditions were numerically studied. We also repeated the same numerical examples using conventional IMPES (not shown in this paper for brevity) and enhanced numerical stability has been observed in our scheme as compared to the conventional one. In numerical example 4, we also considered the effects of capillary and gravity for three-phase flow in heterogeneous porous media. Since the gravity segregation among three phases, the counter-current flow is appeared in vertical orientation of the reservoir. We did not consider the compressibility and miscibility of gas for three-phase flow in this study, which is to consider compressible and partially miscible three phases as our near future work.