\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)$:
\begin{equation}\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),
\end{equation}
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:
\begin{equation}
\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.
\end{equation}
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
\begin{equation}\label{1}
s_{1}(t,x) + s_{2}(t,x) = 1.
\end{equation}
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},
\begin{equation}\label{cpl}
f(s_1(t,x)) = p_{1}(t,x) - p_{2}(t,x).
\end{equation}
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
\begin{equation}\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),
\end{equation}
the global pressure is then defined as in \cite{B-CJ86}
\begin{equation}\label{defp}
p = p_{2}+\tilde{p}(s_1)= p_{1}+\bar{p}(s_1)
\end{equation}
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
\begin{equation} \label{defbeta}
\beta(s)=\int_0^s \alpha(\xi)d\xi.
\end{equation}
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.
\begin{equation} \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}
\end{equation}
where $\mathbf{n}$ is the outward normal to the boundary $\Gamma_{\rm imp}$.
The initial conditions are defined on pressures
\begin{equation}\label{ic}
 p_i(0,x) = p_i^0(x) \quad\mbox{in } \Omega, \quad i=1,2.
\end{equation}

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 $0<s< 1$, and   $\alpha(0)=\alpha(1)=0$.

\end{enumerate}
    We assume that $\beta ^{-1}$, inverse of $\beta(s):=\int_0^s\alpha(z)dz$,  is an  H\"older function
    of order $\theta$, with $0<\theta \le 1$, on $[0,\beta(1)]$.
    Which means there exists a positive $c$ such that for all $s_1,s_2 \in [0,\beta(1)]$, one has
    $|\beta ^{-1}(s_1)-\beta ^{-1}(s_1)|\le c |s_1-s_2|^\theta$.

Assumptions (H1)--(H7) are classical for
porous media.

The main existence result of this paper is given below, for that
let us define the following Sobolev space
$$
H^1_{\Gamma_1}(\Omega)= \{\, u\in H^1(\Omega); u=0 \text{ on } \Gamma_1\},
$$
this is an Hilbert space when equipped with the norm
$\| u\|_{H^1_{\Gamma_1}(\Omega)} =\|\nabla u\|_{(L^2(\Omega))^d}$.
Let us state the main results of this paper.


