Implicit linearization scheme for nonstandard two‐phase flow in porous media

In this article, we consider a nonlocal (in time) two‐phase flow model. The nonlocality is introduced through the wettability alteration induced dynamic capillary pressure function. We present a monotone fixed‐point iterative linearization scheme for the resulting nonstandard model. The scheme treats the dynamic capillary pressure functions semiimplicitly and introduces an L‐scheme type stabilization term in the pressure as well as the transport equations. We prove the convergence of the proposed scheme theoretically under physically acceptable assumptions, and verify the theoretical analysis with numerical simulations. The scheme is implemented and tested for a variety of reservoir heterogeneities in addition to the dynamic change of the capillary pressure function. The proposed scheme satisfies the predefined stopping criterion within a few number of iterations. We also compared the performance of the proposed scheme against the iterative implicit pressure explicit saturation scheme.

wetting phase saturation. This implies that a nonlocal capillary diffusion term in time is introduced in a two-phase flow model. These all impose an additional complexity onto the standard two-phase porous media flow model.
Due to the nonlinearity and dynamic heterogeneity of the designed model, it is impossible to derive analytic solutions. As a consequence, a numerical approach is the only option to predict such flow dynamics. However, developing efficient algorithms for finding numerical solutions is also a challenge in itself even for standard models. 15 Besides the nonlinearity and heterogeneity of the designed model, long-term temporal dynamics adds an extra difficulty for proposing a reliable numerical model. Implicit discretization in time has been employed to handle long-term subsurface evolution as it allows large time step sizes. Newton-type methods are usually applied to solve the resulting nonlinear system of equations. These approaches are second-order convergent. However, this order of convergence comes at a price of a costly computation of the Jacobian of a system at each time step. [16][17][18][19] In addition, these methods are only locally convergent. 19,20 However, the Newton method can be improved by line search strategies as, for example, Armijo's rule.
The other alternative approach is the splitting and then coupling (splitting-coupling) scheme. It splits the entire system into subsystems. The decomposed subproblems are then solved sequentially and are coupled by data exchanges at each time step. The implicit pressure explicit saturation (IMPES) scheme is a widely used splitting-coupling approach to model two-phase flow and component transport processes. 5,[17][18][19][20][21][22] IMPES solves the pressure equation implicitly and updates the saturation explicitly. This approach eliminates the nonlinear terms in the pressure and saturation equations by evaluating them at saturation and fluid properties up-winded from the previous time step. As a consequence, the scheme is conditionally stable, and hence it requires a sufficiently small time step size to approximate the solution.
Several techniques can be implemented to improve the IMPES approach. A very straightforward approach imposes a large time step for the pressure and then subdivides the time step size for the transport equation. 17,22 This approach relies on the assumption that the reservoir pressure changes slowly in time compared with saturation evolution. The other approach solves the transport equation implicitly using a Newton method while the pressure is treated in the same way as the classical IMPES. 23,24 In Kou and Sun, 20 the capillary pressure function in the pressure equation is approximated by a linear function. This helps to couple the pressure and saturation equations at the current time step. However, the scheme involves the calculation of matrix inverse and multiple numbers of matrix multiplications, which greatly increases the computational cost of the scheme. Furthermore, the transport equation is still solved explicitly in time, and the scheme is reduced to the classical IMPES when the capillary pressure is neglected.
Iterative coupling techniques are also applied to improve the classical IMPES scheme. For instance, in Reference 21 an iteration between the pressure and saturation equation is introduced. This iterative scheme is based on their previous work. 20 Radu et al, 19 have proposed a fixed-point iterative scheme for two-phase flow model (in global pressure formulation).
Recently, Kvashchuk and Radu, 18 have proposed an iterative linearization scheme for two-phase flow (in average pressure formulation) following IMPES. The scheme approximates the capillary pressure function by applying a chain rule and evaluating the nonlinear terms at the previous iteration. This approximates the transport equation semiimplicitly. However, the pressure equation was evaluated at the previous iteration saturation profile. This implies that the scheme lacks a coupling term at the current time step. As a consequence, the scheme might be challenged by dynamic capillary pressure forces that change the saturation distributions in a very short time.
In this article, we propose and analyse an iterative linearization scheme for the designed nonstandard model above based on an iterative IMPES approach, typically we followed the work of Kvashchuk and Radu. 18 We discretize the dynamic capillary pressure functions semiimplicitly in time, where the gradient of the dynamic capillary pressure function (in the pressure and saturation equations) is reformulated by applying the chain rule (see Equation (14) in Section 3.2). We then introduce an iteration step and evaluate the nonlinear terms at the previous iteration. We further introduce an L-scheme type 19,25 stabilization term in the pressure and transport equations. We prove the convergence and robustness of the proposed scheme under natural assumptions. The convergence proof shows that the linearization technique and the introduced stabilization terms allowed the scheme to take a large time step size. By contrast to the classical Newton method, the proposed scheme can be seen as an inexact Newton method which has the advantage of not computing the Jacobian of the system. However, this article is not intended to compare the proposed scheme with the existing linearization schemes including the Newton method. Our colleagues [25][26][27] have done comparison studies on the performance of linearization techniques, and concluded that the fixed-point methods are slower but robust than the Newton method.
This article is organized as follows. Section 2 describes the mathematical model of nonstandard immiscible incompressible two-phase flow in porous media. In Section 3, we introduce a linearization scheme for the resulting model, and prove the convergence of the proposed scheme. We further discuss the choice of a relaxation factor in this section. Numerical simulations in 2D and 3D models are presented in Section 4. This section shows the performance of the proposed scheme and compares it with iterative IMPES. The article ends with a conclusive remark in Section 5.

