Application of the Picard–Newton Method to Simulate Enhanced Oil Recovery in Heterogeneous Reservoirs: AWAG Injection Study Using Sequential Gaussian Simulation for Permeability Modeling

1. Introduction

Enhanced Oil Recovery (EOR) techniques have become increasingly critical for maximizing hydrocarbon production from mature and unconventional reservoirs, where primary and secondary recovery methods have reached their economic limits. As global energy demands continue to rise while easily accessible oil reserves decline, the petroleum industry faces mounting pressure to develop more efficient and cost-effective recovery strategies [1,2,3]. Among the various EOR approaches, including chemical, thermal, and miscible gas injection methods, Water Alternating Gas (WAG) injection has emerged as one of the most widely implemented and technically proven techniques for enhancing oil recovery in both onshore and offshore operations [4,5].

The WAG process addresses fundamental limitations inherent in conventional water flooding and gas injection schemes. Pure water injection, while effective in maintaining reservoir pressure and displacing oil, often suffers from viscous fingering and poor sweep efficiency due to unfavorable mobility ratios between water and oil phases. Conversely, gas injection, particularly CO₂ flooding, can achieve excellent displacement efficiency through miscible or near-miscible interactions with reservoir fluids, but is prone to gravity override and channeling effects that significantly reduce volumetric sweep efficiency [6,7]. By alternating the injection of water and gas slugs, WAG injection combines the mobility control benefits of water flooding with the displacement efficiency advantages of gas injection, resulting in improved overall recovery performance [8,9].

The effectiveness of WAG processes depends on numerous factors, including reservoir heterogeneity, fluid properties, injection strategy parameters (such as WAG ratio, slug size, and injection rates), and operational constraints. Reservoir heterogeneity, in particular, plays a crucial role in determining the success of WAG operations, as permeability variations can lead to preferential flow paths that compromise sweep efficiency and ultimate recovery [10,11]. Understanding and predicting the complex interactions between these factors requires sophisticated numerical modeling capabilities that can accurately capture the physics of three-phase flow in heterogeneous porous media.

Indeed, numerical simulation has become an indispensable tool for the design, optimization, and risk assessment of EOR strategies, providing a cost-effective and safe alternative to extensive field trials [12,13]. However, modeling three-phase flow in heterogeneous reservoirs presents significant computational challenges due to the strongly nonlinear nature of the governing equations, the presence of sharp saturation fronts, and the need to accurately represent complex geological features across multiple length scales. The coupled system of mass conservation equations for oil, water, and gas phases, combined with modified Darcy’s law and appropriate constitutive relationships, results in a highly nonlinear system of partial differential equations that requires robust and efficient numerical solution strategies [14,15].

Traditional numerical approaches for reservoir simulation can be broadly classified into two categories: fully implicit methods (FIM) and sequential methods such as IMPES (Implicit Pressure Explicit Saturation). While FIM offers unconditional stability and the ability to handle large time steps, making it suitable for strongly nonlinear problems, it requires the simultaneous solution of large, nonlinear systems that become computationally prohibitive for fine-grid simulations in heterogeneous domains [16,17]. IMPES formulations, on the other hand, provide computational efficiency by decoupling pressure and saturation calculations, but suffer from severe stability constraints and convergence difficulties in the presence of strong nonlinearities, particularly those arising from three-phase flow with unfavorable mobility ratios [18,19].

To address these limitations, significant research efforts have been directed toward developing alternative linearization strategies that can maintain the numerical robustness of fully implicit schemes while achieving the computational efficiency of sequential approaches. Among these emerging methodologies, the Picard–Newton sequential methods have demonstrated considerable promise for simulating complex multiphase flow systems [20,21]. These hybrid approaches strategically combine the convergence stability of Picard iterations with the rapid convergence characteristics of Newton’s method by applying them selectively to different components of the coupled system. As an example, in [20], the authors treat the oil-pressure equation implicitly using Picard iterations, while employing Newton-type schemes for the water-saturation updates, thereby achieving an optimal balance between computational efficiency and numerical stability.

Recent developments in this field have shown the potential of Picard–Newton methodologies for handling highly nonlinear reservoir simulation problems. Debossam et al. [22] successfully applied these methods to simulate three-phase flow in petroleum reservoirs under various production and injection scenarios, including WAG processes, demonstrating agreement with conventional approaches while achieving substantial reductions in computational time. Their work highlighted the particular advantages of these methods for handling the complex phase behavior and saturation-dependent properties that characterize EOR processes.

Building upon these previous contributions, the present study focuses specifically on a comprehensive evaluation of the Picard–Newton sequential method for simulating WAG injection in heterogeneous petroleum reservoirs. The reservoir heterogeneity is represented through permeability fields generated by the Sequential Gaussian Simulation Method (SGSIM), providing a realistic description of geological variability within the simulation domain. We consider the classical quarter five-spot configuration, which represents a fundamental flow pattern widely used in both academic benchmarks and industrial applications due to its symmetric structure and practical relevance for injection-production well arrangements [23,24]. The mathematical formulation encompasses compressible, immiscible three-phase flow (oil, water, and gas) governed by mass conservation principles and modified Darcy’s law, with full consideration of capillary pressure effects and gravitational forces.

The primary objectives of this research are to demonstrate the practical applicability of this approach for large-scale EOR simulations involving complex geological settings. Through comprehensive numerical experiments considering realistic reservoir properties and injection scenarios, this study aims to establish the Picard–Newton methodology as an efficient alternative for WAG simulation in the petroleum industry.

The findings presented herein confirm that the Picard–Newton method is capable of reproducing multiphase flow behavior characteristic of WAG processes. These results suggest that such method represents a valuable tool for the design and optimization of EOR strategies in heterogeneous reservoirs, potentially enabling more detailed and comprehensive reservoir studies.

2. Materials and Methods

This section presents the computational methodology employed to investigate fluid flow behavior in heterogeneous petroleum reservoirs under enhanced oil recovery conditions. The study combines geostatistical field generation techniques with numerical flow simulation to analyze the influence of spatial permeability variations on recovery performance and breakthrough characteristics.

The quarter five-spot well pattern is adopted as the fundamental study configuration, representing a common enhanced oil recovery scenario that facilitates systematic analysis of sweep efficiency and fluid displacement mechanisms in heterogeneous porous media. This configuration provides a controlled framework for investigating the interactions between geological heterogeneity and fluid flow dynamics in petroleum reservoirs.

2.1. Physical-Mathematical Model

This section presents the physical-mathematical model governing multiphase flow in porous media. The model considers three immiscible phases (water, oil, and gas) in an isothermal reservoir. The governing equations comprise mass conservation and Darcy’s law modified for multiphase flow [25]. Gas dissolution in oil is not considered in this work [26].

2.1.1. Conservation of Mass

The mass conservation equation (continuity equation) for flow in porous media is expressed in differential form as [14]:

$${\displaystyle \frac{\partial }{\partial t}}\left({\displaystyle \frac{\varphi S_{\alpha }}{B_{\alpha }}}\right)=-\nabla ·\left({\displaystyle \frac{{\mathbf{u}}_{\alpha }}{B_{\alpha }}}\right)+{\dot{q}}_{{\alpha }_{sc}},$$

(1)

where $\varphi$ is porosity, $B_{\alpha }={\rho }_{{\alpha }_{sc}}/{\rho }_{\alpha }$ is the Formation Volume Factor (FVF), and $S_{\alpha }$, ${\mathbf{u}}_{\alpha }$, ${\rho }_{\alpha }$, and ${\dot{q}}_{\alpha }$ represent saturation, superficial velocity vector, density, and source term for each phase $\alpha$ = w (water), o (oil), and g (gas), respectively. The subscript $sc$ denotes standard conditions.

2.1.2. Darcy’s Law

Darcy’s law, extended to multiphase flow through relative permeabilities, is written as [14]:

$${\mathbf{u}}_{\alpha }=-{\displaystyle \frac{k_{r\alpha }}{{\mu }_{\alpha }}}\mathbf{k}\left(\nabla p_{\alpha }-{\gamma }_{\alpha }\nabla Z\right)$$

(2)

where $\mathbf{k}$ is the absolute permeability tensor, and ${\mu }_{\alpha }$, $k_{r\alpha }$, $p_{\alpha }$, ${\gamma }_{\alpha }={\rho }_{\alpha }g$, and Z are viscosity, relative permeability, pressure, specific weight of phase $\alpha$, and depth, respectively. Here, g represents gravitational acceleration magnitude.

2.1.3. Governing Equations for Three-Phase Flow

Substituting velocity from Darcy’s law into the continuity equation yields three equations with six unknowns: $p_o$, $p_w$, $p_g$, $S_o$, $S_w$, and $S_g$. Therefore, three additional relationships are required: full saturation constraint $S_o+S_w+S_g=1$, and capillary pressure definitions $P_{cow}=p_o-p_w$ and $P_{cgo}=p_g-p_o$. This work employs the $p_o$–$S_w$–$S_g$ formulation [25].

