Buckley–Leverett Theory for a Forchheimer–Darcy Multiphase Flow Model with Phase Coupling

This paper is dedicated to the modeling, analysis, and numerical simulation of a two-phase non-Darcian flow through a porous medium with phase-coupling. Specifically, we introduce an extended Forchheimer–Darcy model where the interaction between phases is taken into consideration. From the modeling point of view, the extension consists of the addition to each phase equation of a term depending on the gradient of the pressure of the other phase, leading to a coupled system of differential equations. The obtained system is much more involved than the classical Darcy system since it involves the Forchheimer equation in addition to the Darcy one. This model is more appropriate when there is a substantial difference between the phases’ velocities, for instance in the case of gas/water phases, and applications in oil recovery using gas flooding. Based on the Buckley–Leverett theory, including capillary pressure, we derive an explicit expression of the phases’ velocities and fractional water flows in terms of the gradient of the capillary pressure, and the total constant velocity. Various scenarios are considered, and the respective numerical simulations are presented. In particular, comparisons with the classical models (without phase coupling) are provided in terms of breakthrough time among others. Eventually, we provide a post-processing method for the derivation of the solution of the new coupled system using the classical non-coupled system. This method is of interest for industry since it allows for including the phase coupling approach in existing numerical codes and software (designed for solving classical models) without major technical changes.


Introduction
Most applications involving fluid dynamics in laminar flow regime through porous media are described using Darcy's equation. Forchheimer's equation suggests the addition of a quadratic correction to Darcy's equation for high flow values, and it was shown that it provides an alternative model to the basic Darcy's law. In literature, the upper limit of Darcy's law is characterized by a limit of the Reynolds number R e (for instance, see Equation (10) in [1]) that is limited between R e = 1 and R e = 10 (cf., e.g., Bear [2], and Chapman [3]). In this paper, we consider the evolution of a two-phase system with significantly different Reynolds numbers and considerable phase-interface interaction, water/gas for instance. We describe the high flow rate phase velocity and pressure using the Forchheimer equation and the relatively low flow rate phase using Darcy's equation. Therefore, we consider a coupled Frochhiemer-Darcy's system that governs the water and gas phases flow displacement in a porous medium. The modeling novelty here consists of taking into account the effect of one phase on the other by introducing a coupling of the Forchheimer equation and the Darcy one through pressure cross-terms in addition to the closure of the system through the capillary pressure and the saturation. The idea is inspired from the results of Kalaydjian et al. [4] and Guérillot et al. [5] in the Darcy scenario. Thus, our main objective in this paper is to extend the pressure cross-terms classical (or fully Darcian) model to the non-Darcy case. Specifically, we generalize the classical Darcy's equation for the slow phase by adding a pressure cross-term (of the fast phase) and coupling it to a generalized Forchheimer equation for the fast phase, including a pressure cross-term (of the slow phase).
The application of our model ranges from a formation near high-rate oil or gas production, groundwater pumping to liquid-waste injection wells. The literature dedicated to the analysis of single-phase high flow rate is rich, and we refer, non-exhaustively, to Tek et al. [6], Swift and Kiel [7], Lee et al. [8], Skjetne et al. [9], Yu-Shu Wu in [10], Guppy et al. [11,12], Evans et al. [13], and Evans and Evans [14]. In particular, Tek et al. [6] derived Forchheimer's equation to generalize Darcy's law and extend its application to various fluid flow rates. Swift and Kiel [7] studied a gas-well performance based on a Darcy and non-Darcy model for the flow of gas through a porous medium. Lee et al. [8] quantified the "turbulence intensity" as related to deliverability front gas wells based on a dimensionless number called the Forchheimer number. Yu-Shu Wu in [10] showed that the values of the non-Darcy coefficients play a crucial role in the derivation of a non-Darcy model by studying the impact of permeabilities on the flow dynamics. For fractured reservoirs, Guppy et al. [11,12] studied the dynamics of a non-Darcy flow in the case of high-permeable fractured wells. Evans et al. [13,14] analyzed the effect of an immobile/mobile liquid saturation on a non-Darcy flow for single-phase displacement. Other works concerning single and multiphase flow with the non-Darcy regime are discussed in Ahmadi et al. [1], and Lasseux et al. [15].
In summary, this work introduces a new coupled Forchheimer-Darcy's system to model the evolution, behavior, and interaction of a Darcian flow interacting with a non-Darcian one in a porous reservoir. The model being entirely new, we proceed to its full rigorous mathematical analysis to prove its well-posedness that is the existence and uniqueness of solutions in the adequate functional spaces. From a mathematical point of view, the main difficulty consists of dealing with the coupling terms. We overcome this difficulty by recasting the system into a vectorial form. From the applications point of view, we investigate the effect of the introduced coupling numerically by considering several values for the amplitudes of the coupling coefficients, and provide comparisons of the results with the classical Forchheimer-Darcy case (without coupling). Moreover, the numerical analysis of the impact of the Forchheimer coefficient on the fractional flow and the water saturation profile is provided. Eventually, application of our model to the estimate of the breakthrough time is presented emphasizing the potential of our model. Since all the models being used on a daily basis in the oil and gas industry ignore the phase coupling, and therefore developing new codes and softwares to include the coupling is rather heavy from the technical point of view, we propose a post-processing scheme allowing for including the coupling using the latter mentioned and existing codes. Our idea consists of using the classical codes to solve the classical Darcy-Forchheimer system with modified phase mobilities, and then post-processing the solution to take into account the coupling. This paper is organized as follows: Section 2 is devoted to the physical model describing the two-phase fluid flow model and fluids' parameters. Section 3 is dedicated to the rigorous mathematical analysis of the coupled set of differential equations at hand, namely the coupled Darcy-Forchheimer system. In addition, in this section, we derive explicit mathematical expressions of the velocity solutions of the coupled system in addition to the phases' fractional flows to be used for the development of a Buckley-Leverett theory. Section 4 is dedicated to the numerical resolution of the system by solving the associated Buckley-Leverett equation using the Euler scheme and finite volume method. We investigate numerically several scenarios: different amplitudes of the coupling amplitudes, different values of the Forchheimer coefficient, comparison with the classical system (without coupling), etc. Eventually, the last section is dedicated to the presentation of our post-processing method and its theoretical and numerical validation.