NONLOCAL TWO-PHASE FLOW MODEL
where ∶ Ω × [0, T] → R is phase mobility, ∶ Ω × [0, T] → R is phase density that controls the buoyancy force, and g is the gravitational vector. The phase mobility is defined as = Kk r , where K ∶ Ω → R d×d is the absolute permeability of the rock, k r is phase relative permeability, and is phase viscosity. For each phase ∈ {w, o}, the balance of mass for the incompressible immiscible fluids yield the transport equations, where is the porosity of the medium Ω, and f is source or sink term in each phase. From model (1) and (2), we obtained two equations with four unknown variables. To close the system the following constraints must also be satisfied: where P c is the capillary pressure that relates the phase saturation to the phase pressures difference. Equations (1) to (3) with appropriate initial and boundary conditions are used to describe two-phase flow dynamics in a porous medium.

Model reformulation
Since we are dealing with incompressible fluids and matrix, we can sum up the mass balance models in Equation (2) to get the pressure equation, where tot = w + o is the total mobility. In Equation (4), we have one equation and two unknowns, namely, P o and S w . As a consequence, the transport equation for the wetting or nonwetting phase saturation should be coupled with Equation (4) in order to close the system. Therefore, we get a system of two equations with two unknowns, where, f t = f w + f n is the total source. In order to solve the two Equations (5a) and (5b), one needs to impose appropriate initial and boundary conditions, such as Neumann and Dirichlet conditions. Thus, we assume that the boundary of the system is divided into disjoint sets such that Ω = Γ D ∪ Γ N . We denote by the outward unit vector normal to Ω, and set where J ∈ R d is phase inflow rate. In order to make the model uniquely determined, it is required that Γ d ≠ ∅.