\begin{theorem}\label{main}
Let {\rm (H1)-(H7)} hold.
Let $(p_1^0,\, p_2^0)$ belongs to $L^2(\Omega)\times L^2(\Omega)$.
Then  there exists $(p_1,p_2)$ satisfying
\begin{gather}
 p_i \in L^2(0,T;H^1_{\Gamma_1}(\Omega)),
 \phi\partial_{t}(\rho_i(p_i)s_i)\in L^2(0,T;(H^1_{\Gamma_1}(\Omega))'),\;
  i=1,2,\label{mainc1}\\
 0\leq s_i(t,x)\leq1 \quad \text{ a.e. in } Q_T,\;
 i=1,2,; \beta(s_1) \in L^2(0,T;H^1(\Omega)) \label{mainc2}
\end{gather}
such that for all $\varphi_i \in C^1(0,T;H^1_{\Gamma_1}(\Omega))$
with $\varphi_i(T)=0$,
\begin{equation} \label{degv1}
\begin{aligned}
& -\int_{Q_T}  \phi \rho_i(p_i) s_i \partial_t\varphi_i\,dx\,dt
-\int_\Omega \phi (x) \rho_i(p_i^0(x)) s_i^0(x)\varphi_i(0,x)\,dx \\
&+\int_{Q_T}\mathbf{K}M_i(s_i)\rho_i(p_2) { \nabla} p_i
\cdot\nabla \varphi_i \,dx\,dt
-\int_{Q_T}\mathbf{K}M_i(s_i)\rho_i^2(p_i)\mathbf{g} \cdot\nabla \varphi_i \,dx\,dt\\
&+\int_{Q_T}\rho_i(p_i)s_i f_{P}\varphi_i \,dx\,dt\\
&=\int_{Q_T}\rho_i(p_i)s_i^If_{I}\varphi_i \,dx\,dt,
\end{aligned}
\end{equation}
and finally the initial conditions are satisfied in a weak sense
as follows: For all $\psi\in H^1_{\Gamma_1}(\Omega)$  the function
\begin{equation} \label{intlcdwksense}
t\to  \int_\Omega\phi \rho_i(p_i) s_i \psi\,dx \in {\mathcal C}^0([0,T]),
\end{equation}
furthermore,
\begin{equation}
\label{intlcdwksense1}
\Big(\int_\Omega \phi \rho_i(p_i) s_i \psi\,dx\Big)(0)
 = \int_\Omega\phi \rho_i(p_i^0) s_i^0\psi\,dx.
\end{equation}
\end{theorem}

As we can see, the above notion of weak solutions is very natural
provided that we explain the origin of the requirements
\eqref{mainc1}--\eqref{mainc2}.
Obviously, they correspond to {\it a priori} estimates.
Indeed, the equations \eqref{degv1} ensure that $s_i\geq 0$ ($i=1,2$)
which is equivalent to $0\le s_i\le 1$ (the proof is detailed in
lemma \ref{max}.
The key point is to obtain the estimates on  $\nabla p$ and
$\nabla \beta (s_1)$.
For that, define
\begin{gather}
g_i(p_i):=\int_0^{p_i} \frac{1}{\rho_i(\xi)} \, d \xi,\quad
 i=1,2, \label{gipi}\\
\mathcal{H}_i(p_i):=\rho_i(p_i)g_i(p_i)-p_i, \quad i=1,2\label{hipi},
\end{gather}
then $\mathcal{H}_i'(p_i)=\rho_i'(p_i) g_i(p_i)$, $\mathcal{H}_i(0)=0$,
$\mathcal{H}_i(p_i)\ge 0$ for all $p_i$, and $\mathcal{H}_i$ is sublinear with respect
to $p_i$.
Multiplying \eqref{conslaw} by $g_1(p_1)$ for $i=1$ and
\eqref{conslaw} by $g_2(p_2)$ for $i=2$ then integrate the
equations with respect to $x$ and adding them, we deduce at least
formally,
\begin{equation}
\begin{aligned}
&\frac{d}{dt}\int_\Omega  \phi \Big( s_1 \mathcal{H}_1(p_1) +s_2 \mathcal{H}_2(p_2)
 + \int_{0}^{s_1}f(\xi)\, d\xi \Big)dx \\
&+\int_\Omega \mathbf{K} M_1(s_1) \nabla p_1 \cdot \nabla p_1 \, dx
+\int_\Omega \mathbf{K} M_2(s_2) \nabla p_2 \cdot \nabla p_2 \,dx\\
&-\int_\Omega \mathbf{K} M_1(s_1) \rho_1(p_1)\mathbf{g} \cdot \nabla
p_1 \, dx
-\int_\Omega \mathbf{K} M_2(s_2) \rho_2(p_2)\mathbf{g} \cdot \nabla p_2 \, dx
\\
&+ \int_\Omega \rho_1(p_1) s_1f_p g_1(p_1) \, dx
+\int_\Omega \rho_2(p_2) s_2f_p g_2(p_2) \, dx\\
&=\int_\Omega \rho_1(p_1) s_1^If_I g_1(p_1) \, dx
+\int_\Omega \rho_2(p_2) s_2^If_I g_2(p_2) \, dx.
\end{aligned} \label{theesti}
\end{equation}
A key point  is to obtain formally the first term in the above equality,
for that let
\begin{align*}
D&=\partial_t(\rho_1(p_1)s_1)g_1(p_1)+\partial_t(\rho_2(p_2)s_2)g_2(p_2)\\
&= \partial_t(\rho_1(p_1)s_1g_1(p_1))+\partial_t(\rho_2(p_2)s_2g_2(p_2)) -s_1\partial_tp_1 -s_2\partial_t p_2.
 \end{align*}
We have $s_1+s_2=1$, then
$s_1\partial_tp_1 +s_2\partial_tp_2 = s_1\partial_t f(s_1)+\partial_tp_2
= \partial_tG(s_1)+\partial_tp_2$, where $G$
is a primitive of $s_1f'(s_1)$. We can write $D$ as
$D=\partial_tE$ where $E$ is defined by
\begin{align*}
E&= \rho_1(p_1)s_1g_1(p_1)+\rho_2(p_2)s_2g_2(p_2) -G(s_1)-p_2\\
&=s_1(\rho_1(p_1)g_1(p_1)-p_1)+s_2(\rho_2(p_2)s_2g_2(p_2)-p_2)-G(s_1)
+s_1f(s_1),
 \end{align*}
from the definition of the functions $\mathcal{H}_i$ ($i=1,2)$ and $G$,
the expression of $E$ is equivalent to:
$$
E=s_1\mathcal{H}_1(p_1)+s_2\mathcal{H}_2(p_2)+\int_0^{s_1}f(\xi)\,d\xi.
$$

Using the assumptions (H1)--(H6) and the
fact that $\mathcal{H}_i \ge 0$, $g_i(p_i)$ is sublinear with respect to
$p_i$ we deduce from \eqref{theesti} that
\begin{gather}\label{magicineq}
\int_{Q_T}  M_1(s_1) \nabla p_1 \cdot \nabla p_1 \, dx
+\int_{Q_T}  M_2(s_2) \nabla p_2 \cdot \nabla p_2 \,dx <\infty,\\
\nabla p=\nabla p_2 + \frac{M_1}{M} \nabla f(s_1)
 =\nabla p_1 - \frac{M_2}{M} \nabla f(s_1), \label{global1}
\end{gather}
 then, we deduce a magic equality
\begin{equation}\label{magic11}
\begin{aligned}
&\int_{Q_T} M |\nabla p|^2\,dx
 + \int_{Q_T} \frac{M_1 M_2}{M}|\nabla f(s_1)|^2 \, dx\\
&=\int_{Q_T}  M_1(s_1) |\nabla p_1|^2 \, dx
+\int_{Q_T}  M_2(s_2) |\nabla p_2|^2 \,dx,
\end{aligned}
\end{equation}
thus,  the equality \eqref{magic11} and the assumption
(H3)  ensure that $p\in L^2(0,T;H^1_{\Gamma_1}(\Omega))$ and
$\beta (s_1) \in L^2(0,T;H^1(\Omega))$.

Before establishing theorem \ref{main}, we introduce the existence
of solutions  to system \eqref{conslaw} under the assumptions
(H1)--(H7), with the addition of some
terms on each equation to save a maximum principle on saturations,
to conserve the  existence of solutions of a time discretization,
and to insure a compactness {\it lemma} which is necessary to
pass from an elliptic approximation to the original
parabolic systemfrom, after that we get rid of these terms by
a limit process which is also conserved. Thus we consider the
non-degenerate system:
\begin{gather} \label{conslawnondeg1}
\begin{aligned}
&\phi\partial_{t}( \rho_{1}({p_1^\eta}){s_1^\eta}) - \operatorname{div}
(\mathbf{K}\rho_{1}({p_1^\eta})M_1({s_1^\eta})\nabla {p_1^\eta})
+\operatorname{div} (\mathbf{K}\rho_{1}^2({p_1^\eta})M_1({s_1^\eta})\mathbf{g})\\
&- \eta \operatorname{div} (\rho_1({p_1^\eta}) \nabla({p_1^\eta}-{p_2^\eta}))
+\rho_1({p_1^\eta}) {s_1^\eta}  f_{P} = \rho_1({p_1^\eta}) s^I_1 f_{I},
\end{aligned}\\
\label{conslawnondeg2}
\begin{aligned}
&\phi\partial_{t}( \rho_{2}({p_2^\eta}){s_2^\eta}) - \operatorname{div}
(\mathbf{K}\rho_{2}({p_2^\eta})M_2({s_2^\eta})\nabla {p_2^\eta})
+\operatorname{div} (\mathbf{K}\rho_{2}^2({p_2^\eta})M_2({s_2^\eta})\mathbf{g})\\
&- \eta \operatorname{div} (\rho_2({p_2^\eta})\nabla({p_2^\eta}-{p_1^\eta}))
+\rho_2({p_2^\eta}) {s_2^\eta}  f_{P} = \rho_2({p_2^\eta}) s^I_2 f_{I},
\end{aligned}
\end{gather}
completed with the initial conditions \eqref{ic}, and the following
mixed boundary conditions, for $i=1,2$,
\begin{equation} \label{bcnondeg}
\begin{gathered}
 {p_1^\eta}(t,x) = 0, \quad {p_2^\eta}(t,x) = 0 \quad \text{on } \Gamma_{1},\\
 \Big(  \mathbf{K}M_i({s_i^\eta})( \nabla
{p_i^\eta}-\rho_i({p_i^\eta})\mathbf{g})+(-1)^{i+1}\eta \nabla({p_1^\eta}
-{p_2^\eta})\Big)\cdot\mathbf{n} =  0 \quad \text{on }\Gamma_{\rm imp},
\end{gathered}
\end{equation}
where $\mathbf{n}$ is the outward normal to the boundary
 $\Gamma_{\rm imp}$.

Now, we state existence of solutions of the above system  by
the following theorem.

\begin{theorem}[Non-degenerate system] \label{non-degt}
Let {\rm (H1)--(H6)} hold.
 Let $(p_1^0,\, p_2^0)$ belongs to $L^2(\Omega)\times L^2(\Omega)$.
Then for all $\eta>0$, 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$,
\begin{equation}\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}
\end{equation}
 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$,
\begin{equation}\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{equation}
\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,
\begin{equation} \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}
\end{equation}
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,
\begin{equation}
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}
\end{equation}


 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:
