COMPRESSIBLE AND VISCOUS TWO-PHASE FLOW IN POROUS MEDIA BASED ON MIXTURE THEORY FORMULATION

. The purpose of this work is to carry out investigations of a generalized two-phase model for porous media ﬂow. The momentum balance equations account for ﬂuid-rock resistance forces as well as ﬂuid-ﬂuid drag force eﬀects, in addition, to internal viscosity through a Brinkmann type viscous term. We carry out detailed investigations of a one-dimensional version of the general model. Various a priori estimates are derived that give rise to an existence result. More precisely, we rely on the energy method and use compressibility in combination with the structure of the viscous term to obtain H 1 -estimates as well upper and lower uniform bounds of mass variables. These a priori es- timates imply existence of solutions in a suitable functional space for a global time T > 0. We also derive discrete schemes both for the incompressible and compressible case to explore the role of the viscosity term (Brinkmann type) as well as the incompressible versus the compressible case. We demonstrate similarities and diﬀerences between a formulation that is based, respectively, on interstitial velocity and Darcy velocity in the viscous term. The investiga- tions may suggest that interstitial velocity seems more natural to use in the formulation of momentum balance than Darcy velocity.


(Communicated by Paola Goatin)
Abstract. The purpose of this work is to carry out investigations of a generalized two-phase model for porous media flow. The momentum balance equations account for fluid-rock resistance forces as well as fluid-fluid drag force effects, in addition, to internal viscosity through a Brinkmann type viscous term. We carry out detailed investigations of a one-dimensional version of the general model. Various a priori estimates are derived that give rise to an existence result. More precisely, we rely on the energy method and use compressibility in combination with the structure of the viscous term to obtain H 1 -estimates as well upper and lower uniform bounds of mass variables. These a priori estimates imply existence of solutions in a suitable functional space for a global time T > 0. We also derive discrete schemes both for the incompressible and compressible case to explore the role of the viscosity term (Brinkmann type) as well as the incompressible versus the compressible case. We demonstrate similarities and differences between a formulation that is based, respectively, on interstitial velocity and Darcy velocity in the viscous term. The investigations may suggest that interstitial velocity seems more natural to use in the formulation of momentum balance than Darcy velocity.