Relative permeability and dynamic capillary pressure functions
Commonly, the Brooks and Corey 7 and van Genuchten 6 models are used to represent the capillary pressure and relative permeabilities for equilibrium system. For nonequilibrium systems, explicit time-dependency of P c -S w curves have been developed (eg, see References 28,29) to capture changes in capillary pressure induced by dynamic flow conditions. These models are developed under a static wettability condition.
In this article, we consider an extended capillary pressure model that captures the dynamic change of rock wettability at pore-scale. Kassa et al 14 have introduced the dynamic term as an interpolation between the end wetting state curves. This can be described mathematically as follows, where, P ww c and P ow c are end wetting (respectively, the water-wet and oil-wet) capillary pressure functions. Here, the water-wet and oil-wet capillary pressure functions are represented, respectively, with large and small (possibly negative) entry pressures. The dynamic coefficient (⋅) is designed to upscale the dynamics of (pore-scale) time-dependent WA mechanism. In Reference 14 a fluid-fluid contact angle (CA) change model (that changes the wettability from an arbitrary initial wetting state to the final wetting state) was introduced at the pore-level. In Kassa et al, 14 two approaches were considered, namely, uniform and nonuniform WA. The first assumes all pores in the REV has been altered identically through exposure time to the WA agent, whereas the nonuniform WA case considers a CA change in a pore only if that particular pore is exposed to the WA agent. These CA change models were coupled with a bundle-of-tubes model to simulate WA induced capillary pressure curves. The obtained curves and the interpolation model were combined to propose the dynamic coefficient . According to the results in Reference 14, can be represented as: where 1 and 2 are fitting parameters that have a clear relation with the pore-scale CA change model parameter (see the details in Reference 14).
The variable is defined as where T is a predetermined scaling factor, and we recommend to choose T such that ∈ [0, 1], that is, T should be above or equivalent to the life time of the exposure period. The variable is used to measure the averaged exposure time of the REV to the WA agent, and it is an increasing function of exposure time. Thus, the models (7) and (8) describe time-dependent WA induced dynamic capillary pressure model. The change of the capillary pressure in time continues even for constant water saturation. However, keeps constant for pores that are fully occupied with water, that is, S w = 1. In this case, the capillary pressure is only dependent on the current water saturation path. This may lead to a discontinuity of the capillary pressure function at the interface of grid blocks. Thus, the continuity of capillary pressure results in a saturation discontinuity. For end wetting (water-wet and oil-wet) conditions, we considered two consistent sets of capillary pressure functions. Qualitatively, these curves represent either water-wet (ww) or oil-wet (ow) conditions. We adopted the van Genuchten constitutive model for these conditions and can be read as, where P e is phase wetting condition entry pressure, m is pore volume distribution of the porous domain for the 's wetting condition and can be related to n as m = 1∕n . In this study, only the standard van Genuchten relative permeability functions are considered to describe the relative movement of fluids, Coupling the relations k ww r -S w and P c -S w -into the flow model (5a) and (5b) will give nonstandard dynamic two-phase flow model in a porous medium. The goal of this study is to propose a stable and flexible scheme that handles such dynamics efficiently for simulations that consider long-term time evolution.

DISCRETIZATION, LINEARIZATION, AND ITERATIVE COUPLING TECHNIQUE
Let the total simulation time interval [0,T] be divided into N time steps in such a way that 0 = t 0 < t 1 < … < t N = T, and define the time step t = T∕N, as well as t n = n t, n ∈ {1, 2, … , N}.
The backward Euler method is applied to discretize the resulting nonlocal two-phase flow model in time and the semidiscretized model can be read as The superscripts (n+1) and n represent the current and previous time steps, respectively. Above, we omitted the gravity term and the analysis will continue in this form for the sake of clarity and brevity of the presentation.
The above system is fully coupled and challenging to solve directly because of its nonlinearity. Due to the nonlinearity, iterative linearization and sequential coupling methods such as iterative IMPES are needed to solve such systems.