The resulting governing equations are:

$${\displaystyle \frac{\partial }{\partial t}}\left[{\displaystyle \frac{\varphi (1-S_w-S_g)}{B_o}}\right]=\nabla ·\left({\displaystyle \frac{k_{ro}}{{\mu }_o{B}_o}}\mathbf{k}\nabla p_o\right)-\nabla ·\left[{\displaystyle \frac{k_{ro}}{{\mu }_o{B}_o}}\mathbf{k}{\gamma }_o\nabla Z\right]+{\dot{q}}_{osc}$$

(3)

for oil,

$${\displaystyle \frac{\partial }{\partial t}}\left({\displaystyle \frac{\varphi S_w}{B_w}}\right)=\nabla ·\left[{\displaystyle \frac{k_{rw}}{{\mu }_w{B}_w}}\mathbf{k}\left(\nabla p_o-\nabla P_{cow}\right)\right]-\nabla ·\left[{\displaystyle \frac{k_{rw}}{{\mu }_w{B}_w}}\mathbf{k}{\gamma }_w\nabla Z\right]+{\dot{q}}_{wsc}$$

(4)

for water, and

$${\displaystyle \frac{\partial }{\partial t}}\left({\displaystyle \frac{\varphi S_g}{B_g}}\right)=\nabla ·\left[{\displaystyle \frac{k_{rg}}{{\mu }_g{B}_g}}\mathbf{k}\left(\nabla p_o+\nabla P_{cgo}\right)\right]-\nabla ·\left[{\displaystyle \frac{k_{rg}}{{\mu }_g{B}_g}}\mathbf{k}{\gamma }_g\nabla Z\right]+{\dot{q}}_{gsc}$$

(5)

for gas.

2.1.4. Initial and Boundary Conditions

Initial conditions require specification of pressure and saturation fields throughout the domain at $t=t_0$: $p(x,y,z,t_0)=p_0$, $S_w(x,y,z,t_0)=S_{w_0}$, and $S_g(x,y,z,t_0)=S_{g_0}$. When gravitational and capillary effects are considered, fluids distribute in distinct vertical regions with pressures determined by hydrostatic relations and saturations by capillary pressure curves [12,27].

In this study, only no-flow boundary conditions are considered, representing impermeable reservoir limits. These boundaries are implemented as zero normal pressure gradient conditions, ensuring that no mass flux occurs across the domain boundaries [13,27].

2.2. Finite Volume Method Applied to the Three-Phase Model

In the context of the Finite Volume Method, discretization means replacing the continuous domain by partitioning it into a finite number of small blocks, also called cells or finite volumes, with well-defined dimensions and positions so that the properties (pressure and saturation) are evaluated at the center of these blocks [22]. These pressure and saturation values at the center of the blocks represent the approximate average value of these variables and are considered constant within each respective finite volume.

The computational mesh has $n_x$, $n_y$, and $n_z$ cells along the directions x, y, and z, respectively. The dimensions of the domain in the three spatial directions are $L_x$, $L_y$, and $L_z$, and the block sizes satisfy

$$\sum _{i=1}^{n_x}\Delta x_i=L_x,\sum _{j=1}^{n_y}\Delta y_j=L_y,\sum _{k=1}^{n_z}\Delta z_k=L_z.$$

(6)

After establishing the computational mesh, the discretization of the partial differential equations that govern the three-phase flow occurs through integration over each finite volume and over a finite time interval [22]. From Eqs. (3), (4) and (5) we have

$$\begin{array}{cc}\hfill {\int }_{t^n}^{t^{n+1}}{\int }_V{\displaystyle \frac{\partial }{\partial t}}\left[{\displaystyle \frac{\varphi (1-S_w-S_g)}{B_o}}\right]\mathrm{d}V\mathrm{d}t& ={\int }_{t^n}^{t^{n+1}}{\int }_V\nabla ·\left({\displaystyle \frac{k_{ro}}{{\mu }_o{B}_o}}\mathbf{k}\nabla p_o\right)\mathrm{d}V\mathrm{d}t\hfill \\ \hfill & -{\int }_{t^n}^{t^{n+1}}{\int }_V\nabla ·\left[{\displaystyle \frac{k_{ro}}{{\mu }_o{B}_o}}\mathbf{k}\left({\gamma }_o\nabla Z\right)\right]\mathrm{d}V\mathrm{d}t\hfill \\ \hfill & +{\int }_{t^n}^{t^{n+1}}{\int }_V{\dot{q}}_{osc}\mathrm{d}V\mathrm{d}t\hfill \end{array}$$

(7)

for oil,

$$\begin{array}{cc}\hfill {\int }_{t^n}^{t^{n+1}}{\int }_V{\displaystyle \frac{\partial }{\partial t}}\left({\displaystyle \frac{\varphi S_w}{B_w}}\right)\mathrm{d}V\mathrm{d}t& ={\int }_{t^n}^{t^{n+1}}{\int }_V\nabla ·\left({\displaystyle \frac{k_{rw}}{{\mu }_w{B}_w}}\mathbf{k}\nabla p_o\right)\mathrm{d}V\mathrm{d}t\hfill \\ \hfill & -{\int }_{t^n}^{t^{n+1}}{\int }_V\nabla ·\left({\displaystyle \frac{k_{rw}}{{\mu }_w{B}_w}}\mathbf{k}\nabla P_{cow}\right)\mathrm{d}V\mathrm{d}t\hfill \\ \hfill & -{\int }_{t^n}^{t^{n+1}}{\int }_V\nabla ·\left[{\displaystyle \frac{k_{rw}}{{\mu }_w{B}_w}}\mathbf{k}\left({\gamma }_w\nabla Z\right)\right]\mathrm{d}V\mathrm{d}t\hfill \\ \hfill & +{\int }_{t^n}^{t^{n+1}}{\int }_V{\dot{q}}_{wsc}\mathrm{d}V\mathrm{d}t\hfill \end{array}$$

(8)

for water, and

$$\begin{array}{cc}\hfill {\int }_{t^n}^{t^{n+1}}{\int }_V{\displaystyle \frac{\partial }{\partial t}}\left({\displaystyle \frac{\varphi S_g}{B_g}}\right)\mathrm{d}V\mathrm{d}t& ={\int }_{t^n}^{t^{n+1}}{\int }_V\nabla ·\left({\displaystyle \frac{k_{rg}}{{\mu }_g{B}_g}}\mathbf{k}\nabla p_o\right)\mathrm{d}V\mathrm{d}t\hfill \\ \hfill & +{\int }_{t^n}^{t^{n+1}}{\int }_V\nabla ·\left({\displaystyle \frac{k_{rg}}{{\mu }_g{B}_g}}\mathbf{k}\nabla P_{cgo}\right)\mathrm{d}V\mathrm{d}t\hfill \\ \hfill & -{\int }_{t^n}^{t^{n+1}}{\int }_V\nabla ·\left[{\displaystyle \frac{k_{rg}}{{\mu }_g{B}_g}}\mathbf{k}\left({\gamma }_g\nabla Z\right)\right]\mathrm{d}V\mathrm{d}t\hfill \\ \hfill & +{\int }_{t^n}^{t^{n+1}}{\int }_V{\dot{q}}_{gsc}\mathrm{d}V\mathrm{d}t\hfill \end{array}$$

(9)

for gas, where V denotes the volume of a finite-volume cell (control volume) and $t^{n+1}=t^n+\Delta t$.

2.2.1. Spatial Discretization

In Eqs. (7), (8), and (9), we distinguish three types of terms: the accumulation term (partial derivative with respect to time), the terms with the divergence operator, and the source term. A simplified notation is used [22], where the generic finite volume indexed by $(i,j,k)$ is represented by P. For a Cartesian coordinate system with parallelepiped volumes, the direct neighbors are indicated by: W $(i-1,j,k)$, E $(i+1,j,k)$, N $(i,j+1,k)$, S $(i,j-1,k)$, A $(i,j,k-1)$ and B $(i,j,k+1)$.

The interfaces between finite volumes are indicated by lowercase letters: w $(i-{\displaystyle \frac{1}{2}},j,k)$, e $(i+{\displaystyle \frac{1}{2}},j,k)$, n $(i,j+{\displaystyle \frac{1}{2}},k)$, s $(i,j-{\displaystyle \frac{1}{2}},k)$, a $(i,j,k-{\displaystyle \frac{1}{2}})$ and b $(i,j,k+{\displaystyle \frac{1}{2}})$.

In the finite volume method, we use the divergence theorem and Two Points Flux Approximation (TPFA) [22] to approximate the spatial derivatives:

