\documentclass[reqno]{amsart} \usepackage{hyperref} \AtBeginDocument{{\noindent\small \emph{Electronic Journal of Differential Equations}, Vol. 2010(2010), No. 122, pp. 1--33.\newline ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu \newline ftp ejde.math.txstate.edu} \thanks{\copyright 2010 Texas State University - San Marcos.} \vspace{9mm}} \begin{document} \title[\hfilneg EJDE-2010/122\hfil Compressible immiscible two phase flow] {Solutions to a model for compressible immiscible two phase flow in porous media} \author[Z. Khalil, M. Saad\hfil EJDE-2010/122\hfilneg] {Ziad Khalil, Mazen Saad} % in alphabetical order \address{Ziad Khalil \newline Ecole Centrale de Nantes, Laboratoire de Math\'ematiques Jean Leray, UMR CNRS 6629, 1, rue de la No\'e, 44321 Nantes, France} \email{Ziad.Khalil@ec-nantes.fr} \address{Mazen Saad \newline Ecole Centrale de Nantes, Laboratoire de Math\'ematiques Jean Leray, UMR CNRS 6629, 1, rue de la No\'e, 44321 Nantes, France} \email{Mazen.Saad@ec-nantes.fr} \thanks{Submitted May 7, 2010. Published August 30, 2010.} \thanks{Partially supported by GNR MoMaS} \subjclass[2000]{35K57, 35K55, 92D30} \keywords{Degenerate system; nonlinear parabolic system; compressible flow; \hfill\break\indent porous media} \begin{abstract} In this article, we study the existence of solutions to a nonlinear degenerate system modelling the displacement of two-phase compressible immiscible flow in a three dimensional porous media. The aim of this work is to treat the model with its general form with the whole nonlinear terms. Especially, we consider the case where the density of each phase depends on its corresponding pressure. We derive new energy estimates on velocities, saturations and pressures to treat the degeneracy of the system. A compactness result is shown for degenerate systems. \end{abstract} \maketitle \numberwithin{equation}{section} \newtheorem{theorem}{Theorem}[section] \newtheorem{proposition}[theorem]{Proposition} \newtheorem{lemma}[theorem]{Lemma} \allowdisplaybreaks \section{Introduction, assumptions and main results} \label{intro} The mathematical and numerical study of the miscible flow models has been investigated in \cite{A-AHZ91,A-AHZ96,A-FL92} and recently in \cite{amaziane2,choquet1-08,choquet2-08,choquet3-08}. The immiscible and incompressible flows have been treated by many authors \cite{A-AHZ91,arbogast,B-CJ86,A-FG00,A-DEH06,A-FLT,A-FS93}. For two immiscible compressible flows, we refer to \cite{A-GS04,A-GS08-2}, and recently \cite{A-GS08} and \cite{brull}. The immiscible flow models developed by \cite{A-GS04,A-GS08,A-GS08-2} use the feature of global pressure even if the density of each phase depends on its own pressure, then the context was to assume small capillary pressure so that the densities are assumed to depend on the global pressure, recently and under that context Galusinski, Saad \cite{A-GS08} obtained an existence result of solutions. The employed global pressure is that defined by Chavent et al. \cite{B-CJ86} for incompressible immiscible two phases, there is no assumptions to define this pressure. In \cite{amaziane1}, a new notion of global pressure is introduced especially for two compressible immiscible fluids, the new global pressure is defined implicitly and depends on state law of density. In this work, we consider the two compressible immiscible flows model studied in \cite{A-GS08}. The model is treated in its general form under the physical assumption that the density of each phase depends on its own pressure. The mathematical study of this model is based on new energy estimates on the velocity of each phase. The main idea consists in deriving from degenerate estimates on pressure of each phase, which not allowed straight bound on pressures, an estimate on global pressure and degenerate capillary term. An appropriate compactness lemma is shown with the help of the feature of global pressure to pass from non-degenerate case to the degenerate case. We give below the basic model written in variables pressures and saturations. The equations describing the immiscible displacement of two compressible fluids are given by the following mass conservation of each phase $(i=1,2)$: $$\label{conslaw} \phi(x)\partial_{t}( \rho_{i}(p_i)s_i)(t,x) + \operatorname{div} (\rho_{i}(p_i) \mathbf{V}_{i})(t,x) +\rho_i(p_i) s_i f_{P}(t,x) = \rho_i(p_i) s^I_i f_{I}(t,x),$$ where $\phi$ is the porosity of the medium, $\rho_i$ and $s_i$ are respectively the density and the saturation of the $i ^{th}$ fluid. The velocity of each fluid $\mathbf{V}_i$ is given by the Darcy law: $$\mathbf{V}_{i}(t,x) = - \mathbf{K}(x) \frac{k_{i}(s_i(t,x))}{\mu_{i}}\big( \nabla p_{i}(t,x)-\rho_i(p_i)\mathbf{g}\big),\quad i = 1, 2.$$ where $\mathbf{K}(x)$ is the permeability tensor of the porous medium at point $x$ to the fluid under consideration, $k_i$ the relative permeability of the $i ^{th}$ phase, $\mu_i$ the constant $i$-phase's viscosity, $p_i$ the $i$-phase's pressure and ${\bf g }$ is the gravity term. Here the functions $f_{I}$ and $f_{P}$ are respectively the injection and production terms. Note that in equation \eqref{conslaw} the injection term is multiplied by a known saturation $s^I_i$ corresponding to the known injected fluid, whereas the production term is multiplied by the unknown saturation $s_i$ corresponding to the produced fluid. By definition of saturations, one has $$\label{1} s_{1}(t,x) + s_{2}(t,x) = 1.$$ The curvature of the contact surface between the two fluids links the jump of pressure of the two phases to the saturation by the capillary pressure law in order to close the system (\eqref {conslaw}-\eqref{1}, $$\label{cpl} f(s_1(t,x)) = p_{1}(t,x) - p_{2}(t,x).$$ the application $s_1\mapsto f(s_1)$ is non-decreasing, $(\frac{df}{ds_1}(s_1) > 0$, for all $s_1 \in[0,1])$, and usually $f(s_1=1)=0$, in the case of two phases, when the wetting fluid is at its maximum saturation. In this study we consider that the index $i=1$ represents the wetting fluid, and for this choice capillary pressure vanishes when $s_1=1$. This point is crucial in determining the spaces that the saturation of each phase belongs to. We take the capillary pressure function $f$ as considered in \cite{B-CJ86}, defined on $[0,1]$, increasing and $f(1)=0$. In section \ref{sec-maint} we will use the feature of global pressure. For that let us denote, \begin{gather*} M_i(s_i) = k_{i}(s_i)/\mu_{i} \quad i-\text{phase's mobility,} \\ M(s_1) =M_1(s_1)+M_2(1-s_1) \quad \text{the total mobility},\\ \mathbf{V} =\mathbf{V}_1 + \mathbf{V}_2 \quad \text{the total velocity.} \end{gather*} As in \cite{B-CJ86,A-SA97,A-GS08} we can express the total velocity in terms of $p_2$ and $f(s_1)$. We have $$\mathbf{V} =- \mathbf{K}M(s_1) \Big( \nabla p_{2}+\frac{M_1(s_1)}{M(s_1)}\nabla f(s_1)\Big) + \mathbf{K}(M_1(s_1)\rho_1(p_1)+ M_2(s_2)\rho_2(p_2))\mathbf{g}.$$ Defining the functions $\tilde{p}(s_1)$, $\bar{p}(s_1)$ such that $$\label{dptilde} \tilde{p}'(s_1) = \frac{M_1(s_1)}{M(s_1)}f'(s_1),\quad \bar{p}'(s_1) = -\frac{M_2(s_2)}{M(s_1)}f'(s_1),$$ the global pressure is then defined as in \cite{B-CJ86} $$\label{defp} p = p_{2}+\tilde{p}(s_1)= p_{1}+\bar{p}(s_1)$$ Thus, the total velocity can be expressed as \begin{gather*} \mathbf{V} = - \mathbf{K}M(s_1)\nabla p + \mathbf{K}(M_1(s_1)\rho_1(p_1)+ M_2(s_2)\rho_2(p_2))\mathbf{g},\\ \mathbf{V}_{i} =-\mathbf{K}M_i(s_i)\nabla p - \mathbf{K} \alpha(s_1) \nabla s_i +\mathbf{K}M_i(s_i)\rho_i(p_i)\mathbf{g}. \end{gather*} where $\alpha (s_1)= \frac{M_1(s_1)M_2(s_2)}{M(s_1)} \frac{df}{ds}(s_1)\geq 0$. Define $$\label{defbeta} \beta(s)=\int_0^s \alpha(\xi)d\xi.$$ In this paper we do not use this concept of writing the total and each velocity in terms of the global pressure and one saturation, but just to show the source of definitions of some functions. We detail the physical context by introducing the boundary conditions, the initial conditions and some assumptions on the data of the problem. Let $T>0$, fixed and let $\Omega$ be a bounded open set of $\mathbb{R}^d$ ($d\ge 1$), with Lipschitz boundary. We set $Q_T=(0,T)\times\Omega$, $\Sigma_T=(0,T)\times\partial\Omega$. To the system \eqref{conslaw}-\eqref{1}-\eqref{cpl} ($i=1,2$), we add the following mixed boundary conditions and initial conditions. We consider the boundary $\partial \Omega=\Gamma_{1} \cup \Gamma_{\rm imp}$, where $\Gamma_{1}\ne \emptyset$ denotes the injection boundary of the first phase and $\Gamma_{\rm imp}$ the impervious one. $$\label{bc} \begin{gathered} p_1(t,x) = 0, \quad p_2(t,x) = 0 \quad \text{on } \Gamma_{1}\\ \mathbf{V}_1\cdot\mathbf{n} = \mathbf{V}_2\cdot\mathbf{n} = 0\quad \text{on }\Gamma_{\rm imp}, \end{gathered}$$ where $\mathbf{n}$ is the outward normal to the boundary $\Gamma_{\rm imp}$. The initial conditions are defined on pressures $$\label{ic} p_i(0,x) = p_i^0(x) \quad\mbox{in } \Omega, \quad i=1,2.$$ Next, we introduce some physically relevant assumptions on the coefficients of the system. \begin{enumerate} \item[(H1)] The porosity $\phi \in L^{\infty}(\Omega)$ and there is two positive constants $\phi_{0}$ and $\phi_{1}$ such that $\phi_0 \le \phi(x)\le \phi_1$ almost everywhere $x\in\Omega$. \item[(H2)] The tensor $\mathbf{K}$ belongs to $(L^{\infty}(\Omega))^{d\times d}$. Moreover, there exist two positive constants $k_0$ and $k_\infty$ such that \begin{equation*} \|\mathbf{K}\|_{(L^\infty(\Omega))^{d\times d}}\le k_\infty,\quad (\mathbf{K}(x)\xi,\xi)\ge k_0|\xi|^2, \quad \text{for all } \xi\in \mathbb{R}^d, \text{ a.e. } x \in \Omega. \end{equation*} \item[(H3)] The functions $M_1$ and $M_2$ belong to ${\mathcal C}^0([0,1]; \mathbb{R}^+)$, $M_1(s_1=0)=0$ and $M_2(s_2=0)=0$. In addition, there is a positive constant $m_0$, such that, for all $s_1\in [0,1]$, $M_1(s_1)+M_2(s_2) \ge m_0;\quad \text{with }s_2=1-s_1.$ \item[(H4)] $(f_{P},f_{I})\in (L^2(Q_T))^2$, $f_{P}(t,x)$, $f_{I}(t,x) \ge 0$ almost everywhere $(t,x)\in Q_T$, $s^I_i (t,x) \ge 0$ $(i=1,2)$ and $s^I_1 (t,x)+s^I_2 (t,x) =1$ almost everywhere in $(t,x)\in Q_T$. \item[(H5)] The densities $\rho_i$ ($i=1,2$) are ${\mathcal C}^2(\mathbb{R})$, increasing and there exist two positive constants $\rho_m$ and $\rho_M$ such that $0<\rho_m \leq \rho_i(p_i)\leq \rho_M$, for $i=1,2$. \item[(H6)] The capillary pressure function $f\in {\mathcal C}^1([0,1]; \mathbb{R}^-)$ and $0< \underline{f}\leq\frac{df}{ds}$. \item[(H7)] The function $\alpha \in {\mathcal C}^1([0,1];\mathbb{R}^+)$ satisfies $\alpha(s)>0$ for $00$, there exists $({p_1^\eta},{p_2^\eta})$ satisfying \begin{gather*} {p_i^\eta}\in L^2(0,T;H^1_{\Gamma_1}(\Omega)),\quad {s_1^\eta} \in L^2(0,T;H^1(\Omega)),\quad {s_2^\eta} \in L^2(0,T;H^1_{\Gamma_1}(\Omega)),\\ \phi\partial_{t}(\rho_i({p_i^\eta}){s_i^\eta})\in L^2(0,T;(H^1_{\Gamma_1}(\Omega))'),\quad \rho_i({p_i^\eta}){s_i^\eta})\in C^0(0,T;L^2(\Omega)),\\ 0\leq {s_i^\eta}(t,x)\leq1 \quad\text{ a.e. in } Q_T,\; i=1,2, \end{gather*} for all $\varphi_i \in L^2(0,T;H^1_{\Gamma_1}(\Omega))$, $i=1, 2$, \label{non-degv1} \begin{aligned} &\langle\phi \partial_t(\rho_i({p_i^\eta}) {s_i^\eta}), \varphi_i\rangle +\int_{Q_T}\mathbf{K} M_i({s_i^\eta})\rho_i({p_i^\eta}) { \nabla} {p_i^\eta} \cdot\nabla \varphi_i \,dx\,dt \\ &-\int_{Q_T}\mathbf{K} M_i({s_i^\eta})\rho_i^2({p_i^\eta})\mathbf{g} \cdot\nabla \varphi_i \,dx\,dt +(-1)^{i+1}\eta \int_{Q_T}\rho_i({p_i^\eta})\nabla ({p_1^\eta} -{p_2^\eta}) \cdot \nabla \varphi \,dx\,dt \\ &+\int_{Q_T}\rho_i({p_i^\eta}){s_i^\eta} f_{P}\varphi_i \,dx\,dt =\int_{Q_T}\rho_i({p_i^\eta})s_i^If_{I}\varphi_i \,dx\,dt \end{aligned} where the symbol $\langle\cdot,\cdot\rangle$ represents the duality product between $L^2(0,T;(H^1_{\Gamma_1}(\Omega))')$ and $L^2(0,T;H^1_{\Gamma_1}(\Omega))$. \end{theorem} The first step to establish theorem \ref{non-degt} is based on a time discretization scheme of \eqref{conslawnondeg1}-\eqref{conslawnondeg2}. For that, let $\rho ^\star_i$ and $s^\star_i$ be the values of the $h$-translated in time of $\rho_i(p_i)$ and $s_i$, respectively, $i=1,2$, we state the following theorem. \begin{theorem}[Non-degenerate elliptic system] \label{deg-ell-t} Let {\rm (H1)--(H6)} hold. Let $(p_1^0,\, p_2^0)$ belongs to $L^2(\Omega)\times L^2(\Omega)$. Then for all $h>0$, there exists $(p_1^h,p_2^h)=(p_1^{\eta,h},p_2^{\eta,h})$ satisfying \begin{gather*} p_1^h\in H^1_{\Gamma_1}(\Omega),\quad p_2^h\in H^1_{\Gamma_1}(\Omega),\quad s_1^h \in H^1(\Omega),\quad s_2^h \in H^1_{\Gamma_1}(\Omega),\\ 0\leq s_i^h(t,x)\leq 1 \quad\text{a.e. in } Q_T, \; i=1,2, \end{gather*} for all $\varphi_i \in H^1_{\Gamma_1}(\Omega)$, $i=1,2$, \label{degv-ell1} \begin{aligned} &\int_\Omega \phi\frac{ \rho_i({p_i^h}){s_i^h}-\rho ^\star_i s^\star_i}{h} \varphi_i \,dx +\int_{\Omega}\mathbf{K} M_i({s_i^h} )\rho_i({p_i^h} ) { \nabla} {p_i^h} \cdot\nabla \varphi_i \,dx \\ &-\int_{\Omega}\mathbf{K} M_i({s_i^h} )\rho_i^2 ({p_i^h} )\mathbf{g} \cdot\nabla \varphi_i \,dx \\ &+(-1)^{i+1}\eta \int_{\Omega}\rho_i({p_i^h})\nabla ({p_1^h} -{p_2^h}) \cdot \nabla \varphi_i\,dx +\int_{\Omega}\rho_i({p_i^h}){s_i^h} f_{P}\varphi_i \,dx\\ &=\int_{\Omega}\rho_i({p_i^h})s_i^I f_{I}\varphi \,dx, \end{aligned} \end{theorem} The rest of the paper is organized as follows. In the next section we deal with the time discrete model to prove Theorem \ref{deg-ell-t} in two steps. The first step deals with an elliptic system with non degenerate mobilities, $M_i^\epsilon=M_i+\epsilon$ with $\epsilon>0$, in this step we apply a suitable fixed point theorem, Leray-Schauder, to get weak solution. The second step is to pass to the limit as $\epsilon$ goes to zero depending on a suitable uniform estimate (w. r. to $\epsilon$), and a maximum principle ensures the positivity of saturations which achieves the proof of theorem \ref{deg-ell-t}. In the third section we introduce a sequence of solutions solving \eqref{degv-ell1}. This choice is motivated by the fact that no evolution have to be considered in a first step. The problem of degeneracy of evolution term is temporarily sat aside. Furthermore, the maximum principle is conserved on saturation after the passage to the limit on in the non linear variational elliptic system. The last section is devoted to pass from non-degenerate case to degenerate case through a compactness {\it lemma} which allow us with the help of some estimates to pass the limit and end the proof of existence of weak solutions of the system under consideration. The next section is devoted to the analysis of the elliptic problem. \section{Study of a nonlinear elliptic system (proof of theorem \ref{deg-ell-t})} Having in mind a time discretization of \eqref{conslawnondeg1}-\eqref{conslawnondeg2}, we are concerned with the following system, \label{nles01} \begin{aligned} &\phi\frac{ \rho_i(p_i)s_i-\rho ^\star_i s^\star_i}{h} - \operatorname{div} (\mathbf{K}\rho_i(p_i)M_i(s_i) \nabla p_i) +\operatorname{div}(\mathbf{K}M_i(s_i)\rho^2_i \mathbf{g})\\ &- (-1)^{i+1}\eta \operatorname{div} (\rho_i(p_i) \nabla(p_1-p_2)) +\rho_i(p_i)s_if_{P}\\ &= \rho_i(p_i)s_i^If_{I}\;\text{ in }\Omega. \end{aligned} Before establishing theorem \ref{deg-ell-t} which is the main purpose of this section, we introduce the existence of solutions of system \eqref{nles01}, when the mobilities $M_i$, $(i=1,2)$ are replaced by a non-degenerate positive functions, $$M_i^\epsilon=M_i+\epsilon,\quad i=1,2,\quad \text{and } \epsilon>0,$$ which reinforce the passage to the limit in another regularization which is the trunk high frequencies of nonlinear elliptic term in pressure $p_2$ in the equation \eqref{nles01}. Let $\mathcal{P}_N$ be the orthogonal projector of $L^2(\Omega)$ on the first $N$ eigenvectors of the operator $$p\to -\Delta p$$ with homogeneous Dirichlet boundary conditions. The projector $\mathcal{P}_N$ appears in \eqref{ellvNeps1} to make regular the implied term. The necessity of this regularization appears in the coming proposition in order to define the operator which we apply on the Leray-Schauder fixed point theorem. The addition of such $\epsilon$ to the mobilities lead to the loss of maximum principle on the saturations $s_i$ $(i=1,2)$. So the functions $M_1$ and $M_2$ are extended on $\mathbb{R}$ by continuous constant functions outside $[0,1]$ and then are bounded on $\mathbb{R}$. For the same reason we denote, $$Z(s_i)= \begin{cases} 0 &\text{for } s_i \le 0\\ s_i &\text{for } s_i\in [0;1]\\ 1 &\text{for } s_i \ge 1. \end{cases}$$ In the same spirit and in order to write the saturations $s_i$ $(i=1,2.)$ as functions of the principle unknowns $p_1$ and $p_2$ of the system, we extend the capillary pressure function $f$ by continuity and strict monotony outside $[0,1]$ in to $\bar{f}$, this is possible in the case when the capillary function $f$ is bounded, in other words when $\mid f(0)\mid< \infty$, and denote by $s_1={{\bar{f}}}^{-1}(p_1-p_2)$ and $s_2=1-\bar{f}^{-1}(p_1-p_2)$. Existence of solution to \eqref{nles01} is constructed in three steps. The first one consists in studying the following problem for fixed parameters $\epsilon>0$, $N>0$ and $\eta >0$. Then, we are concerned with the regularized elliptic system: \label{ellvNeps1} \begin{aligned} &\int_\Omega \phi\frac{ \rho_1(p_1^{\epsilon,N})Z(s_1^{\epsilon,N})-\rho ^\star_1 s^\star_1}{h} \varphi \,dx +\int_{\Omega}\mathbf{K} M_1^\epsilon(s_1^{\epsilon,N} )\rho_1(p_1^{\epsilon,N} ) { \nabla} p_1^{\epsilon,N} \cdot\nabla \varphi \,dx \\ &-\int_{\Omega}\mathbf{K} M_1(s_1^{\epsilon,N} )\rho_1^2 (p_1^{\epsilon,N} )\mathbf{g} \cdot\nabla \varphi \,dx\\ &+\eta \int_{\Omega}\rho_1(p_1^{\epsilon,N})\nabla (\mathcal{P}_N p_1^{\epsilon,N} - \mathcal{P}_N p_2^{\epsilon,N}) \cdot \nabla \varphi \,dx +\int_{\Omega}\rho_1(p_1^{\epsilon,N})Z(s_1^{\epsilon,N}) f_{P}\varphi \,dx\\ &=\int_{\Omega}\rho_1(p_1^{\epsilon,N})s_1^I f_{I}\varphi \,dx, \end{aligned} \label{ellvNeps2} \begin{aligned} &\int_\Omega \phi\frac{ \rho_2(p_2^{\epsilon,N})Z(s_2^{\epsilon,N})-\rho ^\star_2 s^\star_2}{h} \xi \,dx +\int_{\Omega}\mathbf{K} M_2^\epsilon(s_2^{\epsilon,N} )\rho_2(p_2^{\epsilon,N}) { \nabla} p_2^{\epsilon,N} \cdot\nabla \xi \,dx \\ &-\int_{\Omega}\mathbf{K} M_2(s_2^{\epsilon,N} )\rho_2^2 (p_2^{\epsilon,N})\mathbf{g} \cdot\nabla \xi \,dx\\ &-\eta \int_{\Omega}\rho_2(p_2^{\epsilon,N})\nabla (\mathcal{P}_N p_1^{\epsilon,N} -\mathcal{P}_N p_2^{\epsilon,N}) \cdot \nabla \xi \,dx +\int_{\Omega}\rho_2(p_2^{\epsilon,N})Z(s_2^{\epsilon,N}) f_{P}\xi \,dx\\ &=\int_{\Omega}\rho_2(p_2^{\epsilon,N})s_2^I f_{I}\xi \,dx, \end{aligned} for all $(\varphi,\xi)$ belonging to $H^1_{\Gamma_1}(\Omega)\times H^1_{\Gamma_1}(\Omega)$; with $s_1^{\epsilon,N}={{\bar{f}}}^{-1}(p_1^{\epsilon,N}-p_2^{\epsilon,N})$ and $s_2^{\epsilon,N}=1-s_1^{\epsilon,N}$.\\ The second step concerns the passage to the limit as $N$ goes to infinity in order to recover the full physical diffusion on pressures $p_1$ and $p_2$, while the third one is the passage to the limit as $\epsilon$ approaches zero. \smallskip \noindent{\bf Step 1.} We show for fixed $N>0$ and $\epsilon >0$ existence of solutions to \eqref{ellvNeps1}-\eqref{ellvNeps2}. We omit for the time being the dependence of solutions on parameter $N>0$ and $\epsilon$. \begin{proposition} \label{pleray-sch} Assume $\rho_i^\star s_i^\star$ belongs to $L^2(\Omega)$ and $\rho_i^\star s_i^\star \geq 0$. Then there exists $(p_1,p_2)$ belonging to $H^1_{\Gamma_1}(\Omega)\times H^1_{\Gamma_1}(\Omega)$, solution of \eqref{ellvNeps1}-\eqref{ellvNeps2}. \end{proposition} \begin{proof} The proof is based on the Leray-Schauder fixed point theorem. Let $\mathcal{T}$ be a map from $L^2(\Omega)\times L^2(\Omega)$ to $L^2(\Omega)\times L^2(\Omega)$ defined by $$\mathcal{T}(\overline{p}_1,\overline{p}_2) = (p_1,p_2),$$ where the pair $(p_1,p_2)$ is the unique solution of the system \eqref{nles1b}-\eqref{nles2b} \begin{align} &\int_{\Omega}\phi\frac{\rho_1 (\overline{p}_1)Z( \overline{s}_1)- \rho_1^\star s_1^\star}{h} \varphi \, dx +\int_{\Omega}\mathbf{K} M_1^\epsilon( \overline{s}_1)\rho_1(\overline{p}_1) { \nabla} p_1 \cdot\nabla \varphi \,dx \nonumber\\ &-\int_{\Omega}\mathbf{K} M_1( \overline{s}_1)\rho_1^2(\overline{p}_1)\mathbf{g} \cdot\nabla \varphi \,dx+\eta \int_{\Omega}\rho_1(\overline{p}_1)\nabla (\mathcal{P}_N \overline{p}_1 - \mathcal{P}_N \overline{p}_2) \cdot \nabla \varphi \,dx \nonumber\\ &+\int_{\Omega}\rho_1(\overline{p}_1)Z( \overline{s}_1) f_{P}\varphi \,dx =\int_{\Omega}\rho_1(\overline{p}_1)s_1^If_{I}\varphi \,dx, \label{nles1b} \end{align} \begin{align} &\int_{\Omega}\phi\frac{\rho_2 (\overline{p}_2)Z(\overline{s}_2)- \rho_2^\star s_2^\star}{h} \xi \, dx +\int_{\Omega}\mathbf{K} M_2^\epsilon(\overline{s}_2)\rho_2(\overline{p}_2) { \nabla} p_2 \cdot\nabla \xi \,dx \nonumber \\ &-\int_{\Omega}\mathbf{K} M_2(\overline{s}_2)\rho_2^2(\overline{p}_2)\mathbf{g} \cdot\nabla \xi \,dx-\eta \int_{\Omega}\rho_2(\overline{p}_2)\nabla (\mathcal{P}_N \overline{p}_1 -\mathcal{P}_N \overline{p}_2) \cdot \nabla \xi \,dx \nonumber \\ &+\int_{\Omega}\rho_2(\overline{p}_2)Z(\overline{s}_2) f_{P}\xi \,dx =\int_{\Omega}\rho_2(\overline{p}_2)s_2^If_{I}\xi \,dx, \label{nles2b} \end{align} for all $(\varphi,\xi)$ belonging to $H^1_{\Gamma_1}(\Omega)\times H^1_{\Gamma_1}(\Omega)$, $\overline{s}_1=\bar{f}^{-1}(\overline{p}_1-\overline{p}_2)$ and $\overline{s}_2=1-\bar{f}^{-1}(\overline{p}_1-\overline{p}_2)$. The functions $M_1$ and $M_2$ are the extended mobilities which operates on $\mathbb{R}$. Such extensions of the mobilities $M_i$ $(i=1,2)$, the capillary function $f$ and such bound of the saturations $s_i$ $(i=1,2)$ by introducing the map $Z$ are temporary; we deal it at the end of this section after the passage to the limit in $\epsilon$ by a maximum principle on saturations and then the mobilities, the map $Z$ and the extended capillary function $\bar{f}$ operates only on $[0,1]$ where they have a physical meaning. The system $\eqref{nles1b}-\eqref{nles2b}$ can be written in the form $B_1(p_1,\varphi)=f_1(\varphi)$, $B_2(p_2,\xi)=f_2(\xi)$, where $f_1(\cdot)$, $f_2(\cdot)$ are linear continuous mappings on $H^1_{\Gamma_1}(\Omega)$. Then, apply Lax-Milgram theorem to get the existence of the unique pair $(p_1,p_2)$ in $H^1_{\Gamma_1}(\Omega)\times H^1_{\Gamma_1}(\Omega)$ which ensures that the map $\mathcal{T}$ is well defined on $L^2(\Omega)\times L^2(\Omega)$. \begin{lemma}\label{concomp} The map $\mathcal{T}$ is a continuous operator which maps every bounded subsets of $(L^2(\Omega))^2$ into a relatively compact set $(L^2(\Omega))^2$. \end{lemma} \begin{proof} Consider a sequence $(\overline{p}_{1,n},\overline{p}_{2,n})$ of a bounded set of $L^2(\Omega)\times L^2(\Omega)$ which converges to $(\overline{p}_1,\overline{p}_2)\in L^2(\Omega)\times L^2(\Omega)$, and let us prove that $(p_{1,n},p_{2,n})=\mathcal{T}(\overline{p}_{1,n},\overline{p}_{2,n})$ is bounded in $H^1_{\Gamma_1}(\Omega)\times H^1_{\Gamma_1}(\Omega)$ which converges to $(p_1,p_2)=\mathcal{T}(\overline{p}_1,\overline{p}_2)$. The sequences $p_{1,n},~p_{2,n}$ verify respectively \begin{align}\label{nles1bn} &\int_{\Omega}\phi\frac{\rho_1 (\overline{p}_{1,n})Z(\overline{s}_{1,n})- \rho_1^\star s_1^\star}{h} \varphi \, dx +\int_{\Omega}\mathbf{K} M_1^\epsilon(\overline{s}_{1,n})\rho_1(\overline{p}_{1,n}) { \nabla} p_{1,n} \cdot\nabla \varphi \,dx \nonumber\\ &-\int_{\Omega}\mathbf{K} M_1(\overline{s}_{1,n})\rho_1^2(\overline{p}_{1,n})\mathbf{g} \cdot\nabla \varphi \,dx +\eta \int_{\Omega}\rho_1(\overline{p}_{1,n})\nabla (\mathcal{P}_N \overline{p}_{1,n} - \mathcal{P}_N \overline{p}_{2,n}) \cdot \nabla \varphi \,dx \nonumber\\ &+\int_{\Omega}\rho_1(\overline{p}_{1,n})Z(\overline{s}_{1,n}) f_{P}\varphi \,dx =\int_{\Omega}\rho_1(\overline{p}_{1,n})s_1^If_{I}\varphi \,dx, \end{align} \begin{align}\label{nles2bn} &\int_{\Omega}\phi\frac{\rho_2 (\overline{p}_{2,n})Z(\overline{s}_{2,n})- \rho_2^\star s_2^\star}{h} \xi \, dx +\int_{\Omega}\mathbf{K} M_2^\epsilon(\overline{s}_{2,n})\rho_2(\overline{p}_{2,n}) { \nabla} p_{2,n} \cdot\nabla \xi \,dx \nonumber\\ &-\int_{\Omega}\mathbf{K} M_2(\overline{s}_{2,n})\rho_2^2(\overline{p}_{2,n})\mathbf{g} \cdot\nabla \xi \,dx -\eta \int_{\Omega}\rho_2(\overline{p}_{2,n})\nabla (\mathcal{P}_N \overline{p}_{1,n} -\mathcal{P}_N \overline{p}_{2,n}) \cdot \nabla \xi \,dx \nonumber\\ &+\int_{\Omega}\rho_2(\overline{p}_{2,n})Z(\overline{s}_{2,n}) f_{P}\xi \,dx =\int_{\Omega}\rho_2(\overline{p}_{2,n})s_2^If_{I}\xi \,dx, \end{align} for all $(\varphi,\xi)$ belonging to $H^1_{\Gamma_1}(\Omega)\times H^1_{\Gamma_1}(\Omega)$. Let us take $\varphi=p_{1,n}$ in \eqref{nles1bn}, \label{ad} \begin{aligned} &\int_{\Omega}\mathbf{K} M_1^\epsilon(\overline{s}_{1,n})\rho_1(\overline{p}_{1,n}) { \nabla} p_{1,n} \cdot\nabla p_{1,n} \,dx\\ &= \int_{\Omega}\mathbf{K} M_1(\overline{s}_{1,n})\rho_1^2(\overline{p}_{1,n})\mathbf{g} \cdot\nabla p_{1,n} \,dx\\ &\quad -\eta \int_{\Omega}\rho_1(\overline{p}_{1,n})\nabla (\mathcal{P}_N \overline{p}_{1,n} - \mathcal{P}_N \overline{p}_{2,n}) \cdot \nabla p_{1,n} \,dx\\ &\quad -\int_{\Omega}\phi\frac{\rho_1 (\overline{p}_{1,n})Z(\overline{s}_{1,n}) - \rho_1^\star s_1^\star}{h} p_{1,n} \, dx\\ &\quad -\int_{\Omega}\rho_1(\overline{p}_{1,n})Z(\overline{s}_{1,n}) f_{P}p_{1,n} \,dx \int_{\Omega}\rho_1(\overline{p}_{1,n})s_1^If_{I}p_{1,n} \,dx, \end{aligned} we deduce from the Cauchy-Schwarz inequality that \eqref{ad} reduces to, \label{nabs1n} \begin{aligned} \epsilon k_0 \rho_m \int_\Omega \vert \nabla p_{1,n}\vert^2 \,dx &\le C\Big(1 + \Vert p_{1,n}\Vert_{L^2(\Omega)} +\Vert \nabla p_{1,n}\Vert_{L^2(\Omega)}\\ &\quad +\Vert \nabla \mathcal{P}_N\overline{p}_{1,n} \Vert_{L^2(\Omega)} +\Vert \nabla \mathcal{P}_N\overline{p}_{2,n} \Vert_{L^2(\Omega)}\Big), \end{aligned} where $C$ depends on $\Omega$, $\eta$, $h$, $\phi_1$, $\Vert f_{P} \Vert_{L^2(\Omega)}$, $\Vert f_{I} \Vert_{L^2(\Omega)}$, $\rho_M$, $k_\infty$ and $\Vert\rho_1^\star s_1^\star\Vert_{L^2(\Omega)}$. As, $$\Vert \nabla\mathcal{P}_N \overline{p}_{i,n}\Vert_{L^2(\Omega)} \le c_N \Vert \overline{p}_{i,n} \Vert_{L^2(\Omega)},\quad (i=1,2)$$ where $c_N$ is the square root of the $N${th} eigenvalue of the laplace operator (by considering the set of eigenvalues as increasing sequence), the Poincar\'e and Young inequalities and the estimate \eqref{nabs1n} ensure that the sequence $(P_{1,n})_n$ is uniformly bounded in $H^1_{\Gamma_1}(\Omega)$. Then, taking $\xi=p_{2,n}$ in \eqref{nles2bn}, we deduce similarly that \label{nabs2n} \begin{aligned} \epsilon k_0 \rho_m \int_\Omega \vert \nabla p_{2,n}\vert^2 \,dx &\le C \Big(1+ \Vert p_{2,n}\Vert_{L^2(\Omega)}+ \Vert \nabla p_{2,n}\Vert_{L^2(\Omega)}\\ &\quad +\Vert \nabla \mathcal{P}_N\overline{p}_{1,n} \Vert_{L^2(\Omega)} +\Vert \nabla \mathcal{P}_N\overline{p}_{2,n} \Vert_{L^2(\Omega)}\Big), \end{aligned} where $C$ depends on $\Omega$, $\eta$, $h$, $\phi_1$, $\Vert f_{P} \Vert_{L^2(\Omega)}$, $\Vert f_{I} \Vert_{L^2(\Omega)}$, $\rho_M$, $k_\infty$ and $\Vert\rho_2^\star s_2^\star\Vert_{L^2(\Omega)}$. Then the sequence $(P_{2,n})_n$ is uniformly bounded in $H^1_{\Gamma_1}(\Omega)$. This establishes the relative compactness property of the map $\mathcal{T}$. Furthermore, up to a subsequence, we have the convergence \begin{gather} p_{1,n} \to p_1 \quad\text{ weakly in } H^1_{\Gamma_1}(\Omega),\label{cvp1nw}\\ p_{2,n} \to p_2 \quad\text{ weakly in } H^1_{\Gamma_1}(\Omega),\label{cvp2nw}\\ p_{1,n} \to p_1 \quad\text{ strongly in } L^2 (\Omega)\text{ and a.e. in }\Omega,\label{cvp1ns}\\ p_{2,n} \to p_2 \quad\text{ strongly in } L^2 (\Omega)\text{ and a.e. in }\Omega.\label{cvp2ns} \end{gather} To complete the proof of continuity of the operator $\mathcal T$, it is sufficient to show that $(p_1,p_2)$ is the unique adherent value of the sequence $(p_{1,n},p_{2,n})$, for that let us show $(p_1,p_2)$ is the unique solution of \eqref{nles1b}-\eqref{nles2b} by passing the limit in \eqref{nles1bn}-\eqref{nles2bn}. Passing to the limit in \eqref{nles1bn}: \begin{align*} &\int_{\Omega}\phi\frac{\rho_1 (\overline{p}_{1,n})Z(\overline{s}_{1,n})- \rho_1^\star s_1^\star}{h} \varphi \, dx +\int_{\Omega}\mathbf{K} M_1^\epsilon(\overline{s}_{1,n})\rho_1(\overline{p}_{1,n}) { \nabla} p_{1,n} \cdot\nabla \varphi \,dx\\ &-\int_{\Omega}\mathbf{K} M_1(\overline{s}_{1,n})\rho_1^2(\overline{p}_{1,n})\mathbf{g} \cdot\nabla \varphi \,dx +\eta \int_{\Omega}\rho_1(\overline{p}_{1,n})\nabla (\mathcal{P}_N \overline{p}_{1,n} - \mathcal{P}_N \overline{p}_{2,n}) \cdot \nabla \varphi \,dx\\ &+\int_{\Omega}\rho_1(\overline{p}_{1,n})Z(\overline{s}_{1,n}) f_{P}\varphi \,dx =\int_{\Omega}\rho_1(\overline{p}_{1,n})s_1^If_{I}\varphi \,dx, \end{align*} where $\overline{s}_{1,n}=\bar{f}^{-1}(\overline{p}_{1,n}-\overline{p}_{2,n})$. The passage to the limit in the first term is due to the continuity of $Z,$ $\bar{f}^{-1}$ and $\rho_1,$ the convergence \eqref{cvp1ns} and \eqref{cvp2ns}, and the domination of $\rho_1 (\overline{p}_{1,n})Z(\overline{s}_{1,n})\varphi$ by $\rho_M|\varphi|$, which allow us to apply the Lebesgue theorem. The second term is treated as follows, the sequence $\bigl (KM_1^\varepsilon( {\overline{s}_1}_{n})\rho_1({\overline{p}_1}_n) \nabla\varphi\bigr )_n$ is dominated and converges a.e. as $n$ goes to infinity. Then, by Lebesgue theorem, we have the following strong convergence in $L^2(\Omega)$, $$KM_1^\varepsilon(\overline{s}_{1,n})\rho_1({\overline{p}_1}_n) \nabla \varphi \to KM_1^\varepsilon( \overline{s}_1)\rho_1(\overline{p}_1) \nabla \varphi. \label{cvterm2}$$ Furthermore and due to the convergence \eqref{cvp1nw}, it follows that $$\nabla p_1^n \to \nabla p_1 \quad \text{weakly in } L^2(\Omega), \label{cvterm2bis}$$ then, the convergence \eqref{cvterm2} and \eqref{cvterm2bis} establish the limit for the second term. The fourth term $$\eta \int_{\Omega}\rho_1(\overline{p}_{1,n})\nabla (\mathcal{P}_N \overline{p}_{1,n} - \mathcal{P}_N \overline{p}_{2,n}) \cdot \nabla \varphi \,dx,$$ is treated as follows, $$\rho_i(\overline{p}_{i,n}) \nabla \varphi \to \rho_i(\overline{p}_i) \nabla \varphi \quad \text{strongly in } (L^2(\Omega))^d \; (i=1,2). \label{cvterm2b}$$ Furthermore $\overline{p}_{i,n}$ converges in $L^2(\Omega)$, it follows that $$\nabla \mathcal{P}_N\overline{p}_{i,n} \rightharpoonup\nabla \mathcal{P}_N\overline{p}_i \quad \text{weakly in } (L^2(\Omega))^d \; (i=1,2). \label{cvterm2bisb}$$ Then, the convergence \eqref{cvterm2b}-\eqref{cvterm2bisb} allow us to pass the limit in the fourth term. The convergence of the other terms are also an application of the Lebesgue convergence theorem. The passage to the limit on \eqref{nles2bn} is obtained in the same manner. Thus $(p_1,p_2)$ is a solution of \eqref{nles1b}-\eqref{nles2b}, which establishes the continuity and completes the proof of the lemma. \end{proof} \begin{lemma}[A priori estimate] \label{apriori} There exists a positive constant $r$ such that, if $(p_1,p_2)=\lambda \mathcal{T}(p_1,p_2)$ with $\lambda\in (0,1)$, then $$\Vert (p_1,p_2)\Vert_{L^2(\Omega)\times L^2(\Omega)} \le r,$$ where $r$ is independent of $\lambda$. \end{lemma} \begin{proof} Assume $(p_1,p_2)=\lambda \mathcal{T}(p_1,p_2)$ holds, recall that $s_1=\overline{f}^{-1}(p_1-p_2)$ and $s_2=1-s_1$, then $(p_1,p_2)$ satisfies \label{apriori1} \begin{aligned} &\int_{\Omega}\mathbf{K} M_1^\epsilon(s_1)\rho_1(p_1) { \nabla} p_1 \cdot\nabla \varphi \,dx\\ &=-\lambda \int_{\Omega}\phi\frac{\rho_1 (p_1)Z(s_1) - \rho_1^\star s_1^\star}{h} \varphi \, dx +\lambda \int_{\Omega}\mathbf{K} M_1(s_1)\rho_1^2(p_1)\mathbf{g} \cdot\nabla \varphi \,dx \\ &\quad -\lambda\int_{\Omega}\rho_1(p_1)Z(s_1) f_{P}\varphi \,dx +\lambda\int_{\Omega}\rho_1(p_1)s_1^If_{I}\varphi \,dx\\ &\quad -\lambda \eta \int_{\Omega}\rho_1(p_1)\nabla(\mathcal{P}_N p_1- \mathcal{P}_N p_2) \cdot \nabla \varphi \,dx, \end{aligned} \label{apriori2} \begin{aligned} &\int_{\Omega}\mathbf{K} M_2^\epsilon(s_2)\rho_2(p_2) { \nabla} p_2 \cdot\nabla \xi \,dx\\ &=-\lambda \int_{\Omega}\phi\frac{\rho_2 (p_2)Z(s_2) - \rho_2^\star s_2^\star}{h} \xi \, dx +\lambda \int_{\Omega}\mathbf{K} M_2(s_2)\rho_2^2(p_2)\mathbf{g} \cdot\nabla \xi \,dx \\ &\quad -\lambda \int_{\Omega}\rho_2(p_2)Z(s_2)f_{P}\xi \,dx +\lambda \int_{\Omega}\rho_2(p_2)s_2^If_{I}\xi \,dx\\ &\quad +\lambda\eta \int_{\Omega}\rho_2(p_2)\nabla (\mathcal{P}_N p_1- \mathcal{P}_N p_2) \cdot \nabla \xi \,dx. \end{aligned} for all $(\varphi,\xi)$ belonging to $H^1_{\Gamma_1}(\Omega)\times H^1_{\Gamma_1}(\Omega)$. Consider $\varphi=g_1(p_1):=\int_0^{p_1}\frac{1}{\rho_1(\zeta)} \,d\zeta\in H^1_{\Gamma_1}(\Omega)$ in \eqref{apriori1} and $\xi=g_2(p_2):=\int_0^{p_2}\frac{1}{\rho_2(\zeta)}\,d\zeta \in H^1_{\Gamma_1}(\Omega)$ in \eqref{apriori2}. Summing up these quantities, we obtain \begin{aligned} &\lambda \int_\Omega \frac \phi h \Bigl( \big(\rho_1(p_1)Z(s_1) -\rho ^\star_1 s^\star_1\big) g_1(p_1) + \big(\rho_2(p_2)Z(s_2) -\rho ^\star_2 s^\star_2\big)g_2(p_2) \Bigr) \,dx\\ &+\int_\Omega \mathbf{K}M_1^\varepsilon \nabla p_1\cdot \nabla p_1\,dx +\lambda\eta \int_{\Omega}\nabla (\mathcal{P}_N p_1- \mathcal{P}_N p_2) \cdot \nabla(p_1-p_2) \,dx \\ &-\lambda \int_\Omega \mathbf{K}\rho_1(p_1)M_1(s_1)g\cdot\nabla p_1\,dx +\int_\Omega \mathbf{K}M_2^\varepsilon \nabla p_2\cdot \nabla p_2\,dx\\ &-\lambda\int_\Omega \mathbf{K}\rho_2(p_2)M_2(s_2)g\cdot\nabla p_2\,dx\\ &+\lambda \int_{\Omega}\big(\rho_1(p_1)Z(s_1) g_1(p_1)+ \rho_2(p_2)Z(s_2) g_2(p_2)\big)f_{P}\,dx\\ &= \lambda \int_{\Omega}\big((\rho_1(p_1)s_1^Ig_1(p_1) + \rho_2(p_2)s_2^I g_2(p_2)\big)f_{I}\,dx. \end{aligned}\label{iiplus} Remark that the functions $p_i\to g_i(p_i)$ is sub-linear, we deduce from Cauchy-Schwarz and Poincar\'e inequalities that \eqref{iiplus} reduces to \label{iplus} \begin{aligned} &\epsilon \int_\Omega \vert{ \nabla} p_1\vert^2 \,dx +\epsilon \int_\Omega \vert{ \nabla} p_2\vert^2 \,dx +\lambda \eta\int_\Omega \vert{ \nabla}(\mathcal{P}_N p_1- \mathcal{P}_N p_2)\vert^2 \,dx\\ &\le C_1 (1+\Vert f_{P} \Vert^2_{L^2(\Omega)}+\Vert f_{I} \Vert^2_{L^2(\Omega)}+\Vert \rho ^\star_1 s^\star_1 \Vert^2_{L^2(\Omega)} + \Vert \rho ^\star_2 s^\star_2 \Vert^2_{L^2(\Omega)}), \end{aligned} where $C_1$ depends on $\epsilon$ and not on $\lambda$. It is important in \eqref{i+N} to ensure that $C_1$ does not depend on $N$. Lemma \ref{concomp}, Lemma \ref{apriori} allow to apply the Leray-Schauder fixed point theorem \cite{Zeidler}, thus the proof of proposition \ref{pleray-sch} is completed. \end{proof} \noindent{\bf Step 2.} Now we are concerned with the limit $N$ goes to infinity (we omit the dependence of solutions on $\epsilon$). For all $N$, we have established a solution $(p_{1,N},p_{2,N})\in H^1_{\Gamma_1}(\Omega)\times H^1_{\Gamma_1}(\Omega)$ to \eqref{ellvNeps1} \eqref{ellvNeps2} satisfying \begin{align}\label{ellvN1} &\int_\Omega \phi\frac{ \rho_1(p_1^{N})Z(s_1^{N})-\rho ^\star_1 s^\star_1}{h} \varphi \,dx +\int_{\Omega}\mathbf{K} M_1^\epsilon(s_1^{N} )\rho_1(p_1^{N} ) { \nabla} p_1^{N} \cdot\nabla \varphi \,dx \nonumber\\ &-\int_{\Omega}\mathbf{K} M_1(s_1^{N} )\rho_1^2 (p_1^{N} )\mathbf{g} \cdot\nabla \varphi \,dx +\eta \int_{\Omega}\rho_1(p_1^{N})\nabla (\mathcal{P}_N p_1^{N} - \mathcal{P}_N p_2^{N}) \cdot \nabla \varphi \,dx \nonumber\\ &+\int_{\Omega}\rho_1(p_1^{N})Z(s_1^{N}) f_{P}\varphi \,dx =\int_{\Omega}\rho_1(p_1^{N})s_1^I f_{I}\varphi \,dx, \end{align} \begin{align}\label{ellvN2} &\int_\Omega \phi\frac{ \rho_2(p_2^{N})Z(s_2^{N})-\rho ^\star_2 s^\star_2}{h} \xi \,dx +\int_{\Omega}\mathbf{K} M_2^\epsilon(s_2^{N} )\rho_2(p_2^{N}) { \nabla} p_2^{N} \cdot\nabla \xi \,dx \nonumber\\ &-\int_{\Omega}\mathbf{K} M_2(s_2^{N} )\rho_2^2 (p_2^{N})\mathbf{g} \cdot\nabla \xi \,dx -\eta \int_{\Omega}\rho_2(p_2^{N})\nabla (\mathcal{P}_N p_1^{N} -\mathcal{P}_N p_2^{N}) \cdot \nabla \xi \,dx \nonumber\\ &+\int_{\Omega}\rho_2(p_2^{N})Z(s_2^{N}) f_{P}\xi \,dx =\int_{\Omega}\rho_2(p_2^{N})s_2^I f_{I}\xi \,dx, \end{align} for all $(\varphi,\xi)$ belonging to $H^1_{\Gamma_1}(\Omega)\times H^1_{\Gamma_1}(\Omega)$. Reproducing the estimate \eqref{iplus} with $\lambda=1$, we get \label{i+N} \begin{aligned} &\epsilon \int_\Omega \vert{ \nabla} p_1\vert^2 \,dx +\epsilon \int_\Omega \vert{ \nabla} p_2\vert^2 \,dx + \eta\int_\Omega \vert{ \nabla}(\mathcal{P}_N p_1- \mathcal{P}_N p_2)\vert^2 \,dx\\ &\le C_1 (1+\Vert f_{P} \Vert^2_{L^2(\Omega)}+\Vert f_{I} \Vert^2_{L^2(\Omega)}+\Vert \rho ^\star_1 s^\star_1 \Vert^2_{L^2(\Omega)} + \Vert \rho ^\star_2 s^\star_2 \Vert^2_{L^2(\Omega)}), \end{aligned} where $C_1$ depends on $\epsilon$ and not on $N$. Then, up to a subsequence, we have the convergence, \begin{align} p_{1,N} \to p_1 \quad\text{weakly in } H^1_{\Gamma_1}(\Omega),\text{ strongly in } L^2 (\Omega)\text{ and a.e. in }\Omega\\ p_{2,N} \to p_2 \quad\text{weakly in } H^1_{\Gamma_1}(\Omega),\text{ strongly in } L^2 (\Omega)\text{ and a.e. in }\Omega. \end{align} The convergence in \eqref{ellvN1}-\eqref{ellvN2} with respect to $N$ are obtained in the same manner as for the convergence with respect to $n$ in \eqref{nles1bn} \eqref{nles2bn}. \medskip \noindent{\bf Step 3.} Passage to the limit as $\epsilon$ approaches zero. For all $\epsilon >0$, we have shown that there exists $(p_{1,\epsilon},p_{2,\epsilon})\in H^1_{\Gamma_1}(\Omega)\times H^1_{\Gamma_1}(\Omega)$, satisfying \label{ellveps1} \begin{aligned} &\int_\Omega \phi\frac{ \rho_1(p_1^{\epsilon})Z(s_1^{\epsilon})-\rho ^\star_1 s^\star_1}{h} \varphi \,dx +\int_{\Omega}\mathbf{K} M_1^\epsilon(s_1^{\epsilon} )\rho_1(p_1^{\epsilon} ) { \nabla} p_1^{\epsilon} \cdot\nabla \varphi \,dx \\ &-\int_{\Omega}\mathbf{K} M_1(s_1^{\epsilon} )\rho_1^2 (p_1^{\epsilon} )\mathbf{g} \cdot\nabla \varphi \,dx +\eta \int_{\Omega}\rho_1(p_1^{\epsilon})\nabla ( p_1^{\epsilon} - p_2^{\epsilon}) \cdot \nabla \varphi \,dx \\ &+\int_{\Omega}\rho_1(p_1^{\epsilon})Z(s_1^{\epsilon}) f_{P}\varphi \,dx =\int_{\Omega}\rho_1(p_1^{\epsilon})s_1^I f_{I}\varphi \,dx, \end{aligned} \label{ellveps2} \begin{aligned} &\int_\Omega \phi\frac{ \rho_2(p_2^{\epsilon})Z(s_2^{\epsilon})-\rho ^\star_2 s^\star_2}{h} \xi \,dx +\int_{\Omega}\mathbf{K} M_2^\epsilon(s_2^{\epsilon} )\rho_2(p_2^{\epsilon}) { \nabla} p_2^{\epsilon} \cdot\nabla \xi \,dx \\ &-\int_{\Omega}\mathbf{K} M_2(s_2^{\epsilon} )\rho_2^2 (p_2^{\epsilon})\mathbf{g} \cdot\nabla \xi \,dx -\eta \int_{\Omega}\rho_2(p_2^{\epsilon})\nabla (p_1^{\epsilon} - p_2^{\epsilon}) \cdot \nabla \xi \,dx\\ &+\int_{\Omega}\rho_2(p_2^{\epsilon})Z(s_2^{\epsilon}) f_{P}\xi \,dx =\int_{\Omega}\rho_2(p_2^{\epsilon})s_2^I f_{I}\xi \,dx, \end{aligned} for all $(\varphi,\xi)$ belonging to $H^1_{\Gamma_1}(\Omega)\times H^1_{\Gamma_1}(\Omega)$, with $s_1^{\epsilon}=\overline{f}^{-1}(p_1^{\epsilon}-p_2^{\epsilon})$ and $s_2^{\epsilon}=1-s_1^{\epsilon}$. We need uniform estimates on the solutions independent of the regularization $\epsilon$ in order to pass to the limit in $\epsilon$. For that, we are going to use the feature of global pressure. After the passage to the limit in $\epsilon$, a maximum principle on saturations is possible. We are now concerned with a uniform estimate on the gradient of $\beta ({s_1^\epsilon})$, and on the global pressure $p^\epsilon$. We state the following two lemmas. \begin{lemma}\label{lemunifeps} The sequences $({s_i^\epsilon})_\epsilon$, $(p^{\epsilon}:=p_2^{\epsilon}+\tilde{p}({s_1^\epsilon}))_\epsilon$ defined by Proposition \ref{pleray-sch} satisfy \begin{gather} (p^{\epsilon})_\epsilon \text{ is uniformly bounded in } H^1_{\Gamma _1}({\Omega})\label{ppc}\\ (\sqrt{\epsilon}\,\, \nabla{p_i^\epsilon})_\epsilon \text{is uniformly bounded in } L^2(\Omega)\label{epsestimatec}\\ (\beta({s_1^\epsilon}))_\epsilon \text{ is uniformly bounded in } H^1({\Omega})\label{betabetac}\\ \nabla f({s_1^\epsilon}))_\epsilon \text{ is uniformly bounded in } L^2({\Omega})\label{ffc} \end{gather} \end{lemma} \begin{proof} Consider $\varphi=g_1(p_1^\epsilon):=\int_0^{p_1^\epsilon}\frac{1}{\rho_1(\zeta)} \,d\zeta\in H^1_{\Gamma_1}(\Omega)$ in \eqref{ellveps1} and $\xi=g_2(p_2^\epsilon):=\int_0^{p_2^\epsilon}\frac{1}{\rho_2(\zeta)}\,d\zeta \in H^1_{\Gamma_1}(\Omega)$ in \eqref{ellveps2}. Summing these quantities, we obtain \begin{align*} & \int_\Omega \frac \phi h \Bigl( \big(\rho_1(p_1^\epsilon)Z(s_1^\epsilon)-\rho ^\star_1 s^\star_1\big) g_1(p_1^\epsilon) + \big(\rho_2(p_2^\epsilon)Z(s_2^\epsilon)-\rho ^\star_2 s^\star_2\big)g_2(p_2^\epsilon) \Bigr) \,dx \\ &+\int_\Omega \mathbf{K}M_1^\varepsilon \nabla p_1^\epsilon\cdot \nabla p_1^\epsilon\,dx +\eta \int_{\Omega}\nabla ( p_1^\epsilon- p_2^\epsilon) \cdot \nabla(p_1^\epsilon-p_2^\epsilon) \,dx\\ &- \int_\Omega \mathbf{K}\rho_1(p_1^\epsilon)M_1(s_1^\epsilon)\mathbf{g}\cdot\nabla p_1^\epsilon\,dx +\int_\Omega \mathbf{K}M_2^\varepsilon \nabla p_2^\epsilon\cdot \nabla p_2^\epsilon\,dx\\ &-\int_\Omega \mathbf{K}\rho_2(p_2^\epsilon)M_2(s_2^\epsilon)\mathbf{g}\cdot\nabla p_2^\epsilon\,dx \\ &+\int_{\Omega}\big(\rho_1(p_1^\epsilon)Z(s_1^\epsilon) g_1(p_1^\epsilon)+ \rho_2(p_2^\epsilon)Z(s_2^\epsilon) g_2(p_2^\epsilon)\big)f_{P}\,dx \\ &= \int_{\Omega}\big((\rho_1(p_1^\epsilon)s_1^Ig_1(p_1^\epsilon)+ \rho_2(p_2^\epsilon)s_2^I g_2(p_2^\epsilon)\big)f_{I}\,dx, %\label{unifeps} \end{align*} then \label{unifepsf} \begin{aligned} &\int_\Omega \mathbf{K}M_1 \nabla p_1^\epsilon\cdot \nabla p_1^\epsilon\,dx +\int_\Omega \mathbf{K}M_2 \nabla p_2^\epsilon\cdot \nabla p_2^\epsilon\,dx\\ &+\epsilon \int_\Omega \mathbf{K} \nabla p_1^\epsilon\cdot \nabla p_1^\epsilon\,dx +\epsilon \int_\Omega \mathbf{K} \nabla p_2^\epsilon\cdot \nabla p_2^\epsilon\,dx +\eta \int_{\Omega}\nabla f(s_1) \cdot \nabla f(s_1) \,dx \\ &= \int_\Omega \mathbf{K}\rho_1(p_1^\epsilon)M_1(s_1^\epsilon)\mathbf{g}\cdot\nabla p_1^\epsilon\,dx +\int_\Omega \mathbf{K}\rho_2(p_2^\epsilon)M_2(s_2^\epsilon) \mathbf{g}\cdot\nabla p_2^\epsilon\,dx \\ &\quad-\int_{\Omega}\big(\rho_1(p_1^\epsilon)Z(s_1^\epsilon) g_1(p_1^\epsilon)+ \rho_2(p_2^\epsilon)Z(s_2^\epsilon) g_2(p_2^\epsilon)\big)f_{P}\,dx\\ &\quad +\int_{\Omega}\big((\rho_1(p_1^\epsilon)s_1^Ig_1(p_1^\epsilon) + \rho_2(p_2^\epsilon)s_2^I g_2(p_2^\epsilon)\big)f_{I}\,dx\\ &\quad - \int_\Omega \frac \phi h \Bigl( \big(\rho_1(p_1^\epsilon) Z(s_1^\epsilon)-\rho ^\star_1 s^\star_1\big) g_1(p_1^\epsilon) + \big(\rho_2(p_2^\epsilon) Z(s_2^\epsilon)-\rho ^\star_2 s^\star_2\big)g_2(p_2^\epsilon) \Bigr) \,dx . \end{aligned} The hypothesis $(H2)$ and with the help of Cauchy-Schwarz inequality, we have \begin{gather*} \big|\int_\Omega \mathbf{K}\rho_1(p_1^\epsilon)M_1(s_1^\epsilon)\mathbf{g} \cdot\nabla p_1^\epsilon\,dx\big| \le C+\frac{k_0}{2}\int_\Omega M_1 \nabla p_1^\epsilon\cdot \nabla p_1^\epsilon\,dx,\\ \big|\int_\Omega \mathbf{K}\rho_2(p_2^\epsilon)M_2(s_2^\epsilon)\mathbf{g} \cdot\nabla p_2^\epsilon\,dx\big| \le C+\frac{k_0}{2}\int_\Omega M_2 \nabla p_2^\epsilon\cdot \nabla p_2^\epsilon\,dx, \end{gather*} then the gravity terms in \eqref{unifepsf} on the right hand side are absorbed by pressures dissipative terms. Recall that, the functions $p_i\to g_i(p_i)$ is sub-linear (i.e $|g_i(p_i)|\le\frac{1}{\rho_m}|p_i|$ ), then from \eqref{unifepsf}, one obtains \begin{aligned} &\int_\Omega M_1({s_1^\epsilon}) |\nabla {p_1^\epsilon}|^2 \, dx +\int_\Omega M_2({s_2^\epsilon}) |\nabla {p_2^\epsilon}|^2 \,dx +\eta \int_{\Omega}|\nabla f(s_1)|^2 \,dx\\ &+\epsilon\int_\Omega |\nabla {p_1^\epsilon}|^2 \, dx +\epsilon\int_\Omega |\nabla {p_2^\epsilon}|^2 \,dx\\ &\le C(1+\| {p_1^\epsilon}\|_{L^2(\Omega)}+\| {p_2^\epsilon}\|_{L^2(\Omega)}). \end{aligned} \label{estip1p2} Return now to the relationship between pressures and global pressure. From \eqref{defp}, we have $p^\epsilon=p_2^\epsilon+\tilde{p}(s_1^\epsilon)=p_1^\epsilon+\bar{p}(s_1^\epsilon)$, and \begin{align}\label{graddefp} \nabla p^\epsilon=\nabla {p_2^\epsilon} + \frac{M_1({s_1^\epsilon})}{M({s_1^\epsilon})} \nabla f({s_1^\epsilon}) =\nabla {p_1^\epsilon} - \frac{M_2({s_2^\epsilon})}{M({s_1^\epsilon})} \nabla f({s_1^\epsilon}), \end{align} which imply that \label{magic1eps} \begin{aligned} &\int_{\Omega} M({s_1^\epsilon}) |\nabla p^\epsilon|^2\,dx + \int_{\Omega} \frac{M_1({s_1^\epsilon}) M_2({s_2^\epsilon})}{M({s_1^\epsilon})} |\nabla f({s_1^\epsilon})|^2 \, dx\\ &=\int_{\Omega} M_1({s_1^\epsilon}) \nabla {p_1^\epsilon} \cdot \nabla {p_1^\epsilon} \, dx +\int_{\Omega} M_2({s_2^\epsilon}) \nabla {p_2^\epsilon} \cdot \nabla {p_2^\epsilon} \,dx. \end{aligned} The estimate \eqref{estip1p2} is equivalent to \begin{align*} &\int_{\Omega} M({s_1^\epsilon}) |\nabla p^\epsilon|^2\,dx + \int_{\Omega} \frac{M_1({s_1^\epsilon}) M_2({s_2^\epsilon})}{M({s_1^\epsilon})} |\nabla f({s_1^\epsilon})|^2 \, dx\\ &+\eta \int_{\Omega}|\nabla f(s_1)|^2 \,dx+ \epsilon\int_\Omega |\nabla {p_1^\epsilon}|^2 \, dx +\epsilon\int_\Omega |\nabla {p_2^\epsilon}|^2 \,dx\\ &\le C(1+\| {p_1^\epsilon}\|_{L^2(\Omega)}+\| {p_2^\epsilon}\|_{L^2(\Omega)})\\ &\le C(1+\| p^\epsilon\|_{L^2(\Omega)}+\| \overline{p}({s_1^\epsilon})\|_{L^2(\Omega)} +\| \tilde{p}({s_1^\epsilon})\|_{L^2(\Omega)})\\ &\le C(1+\| \nabla p^\epsilon\|_{L^2(\Omega)}+\| \overline{p}({s_1^\epsilon})\|_{L^2(\Omega)} +\| \tilde{p}({s_1^\epsilon})\|_{L^2(\Omega)}), \end{align*} due to the Poincar\'e's inequality. Finally, using the fact that the function $\tilde{p}$ and $\bar{p}$ are bounded, and the global pressure term on the right hand side in the above inequality can be absorbed by the dissipative term in global pressure, on the left hand side, we deduce that there exists a constant $C_1$ independent of $\epsilon$, $C_1=C_1(h, \rho_m, M_i, \mathbf{g},f_p,f_I,s_1^I,s_2^I,\rho ^\star_1 s^\star_1,\rho ^\star_2s^\star_2,h,\phi,k_\infty,k_0)$ such that \label{pffp} \begin{aligned} &\int_{\Omega} M({s_1^\epsilon}) |\nabla p^\epsilon|^2\,dx + \int_{\Omega} \frac{M_1({s_1^\epsilon}) M_2({s_2^\epsilon})}{M({s_1^\epsilon})} |\nabla f({s_1^\epsilon})|^2 \, dx\\ &+\eta \int_{\Omega}|\nabla f(s_1)|^2 \,dx+ \epsilon\int_\Omega |\nabla {p_1^\epsilon}|^2 \, dx +\epsilon\int_\Omega |\nabla {p_2^\epsilon}|^2 \,dx \le C_1, \end{aligned} which establish the estimates \eqref{ppc}, \eqref{epsestimatec} and \eqref{ffc}. For the estimate \eqref{betabetac}, we use the fact that the second term on the left hand side in \eqref{pffp} is bounded and the total mobility is bounded below due to the assumption (H3), we have $$\int_{\Omega} M_1({s_1^\epsilon}) M_2({s_2^\epsilon})|\nabla f({s_1^\epsilon})|^2 \, dx \le m_0C_1,$$ which implies \begin{align*} \int_{\Omega} |\nabla \overline{\beta} ({s_1^\epsilon}) |^2\, dx &=\int_{\Omega} \frac{M_1^2({s_1^\epsilon}) M_2^2({s_2^\epsilon})}{M^2({s_1^\epsilon})} |\nabla f({s_1^\epsilon})|^2 \, dx\\ &\le \int_{\Omega} M_1({s_1^\epsilon}) M_2({s_2^\epsilon})|\nabla f({s_1^\epsilon})|^2 \, dx \le m_0C_1, \end{align*} and completes the proof of lemma. \end{proof} From the previous lemma, we deduce the following convergence. \begin{lemma}[Strong and weak convergence] Up to a subsequence the sequences $({s_i^\epsilon})_\epsilon$, $(p^\epsilon)_\epsilon$, $({p_i^\epsilon})_\epsilon$ have the following convergence \begin{gather} p^\epsilon \rightharpoonup p \quad\text{weakly in } H^1_{\Gamma_1}(\Omega) \label{ph11c}\\ \beta({s_1^\epsilon}) \rightharpoonup \beta(s_1) \quad\text{weakly in } H^1(\Omega),\label{betah11c}\\ p^\epsilon \to p \quad \text{almost everywhere in } \Omega \label{ppppc}\\ \beta({s_1^\epsilon}) \to \beta(s_1) \quad\text{ almost everywhere in } \Omega \label{betaahc}\\ Z({s_1^\epsilon}) \to Z(s_1)\quad\text{almost everywhere in } \Omega \label{beta1123c}\\ Z({s_1^\epsilon}) \to Z(s_1) \quad\text{strongly in } L^2(\Omega) \label{beta121c}\\ {p_i^\epsilon} \to p_i \quad \text{almost everywhere in } \Omega \label{ppppic}. \end{gather} \end{lemma} \begin{proof} The weak convergence \eqref{ph11c}--\eqref{betah11c} follows from the uniform estimates \eqref{ppc} and \eqref{betabetac} of lemma \ref{lemunifeps}, while \begin{gather*} p^\epsilon\to p \quad \text{strongly in $L^2(\Omega)$ and a.e. in } \Omega,\\ \beta({s_1^\epsilon})\to \beta^\star \quad\text{strongly in $L^2(\Omega)$ and a. e. in } \Omega \end{gather*} is due to the compact injection of $H^1_{\Gamma_1}$ in to $L^2(\Omega)$. As $\beta(s_1):=\beta(Z(s_1))$ and $\beta ^{-1}$ is continuous, $$Z({s_1^\epsilon}) \to Z(s_1) \text{ a. e. in } \Omega,$$ while the Lebesgue theorem ensures the strong convergence \eqref{beta121c}. The convergence \eqref{ppppic} is a consequence of \eqref{ppppc}--\eqref{beta1123c} and the fact that $p_1^\epsilon:=p^\epsilon-\bar{p}(Z(s_1^\epsilon))$, $p_2^\epsilon:=p^\epsilon-\tilde{p}(Z(s_1^\epsilon))$. \end{proof} To achieve the proof of Theorem \ref{deg-ell-t}, it remains to pass to the limit as $\epsilon$ goes to zero in the formulations \eqref{ellveps1}\eqref{ellveps2} and a proof of a maximum principle on saturations. For all test functions $(\varphi, \xi) \in H^1_{\Gamma_1}(\Omega) \times H^1_{\Gamma_1}(\Omega)$, \begin{align*}%\label{ellveps1} &\int_\Omega \phi\frac{ \rho_1(p_1^{\epsilon})Z(s_1^{\epsilon})-\rho ^\star_1 s^\star_1}{h} \varphi \,dx +\int_{\Omega}\mathbf{K} M_1^\epsilon(s_1^{\epsilon} )\rho_1(p_1^{\epsilon} ) { \nabla} p_1^{\epsilon} \cdot\nabla \varphi \,dx\\ &-\int_{\Omega}\mathbf{K} M_1(s_1^{\epsilon} )\rho_1^2 (p_1^{\epsilon} )\mathbf{g} \cdot\nabla \varphi \,dx +\eta \int_{\Omega}\rho_1(p_1^{\epsilon})\nabla ( p_1^{\epsilon} - p_2^{\epsilon}) \cdot \nabla \varphi \,dx\\ &+\int_{\Omega}\rho_1(p_1^{\epsilon})Z(s_1^{\epsilon}) f_{P}\varphi \,dx =\int_{\Omega}\rho_1(p_1^{\epsilon})s_1^I f_{I}\varphi \,dx, \end{align*} \begin{align*}%\label{ellveps2} &\int_\Omega \phi\frac{ \rho_2(p_2^{\epsilon})Z(s_2^{\epsilon})-\rho ^\star_2 s^\star_2}{h} \xi \,dx +\int_{\Omega}\mathbf{K} M_2^\epsilon(s_2^{\epsilon} )\rho_2(p_2^{\epsilon}) { \nabla} p_2^{\epsilon} \cdot\nabla \xi \,dx\\ &-\int_{\Omega}\mathbf{K} M_2(s_2^{\epsilon} )\rho_2^2 (p_2^{\epsilon})\mathbf{g} \cdot\nabla \xi \,dx -\eta \int_{\Omega}\rho_2(p_2^{\epsilon})\nabla (p_1^{\epsilon} - p_2^{\epsilon}) \cdot \nabla \xi \,dx\\ &+\int_{\Omega}\rho_2(p_2^{\epsilon})Z(s_2^{\epsilon}) f_{P}\xi \,dx =\int_{\Omega}\rho_2(p_2^{\epsilon})s_2^I f_{I}\xi \,dx, \end{align*} The first terms of the above equalities converge due to the strong convergence of $\rho_i({p_i^\epsilon}) Z({s_i^\epsilon})$ to $\rho_i(p_i) Z(s_i)$ in $L^2(\Omega)$. The second terms can be written as, \label{lovec} \begin{aligned} &\int_{\Omega}\mathbf{K}M_i^\epsilon({s_i^\epsilon})\rho_i({p_i^\epsilon}) { \nabla} {p_i^\epsilon} \cdot\nabla \varphi \,dx\\ &= \int_{\Omega}\mathbf{K}M_i({s_i^\epsilon})\rho_i({p_i^\epsilon}) { \nabla}p^\epsilon \cdot\nabla \varphi \,dx\\ &\quad + \int_{\Omega}\mathbf{K}\rho_i({p_i^\epsilon}) { \nabla} \beta({s_i^\epsilon}) \cdot\nabla \varphi \,dx +\sqrt{\varepsilon} \int_{\Omega}\mathbf{K}\rho_i({p_i^\epsilon}) (\sqrt{\varepsilon}\, \, { \nabla} {p_i^\epsilon}) \cdot\nabla \varphi \,dx. \end{aligned} The first two terms on the right hand side of the equation converge arguing in two steps. Firstly, the Lebsgue theorem and the convergence \eqref{beta1123c}\eqref{ppppic} establish \begin{gather*} \rho_i({p_i^\epsilon}) M_i({s_i^\epsilon})\nabla \varphi \to \rho_i(p_i)M_i(s_i)\nabla \varphi \quad \text{strongly in } (L^2(Q_T))^d, \\ \rho_i({p_i^\epsilon}) \nabla \varphi \to \rho_i(p_i)\nabla \varphi \quad \text{strongly in } (L^2(Q_T))^d. \end{gather*} Secondly, the weak convergence on pressure \eqref{ph11c} combined to the above strong convergence validate the convergence for the first term of the right hand side of \eqref{lovec}, and the weak convergence \eqref{betah11c} combined to the above strong convergence validate the convergence for the second term of the right hand side of \eqref{lovec}. The third term converges to zero due to the uniform estimate \eqref{epsestimatec}, and this achieves the passage to the limit on the second terms. The convergence of the fourth terms of the above equations are due to the uniform estimate \eqref{ffc}. The other terms converge using \eqref{beta1123c}\eqref{ppppic} and the Lebesgue dominated convergence theorem. In summarize, we have shown, there exists $(p_1^{h},p_1^{h})\in H^1_{\Gamma_1}(\Omega)\times H^1_{\Gamma_1}(\Omega)$ solution of \label{maxv1} \begin{aligned} &\int_\Omega \phi\frac{ \rho_1({p_1^h})Z({s_1^h})-\rho ^\star_1 s^\star_1}{h} \varphi \,dx +\int_{\Omega}\mathbf{K} M_1({s_1^h} )\rho_1({p_1^h} ) { \nabla} {p_1^h} \cdot\nabla \varphi \,dx \\ &-\int_{\Omega}\mathbf{K} M_1({s_1^h} )\rho_1^2 ({p_1^h} )\mathbf{g} \cdot\nabla \varphi \,dx +\eta \int_{\Omega}\rho_1(p_1^{h})\nabla ( p_1^{h} - p_2^{h}) \cdot \nabla \varphi \,dx\\ &+\int_{\Omega}\rho_1({p_1^h})Z({s_1^h}) f_{P}\varphi \,dx =\int_{\Omega}\rho_1({p_1^h})s_1^I f_{I}\varphi \,dx, \end{aligned} \label{maxv2} \begin{aligned} &\int_\Omega \phi\frac{ \rho_2({p_2^h})Z({s_2^h})-\rho ^\star_2 s^\star_2}{h} \xi \,dx +\int_{\Omega}\mathbf{K} M_2({s_2^h} )\rho_2({p_2^h} ) { \nabla} {p_2^h} \cdot\nabla \xi \,dx\\ &-\int_{\Omega}\mathbf{K} M_2({s_2^h} )\rho_2^2 ({p_2^h})\mathbf{g} \cdot\nabla \xi \,dx -\eta \int_{\Omega}\rho_2(p_2^{h})\nabla ( p_1^{h} - p_2^{h}) \cdot \nabla \xi \,dx \\ &+\int_{\Omega}\rho_2({p_2^h})Z({s_2^h}) f_{P}\xi \,dx =\int_{\Omega}\rho_2({p_2^h})s_2^I f_{I}\xi \,dx, \end{aligned} for all $\varphi, \xi \in H^1_{\Gamma_1}(\Omega)$, with ${s_1^h}=\overline{f}^{-1}(p_1^{h}-p_2^{h})$ and ${s_2^h}=1-{s_1^h}$. \begin{lemma}[Maximum principle] \label{max} Under the conditions of Theorem \ref{deg-ell-t}, the saturation functions ${s_1^h}$ and ${s_2^h}$ which verify \eqref{maxv1}-\eqref{maxv2} are between zero and one a.e. in $\Omega$. \end{lemma} \begin{proof} It is sufficient to show that $s_i^h\ge 0$ a.e. in $\Omega$. For that, consider $\varphi=-(s_1)^-, \xi=-(s_2)^-$ respectively in \eqref{maxv1} and \eqref{maxv2} and by taking into consideration the definition of the map $Z$, and according to the extension of the mobility of each phase, $M_i(s_i^h)(s_i^h)^-=0$ $(i=1,2)$ we obtain \begin{gather*}%\label{maxv11} \int_\Omega \phi\frac{\rho ^\star_1 s^\star_1}{h} ({s_1^h})^- \,dx +\eta \int_\Omega \bar{f}' ({s_1^h})\nabla ({s_1^h})^- \cdot \nabla ({s_1^h})^- \,dx =-\int_{\Omega}\rho_1({p_1^h})s_1^I f_{I}(s_1^h)^- \,dx, \\ %\label{maxv22} \int_\Omega \phi\frac{ \rho ^\star_2 s^\star_2}{h} ({s_2^h})^- \,dx +\eta \int_\Omega \bar{f}' ({s_1^h})\nabla ({s_2^h})^- \cdot \nabla ({s_2^h})^- \,dx =-\int_{\Omega}\rho_2({p_2^h})s_2^I f_{I}(s_2^h)^- \,dx. \end{gather*} % Since it is possible to choose an extension $\bar{f}$ of $f$ outside $[0,1]$ in a way that ensures $\bar{f}' (s_1)$ different from zero outside $[0,1]$, we get $$\eta\int_\Omega |\nabla ({s_i^h})^-|^2\,dx\le 0 \quad (i=1,2)$$ which proves the maximum principle since $s_i^-$ vanishes on $\Gamma_1$, $i=1,2$. \end{proof} After this maximum principle, the weak formulations \eqref{degv-ell1} are established, and thus the theorem \ref{deg-ell-t} is then established. \section{Proof of Theorem \ref{non-degt}} \label{secnon-degt} The proof is based on a semi-discretization method in time \cite{A-AL83}. Let be $T>0$, $N\in \mathbb{N}^*$ and $h=\frac{T}{N}$. We define the following sequence parameterized by $h$: $p_{i,h}^0(x)=p_i^0(x)\quad \text{a.e. in} \Omega \quad i=1,2,$ for all $n\in [0,N-1]$, consider $(p_{1,h}^{n},\,p_{2,h}^{n}) \in L^2(\Omega)\times L^2(\Omega)$ with $\rho_1(p_{i,h}^{n}) s_{i,h}^{n}\geq 0$ for $i=1,2$, denote by $(f_{P})_h^{n+1}= \frac 1 h \int_{nh}^{(n+1)h} f_{P}(\tau)\,d\tau$, $(f_{I})_h^{n+1} = \frac 1 h \int_{nh}^{(n+1)h} f_{I}(\tau)\,d\tau$ and $(s_{i}^{I})_h^{n+1}= \frac 1 h \int_{nh}^{(n+1)h} s_{i}^{I}(\tau)\,d\tau$ for $i=1,2$, then define $(p_{1,h}^{n+1},\,p_{2,h}^{n+1})$ solution of \begin{aligned} &\phi\frac{ \rho_1(p_{1,h}^{n+1})s_{1,h}^{n+1}-\rho_1(p_{1,h}^{n}) s_{1,h}^{n}}{h} - \operatorname{div} (\mathbf{K}M_1 (s_{1,h}^{n+1})\rho_1(p_{1,h}^{n+1}) \nabla p_{1,h}^{n+1}) \\ &+ \operatorname{div} (\mathbf{K}\rho_1^2(p_{1,h}^{n+1}) M_1(s_{1,h}^{n+1}) \mathbf{g}) -\eta\operatorname{div} (\rho_1(p_{1,h}^{n+1})\nabla(p_{1,h}^{n+1}-p_{2,h}^{n+1}) ) \\ &+\rho_1(p_{1,h}^{n+1})s_{1,h}^{n+1}(f_{P})_h^{n+1}= \rho_1(p_{1,h}^{n+1})(s_{1}^{I})_h^{n+1}(f_{I})_h^{n+1}, \end{aligned}\label{gasdhn} \begin{aligned} &\phi\frac{ \rho_2(p_{2,h}^{n+1})s_{2,h}^{n+1}-\rho_2(p_{2,h}^{n}) s_{2,h}^{n}}{h} - \operatorname{div} (\mathbf{K}M_2 (s_{2,h}^{n+1})\rho_2(p_{2,h}^{n+1}) \nabla p_{2,h}^{n+1}) \\ &+ \operatorname{div} (\mathbf{K}\rho_2^2(p_{2,h}^{n+1}) M_2(s_{2,h}^{n+1}) \mathbf{g}) +\eta\operatorname{div} (\rho_2(p_{2,h}^{n+1})\nabla(p_{1,h}^{n+1}-p_{2,h}^{n+1}) ) \\ &+\rho_2(p_{2,h}^{n+1})s_{2,h}^{n+1}(f_{P})_h^{n+1}= \rho_2(p_{2,h}^{n+1})(s_{2}^{I})_h^{n+1}(f_{I})_h^{n+1}, \end{aligned}\label{waterdhn} with the boundary conditions \eqref{bcnondeg}. This sequence is well defined for all $n\in [0,N-1]$ by virtue of theorem \ref{deg-ell-t}. As a matter of fact, for given $s_{i,h}^{n}\rho_i(p_{i,h}^{n})\geq 0$ and $\rho_i(p_{i,h}^{n}) s_{i,h}^{n}\in L^2(\Omega)$, $i=1,2$, we construct $(p_{1,h}^{n+1},p_{2,h}^{n+1})\in H^1_{\Gamma_1}(\Omega)\times H^1_{\Gamma_1}(\Omega)$ so that $s_{i,h}^{n+1}\in [0,1]$. Now, we are concerned with uniform estimates with respect to $h$. We state the following lemma. \begin{lemma}[Uniform estimates with respect to $h$] \label{lemunifh} The solutions of \eqref{gasdhn}-\eqref{waterdhn} satisfy \begin{aligned} &\frac 1 h \int_\Omega \phi\Big( \mathcal{H}_1(p_{1,h}^{n+1})s_{1,h}^{n+1} - \mathcal{H}_1(p_{1,h}^{n})s_{1,h}^{n}\Big) \, dx\\ &+\frac 1 h \int_\Omega \phi\Big( \mathcal{H}_2(p_{2,h}^{n+1})s_{2,h}^{n+1} - \mathcal{H}_2(p_{2,h}^{n})s_{2,h}^{n}\Big) \, dx\\ &+\frac 1 h \int_\Omega \phi\Big( \mathcal{F}(s_{1,h}^{n+1})-\mathcal{F}(s_{1,h}^{n})\Big)\, dx + \eta \int_\Omega \vert{ \nabla} (p_{1,h}^{n+1}-p_{2,h}^{n+1})\vert^2 \,dx\\ &+k_0\int_\Omega M_1 (s_{1,h}^{n+1}) \nabla p_{1,h}^{n+1} \cdot \nabla p_{1,h}^{n+1} \, dx +k_0\int_\Omega M_2 (s_{2,h}^{n+1}) \nabla p_{2,h}^{n+1} \cdot \nabla p_{2,h}^{n+1} \, dx \\ &\le C(1+\Vert (f_{P})_h^{n+1} \Vert^2_{L^2(\Omega)} +\Vert (f_{I})_h^{n+1} \Vert^2_{L^2(\Omega)}) \end{aligned}\label{nablapn} where $C$ does not depend on $h$. For $i=1,2$, $\mathcal{H}_i(p_i):=\rho_i(p_i)g_i(p_i)-p_i,\quad \mathcal{F}(s):=\int_0^s f(\zeta)\, d\zeta, \quad g_i(p_i)=\int_0^{p_i}\frac {1} {\rho_i({\zeta})} \,d\zeta.$ \end{lemma} \begin{proof} First of all, let us prove that: for all $s_i\ge 0$ and $s^\star_i\ge 0$ such that $s_1+s_2=s^\star_1 +s^\star_2=1$, \label{magic} \begin{aligned} &\bigl(\rho_1(p_1)s_1-\rho_1(p_1^\star) s^\star_1\bigr)g_1(p_1)+ \bigl(\rho_2(p_2)s_2-\rho_2(p_2^\star) s^\star_2\bigr)g_2(p_2)\\ &\ge \mathcal{H}_1(p_1)s_1-\mathcal{H}_1(p_1^\star)s^\star_1 + \mathcal{H}_2(p_2)s_2-\mathcal{H}_2(p_2^\star)s^\star_2+ \mathcal{F}(s_1)-\mathcal{F}(s_1^\star). \end{aligned} Let us denote by $\mathcal{J}$ the left hand side of \eqref{magic}, $\mathcal{J} = \bigl(\rho_1(p_1)s_1-\rho_1(p_1^\star)s^\star_1\bigr)g_1(p_1)+ \bigl(\rho_2(p_2)s_2-\rho_2(p_2^\star) s^\star_2\bigr)g_2(p_2).$ Since the function $g_i$ is concave, we have $$\label{gconcave} g_i(p_i)\le g_i(p_i^\star) + g' _i(p_i^\star)(p_i-p_i^\star)=g_i(p_i^\star) + \frac{1}{\rho_i(p_i^\star)}(p_i-p_i^\star).$$ From the definition of $\mathcal{H}_i$, we have \begin{align*} \mathcal{J} &=\bigl[ \bigl(\rho_1(p_1)s_1g_1(p_1)-s_1 p_1\bigr)+ s_1 p_1 -\rho_1(p_1^\star)s^\star_1g_1(p_1)\bigr]\\ &\quad+\bigl[ \bigl(\rho_2(p_2)s_2g_2(p_2)-s_2 p_2\bigr) + s_2 p_2 -\rho_2(p_2^\star)s^\star_2g_2(p_2)\bigr]\\ &=s_1\mathcal{H}_1(p_1)+ s_1 p_1 -\rho_1(p_1^\star)s^\star_1g_1(p_1)+ s_2\mathcal{H}_2(p_2)+ s_2 p_2 -\rho_2(p_2^\star)s^\star_2g_2(p_2) \end{align*} and the concavity property of $g_i$ leads to \label{magic1} \begin{aligned} \mathcal{J} &\ge s_1 \mathcal{H}_1(p_1)- s_1^\star \mathcal{H}_1(p_1^\star) +s_2 \mathcal{H}_2(p_2)- s_2^\star \mathcal{H}_2(p_2^\star) +s_1p_1 -s_1^\star p_1 +s_2p_2 -s_2^\star p_2\\ &\ge s_1 \mathcal{H}_1(p_1)- s_1^\star \mathcal{H}_1(p_1^\star) +s_2 \mathcal{H}_2(p_2)- s_2^\star \mathcal{H}_2(p_2^\star) +s_1 \bigl(p_1 -p_2 \bigr) - s_1^\star \bigl(p_1 -p_2 \bigr) \\ &= s_1 \mathcal{H}_1(p_1)- s_1^\star \mathcal{H}_1(p_1^\star) +s_2 \mathcal{H}_2(p_2)- s_2^\star \mathcal{H}_2(p_2^\star) +\bigl( s_1 - s_1^\star \bigr) f(s_1). \end{aligned} Since the function $\mathcal{F}$ is convex, $$\label{convexeq} (s_1-s_1^\star)f(s_1)\ge \mathcal{F} (s_1)-\mathcal{F} (s_1^\star).$$ The above inequalities \eqref{magic1} and \eqref{convexeq} ensure that the assertion \eqref{magic} is satisfied. Let us multiply scalarly \eqref{gasdhn} with $g_1(p_{1,h}^{n+1})$ and add the scalar product of \eqref{waterdhn} with $g_2(p_{2,h}^{n+1})$, we have \label{long} \begin{aligned} &\frac 1 h \int_\Omega \phi\Big( \bigl( \rho_1(p_{1,h}^{n+1})s_{1,h}^{n+1}-\rho_1(p_{1,h}^{n}) s_{1,h}^{n} \bigr) g_1(p_{1,h}^{n+1}) \\ &+\bigl( \rho_2(p_{2,h}^{n+1})s_{2,h}^{n+1}-\rho_2(p_{2,h}^{n}) s_{2,h}^{n} \bigr) g_2(p_{2,h}^{n+1}) \Big)\, dx\\ &+\int_\Omega\mathbf{K}M_1 (s_{1,h}^{n+1}) \nabla p_{1,h}^{n+1} \cdot \nabla p_{1,h}^{n+1} \, dx +\int_\Omega\mathbf{K}M_2 (s_{2,h}^{n+1}) \nabla p_{2,h}^{n+1} \cdot \nabla p_{2,h}^{n+1} \, dx\\ &+\eta\int_\Omega |\nabla f(s_{1,h}^{n+1})|^2\,dx\\ &=\int_\Omega\mathbf{K}M_1(s_{1,h}^{n+1}) \rho_1(p_{1,h}^{n+1})\mathbf{g}\cdot \nabla p_{1,h}^{n+1}\,dx\\ &\quad +\int_\Omega\mathbf{K}M_2(s_{2,h}^{n+1}) \rho_2(p_{2,h}^{n+1})\mathbf{g}\cdot \nabla p_{1,h}^{n+1}\,dx -\int_\Omega \rho_1(p_{1,h}^{n+1})s_{1,h}^{n+1}(f_{P})_h^{n+1}g_1(p_{1,h}^{n+1}) \,dx\\ &\quad -\int_\Omega\rho_2(p_{2,h}^{n+1})s_{2,h}^{n+1}(f_{P})_h^{n+1}g_2(p_{2,h}^{n+1})\,dx+ \int_\Omega \rho_1(p_{1,h}^{n+1})(s_{1}^{I})_h^{n+1}(f_{I})_h^{n+1} g_1(p_{1,h}^{n+1})\, dx\\ &\quad +\int_\Omega \rho_2(p_{2,h}^{n+1})(s_{2}^{I})_h^{n+1}(f_{I})_h^{n+1} g_2(p_{2,h}^{n+1})\, dx. \end{aligned} Using \eqref{magic} and following the proof of Lemma \ref{lemunifeps}, one gets \eqref{nablapn}. \end{proof} For a given sequence $(u_h^n)_n$, let us denote \label{defraccord} \begin{aligned} u_h(0) = u_h^0,\\ u_h(t) =\sum_{n=0}^{N-1} u_h^{n+1} \chi_{]nh,(n+1)h]}(t),\quad \forall t\in ]0,T] \end{aligned} and $$\label{defraccordtilde} \tilde u_h(t) = \sum_{n=0}^{N-1} \bigl ((1+n-\frac t h )u_h^{n} + (\frac t h -n) u_h^{n+1}\bigr )\chi_{[nh,(n+1)h]}(t),\quad \forall t\in [0,T].$$ Then $$\partial_t \tilde u_h(t) =\frac 1 h \sum_{n=0}^{N-1} ((u_h^{n+1} - u_h^{n})\chi_{]nh,(n+1)h[}(t),\quad \forall t\in [0,T]\backslash \{\cup_{n=0}^{N} nh\}$$ Let the functions $p_{i,h}$ and $s_{i,h}$ be defined as in \eqref{defraccord}. For $i=1,2$, we denote by $r_{i,h}$ the function defined similarly to \eqref{defraccord} corresponding to $r_{i,h}^n=\rho_i(p_{i,h}^n)s_{i,h}^n$ and $\tilde r_{i,h}$ the function defined similarly to \eqref{defraccordtilde} corresponding to $r^n_{i,h}$. In the same way, we denote by $f_{P,h}$, $f_{I,h}$ and $(s_i ^I)_h$ the functions corresponding to $(f_{P})_h^{n+1}$, $(f_{I})_h^{n+1}$ and $(s_i ^I)_h^{n+1}$ respectively. \begin{proposition} \label{prop3.1} We have \begin{gather}\label{esthsh1} (s_{2,h})_h \quad \text{is uniformly bounded in } L^2(0,T;H_{\Gamma_1}^1(\Omega)),\\ \label{esthph1} (p_{i,h})_h \quad \text{is uniformly bounded in } L^2(0,T;H_{\Gamma_1}^1(\Omega)), i=1,2\\ \label{esthrh1} (r_{i,h})_h \quad\text{is uniformly bounded in } L^2(0,T;H^1(\Omega)), i=1,2\\ \label{estrtldh} (\tilde{r}_{i,h})_h \quad\text{ is uniformly bounded in } L^2(0,T;H^1(\Omega)), i=1,2\\ \label{esthdtr} (\phi \partial_t \tilde r_{i,h})_h \quad \text{is uniformly bounded in } L^2(0,T;(H_{\Gamma_1}^1(\Omega))'), i=1,2. \end{gather} \end{proposition} \begin{proof} At the beginning of this proof, we indicate some useful remarks which can be established by a classical calculations, \begin{gather} \int_{Q_T} M_i(s_{i,h})|\nabla p_{i,h}|^2\,dx\,dt =h\sum_{n=0}^{N-1}\int_\Omega M_i(s_{i,h}^{n+1})|\nabla p_{i,h}^{n+1}|^2\,dx ~~(i=1,2.),\\ \int_{Q_T}|\nabla f(s_{1,h})|^2\,dx\,dt =h\sum_{n=0}^{N-1}\int_\Omega |\nabla f(s_{1,h}^{n+1})|^2\,dx,\\ \int_{Q_T}|f_p(t,x)|^2\,dt\,dx \ge h\sum_{n=0}^{N-1}\|(f_p)_h^{n+1}\|_{L^2(\Omega)},\\ \int_{Q_T}|f_I(t,x)|^2\,dt\,dx \ge h\sum_{n=0}^{N-1}\|(f_I)_h^{n+1}\|_{L^2(\Omega)}. \end{gather} Now, multiply \eqref{nablapn} by $h$ and summing it from $n=0$ to $n=N-1$, \label{esthph1bis} \begin{aligned} &\int_\Omega \phi\mathcal{H}_1(p_{1,h}(T))s_{1,h}(T)+\phi\mathcal{H}_2(p_{2,h}(T))s_{2,h}(T)\, dx \\ &+k_0 \int_{Q_T}M_1(s_{1,h}) \vert{ \nabla} p_{1,h}\vert^2 \,dx\,dt +k_0 \int_{Q_T}M_2(s_{2,h})\vert{ \nabla} p_{2,h}\vert^2 \,dx\,dt \\ &+\eta\int_{Q_T}|\nabla f(s_{1,h})|^2\,dx\,dt\\ & \le \int_\Omega \bigl( \phi\mathcal{H}_1(p_{1,h}(0))s_{1,h}(0)+\phi\mathcal{H}_2(p_{2,h}(0))s_{2,h}(0)\,\bigr) dx\\ &\quad +\mathcal{F}(s_{1,h}(0))-\mathcal{F}(s_{1,h}(T))+ C \bigl(1+\Vert f_{P} \Vert^2_{L^2(Q_T)} +\Vert f_{I} \Vert^2_{L^2(Q_T)}\bigr), \end{aligned} where $C$ is a constant independent of $h$. The positivity of the first term on the left hand side of \eqref{esthph1bis} ensures that there exists a positive constant $C$ independent of $h$ such that \begin{align*} &k_0 \int_{Q_T}M_1(s_{1,h}) \vert{ \nabla} p_{1,h}\vert^2 \,dx\,dt +k_0 \int_{Q_T}M_2(s_{2,h})\vert{ \nabla} p_{2,h}\vert^2 \,dx\,dt\\ &+\eta\int_{Q_T}|\nabla f(s_{1,h})|^2\,dx\,dt\le C, \end{align*} since we have, \begin{align*} &\int_{Q_T}M_1(s_{1,h}) \vert{ \nabla} p_{1,h}\vert^2 \,dx\,dt + \int_{Q_T}M_2(s_{2,h})\vert{ \nabla} p_{2,h}\vert^2 \,dx\,dt\\ &=\int_{Q_T}M(s_{1,h})\vert{ \nabla} p_{h}\vert^2 \,dx\,dt +\int_{Q_T}\frac{M_1(s_{1,h})M_2(s_{2,h})}{M(s_{1,h})}\vert \nabla f(s_{1,h})\vert^2\,dx\,dt, \end{align*} we deduce $$\label{love1} \int_{Q_T}M(s_{1,h})\vert{ \nabla} p_{h}\vert^2 \,dx\,dt +\eta\int_{Q_T}|\nabla f(s_{1,h})|^2\,dx\,dt\le C.$$ For the first estimate \eqref{esthsh1} and first of all, let us indicate to the fact that, $$p_{1,h}(t,x)-p_{2,h}(t,x)=0=f(s_{1,h}(t,x))\quad\text{for } x\in\Gamma_{1}$$ which gives ${s_{2,h}}_{|\Gamma_{1}}=0$. The assumption (H6) on the capillary function $f$ with the second term on the right hand side of \eqref{love1} lead to $$\int_{Q_T}|\nabla s_{1,h}|^2\,dx\,dt\le C,$$ where $C$ is a constant independent of $h$, which establishes \eqref{esthsh1}. Since we have $\nabla p_{1,h}=\nabla p_{h}+\frac{M_2}{M} \nabla f(s_{1,h}) \quad\text{and}\quad \nabla p_{2,h}=\nabla p_{h}-\frac{M_1}{M} \nabla f(s_{1,h}),$ the estimate \eqref{esthph1} becomes a consequence of \eqref{love1}. The uniform estimate \eqref{esthrh1} is a consequence of the two previous ones since the densities $\rho_i$ are bounded and of class $\mathcal{C}^1$functions as well as the saturations $0\le s_{i,h}\le 1$, $$\nabla r_{i,h}= \sum_{n=0}^{N-1} \bigl (\rho_i'(p_{i,h}^{n+1})s_{i,h}^{n+1} \nabla p_{i,h}^{n+1} +\rho_i(p_{i,h}^{n+1}) \nabla s_{i,h}^{n+1}\bigr) \chi_{]nh,(n+1)h]}(t).$$ Now, for estimate \eqref{estrtldh} we have \begin{aligned} \nabla \tilde{r}_{i,h} &= \sum_{n=0}^{N-1} \bigl ((1+n-\frac{t}{h}) [\rho_i'(p_{i,h}^{n})s_{i,h}^{n}\nabla p_{i,h}^{n} +\rho_i(p_{i,h}^{n})\nabla s_{i,h}^{n}] \\ &\quad +(\frac{t}{h}-n) [\rho_i'(p_{i,h}^{n+1}) s_{i,h}^{n+1}\nabla p_{i,h}^{n+1} +\rho_i(p_{i,h}^{n+1}) \nabla s_{i,h}^{n+1}]\bigr) \chi_{]nh,(n+1)h]}(t). \end{aligned} since the densities $\rho_i$ are bounded and of class $\mathcal{C}^1$ functions as well as the saturations $0\le s_{i,h}^n\le 1$, $|\nabla \tilde{r}_{i,h}|^2\le C \sum_{n=0}^{N-1} \bigl( |\nabla p_{i,h}^{n}|^2+ |\nabla s_{i,h}^{n}|^2 + |\nabla p_{i,h}^{n+1}|^2+ |\nabla s_{i,h}^{n+1}|^2\bigr) \chi_{]nh,(n+1)h]}(t),$ and this implies $||\nabla \tilde{r}_{i,h}||^2_{L^2(Q_T)} \le C (\|\nabla p_{i,h}^{0}\|^2_{L^2(\Omega)} +\|\nabla s_{i,h}^{0}\|^2_{L^2(\Omega)} +\|\nabla p_{i,h}\|^2_{L^2(Q_T)} +\|\nabla s_{i,h}\|^2_{L^2(Q_T)}),$ where $C$ is a constant independent of $h$, and the estimate \eqref{estrtldh} is established.\\ From equations \eqref{gasdhn} and \eqref{waterdhn}, we have for all $\varphi\in L^2(0,T;H^1_{\Gamma_1}(\Omega))$, \begin{align*} &\langle \phi \partial_t \tilde r_{i,h}, \varphi \rangle\\ &=-\int_{Q_T}\mathbf{K} M_i(s_{i,h})\rho_i(p_{i,h}) \nabla p_{i,h}\cdot\nabla \varphi \,dx\,dt \\ &\quad+\int_{Q_T}\mathbf{K}\rho_i^2(p_{i,h}) M_i(s_{i,h})\mathbf{g}\cdot\nabla \varphi \,dx\,dt +\eta (-1)^i\int_{Q_T} \nabla (p_{1,h}-p_{2,h})\cdot\nabla \varphi \,dx\,dt \\ &\quad-\int_{Q_T}\rho_i(p_{i,h})s_{i,h} f_{P,h} \varphi\,dx\,dt +\int_{Q_T}\rho_i(p_h)s^I_{i,h} f_{I,h} \varphi\,dx\,dt. \end{align*} The above estimates \eqref{esthsh1}--\eqref{esthph1} with \eqref{love1} ensure that $(\phi \partial_t \tilde r_{i,h})_h$ is uniformly bounded in $L^2(0,T;(H_{\Gamma_1}^1(\Omega))').$ \end{proof} The next step is to pass from an elliptic problem to a parabolic one. Then, we pass to the limit on $h$, using some compactness theorems. \begin{proposition}[Convergence with respect to $h$] \label{propconvh} We have the following convergence as $h$ goes to zero, \begin{gather} \label{cvhrmr} \Vert r_{i,h} -\tilde r_{i,h}\Vert_{L^2(Q_T)}\to 0,\\ \label{cvhs} s_{2,h} \rightharpoonup s_2 \quad \text{ weakly in } L^2(0,T;H^1_{\Gamma_1}(\Omega)),\\ \label{cvhp} p_{i,h} \rightharpoonup p_i \quad \text{ weakly in } L^2(0,T;H^1_{\Gamma_1}(\Omega)),\\ \label{cvhr} r_{i,h} \to r_i \quad \text{ strongly in } L^2(Q_T). \end{gather} Furthermore, \begin{gather} \label{cvpps} s_{i,h}\to s_i \quad\text{almost everywhere in } Q_T,\\ 0\le s_i\le 1 \quad \text{almost everywhere in } Q_T, \label{spmh}\\ \label{cvppp} p_{i,h}\to p_i \quad \text{almost everywhere in } Q_T, \\ \label{r=rhos} r_i=\rho_i(p_i)s_i \text{ almost everywhere in } Q_T. \end{gather} Finally, we have $$\label{cvhdtr} \phi\partial_t\tilde r_{i,h} \rightharpoonup \phi\partial_t (\rho_i(p_i)s_i) \quad \text{weakly in } L^2(0,T;(H^1_{\Gamma_1}(\Omega))').$$ \end{proposition} \begin{proof} Note that \begin{align*} \Vert r_{i,h} -\tilde r_{i,h}\Vert^2_{L^2(Q_T)} &=\sum_{n=0}^{N-1} \int_{nh}^{(n+1)h} \Vert((1+n-\frac t h )(r_{i,h}^{n+1}-r_{i,h}^{n})\Vert^2_{L^2(\Omega)} \,dt\\ &=\frac h 3 \sum_{n=0}^{N-1}\Vert r_{i,h}^{n+1}-r_{i,h}^{n}\Vert^2_{L^2(\Omega)}. \end{align*} We multiply scalarly \eqref{gasdhn} and \eqref{waterdhn} respectively with $r_{1,h}^{n+1}-r_{1,h}^{n}$ and $r_{2,h}^{n+1}-r_{2,h}^{n}$. Then, summing for $n=0$ to $N-1$, we obtain, for $i=2$, \begin{align*} &\frac {\phi_0} h \sum_{n=0}^{N-1}\Vert r_{2,h}^{n+1}-r_{2,h}^{n} \Vert^2_{L^2(\Omega)}\\ &\le C\sum_{n=0}^{N-1}\bigl (\Vert\nabla r_{2,h}^{n} \Vert^2_{L^2(\Omega)}+\Vert\nabla r_{2,h}^{n+1}\Vert^2_{L^2(\Omega)} +\Vert\nabla s_{2,h}^{n+1}\Vert^2_{L^2(\Omega)}+\Vert\nabla p_{2,h}^{n+1} \Vert^2_{L^2(\Omega)}\\ &\quad +\Vert (f_{P})^{n+1}_h\Vert^2_{L^2(\Omega)}+\Vert (f_{I})^{n+1}_h \Vert^2_{L^2(\Omega)}\bigr ). \end{align*} This yields \begin{align*} &\sum_{n=0}^{N-1}\Vert r_{2,h}^{n+1}-r_{2,h}^{n}\Vert^2_{L^2(\Omega)}\\ &\le C\bigl (1+\Vert\nabla r_{2,h}\Vert^2 +\Vert\nabla s_{2,h}\Vert^2_{L^2(Q_T)}+\Vert\nabla p_{2,h} \Vert^2_{L^2(Q_T)}\\ &\quad +\Vert f_{P}\Vert^2_{L^2(Q_T)}+\Vert f_{I}\Vert^2_{L^2(Q_T)}\bigr ). \end{align*} From \eqref{esthsh1},\eqref{esthph1}, and \eqref{esthrh1}, we conclude that $\Vert r_{2,h} -\tilde r_{2,h}\Vert_{L^2(Q_T)}\to 0$. For $i=1$, \begin{align*} &\frac {\phi_0} h \sum_{n=0}^{N-1}\Vert r_{1,h}^{n+1}-r_{1,h}^{n} \Vert^2_{L^2(\Omega)}\\ &\le C\sum_{n=0}^{N-1}\bigl (\Vert\nabla r_{1,h}^{n}\Vert^2_{L^2(\Omega)} +\Vert\nabla r_{1,h}^{n+1}\Vert^2_{L^2(\Omega)}\\ &\quad +\Big| \int_{\Gamma_1} {\bf{K}} \rho_1^2(p^{n+1}_{1,h}) M_1(s_{1,h}^{n+1}) {\bf{g}}.\nu (r_{1,h}^{n+1}-r_{1,h}^{n})\, d\gamma \Big| \\ &\quad +\Vert\nabla s_{2,h}^{n+1}\Vert^2_{L^2(\Omega)} +\Vert\nabla p_{1,h}^{n+1}\Vert^2_{L^2(\Omega)} +\Vert (f_{P})^{n+1}_h\Vert^2_{L^2(\Omega)} +\Vert (f_{I})^{n+1}_h\Vert^2_{L^2(\Omega)}\bigr ). \end{align*} where $\nu$ is the outward normal to the injection boundary. This yields, with the help of trace theory to handle with the third term in the right hand side namely ($\Vert r_{1,h}^{n+1}-r_{1,h}^{n}\Vert_{L^2(\Gamma_1)} \le C \Vert \nabla (r_{1,h}^{n+1}-r_{1,h}^{n})\Vert_{L^2(\Omega)}$), that \begin{align*} \sum_{n=0}^{N-1}\Vert r_{1,h}^{n+1}-r_{1,h}^{n}\Vert^2_{L^2(\Omega)} & \le C\bigl (1+\Vert\nabla r_{1,h}\Vert^2+\Vert\nabla s_{2,h} \Vert^2_{L^2(Q_T)}+\Vert\nabla p_{1,h}\Vert^2_{L^2(Q_T)}\\ &\quad +\Vert f_{P}\Vert^2_{L^2(Q_T)} +\Vert f_{I} \Vert^2_{L^2(Q_T)}\bigr). \end{align*} From \eqref{esthsh1},\eqref{esthph1}, and \eqref{esthrh1}, we conclude that $%\label{difrrtld} \Vert r_{1,h} -\tilde r_{1,h}\Vert_{L^2(Q_T)}\to 0,$ and this achieves \eqref{cvhrmr}. From \eqref{esthph1} \eqref{esthsh1}, the sequences $(p_{i,h})_h$, $(s_{2,h})_h$ are uniformly bounded in $L^2(0,T;H_{\Gamma_1}^1(\Omega))$, we have up to a subsequence the convergence results \eqref{cvhs}, \eqref{cvhp}.\\ The sequences $(\tilde r_{i,h})_h$ are uniformly bounded in $L^2(0,T;H^1(\Omega))$. In light of \eqref{esthdtr} we have the strong convergence $$\label{cvrtld} \tilde r_{i,h} \to r_i \text{ strongly in } L^2(Q_T).$$ This compactness result is classical and can be found in \cite{A-SI87}, \cite{B-CJ86} when the porosity is constant, and under the assumption (H1) (the porosity belongs to $W^{1,\infty}(\Omega)$), the proof can be adapted with minor modifications. The convergence \eqref{cvrtld} with \eqref{cvhrmr} ensures the following strong convergence \begin{gather} \label{cvrho1p1h1} \rho_1(p_{1,h})s_{1,h} \to r_1 \quad \text{strongly in $L^2(Q_T)$ and a.e. in } Q_T,\\ \label{cvrho2p2h2} \rho_2(p_{2,h})s_{2,h} \to r_2 \quad \text{strongly in $L^2(Q_T)$ and a.e. in } Q_T, \end{gather} and this achieves \eqref{cvhr}. We are now concerned with almost everywhere convergence on pressures $p_{i,h}$ and saturations $s_{i,h}$. Denote $$u= \rho_1(p_{1,h})s_{1,h},\quad v= \rho_2(p_{1,h}-f(s_{1,h}))(1-s_{1,h}).$$ Define the map $H : \mathbb{R}^+ \times \mathbb{R}^+ \to \mathbb{R}^+ \times [0,1]$ defined by $$H(u,v) = (p_{1,h},s_{1,h}) \label{defHb}$$ where $u$ and $v$ are solutions of the system \begin{gather*} u(p_{1,h},s_{1,h}) = \rho_1(p_{1,h})s_{1,h},\\ v(p_{1,h},s_{1,h}) = \rho_2(p_{1,h}-f(s_{1,h}))(1-s_{1,h}). \end{gather*} Note that $H$ is well defined as a diffeomorphism, \begin{align*} &\left| \begin{matrix} \frac{\partial u }{\partial p_{1,h}} & \frac{\partial u }{\partial s_{1,h}} \\ \frac{\partial v }{\partial p_{1,h}} & \frac{\partial v }{\partial s_{1,h}} \end{matrix} \right|\\ &= -\rho_1'(p_{1,h})\rho_2(p_{1,h}-f(s_{1,h}))s_{1,h} -\rho_1(p_{1,h})\rho_2'(p_{1,h}-f(s_{1,h}))(1-s_{1,h}) \\ &\quad -\rho_1'(p_{1,h})s_{1,h}(1- s_{1,h}) \rho_2'(p_{1,h}-f(s_{1,h}))f'(s_{1,h})<0. \end{align*} As we have the almost everywhere convergence \eqref{cvrho1p1h1}, \eqref{cvrho2p2h2} and the map $H$ defined in \eqref{defH} is continuous, we deduce \begin{gather*} p_{1,h} \to p_1 \quad \text{a.e. in } Q_T.\\ s_{1,h} \to s_1 \quad \text{a.e. in } Q_T. \end{gather*} The identification of the limit is due to \eqref{esthph1}, \eqref{esthsh1}. The continuity of the capillary pressure function ensures that $p_{2,h} \to p_2 \text{ a.e. in } Q_T,$ the saturation equation ensures also, $s_{2,h} \to s_2 \text{ a.e. in } Q_T,$ and this achieve \eqref{cvpps}, \eqref{cvppp}. The maximum principle \eqref{spmh} and the identification \eqref{r=rhos} are conserved through a limit process. Finally the weak convergence \eqref{cvhdtr} is a consequence of \eqref{esthdtr}, and the identification of the limit is due to \eqref{r=rhos}. \end{proof} The technique for obtaining solutions of the system \eqref{conslawnondeg1}--\eqref{conslawnondeg2} is to pass to the limit as h goes to zero on the solutions of \label{non-degttld} \begin{aligned} &\phi\partial_t( \tilde{r}_{i,h}) - \operatorname{div} (\mathbf{K}M_i(s_{i,h})\rho_i(p_{i,h}) \nabla p_{i,h}) +\operatorname{div} (\mathbf{K}M_i(s_{i,h})\rho^2_i(p_{i,h}) \mathbf{g})\\ &+ (-1)^i \eta\operatorname{div} (\rho_i(p_{i,h})\nabla(p_{1,h}-p_{2,h})) +\rho_i(p_{i,h})s_{i,h} f_{P,h}\\ &= \rho_i(p_{i,h})s_i^If_{I,h} \end{aligned} We remark that this system ($i=1,2$) is nothing else than \eqref{gasdhn}-\eqref{waterdhn}, written for $n=0$ to $N-1$ by using the definition \eqref{defraccord} and \eqref{defraccordtilde}. Let us consider the weak formulations ($i=1,2$) on which we have to pass to the limit \label{weakfh} \begin{aligned} &\langle \phi \partial_t \tilde r_{i,h}, \varphi_i \rangle +\int_{Q_T}\mathbf{K} M_i(s_{i,h})\rho_i(p_{i,h}) \nabla p_{i,h} \cdot\nabla \varphi_i \,dx\,dt \\ &-\int_{Q_T}\mathbf{K}\rho_i^2(p_{i,h}) M_i(s_{i,h})\mathbf{g}\cdot\nabla \varphi_i \,dx\,dt\\ &-(-1)^i \eta\int_{Q_T} \rho_i(p_{i,h})\nabla(p_{1,h}-p_{2,h}) \cdot\nabla \varphi_i \,dx\,dt\\ &+\int_{Q_T}\rho_i(p_{i,h})s_{i,h} f_{P,h} \varphi_i\,dx\,dt\\ &=\int_{Q_T}\rho_i(p_h)s^I_{i,h} f_{I,h} \varphi_i\,dx\,dt. \end{aligned} where $\varphi_i$ ($i=1,2$) belongs to $L^2(0,T;H^1_{\Gamma_1}(\Omega))$. Next, we pass to the limit on each term of \eqref{weakfh} which is conserved by the previous proposition. The passage to the limit on the first term is due to \eqref{cvhdtr}, for the second term we have $M_i(s_{i,h})\rho_i(p_{i,h}) \nabla \varphi_i$ converges almost everywhere in $Q_T$ and dominated which leads by Lebesgue theorem to a strong convergence in $L^2(Q_T)$ and by virtue of the weak convergence \eqref{cvhp} we establish the convergence of the second term of \eqref{weakfh} to the desired term. The last three terms converge obviously to the wanted limit due to the previous proposition and Lebesgue theorem. We then have established the weak formulation \eqref{non-degv1} of theorem \ref{non-degt}. Furthermore, we have well obtained by proposition \ref{propconvh}, \begin{gather*} 0\le s_i(t,x)\le 1 \quad \text{a.e. in } Q_T, \quad s_2\in L^2(0,T;H^1_{\Gamma_1}(\Omega)),\\ p_i\in L^2(0,T;H^1_{\Gamma_1}(\Omega)), \quad \phi\partial_t (\rho_i(p_i)s_i) \in L^2(0,T;(H^1_{\Gamma_1}(\Omega))'), \quad i=1,2. \end{gather*} We recall that $s_i$ is a given function of $p_1$ and $p_2$. The compactness property on $\rho_i(p_{i,h})s_{i,h}$ implies $\rho_i(p_i)s_{i}\in C^0([0,T];L^2(\Omega))$, for $i=1,2$. Theorem \ref{non-degt} is then proved. \section{Proof of Theorem \ref{main} (Degenerate case)} \label{sec-maint} The proof is based on the existence result established for the non-degenerate case and the compactness lemma \ref{lemcompact}. \begin{lemma}\label{lemunifeta} The sequences $({s_i^\eta})_\eta$, $(p^{\eta}:=p_2^{\eta}+\tilde{p}({s_1^\eta}))_\eta$ defined by Theorem \ref{non-degt} satisfy \begin{gather} 0\le {s_i^\eta} (t,x) \le 1 \quad \text{a.e. in $x\in \Omega$ for all } t\in[0,T], \label{maxprin1}\\ (p^{\eta})_\eta \quad\text{is uniformly bounded in } L^2(0,T;H^1_{\Gamma _1}({\Omega})), \label{pp}\\ (\sqrt{\eta} \nabla f({s_1^\eta}))_\eta \text{is uniformly bounded in}~ L^2(Q_T)\label{epsestimate}\\ (\sqrt{M_i({s_i^\eta})}\nabla {p_i^\eta})_\eta \quad \text{is uniformly bounded in } L^2(Q_T), \label{mpeta}\\ (\beta({s_1^\eta}))_\eta \quad \text{is uniformly bounded in } L^2(0,T;H^1({\Omega})), \label{betabeta}\\ (\phi\partial_t(\rho_i(p_i^\eta){s_i^\eta}))_\eta \quad \text{is uniformly bounded in } L^2(0,T;H^1_{\Gamma _1}({\Omega})').\label{drhoi} \end{gather} \end{lemma} \begin{proof} The maximum principle \eqref{maxprin1} is conserved through the limit process. For the next three estimates, consider the $L^2(\Omega)$ scalar product of $\eqref{conslawnondeg1}$ by $g_1({p_1^\eta})=\int_0^{{p_1^\eta}}\frac{1}{\rho_1(\xi)} \, d\xi$ and \eqref{conslawnondeg2} by $g_2({p_2^\eta})=\int_0^{{p_2^\eta}}\frac{1}{\rho_2(\xi)} \, d\xi$ and adding them after denoting by $\mathcal{H}_i({p_i^\eta})=\rho_i({p_i^\eta}) g_i({p_i^\eta})-{p_i^\eta}$ $(i=1,2)$, then we have \label{estequ} \begin{aligned} &\frac{d}{dt}\int_\Omega \phi \Big( {s_1^\eta} \mathcal{H}_1({p_1^\eta}) +{s_2^\eta} \mathcal{H}_2({p_2^\eta}) + \int_{0}^{{s_1^\eta}}f(\xi)\, d\xi \Big)dx +\int_\Omega \mathbf{K} M_1({s_1^\eta}) \nabla {p_1^\eta} \cdot \nabla {p_1^\eta} \, dx \\ &+\eta\int_\Omega| \nabla f({s_1^\eta})|^2\,dx +\int_\Omega \mathbf{K} M_2({s_2^\eta}) \nabla {p_2^\eta} \cdot \nabla {p_2^\eta} \,dx \\ &=\int_\Omega \mathbf{K} M_1({s_1^\eta}) \rho_1({p_1^\eta})\mathbf{g} \cdot \nabla {p_1^\eta} \, dx +\int_\Omega \mathbf{K} M_2({s_2^\eta}) \rho_2({p_2^\eta})\mathbf{g} \cdot \nabla {p_2^\eta} \, dx \\ &\quad + \int_\Omega \rho_1({p_1^\eta}) s_1^If_I g_1({p_1^\eta}) \, dx- \int_\Omega \rho_1({p_1^\eta}) {s_1^\eta} f_p g_1({p_1^\eta}) \, dx \\ &\quad -\int_\Omega \rho_2({p_2^\eta}) {s_2^\eta} f_p g_2({p_2^\eta}) \, dx +\int_\Omega \rho_2({p_2^\eta}) s_2^If_I g_2({p_2^\eta}) \, dx. \end{aligned} Integrate \eqref{estequ} over$(0,T)$ to obtain \label{estinequT} \begin{aligned} &\int_\Omega\phi \big( {s_1^\eta} \mathcal{H}_1({p_1^\eta}) +{s_2^\eta} \mathcal{H}_2({p_2^\eta}) \big)(x,T)\,dx+\int_{Q_T} \mathbf{K} M_1({s_1^\eta}) \nabla {p_1^\eta} \cdot \nabla {p_1^\eta} \, dx\,dt \\ &+\eta\int_{Q_T}| \nabla f({s_1^\eta})|^2\,dx\,dt +\int_{Q_T} \mathbf{K} M_2({s_2^\eta}) \nabla {p_2^\eta} \cdot \nabla {p_2^\eta} \,dx\,dt \\ &=\int_\Omega\phi \big( s_1^0 \mathcal{H}_1(p_1^0) +s_2^0 \mathcal{H}_2(p_2^0) \big)\,dx -\int_\Omega\int_{s_1^0}^{{s_1^\eta}(x,T)}f(\xi)\, d\xi dx \\ &\quad +\int_{Q_T} \mathbf{K} M_1({s_1^\eta}) \rho_1({p_1^\eta})\mathbf{g} \cdot \nabla {p_1^\eta} \, dx\,dt +\int_{Q_T} \mathbf{K} M_2({s_2^\eta}) \rho_2({p_2^\eta})\mathbf{g} \cdot \nabla {p_2^\eta} \, dx\,dt \\ &\quad + \int_{Q_T} \rho_1({p_1^\eta}) s_1^If_I g_1({p_1^\eta}) \, dx\,dt - \int_{Q_T} \rho_1({p_1^\eta}) {s_1^\eta} f_p g_1({p_1^\eta}) \, dx\,dt \\ &\quad -\int_{Q_T} \rho_2({p_2^\eta}) {s_2^\eta} f_p g_2({p_2^\eta}) \, dx\,dt +\int_{Q_T} \rho_2({p_2^\eta}) s_2^If_I g_2({p_2^\eta}) \, dx\,dt. \end{aligned} The first term on the left hand side of \eqref{estinequT} is positive and the two first terms on the right hand side are bounded since $p_i^0\in L^2(\Omega)$ and $0\le s_i^0\le 1$. The third and the fourth terms on the right hand side, corresponding to gravity term, can be absorbed by the degenerate dissipative term on pressures (namely the second and fourth terms on the left hand side of \eqref{estinequT}) since: $$\big|\int_{Q_T} \mathbf{K} M_i(s_i^\eta) \rho_i(p_i^\eta)\mathbf{g} \cdot \nabla p_i^\eta\, dx\,dt\big| \le C +\frac{k_0}{2} \int_{Q_T} M_i(s_i^\eta) | \nabla p_i^\eta |^2\; dx\,dt, \quad i=1,2.$$ Finally, using the fact that the functions $g_i$ $(i=1,2)$ are sublinear, we deduce from \eqref{estinequT} that \label{magicineqeps} \begin{aligned} & \int_{Q_T} M_1({s_1^\eta}) |\nabla {p_1^\eta}|^2 \, dx\,dt +\int_{Q_T} M_2({s_2^\eta}) |\nabla {p_2^\eta} |^2 \,dx\,dt\\ &+\eta\int_{Q_T} \nabla f({s_1^\eta}) \cdot \nabla f({s_1^\eta}) \, dx\,dt \\ &\le C(1+ \| {p_1^\eta}\|_{L^2(Q_T)}+\| {p_2^\eta}\|_{L^2(Q_T)}), \end{aligned} where $C$ is a constant independent of $\eta$. From the definition of the global pressure, we have \begin{align}\label{defpeta} \nabla p^\eta=\nabla {p_2^\eta} + \frac{M_1({s_1^\eta})}{M({s_1^\eta})} \nabla f({s_1^\eta}) =\nabla {p_1^\eta} - \frac{M_2({s_2^\eta})}{M({s_1^\eta})} \nabla f({s_1^\eta}), \end{align} and consequently, \label{magic1eta} \begin{aligned} &\int_{Q_T} M({s_1^\eta}) |\nabla p^\eta |^2\,dx\,dt + \int_{Q_T} \frac{M_1({s_1^\eta}) M_2({s_2^\eta})}{M({s_1^\eta})} |\nabla f({s_1^\eta})|^2 \, dx\,dt\\ &=\int_{Q_T} M_1({s_1^\eta}) \nabla {p_1^\eta} \cdot \nabla {p_1^\eta} \, dx\,dt +\int_{Q_T} M_2({s_2^\eta}) \nabla {p_2^\eta} \cdot \nabla {p_2^\eta} \,dx\,dt. \end{aligned} On the other hand, $\| {p_1^\eta}\|_{L^2(Q_T)}\le \| p^\eta\|_{L^2(Q_T)} +\| {\bar p}(s_1^\eta)\|_{L^2(Q_T)} \le C\|\nabla p^\eta\|_{L^2(Q_T)} +\| {\bar p}(s_1^\eta)\|_{L^2(Q_T)},$ due to Poincar\'e's inequality, in the same way we have $$\| {p_2^\eta}\|_{L^2(Q_T)}\le C\|\nabla p^\eta\|_{L^2(Q_T)} +\| {\tilde p}(s_1^\eta)\|_{L^2(Q_T)}.$$ From the above estimates and \eqref{magic1eta}, the estimate \eqref{magicineqeps} yields \label{magic2eta} \begin{aligned} &\int_{Q_T} M({s_1^\eta}) |\nabla p^\eta |^2\,dx\,dt + \int_{Q_T} \frac{M_1({s_1^\eta}) M_2({s_2^\eta})}{M({s_1^\eta})} |\nabla f({s_1^\eta})|^2 \, dx\,dt\\ &+\eta\int_{Q_T} \nabla f({s_1^\eta}) \cdot \nabla f({s_1^\eta}) \, dx\,dt\\ &\le C(1+\|\nabla p^\eta\|_{L^2(Q_T)}). \end{aligned} The Young inequality permits to absorb the last term by the first term on the left hand side of \eqref{magic2eta} to obtain \label{magicineqepsb} \begin{aligned} &\int_{Q_T} M({s_1^\eta}) |\nabla p^\eta|^2\,dx\,dt + \int_{Q_T} \frac{M_1({s_1^\eta}) M_2({s_2^\eta})}{M({s_1^\eta})} |\nabla f({s_1^\eta})|^2 \, dx\,dt \\ &+ \int_{Q_T} M_1({s_1^\eta}) \nabla {p_1^\eta} \cdot \nabla {p_1^\eta} \, dx\,dt +\int_{Q_T} M_2({s_2^\eta}) \nabla {p_2^\eta} \cdot \nabla {p_2^\eta} \,dx\,dt \\ &+\eta\int_{Q_T} \nabla f({s_1^\eta}) \cdot \nabla f({s_1^\eta}) \, dx\,dt \le C, \end{aligned} where $C$ is a constant independent of $\eta$. Thus, the estimates \eqref{epsestimate}--\eqref{mpeta} are established. The estimate \eqref{betabeta} is also a consequence of \eqref{magicineqeps} since \begin{align*} \int_{Q_T} |\nabla \beta ({s_1^\eta}) |^2\, dx\,dt &=\int_{Q_T} \frac{M_1^2({s_1^\eta}) M_2^2({s_2^\eta})}{M^2({s_1^\eta})} |\nabla f({s_1^\eta})|^2 \, dx\,dt\\ &\le \int_{Q_T} M_1({s_1^\eta}) M_2({s_2^\eta})|\nabla f({s_1^\eta})|^2 \, dx\,dt \le C. \end{align*} For the last estimate \eqref{drhoi}, let $\varphi_i \in L^2(0,T;H^1_{\Gamma_1}(\Omega))$ and denote the bracket $\langle\cdot,\cdot\rangle$ the duality product between $L^2(0,T;(H^1_{\Gamma_1}(\Omega))')$ and $L^2(0,T;H^1_{\Gamma_1}(\Omega))$, using \eqref{defp}, one gets \begin{align*}%\label{weakfeps} &|\langle \phi\partial_t(\rho_i(p_i^\eta){s_i^\eta}), \varphi_i \rangle|\\ &\le \big| \eta \int_{Q_T} \rho_i({p_i^\eta}) \nabla f({s_i^\eta})\cdot\nabla \varphi_i \,dx\,dt \big|\\ &\quad + \big|\int_{Q_T}\mathbf{K} \rho_i({p_i^\eta})(M_i({s_i^\eta}) \nabla p^\eta + \nabla \beta ({s_1^\eta}))\cdot\nabla \varphi_i \,dx\,dt\big| \\ &\quad +\big|\int_{Q_T}\mathbf{K}\rho_i^2({p_i^\eta}) M_i({s_i^\eta})\mathbf{g}\cdot\nabla \varphi_i \,dx\,dt\big| +\big|\int_{Q_T}\rho_i({p_i^\eta}){s_i^\eta} f_{P} \varphi_i\,dx\,dt\big| \\ &\quad + \big|\int_{Q_T}\rho_i({p_i^\eta})s^I_{i} f_{I} \varphi_i\,dx\,dt \big|, \end{align*} and from the estimates \eqref{pp}--\eqref{betabeta}, we deduce $$|\langle \phi\partial_t(\rho_i(p_i^\eta){s_i^\eta}), \varphi_i \rangle| \le C \|\varphi_i\|_{L^2(0,T;H^1_{\Gamma _1} ({\Omega}))},$$ which establishes \eqref{drhoi} and completes the proof of the lemma. \end{proof} \begin{lemma}[Compactness result for degenerate case] \label{lemcompact} For every $M$, the following set \begin{align*} E_M = \Big\{ &(\rho_1(p_1)s_1,\rho_2(p_2)s_2) \in L^2(Q_T)\times L^2(Q_T): \|\beta(s_1)\|_{L^2(0,T;H^1(\Omega))}\le M,\,\\ &\|\sqrt{M_1(s_1)}\nabla p_1\|_{L^2(Q_T)} +\|\sqrt{M_2(s_2)}\nabla p_2\|_{L^2(Q_T)}\le M,\,\\ &\|\phi\partial_t (\rho_1(p_1)s_1)\|_{L^2(0,T;(H_{\Gamma_1}^1(\Omega))')} + \|\phi\partial_t (\rho_2(p_2)s_2)\|_{L^2(0,T;(H_{\Gamma_1}^1(\Omega))')} \le M \Big\} \end{align*} is relatively compact in $L^2(Q_T)\times L^2(Q_T)$, and $\gamma(E_M)$ is relatively compact in $L^2(\Sigma_T)\times L^2(\Sigma_T)$, ($\gamma$ denotes the trace on $\Sigma_T$ operator). \end{lemma} \begin{proof} The proof is inspired by the compactness lemma 4.3 \cite[p. 37]{A-GS08}, which introduced for compressible degenerate model. We generalize this result for our compressible degenerate model. Denote by $$u= \rho_1(p_1)s_1,\quad v= \rho_2(p_2)(1-s_1).$$ Define the map $\mathbb{H} : \mathbb{R}^+ \times \mathbb{R}^+ \mapsto \mathbb{R}^+ \times [0,\beta(1)]$ defined by $$\mathbb{H}(u,v) = (p,\beta(s_1)) \label{defH}$$ where $u$ and $v$ are solutions of the system \begin{gather*} u(p,\beta(s_1)) = \rho_1(p-\bar{p}(\beta ^{-1}(\beta(s_1)))) \beta ^{-1}(\beta(s_1)), \\ v(p,\beta(s_1)) = \rho_2(p-\tilde{p}(\beta ^{-1}(\beta(s_1)))) (1-\beta ^{-1}(\beta(s_1)). \end{gather*} Note that $\mathbb{H}$ is well defined as a diffeomorphism, since \begin{align*} \frac{\partial u }{\partial p} &=\rho_1'(p-\bar{p}(\beta ^{-1}(\beta(s_1))))\beta ^{-1}(\beta(s_1))\ge0\\ \frac{\partial u }{\partial \beta} &= \rho_1'(p-\bar{p}(\beta ^{-1}(\beta(s_1)))) [ -\bar{p}'(\beta ^{-1}(\beta(s_1))) ({\beta ^{-1}}' (\beta(s_1)))]\beta ^{-1}(\beta(s_1)) \\ &\quad + \rho_1(p-\bar{p}(\beta ^{-1}(\beta(s_1)))) {\beta ^{-1}}' (\beta(s_1))\ge0\\ \frac{\partial v }{\partial p} &=-\rho_2'(p-\tilde{p}(\beta ^{-1}(\beta(s_1))))(1-\beta ^{-1} (\beta(s_1)))\ge0\\ \frac{\partial v }{\partial \beta} &= \rho_2'(p-\tilde{p}(\beta ^{-1}(\beta(s_1)))) [-\tilde{p}'(\beta ^{-1}(\beta(s_1))) ({\beta ^{-1}}' (\beta(s_1)))] [ 1-\beta ^{-1}(\beta(s_1))]\\ &\quad - \rho_2(p-\tilde{p}(\beta ^{-1}(\beta(s_1)))) {\beta ^{-1}}' (\beta(s_1))\le0,\\ \end{align*} and if one of the saturations is zero the other one is one, this conserves that the jacobian determinant of the map $\mathbb{H}^{-1}$ is strictly negative. Furthermore, $\mathbb{H}^{-1}$ is an H\"older function, in the sense that $u$ and $v$ are H\"older functions of order $\theta$ with $0<\theta \le 1$ . For that, let $(q_1,\sigma_1)$ and $(q_2,\sigma_2)$ in $\mathbb{R}^+ \times [0,\beta(1)]$, we have \begin{align*} &|u(q_1,\sigma_1)-u(q_2,\sigma_2)| \\ &= |\rho_1(q_1-\bar{p}(\beta ^{-1}(\sigma_1)))\beta ^{-1}(\sigma_1) - \rho_1(q_2-\bar{p}(\beta ^{-1}(\sigma_2)))\beta ^{-1}(\sigma_2)| \\ &\le |\rho_1(q_1-\bar{p}(\beta ^{-1}(\sigma_1)))-\rho_1(q_2-\bar{p}(\beta ^{-1}(\sigma_2)))| + \rho_M|\beta ^{-1}(\sigma_1)-\beta ^{-1}(\sigma_2)|, \end{align*} since $\beta ^{-1}$ is an H\"older function of order $\theta$, $0<\theta \le 1$, and the map $\rho_1$ is bounded and of class $C^1$, we deduce up to two cases: The first case $|q_1-q_2|\ge1$: \begin{align*} &|u(q_1,\sigma_1)-u(q_2,\sigma_2)|\\ &\le|\rho_1(q_1-\bar{p}(\beta ^{-1}(\sigma_1)))-\rho_1(q_2-\bar{p}(\beta ^{-1}(\sigma_2)))| + \rho_M|\beta ^{-1}(\sigma_1)-\beta ^{-1}(\sigma_2)|\\ &\le \rho_M + \rho_M|\beta ^{-1}(\sigma_1)-\beta ^{-1}(\sigma_2)|\\ &\le \rho_M |q_1-q_2|^\theta +\rho_M C_\beta|\sigma_1-\sigma_2|^\theta, \end{align*} for the other case $|q_1-q_2|<1$, we have \begin{align*} &|u(q_1,\sigma_1)-u(q_2,\sigma_2)| \\ &\le |\rho_1(q_1-\bar{p}(\beta ^{-1}(\sigma_1)))-\rho_1(q_2-\bar{p}(\beta ^{-1}(\sigma_2)))| + \rho_M|\beta ^{-1}(\sigma_1)-\beta ^{-1}(\sigma_2)|\\ &\le C (|q_1-q_2|+|\bar{p}(\beta ^{-1}(\sigma_1))-\bar{p}(\beta ^{-1}(\sigma_2))|) +\rho_M C_\beta|\sigma_1-\sigma_2|^\theta \\ &\le C |q_1-q_2|^\theta + C |\bar{p}(\beta ^{-1}(\sigma_1))-\bar{p}(\beta ^{-1}(\sigma_2))|+\rho_M C_\beta|\sigma_1-\sigma_2|^\theta \end{align*} further more one can easily show that $\bar{p}$ is a $C^1([0,1];\mathbb{R})$, it follows that \begin{align} |u(q_1,\sigma_1)-u(q_2,\sigma_2)| &\le C|q_1-q_2|^\theta + C|\sigma_1-\sigma_2|^\theta. \end{align} In the same way, we have \begin{align} |v(q_1,\sigma_1)-v(q_2,\sigma_2)| \le c_1 |q_1 - q_2|^\theta + c_2 |\sigma_1 - \sigma_2|^\theta. \end{align} For $0<\tau <1$, and $1 \theta/2$ (see \cite{lionsmagenes} for more details). Choosing for example $\tau'=\frac{3\theta}{4}$, we deduce the relative compactness of $\gamma(E_M)$ into $L^2(\Sigma_T)\times L^2(\Sigma_T)$. This completes the proof. \end{proof} From the previous two lemmas, we deduce the following convergence. \begin{lemma}[Strong and weak convergence] \label{lem4.3} Up to a subsequence, the sequences $({s_i^\eta})_\eta$, $(p^\eta)_\eta$, $({p_i^\eta})_\eta$ verify the following convergence \begin{gather} p^\eta \rightharpoonup p \quad \text{weakly in } L^2(0,T;H^1_{\Gamma_1}(\Omega)),\label{ph11}\\ \beta({s_1^\eta}) \rightharpoonup \beta(s_1) \quad \text{weakly in } L^2(0,T;H^1(\Omega)),\label{betah11}\\ p^\eta \to p \quad\text{almost everywhere in } Q_T, \label{pppp}\\ {s_1^\eta} \to s_1 \quad \text{almost everywhere in $Q_T$ and }\Sigma_T, \label{beta1123}\\ {s_1^\eta} \to s_1 \quad \text{strongly in $L^2(Q_T)$ and } L^2(\Sigma_T), \label{beta121}\\ 0\le s_i(t,x)\le 1 \quad \text{almost everywhere in } (t,x) \in Q_T,\label{s111}\\ {p_i^\eta} \to p_i \quad \text{almost everywhere in } Q_T, \label{ppppi}\\ \phi\partial_t (\rho_i({p_i^\eta}){s_i^\eta}) \rightharpoonup \phi\partial_t (\rho_i(p_i)s_i) \quad \text{weakly in } L^2(0,T;(H^1_{\Gamma_1}(\Omega))'),\label{dts11}. \end{gather} \end{lemma} \begin{proof} The weak convergence \eqref{ph11}--\eqref{betah11} follow from the uniform estimates \eqref{pp} and \eqref{betabeta} of lemma \ref{lemunifeta}. The lemma \ref{lemcompact} ensures the following strong convergence \begin{gather*} \rho_i({p_i^\eta}){s_i^\eta}\to l_i \text{ in $L^2(Q_T)$ and a. e. in } Q_T, \\ \rho_i({p_i^\eta}){s_i^\eta}\to l_i\text{ in $L^2(\Sigma_T)$ and a. e. in } \Sigma_T, \end{gather*} As the map $\mathbb{H}$ defined in \eqref{defH} is continuous, we deduce \begin{gather*} p^\eta\to p \text{ a. e. in $Q_T$ and a. e. in } \Sigma_T, \\ \beta({s_1^\eta})\to \beta^\star \text{ a. e. in $Q_T$ and a. e. in } \Sigma_T. \end{gather*} The convergence \eqref{pppp} is then established and as $\beta ^{-1}$ is continuous, $${s_1^\eta} \to s_1=\beta^ {-1}(\beta^\star)\text{ a. e. in Q_T and a. e. in } \Sigma_T.$$ From \eqref{maxprin1}, the estimate \eqref{s111} holds and the Lebesgue theorem ensures the strong convergence \eqref{beta121}. The convergence \eqref{ppppi} is a consequence of \eqref{pppp}--\eqref{beta1123}.\\ At last, the weak convergence \eqref{dts11} is a consequence of the estimate \eqref{drhoi}, and the identification of the limit follows from the previous convergence. \end{proof} To achieve the proof of Theorem \ref{main}, it remains to pass to the limit as $\eta$ goes to zero in the formulations \eqref{non-degv1}, for all smooth test functions $\varphi \in C^1([0,T];H^1_{\Gamma_1}(\Omega)) \cap L^2(0,T;H^2(\Omega))$ such that $\varphi(T)=0$ \label{non-degveps} \begin{aligned} &-\int_{Q_T}\phi \rho_i({p_i^\eta}) {s_i^\eta} \partial_t\varphi\,dx\,dt +\int_{Q_T}\mathbf{K} M_i({s_i^\eta})\rho_i({p_i^\eta}) { \nabla} {p_i^\eta} \cdot\nabla \varphi \,dx\,dt \\ &-\int_{Q_T}\mathbf{K} M_i({s_i^\eta})\rho_i^2({p_i^\eta})\mathbf{g} \cdot\nabla \varphi \,dx\,dt +\int_{Q_T}\rho_i({p_i^\eta}){s_i^\eta} f_{P}\varphi \,dx\,dt\\ &- (-1)^i~\eta \int_\Omega \rho_i({p_i^\eta})\nabla({p_1^\eta}-{p_2^\eta}) \cdot\nabla \varphi\,dx\,dt \\ &=\int_{Q_T}\rho_i({p_i^\eta})s_i^If_{I}\varphi \,dx\,dt +\int_{\Omega}\phi \rho_i(p_i^0) s_i^0 \varphi(0,x) \,dx\,dt, \quad i=1,2. \end{aligned} The first term converges due to the strong convergence of $\rho_i({p_i^\eta}) {s_i^\eta}$ to $\rho_i(p_i) s_i$ in $L^2(Q_T)$. The second term can be written, with the help of global pressure, as \label{love} \begin{aligned} &\int_{Q_T}\mathbf{K}M_i({s_i^\eta})\rho_i({p_i^\eta}) { \nabla} {p_i^\eta} \cdot\nabla \varphi \,dx\,dt\\ &= \int_{Q_T}\mathbf{K}M_i({s_i^\eta})\rho_i({p_i^\eta}) { \nabla}p^\eta \cdot\nabla \varphi \,dx\,dt + \int_{Q_T}\mathbf{K}\rho_i({p_i^\eta}) { \nabla} \beta({s_i^\eta}) \cdot\nabla \varphi \,dx\,dt. \end{aligned} The two terms on the right hand side of the equation \eqref{love} converge arguing in two steps. Firstly, the Lebesgue theorem and the convergence \eqref{beta1123}\eqref{ppppi} establish \begin{gather*} \rho_i({p_i^\eta}) M_i({s_i^\eta})\nabla \varphi \to \rho_i(p_i)M_i(s_i)\nabla \varphi \quad \text{strongly in } (L^2(Q_T))^d,\\ \rho_i({p_i^\eta}) \nabla \varphi \to \rho_i(p_i)\nabla \varphi \quad \text{ strongly in } (L^2(Q_T))^d. \end{gather*} Secondly, the weak convergence on global pressure \eqref{ph11} and the weak convergence \eqref{betah11} combined to the above strong convergence allow the convergence for the terms of the right hand side of \eqref{love}. The fifth term can be written as $\eta \int_\Omega \rho_i({p_i^\eta})\nabla({p_1^\eta}-{p_2^\eta}) \cdot\nabla \varphi\,dx\,dt=\sqrt{\eta} \int_\Omega \rho_i({p_i^\eta})(\sqrt{\eta}\nabla f({s_1^\eta})) \cdot\nabla \varphi\,dx\,dt,$ The Cauchy-Schwarz inequality and the uniform estimate \eqref{epsestimate} ensure the convergence of this term to zero as $\eta$ goes to zero. The other terms converge using \eqref{beta1123}\eqref{ppppi} and the Lebesgue dominated convergence theorem. The weak formulations \eqref{degv1} are then established. The main theorem \ref{main} is then established. \end{proof} \begin{thebibliography}{00} \bibitem{A-AD85} H. W. Alt, E. Di Benedetto; \emph{Nonsteady flow of water and oil through inhomogeneous porous media}, Annali della Scuola Normale Superiore di Piza, S\'eries IV, {\textbf XII}, 3, (1985). \bibitem{A-AL83} H. W. Alt, S. Luckhaus; \emph{Quasilinear elliptic-parabolic differential equations}, Math. Z., 3, (1983), pp. 311--341. \bibitem{amaziane1} B. Amaziane, M. Jurak; \emph{A new formulation of immiscible compressible two-phase flow in porous media}, C.R.A.S. Mecanique 336 (7), 600-605, 2008. \bibitem{amaziane2} B. Amaziane, M. El Ossmani; \emph{Convergence analysis of an approximation to miscible fluid flows in porous media by combining mixed finite element and finite volume methods}, Numerical Methods for Partial Differential Equations 24 (3), 799-832, 2008. \bibitem{A-AHZ91} Y. Amirat, K. Hamdache, A. Ziani; \emph{Homogenization of a model of compressible miscible flow in porous media}, Boll. Unione Mat. Ital., VII. Ser., B (1991), pp. 463--487. \bibitem{A-AHZ96} Y. Amirat, K. Hamdache, A. Ziani; \emph{Mathematical analysis for compressible miscible displacement models in porous media,} Math. Models Methods Appl. Sci., 6, no. 6 (1996), pp. 729--747. \bibitem{arbogast} T. Arbogast; \emph{Two-phase incompressible flow in a porous medium with various non homogeneous boundary conditions,} IMA Preprint series 606, February 1990. \bibitem{B-AS79} K. Aziz, A. Settari; \emph{Petroleum reservoir simulation}, Applied Science Publisher LTD, London, 1979. \bibitem{brull} S. Brull; \emph{Two compressible immiscible fluids in porous media: The case where the porosity depends on the pressure.} Advances in Differential equations (2008), Vol 13, No 7-8, 781-800. \bibitem{B-CJ86} G. Chavent, J. Jaffre; \emph{Mathematical models and finite elements for reservoir simulation. Single phase, multiphase and multicomponent flows through porous media,} Studies in Mathematics and its Applications; 17, North-Holland Publishing Comp., 1986. \bibitem{choquet1-08} C. Choquet; \emph{Asymptotic analysis of a nonlinear parabolic problem modelling miscible compressible displacement in porous media,} NoDEA, Nonlinear Differential Equations and Appl. (2008), vol. 15, no. 6, 757--782. \bibitem{choquet2-08} C. Choquet; \emph{Nuclear contamination on a naturally fractured porous medium}, Int. J. of Pure and Applied Math. (2008), vol. 42, no. 2, 281-288. \bibitem{choquet3-08} C. Choquet; \emph{On a fully nonlinear parabolic problem modelling miscible compressible displacement in porous media,} Journal of Mathematical Analysis and Applications, (2008), no. 339, 1112--1133. \bibitem{A-DEH06} F.Z. Da\" \i m, , R. Eymard, D. Hilhorst; \emph{Existence of a solution for two phase flow in porous media : The case that the porosity depends on pressure,} Journal of Mathematical Analysis and Applications, 326(1):332-351, 2007. \bibitem{A-FG00} P. Fabrie, T. Gallou\"et; \emph{Modeling wells in porous media flow,} Math. Models Methods Appl., 10, no 5 (2000), pp. 673--709. \bibitem{A-FL92} P. Fabrie, M. Langlais; \emph{Mathematical analysis of miscible displacement in porous medium}, {SIAM J. Math. Anal.}, 23, no 6 (1992), pp. 1375--1392. \bibitem{A-FLT} P. Fabrie, P. Le Thiez, P. Tardy; \emph{On a system of nonlinear elliptic and degenerate parabolic equations describing compositional water-oil flows in porous media,} Nonlinear Anal., 28, no 9 (1997), pp. 1565--1600. \bibitem{A-FS93} P. Fabrie, M. Saad; \emph{Existence de solutions faibles pour un mod{\'e}le d'{\'e}coulement triphasique en milieu poreux,} Ann. Fac. Sci. Toulouse, 2(3) (1993), pp. 337--373. \bibitem{B-GM95} G. Gagneux, M. Madaune-Tort; \emph{Analyse math\'ematique de mod\eles non lin\'eaires de l'ing\'en\'erie p\'etroli\ere,} Mathematiques and Applications , 22, Springer Verlag, 1995. \bibitem{A-GS04} C. Galusinski, M. Saad; \emph{On a degenerate parabolic system for compressible immiscible two-phase flows in porous media}, Adv. in Diff. Equ., 9(11-12) (2004), pp. 1235--1278. \bibitem{A-GS08} C. Galusinski, M. Saad; \emph{Two compressible immiscible fluids in porous media }, J. Differential Equations 244 (2008), pp. 1741–-1783 \bibitem{P-GS05} C. Galusinski, M. Saad; \emph{Water gas flows in porous media}, Proceedings in fifth AIMS Meeting (Pomona) (2004). \bibitem{A-GS08-2} C. Galusinski, M. Saad; \emph{A nonlinear degenerate system modeling water-gas in reservoir flow}, Discrete and continuous dynamical systems series B, Vol. 9, Num. 2, pp. 281--308, March 2008. \bibitem{B-LS65} O. Ladyzenskaia, V. Solonnikov, N. Uraltseva; \emph{Linear and quasi-linear equations of parabolic type}, Amer. Math. Soc., 1968. \bibitem{lionsj} J. L. Lions; \emph{Quelques m\'ethodes de r\'esolution des probl\emes aux limites non-lin\'eaires}Dunod-Gauthier-Villars, Paris, 1969. \bibitem{lionsmagenes} J. L. Lions, E. Magenes; \emph{Probl\emes aux limites non homog\`nes et applications.} Dunid Paris, 1968. \bibitem{lions} P. L. Lions; \emph{Mathematical topics in fluid mechanics. Vol. 1. Incompressible models.} Oxford Lecture Series in Mathematics and its Applications, 3. Oxford Science Publications. The Clarendon Press, Oxford University Press, New York, 1996. \bibitem{A-SA97} M. Saad; \emph{An accurate numerical algorithm for solving three-phase flow in porous media}, Applicable Analysis, 66 (1997), pp. 57--88. \bibitem{A-SI87} J. Simon; \emph{Compact sets in ${L^{p}(0,T;B)}$}, Ann. Mat. Pura Appl., IV(146) (1987), pp. 65--96. \bibitem{simon} J. Simon; \emph{Nonhomogeneous viscous incompressible fluids: existence of velocity, density, and pressure.} SIAM J. Math. Anal. 21 (1990), no. 5, pp. 1093--1117. \bibitem{Zeidler} E. Zeidler; \emph{Nonlinear Analysis and Fixed-Point Theorems,} Berlin, Springer-Verlag, 1993. \end{thebibliography} \end{document}