Iterative IMPES
The iterative IMPES linearizes the given two-phase flow problem by evaluating the nonlinear terms from the previous iteration step. 21 Thus, the nonlinear model (12a) and (12b) can be reduced to The iterative IMPES solver starts with S n+1,i w = S n w and thus, the system above is linear and decoupled. Usually the pressure equation (13a) is solved for P n+1,i+1 o first. The computed pressure and the previous iteration saturation profile are used to update the current iteration saturation profile explicitly from Equation (13b). The iteration will continue until the convergence criterion has been satisfied.

Semiimplicit time discretization
The iterative IMPES formulation above splits the pressure and saturation equations in each iteration step. Hence, the approach has missed the inherent coupled nature of the original problem (12a) and (12b). This may cause instability on the convergence of the method in particular for long-term reservoir processes.
In this article, we propose a scheme that couples the pressure, and saturation equations at the (n+1)th time step in addition to the current iteration step. The scheme treats the dynamic capillary pressure function (in the pressure and saturation equations) semiimplicitly in time. We then introduce a monotone fixed-point iteration. 18,19,21,25,30 The development of the scheme is discussed below.
The scheme starts with approximating the capillary pressure function at the current time step (in the pressure (12a) and saturation (12b) equations) by applying chain rule and semibackward Euler discretization in time. The resulting approximation is read as The obtained approximate capillary pressure is substituted back to the two-phase flow model to give the following linear system (we call this linearization technique pseudo-monolithic scheme), The above approach (15a) and (15b) couples the pressure and saturation equations at the current time step weakly. But, importantly, the saturation and pressure state variables communicate each other at the same degree of decision making level. Recall that the variable is also a function of saturation, and thus, the number of equations and unknowns are compatible.
Then stability and accuracy of the pseudo-monolithic scheme (15a) and (15b) is improved further by introducing outer iteration steps (ie, (i + 1) and i), and evaluating the nonlinear terms at the current time step (n+1 instead of n) but at the previous iteration i. We controlled the convergence of the proposed fixed-point iteration by adding an L-scheme type 19,21,25 stabilization term. We named this linearization technique as iterative pseudo-monolithic scheme, and read as where L i + 1 ∈ (0,1] is a stabilization constant that has an important role on the convergence of the proposed scheme. The choice of L i + 1 in each iteration will be discussed later in this article. Equations (16c) and (16d) can be substituted into Equations (16a) and (16b) directly during the solution processes. Here, we note that S n+1,i+1 w and n+1,i+1 are used as a previous iteration values for the next iteration and we set S n+1,0 w = S n+1 w for the first iteration step. Remark 1. The pseudo-monolithic and the iterative pseudo-monolithic schemes reduced to IMPES and iterative L-scheme, respectively, if the capillary pressure is neglected.
Below, we demonstrate the convergence of the iterative pseudo-monolithic scheme, and in Section 4, we compare its performance against the pseudo-monolithic scheme (15a) and (15b) and the iterative IMPES.

Convergence analysis of the iterative pseudo-monolithic scheme
We denote by L 2 (Ω) the space of real valued square integrable functions, and by H 1 (Ω) its subspace containing functions having also the first-order derivatives in L 2 (Ω). Let H 1 0 (Ω) be the space of functions in H 1 (Ω) which vanish on the boundary. Furthermore, we denote by ⟨⋅ , ⋅⟩ the inner product on L 2 (Ω), and by || ⋅ || the norm of L 2 (Ω). L f stays for the Lipschitz constant of a Lipschitz continuous function f (⋅).
Let T h is a regular decomposition of Ω, which decomposes Ω into closed d-simplices; h stands for the mesh diameter. Here, we assume Ω = ∪  ∈ h  . The Galerkin finite element space is given by where P 1 ( ) denotes the space of linear polynomials on any simplex T.
We use the definition of spaces and notations above to write the variational form of Equations (12a) to (12b) which finds P n+1 o , S n+1 w ∈ V h for a given S n w such that the following holds, for all v h ∈ V h . Similarly, we can also write the variational form of the iterative pseudo-monolithic method (16a) and (16b) that find P n+1,i+1 holds for all v h ∈ V h . The aim is to show that the linearized model (19a) and (19b) converges to the nonlinear problem (18a) and (18b) within few outer iteration steps in each time step. The convergence analysis of the scheme is proved theoretically by assuming that the continuous model has a solution. Furthermore, the following assumptions on the coefficient functions and the discrete solutions are defining the framework in which we can prove the convergence of the proposed scheme.

