Iterative Numerical Scheme for Non-Isothermal Two-Phase Flow in Heterogeneous Porous Media

In the current paper, an iterative algorithm is developed to simulate the problem of two-phase flow with heat transfer in porous media. The convective body force caused by heat transfer is described by Boussinesq approximation throughout with the governing equations, namely, pressure, saturation, and energy. The two coupled equations of pressure and saturation are solved using the implicit pressure-explicit saturation (IMPES) scheme, while the energy equation is treated implicitly, and the scheme is called iterative implicit pressure, explicit saturation, implicit temperature (I-IMPES-IMT). In order to calculate the pressure implicitly, the equations of pressure and saturation are coupled by linearizing the capillary pressure which is a function of saturation. After that, the equation of saturation is solved explicitly. Then, the velocity is computed which is used in the energy equation to calculate the temperature implicitly. The cell-centered finite difference (CCFD) method is utilized for spatial discretization. Furthermore, a relaxation factor along is used with the Courant–Friedrichs–Lewy (CFL) condition. Finally, in order to illustrate the efficiency of the developed algorithm, error estimates for saturation and temperature for different values of time steps and number of iterations are presented. Moreover, numerical examples of different physical scenarios of heterogamous media are presented.


Introduction
In the reservoir simulation, there are two numerical methods for solving the governing equations, namely, implicit pressure, explicit saturation (IMPES) method, and the fully implicit method. For each cell, the solution of these equations provides a way to estimate specific time-dependent changes in pressure, saturation, and other key parameters. Although IMPES is much faster than the fully implicit method, it is much less stable because of the explicit evaluation of mobility terms and capillary pressure at the previous time step n instead of the current one n + 1 which requires smaller time steps. On the other hand, as more equations have to be solved in the fully implicit scheme, it requires more CPU memory. Since the fully implicit method is more stable than IMPES, the simulator can take bigger time steps. The ongoing research into IMPES mainly focuses on the improvement of its stability (such as using iterative schemes) and so keeping it fast and robust. The model of two-phase flow in porous media can be solved numerically using a fully implicit scheme [1][2][3], or using an implicit-explicit (IMEX) one [4,5]. The implicit schemes are unconditionally stable but they are time-consuming and expensive from a computational point of view. Therefore, despite IMEX schemes being conditionally stable, they are usually chosen to solve large-scale systems such as oil reservoir models using the IMPES approach. The IMPES schemes have been improved many times (e.g., [6,7]).
The iterative IMPES method is considered to be one of the improved versions of IMPES for two-phase flow [8][9][10]. The IMPES approach divides the system of the equation into the pressure equation and saturation equation which are solved sequentially. Iterative IMPES uses IMPES instead of Newton iteration as an iterative scheme for full implicit systems. In some iterative schemes [11,12], the pressure equation is implicitly solved and for each iteration, an implicit saturation equation with the implicit capillary pressure is constructed and solved. The iterative method was used along with the IMPES scheme by Kou and Sun [13] and El-Amin et al. [14]. In [13], an iterative version of the IMPES scheme was introduced. Kou and Sun [13] linearized the capillary pressure function to be used to couple the implicit saturation equation into the implicit pressure equation. In [14], an iterative implicit pressure, explicit saturation, implicit concentration scheme was developed to solve the model of nanoparticle transport with two-phase flow in porous media.The linearized capillary pressure function was also used to couple the implicit saturation and pressure equations, while the concentration equations were treated implicitly. Recently, El-Amin et al. [15] introduced a convergence analysis of the nonlinear iterative method for two-phase flow and nanoparticle transport in porous media.
The non-isothermal two-phase flow in porous media appears in many natural and industrial applications such as steam injection in oil recovery, geothermal reservoirs, and other environmental and engineering applications [16][17][18].
In the current work, we develop an iterative implicit pressure, explicit saturation, implicit temperature (I-IMPES-IMT) scheme to simulate the problem of non-isothermal two-phase flow in porous media. To the best of the author's knowledge, no research develops a similar scheme to simulate the non-isothermal two-phase flow in porous media. The remaining sections of this paper are arranged as follows: in Section 2, we develop a mathematical model consisting of a group of differential equations, algebraic constraints, and initial/boundary conditions. Then, in Section 3, we introduce an iterative scheme for the governing equation. The spatial discretization is presented in Section 4. In Section 5, error estimation and solution verification are discussed and some numerical examples are listed. Finally, the conclusion is drawn in Section 6.