$$\begin{array}{cc}\hfill \Xi {\left(\Gamma \Delta \phi \right)}_P& ={\Gamma }_w\left({\phi }_W-{\phi }_P\right)+{\Gamma }_e\left({\phi }_E-{\phi }_P\right)+{\Gamma }_n\left({\phi }_N-{\phi }_P\right)\hfill \\ \hfill & +{\Gamma }_s\left({\phi }_S-{\phi }_P\right)+{\Gamma }_a\left({\phi }_A-{\phi }_P\right)+{\Gamma }_b\left({\phi }_B-{\phi }_P\right).\hfill \end{array}$$

(10)

After evaluating the integrations and introducing appropriate approximations, Eqs. (7), (8) and (9) become [22]

$${\int }_{t^n}^{t^{n+1}}{\displaystyle \frac{\partial }{\partial t}}{\left[{\displaystyle \frac{\varphi (1-S_w-S_g)}{B_o}}\right]}_P{V}_P\mathrm{d}t={\int }_{t^n}^{t^{n+1}}\Xi {\left[T_o\left(\Delta p_o-{\gamma }_o\Delta Z\right)\right]}_P\mathrm{d}t+{\int }_{t^n}^{t^{n+1}}{q}_{os{c}_P}\mathrm{d}t$$

(11)

for oil,

$$\begin{array}{cc}\hfill {\int }_{t^n}^{t^{n+1}}{\displaystyle \frac{\partial }{\partial t}}{\left({\displaystyle \frac{\varphi S_w}{B_w}}\right)}_P{V}_P\mathrm{d}t& ={\int }_{t^n}^{t^{n+1}}\Xi {\left[T_w\left(\Delta p_o-\Delta P_{cow}\right)\right]}_P\mathrm{d}t-{\int }_{t^n}^{t^{n+1}}\Xi {\left[T_w\left({\gamma }_w\Delta Z\right)\right]}_P\mathrm{d}t\hfill \\ \hfill & +{\int }_{t^n}^{t^{n+1}}{q}_{ws{c}_P}\mathrm{d}t\hfill \end{array}$$

(12)

for water, and

$$\begin{array}{cc}\hfill {\int }_{t^n}^{t^{n+1}}{\displaystyle \frac{\partial }{\partial t}}{\left({\displaystyle \frac{\varphi S_g}{B_g}}\right)}_P{V}_P\mathrm{d}t& ={\int }_{t^n}^{t^{n+1}}\Xi {\left[T_g\left(\Delta p_o+\Delta P_{cgo}\right)\right]}_P\mathrm{d}t-{\int }_{t^n}^{t^{n+1}}\Xi {\left[T_g\left({\gamma }_g\Delta Z\right)\right]}_P\mathrm{d}t\hfill \\ \hfill & +{\int }_{t^n}^{t^{n+1}}{q}_{gs{c}_P}\mathrm{d}t\hfill \end{array}$$

(13)

for gas, where $q_{\alpha s{c}_P}={\dot{q}}_{\alpha s{c}_P}{V}_P$ and $V_P=\Delta x_i\Delta y_j\Delta z_k$.

The transmissibilities $T_o$, $T_w$, and $T_g$ are defined as [12]

$$T_{{\alpha }_f}={\left({\displaystyle \frac{k_x{A}_x{k}_{r\alpha }}{{\mu }_{\alpha }{B}_{\alpha }\Delta x}}\right)}_{f}f=w,e$$

(14)

for the x-direction, where $A_x=\Delta y_j\Delta z_k$ is the interface area. Transmissibilities for other spatial directions are defined analogously.

2.2.2. Time Discretization

For temporal discretization, we evaluate the integrals in Eqs. (11), (12) and (13). The solution advances one time increment per iteration, obtaining unknowns at time $t^{n+1}=t^n+\Delta t$, while all information at time $t^n$ and earlier is known.

For a fully implicit scheme, the approximation of integrals containing transmissibilities is approximated by

$${\int }_{t^n}^{t^{n+1}}\Xi {\left(\Gamma \Delta \phi \right)}_P\mathrm{d}t=\Xi {\left(\Gamma \Delta \phi \right)}_P^{n+1}\Delta t,$$

(15)

and for the source term,

$${\int }_{t^n}^{t^{n+1}}{q}_{\alpha s{c}_P}\mathrm{d}t=q_{\alpha s{c}_P}^{n+1}\Delta t.$$

(16)

After spatial and temporal discretization, Eqs. (11), (12), and (13) result in [22]

$${\displaystyle \frac{V_P}{\Delta t}}{\left[{\displaystyle \frac{\varphi (1-S_w-S_g)}{B_o}}\right]}_P^{n+1}-{\displaystyle \frac{V_P}{\Delta t}}{\left[{\displaystyle \frac{\varphi (1-S_w-S_g)}{B_o}}\right]}_P^n=\Xi {\left[T_o\left(\Delta p_o-{\gamma }_o\Delta Z\right)\right]}_P^{n+1}+q_{os{c}_P}^{n+1}$$

(17)

for oil,

$${\displaystyle \frac{V_P}{\Delta t}}\left[{\left({\displaystyle \frac{\varphi S_w}{B_w}}\right)}_P^{n+1}-{\left({\displaystyle \frac{\varphi S_w}{B_w}}\right)}_P^n\right]=\Xi {\left[T_w\left(\Delta p_o-\Delta P_{cow}-{\gamma }_w\Delta Z\right)\right]}_P^{n+1}+q_{ws{c}_P}^{n+1}$$

(18)

for water, and

$${\displaystyle \frac{V_P}{\Delta t}}\left[{\left({\displaystyle \frac{\varphi S_g}{B_g}}\right)}_P^{n+1}-{\left({\displaystyle \frac{\varphi S_g}{B_g}}\right)}_P^n\right]=\Xi {\left[T_g\left(\Delta p_o+\Delta P_{cgo}-{\gamma }_g\Delta Z\right)\right]}_P^{n+1}+q_{gs{c}_P}^{n+1}$$

(19)

for gas.

Since these equations are nonlinear, a conservative scheme is used in the expansions to avoid material balance errors and numerical instabilities [14]. Following Ertekin et al. [27], after conservative expansions, Eqs. (17), (18), and (19) become [22]

$$C_{op}{\Delta }_t{p}_o+C_{ow}{\Delta }_t{S}_w+C_{og}{\Delta }_t{S}_g=\Xi {\left[T_o\left(\Delta p_o-{\gamma }_o\Delta Z\right)\right]}_P^{n+1}+q_{os{c}_P}^{n+1}$$

(20)

for oil,

$$C_{wp}{\Delta }_t{p}_o+C_{ww}{\Delta }_t{S}_w=\Xi {\left[T_w\left(\Delta p_o-\Delta P_{cow}-{\gamma }_w\Delta Z\right)\right]}_P^{n+1}+q_{ws{c}_P}^{n+1}$$

(21)

for water, and

$$C_{gp}{\Delta }_t{p}_o+C_{gg}{\Delta }_t{S}_g=\Xi {\left[T_g\left(\Delta p_o+\Delta P_{cgo}-{\gamma }_g\Delta Z\right)\right]}_P^{n+1}+q_{gs{c}_P}^{n+1}$$

(22)

for gas, where ${\Delta }_t\phi ={\phi }^{n+1}-{\phi }^n$, and the coefficients are defined as

$$\begin{array}{cc}\hfill C_{op}& ={\displaystyle \frac{V_P}{\Delta t}}\left({\displaystyle \frac{{\varphi }_P^{\prime }}{B_{o_P}^n}}\right)\left(1-S_{w_P}^n-S_{g_P}^n\right)+{\displaystyle \frac{V_P}{\Delta t}}\left[{\varphi }_P^{n+1}{\left({\displaystyle \frac{1}{B_o}}\right)}_P^{\prime }\right]\left(1-S_{w_P}^n-S_{g_P}^n\right),\hfill \end{array}$$

(23)

$$C_{ow}=C_{og}=-{\displaystyle \frac{V_P}{\Delta t}}{\left({\displaystyle \frac{\varphi }{B_o}}\right)}_P^{n+1},$$

(24)

$$C_{wp}={\displaystyle \frac{V_P}{\Delta t}}\left({\displaystyle \frac{{\varphi }_P^{\prime }}{B_{w_P}^n}}\right)S_w^n+{\displaystyle \frac{V_P}{\Delta t}}\left[{\varphi }_P^{n+1}{\left({\displaystyle \frac{1}{B_w}}\right)}_P^{\prime }\right]S_w^n,$$

(25)

$$C_{ww}={\displaystyle \frac{V_P}{\Delta t}}{\left({\displaystyle \frac{\varphi }{B_w}}\right)}_P^{n+1},$$

(26)