1.
Introduction. The importance of multiphase flow in porous media has long been recognized in many fields. Mathematical modelling of multiphase flow is essential in practical applications like enhanced oil recovery and geological CO 2 storage in depleted oil and gas reservoirs [25,42] as well as biological processes [29,14,18,39,40,37]. Traditional formulations of multiphase flow describe macroscopic fluid fluxes with a straightforward extension, first proposed by Muskat [30,3], of Darcy's equation for single-phase flow. Unlike in the single-phase case, this extension cannot be rigorously obtained from first principles [22,23]. The multiphase extension of Darcy's equation may be described as a quasi-linear relation, because the fluid flux depends linearly on the "driving force", which includes viscous, capillary, and gravity forces, and all the nonlinearity is agglutinated in the relative permeability and capillary pressure functions [25]. An instructive overview is given in [32] of how generalizations of the standard Darcy's law for single phase flow can be derived within the context of mixture theory. Starting with more general momentum balance equations and using different sets of assumptions leads to a hierarchy of mathematical models. In particular, it can be shown that popular models due to Brinkman, Biot and many others can be obtained via various approximations. Interesting extensions of the classical multiphase formulation are also discussed by Wu [42].
1.1. A compressible and viscous two-fluid model for porous media flow. The model we are interested in describes flow of two compressible immiscible fluids, e.g., water (w), oil (o), or gas (g), moving in a porous media and takes the following form (we use "w" and "o" in the following as index): with capillary pressure P c defined as the difference between the non-wetting fluid (oil) pressure P o and wetting fluid (water) pressure P w P c = P o − P w = P c (s w ), P c (s w ) < 0. (1.2) Herein, φ is the porosity of the medium, ρ i represents density and s i the volume fraction (saturation) where i = w, o. In addition, we have the fundamental relation that expresses that the water and oil occupy the pore space Furthermore, ε w , ε o (assumed to be constant) characterize the magnitude of the viscous terms. The model can be derived (or at least motivated) from general mixture theory [11,32] where we study a continuum composed of matrix occupying a volume fraction (1 − φ) and a pore space of volume φ that is filled with a mixture of water and oil represented, respectively, by φs w and φs o such that (1 − φ) + φs w + φs o = 1. The matrix is stagnant whereas the two fluids move with (locally) averaged interstitial velocities u w and u o . We refer to the recent work [31] for more details leading to (1.1). See also [1,29,34] and references therein for interesting examples of similar models developed in the mixture theory framework. Note that the viscous terms ε o ∇ · (n∇u o ) and ε w ∇ · (m∇u w ) (Brinkman type of term) included in (1.1) 3,4 involve a mass dependent coefficient whose magnitude is governed by the parameter ε i . This reflects that we have introduced kinematic viscosity ε that is related to dynamic viscosity µ by ερ = µ for single-phase flow of a fluid with density ρ [28]. Combined with the two-phase momentum balance for water and oil this gives rise to mass dependent viscosity coefficients of the form ε o n and ε w m. We refer to [14] (and references) therein for more details. More generally, we may think of ε o n and ε w m as "effective" viscosities since the model (1.1) must be understood as the result of a volume averaging process where variables have been obtained through averaging over a small representative volume element, implying that detailed information about complex interfaces between the two phases have been lost and are represented only in an averaged sense [11,32]. Some authors also denote this viscosity as the "Brinkman viscosity" whereas the viscosity associated with the rock-fluid friction termk i (i = o, w) is denoted the "Darcy viscosity" [27]. This issue is also discussed in [38] where it is observed by means of an up-scaling procedure based on volume averaging methods, that the use of a slip boundary condition gives rise to an effective viscosity different from the one corresponding to the fluid phase.
In the following we will focus on nonlinear coupling mechanisms and we therefore assume physical parameters like porosity φ, absolute permeability K, Darcy viscosity µ i (will be introduced later), and Brinkman viscosity ε i to be constant. Generally,k o ,k w , andk ow depend on the fluid composition through s i and ρ i .
1.1.1. Closure relations. The above model must be endowed with appropriate closure relations for densities ρ i = ρ i (P i ). The two phases will be treated as weakly compressible fluids. More precisely, we represent the water and the oil by linear pressure-density relations of the form where C w and C o reflect the compressibility of water and oil, respectively. An essential role is played by the interaction coefficientsk ow ,k w , andk o . We will come back with more details about the choice of these. In addition, a functional form of the capillary pressure P c (s w ) must also be specified. Combining where ν is the outward normal on ∂Ω. The corresponding initial data is n(x, t = 0) = n 0 (x), m(x, t = 0) = m 0 (x), x ∈ Ω. (1.6) 1.2. The model (1.1) as a generalization of Darcy's equation based formulation. We may ignore the effects from the viscous terms in (1.1) 3,4 by setting ε o = ε w = 0. In addition, we neglect the fluid-fluid interaction effect by settinĝ k ow = 0, combined with the assumption that fluid-pore resistance force coefficient k i takes the formk where K is the absolute permeability (assumed here to be a scalar, i.e., we assume a homogeneous media), k ri is relative permeability, and µ i viscosity. This gives from (1.1) 3,4 the following reduced momentum equations for i = w, o which is nothing but the classical Darcy law where U i is the Darcy velocity. When combined with (1.1) 1,2 we arrive at the classical two-phase formulation [42]. The model (1.1) involves two main extensions from classical formulation based on two-phase Darcy's equation, as elaborated upon in the following: • The interaction forces on the RHS of (1.1) 3,4 involve a fluid-fluid drag force effectk ow (u w − u o ) in addition to fluid-rock drag force −k o u o and −k w u w . The only interaction force in Darcy's equation is the latter one representing friction between fluid and boundaries of the pores [32]. Moreover, while the drag force depends on the velocity, it is by no means necessary that it in general depends linearly on the relative velocity. See for example [42] for extensions that include nonlinear dependence on fluid velocity (Forchheimer). Moreover, it has been observed that inclusion of the fluid-fluid interaction termk ow (u w − u o ) can give improvements over standard Darcy's equation based formulation for water-oil flow in porous media. We refer to [35,31] for a first discussion of this in the context of core scale modelling and generalized permeability functions and [36] for a discussion of this generalized two-phase flow in the context of imbibition (i.e., capillary pressure driven counter-current flow). • The viscous terms ε o ∇ · (n∇u o ) and ε w ∇ · (m∇u w ) in (1.1) 3,4 can account for frictional forces within the fluid due to its viscosity. Ignoring these terms essentially imply that the viscosity of the fluid and the roughness of the solid surface lead to far greater frictional resistance (and hence dissipation) at the porous boundaries of the solid in comparison to the frictional resistance in the fluid [32]. Note that (1.1) 3,4 can naturally be interpreted as a two-phase version of Brinkman's equation. Brinkman's equation amounts to using a momentum balance equation that takes the following form for a fluid-matrix system [32]: ∇p + ρg = −α 0 u + µ∇ · (∇u), (1.9) where the quantities refer to the fluid phase, i.e., p is fluid pressure and u is pore velocity, ρ phase density and µ phase viscosity. g denotes the external gravity force and α 0 the magnitude of the fluid-matrix drag force. The final term, involving second order derivatives of velocity, marks a clear distinction from Darcy's law for single phase flow. This may be relevant for porous media having large permeability and/or being dominated by a network of fractures [27]. Some more precise remarks seem useful in order to set the model (1.1) into a broader context. Remark 1.1. We may replace the viscous term in (1.1) 3,4 that accounts for the fluid viscous shear effects that oppose the flow through the porous structure by a more general term ∇ · (ε i,eff ∇u i ), where the effective viscosity coefficient ε i,eff depends on other variables than the mass. We refer to [38] for a discussion of this, both from theoretical and numerical investigations. For example, from physical considerations and experimental investigations it seems clear that it could depend on pressure [32]. It is concluded that it might be reasonable to include dependence on pressure both in modelling of fluid-pore friction force as well as frictional effects within the fluid itself [32]. More generally, one should also account for the possibility that the flow may not be steady implying that one should add a term ρu t to (1.9) in the single-phase case or (φs i ρ i u i ) t in the more general two-phase model (1.1). Moreover, nonlinear inertial effects in the fluid cannot be ignored if the flow is not sufficiently slow. Remark 1.2. An interesting study of a two-phase Brinkman type of model is found in [9]. A two-phase formulation based on Darcy's equations of the following form is studied with s = s w is total mobility and f (s) = λ w (s)/λ T (s) is the fractional flow function for water phase. U T is the total Darcy velocity U T = U w + U o . The following two-phase version of Brinkman's equation based on Darcy's velocity U i is proposed as a generalization of Darcy's equation: Summing the two momentum equations in (1.11) gives rise to a Brinkman type of momentum equation for the mixture of the two phases expressed in terms of the total Darcy velocity U T = U w + U o : Taking the divergence of (1.12) gives in light of the total mass balance equation (1.10) 2 the following Brinkman based approximation of (1.10)  (1.11). In [41] it is argued that the most natural choice, at least for the flow system they consider with creeping flow inside moving permeable particles, is to use interstitial (intrinsic) phase velocity u i in the macroscopic equations. Their conclusion is based on numerical computations and comparison of the model based on, respectively, phase velocity U i = s i u i and interstitial phase velocity u i .