Physical Model
In this work, the gravity effects are neglected. We assume that the flow in the fluid domain is governed by the following system (for instance, see [16,17] where µ g and µ w are the viscosities of gas-phase and water-phase, respectively, and k is the absolute permeability tensor. kr g and kr w denote the relative permeability of gas-phase and water-phase, respectively. (p g , v g ), and (p w , v w ) are the pressures and velocities of gas-phase and water-phase, respectively. The parameter β g denotes the Forchheimer's coefficient flow, and ρ g is the density of the gas phase. This study introduces a coupled Forchheimer-Darcy's system to justify the fact that each phase acts on the other one (cf., e.g., [4,5]). First, we introduce the mathematical formulation of the coupled Forchheimer-Darcy system that governs the fluid evolution (cf., e.g., [4,16,17]) where kr g,w and kr w,g denote the pseudo-permeabilities of gas and water phases, respectively. Let us mention that all the permeabilities and pseudo permeabilities are assumed to be functions depending on water saturation only. Now, we introduce the capillary pressure p c,g,w defined as the difference between the gas-phase and the water-phase pressures. That is, p c,g,w = p g − p w .
To compute the fractional flow of the water-phase and derive a mathematical formula of the velocity vectors, we denote by M g and M w the gas and water mobilities. In addition, we introduce the gas and water pseudo-mobilities M g,w and M w,g , respectively, depending on the cross-terms amplitudes of their respective gas and water phases in (2), as We introduce the non-Darcy flow coefficient following [18] and given by where C g is a non-Darcy flow constant, S w is the water-phase saturation, and S wr is the connate saturation of the water-phase. φ is the porosity of the porous medium. Now, we assume that the permeabilities kr g and kr w are as in the Brooks and Corey model [19], that is, for all S wr ≤ S w ≤ 1 − S nr , and we assume that the cross-terms in system (2) are given by kr g,w = (S w − S wr ) n g,1 (1 − S w − S nr ) n g,2 , and kr w,g = (S w − S wr ) n w,1 (1 − S w − S nr ) n w,2 , where n g,1 , n g,2 , n w,1 , and n w,2 are real numbers that can be obtained from laboratory experiments and measurements, and S nr and S wr denote the residual gas-phase saturation and the connate water saturation, respectively. In addition, k max rw and k max rg denote the maximal relative permeabilities of water and gas phases, respectively. We consider the Van Genuchten capillary pressure [20] (assumed function of the water saturation only) whose expression is given by where τ denotes a displacement pressure, and m is a prescribed empirical real number between 0 and 1.

Mathematical Solutions
We are interested in the description of incompressible fluids displacement in a porous medium reservoir. The mass conservation equations for the gas and the water phases, denoted with subscripts g and w respectively, read where ρ g , S g , and ρ w , S w denote the density and saturations of gas and water, respectively. Moreover, v g and v w denote the superficial velocity of the gas and the water phases, respectively.

Well-Posedness Result
In this section, we develop an existence and uniqueness theory of the solutions to the coupled system (2) in general d-dimensional flow. For this purpose, we recast the system (2) as follows: where U is the vector solution, and A and B symbolize the following matrices given by and We complete that matrix system (9) with the following boundary conditions: in Ω, where L is the reservoir length, and n is the unit outer normal vector to the fluid domain boundary denoted by Γ. Moreover, we assume that the terms in the boundary conditions (12) satisfy Now, we denote by V g , V w , W g and W w the following vectors: Thus, we can write the conditions (12) and (13) as follows: In addition, we denote by V g,w and W g,w the following vectors: Before going further, we shall need the following preliminaries. First, we give the following classical inequalities (for instance, see [21,22]): There exists a constant α > 0 such that 3. B * is an isomorphism from X/V onto M * and ||Ṁu|| M * ≥ α||u|| X/V , ∀u ∈ X/V. Now, let M p,q (X) be the space of matrices of p rows and q columns in X. We define the following spaces: and The weak formulation of system (9)-(13) reads and for all q = q 1 g 2 ∈ Y, and for all Φ ∈ X, and for every satisfying (15). Next, define the spaces Proposition 1. Assuming that V g,w ∈ X g,w , W g,w ∈ Y g,w satisfies (15), then there is unique U g,w ∈ X/R g such that The proof of this proposition is similar to the proof of Proposition 3.1 in [24]. Inspired by Proposition 1, we write U = U g,w,0 + U g,w where U g,w,0 ∈ R g . The weak formulation (23) and (24) reads where the operator A is defined as Proposition 2 (see Proposition 3.2 in Audu et al. [24]). Problem (31) is equivalent to problem (23) and (24).
Lemma 2 (see Lemma 4.1 in Audu et al. [24]). The operator A : X −→ Y satisfies the following inequalities: and Lemma 3 (see Lemma 4.2 in Audu et al. [24]). The function U → A(U + U g,w ) defined in (32) is strongly monotone, that is, for all U, V ∈ X, we have Proof. Define the map F as F : X −→ R given by Let U, V ∈ X, for all h ∈ R, we write Consequently, Using Lebesgue convergence theorem (for details, see [25]), we have Equivalently, we have Thus, Setting U = 0, we infer Noticing that F is positive definite and symmetric, we have Let U, V ∈ X, and setÛ = U + U g,w andV = V + U g,w with U g,w ∈ X fixed. Therefore, we have The mean value theorem allows for writing Lemma 4. The function U → A(U + U g,w ) in (32) is coercive in X, for any fixed U g,w ∈ X. Moreover, Proof. Let U ∈ X be arbitrarily chosen, and settingÛ = U + U g,w equivalent to (39), we have In addition, we have Now, let λ 1 = 0 be the least eigenvalue of the matrix A that leads to Thus, Proposition 3. Let U g,w be fixed in X. Then, the function is continuous on R for all U, V, W ∈ X.
We refer to Proposition 4.1 in [24] for a proof. Indeed, summing-up and combining all the previous results, using the inf-sup condition in Lemma 1, and taking into consideration the linearity and continuity of the operator G, one establishes the existence of (U, P) solution of the system (9).

Explicit Solutions
In this section, we provide a mathematical explicit expressions of the velocities solutions to the coupled system (2). We assume that the reservoir is a horizontal linear reservoir with a uniform cross-section A, and we consider only the x-direction components of the velocities for the rest of this work. Thus, system (2) now reads v g + β g ρ g k kr g The following lemma gives an explicit formulation of the velocity solutions of system (54) in terms of the capillary pressure. We have Lemma 5. Denote by (v w , p w ) and by (v g , p g ) the velocities and pressures of water and gas phases, respectively, solutions of the coupled Forchheimer-Darcy's system (54). Then, the velocity v g is given with and To deduce the velocity of water v w , we apply the relation v t = v g + v w , where v t is the total constant velocity (i.e., A v t = q t ).
Proof. Using the definition of the capillary pressure, we can write Replacing the terms in (58), we get 1 + 2β g ρ g kM w,g ∂ x p c,g,w ṽ g + β g ρ gṽ Set the term and ξ is given by the relation (56). We get an equation on the pressure gradient ∂ x p w , by substituting the second equation in the first equation in (59), Plugging the expression of the capillary pressure gradient (63) into system (58), we infer (55) when ξ * is given by (56).
Next, we introduce the concept of fractional flow as proposed in Willhite [26]. The fractional flow is defined as the volumetric flux fraction of the water-phase flowing position x and time t. For instance, the water-phase fractional flow is given by To compute the fractional flow of gas-phase, we use the closure relation f g + f w = 1. Thanks to Lemma 5, we infer The water-phase fractional flow is a function of water saturation S w and capillary pressure gradient to x.

Numerical Simulation
This section is dedicated to the presentation of numerical simulations for comparison purposes. Indeed, we compare the water saturation and the fractional flow profiles, and the front position in several scenarios including the non-coupled situation. In addition, we analyze numerically the impact of a variation of the Forchheimer coefficient on the water saturation profile and the waterfront position. Eventually, as an application of our model, we provide estimates of the breakthrough time and compare it to the classical case among other cases.

Numerical Scheme
To study the water-phase saturation profile, we introduce the Buckley-Leverett equation [27] (also referred as the front advance equation) where q t is the total constant injection rate of the water-phase. Equation (66) is discretized in a one-dimensional space of length L through a finite volume method. We consider a mesh of N + 1 cells x i , of length ∆x = L/N (see Figure 1). First, we integrate the Buckley-Leverett Equation (66) over a grid cell x i , and rearranging the derivation terms, we obtain Denoting S w,i the average water-phase saturation in a single cell x i , (67) leads to The term on the left-hand-side in (68) is discretized using the 1st order accurate Euler forward method: with ∆t being the time step. Equation (69) is then cast to For the first order upwind, the terms f w,i+ The matrix form of Equation (72) reads where and F w is a N × N matrix defined by where f w,0 corresponds to the flow at the virtual cell, which is at distance ∆x to the first grid cell x 1 (see Figure 1), and extrapolated and equals f w,0 = 1. From the numerical point of view, the results of the first-order upwind present some small oscillations, which is not satisfactory. An alternative modified Euler method (cf. e.g., Koenig [28]) is used. This method explores both water-phase saturation and fractional flow of water at the end of each time step. Specifically, wheref n+1 w (S n+1 w ) is computed using the first order upwind given by Equation (73) as follows:S n+1 w The numerical scheme in (76) is stable under the following CFL condition: The convergence of our scheme can be seen on the graphs of the water saturation error curves in Figures 2, 6c and 7c.

Model Parameters
In this short paragraph, we present the numerical parameters we shall use for the numerical simulations below. These numerical values are extracted from the book [10] dedicated to the analysis and numerical simulation of similar models describing the same physical context but without coupling. The initial saturation S 0 wr is set to 0.2, and the outlet pressure denoted by p out w is set to p out w = 80.10 5 Pa. We assume that our reservoir is horizontal to neglect the gravity for simplicity (our model is applicable in case of a non horizontal reservoir) with length L = 3280 ft and cross section A = 10.764 × 10 4 ft 2 . The absolute permeability is set to 300 mD. The viscosity of water-phase and gas phase are equal to µ w = 1 cP, and µ g = 5 cP respectively. The density of water phase, and gas phase are equal to ρ w = 22.653 kg/ft 3 , and ρ g = 28.316 kg/ft 3 , respectively. We set the initial injection rate to q t = 1259.96 bll/s. The Forchheimer coefficient of gas phase is set to C g = 3.2 × 10 −6 m 3/2 (other values will be considered below for comparison purposes). The porosity is equal to φ = 0.25. The connate saturation and the residual saturation of gas-phase are equal to S wr = 0.20 and S nr = 0.30, respectively. The maximum of the relative permeability of water-phase is assumed to be kr max w = 0.5. The maximum of the relative permeability of gas-phase is set to kr max w = 0.8. The parameters n g , n w , n g,1 , n g,2 , n w,1 , n w,2 are set to 2, and the displacement pressure τ is equal to 10 6 Pa. Eventually, the real m is set to 0.5.

Numerical Results
This section is dedicated to the presentation of numerical simulations of our model in several scenarios, using the numerical parameters given above. In particular, we shall consider different values of the coupling cross-terms amplitudes to track their impact on the immiscible fluid displacement and the waterfront position. More precisely, the following scenarios will be investigated: first, we compare the solution of the coupled system (54) to the solution of the classical system (1) in the particular setting of kr g,w = kr w,g . Next, we consider equal amplitudes of coupling cross-terms, that is, M g,w = M w,g . In addition, we shall introduce a tolerance parameter denoted ε and defined as kr w,g = ε kr g,w allowing for relative variation of the coupling cross-terms, and we shall consider the particular values ε ∈ [1/2, 3/2, 5/2]. Eventually, the impact of variations of the Forchheimer coefficient on the solutions will be shown numerically. Moreover, we shall estimate the breakthrough time for all these scenarios as an application of our model.

Particular Coupling Cross-Terms
We consider the particular case of the coupled system (54) with equal coupling crossterms, which is kr g,w = kr w,g , where Our simulations show that, at the initial stage, there is no notable difference in the water saturation profile solution of the coupled Forchheimer-Darcy one (54) compared to the solution of the classical Forchheimer-Darcy system (1) (see Figure 3b), the difference becomes important at the front where the position of the waterfront is different by about 5.66% (see the zoom in Figure 3d). This observation is confirmed by the water fractional flow graph corresponding to both systems, where we see a clear difference in the value of the water saturation at the front, S w f (see Figure 3a). The graphs of the different velocities are given in Figure 3d, showing in particular that the total velocity remains absolutely constant, which is the case theoretically. The water and gas pressure graphs are given in Figure 4a,b, presenting an expected linear behavior of water-phase and gas-phase pressures after the front. The pressure before the front depends on the water saturation nature, which justifies its curved profile.
Next, we consider a second particular case, namely the case of equal coupling crossterms amplitudes. More precisely, we assume that M g,w = M w,g , thus kr w,g = µ g µ w kr g,w and we refer to (79) for detailed expressions (observe that this case is different from the case presented above since µ g = µ w ). We recall that the coupling terms are given by the expression (79). In this scenario, we observe that the impact of the coupling is rather negligible, as it can be seen in the graphs of the fractional water flow and the water saturation profiles (see Figure 5). Indeed, from Figure 5d, one can estimate that the waterfront positions are different by only 0.05%.  (c) water velocity v w , gas velocity v g , the total constant velocity v t ; and water-phase pressure and gas-phase pressure after 1000 days.

Different Coupling Cross-Terms' Amplitudes
In this scenario, we would like to emphasize the impact of a relative variation of the coupling cross-terms amplitudes that is when a notable difference of their amplitudes is chosen, on the waterfront position and the solution of system (54). For this purpose, we assume that the pseudo-permeabilities kr g,w and kr w,g are linked via a tolerance parameter ε that is kr w,g = ε kr g,w , and we shall vary the tolerance ε as follows: ε = {1/2, 3/2, 5/2} such that kr w,g = ε kr g,w , and we remind readers that kr g,w given by (79). Figure 6b, and particularly the zoom in Figure 6d, shows numerically that the larger the tolerance parameter, the more advanced the waterfront position in the reservoir. It can be computed that the waterfront position is in advance of about 4.95%, 5.35% and 5.73% after 1000 days for tolerance ε values 1/2, 3/2 and 5/2, respectively.

Variation of the Forchheimer Coefficient
The aim of this section is to numerically investigate the behavior of the solution of the coupled system with respect to variation of the Forchheimer coefficient. From the theoretical point of view, we believe one can exhibit continuous dependence of the solution on the Forchheimer coefficient. Thus, we vary the values of the coefficient C g and plot the associated fractional flow and the water saturation profiles describing the two-phase flow displacement in the reservoir. We shall consider the following values (we refer to [10] for more details): C g = [3.23 × 10 −6 , 3.43 × 10 −6 , 3.54 × 10 −6 ]. As an application of the previous results, one can compute the breakthrough time associated with the previous considered scenarios using the following formula and compare them to the classical system: For the first scenario, the breakthrough in terms of days: 1.
For the first scenario of kr g,w = kr w,g , we obtain Nature of the system Breakthrough t L Classical system t L = 1224.5 days Coupled system t L = 1163.8 days

2.
For the second scenario M g,w = M w,g Nature of the system Breakthrough t L Classical system t L = 1224.5 days Coupled system t L = 1166.2 days 3.
For the third scenario, where variations of the Forchheimer coefficient are considered: Tolerance ε Breakthrough t L ε = 0.5 t L = 1164.5 days ε = 1.5 t L = 1169.1 days ε = 2.5 t L = 1178.5 days

Post-Processing: Solution of Coupled Forchheimer-Darcy's System via a Decoupled System
In this section, we develop a post-processing algorithm for solving the physical coupled Forchheimer-Darcy system (54) using a decoupled non-physical system since the idea is based on the introduction of modified mobilities, that is, by solving first a classical Forchheimer-Darcy system (1) with modified mobilities and next post-processing the solution to obtain the solution of the coupled system (54). The advantage of this postprocessing algorithm is that it allows for solving coupled systems (and therefore taking into account the phases coupling) using any existing code or industrial software designed to solve only classical (that is not coupled, i.e., (1)) without any technical modification. More precisely, with the post-processing algorithm, there is no need to change the software's source code that solves classical systems, but one only solves a classical system with the modified mobilities and then post-processes the solution to take into account the coupling. The modified mobilities are given by 1.
For the water phase,M 2.
For the gas phase, we define the following termsM g,w andM g,w such as In addition, we introduce the pressuresp g andp w related to the pressure p g and p w through the following expressions: Let us mention that, given the numerical values of the mobilities, modified mobilities, and the capillary pressure,p g andp w can be obtained by integrating the previous expressions with respect to x. In addition, observe that the new capillary pressure is now related to the original capillary pressure as follows: g,wMg,w +M w,gMw,g 2M g,wMw,g ∂ x p c,g,w .
Regarding the post-processing approach, we have the following: Lemma 6. Let us assume M w,g = M g,w , and let (v g ,p g ) and (v w ,p w ) be solutions of the modified system: v g + k β g ρ g M g v 2 g = − kM g,w ∂ xpg , v w = − kM w,g ∂ xpw . (86) Then, (v g , p g ) and (v w , p w ) are solutions of system (54).

Proof.
Using the pressures' definition (85), system (86) reads v g + k β g ρ g M g |v g |v g = − k M g ∂ x p g − kM g,w ∂ x p g + k M w,g + M g,w 2 ∂ x p c,g,w , Using the capillary pressure definition, we infer Now, using the fact that M g,w = M w,g , we infer that (v g , p g ) and (v w , p w ) are solutions of the coupled Forchheimer-Darcy system (54).
Let us mention that Lemma 6 shows the equivalence of both systems in the very specific case of M g,w = M w,g , which is clearly restrictive in terms of applications. This is obviously due to the fact that we were not able to prove theoretically the equivalence of both systems theoretically in more general cases. However, we believe that a more involved post-processing approach might be possible for the general case, at least if applied at the level of the numerical scheme and not the mathematical systems.

Conclusions
In this paper, we introduced a mathematical model for two-phase non-Darcian flow where the phases interaction is modeled through coupling pressure cross-terms depending on pseudo-permeabilities and the fluids' viscosities, namely a coupled Darcy-Forchheimer model. The model being new to the literature (the novelty resides in the modeling of the coupling), a rigorous mathematical analysis of the coupled system of differential equations was provided leading to the existence and uniqueness of solutions. In addition, we developed the associated Buckley-Leverett theory by providing explicit expression of the phases' velocities and the fractional flows depending only on the gradient of the capillary pressure and the water saturation. We provided numerical simulations solving the front advance Equation (66) using a finite volume method. Several scenarios were considered: first, we considered the case of equal pseudo permeabilities, kr g,w = kr w,g , where we observed a clear difference in the water saturation profiles compared to the classical (non coupled) system; the difference is about 5.66%. In the second scenario, we set equal amplitudes for the coupling pressure cross-terms, that is, we set M g,w = M w,g . In this situation, we showed numerically that the impact of the coupling is rather negligible; indeed, the difference of the waterfront positions of the coupled system and the non coupled one differ only by 0.05%. In addition, we analyzed numerically the change in the front position when a relative variation of the coupling cross-terms amplitudes is introduced. We observed a higher tolerance (and therefore difference between the amplitudes), the more advanced the waterfront position is. Eventually, we provided numerical illustrations of the impact of variation of the Forchheimer coefficient and observed a clear impact on the waterfront position. Eventually, as application of our model, we computed the breakthrough time for all previous scenarios. In particular, we observe a rather important gain compared to the classical (non coupled) system. Eventually, a post-processing algorithm for solving the coupled Forchheimer-Darcy system was introduced. The idea is to use existing codes and industrial software designed for solving non coupled systems without any technical change to solve our model, and therefore take into account the phase coupling we proposed. Our method is based on the introduction of modified mobilities. Funding: The work of Mostafa Kadiri was funded by the TAMUQ start-up fund of Saber Trabelsi No.482202-50120. The authors would like to acknowledge Qatar Foundation and Qatar National Library for their financial support.

Conflicts of Interest:
The authors declare no conflict of interest.