$$C_{gp}={\displaystyle \frac{V_P}{\Delta t}}\left({\displaystyle \frac{{\varphi }_P^{\prime }}{B_{g_P}^n}}\right)S_{g_P}^n+{\displaystyle \frac{V_P}{\Delta t}}\left[{\varphi }_P^{n+1}{\left({\displaystyle \frac{1}{B_g}}\right)}_P^{\prime }\right]S_{g_P}^n,$$

(27)

$$C_{gg}={\displaystyle \frac{V_P}{\Delta t}}{\left({\displaystyle \frac{\varphi }{B_g}}\right)}_P^{n+1},$$

(28)

and the terms in Eqs. (20), (21), and (22), evaluated at time $n+1$, still require linearization.

Boundary conditions are implemented using forward or backward difference approximations on finite volumes with half the length of inner volumes in the direction perpendicular to the boundary [22]. The mathematical formulation for evaluating transmissibilities at finite volume interfaces, including the decomposition into geometric, pressure-dependent, and saturation-dependent terms, as well as the use of harmonic and arithmetic averaging and upwind schemes for relative permeabilities, is detailed in Debossam et al. [22]. The same reference also presents the methodology for variable time step determination, covering the adaptive iteration-based strategy with increment and decrement factors, time step adjustment based on convergence behavior, and the use of maximum time step limits to control truncation errors.

2.3. Well-Reservoir Coupling

This section details the implementation of the source term, which represents producing and injecting wells in reservoir simulation. The objective is to derive relationships to evaluate pressure or flow rate within the well [22]. The region where the largest pressure gradients typically occur is located in the vicinity of the well, whose dimensions are considerably smaller than those of a typical cell in the computational grid. Therefore, the main challenge in well–reservoir coupling lies in relating the pressure of the grid block containing the well to the actual wellbore pressure.

The simplest and most used approach relates the pressure in finite volumes containing wells to the well pressure through the flow rate expressed in terms of a productivity index [22]:

$$q_{\alpha s{c}_m}=-J_{{\alpha }_m}\left(p_{{\alpha }_m}-p_{w{f}_m}\right)$$

(29)

where $J_{\alpha }$ is the productivity index, $p_{wf}$ is the well pressure, and m is the index of cells containing at least one well section. For each phase, the flow rate is:

$$q_{os{c}_m}=-J_{o_m}\left(p_{o_m}-p_{w{f}_m}\right)$$

(30)

for oil,

$$q_{ws{c}_m}=-J_{w_m}\left(p_{o_m}-P_{co{w}_m}-p_{w{f}_m}\right)$$

(31)

for water, and

$$q_{gs{c}_m}=-J_{g_m}\left(p_{o_m}+P_{cg{o}_m}-p_{w{f}_m}\right)$$

(32)

for gas.

The productivity index for the multiphase case can be written as [27]

$$J_{{\alpha }_m}=G_{w{f}_m}\left({\displaystyle \frac{k_{{r\alpha }_m}}{{\mu }_{{\alpha }_m}{B}_{{\alpha }_m}}}\right)$$

(33)

where $G_{wf}$ is the geometric factor of the well. For vertical wells [27]:

$$G_{w{f}_m}={\displaystyle \frac{2\pi \sqrt{k_{x_m}{k}_{y_m}}\Delta z_m}{ln\left(r_{e{q}_m}/r_{wf}\right)}}$$

(34)

where $r_{eq}$ is the equivalent radius and $r_{wf}$ is the well radius. The equivalent radius is evaluated using the generalized equation by Peaceman [28], which accounts for medium anisotropy and non-uniform mesh effects:

$$r_{e{q}_m}=0.28\left({\displaystyle \frac{\sqrt{{\chi }_1}}{{\chi }_2}}\right),$$

(35)

where

$$\begin{array}{cc}\hfill {\chi }_1& ={\left({\displaystyle \frac{k_{y_m}}{k_{x_m}}}\right)}^{1/2}\Delta x_m^2+{\left({\displaystyle \frac{k_{x_m}}{k_{y_m}}}\right)}^{1/2}\Delta y_m^2,\hfill \\ \hfill \end{array}$$

(36)

$$\begin{array}{cc}\hfill {\chi }_2& ={\left({\displaystyle \frac{k_{y_m}}{k_{x_m}}}\right)}^{1/4}+{\left({\displaystyle \frac{k_{x_m}}{k_{y_m}}}\right)}^{1/4}.\hfill \end{array}$$

(37)

Since the well typically crosses multiple cells, the production flow of each phase is the sum of respective flows in each layer. The total flow (without wellbore storage effect) is:

$$q_{sc}=-\sum _m\left(J_{o_m}+J_{w_m}+J_{g_m}\right)\left(p_{o_m}-p_{w{f}_m}\right)+\sum _m\left(J_{w_m}{P}_{co{w}_m}-J_{g_m}{P}_{cg{o}_m}\right).$$

(38)

The well pressure in each cell, considering the hydrostatic gradient inside the well (neglecting friction losses), can be expressed in terms of a reference pressure [27]:

$$p_{w{f}_m}=p_{w{f}_{ref}}+{\overline{\gamma }}_{wf}\left(Z_m-Z_{ref}\right)$$

(39)

where $p_{w{f}_{ref}}$ is the reference well pressure at depth $Z_{ref}$ and ${\overline{\gamma }}_{wf}$ is the average specific weight of fluids inside the wellbore. Substituting Eq. (39) into Eqs. (30)–(32):

$$q_{os{c}_m}=J_{o_m}{p}_{w{f}_{ref}}-J_{o_m}{p}_{o_m}+J_{o_m}{\overline{\gamma }}_{wf}\left(Z_m-Z_{ref}\right)$$

(40)

for oil,

$$q_{ws{c}_m}=J_{w_m}{p}_{w{f}_{ref}}-J_{w_m}{p}_{o_m}+J_{w_m}{P}_{co{w}_m}+J_{w_m}{\overline{\gamma }}_{wf}\left(Z_m-Z_{ref}\right)$$

(41)

for water, and

$$q_{gs{c}_m}=J_{g_m}{p}_{w{f}_{ref}}-J_{g_m}{p}_{o_m}-J_{g_m}{P}_{cg{o}_m}+J_{g_m}{\overline{\gamma }}_{wf}\left(Z_m-Z_{ref}\right)$$

(42)

for gas.

Therefore, Eq. (38) becomes:

$$\begin{array}{cc}\hfill q_{sc}& =-\sum _m\left(J_{o_m}+J_{w_m}+J_{g_m}\right)\left(p_{o_m}-p_{w{f}_{ref}}\right)\hfill \\ \hfill & +\sum _m\left(J_{o_m}+J_{w_m}+J_{g_m}\right){\overline{\gamma }}_{wf}\left(Z_m-Z_{ref}\right)\hfill \\ \hfill & +\sum _m\left(J_{w_m}{P}_{co{w}_m}-J_{g_m}{P}_{cg{o}_m}\right).\hfill \end{array}$$

(43)

The source terms in Eqs. (20), (21), and (22) are replaced by Eqs. (40), (41), and (42), respectively, for simulations with specified flow rate ($p_{w{f}_{ref}}$ unknown) or prescribed pressure ($p_{w{f}_{ref}}$ known). For specified flow rate, Eq. (43) becomes an additional equation in the system.

For injection wells, only water or gas is injected, and the well-reservoir coupling equation is similar to Eq. (29):

$$q_{\alpha s{c}_m}=-J_{i_m}\left(p_{{\alpha }_m}-p_{w{f}_m}\right)$$

(44)

where $\alpha$ represents water or gas phase, and $J_{i_m}$ is the injectivity index [29]. The mobility of the injected fluid equals the total mobility in the cell containing the well [27]:

$$J_{i_m}={\displaystyle \frac{G_{w{f}_m}}{B_{w_m}}}\left({\displaystyle \frac{k_{r{w}_m}}{{\mu }_{w_m}}}+{\displaystyle \frac{k_{r{o}_m}}{{\mu }_{o_m}}}+{\displaystyle \frac{k_{r{g}_m}}{{\mu }_{g_m}}}\right)$$

(45)

for water injection, and

$$J_{i_m}={\displaystyle \frac{G_{w{f}_m}}{B_{g_m}}}\left({\displaystyle \frac{k_{r{w}_m}}{{\mu }_{w_m}}}+{\displaystyle \frac{k_{r{o}_m}}{{\mu }_{o_m}}}+{\displaystyle \frac{k_{r{g}_m}}{{\mu }_{g_m}}}\right)$$

(46)

for gas injection.

Following the same development as for producing wells, the source terms for injection wells are:

$$q_{ws{c}_m}=-J_{i_m}\left(p_{o_m}-P_{co{w}_m}-p_{w{f}_{ref}}\right)+J_{i_m}{\overline{\gamma }}_{wf}\left(Z_m-Z_{ref}\right)$$

(47)

for water injection, and