A1:
The mobilities satisfy the Lipschitz continuity condition in the wetting phase saturation, that is, there exist constants L such that This implies that any linear combination of is also Lipschitz continuous. A2: The dynamic capillary pressure function P c , and its partial derivatives Furthermore, we assume that the dynamic capillary pressure P c ( , S w ) is decreasing function, that is, and P c ( ,S w ) < 0 A3: We assumed that the initial wetting phase saturation satisfies ||∇S n w || ∞ ≤ M s with || ⋅ || denoting the L ∞ (Ω)-norm. This implies also ||∇S n+1 w || ∞ ≤ M s and ||∇P n+1 n || ∞ ≤ M p . A4: The total derivative of P c with respect to S w is bounded above by zero. A5: Assume that for any time step (n+1) with n ≥ 0, there exist a solution for saturation S n+1 w and pressure P n+1 o such that the Equations (18a) to (18b) are satisfied.
From now on, we denote by for any v h ∈ V h . Applying the Cauchy-Schwartz inequality followed by the assumptions (A1) to (A4), and testing with v h = e i+1 p , we get the following estimate, Applying (A4) once more, we obtain an estimate as follows, Similarly, we subtract Equation (19a) from Equation (18b) to get, Now by taking the advantage of assumptions (A1) to (A3) and applying the Cauchy-Schwartz inequality with the definition of , Equation (28) can be estimated as, Considering all these and after some algebraic manipulation, the inequality (29) can be rewritten as where s = p + L w M p . Substitute the pressure estimate (27) into (30) to give an estimate for the saturation error: Let us define and apply Young's inequality for > 0 to the inequality (31), and choosing the parameter to be = C At this stage, we can substitute the stabilization term from Equation (16c) into Equation (33) to get the following estimate, For any choice of L i + 1 ∈ (0,1], 1 − 1 L i+1 ≤ 0, and thus, by applying the reverse triangle inequality, we can obtain, Thus, the scheme converges linearly for the designed nonlocal two-phase flow model when is satisfied. ▪ Remark 2. If we choose a small L, convergence of the scheme is guaranteed for large time step. However, the rate of convergence may be slow and thus, we may encounter large number of iterations.

Choice of the relaxation factor
Above we observed that the choice of the relaxation factor plays an important role on the convergence of the scheme.
Here, we introduce a choice strategy for the relaxation factor based on the history of the errors at previous and current iterations. Following Reference 21, we define the length of the residual of the transport equation at the current iteration by The aim is finding a relaxation factor that makes (37) sufficiently small. However, this problem is highly nonlinear optimization problem and thus, challenging to come up with optimal global solution. As a consequence, we compute and bound the relaxation factor adaptively in each iteration.
To support the convergence of the iterative pseudo-monolithic scheme, the relaxation factor L should be chosen such that the residual defined by (37) is decreasing with each successive iteration, that is, From (16c) and (38), there exists a constant L such that, We denote the relaxation factor at the ith iteration step by L i , and thus the relaxation equation for wetting phase saturation can be rewritten as, where L i + 1 ∈ (0,1]. Substituting Equation (40) into Equation (39) and rearranging will give, Recall that L i + 1 ∈ (0,1] and from Equation (38), and thus the choice of L i + 1 should satisfy instead, where L min ,L max ∈ (0,1] and L are specified a priori.