Modeling and Mathematical Formulation
Consider the problem of two-phase flow with convective heat transfer in heterogeneous porous media. The flow is assumed to be incompressible and immiscible without phase transition. The variation in density is caused by differences in temperature and is related to the buoyancy effect which is described by the Boussinesq approximation. The basis of this approximation is that there are flows in which the temperature varies slightly, in these, the density varies slightly, and this buoyancy contributes the motion. Therefore, the variation in density is neglected anywhere except in the buoyancy term. Let ρ α,0 be a reference density of the phase α, the temperature is T r . For small temperature difference we can write where ρ α,0 is a reference density of the phase α. β α is the coefficient of volume expansion. β α is in the range 10 −3 -10 −4 . For a temperature variation of moderate amount, we have Thus, the variation in the density and therefore in the specific volume is negligible. The velocity equation of the phase α is obtained by Darcy's law, where k rα is the relative permeability of the phase α; K is the permeability tensor; µ α is the α-fluid viscosity; p α is the phase-α pressure; g is the gravitational acceleration; z is the depth; q α is the external mass flow rate; w refers to the wetting phase; and n refers to the non-wetting phase. The temperature difference may cause a convective flow, which is represented by the Boussinesq approximation in the body force term of the previous Equation (1). For measures of central tendency (mean, mode, median), pore-throat sizes (diameters) are generally greater than 2 mm in conventional reservoir rocks, they range from about 2 to 0.03 mm in tight-gas sandstones, and range from 0.1 to 0.005 mm in shales [19]. On the other hand, the ratio between the inertial forces and the viscous forces driving the flow is computed by the Reynolds number, which is used as a criterion to distinguish between the laminar flow, the turbulent flow, and the transition zone. For porous media, the Reynolds number is defined as Re = ρu i d i /µ, where u i is the injection velocity and d i is the injection width of the inlet. According to Bear [20], Darcy's law, which supposes a laminar flow, is valid for Reynolds number less than 1, but the upper limit can be extended up to 10. This may be because the Reynolds number is calculated based on the discharge (inlet) flow rate, however, the velocity reduces significantly in the porous media as it depends on the permeability. For the current problem, the Reynolds number is about 0.6-1.8 for the injected water depending on the discharge inlet diameter 0.01-0.03 m. Therefore, based on Bear's classification, Darcy's law is valid for our study. Moreover, according to Bear's classification , the nonlinear effects are taken into consideration for 10 < Re < 100 [20]. The flow in the current problem is linear, therefore, no need to consider the Brinkman and Forchheimer effects.
The saturation equation which represents the mass conservation law can be written as where φ [-] is the porosity, and s α is the saturation of the phase α. If we sum the two equations of saturation, (2), namely, α = w and α = n, whilst implementing the constraint s w + s n = 1, we reach the following pressure equation where u t = u w + u n is the total velocity; p c = p n − p w is the capillary pressure; ∆ρ 0 = ρ n,0 − ρ w,0 is the density difference; λ α = k rα /µ α is the phase mobility; λ t = λ w + λ n is the total mobility; and q t = q w + q n is the total source mass transfer. Equation (3) with the following saturation equation are solved simultaneously to compute the pressure of water phase: Therefore, we can update the saturation using the following equation: where The relative permeability of the two phases is defined as where S = s w −s wr 1−s nr −s wr is the normalized wetting phase saturation, such that, s nr is the residual non-wetting saturation, s wr is the irreducible saturation, 0 ≤ S ≤ 1. In addition, k 0 rw = k rw (S = 1) , k 0 rw = k rw (S = 1).
Assume that the two fluid phases and the solid phase are of local thermal equilibrium of low flow velocities and relatively small grain sizes (cf. Hassanizadeh [21]): This yields the conservation equation of energy for a two-phase system, consisting of a solid phase plus fluid phases [16,22]: where h α,d is the deterministic thermal conductivity of the phase α, while the mechanical thermal-dispersion tensor, h disp , may be written as Therefore, where d l,w is the longitudinal dispersion coefficient and d t,w is the transverse dispersion coefficient. The capillary pressure is defined by where p e (T) is the entry pressure which is a function of temperature [23], is the surface tension.
The initial and boundary conditions are given by u t · n = q N , s w (or S n ) = S N , where Ω is the computational domain; Γ D is the Dirichlet boundary; Γ N is the Neumann boundary, such that Γ D ∩ Γ N = Φ. n is the outward unit normal vector to the boundary; p D is the pressure on Γ D ; and q N the inflow rate on Γ N .