$$q_{gs{c}_m}=-J_{i_m}\left(p_{o_m}+P_{cg{o}_m}-p_{w{f}_{ref}}\right)+J_{i_m}{\overline{\gamma }}_{wf}\left(Z_m-Z_{ref}\right)$$

(48)

for gas injection.

The total injection flow consists of only one phase:

$$q_{wsc}=-\sum _m{J}_{i_m}\left(p_{o_m}-P_{co{w}_m}-p_{w{f}_{ref}}\right)+\sum _m{J}_{i_m}{\overline{\gamma }}_{wf}\left(Z_m-Z_{ref}\right)$$

(49)

for water injection, and

$$q_{gsc}=-\sum _m{J}_{i_m}\left(p_{o_m}+P_{cg{o}_m}-p_{w{f}_{ref}}\right)+\sum _m{J}_{i_m}{\overline{\gamma }}_{wf}\left(Z_m-Z_{ref}\right)$$

(50)

for gas injection.

Analogously, for injection wells, the source terms in Eqs. (21) and (22) are replaced by those in Eqs. (47) and (48), respectively, depending on whether the simulation specifies a flow rate (unknown $p_{w{f}_{ref}}$) or a pressure (known $p_{w{f}_{ref}}$). In the specified flow rate case, Eqs. (49) and (50) become additional equations in the system, accounting for water and gas injection, respectively.

2.4. Picard-Newton Sequential Method

This section investigates a solution methodology that combines the ease of implementation of the IMPES method with the robustness and stability of the fully implicit method for three-phase flow simulations [22]. The starting point is the pressure equation, identical to that used in the IMPES method, which is linearized using Picard iteration:

$$\begin{array}{cc}\hfill & -\sum _{nb}\left(B_{oP}^v{T}_{o_f}^v+B_{w_P}^v{T}_{w_f}^v+B_{g_P}^v{T}_{g_f}^v\right)p_{o_{nb}}^{v+1}+\left(B_{oP}^v{C}_{op}+B_{w_P}^v{C}_{wp}+B_{g_P}^v{C}_{gp}\right)p_{o_P}^{v+1}\hfill \\ \hfill & +\sum _{nb}\left(B_{oP}^v{T}_{o_f}^v+B_{w_P}^v{T}_{w_f}^v+B_{g_P}^v{T}_{g_f}^v\right)p_{o_P}^{v+1}=\left(B_{oP}^v{C}_{op}+B_{w_P}^v{C}_{wp}+B_{g_P}^v{C}_{gp}\right)p_{o_P}^n\hfill \\ \hfill & -\sum _{nb}{B}_{w_P}^v{T}_{w_f}^v\left(P_{co{w}_{nb}}^v-P_{co{w}_P}^v\right)+\sum _{nb}{B}_{g_P}^v{T}_{g_f}^v\left(P_{cg{o}_{nb}}^v-P_{cg{o}_P}^v\right)+B_{oP}^v{q}_{os{c}_P}^{v*}+B_{w_P}^v{q}_{sc{w}_P}^{v*}+B_{g_P}^v{q}_{gs{c}_P}^{v*}\hfill \\ \hfill & -\sum _{nb}{B}_{oP}^v{T}_{o_f}^v{\gamma }_{o_f}^v\left(Z_{nb}-Z_P\right)-\sum _{nb}{B}_{w_P}^v{T}_{w_f}^v{\gamma }_{w_f}^v\left(Z_{nb}-Z_P\right)-\sum _{nb}{B}_{g_P}^v{T}_{g_f}^v{\gamma }_{g_f}^v\left(Z_{nb}-Z_P\right),\hfill \end{array}$$

(51)

where $nb=\{A,N,E,W,S,B\}$ represents the neighboring volumes, $f=\{a,n,e,w,s,b\}$ indicates the corresponding interfaces, and the superscript $v^{*}$ denotes the flow rates whose coefficients are linearized through Picard iterations (superscript v), whereas the pressures, both inside the control volume and in the well, are variables to be computed (superscript $v+1$). In this linearization, according to the adopted notation, the superscript $n+1$ has been omitted from the terms containing v, $v^{*}$, and $v+1$. Equation (51) forms a system $\mathbf{A}\mathbf{x}=\mathbf{b}$, which, when solved, yields the pressure field $p_o^{n+1,v+1}$.

Although the transmissibilities in Eq. (51) vary significantly with saturation, the coefficients associated with them contain the sum of the oil, water, and gas transmissibilities. While the individual transmissibilities can vary greatly, their sum exhibits a smoother variation within a time step, making the method’s stability primarily dependent on saturation evaluations. In the fully implicit method, stability is independent of the time step owing to the implicit linearization of the strong nonlinearities in the saturation-dependent terms.

Since the time step restriction in the IMPES method arises from explicit saturation calculations, improved stability requires fully implicit linearization applied to previously explicit equations, suggesting a sequential solution based on operator splitting. This linearization involves derivatives of residuals with respect to water and gas saturations, which are responsible for the strong nonlinearities.

The implicit linearization for saturation calculations approximates the IMPES method to the fully implicit method, yielding the sequential Picard-Newton method.

2.4.1. Sequential Picard-Newton Method: Simultaneous Saturation Solution

After solving the system in Eq. (51), the pressure is obtained at iterative level $v+1$ (at time $n+1$). At this level, pressure is no longer an unknown variable, and the Jacobian matrix no longer contains partial derivatives with respect to pressure. The system becomes:

$${\mathbf{J}}^v\delta {\mathbf{S}}^{v+1}=-{\mathbf{R}}^v$$

(52)

with $\mathbf{S}={\left({\mathbf{S}}_1,{\mathbf{S}}_2,\dots ,{\mathbf{S}}_N\right)}^T$ and ${\mathbf{S}}_m={\left(S_{w_m},S_{g_m}\right)}^T$. The residue vector $\mathbf{R}$ is similar to that of the fully implicit method but now comprises only two residues chosen from $R_o$, $R_w$, or $R_g$. For simplicity, we select ${\mathbf{R}}_m={\left(R_{w_m},R_{o_m}\right)}^T$. With pressure at $v+1$ known, the residuals are:

$$\begin{array}{cc}\hfill R_{w_P}^v& =C_{wp}\left(p_{o_P}^{v+1}-p_{o_P}^n\right)+C_{ww}\left(S_{w_P}^v-S_{w_P}^n\right)-\sum _{nb}{T}_{w_f}^v\left(p_{o_{nb}}^{v+1}-p_{o_P}^{v+1}\right)\hfill \\ \hfill & +\sum _{nb}{T}_{w_f}^v\left(P_{co{w}_{nb}}^v-P_{co{w}_P}^v\right)+\sum _{nb}\left[T_{w_f}^v{\gamma }_{w_f}^v\left(Z_{nb}-Z_P\right)\right]-q_{w_{scP}}^{v*}\hfill \end{array}$$

(53)

for water, and

$$\begin{array}{cc}\hfill R_{o_P}^v& =C_{op}\left(p_{o_P}^{v+1}-p_{o_P}^n\right)+C_{ow}\left(S_{w_P}^v-S_{w_P}^n\right)+C_{og}\left(S_{g_P}^v-S_{g_P}^n\right)\hfill \\ \hfill & -\sum _{nb}{T}_{o_f}^v\left(p_{o_{nb}}^{v+1}-p_{o_P}^{v+1}\right)+\sum _{nb}\left[T_{o_f}^v{\gamma }_{o_f}^v\left(Z_{nb}-Z_P\right)\right]-q_{o_{scP}}^{v*}\hfill \end{array}$$

(54)

for the oil residue, used for gas saturation calculation.

The submatrices composing the Jacobian matrix are:

$${\mathbf{J}}_{P,m}=\left[\begin{array}{cc}{\displaystyle \frac{\partial R_{w_P}}{\partial S_{w_m}}}& {\displaystyle \frac{\partial R_{w_P}}{\partial S_{g_m}}}\\ \\ {\displaystyle \frac{\partial R_{o_P}}{\partial S_{w_m}}}& {\displaystyle \frac{\partial R_{o_P}}{\partial S_{g_m}}}\end{array}\right],$$

(55)

where $m=\{A,N,E,P,W,S,B\}$.

The system in Eq. (52) represents two equations for each finite volume:

$$-R_{{\alpha }_P}^v={\displaystyle \frac{\partial R_{{\alpha }_P}}{\partial S_{w_P}}}{|}^v\delta S_{w_P}^{v+1}+{\displaystyle \frac{\partial R_{{\alpha }_P}}{\partial S_{g_P}}}{|}^v\delta S_{g_P}^{v+1}+\sum _{nb}{\displaystyle \frac{\partial R_{{\alpha }_P}}{\partial S_{w_{nb}}}}{|}^v\delta S_{w_{nb}}^{v+1}+\sum _{nb}{\displaystyle \frac{\partial R_{{\alpha }_P}}{\partial S_{g_{nb}}}}{|}^v\delta S_{g_{nb}}^{v+1},$$