NUMERICAL RESULTS
In this section, we examine the convergence and accuracy of the iterative pseudo-monolithic scheme presented in this work. Section 4.1 presents a comparison between the pseudo-monolithic (15a) and (15b) scheme and the iterative pseudo-monolithic (16a) to (16d) scheme. We also carry out comparisons between iterative IMPES, and iterative pseudo-monolithic scheme in Section 4.2. All the schemes are implemented in the open source software package MRST. 31 Here, we applied two point flux approximation to discretize the models designed below. However, we recall that we applied a Galerkin finite elements to show the convergence of the scheme theoretically. This is to show that the scheme is independent of space discretization methods.

Analytic example
In this subsection, a porous medium flow model is designed by choosing exact solutions followed by constructing source terms and boundary conditions. For this particular example, we set t ∈ [0,1] and Ω = (0, 1) × (0, 1). Furthermore, we consider unit magnitude for rock as well as fluid properties in order to ease the construction of the source terms. We applied van Genuchten relative permeability relations (11) and dynamic capillary pressure We experimented a convergence test considering the inputs above, and Figure 1A presents the number of iterations of the iterative pseudo-monolithic scheme for different time step sizes. Furthermore, we also plotted the number of iterations of the pseudo-monolithic scheme just as a reference.
Obviously, the proposed pseudo-monolithic scheme exits the iteration steps at the first iteration for all time steps. On the other hand, the proposed iterative pseudo-monolithic scheme converges to the solution within two iterations for all time steps except for the larger time step t = 0.2 which needs one extra iteration. Figure 1B shows the associated L 2 error ||S n+1 w − S an w || L 2 (Ω) for the saturation and pressure profiles. S an w represents the analytical solution of the saturation. The pseudo-monolithic method approximates the exact solution efficiently. The iterative pseudo-monolithic scheme has improved the accuracy of the pseudo-monolithic scheme as proposed in Section 3. This explains that the efficiency and accuracy of the iterative scheme can be gained with only the cost of a few extra iterations. Note that the number of iterations can be reduced by considering larger stopping criteria for outer iterations without affecting the accuracy.
We have also performed a numerical experiment to analyze the convergence of the proposed iterative method by fixing the time step size t for a different number of grid cells with where h is the side length of a uniform grid cell. The obtained results are listed in Table 1. The iterative pseudo-monolithic method converges with a maximum iterations of three. This maximum number of iteration was needed for the largest time step size t = 0.2. From Table 1, we observe that the number of iterations keeps the same while the grid size varies. This implies that the proposed iterative pseudo-monolithic scheme is not dependent on the mesh size.

Physical test
Above, we considered an academic example and studied the accuracy of the iterative pseudo-monolithic method over the pseudo-monolithic scheme. In the following, we will compare the iterative pseudo-monolithic scheme and IMPES by considering complex porous media geometries with fluid and rock properties given in Table 2. Note that fluid properties and model parameters given in Table 2 are applied for both iterative IMPES and iterative pseudo-monolithic schemes. The outer iteration loop is allowed to continue until ||S n+1,i+1 w − S n+1,i w || ≤ 2.5 × 10 −5 is satisfied. The relaxation factor choice strategy starts with computing L 1 from (37), where we take ||S n+1,0 w − S n+1,−1 w || = 1 in each time step. The relative permeabilities and capillary pressure models in Examples 1 and 2 below are considering a zero residual saturations for the wetting and nonwetting fluids.