\begin{equation}\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}
\end{equation}
\begin{equation}\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}
\end{equation}
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},
\begin{equation} \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}
\end{equation}
we deduce from the Cauchy-Schwarz inequality that \eqref{ad} reduces to,
\begin{equation} \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}
\end{equation}
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
\begin{equation} \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}
\end{equation}
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)$,
\begin{equation}
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}
\end{equation}
 Furthermore and due to the convergence \eqref{cvp1nw},
it follows that
\begin{equation}
\nabla p_1^n \to  \nabla p_1 \quad \text{weakly in }
L^2(\Omega), \label{cvterm2bis}
\end{equation}
 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,
\begin{equation}
\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}
\end{equation}
 Furthermore $\overline{p}_{i,n}$ converges in $L^2(\Omega)$, it follows that
\begin{equation}
\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}
\end{equation}

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
\begin{equation} \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}
\end{equation}
\begin{equation} \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}
\end{equation}
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{equation}
\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}
\end{equation}
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
\begin{equation} \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}
\end{equation}
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
\begin{equation}\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}
\end{equation}
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
\begin{equation}\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}
\end{equation}
\begin{equation}
\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}
\end{equation}
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
\begin{equation} \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}
\end{equation}
 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{equation}
\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}
\end{equation}
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
\begin{equation} \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}
\end{equation}
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
\begin{equation} \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}
\end{equation}
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,
\begin{equation}\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}
\end{equation}
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
\begin{equation} \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}
\end{equation}
\begin{equation} \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}
\end{equation}
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{equation}
\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}
\end{equation}
\begin{equation}
\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}
\end{equation}
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{equation}
\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}
\end{equation}
 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$,