(56)

where $\alpha =w,o$ represents the water and oil phases, and $nb=\{A,N,E,W,S,B\}$ represents neighboring volumes.

The Picard-Newton sequential method with simultaneous saturation solution follows these steps: solve the system in Eq. (51), then solve the $2N\times 2N$ system in Eq. (52). The saturation fields are updated by:

$${\mathbf{S}}^{v+1}={\mathbf{S}}^v+\delta {\mathbf{S}}^{v+1}.$$

(57)

All properties are then recalculated, and the iterative process continues until convergence.

2.4.2. Sequential Picard-Newton Method: Well-Reservoir Coupling

The pressure calculation in the sequential Picard–Newton method is performed in the same manner as in the IMPES method. When wells with unknown well pressure $p_{w{f}_{ref}}$ are present, additional elements are incorporated into the coefficient matrix. Furthermore, an additional equation is included in the system to compute this unknown variable, and well pressures are obtained simultaneously with reservoir pressures [30]:

$$\begin{array}{cc}\hfill & -\sum _m\left[\left(J_{o_m}^v+J_{w_m}^v+J_{g_m}^v\right)\left(p_{o_m}^{v+1}-p_{w{f}_{ref}}^{v+1}\right)\right]=q_{sc}\hfill \\ \hfill \\ \hfill & -\sum _m\left\{\left(J_{o_m}^v+J_{w_m}^v+J_{g_m}^v\right){\overline{\gamma }}_{wf}^v\left(Z_m-Z_{ref}\right)+J_{w_m}^v{P}_{co{w}_m}^v-J_{g_m}^v{P}_{cg{o}_m}^v\right\}\hfill \end{array}$$

(58)

for production wells,

$$-\sum _m\left[J_{i_m}^v\left(p_{o_m}^{v+1}-p_{w{f}_{ref}}^{v+1}\right)\right]=q_{wsc}-\sum _m{J}_{i_m}^v\left[P_{co{w}_m}^v+{\overline{\gamma }}_{wf}^v\left(Z_m-Z_{ref}\right)\right]$$

(59)

for water injection, and

$$-\sum _m\left[J_{i_m}^v\left(p_{o_m}^{v+1}-p_{w{f}_{ref}}^{v+1}\right)\right]=q_{gsc}+\sum _m{J}_{i_m}^v\left[P_{cg{o}_m}^v-{\overline{\gamma }}_{wf}^v\left(Z_m-Z_{ref}\right)\right]$$

(60)

for gas injection.

After calculating the pressures, both $p_o^{v+1}$ and $p_{w{f}_{ref}}^{v+1}$ are known, and the flow rates evaluated at the iterative level $v^{*}$, which appear in Eqs. (53) and (54), can be computed for both production and injection wells.

The derivatives with respect to saturations in the flow rate terms are:

$$\begin{array}{cc}\hfill {\displaystyle \frac{\partial q_{sc{w}_P}}{\partial S_{w_P}}}{|}^v& =-{\displaystyle \frac{\partial J_{w_P}}{\partial S_{w_P}}}{|}^v\left(p_{o_P}^{v+1}-P_{co{w}_P}^v\right)+{\displaystyle \frac{\partial J_{w_P}}{\partial S_{w_P}}}{|}^v{p}_{w{f}_{ref}}^{v+1}\hfill \\ \hfill & +{\displaystyle \frac{\partial J_{w_P}}{\partial S_{w_P}}}{|}^v{\overline{\gamma }}_{wf}^v\left(Z_P-Z_{ref}\right)+J_{w_P}^v{\displaystyle \frac{\partial P_{co{w}_P}}{\partial S_{w_P}}}{|}^v\hfill \end{array}$$

(61)

for $q_{sc{w}_P}^{v*}$ with respect to $S_w$,

$${\displaystyle \frac{\partial q_{sc{w}_P}}{\partial S_{g_P}}}{|}^v=0$$

(62)

for $q_{sc{w}_P}^{v*}$ with respect to $S_g$,

$${\displaystyle \frac{\partial q_{sc{o}_P}}{\partial S_{w_P}}}{|}^v=-{\displaystyle \frac{\partial J_{o_P}}{\partial S_{w_P}}}{|}^v\left(p_{o_P}^{v+1}-p_{w{f}_{ref}}^{v+1}\right)+{\displaystyle \frac{\partial J_{o_P}}{\partial S_{w_P}}}{|}^v{\overline{\gamma }}_{wf}^v\left(Z_P-Z_{ref}\right)$$

(63)

for $q_{sc{o}_P}^{v*}$ with respect to $S_w$, and

$$\begin{array}{c}\hfill {\displaystyle \frac{\partial q_{sc{o}_P}}{\partial S_{g_P}}}{|}^v=-{\displaystyle \frac{\partial J_{o_P}}{\partial S_{g_P}}}{|}^v\left(p_{o_P}^{v+1}-p_{w{f}_{ref}}^{v+1}\right)+{\displaystyle \frac{\partial J_{o_P}}{\partial S_{g_P}}}{|}^v{\overline{\gamma }}_{wf}^v\left(Z_P-Z_{ref}\right)\end{array}$$

(64)

for $q_{sc{o}_P}^{v*}$ with respect to $S_g$.

For injection wells, the term $q_{sco}$ is not considered since oil injection is not contemplated, only water or gas injection. Therefore, only the following terms need evaluation:

$${\displaystyle \frac{\partial q_{sc{w}_P}}{\partial S_{w_P}}}{|}^v=-{\displaystyle \frac{\partial J_{i_P}}{\partial S_{w_P}}}{|}^v\left(p_{o_P}^{v+1}-P_{co{w}_P}^v\right)+{\displaystyle \frac{\partial J_{i_P}}{\partial S_{w_P}}}{|}^v{p}_{w{f}_{ref}}^{v+1}+{\displaystyle \frac{\partial J_{i_P}}{\partial S_{w_P}}}{|}^v{\overline{\gamma }}_{wf}^v\left(Z_P-Z_{ref}\right)+J_{i_P}^v{\displaystyle \frac{\partial P_{co{w}_P}}{\partial S_{w_P}}}{|}^v$$

(65)

for $q_{sc{w}_P}^{v*}$ with respect to $S_w$, and

$${\displaystyle \frac{\partial q_{sc{w}_P}}{\partial S_{g_P}}}{|}^v=-{\displaystyle \frac{\partial J_{i_P}}{\partial S_{g_P}}}{|}^v\left(p_{o_P}^{v+1}-P_{co{w}_P}^v\right)+{\displaystyle \frac{\partial J_{i_P}}{\partial S_{g_P}}}{|}^v{p}_{w{f}_{ref}}^{v+1}+{\displaystyle \frac{\partial J_{i_P}}{\partial S_{g_P}}}{|}^v{\overline{\gamma }}_{wf}^v\left(Z_P-Z_{ref}\right)$$

(66)

for $q_{sc{w}_P}^{v*}$ with respect to $S_g$.

2.5. Fundamental Properties

The fundamental properties used in this work vary only as a function of pressure in an isothermal flow. The detailed formulations and correlations for these properties are presented in our previous work [22]. Here, we provide a summary of the key relationships.

For slightly compressible fluids (water and oil without dissolved gas), the FVF is calculated as [27]:

$$B_{\alpha }={\displaystyle \frac{B_{\alpha }^0}{1+c_{\alpha }\left(p_{\alpha }-p_{\alpha }^0\right)}},$$

(67)

where $c_{\alpha }$ is the fluid compressibility coefficient and $B^0$ is the FVF at reference pressure $p_{\alpha }^0$. Values provided in Table 1.

For gas, which exhibits higher compressibility than liquid phases, the real gas law is applied. Assuming $Z_{sc}\approx 1$ for standard conditions [31]:

$$B_g={\displaystyle \frac{p_{sc}}{p_g}}{\displaystyle \frac{Z_{g}T}{T_{sc}}},$$

(68)

where $Z_g(p,T,S{G}_g)$ is the gas compressibility factor obtained using the Dranchuk-Abou-Kassem correlation [32], and $S{G}_g$ denotes the specific gravity of gas.

Rock porosity variation with pressure is described by [13]:

$$\varphi ={\varphi }^0\left[1+c_{\varphi }\left(p_o-p^0\right)\right],$$

(69)

where $c_{\varphi }$ is the rock compressibility coefficient, as listed in Table 2.