Example 1
The computational domain with 300m × 150m dimensions, consisting of different subdomains for the distribution of permeability, is considered. This porous medium model is shown in Figure 2. Note: The "-" sign stands for the scheme is not convergent, "Total iteration" stands for the overall number of iterations to complete the simulation, and "Average iteration" represents the average iteration number per time step. Abbreviation: IMPES, implicit pressure explicit saturation.
We applied the water-wet van Genuchten relative permeabilities (11) and a capillary pressure function as given below, where P ww c and P ow c are as described in Equation (10). The capillary pressure (44) is changing from P ww c to P ow c dynamically for 7.5 years. Here 7.5 years represent the life of injection for this particular simulation. Further data on the model parameters are given in Table 2 above. We complete the model by injecting the nonwetting fluid to the left-bottom corner of the domain with an injection rate of 0.35m 3 per day and we impose a zero Dirichlet condition at the middle of the right side of the domain. The rest of the boundaries are considered impermeable.
We discretized the above model with 2500 regular grid cells and performed numerical experiments to evaluate the convergence behavior of the iterative IMPES and pseudo-monolithic scheme for different time step sizes. Table 3 shows the convergence results of the two methods. As shown in Table 3, the iterative IMPES only converges if the time step size t ≤ 0.85 day, and the iterative pseudo-monolithic scheme converges for all time step sizes. This shows that the iterative IMPES is subject to strong restrictions with respect to the time step size choice in this example. Usually, IMPES encountered a difficulty regarding the choice of time step size even for standard multiphase flow models. [20][21][22]32 The dynamic nature of the capillary pressure function further worsens the flexibility of iterative IMPES on the choice of the time step size in this example. By contrast, the iterative pseudo-monolithic scheme shows its strength allowing for the relaxed choice of time step sizes. The scheme is capable of taking a large time step size. Furthermore, the total and average number of iterations are increasing and decreasing, respectively, when the scheme considers smaller time step sizes. The decreasing number of average iterations per time step size is a positive sign toward the stability of the iterative pseudo-monolithic scheme.
We further studied the convergence stability of the iterative pseudo-monolithic scheme by controlling the speed of capillary pressure alteration. To do so, we vary the dynamic coefficient parameter 1 from Equation (44). For this numerical experiment, we used the same mesh size as before and chose a larger time step size t = 30.5 days. Table 4 shows the convergence behavior of the iterative pseudo-monolithic scheme for different dynamic coefficient parameter, 1 , values. From Table 4, we observe that the scheme requires more iterations as 1 increases. That means the scheme needs a few extra iterations to converge as the alteration speed of the capillarity becomes faster. Furthermore, the scheme may fail to converge for this model if we choose sufficiently large 1 > 2 × 10 4 (not shown here). For such case, the proposed scheme is enforced to choose a relatively small time step size. Nevertheless, the results above show that the scheme converges successfully for nonlocal (in time) two-phase flow model that considers physically reasonable dynamic capillary pressure alteration (even with capillary pressure jumps). These all support the theoretical convergence analysis of the proposed scheme as discussed in Section 3. In general, the reliability of the scheme to handle the dynamic alteration of the capillary pressure and nonlocality of the problem has been successfully demonstrated.
After a successful convergence stability experiment, we also studied the impact of the dynamic coefficient parameter, As shown in Figure 3, the effect of the dynamic coefficient on fluid displacement is distinct. The movement of the nonwetting fluid is restricted for fast dynamic alteration processes. In other words, the displacing fluid remains the resident fluid for slow dynamic capillary pressure alteration, whereas it swipes the resident fluid when we consider a fast capillary pressure alteration. This happens because the wetting property of the volumes (occupied by the displacing fluid) is altered to the intermediate-wet system before it leaves the volume, and thus, the nonwetting fluid preferred to be in contact with the solids when we consider fast capillary pressure alteration.