\begin{equation} \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}
\end{equation}
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
\begin{equation}\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).
\end{equation}
 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
\begin{equation} \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}
\end{equation}
Since the function $\mathcal{F}$ is convex,
\begin{equation}\label{convexeq}
(s_1-s_1^\star)f(s_1)\ge \mathcal{F} (s_1)-\mathcal{F} (s_1^\star).
\end{equation}
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
\begin{equation} \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}
\end{equation}
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
\begin{equation} \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}
\end{equation}
and
\begin{equation}
\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].
\end{equation}
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$,
\begin{equation} \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}
\end{equation}
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
\begin{equation}\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.
\end{equation}

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{equation}
\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}
\end{equation}
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
\begin{equation} \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{equation}
\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
\begin{equation} \label{cvrtld}
\tilde r_{i,h} \to  r_i \text{ strongly in } L^2(Q_T).
\end{equation}
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
\begin{equation}
H(u,v) = (p_{1,h},s_{1,h})
\label{defHb}
\end{equation}
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
\begin{equation} \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}
\end{equation}
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
\begin{equation} \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}
\end{equation}
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
\begin{equation} \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}
\end{equation}
Integrate \eqref{estequ} over$(0,T)$ to obtain
\begin{equation} \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}
\end{equation}
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
\begin{equation}\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}
\end{equation}
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,
\begin{equation} \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}
\end{equation}
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
\begin{equation} \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}
\end{equation}
The Young inequality permits to absorb the last term  by the
 first term on the left hand side of \eqref{magic2eta} to  obtain