In this work, isothermal flow is assumed, with reservoir fluids exhibiting distinct viscosities: gas is the least viscous, water has nearly constant viscosity, and oil ranges from water-like to very high viscosities for heavy oils. Viscosity determinations follow Debossam [30] and Debossam et al. [22]. In the absence of experimental data, correlations are used: McCain [33] for water (as a function of temperature, pressure, and salinity), Petrosky and Farshad [34] for oil (accounting for dead, saturated, and undersaturated conditions), and Lee et al. [35] for gas (considering temperature, density, and molar mass). Detailed formulations are available in the original sources.

Relative permeabilities follow power-law correlations [31]:

$$k_{rw}\left(S_w\right)=k_{r{w}_{max}}{\left({\displaystyle \frac{S_w-S_{wi}}{1-S_{or}-S_{wi}}}\right)}^{ew},$$

(70)

$$k_{rg}\left(S_g\right)=k_{r{g}_{max}}{\left({\displaystyle \frac{S_g-S_{gr}}{1-S_{gr}-S_{or}-S_{wi}}}\right)}^{eg}$$

(71)

where $k_{r{w}_{max}}$ and $k_{r{g}_{max}}$ define their respective maximum values, the subscripts i and r indicate the initial and residual values, and the exponents $ew$ and $enw$ must be provided.

For three-phase oil relative permeability, Stone’s first model is employed [36]:

$$k_{ro}=\left[{\displaystyle \frac{k_{row}{k}_{rog}}{k_{row}^0\left(1-{\overline{S}}_w\right)\left(1-{\overline{S}}_g\right)}}\right]{\overline{S}}_o,$$

(72)

where $k_{row}(S_w,S_g)$ and $k_{rog}(S_g,S_w)$ denote the relative permeabilities of oil in the oil–water and oil–gas systems of the modified Corey model [30], respectively; $k_{row}^0$ is the relative permeability from the same model, evaluated using $S_w=S_{iw}$ and $S_g=0$; and ${\overline{S}}_{\alpha }$ denotes the normalized saturation [30].

Capillary pressures are also represented by power laws [31]:

$$P_{cow}=P_{co{w}_{max}}{\left({\displaystyle \frac{1-S_{or}-S_w}{1-S_{or}-S_{wi}}}\right)}^{epcow},$$

(73)

$$P_{cgo}=P_{cg{o}_{max}}{\left({\displaystyle \frac{1-S_o-S_w}{1-S_{or}-S_{wi}}}\right)}^{epcgo}$$

(74)

where the subscript $max$ refers to the maximum values and the exponents are obtained through historical fitting.

The parameters utilized in determining relative permeabilities and capillary pressures are shown in Table 3.

2.6. Heterogeneous Permeability Field Generation

The generation of heterogeneous permeability fields is essential for realistic numerical simulations of petroleum reservoirs, as it captures the inherent spatial variability of subsurface formations. In this study, the Sequential Gaussian Simulation Method (SGSIM) [37] is employed, a well-established geostatistical technique that produces stochastic realizations of reservoir property fields while honoring prescribed spatial statistics represented by variogram models.

The computational domain is discretized using a regular three-dimensional grid of $n_x\times n_y\times n_z$ cells, corresponding to physical dimensions of $L_x\times L_y\times L_z$. The statistical parameters governing the permeability field (in mD) include a mean $k_m$, a standard deviation $\sigma$, and bounds $k_{min}$ and $k_{max}$ to ensure physically realistic limits.

Spatial correlation is described by an anisotropic exponential variogram with range parameters $R_x$, $R_y$, and $R_z$ m in the x-, y-, and z-directions, respectively, representing typical geological layering. The variogram includes a nugget effect $C_0$ and a sill $C_1$, corresponding to micro-scale variability and the total variance of the field.

The SGSIM algorithm operates through a sequential simulation procedure in which values are generated following a random visitation path. Each simulated point uses previously simulated values as conditioning data, ensuring that spatial correlations are preserved across the domain. The method employs a Gaussian transformation, performing computations in the normal-score (Gaussian) space and subsequently applying an exponential back-transformation to the permeability domain, consistent with the log-normal behavior typically observed in permeability data. The procedure is summarized in Algorithm 1. The SGSIM algorithm was implemented in Python with support for multiprocessing, while the main simulations employed a fast Fourier transform–based approach, which efficiently generated large-scale heterogeneous permeability fields.

Algorithm 1 Sequential Gaussian Simulation Method (SGSIM)

1: Initialize permeability field: $k\left(\mathbf{x}\right)=0$ 2: Define random visitation order for all grid points 3: for each point ${\mathbf{x}}_i$ in random order do 4:    Identify conditional points $\mathbf{x}j,k\left(\mathbf{x}j\right)$ within the search neighborhood 5:    if conditional points exist then 6:       Apply simple kriging in log-space to estimate local mean $k_m$ and variance ${\sigma }_i^2$ 7:    else 8:      Set $k_m=km,global$, ${\sigma }_i^2=\sigma {global}^2$ 9:    end if 10:    Draw a Gaussian value: $G_i\sim \mathcal{N}(k_m,{\sigma }_i^2)$ 11:    Back-transform to permeability space: $k\left({\mathbf{x}}_i\right)=exp\left(G_i\right)$ 12:    Enforce bounds: $k\left(\mathbf{x}i\right)=max(min(k\left(\mathbf{x}i\right),kmax),kmin)$ 13:    if point lies near a well location then 14:      Apply well enhancement: $k\left({\mathbf{x}}_i\right)=k\left(\mathbf{x}i\right)·[1+\alpha ·exp(-d/rwell)]$ 15:    end if 16: end for

For each simulated point, the algorithm searches for conditioning data within an ellipsoidal neighborhood defined by the variogram ranges, considering up to $N_{max}$ conditioning points to ensure computational efficiency. When sufficient local data are available, simple kriging is applied to estimate local statistics; otherwise, global parameters are used.

To represent enhanced permeability zones around wells due to stimulation effects, local permeability values are modified using an exponential decay function within a radius of $r_{well}$ cells surrounding each well. Enhancement factors ${\alpha }_{inj}$ and ${\alpha }_{prod}$ are applied for injection and production wells, respectively.

The generated permeability fields are subsequently used as input for flow simulation studies, providing realistic representations of reservoir heterogeneity. A quarter five-spot well configuration is employed, with an injection well at $(x_{iw},y_{iw},z_{iw})$ and a production well at $(x_{pw},y_{pw},z_{pw})$.

For simulation purposes, two media with low permeability field heterogeneity were generated and employed. The parameters used for their generation through the SGSIM algorithm are presented in Table 4, Table 5 and Table 6, grouped by category. Permeability field 1 uses the base parameters as shown, while permeability field 2 employs doubled statistical parameters (Table 5) to investigate the influence of mean permeability on flow dynamics. Well enhancement in the vicinity of the wells was not considered; therefore, the factor $\alpha$ in Algorithm 1 was set equal to zero.

As an illustrative example, Figure 1 shows three-dimensional visualizations of the two permeability fields obtained through the Sequential Gaussian Simulation Method. A threshold corresponding to the 75th percentile was applied to simplify their representation by reducing the number of points included.

3. Numerical Results

In this section, we present numerical results demonstrating the performance of the proposed Picard-Newton sequential method for simulating WAG injection processes in heterogeneous petroleum reservoirs. The investigations focus on the classical quarter five-spot configuration.

The numerical experiments address the primary objective of demonstrating the practical applicability of the proposed method for large-scale EOR simulations. The analysis encompasses scenarios that capture the essential physics of three-phase flow during WAG operations.

3.1. Numerical Verification

The numerical simulator developed in this research has been thoroughly verified and validated against analytical solutions and established benchmarks for three-phase flow in porous media. The verification process included comparisons with analytical solutions for incompressible three-phase flow presented by Barros et al. [38] and Guzmán and Fayers [39], covering both concave and convex relative permeability functions under different injection conditions.

The validation results demonstrate good agreement between the numerical solutions obtained with our simulator and the reference analytical solutions, confirming the accuracy and reliability of the implemented method. The detailed verification and validation analyses, including saturation profiles and comprehensive comparisons, can be found in Debossam’s doctoral thesis [30] and in Debossam et al. [22].

3.2. Mesh Selection

The mesh selection employed for the simulations and permeability field generation with low heterogeneity was based on the mesh refinement studies conducted in Debossam [30] and Debossam et al. [22]. According to the obtained results, a grid with 200 × 200 × 30 finite volumes was determined to be sufficiently refined for three-phase flow simulation using the WAG strategy.

3.3. Water-Alternating-Gas Strategy

The Water-Alternating-Gas injection strategy represents a crucial component of the enhanced oil recovery process, requiring careful consideration of operational parameters to optimize sweep efficiency and oil displacement. The implementation of WAG involves the sequential injection of water and gas slugs, with the alternating pattern designed to combine the favorable aspects of both waterflooding and gas injection while mitigating their individual limitations.