Example 2
We considered 50m × 50m × 10m dimensional heterogeneous medium, with permeability distribution, We employed the same relative permeabilities functions as above and dynamic capillary pressure given as where P ww c and P ow c are as described in Equation (10). The capillary pressure is allowed to change from P ww c to P ow c dynamically according to model (46) in each subdomain for 2.5 years. Additional data on model parameters are listed in Table 2. We inject the nonwetting fluid to the west particularly at (y,z) ∈ (10m,15m) × (2m,8m) with an injection rate of 0.15m 3 /day for 2.5 years, and impose a zero Dirichlet boundary condition to the east side of the domain, particularly at (y,z) ∈ (10m,15m) × (2m,8m). The rest of the boundaries are considered to be impermeable. We discretized the designed model above with 3125 grid cells and we did simulations to examine the convergence behavior of the iterative IMPES and pseudo-monolithic linearization techniques. The obtained results are published in Table 5. In Table 5, we noticed that the iterative IMPES fails to converge for time step sizes bigger than 0.26 day. This shows that the choice of a time step size is strongly restricted for iterative IMPES linearization which confirms the results reported in References 20,21,32. Unlike the iterative IMPES, iterative pseudo-monolithic scheme relaxes the choice of the Note: Here "-" represents that the scheme is not convergent. Abbreviation: IMPES, implicit pressure explicit saturation.

TA B L E 6
The impact of dynamic coefficient parameter on the stability of the iterative pseudo-monolithic scheme time step size. The scheme return with an approximate solution for a relatively large time step size compared with the iterative IMPES scheme. The total and average number of iterations, respectively, are increasing and decreasing when we consider smaller time step sizes, see Table 5. We further investigate the sensitivity of the dynamic coefficient parameter 2 . For this, we keep the number of grid elements as before and the time step size to be the larger one in Table 5, that is, t = 90.25 days. Then, we vary 2 , and observe its impact on the convergence of the scheme. Table 6 shows the convergence results for different values of 2 . As shown in Table 6, the scheme converges with the same number of iterations for all values of 2 . This implies that the proposed scheme is not affected by the speed of the capillary pressure alteration dynamics for this particular example.
Above, we studied the convergence of the iterative pseudo-monolithic scheme for the nonlocal two-phase flow model. Below, we investigate the impact of the dynamic capillary pressure model on the injected fluid distribution. Figure 4 compares the injected fluid distribution for the initial wetting condition capillary pressure (ie, 2 = 0 in Equation (46)), and dynamic capillary pressure model (46) with 2 = 1 × 10 5 .
In Figure 4, we observe that the displacing fluid leaves the resident fluid behind when we consider 2 = 0 in Equation (46). This is due to the fact that the rock surfaces are water-wet in this case, that is, no WA, and thus, the resident fluid prefers to remain in the pores. On the other hand, the nonwetting fluid has displaced the resident fluid and concentrated near the injection area when we employed the dynamic capillary pressure model (46) with 2 = 1 × 10 5 . In this case, the wettability of the rock surfaces near to the injection area has been changed (in time) to the intermediate-wet system before the displacing fluid leaves the volume, and thus, the displacing fluid preferred to occupy these pores. That means the injected fluid pressure is able to displace the resident fluid with relatively small pressure. This shows that the dynamic capillary pressure results in a large change of fluid saturation as compared with the standard capillary pressure model, P ww c , in Equation (10). This might be one of the reasons that restrict the time step size choice of the iterative IMPES.

CONCLUSION
In this article, we introduced fluid history and time-dependent dynamic capillary pressure model in a two-phase immiscible incompressible porous media flow model. We developed a linearization scheme for the resulting nonstandard two-phase flow model by treating the capillary pressure implicitly and adding stabilization terms. This implicit treatment of the dynamic capillary pressure model couples the pressure and saturation equations strongly, and makes the scheme stable. We gave a theoretical convergence analysis of the scheme under some meaningful assumptions. The scheme has been successfully implemented and tested for different illustrative examples. We found that the proposed scheme is efficient to approximate the solution of the resulting nonstandard two-phase flow model. Most importantly, the scheme demonstrates flexibility regarding the choice of time step size for dynamic capillary pressure alteration (possibly with capillary jumps). Thus, combining the scheme with a Newton method is a straightforward application. This implies that one can alternate between the iterative pseudo-monolithic scheme and the Newton method as mentioned in Reference 25. This may further improve the convergence speed and accuracy of the approximation to simulate such complex models. Furthermore, the convergence of the scheme shall be investigated for degenerate cases in the future.