Compressible and viscous two-phase flow in porous media based on mixture theory formulation

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 $

Herein, $ \phi $ is the porosity of the medium, $ \rho_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, $ \varepsilon_w, \varepsilon_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-\phi) $ and a pore space of volume $ \phi $ that is filled with a mixture of water and oil represented, respectively, by $ \phi s_w $ and $ \phi s_o $ such that $ (1-\phi) + \phi s_w + \phi s_o = 1 $. The matrix is stagnant whereas the two fluids move with (locally) averaged interstitial velocities $ \mathbf{u}_w $ and $ \mathbf{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 $ \varepsilon_o\nabla\cdot ( n \nabla \mathbf{u}_o) $ and $ \varepsilon_w \nabla\cdot ( m \nabla \mathbf{u}_w) $ (Brinkman type of term) included in (1.1)$ _{3, 4} $ involve a mass dependent coefficient whose magnitude is governed by the parameter $ \varepsilon_i $. This reflects that we have introduced kinematic viscosity $ \varepsilon $ that is related to dynamic viscosity $ \mu $ by $ \varepsilon\rho = \mu $ for single-phase flow of a fluid with density $ \rho $ [28]. Combined with the two-phase momentum balance for water and oil this gives rise to mass dependent viscosity coefficients of the form $ \varepsilon_o n $ and $ \varepsilon_w m $. We refer to [14] (and references) therein for more details. More generally, we may think of $ \varepsilon_on $ and $ \varepsilon_wm $ 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 term $ \hat{k}_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 $ \phi $, absolute permeability $ K $, Darcy viscosity $ \mu_i $ (will be introduced later), and Brinkman viscosity $ \varepsilon_i $ to be constant. Generally, $ \hat{k}_o $, $ \hat{k}_w $, and $ \hat{k}_{ow} $ depend on the fluid composition through $ s_i $ and $ \rho_i $.

1.1.1.   Closure relations

The above model must be endowed with appropriate closure relations for densities $ \rho_i = \rho_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 coefficients $ \hat{k}_{ow} $, $ \hat{k}_w $, and $ \hat{k}_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 (1.2), (1.3), and (1.4) it follows that $ \rho_w = \rho_w(m, n) $ and $ \rho_o = \rho(m, n) $ are well-defined as functions of $ m $ and $ n $ for $ m, n\geq 0 $, from which we also can compute $ s_w = s_w(m, n) $ and $ s_o = s_o(m, n) $, see (2.17)- (2.21) for details.

1.1.2.   Initial and boundary conditions

Boundary conditions are prescribed as no-flux conditions:

where $ \nu $ is the outward normal on $ \partial \Omega $. The corresponding initial data is

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 $ \varepsilon_o = \varepsilon_w = 0 $. In addition, we neglect the fluid-fluid interaction effect by setting $ \hat{k}_{ow} = 0 $, combined with the assumption that fluid-pore resistance force coefficient $ \hat{k}_i $ takes the form

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 $ \mu_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 $ \mathbf{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 effect $ \hat{k}_{ow}( \mathbf{u}_w- \mathbf{u}_o) $ in addition to fluid-rock drag force $ -\hat{k}_o \mathbf{u}_o $ and $ -\hat{k}_w \mathbf{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 term $ \hat{k}_{ow}( \mathbf{u}_w- \mathbf{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 $ \varepsilon_o\nabla\cdot ( n \nabla \mathbf{u}_o) $ and $ \varepsilon_w \nabla\cdot ( m \nabla \mathbf{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]:

where the quantities refer to the fluid phase, i.e., $ p $ is fluid pressure and $ \mathbf{u} $ is pore velocity, $ \rho $ phase density and $ \mu $ phase viscosity. $ \mathbf{g} $ denotes the external gravity force and $ \alpha_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

where the effective viscosity coefficient $ \varepsilon_{i, \text{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 $ \rho \mathbf{u}_t $ to (1.9) in the single-phase case or $ (\phi s_i \rho_i \mathbf{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 $

where $ \phi\equiv 1 $, $ P_w = P_o = P $ (zero capillary pressure), $ \lambda_T(s) = \lambda_w(s)+\lambda_o(s) $ is total mobility and $ f(s) = \lambda_w(s)/\lambda_T(s) $ is the fractional flow function for water phase. $ \mathbf{U}_T $ is the total Darcy velocity $ \mathbf{U}_T = \mathbf{U}_w+ \mathbf{U}_o $. The following two-phase version of Brinkman's equation based on Darcy's velocity $ \mathbf{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 $ \mathbf{U}_T = \mathbf{U}_w+ \mathbf{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)

In [9] a notion of weak solutions to the Brinkman model (1.13) is introduced and convergence of an iterative approximation as well as full numerical scheme is demonstrated. Numerical experiments show that the numerical approximation is quite sensitive to the choice of $ \mu $ and creates oscillatory behavior. Interestingly, numerical experiments also indicate that this solution may not converge to the solution of the model (1.10) corresponding to $ \mu = 0 $.

Remark 1.3. In the literature there seems to be an ongoing interesting discussion of various formulations of two-phase versions of porous media flow based on Brinkman's equation. See for example the work [41] for a discussion of this. Traditionally, the superficial phase velocity (Darcy velocity) $ \mathbf{U}_i $ has been used in two-phase versions of Brinkman's equation similar to (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 $ \mathbf{u}_i $ in the macroscopic equations. Their conclusion is based on numerical computations and comparison of the model based on, respectively, phase velocity $ \mathbf{U}_i = s_i \mathbf{u}_i $ and interstitial phase velocity $ \mathbf{u}_i $.

1.3.   Purpose of this work and brief review of related works

The aim of this paper is three-fold: (ⅰ) 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; (ⅱ) 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; (ⅲ) 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 $ \mathbf{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 + \rho_wQ_I $ for constant $ Q_p, Q_I $. The model takes the following form with $ (n, m, u_w, u_o) $ as the main variables

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 notation $ \hat{k} = \hat{k}_{ow} $ in (2.14), for simplicity reason only.

Definition of $ (\rho_o, \rho_w, s_o, s_w) $. Let us first see how we can obtain $ \rho_w $ and $ \rho_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

Combining (2.18) with (2.17), we get

where we have introduced the function $ F(\rho_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 $ \rho_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 $ \rho_w $ which belongs to $ (m, +\infty) $. Moreover, from (2.19) we observe that (i) $ F(u\rightarrow m^+;m, n) = +\infty $; (ii) $ F(u\rightarrow +\infty;m, n) = -\infty $. Next, we check monotonicity properties of $ F(u;m, n) $ as a function of $ u $.

Since $ F'_{u}(u, m, n)<0 $ in $ (m, +\infty) $ for any given $ m, n>0 $, and $ F:(m, +\infty)\mapsto(-\infty, +\infty) $ as observed above, it follows (from the intermediate value theorem) that there is a unique $ \rho_w = \rho_w(m, n)\in(m, +\infty) $ such that $ F(\rho_w;m, n) = 0 $. In addition, since $ F_u^\prime(\rho_w; m, n)\not = 0 $, it concludes that the function $ \rho_w $ is differentiable with respect to m or n (from the implicit function theorem). Furthermore, $ \rho_o $, $ s_o $, and $ s_w $ are then obtained as follows:

For the limit case when $ m = 0 $, there are two options: (ⅰ) $ s_w = 0 $, which implies that $ \rho_o = n $ and $ \rho_w $ is found from (2.18); (ⅱ) $ s_w>0 $, which implies that $ \rho_w = 0 $ and where $ \rho_o $ is found from (2.18). Similarly, we can compute $ \rho_w $ and $ \rho_o $ for the case $ n = 0 $.

Notation. We first give some notation.

● $ L^p = L^p([0, 1]) $ for $ p\in[1, \infty] $

● We define

Assumptions. The following assumptions are made:

● Capillary pressure $ P_c(s_w) $:

We assume that for $ \Phi(s_w) $ such that $ \Phi'(s_w) = P_c(s_w) $, the following property holds:

where $ \tilde{s}_w = s_w(\tilde{m}, \tilde{n}) $ and $ \tilde{m} $ and $ \tilde{n} $ refer to the total masses given by (2.22), which are constant due to Remark 2.1. Moreover, we assume that

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 term $ \hat{k}_w $, $ \hat{k}_o $, and $ \hat{k} $ 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

where $ \tilde{m}_0, \tilde{n}_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(\delta + \frac{s_w}{a}) $, for some $ \delta, a>0 $ as a typical example of a physical relevant function. Clearly, from the relation $ \Phi'(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\geq e^{-1} $ whereas for $ x\in[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

What remains to show then is that

Clearly, $ (P_c^* - P_c(\tilde{s}_w))s_w \leq C_2 $ for an appropriate choice of the constant $ C_2 = C_2(P_c^*, \tilde{s}_w) $ since $ s_w\in[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.

Theorem 2.1. (Local existence) Assume that $ m_0\in H^1 $, $ n_0\in H^1 $ and $ \inf\limits_{x\in[0, 1]}n_0>0 $, $ \inf\limits_{x\in[0, 1]}m_0>0 $, and that

where $ I_{wo} $ refers to the coefficient in (2.26), $ k_0 = \min\Big\{\frac{\inf n_0}{e}, \frac{\inf m_0}{e}\Big\} $ and $ k_1 = \max\Big\{e\sup m_0, e\sup n_0\Big\} $, 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 $ \varepsilon_o $ and $ \varepsilon_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, u_w, u_o) $ on $ [0, 1]\times[0, T_0] $ in the sense that

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\Big\{\varepsilon_w\tilde{m}, \varepsilon_o\tilde{n}\Big\} $, then the system (2.14) with initial-boundary conditions (2.15) and (2.16) has a unique solution $ (m, n, u_w, u_o) $ on $ [0, 1]\times[0, T] $ in the sense that

where $ K_1 $ is given by (2.45).

Moreover, we have the following estimates:

for any $ t\in[0, T] $.

Remark 2.3. The constraint $ K_1< \min\Big\{\varepsilon_w\tilde{m}, \varepsilon_o\tilde{n}\Big\} $ where $ K_1 $ is given by (2.45), implies smallness of initial data combined with assumption about sufficiently large viscosity $ \varepsilon_w $ and $ \varepsilon_o $. More precisely, from the definition of $ K_1 $ for a fixed $ T>0 $ we may choose $ m_0 $ and $ n_0 $ such that $ K_1\leq 2[g^2( \tilde{m}+\tilde{n} )T + K_0] \exp(T) $ where $ \varepsilon_w, \varepsilon_o $ are chosen sufficiently large to ensure that

(ⅰ) $ \bar{a}\max\{\frac{1}{ \varepsilon_w}, \frac{1}{ \varepsilon_o} \} \leq 1 $; (ⅱ) $ 2[g^2( \tilde{m}+\tilde{n} )T + K_0] \exp(T)< \min\Big\{\varepsilon_w\tilde{m}, \varepsilon_o\tilde{n}\Big\} $. Note that this constraint is only used to get the positive lower bound of $ m $ and $ n $. See Corollary 3 for more details. Hence, the obtained estimates cannot be used to investigate the limit when $ \varepsilon_w, \varepsilon_o\rightarrow 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¹. 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^*\le T $, it will lead to a contradiction based on the following estimates uniformly for $ t $, i.e.,

¹ It means that the solution exists on $ [0, T^*) $ but not on $ [0, T^*] $.

for any $ t\in[0, T^*) $, where $ Q_{T^*} = [0, 1]\times[0, T^*) $. 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)<\infty $. We let $ t<T^*\leq T $ throughout the rest of this section, i.e., in Lemma 2.2–Corollary 2.5. Note that $ C(T)\geq C(T^*) $ and $ K_1 = K_1(T)\ge K_1(T^*) $ in Lemma 2.2.

(a) Energy estimate.

Lemma 2.1. For any $ t\in[0, T^*) $, it holds that

where $ E(t) $ is given by

where $ \tilde{\rho}_w = \rho_w(\tilde{m}, \tilde{n}), \tilde{\rho}_o = \rho_o(\tilde{m}, \tilde{n}) $, and $ \tilde{s}_w = s_w(\tilde{m}, \tilde{n}) $.

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)

For $ I_0 $, integrating the two mass equations (2.14)$ _{1, 2} $ on $ (0, x) $ for any given $ x\in(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

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 $ \tilde{\rho}_w>0 $. Hence, again by using (2.27)

and

Thus, it follows that

Note that in view of (1.4), the last line of (2.39) gives us

where $ \tilde{\rho}_o $, $ \tilde{\rho}_w $, and $ \tilde{s}_w $ are related to each other by common masses $ \tilde{m}, \tilde{n} $, i.e., we have that

(ⅰ) $ \tilde{\rho}_o = \rho_o(\tilde{m}, \tilde{n}) $, $ \tilde{\rho}_w = \rho_w(\tilde{m}, \tilde{n}) $, and $ \tilde{s}_w = s_w(\tilde{m}, \tilde{n}) $;

(ⅱ) $ P_o(\tilde{\rho}_o) = P_w(\tilde{\rho}_w)+P_c(\tilde{s}_w) $.

Hence, it follows from (2.39) that

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 $ \int_0^1 [P_c(\tilde{s}_w)s_w - \Phi(s_w)]\, dx\geq 0. $

Lemma 2.1 implies the following corollary.

Corollary 2.1. For any $ t\in[0, T^*) $, it holds that

where

where $ \rho_{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 coefficients $ \hat{k} $, $ \hat{k}_w $, and $ \hat{k}_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 $ \varepsilon_w $ and $ \varepsilon_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, $ \rho_o $ decreases according to (1.4) and $ s_o $ increases since mass $ m_o $ is conserved and $ m_o = s_o\rho_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).

(b) More regularity estimates.

Lemma 2.2. For any $ t\in[0, T^*) $, it holds that

and

where

and $ \overline{a} = \max\{1+I_{wo}+I_w, 1+I_{wo}+I_o\} $.

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

or

for

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

with

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)

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

In the following we make use of the functional form of the interaction coefficients $ \hat{k}_w $, $ \hat{k}_o $, and $ \hat{k} $ as expressed in (2.26).

and

Furthermore,

and

Putting the estimates (2.54)–(2.58) into (2.53), integrating the result over $ (0, t) $, and using Corollary 1, we have

where $ \overline{a} = \max\{1+I_{wo}+I_w, 1+I_{wo}+I_o\} $.

Using Gronwall's inequality and (2.59), we get (2.43). (2.44) is given by (2.43) and (2.59).

In view of Lemma 2.2 it follows that $ \sqrt{m}, \sqrt{n}\in H^1(0, 1) $. Combined with the Sobolev inequality $ H^1(0, 1)\hookrightarrow C([0, 1]) $ we have the following corollary.

Corollary 2.2. It holds that

for any $ (x, t)\in [0, 1]\times[0, T^*) $.

(c) Upper and lower bounds of density and related quantities.

Lemma 2.3. It holds that

for any $ (x, t)\in [0, 1]\times[0, T^*) $.

Proof. Since

we have

Armed with the upper bounds of $ m $ and $ n $ from Corollary 2.2, we get the upper bounds of $ \rho_o $ and $ \rho_w $ from (2.63) and (2.64). Note that the term $ \rho_i $ appearing on the right-hand-side is grouped with the corresponding $ s_i $ which gives either $ m $ or $ n $. In addition, we make use of the uniform bound on $ P_c(s_w) $ given by (2.24). The lower bounds can be derived from the definitions $ n = s_o\rho_o $ and $ m = s_w\rho_w $ combined with (1.3).

Corollary 2.3. The following uniform lower bound holds

subject to the constraint that $ K_1<\min\{ \varepsilon_w\tilde{m}, \varepsilon_o \tilde{n}\} $ where $ K_1 $ is given by (2.45).

Proof.

Letting $ \left(\frac{\tilde{m}K_1}{\varepsilon_w}\right)^\frac{1}{2}<\tilde{m} $, i.e., $ K_1< \varepsilon_w\tilde{m} $, we have

on $ [0, 1]\times[0, T^*) $ for some positive constant $ C $. Similarly, letting $ \left(\frac{\tilde{n}K_1}{\varepsilon_o}\right)^\frac{1}{2}<\tilde{n} $, i.e., $ K_1< \varepsilon_o\tilde{n} $, we get the positive lower bound of $ n $.

(d) Higher-order estimates.

Lemma 2.4. The following estimate holds

for any $ t\in[0, T^*) $.

Proof. Differentiating (2.64) with respect to $ x $, we have

where we have used $ s_o+s_w = 1 $ such that $ (s_o)_x = -(s_w)_x $.

With (2.67), we proceed to estimate $ (s_w)_x $.

which implies that

By (2.62), we have

Substituting this identity and that $ 1-s_w = s_o $ into (2.68), we have

Since

we have

Combining (2.71) with (2.69), (2.61), (2.65), and the assumptions (2.24) and (2.25), allow us to conclude that

This combined with (2.67) and (2.60) gives

By virtue of (2.70), the fact that $ (s_w)_x = -(s_o)_x $, (2.72), and (2.73), we get the estimates of $ (s_o)_x $ and $ (\rho_o)_x $.

Lemma 2.5. For any $ t\in[0, T^*) $, it holds that

Proof. By virtue of (2.32), we have

where we use $ \int_0^1m\, dx = \tilde{m}, \int_0^1n\, dx = \tilde{n} $ and Hölder inequality in the first inequality, and use the inequality $ \|u_i\|_{\infty}\leq C\|u_{ix}\|_{L^2} $, (2.66), (2.65), and Cauchy inequality in the rest.

Using (2.65) again, and (2.75), we get

Corollary 2.4. For any $ t\in[0, T^*) $, it holds that

Proof. From the equation of $ u_w $, we have

which together with (2.65), (2.60), (2.66), and (2.74) gives

where we use $ W^{1, 1}(0, 1)\hookrightarrow L^\infty(0, 1) $ in the third inequality, and use $ \int ab\, dx \leq C\int a^2\, dx + \varepsilon \int b^2\, dx $ with appropriate choice of $ \varepsilon $ in the last one. This implies

Similarly, we get the estimate of $ (u_o)_{xx} $.

Corollary 2.5. For any $ t\in[0, T^*) $, it holds that

Proof. By using the equation of $ m $, Cauchy inequality, (2.60), and (2.74), we have

Similarly we get

We consider (2.69) with $ \partial_x $ replaced by $ \partial_t $, and use similar analysis as (2.72). Then

where we use (2.65), (2.79), and (2.80). Since $ (s_o)_t = -(s_w)_t $, we get

By virtue of (2.67) where $ \partial_x $ is replaced by $ \partial_t $, (2.79), (2.80), and (2.81), we have

This combined with (2.70) and (2.81) gives

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_iu_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 $ \varepsilon_w, \varepsilon_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 $ \varepsilon_w, \varepsilon_o \rightarrow 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]:

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.004 m^3/d $. The relevant boundary values are: $ s_w(x = 0, t) = 0.8 $ and $ u_w(x = 0, t) = 8.004 m/d $. In addition, 2001 grids are used to simulate this displacement. Water and oil have the same properties such as density and viscosity. Fluid-fluid interaction effect is ignored here by setting $ I = 0 $ in the correlations (4.132). The corresponding initial water saturation profile is shown in Fig. 1 (right figure).

3.1.   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_iu_{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, $ \varepsilon_w = \varepsilon_o = 10^7 $ and $ \varepsilon_w = \varepsilon_o = 10^6 $ (upper row) and compare with the corresponding results from Coclite et al. [9] (lower row). Main observations are:

(ⅰ) 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;

(ⅱ) 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;

(ⅲ) The solution shows a clear sensitivity to the magnitude of $ \varepsilon_o, \varepsilon_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 $ \varepsilon_o, \varepsilon_w $.

The classical Buckley-Leverett model solution (i.e., $ \varepsilon_w = \varepsilon_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_iu_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 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 $ \varepsilon_w, \varepsilon_o $ get smaller. Both schemes seem to reflect convergence toward the solution of the inviscid model (i.e., $ \varepsilon_w = \varepsilon_o = 0 $) with a considerably faster convergence produced by the scheme with viscous term based on interstitial fluid $ u_i $.

Case 2. In Fig. 4 we show simulation results when we vary the internal relation between $ \varepsilon_w $ and $ \varepsilon_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 $ \varepsilon_w $ from $ 10^4 $ to $ 10^6 $ with a constant $ \varepsilon_o $ has a dramatic effect when $ \varepsilon_o $ is small (i.e., $ 10^4 $) whereas the effect is rather small when $ \varepsilon_o $ is large (i.e., $ 10^6 $). We refer to Fig. 4 for simulation results.

Case 3. Now, we move to another case which has a different initial condition but still with the same injection rate of water, $ 8.004 m^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 $ \varepsilon $ 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 pressure-density relation: $ \rho_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).

4.   Concluding remarks

Some main observations from the investigations of the present manuscript are:

● 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 (energy-type 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 $ \varepsilon_w, \varepsilon_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 $ \varepsilon_w, \varepsilon_o $ and cannot be used to ensure convergence to the inviscid model. Hence, this remains an interesting open problem.

Acknowledgments

We are grateful for instructive comments from the anonymous 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

where $ A = \max\{\frac{1}{\epsilon_o}, \frac{1}{\epsilon_w}\} $, $ A_1 $ and $ T_0 $ satisfy (4.101) and (4.87), (4.88), and (4.117), respectively.

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:

and

for $ k = 1, 2, 3, \cdot\cdot\cdot $, where $ (u^0_w, u^0_o) = (0, 0) $, $ s_w^k = s_w(m^k, n^k) $, $ s_o^k = s_o(m^k, n^k) $, $ P_w^k = P_w(m^k, n^k) $, $ P_o^k = P_o(m^k, n^k) $, $ \hat{k}_w^k = \hat{k}_w(m^k, n^k) $, $ \hat{k}_o^k = \hat{k}_o(m^k, n^k) $, $ \hat{k}^k = \hat{k}(m^k, n^k) $, and

Step 2. Boundedness of the sequence.

Assume that $ u^{k-1}_w, u^{k-1}_o\in S $. To prove $ u^{i}_w, u^{i}_o\in S $ for all $ i = 0, 1, 2, 3, ... $, it suffices to prove $ u^{k}_w, u^{k}_o\in S $.

In fact, as a consequence of that $ u^{k-1}_w, u^{k-1}_o\in S $, we have

for $ t\in[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\le T_1 $ for some $ T_1>0 $ which reflects smallness on time, such that

and

In fact, (4.87)$ _1 $ as well as (4.88) are straightforward consequences of the fact that $ m^k $ and $ n^k $ are transported by the smooth vector fields $ u_w^{k-1} $ and $ u_o^{k-1} $. Moreover, (4.87)$ _{2, 3} $ can then be deduced from the regularity of $ (m, n)\mapsto \rho_s(m, n) $ and of $ \rho_s\mapsto P_s(\rho_s) $ for $ s = w, o $. See Step 2, Section 2.2 in [16] for more details.

Multiplying (4.85)$ _3 $ by $ u_{o}^k $ and $ u_{oxx}^k $ 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 $ \delta>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

and

where we take $ \delta = \frac{k_0}{10C} $.

Putting (4.92) into (4.93), we have

where

In view of (4.92) and (4.94), we have

Similar to (4.89) and (4.90), we have

and

which yield

and

where

In view of (4.96) and (4.97), we have

Letting

which can be done by assuming for instance that $ \varepsilon_o $ and $ \varepsilon_w $ are large enough, we obtain from (4.95) and (4.98)

where we choose that

Consequently, we have $ u_w^k, u_o^k\in S $ for all $ k = 1, 2, 3, \cdot\cdot\cdot $, if we assume that (4.99) is satisfied for $ A_1 $ given by (4.101) and for $ T_0\le T_1 $.

Step 3. Compactness arguments.

By virtue of the estimates uniformly for $ k $ in Step 2, there exist a subsequence $ k_i $ ($ i = 1, 2, 3, ... $) and a $ (u_w, u_o, m, n) $, such that

as $ k_i\rightarrow\infty $, where

and $ n_t, m_t\in L^\infty([0, T_0];L^{2}) $. Using the Aubin-Lions' compactness theorem, we can obtain some strong convergence. More precisely, there exists a subsequence still denoted by $ k_i $ without loss of generality ($ i = 1, 2, 3, ... $), such that

as $ k_i\rightarrow\infty $. It follows from (4.103) and (4.87) that

Step 4. Convergence of $ (u^{k_i-1}_w, u^{k_i-1}_o) $.

We are going to investigate the convergence of the neighbor sequence of $ (u^{k_i}_w, u^{k_i}_o) $, i.e., $ (u^{k_i-1}_w, u^{k_i-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. Denote $ \bar{m}^{k+1} = m^{k+1}-m^{k} $ and $ \bar{n}^{k+1} = n^{k+1}-n^{k} $. Then, from (4.85)$ _{1, 2} $ it follows

and

for $ (x, t)\in[0, 1]\times[0, T_0] $.

Using (4.104), we have

where $ \bar{C} $ 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_{w}^{k}\in S $, and the Poincaré inequality. Similarly, we have

Using (4.85)$ _3 $, we have

where $ k = 2, 3, 4, \cdot\cdot\cdot $.

Denote $ \bar{u_o}^{k+1} = u_o^{k+1}-u_o^{k} $ and $ \bar{u_w}^{k+1} = u_w^{k+1}-u_w^{k} $. Then similar to the estimate of $ u_o^k $ (see (4.89)), the $ \|(u^{k}_{o}-u^{k-1}_{o})_x\|_{L^2} $ can be evaluated as follows:

where $ k = 2, 3, 4, \cdot\cdot\cdot $. Note that

This, together with Hölder inequality, (4.109) and the boundedness of $ P^{k}_{ox} $ and $ [s^{k-1}_o]_x $ in $ L^2 $, implies that

where $ k = 2, 3, 4, \cdot\cdot\cdot $.

Similarly, by virtue of (4.85)$ _4 $, we have

and

where $ k = 2, 3, 4, \cdot\cdot\cdot $.

Putting (4.110) into (4.112), we have

for $ k = 2, 3, 4, \cdot\cdot\cdot $, where we have used (4.87) and Sobolev inequality.

Combining (4.106), (4.107), (4.110) and (4.113) with Cauchy inequality, we have

which yields

where

Then, (4.114) together with Gronwall inequality yields

Integrating (4.114) over $ [0, T_0] $, and using (4.115), we have

Letting

for small $ T_0 $ and large $ \varepsilon_o $, $ \varepsilon_w $. We note that (4.116) and (4.115) imply that

Thus, (4.118) combined with (4.102)$ _1 $ implies that

as $ k_i\rightarrow\infty $.

Step 5. Conclusion.

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., $ \rho_w $ and $ \rho_o $ are constant.

● porosity $ \phi $ 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 use $ \hat{k}_{ow} = \hat{k} $. We can solve for $ \mathbf{u}_w $ and $ \mathbf{u}_o $ from the two momentum balance equations in (4.120) and find that

That is, we find for $ \mathbf{U}_i = \phi s_i \mathbf{u}_i $:

where $ \Delta\rho = \rho_w-\rho_o $ and with generalized mobility functions $ \hat{\lambda}_i $ ($ i = w, o, T $) of the form

Note that the resistance force coefficients $ \hat{k}_w(s_w), \hat{k}_o(s_w), \hat{k}(s_w) $ typically are functions of water saturation $ s_w $. From (4.122) we find

Remark 4.1. A fundamental difference between the above expression (4.124) for the total superficial velocity $ \mathbf{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\nabla 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 $ \mathbf{u}_i $ instead of Darcy velocity $ \mathbf{U}_i = \phi s_i \mathbf{u}_i $.

The incompressible, inviscid case $ \varepsilon_w = \varepsilon_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 $ \varepsilon_w = \varepsilon_o = 0 $. We observe that (4.124) gives

where

are the standard fractional flow functions which satisfy that $ \hat{f}_w+\hat{f}_o = 1 $. Using (4.126) in (4.122) we get

where $ \Delta\rho = \rho_w-\rho_o $ and $ \hat{h}(s_w) $ is defined by

From this, it also follows that

Consequently, the model (4.120) has been reduced to solving the following conservation law

where $ \mathbf{U}_w = \mathbf{U}_w(s_w) $ is given by (4.128).

Remark 4.2. The model (4.131) combined with (4.128) (last line), (4.127), (4.129), and the generalized mobility functions as defined by (4.123), recovers the classical incompressible immiscible model except that we now also include for an additional water-oil drag force effect through the term $ \hat{k} $. Setting this term to zero reproduces exactly the classical formulation. To obtain closed expressions for the generalized mobility functions $ \hat{\lambda}_i $ as well as for $ \hat{h}(s_w) $, we must specify appropriate functional correlations for (ⅰ) the water-rock resistance force $ \hat{k}_w $; (ⅱ) the oil-rock resistance force $ \hat{k}_o $; (ⅲ) the water-oil drag force effect $ \hat{k} $. In the recent work [31] we used correlations of the following form:

where $ \mu_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 $ \mathbf{U}_w $ given by (4.128) and $ \mathbf{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 $ \mathbf{U}_T\hat{f}_w $ and counter-current flow (i.e., flow in the opposite direction) represented by $ -\hat{h}(s_w)\Delta \rho \mathbf{g} $ and $ \hat{h}(s_w)\nabla P_{c}(s_w) $ (compare with (4.130)). The fact that the fluid-fluid interaction term $ \hat{k} $ is explicitly accounted for and appears in $ \hat{h} $ (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 $ \rho_w $ and the $ m $ mass balance with $ \rho_o $ and summing the two resulting equations

or

with

Clearly,

Consequently,

Using this in (4.134) we get

that is,

Note that

so that

or

Similarly, we have for $ P_{w} $:

Clearly,

Consequently,

Using this in (4.140) we get

that is,

Note that

so that

or

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\rho_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:

with

We refer to Appendix C and (4.144) which gives the pressure evolution equation (4.145)$ _2 $. Note that $ s_w, s_o, n, P_o $ are determined by

We solve (4.145) on our domain $ \Omega $ with boundary conditions

and initial condition

System of ODEs. We consider the domain $ \Omega = [0, 1] $ and introduce a grid of $ N_x $ cells with nodes $ x_j $ placed at the center of the cells

and cell interfaces $ x_{j+1/2} $ at the cell interfaces

where $ \Delta x = 1/N_x $. We introduce the approximate mass and pressure $ \{m_j(t)\}_{j = 1}^{N_x} $ and $ \{P_{w, j}(t)\}_{j = 1}^{N_x} $ associated with the nodes $ \{x_j\}_{j = 1}^{N_x} $ whereas the approximate velocities $ \{u_{w, j+1/2}\}_{j = 0}^{N_x} $ and $ \{u_{o, j+1/2}\}_{j = 0}^{N_x} $ are associated with the cell interfaces $ \{x_{j+1/2}\}_{j = 0}^{N_x} $. In the following we describe a semi-discrete version of (4.145).

Step 1: Mass transport. We solve for $ m_j(t) $ by considering the following ODE corresponding to (4.145)$ _1 $:

where

This can also be expressed as

Step 2: Computation of velocities and pressure. Next, we solve for $ P_{w, j}(t) $ and $ u_{w, j+1/2}(t) $ and $ u_{o, j+1/2}(t) $ by considering the following ODE system corresponding to (4.145)$ _{2, 3, 4} $:

which is combined with the momentum balance equations

Here we note that the average $ s_{w, j+1/2} $ in (4.153) is based on upwind relatively $ u_{w, j+1/2} $

Similarly, for $ s_{o, j+1/2} $ and for the interaction terms $ \hat{k}_{w, j+1/2} $ and $ \hat{k}_{o, j+1/2} $. In addition, $ \hat{k}_{j+1/2} $ is based on upwind relatively $ u_{w, j+1/2} $ and $ u_{o, j+1/2} $ as follows:

Moreover, $ [nu_o]_{j+1/2} $ and $ [mu_w]_{j+1/2} $ appearing in (4.152) employ upwind as described in (4.151). Now, we are in a position where we can describe a fully discrete model.

A fully discrete scheme. We assume that we have given $ (m_j^{k}, P_{w, j}^k, u_{w, j}^k, u_{o, j}^k) $. We then compute the approximate solution at time $ t^{k+1} $ expressed by $ (m_j^{k+1}, P_{w, j}^{k+1}, $ $ u_{w, j}^{k+1}, u_{o, j}^{k+1}) $ as follows:

Step 1: Mass transport.

where

Having computed $ m_j^{k+1} $ we can compute an updated water saturation $ s_{w, j}^{k+1/2} $ given by

Similarly, we compute updated mass $ n_j^{k+1/2} = (1-s_{w, j}^{k+1/2})\rho_o(P_{o, j}^{k+1/2}) $ and $ P_{o, j}^{k+1/2} = P_{w, j}^k + P_c(s_{w, j}^{k+1/2}) $, according to (4.147), which are needed to evaluate coefficients in the next step.

Step 2: Computation of velocities and pressure. Next, we solve simultaneously for $ P_{w, j}^{k+1} $ and $ u_{w, j+1/2}^{k+1} $ and $ u_{o, j+1/2}^{k+1} $ by considering the following algebraic system

which is combined with the momentum balance equations

Equipped with $ (P_{w, j}^{k+1}, u_{w, j+1/2}^{k+1}, u_{o, j+1/2}^{k+1}) $ we can now update the saturation $ s_{w, j} $ by

from which we also compute the updated oil mass $ n_j^{k+1} $ and pressure $ P_{o, j}^{k+1} $ via (4.147). If necessary, we may repeat step 2 to improve the accuracy before we proceed to next time level.

Remark 4.3. The upwind discretization of $ [n^{k+1/2} u_o^{k+1}]_{j+1/2} $ and $ [m^{k+1} u_w^{k+1}]_{j+1/2} $ appearing in (4.159) are based on "old" velocities $ u_{o, j+1/2}^{k} $ and $ u_{w, j+1/2}^{k} $.

A semidiscrete scheme for the incompressible model. When fluids are incompressible the model (2.14) takes the following form with unknown variables $ (s_w, P_w, u_w, u_o) $

subject to the boundary condition

and initial condition

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.162) is given by (after a summation of the two mass balance equation)

This formulation is consistent with and follows directly from (4.145) by letting $ C_w, C_o\rightarrow \infty $ (i.e., the fluids become incompressible). This is a consequence of the fact that $ \tilde{\eta}\rightarrow \infty $ and $ \tilde{\rho}_o\rightarrow \rho_o $, see (4.146).

Step 1: Mass transport. We solve the following ODE for $ s_{w, j}(t) $ corresponding to (4.165)$ _1 $:

where

Step 2: Computation of velocities and pressure. Next, we solve for $ P_{w, j}(t) $ and $ u_{w, j+1/2}(t) $ and $ u_{o, j+1/2}(t) $ by considering the algebraic system corresponding to (4.165)$ _{2, 3, 4} $:

which is combined with the momentum balance equations

Here we note that the average $ s_{w, j+1/2} $ in (4.169) is based on upwind relatively $ u_{w, j+1/2} $

Similarly, for $ s_{o, j+1/2} $ and for the interaction terms $ \hat{k}_{w, j+1/2} $ and $ \hat{k}_{o, j+1/2} $. In addition, $ \hat{k}_{j+1/2} $ is based on upwind relatively $ u_{w, j+1/2} $ and $ u_{o, j+1/2} $

Moreover, $ [s_ou_o]_{j+1/2} $ and $ [s_wu_w]_{j+1/2} $ appearing in (4.168) employ upwind as described in (4.170).

A fully discrete scheme for the incompressible model. We can now proceed with a description of a full discrete scheme for the incompressible case which bears clear similarity to the scheme for the compressible model.

Step 1: Mass transport.

where

Having computed $ s_{w, j}^{k+1} $ we can compute pressure and velocities simultaneously at time level $ k+1 $.

Step 2: Computation of velocities and pressure. We solve for $ P_{w, j}^{k+1} $ and $ u_{w, j+1/2}^{k+1} $ and $ u_{o, j+1/2}^{k+1} $ by considering the following algebraic system

which is combined with the momentum balance equations

Remark 4.4. The upwind discretization of $ [s_w^{k+1} u_w^{k+1}]_{j+1/2} $ and $ [s_o^{k+1} u_o^{k+1}]_{j+1/2} $ appearing in (4.172) are based on "old" velocities $ u_{w, j+1/2}^{k} $ and $ u_{o, j+1/2}^{k} $.