The fundamental parameters governing WAG performance include the WAG ratio, cycle duration, slug sizes, and injection rates. The WAG ratio, defined as the volume ratio of water to gas injected over a complete cycle, significantly influences the mobility control effectiveness and ultimate recovery efficiency. Typical WAG ratios in field applications range from 1:1 to 3:1, with the optimal value depending on reservoir characteristics, fluid properties, and economic considerations.

As the primary goal of this work is to apply the Picard–Newton method to three-phase flow simulation, a representative WAG strategy with a 1:1 water–gas injection ratio was implemented. The total production period was set to 1,200 days, during which three complete injection cycles were performed. Each cycle comprised 200 days of water injection followed by 200 days of gas injection. For the numerical solution of the resulting linear algebraic systems, the Bi-Conjugate Gradient Stabilized (BiCGSTAB) iterative method with preconditioning was implemented, employing both Jacobi and ILU(0) preconditioners [22]. The remaining parameters and properties used in the simulations are listed in Table 7.

3.4. WAG Injection in Reservoirs with Distinct Low-Heterogeneity Permeability Fields

In this subsection, we investigate how differences in low-level permeability variability affect the performance of a WAG injection scenario. Two reservoir models characterized by low heterogeneity are considered. By keeping all operational parameters identical, the analysis isolates the influence of subtle geological variations on flow behavior, sweep efficiency, and overall oil recovery. This comparative assessment highlights the extent to which even moderate differences in permeability distribution can shape WAG displacement outcomes.

3.4.1. Permeability Field 1

This case corresponds to the reservoir model with the lowest level of permeability variability considered in this study. The permeability distribution in Field 1 is relatively uniform, with small deviations around the mean value of 100 mD, resulting in a smooth spatial variation of flow properties. As such, this model serves as a reference to assess WAG performance under conditions where geological heterogeneity plays a limited role in fluid redistribution and sweep behavior.

The oil production rate and cumulative oil production curves are shown in Figure 2a and Figure 2b, respectively. The production rate curve exhibits a distinctive cyclic behavior directly linked to the alternation between water and gas injection. At the beginning of each water injection stage, the production rate increases as the injected water efficiently displaces the remaining oil in the more permeable zones. As the water front advances, however, the production rate gradually declines, reflecting the onset of water breakthrough and the corresponding decrease in displacement efficiency. The subsequent gas injection phase restores part of the production by re-pressurizing the reservoir and improving microscopic displacement efficiency through oil swelling and viscosity reduction. Over the successive cycles, the amplitude of production oscillations decreases, indicating depletion of the most accessible oil and progressive attenuation of reservoir energy.

The cumulative oil production curve displays a continuously increasing trend throughout the simulation period, confirming the overall effectiveness of the WAG approach in enhancing recovery. The slope of the cumulative curve is steeper during the early stages of water injection, suggesting higher short-term displacement efficiency, and becomes gentler during gas injection, when production stabilizes at lower rates. This pattern highlights the complementary nature of the two injection fluids: water improves macroscopic sweep, while gas enhances microscopic recovery and maintains reservoir pressure.

Figure 3a,b illustrate the temporal evolution of water and gas saturations at the injector well. The results clearly demonstrate the cyclic behavior characteristic of the alternating injection scheme. During the water injection phases, water saturation at the injector well increases progressively from approximately 0.2 to 0.8, while gas saturation decreases to near-zero values. Conversely, during the gas injection phases, gas saturation rises to approximately 0.6, while water saturation declines correspondingly. The saturation profiles exhibit smooth transitions between injection cycles, with well-defined cyclic patterns that repeat consistently across the three WAG cycles, demonstrating the successful implementation of the alternating injection strategy.

Overall, the results for Field 1 demonstrate that the WAG 1:1 injection sequence provides a balanced mechanism between sweep efficiency and displacement efficiency in a moderately permeable heterogeneous system. The alternation between water and gas successfully delays early gas breakthrough and mitigates rapid water channeling, leading to sustained oil recovery over multiple cycles. The gradual reduction in production rates across cycles indicates depletion of the movable oil fraction, suggesting that optimization of the WAG ratio or cycle duration could further enhance recovery performance in low-permeability heterogeneous formations.

3.4.2. Permeability Field 2

Field 2 exhibits a higher average permeability ($k_m$=200 mD) and a slightly greater degree of spatial variability compared to Field 1, although still representative of a relatively low-heterogeneity reservoir. The broader distribution of permeability values introduces more pronounced contrasts in local flow capacity, which influence the multiphase displacement patterns during the alternating water and gas injection stages. This configuration enables assessment of how a modest increase in permeability and heterogeneity affects the sweep process, fluid mobility control, and the resulting oil recovery under identical WAG operating conditions.

Figure 4a shows the oil production rate, while Figure 4b presents the cumulative oil production. The production rate displays a sharp initial increase, reaching a peak during the early water injection stage, followed by a rapid decline. This behavior is consistent with early water breakthrough in systems containing preferential flow channels, where higher-permeability regions promote faster fluid displacement toward the producer. After the first cycle, the production rate stabilizes at lower but more uniform levels, with mild oscillations reflecting the alternation between water and gas injection. Subsequent WAG cycles help mitigate further production decline, but do not recover the initial peak rates observed at the start of the process.

The cumulative oil production increases monotonically throughout the simulation period. The most significant gain occurs during the first cycle, which contributes the largest fraction of recovered oil. Later cycles provide smaller incremental increases, indicating diminishing returns as reservoir depletion progresses and as the remaining oil becomes increasingly trapped in low-permeability zones. The alternating pattern observed in the production curves is reflected in the injector well saturation profiles (Figure 5a,b), which exhibit three clearly defined cycles of water and gas dominance corresponding to the imposed injection schedule.

3.4.3. Comparative Analysis

When the results are compared between Field 1 ($k_m$=100 mD) and Field 2 ($k_m$=200 mD), clear differences emerge in the overall recovery behavior under identical WAG operating conditions. The higher average permeability in Field 2 enhances global injectivity and facilitates the displacement of oil toward the producer, resulting in greater cumulative recovery. However, this improved flow capacity also promotes the formation of preferential high-permeability channels, which accelerate early water breakthrough and locally reduce sweep efficiency. Consequently, while Field 2 achieves higher total production, its recovery dynamics are characterized by stronger initial peaks followed by faster decline rates, reflecting the typical trade-off between injectivity and sweep control in more permeable heterogeneous systems.

Field 1, with its more uniform permeability distribution, exhibits more stable production rates throughout the simulation period, with gentler oscillations between WAG cycles. The cumulative production curve for Field 1 shows a more gradual and sustained increase, indicating better sweep efficiency and delayed breakthrough. In contrast, Field 2 demonstrates higher initial productivity but experiences more pronounced early breakthrough effects, leading to steeper production decline after the first cycle.

These results underscore the sensitivity of the WAG process to permeability distribution and highlight the importance of geological characterization in reservoir management. The findings suggest that optimizing the injection cycle or adjusting the WAG ratio could further improve displacement efficiency across different reservoir conditions, particularly in higher-permeability heterogeneous systems where preferential channeling poses a greater challenge to sweep efficiency.

4. Conclusions

This work presents a sequential Picard–Newton method for simulating three-dimensional, three-phase flows in heterogeneous petroleum reservoirs under isothermal conditions. The formulation accounts for capillary pressure and gravitational effects, employs finite-volume discretization with implicit time integration, and incorporates well–reservoir coupling models for both injection and production operations. The resulting nonlinear systems were solved using the preconditioned BiCGSTAB iterative method.

The proposed methodology was applied to WAG injection scenarios in a quarter five-spot configuration, considering two heterogeneous permeability fields with mean values of 100 mD and 200 mD. The numerical results successfully reproduced the characteristic cyclic behavior of WAG displacement, including alternating saturation profiles at the injector well and oscillatory production responses corresponding to water and gas injection phases. The comparative analysis revealed that reservoir permeability distribution plays a critical role in determining recovery efficiency: while the higher-permeability field (200 mD) achieved greater cumulative oil production due to enhanced injectivity, it also exhibited earlier water breakthrough and steeper production decline, illustrating the trade-off between sweep efficiency and mobility control in heterogeneous systems.

The developed simulator demonstrated robustness in handling the strongly coupled, nonlinear nature of multiphase flow equations in heterogeneous media, accurately capturing the interplay between geological variability and fluid displacement dynamics. The results confirm that even moderate differences in permeability distribution significantly influence WAG performance, suggesting that reservoir characterization and injection strategy optimization are essential for maximizing recovery efficiency.

Future extensions of this framework may include investigation of alternative WAG ratios and cycle durations, incorporation of thermal effects for non-isothermal displacement processes, implementation of compositional formulations for gas injection with miscibility effects, and application to field-scale reservoir models with complex geological features. Such developments would further enhance the predictive capability of the simulator and broaden its applicability to a wider range of enhanced oil recovery scenarios.