Iterative Method
The governing Equations (3)- (6), are strongly coupled and highly nonlinear. In order to treat such a type of complexity, we may use an iterative method [24]. The time interval [0, T] is discretized into a number of N T time steps, ∆t n = t n+1 − t n , n + 1 refers to the current time step, and n refers to the previous time step. Apply iterations of N I iterations for each time step, such that k + 1 refers to the current iteration step and k refers to the previous iteration step. The time discretization of time derivatives is achieved by using the backward Euler. The discretized pressure equation takes the form The capillary pressure p n+1,k+1 c is linearized as Using the saturation equation to find The above three equations are solved to compute the wetting phase pressure. Then, saturation is updated using the explicit equation where u n+1,k+1 a = −λ n+1,k t K∇p n+1,k+1 w + gρ w,0 (1 − β w (T n+1,k − T r ))∇z.
In the above equation, the advection term is multiplied with the relaxation factor, α r . In order to ensure the stability of the proposed scheme, we employ the stability Courant-Friedrichs-Lewy (CFL) condition (CFL < 1) to estimate a relaxation factor. The convection terms in the discretized saturation and temperature equations are multiplied by the following relaxation factor: The energy equation is computed implicitly by the following scheme: The advection term in the energy Equation (17) is also multiplied by the relaxation factor, α r . Finally, the variables, λ w ,λ n ,λ t , f w , and ρ α , are updated. This procedure is repeated until approaching convergence.

Spatial Discretization
The cell-centered finite difference (CCFD) is a type of finite volume method and is locally conservative. In the case of rectangular elements, it was proven early on by Wheeler and Yotov [25] and Arbogast et al. [26] to be equivalent to the mixed finite element method. In the CCFD method, the velocity is located on the edges, while the pressure and permeability in the cell centers. Now, we apply the CCFD scheme to the system of Equations (13)- (17), to obtain the full implicit iterative discretization. The discretization form of the pressure equation may be given as In the above equation, the two matrices, A a and A c , depend on the vectors s w and T from the previous iteration step. On the other hand, the matrix A a depends only on the vector T from the previous iteration step. Moreover, the capillary pressure approximation (14) can be discretized as such that the matrix dP s is diagonally defined, for N c total number of cells, by It is clear that the derivative of p c is a function of p c itself, as the saturation at each spatial point varies with time even though it is discontinuously distributed in space.
The spatial discretization of the saturation equation to be coupled with the pressure equation given as such that M is a diagonal matrix appearing from the porosity coefficient. The coupled pressure equation can be obtained from substituting (20) and (21) into (19), i.e., where The upwind scheme is used to treat the convection term of Equation (16) Then, update f n+1,k+1 w i.e., u n+1,k+1 w to be used in the energy Equation (17), Re-arrange Equation (24) to take the form where