\begin{equation}\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}
\end{equation}
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
\begin{equation}
\mathbb{H}(u,v) = (p,\beta(s_1)) \label{defH}
\end{equation}
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<r<\infty$, let us denote
$$
W^{\tau,r}(\Omega) =\{w \in L^r(\Omega) ;
\int_\Omega \int_\Omega \frac{\vert w(x)-w(y)\vert^r}
{\vert x-y\vert^{d+\tau r}} \,dx dy < +\infty \}
$$
equipped with the norm
$$
\|w\|_{W^{\tau,r}(\Omega)} = \Big( \|w\|_{L^r(\Omega)}^r
+ \int_\Omega \int_\Omega \frac{\vert w(x)-w(y)\vert^r}
{\vert x-y\vert^{d+\tau r}} \,dx dy \Big)^{1/r},
$$
recall that $d$ denotes the space dimension.
Let $q$, $\sigma$ be in  $W^{\tau,r}(\Omega)\times W^{\tau,r}(\Omega)$, then the
H\"older functions $u$ and $v$ belong to
$W^{\theta\tau,r/\theta}(\Omega)$.  In fact,
\begin{align*}
|u(q,\sigma)| \le c_1 |q |^\theta + c_2 |\sigma|^\theta;
\end{align*}
then $u$ belongs to $L^{r/\theta}(\Omega)$. Furthermore,
\begin{align*}
&\int_\Omega \int_\Omega \frac{\vert u(q(x),\sigma(x))-u(q(y),
\sigma(y))\vert^{r/\theta}}
{\vert x-y\vert^{d+\tau r}} \,dx\,dy\\
&\le c_1^{r/\theta}\int_\Omega \int_\Omega \frac{\vert q(x)-q(y) \vert^{r}}
{\vert x-y\vert^{d+\tau r}} \,dx dy
+ c_2^{r/\theta} \int_\Omega \int_\Omega \frac{\vert \sigma(x)
-\sigma(y) \vert^{r}}
{\vert x-y\vert^{d+\tau r}} \,dx dy ,
\end{align*}
which ensures,
$$
\|u(q,\sigma)\|_{W^{\theta\tau,r/\theta}(\Omega)} \le c(
\|q\|_{W^{\tau,r}(\Omega)}^\theta + \|\sigma\|_{W^{\tau,r}(\Omega)}^\theta).
$$
Using the continuity of the injection of
$H^1(\Omega)$ into $W^{\tau,2}(\Omega)$, with $\tau <1$,
\begin{align*}
\|u(p,\beta(s_1))\|_{W^{\theta\tau,2/\theta}(\Omega)}
&\le c( \|p\|_{W^{\tau,2}(\Omega)}^\theta
 + \| \beta(s_1)\|_{W^{\tau,r}(\Omega)}^\theta)\\
&\le  c(\|p\|_{H^1(\Omega)}^\theta + \| \beta(s_1)\|_{H^1(\Omega)}^\theta)
\end{align*}
integrating the above inequality over $(0,T)$,
$$
\|u(p,\beta(s_1))\|_{L^{2/\theta}(0,T;W^{\theta\tau,2/\theta}(\Omega))}
\le c  \|p\|_{L^2(0,T;H^1(\Omega))}^\theta
+  \|\beta(s_1)\|_{L^2(0,T;H^1(\Omega))}^\theta
$$
Furthermore the porosity function $\phi$ belongs to
$W^{1,\infty}(\Omega)$, it follows that
$$
\|\phi u(p,\beta(s_1))\|_{L^{2/\theta}(0,T;W^{\theta\tau,2/\theta}(\Omega))}
 \le C.
$$
As $\Omega$ is bounded and regular, we have, for $\tau'<\theta\tau$,
$$
W^{\theta\tau,2/\theta}(\Omega)\subset W^{\tau',2/\theta}(\Omega)
\subset (H_{\Gamma_1}^1(\Omega))'
$$
with compact injection from
$ W^{\theta\tau,2/\theta}(\Omega)$  into $W^{\tau',2/\theta}(\Omega)$.
 Finally, from a standard compactness argument, we get
$E_M$ is relatively compact in
$L^{2/\theta}(0,T;W^{\tau',2/\theta}(\Omega))\subset L^{2}(0,T;L^2(\Omega))$.
Secondly, the trace operator $\gamma$ maps continuously
$W^{\tau',2/\theta}(\Omega))$ into
$W^{\tau'-\theta/2,2/\theta}(\Gamma))$ as soon as
$\tau' > \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$
\begin{equation} \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}
\end{equation}
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
\begin{equation} \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}
\end{equation}
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}