1.3.
Purpose of this work and brief review of related works. The aim of this paper is three-fold: (i) Present an example of stability analysis motivated by the study of compressible Navier-Stokes equation (and different from traditional two-phase porous media stability analysis) which exploits the structure of the viscous term in the momentum balance and accounts for compressibility; (ii) Present an example of a numerical scheme both for the compressible and incompressible version of (1.1) in a one-dimensional setting and demonstrate similarities and differences through some specific simulations; (iii) Gain some insight into the impact the viscosity terms have on the solution. We end this section by giving a brief review of other works on the two-phase porous media model based on Darcy's law. Most works focus on the incompressible, immiscible two-phase flow case. For example, [2] studied the existence of weak solutions for the incompressible two-phase model in fractured porous media based on a dual-porosity formulation. Regularity and stability results were obtained in [8] when analysing a coupled system involving a saturation and a global pressure. In [6,5] the authors showed existence of a solution for an incompressible two-phase flow within a heterogeneous porous medium made of two rock types. Considering dynamic capillary pressure [24] for the incompressible two-phase flow in porous media, [7] proved the existence of a weak solution to a degenerate elliptic parabolic system whereas in [13,12] existence conditions for the traveling wave solution were derived. In particular, non-monotone weak solutions for the Buckley-Leverett equation were obtained. Interesting contributions have also been made concerning the compressible immiscible two-phase flow in porous media where phase densities are assumed to depend on their own pressure. Without using the feature of global pressure, Khalil and Saad [26] established an existence result for a three-dimensional model. In addition, the implicit finite volume scheme was studied in [33] to obtain convergence to a weak solution.
2. Stability analysis and existence of solution in the one-dimensional setting. The purpose of this section is to derive a priori estimates of the solution of (1.1). The approach is quite different from the approach used for the incompressible model and formulation based on Darcy velocity U i where the first step is to derive estimate for pressure [9]. It is also different from mathematical analysis of compressible two-phase models that are based on global pressure [19,20,21]. We rely on the energy method where we first derive an energy-type of estimate. In addition, the special structure of the viscous terms allows one to obtain estimate of m, n in H 1 along the lines of the Bresch-Desjardin method [4] for a two-phase Navier-Stokes model. For analysis of related models we refer the interesting reader to [15,16,17,18], and references therein. In the following we consider the 1D version of (1.1) where source terms have been set such that water is injected and possibly produced whereas oil is produced only, i.e., Q o = −nQ p whereas Q w = −mQ p + ρ w Q I for constant Q p , Q I . The model takes the following form with (n, m, u w , u o ) as the main variables (2.14) subject to the boundary condition and initial condition Note that the gravity constant g can take both signs depending on the orientation of the x-coordinate axis. Above we assume that positive direction of x-axis points downward and g > 0. We also use the notationk =k ow in (2.14), for simplicity reason only.
Definition of (ρ o , ρ w , s o , s w ). Let us first see how we can obtain ρ w and ρ o as a function of masses m and n. We focus on the situation where m, n > 0. The cases where m = 0 or n = 0 are treated separately. We rewrite s o + s w = 1 as On the other hand, from (1.2) we have (2.19) where we have introduced the function F (ρ w ; m, n) where m, n are thought of as parameters. Clearly, for any choice of m, n > 0, we want to verify that F (u; m, n) (where we use u as the main variable) has a unique zero point which we denote as ρ w (m, n). Let us check some basic properties of F (u; m, n) as a function of u.
By the definition of m, it is natural to look for ρ w which belongs to (m, +∞). Moreover, from (2.19) we observe that (i) F (u → m + ; m, n) = +∞; (ii) F (u → +∞; m, n) = −∞. Next, we check monotonicity properties of F (u; m, n) as a function of u.
Since F u (u, m, n) < 0 in (m, +∞) for any given m, n > 0, and F : (m, +∞) → (−∞, +∞) as observed above, it follows (from the intermediate value theorem) that there is a unique ρ w = ρ w (m, n) ∈ (m, +∞) such that F (ρ w ; m, n) = 0. In addition, since F u (ρ w ; m, n) = 0, it concludes that the function ρ w is differentiable with respect to m or n (from the implicit function theorem). Furthermore, ρ o , s o , and s w are then obtained as follows: For the limit case when m = 0, there are two options: (i) s w = 0, which implies that ρ o = n and ρ w is found from (2.18); (ii) s w > 0, which implies that ρ w = 0 and where ρ o is found from (2.18). Similarly, we can compute ρ w and ρ o for the case n = 0.
Notation. We first give some notation. Assumptions. The following assumptions are made: • Capillary pressure P c (s w ): We assume that for Φ(s w ) such that Φ (s w ) = P c (s w ), the following property holds: Φ(s w ) ≤ P c (s w )s w , 0 ≤ s w ≤ 1 (2.23) wheres w = s w (m,ñ) andm andñ refer to the total masses given by (2.22), which are constant due to Remark 2.1. Moreover, we assume that sup sw∈(0,1) and that Note that these constraints on the capillary pressure P c (s w ) are all mild and physical reasonable conditions. • Source terms in (2.14) 1,2 are ignored by setting Q p = 0 = Q I . • Interaction termk w ,k o , andk are set as follows: Remark 2.1. Clearly, in view of (2.14) 1,2 , the condition (2.15), and assumption Q p = Q I = 0, it follows from (2.22) that wherem 0 ,ñ 0 are constant.