Numerical Investigations
In this section, we present error estimations for the time dependent variables, and introduce two numerical examples to check the validity and the efficiency of the proposed scheme. The rate of injection on the inlet boundary is used to calculate the normal velocity, u i = φd i (PVI), here PVI is the pore volume injection and d i is inlet diameter. Table 1 shows the physical parameters definitions and values.
The error estimates for saturation and temperature are presented in Table 2 for different values of the time steps number (TSN) and iterations number (ITN). A domain of size 1.0 × 0.5 m discretized by 40 × 30 is used in this error estimation. Non-iterative calculations are carried out at 2000, 1000, 500, and 100 time steps. In addition, we carried out some calculations for 10, 30, and 50 iterations with 100 time steps. As there is no exact solution for such a complicated problem, we created a reference case of highly dense mesh and a very small time step which is very close to the exact solution. It is under the same conditions but as mentioned above the non-iterative case is just a reference. Table 2 illustrates that the error deceases as the number of time steps and number of iterations increase. It is also clear that the temperature converges faster than saturation. Moreover, in Table 3  In the following, two examples are provided and discussed for heterogeneous media. The first example is a real field high-heterogeneous permeability media [27]. The spatial discretization is taken as 120 × 50 of a uniform rectangular mesh. A hot wetting fluid is injecting at a flow rate of 0.1 pore volume injection (PVI) into the left bottom corner of the domain. Four cases of injections are taken into consideration, i.e., continue until reaches 0.35, 0.5, 0.75, and 0.85 PV (pore volume). In this example, we used a time step size of 0.01 and a maximum number of iterations of 50. Figure 1 shows the contours of the saturation of the wetting phase after 0.35, 0.5, 0.75, and 0.85 PV. It is interesting to observe discontinuities in the saturation distribution caused by the heterogeneity. Moreover, the contours of temperature at the heterogenous medium are plotted in Figure 2. It seems to be that the heterogeneity of the medium has no effect on the temperature distribution.
In the second example, the medium is artificial and heterogeneous, such that the it consists of two subdomains with different permeability distribution (0.1 md for the barrier and 100 md for the rest of the domain); a similar configuration to the first example. The spatial discretization is 120 × 50 of uniform rectangular mesh. An injection rate of 0.1 PVI of hot water at the left bottom corner is used. Four cases of injections are taken into consideration, i.e., continuous until it reaches 0.35, 0.5, 0.75, and 0.85 PV (pore volume). The time step equals 0.01. The distributions of water saturation is shown in Figure 3, while the temperature distribution is shown in Figure 4. Again, one can make a similar observation that the saturation distribution is clearly affected by the medium heterogeneity which is not the case for temperature distribution.    Finally, let us discuss some results to express the effectiveness of the proposed method (see . The domain of this test with a heterogeneous permeability is of size 0.2 × 0.2 m and a grid map of 20 × 20. The time step is 0.14. The results are plotted at NIT = 1 (non-iterative case) and NIT = 5, 10, 15 (iterative cases). In Figure 5, the saturation profiles at y = 0.01 m are plotted against the x-axis for different number of iterations (NIT). The lower graph is zoomed in. It can be seen from this figure, that there is a difference between the non-iterative case and the iterative cases.
In addition, there are small differences between the iterative cases themselves which vanish as the NIT increases. A similar explanation can be presented for the temperature which is introduced in Figure 6, where the temperature profiles at y = 0.04 m are plotted against the x-axis for different number of iterations. The corresponding contours of saturation and temperature are presented, respectively, in Figures 7 and 8.

Conclusions
This paper dealt with the issue of non-isothermal two-phase flow in porous media. The model consists of an equation of flow that includes a term of body force represented by the approximation of Boussinesq and capillary forces, in addition to saturation and energy equations. We develop an iterative implicit pressure, explicit saturation, implicit temperature scheme to solve the system of the governing equations. In order to control the solution convergence, a relaxation factor is formulated as a function of the CFL stability condition. The pressure and saturation equations are coupled via a linearized capillary pressure. Then, the saturation equation is explicitly solved, while the energy equation is calculated implicitly. Both porosity and permeability are updated after each iteration loop. To demonstrate the performance of the current approach, some error estimates are presented for various values of number of time steps and number of iterations. Two numerical examples are provided. It was found that the medium heterogeneity has a clear effect on the saturation distribution while the opposite is true for the temperature distribution. Finally, the effectiveness of the proposed method has been discussed in terms of the number of iterations.