Remark 2.2.
As far as the condition on capillary pressure P c (s w ) as given by (2.23) is concerned, we may observe that this appears to be a weak structural constraint. Consider for example a capillary pressure curve of the form P c (s w ) = −P * c ln(δ+ sw a ), for some δ, a > 0 as a typical example of a physical relevant function. Clearly, from the relation Φ (s w ) = P c (s w ) we can introduce two positive constants C 1 and C 2 to be determined such that Since x ln(x) is an increasing function for x ≥ e −1 whereas for x ∈ [0, e −1 ) decreases from zero for x = 0 and takes a minimum −e −1 , it is clear that we can secure that for an appropriate choice of C 1 such that we conclude from (2.28) that Φ(s w ) ≤ P * c s w − C 2 . What remains to show then is that Clearly, (P * c − P c (s w ))s w ≤ C 2 for an appropriate choice of the constant C 2 = C 2 (P * c ,s w ) since s w ∈ [0, 1]. 2.1. Main results. First, we present a local (in time) existence result whose proof is presented in Appendix A. Then, we state an (almost) global in time existence result which relies on the local existence result combined with certain a priori estimates, see (2.29). Section 2.2 is devoted to these estimates.
where I wo refers to the coefficient in (2.26), k 0 = min inf n0 e , inf m0 e and k 1 = max e sup m 0 , e sup n 0 , and where I w , I o are coefficients given by (2.26) and C is a positive constant depending on k 0 ,k 1 and some other known data but independent of ε o and ε w (see Step 2 in Appendix A for more details). Then there exists a positive constant T 0 , such that the system (2.14) with initial-boundary conditions (2.15) and (2.16) has a unique solution (m, n, Now we are in the position to state our second result on the almost global existence.
Theorem 2.2. (Almost global existence) In addition to the assumptions of Theorem 2.1, for any given T > 0, if K 1 < min ε wm , ε oñ , then the system (2.14) with initial-boundary conditions (2.15) and (2.16) has a unique solution (m, n, Moreover, we have the following estimates: Remark 2.3. The constraint K 1 < min ε wm , ε oñ where K 1 is given by (2.45), implies smallness of initial data combined with assumption about sufficiently large viscosity ε w and ε o . More precisely, from the definition of K 1 for a fixed T > 0 we may choose m 0 and n 0 such that Note that this constraint is only used to get the positive lower bound of m and n. See Corollary 2.3 for more details. Hence, the obtained estimates cannot be used to investigate the limit when ε w , ε o → 0.

2.2.
Proof of Theorem 2.2. Equipped with Theorem 2.1, we are going to prove Theorem 2.2. More precisely, let T * denote the maximum time for the existence of solutions as in Theorem 2.1 1 . Then Theorem 2.1 implies that T * > 0. To prove the almost global existence, it suffices to show that T * is larger than the given T which can be taken as large as possible. For otherwise, i.e., T * ≤ T , it will lead to a contradiction based on the following estimates uniformly for t, i.e., In fact, (2.29) implies that T * is not the maximum time for the existence which is the desired contradiction.
To get (2.29), we need the following lemmas. To simplify the proof, we let C(T ) denote a generic positive constant depending on the initial data and T. Moreover, for any given T > 0, C(T ) < ∞. We let t < T * ≤ T throughout the rest of this section, i.e., in Lemma 2.2-Corollary 2.5. Note that C(T ) ≥ C(T * ) and (a) Energy estimate.
where E(t) is given by Proof. From the two momentum equations of (2.14) 3,4 we get after a multiplication, respectively, by u o and u w , followed by integration over [0, 1], integration by parts and use of (2.15) (2.32) For I 0 , integrating the two mass equations (2.14) 1,2 on (0, x) for any given x ∈ (0, 1), and using the boundary condition, we have Thus we have As to I 1 , we observe that by using (1.4) and (2.14) 1,2 we can compute as follows and, by similar calculations (2.35)

YANGYANG QIAO, HUANYAO WEN AND STEINAR EVJE
Note that here we have also used (2.27). Consequently, using that P w = P o − P c and (1.3), we find from summing (2.34) and (2.35) That is, Moreover, we see that for some reference densityρ w > 0. Hence, again by using (2.27) Thus, it follows that (2.39) Note that in view of (1.4), the last line of (2.39) gives us whereρ o ,ρ w , ands w are related to each other by common massesm,ñ, i.e., we have that Hence, it follows from (2.39) that (2.40) Inserting (2.40) and (2.33) in (2.32) we get We can rewrite (2.41) to be with E(t) as given by (2.31). Hence, we conclude that (2.30) holds and where it is also clear from (2.23) that Lemma 2.1 implies the following corollary.
where ρ i0 and s i0 refer to initial states.
Remark 2.4. The energy equality as expressed by (2.30) contains several dissipation terms on its left-hand-side. The three last terms reflect that there is a loss of energy (i.e., an energy transformation) through the three different friction terms, respectively, with coefficientsk,k w , andk o and is a consequence of the viscous property of the involved fluids leading to drag force effects. Similarly, the internal viscosity of each fluid also creates resistance to move and is reflected by the two viscous terms, respectively, with coefficients ε w and ε o .
Remark 2.5. The energy E(t) defined in (2.31) contains different terms each having a specific physical meaning. The use of compressible fluids implies that energy can be stored and released in the fluid due to pressure variations. An intuitive example of this is when there is influx of gas (oil) in a low reservoir layer where pressure is high. As this gas migrates towards a higher zone where pressure is lower, the gas (oil) will expand. That is, ρ o decreases according to (1.4) and s o increases since mass m o is conserved and m o = s o ρ o and therfore typically will displace the surrounding water phase represented by s w . This energy exchange is accounted for through the two first terms of E(t). Capillary pressure P c = P o − P w , accounts for the difference between the water and oil pressure P w and P o , and also acts as a driver (an additional pressure effect) for fluid motion. It naturally occurs in the energy functional E(t) similar to the gravitational energy, see the last line of (2.31).

45)
and a = max{1 Proof. Note that from (2.14) 2 we get the following reformulated equation after expanding the advective term and taking a derivative in space: Note the appearance of the viscosity term on the RHS of (2.46). Combining (2.46) with (2.14) 4 we arrive at This is the same as x m which clearly, by using (2.14) 2 , is the same as Now, we test (2.47) with w and combine it with (2.14) 2 and (2.15) which leads us to Similarly, for the oil phase we obtain Next, we focus on the terms appearing on the RHS of (2.48): We note that Similarly, for J o,1 associated with (2.49) (2.51) To conclude, we see that by summing (2.48) and (2.49), using (2.50), and (2.51), we get where we again have used P c (s w ) = P o − P w and s w + s o = 1. That is, For K ow0 , we use Cauchy inequality and the mass equations and have (2.54) In the following we make use of the functional form of the interaction coefficientŝ k w ,k o , andk as expressed in (2.26).
Proof. Differentiating (2.64) with respect to x, we have which implies that (2.68) By (2.62), we have Substituting this identity and that 1 − s w = s o into (2.68), we have we have Proof. By virtue of (2.32), we have where we use W 1,1 (0, 1) → L ∞ (0, 1) in the third inequality, and use ab dx ≤ C a 2 dx + ε b 2 dx with appropriate choice of ε in the last one. This implies Similarly, we get the estimate of (u o ) xx .  We consider (2.69) with ∂ x replaced by ∂ t , and use similar analysis as (2.72). Then With the above estimates, we get (2.29). Thus the proof of Theorem 2.2 is complete.
3. Numerical results. The main objective of this section is to carry out some testing of the numerical schemes presented in Appendix D, respectively, for the compressible and incompressible version of (1.1). First, we want to test general stability properties. Second, we seek some insight into the role played by using Darcy velocity U i = s i u i versus interstitial fluid velocity u i (i = w, o) in the viscous terms. In other words, we modify (4.145) 3,4 and use the following momentum equations for the compressible case and modify (4.165) 3,4 as follows for the incompressible case Third, we also test the behavior of the scheme as the coefficients ε w , ε o are varied to see what kind of impact this term will have on the solution. This also allows us to get some understanding of whether the viscous model seems to converge to the inviscid version obtained by letting ε w , ε o → 0. We conduct numerical tests similar to those reported in [9]. We refer to Remark 1.2 for more details regarding the model they solve. Most importantly, the viscous term in their model depends on the Darcy velocity U w , U o . We apply the scheme for the incompressible model for these investigations, see Section 3.1. However, in Section 3.2 we also include examples where we use the scheme derived for the compressible model (see Appendix D) and do some comparison with the results from the incompressible model. The following input data are chosen for the numerical examples.
We choose parameters as specified in Table 1. In particular, when combined with the relations (4.132) it gives rise to a fractional flow function given by . The function is illustrated in Fig. 1 (left figure). The initial data of water saturation is set to be as in [9]: (3.84) A horizontal reservoir layer is considered in this case and porosity is also assumed to be 1. The whole test is a 10-day flooding process with a constant interstitial water injection rate at left boundary, Q = 8.  Table 1. Input parameters of reservoir and fluid properties used for for the below simulations. Note that P wL is the boundary pressure at left for the incompressible model whereas for the compressible model it represents the initial pressure distribution.  Table 1 (left figure) and initial water saturation (3.84) profile (right), both similar to that used in [9].
correlations (4.132). The corresponding initial water saturation profile is shown in Fig. 1 (right figure).

The incompressible model.
Case 1. First, we want to compare the numerical results obtained by using the scheme from Appendix D (incompressible variant) and compare with similar results presented in [9], which are based on the model (1.13). We also mimic their scheme by slightly modifying the scheme prescribed in Appendix D (incompressible variant) where s i u ix is replaced by U ix in the viscous term, as described by (3.83). Results are illustrated in Fig. 2. We show water saturation profiles after 10 days flooding with, respectively, ε w = ε o = 10 7 and ε w = ε o = 10 6 (upper row) and compare with the corresponding results from Coclite et al. [9] (lower row). Main observations are: The results of two corresponding cases with ε w = ε o = 10 7 and ε w = ε o = 10 6 after a dimensionless time, 0.65, produced by the numerical scheme described in [9] to solve the model (1.13). From these computations we see that the solution is sensitive to whether the interstitial velocity u i or the Darcy velocity U i appear in the viscous term. In particular, the use of Darcy velocity seems to generate considerably more oscillatory behavior behind the "water bank" formed at the front.
(i) A new "water bank" is formed behind the front of the water as a result of the viscous terms. This is a local effect restricted to the region right behind the water front where large gradients in velocity are present; (ii) Internal viscous forces slow down the transport effect, especially at the saturation front where velocity involves dramatic changes. Right behind the water bank, the model with Darcy velocity involved in the viscous term tends to develop oscillatory behavior; (iii) The solution shows a clear sensitivity to the magnitude of ε o , ε w (i.e., 10 7 versus 10 6 ) for the scheme based on Darcy velocity U i in the viscous terms. The scheme with viscous term based on interstitial velocity u i shows less sensitivity to this change in ε o , ε w . The classical Buckley-Leverett model solution (i.e., ε w = ε o = 0) is composed of a sharp water front followed by a rare-faction wave which is due to a viscous effect associated with resistance forces between fluid (water and oil) and walls of the pore space. The new water bank is a consequence of internal viscosity within the fluid felt at the region behind the water front. The difference between solution when viscous term is based on Darcy velocity U i versus solution when viscous term is based on interstitial velocity u i , can be naturally understood in light of the expansion Clearly, the viscous term based on U i = s i u i gives rise to an additional term that naturally can be linked to the observed difference between the two schemes used to produce solutions in Fig. 2. It should be noted that Brinkman equation was developed empirically for single phase flow and afterwards has been extended to the multiphase setting in a heuristic manner. As noted in Remark 1.3 there seems to be an ongoing discussion in the literature whether to base the viscous term on Darcy velocity or interstitial velocity. Finally, in Fig. 3 is shown the result for the two schemes as ε w , ε o get smaller. Both schemes seem to reflect convergence toward the solution of the inviscid model (i.e., ε w = ε o = 0) with a considerably faster convergence produced by the scheme with viscous term based on interstitial fluid u i . Fig. 4 we show simulation results when we vary the internal relation between ε w and ε o . It is intuitively understandable that oil viscous effects can have a strong impact on the constant water level right behind the water front. Apparently, the same change of magnitude of water viscous parameter ε w from 10 4 to 10 6 with a constant ε o has a dramatic effect when ε o is small (i.e., 10 4 ) whereas the effect is rather small when ε o is large (i.e., 10 6 ). We refer to Fig. 4 for simulation results.

Case 2. In
Case 3. Now, we move to another case which has a different initial condition but still with the same injection rate of water, 8.004m 3 /d, as interstitial velocity at left boundary. The initial water saturation is illustrated in Fig. 5. Numerical results are shown in Fig. 6 where we compare the scheme based on interstitial viscosity (right figure) with the simulation result reported in [9] (left figure). In particular, it seems that the numerical solution based on using Darcy velocity produces an unphysical water saturation value near the left boundary. The numerical solution illustrated in Fig. 6 (right) does not contain this "defect". In addition, apparently the solution converges to the classical Buckley-Leverett type result with a small ε such as 10 3 . This behavior seems different from the conclusion in [9].
3.2. The compressible model. Next, we use the numerical scheme described in Appendix D (compressible variant) to compute and illustrate the numerical behavior for the compressible model. Comparison is made with the cases shown in Fig. 2 for the incompressible model. The compressible fluids are assumed to be given by  the pressure-density relation (1.4). The initial water saturation is also the same as shown in Fig. 1. The water saturation s w at left boundary is constrained with 0.8 and the initial water pressure at left boundary is 10 6 Pa. The numerical behavior is shown in Fig. 7. The essential difference is a delay in the solution of the compressible model.
In order to shed light on the difference observed in Fig. 7, we explore the pressure profiles at various times (shown in Fig. 8). It is clear from these plots that water pressure keeps increasing, especially for the water displacing part. Water can be compressed in the compressible model therefore the water density will also increase, which leads to a larger viscous effect since density is included in the coefficient of the viscous term. Water will feel more resistance forces and it is more difficult to displace oil resulting in a delay effect.
Injection of water versus gas. Finally, a numerical example is shown with gas injection to displace oil in the compressible model instead of water. The parameters of gas are the same as for water, as described in Table 1, except using the pressuredensity relation: ρ g = P g /C g where C g = 10 5 . The left figure of Fig. 9 compares the results of gas injection and water injection. As expected, it is a much slower  process for gas to displace oil. This is a natural consequence of the fact that gas is much more compressible than water. The high gas pressure which results from the increased viscous effect will generate a strong compression of gas. It is also interesting to see that the fluctuation in the gas saturation profile becomes stronger as time elapses, which is not observed in the case of water injection. However, the elevated constant water level is almost the same in both cases (compare the left and the right subfigures in Fig. 9). • We have found that exploiting the fact that the viscous term depends on the interstitial fluid flow velocity u w , u o , we can derive stability estimates (energytype estimate and BD-estimate) that also naturally deal with the capillary pressure term P c (s w ). This approach seems strongly linked to the special structure of the viscous coefficients. • We formulated finite difference schemes, both for the incompressible and compressible version of the model. These schemes allow us to systematically gain some insight into the effect of compressibility as well as the effect from the viscous terms that account for the frictional resistance within the fluid. We also observe that by using Darcy velocity in the viscous term, the resulting scheme tends to give more oscillatory behavior similar to that reported in [9]. • In particular, when the viscous coefficients ε w , ε o become small enough, the numerical experiments carried out in a one-dimensional setting indicate that the approximate solution converges to the solution of the inviscid model. The stability estimates for the model based on interstitial fluid velocity, however, do depend on ε w , ε o and cannot be used to ensure convergence to the inviscid model. Hence, this remains an interesting open problem. mous reviewers that helped improving certain parts of a first version of this manuscript.
Appendix A: Proof of Theorem 2.1. We apply a method similar to the one used in our previous work [16] to prove the local existence and uniqueness. Hence, in order to make the proof more compact we will heavily refer to that paper for details and highlight what is different.
First, we consider the solution space Step 1. Construct an iteration sequence.
We define an iteration sequence to approximate (2.14) which takes the following form with the initial-boundary value conditions: Step 2. Boundedness of the sequence. Assume that u k−1 w , u k−1 o ∈ S. To prove u i w , u i o ∈ S for all i = 0, 1, 2, 3, ..., it suffices to prove u k w , u k o ∈ S. In fact, as a consequence of that u k−1 w , u k−1 o ∈ S, we have for t ∈ [0, T 0 ]. By virtue of (4.85) 1 , (4.85) 2 , and (4.86), we can find positive constants C, k 0 and k 1 independent of T 0 , A, and A 1 and where T 0 ≤ T 1 for some T 1 > 0 which reflects smallness on time, such that and . Moreover, (4.87) 2,3 can then be deduced from the regularity of (m, n) → ρ s (m, n) and of ρ s → P s (ρ s ) for s = w, o. See Step 2, Section 2.2 in [16] for more details.
Multiplying (4.85) 3 by u k o and u k oxx respectively, and integrating the results over (0, 1), we have and where we have used Cauchy inequality, (4.87), (4.88). Note that we obtain for any small δ > 0, by using Sobolev inequality, Hölder inequality and Cauchy inequality. By virtue of (4.87), (4.89), (4.90), and (4.91), we have where we take δ = k0 10C . Putting (4.92) into (4.93), we have In view of (4.92) and (4.94), we have εwk0 . In view of (4.96) and (4.97), we have which can be done by assuming for instance that ε o and ε w are large enough, we obtain from (4.95) and (4.98) where we choose that Consequently, we have u k w , u k o ∈ S for all k = 1, 2, 3, · · ·, if we assume that (4.99) is satisfied for A 1 given by (4.101) and for T 0 ≤ T 1 .
as k i → ∞. It follows from (4.103) and (4.87) that Step 4. Convergence of (u ki−1 w , u ki−1 o ). We are going to investigate the convergence of the neighbor sequence of (u ki w , u ki o ), i.e., (u ki−1 w , u ki−1 o ), in order to make sure that their limits are the same, since they both appear in the approximate system (4.85).
To begin with, we need the estimates of the difference between m k+1 (n k+1 ) and m k (n k ), since there is a connection between velocity and mass due to the momentum equation. Denotem k+1 = m k+1 − m k andn k+1 = n k+1 − n k . Then, from (4.85) 1,2 it follows and n k+1 t +n k+1 whereC is a generic positive constant depending only on the initial data, the upper bound of T 0 and other known constants but independent of k , A and T 0 . Here we have used the facts that m k x is bounded in L 2 and that u k w ∈ S, and the Poincaré inequality. Similarly, we have  This, together with Hölder inequality, (4.109) and the boundedness of P k ox and where k = 2, 3, 4, · · ·. Similarly, by virtue of (4.85) 4 , we have where k = 2, 3, 4, · · ·. Putting (4.110) into (4.112), we have which yields where Then, (4.114) together with Gronwall inequality yields for small T 0 and large ε o , ε w . We note that (4.116) and (4.115) imply that (4.118) Thus, (4.118) combined with (4.102) 1 implies that as k i → ∞.
Based on (4.102), (4.103), and (4.119), it can be verified that (m, n, u w , u o ) is a unique solution of (2.14)-(2.16). See [16] for more details. Thus the proof of Theorem 2.1 is complete.
Appendix B: Some issues related to the general model (1.1). The purpose of this section is to elaborate on some of the differences between the model (1.1) and other Brinkman-type of two-phase models like the one mentioned in Remark 1.2.
The incompressible, viscous version of the generalized two-phase model (1.1). It seems useful to impose some simplifying assumptions on the model (1.1) and derive a reduced version of it. We may impose the following assumptions: • incompressible fluid, i.e., ρ w and ρ o are constant.
• porosity φ is constant We may rescale the time such that the porosity disappears in the time derivative terms. The mass and momentum equations in (1.1) now take the form where we have used that P c = P o − P w . Note that in the rest of this paper we will usek ow =k. We can solve for u w and u o from the two momentum balance equations in (4.120) and find that (4.121) That is, we find for U i = φs i u i : where ∆ρ = ρ w − ρ o and with generalized mobility functionsλ i (i = w, o, T ) of the formλ (4.123) Note that the resistance force coefficientsk w (s w ),k o (s w ),k(s w ) typically are functions of water saturation s w . From (4.122) we find (4.124) Remark 4.1. A fundamental difference between the above expression (4.124) for the total superficial velocity U T for the mixture and the model described in Remark 1.2 and expressed by the mixture momentum balance equation (1.12), is the viscous terms. Taking the divergence of (4.124) combined with the sum of the two first equations of (4.120) and (1.3), we arrive at This gives an elliptic pressure equation for P w . However, due to the appearance of complex viscous terms, it seems not easy to obtain H 1 -estimate of P w needed for compactness of the non-conservative pressure term s w ∇P w in the momentum balance equation.
We may conclude that it does not seem straightforward to analyse the incompressible, viscous model (4.120) by relying on an approach similar to the one used in [9]. The main reason is the use of interstitial fluid velocity u i instead of Darcy velocity U i = φs i u i .
The incompressible, inviscid case ε w = ε o = 0. It is instructive to also derive the model where viscous terms are set to zero. Hence, in the following we consider the incompressible model (4.120) where we assume that viscous terms in the momentum equations are ignored by setting ε w = ε o = 0. We observe that (4.124) gives are the standard fractional flow functions which satisfy thatf w +f o = 1. Using (4.126) in (4.122) we get where ∆ρ = ρ w − ρ o andĥ(s w ) is defined bŷ From this, it also follows that Consequently, the model (4.120) has been reduced to solving the following conservation law where U w = U w (s w ) is given by (4.128).  (4.123), recovers the classical incompressible immiscible model except that we now also include for an additional water-oil drag force effect through the termk. Setting this term to zero reproduces exactly the classical formulation. To obtain closed expressions for the generalized mobility functionsλ i as well as forĥ(s w ), we must specify appropriate functional correlations for (i) the water-rock resistance forcek w ; (ii) the oil-rock resistance forcek o ; (iii) the water-oil drag force effectk. In the recent work [31] we used correlations of the following form: where µ i is fluid viscosity, K absolute permeability, I w , I o , I are parameters that characterize the strength of the resistance force (similar to the end points of relative permeability in classical formulation). The generalized mobility functions introduced above account for two different resistance forces, instead of only one, as in standard Darcy's equation-based formulation. Mobility functions that are measured experimentally are known to generally depend on the flow configuration. Two main flow regimes are present in the expression for U w given by (4.128) and U o given by (4.130): Co-current flow (i.e., flow of water and oil in the same direction) represented by the first component U Tfw and counter-current flow (i.e., flow in the opposite direction) represented by −ĥ(s w )∆ρg andĥ(s w )∇P c (s w ) (compare with (4.130)). The fact that the fluid-fluid interaction termk is explicitly accounted for and appears inĥ (4.129) and mobility functions (4.123) implies that a more accurate understanding of water-oil flow mechanisms can be sought [31].
Appendix C: A pressure evolution equation. From the two mass balance equations we get after multiplying the n mass balance with ρ w and the m mass balance with ρ o and summing the two resulting equations with Clearly, Using this in (4.134) we get (4.139) Similarly, we have for P w : Clearly, Consequently, Using this in (4.140) we get (4.144) Appendix D: Discrete schemes for the compressible/incompressible version of (1.1).
A semi-discrete scheme for the compressible model. We consider discrete schemes for both the compressible and incompressible version of (2.14). For that purpose we introduce a reformulation that brings the compressible model closer to the incompressible model. In particular, we solve explicitly only for the mass transport of m = s w ρ w whereas the mass n is computed implicitly. This will be in the spirit of the incompressible approach where we solve the mass balance equation for s w and computes s o = 1 − s w . Details are given below. The starting point is the model (2.14) with (n, m, u w , u o ) as the main (unknown) variables. We rewrite the model in the following equivalent form with (m, P w , u w , u o ) as the main variables: (4.146) We refer to Appendix C and (4.144) which gives the pressure evolution equation and {P w,j (t)} Nx j=1 associated with the nodes {x j } Nx j=1 whereas the approximate velocities {u w,j+1/2 } Nx j=0 and {u o,j+1/2 } Nx j=0 are associated with the cell interfaces {x j+1/2 } Nx j=0 . In the following we describe a semi-discrete version of (4.145).
Step 2: Computation of velocities and pressure. Next, we solve simultaneously for P k+1 w,j and u k+1 w,j+1/2 and u k+1 o,j+1/2 by considering the following algebraic system from which we also compute the updated oil mass n k+1 j and pressure P k+1 o,j via (4.147). If necessary, we may repeat step 2 to improve the accuracy before we proceed to next time level. Note that we can only determine P w up to a constant and a reference pressure P * at some point in the domain may be specified. An equivalent formulation of (4.   ] j+1/2 appearing in (4.172) are based on "old" velocities u k w,j+1/2 and u k o,j+1/2 .