\documentclass[reqno]{amsart} \usepackage{hyperref} \AtBeginDocument{{\noindent\small \emph{Electronic Journal of Differential Equations}, Vol. 2012 (2012), No. 102, pp. 1--22.\newline ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu \newline ftp ejde.math.txstate.edu} \thanks{\copyright 2012 Texas State University - San Marcos.} \vspace{9mm}} \begin{document} \title[\hfilneg EJDE-2012/102\hfil Existence for a global pressure formulation] {Existence for a global pressure formulation of water-gas flow in porous media} \author[B. Amaziane, M. Jurak, A. Vrba\v ski\hfil EJDE-2012/102\hfilneg] {Brahim Amaziane, Mladen Jurak, Anja Vrba\v ski} % in alphabetical order \address{Brahim Amaziane \newline Univ Pau \& Pays Adour, Laboratoire de Math\'ematiques et de leurs Applications, CNRS-UMR 5142, Av. de l'Universit\'e, 64000 Pau, France} \email{brahim.amaziane@univ-pau.fr} \address{Mladen Jurak \newline Faculty of Science, University of Zagreb, Bijenicka 30, 10000 Zagreb, Croatia} \email{jurak@math.hr} \address{Anja Vrba\v ski \newline Faculty of Science, University of Zagreb, Bijenicka 30, 10000 Zagreb, Croatia} \email{avrbaski@math.hr} \thanks{Submitted June 5, 2012. Published June 18, 2012.} \subjclass[2000]{35K55, 35K65, 76S05} \keywords{Degenerate system; global pressure; immiscible compressible; \hfill\break\indent nonlinear parabolic system; nuclear waste; porous media; two-phase flow, water-gas} \begin{abstract} We consider a model of water-gas flow in porous media with an incompressible water phase and a compressible gas phase. Such models appear in gas migration through engineered and geological barriers for a deep repository for radioactive waste. The main feature of this model is the introduction of a new global pressure and it is fully equivalent to the original equations. The system is written in a fractional flow formulation as a degenerate parabolic system with the global pressure and the saturation potential as the main unknowns. The major difficulties related to this model are in the nonlinear degenerate structure of the equations, as well as in the coupling in the system. Under some realistic assumptions on the data, including unbounded capillary pressure function and non-homogeneous boundary conditions, we prove the existence of weak solutions of the system. Furthermore, it is shown that the weak solution has certain desired properties, such as positivity of the saturation. The result is proved with the help of an appropriate regularization and a time discretization of the coupled system. We use suitable test functions to obtain a priori estimates and a compactness result in order to pass to the limit in nonlinear terms. \end{abstract} \maketitle \numberwithin{equation}{section} \newtheorem{theorem}{Theorem}[section] \newtheorem{lemma}[theorem]{Lemma} \newtheorem{proposition}[theorem]{Proposition} \newtheorem{remark}[theorem]{Remark} \allowdisplaybreaks \section{Introduction} Two-phase flow in porous media is important to many practical problems, including those in petroleum reservoir engineering, unsaturated zone hydrology, and soil sciences. Most recently, modeling multiphase flow received an increasing attention in connection with the disposal of radioactive waste and sequestration of CO$_2$. This paper focuses on the modeling and analysis of immiscible compressible two-phase flow through porous media, in the framework of the geological disposal of radioactive waste. As a matter of fact, one of the solutions envisaged for managing waste produced by nuclear industry is to dispose of it in deep geological formations chosen for their ability to prevent and attenuate possible releases of radionuclides in the geosphere. In the frame of designing nuclear waste geological repositories appears a problem of possible two-phase flow of water and gas, mainly hydrogen. For more details, see for instance \cite{Norris,OECD-NEA}. Multiple recent studies have established that in such installations important amounts of gases are expected to be produced in particular due to the corrosion of metallic components used in the repository design, see e.g. \cite{CroiseEtAl,SengerEtAl} and the references therein. The French Agency for the Management of Radioactive Waste (Andra) \cite{ANDRA} is currently investigating the feasibility of a deep geological disposal of radioactive waste in an argillaceous formation. A question related to the long-term performance of the repository concerns the impact of the hydrogen gas generated in the wastes on the pressure and saturation fields in the repository and the host rock. As reported in \cite{AJZ2}, the mathematical analysis of two-phase flow in porous media has been a problem of interest for many years and many methods have been developed. There is an extensive literature on this subject. We will not attempt a literature review here, but merely mention a few references. Here we restrict ourselves to two-phase flow in porous media. We refer for instance to \cite{AB,ant-kaz-mon1990,Arb92,BA-HA-95,GC-JJ,Chen01,DFZ-ER-HD-07,FBK-07,GG-MT-96,KD-LS-84,LS,YLM-02,YLM-06} for more information on the analysis, especially on the existence of solutions, of immiscible incompressible two-phase flow in porous media, and to \cite{AY-HK-ZA-96,AY-MM-95,AY-SV-07,CC-08,FX-94} for miscible compressible flow in porous media. However, the situation is quite different for immiscible compressible two-phase flow in porous media, where, only recently few results have been obtained. In the case of immiscible two-phase flows with one (or more) compressible fluids without any exchange between the phases, some approximate models were studied in \cite{GS1,GS2,GS3,GS4}. Namely, in \cite{GS1} certain terms related to the compressibility are neglected, and in \cite{GS2,GS3,GS4} the mass densities are assumed to depend not on the physical pressure, but on the global pressure. As shown in \cite{AJZ1} the models based on the mass density approximation can be suitable in oil reservoir simulations but are inadequate in many underground gas and water flows where the difference between the phase pressures and the global pressure can be significant. In the articles \cite{GS5,KS1,KS2,KS3}, a more general immiscible compressible two-phase flow model in porous media is considered for homogeneous fields. In these contributions, the models are based on phase formulations, i.e. the main unknowns are the phase pressures and the saturation of one phase, and the feature of the global pressure as introduced in \cite{ant-kaz-mon1990,GC-JJ} for incompressible immiscible flows is used to establish a priori estimates. The obtained results are established under the assumption that the capillary pressure is bounded and no discontinuity of the porosity and the permeability is permitted, which is too restrictive for some realistic problems, such as gas migration through engineered and geological barriers for a deep repository for radioactive waste. In the case of immiscible two-phase flows with one (or more) compressible fluids with exchange between the phases, i.e. a multicomponent model, existence of weak solutions to these equations under some assumptions on the compressibility of the fluids has been recently established in \cite{CF-SB-SM,LB-SW,MA-09,Smai1,Smai2}. For modeling such flow problems, there are two main approaches known as the phase and the global pressure formulations. The phase formulation is based on individual balance equations for each of the fluids. For such formulation, in regions without the wetting fluid, the wetting pressure is physically not well-defined. So the pressures are not mathematically well defined globally. Also, as a consequence of the degeneracy of the relative permeability functions is that no uniform bounds for the pressure gradients in $L^2$-spaces are available. To overcome these difficulties, the global pressure formulation of the original flow equations was introduced for incompressible two-phase flows in \cite{ant-kaz-mon1990,GC-JJ}, and generalized recently to compressible two-phase flows in \cite{AJ,AJZ1,AZK} and for three-phase flows in \cite{CG-09}. In this paper, we focus our attention on the study of immiscible, compressible two-phase flow in porous media under isothermal condition without mass transfer between the phases taking into account gravity, capillary effects, and heterogeneity. The system consists of incompressible wetting phase and compressible non-wetting phase, such as water and hydrogen in the context of gas migration through engineered and geological barriers for a deep repository of nuclear waste. The system is written in a fractional flow formulation with the saturation of one of the phases and the global pressure as independent variables. This new formulation, recently derived in \cite{AJ,AJZ1} without any simplifying assumptions, is fully equivalent to the original phase equations formulation and is more suitable for mathematical and numerical analysis, for more details see \cite{AJZ2,AJZ3}. The fractional flow approach treats the two-phase flow problem as a total fluid flow of a single mixed fluid, and then describes the individual phases as fractions of the total flow. This approach leads to a less strong coupling between the two coupled equations: the global pressure equation and the saturation equation. The mathematical structure is well defined: the system consists of two nonlinear degenerate parabolic equations. The major difficulties related to this model are in the nonlinear degenerate structure of the equations, as well as in the coupling in the system. This formulation leads to a coupled system consisting of a nonlinear parabolic equation for the global pressure and a nonlinear degenerate parabolic diffusion-convection equation for the saturation, subject to appropriate boundary and initial conditions. Our aim is to establish existence of weak solutions for this system of equations under realistic assumptions. Let us mention that the case of two compressible fluids was treated in \cite{AJZ2}. In contrast to this case, the incompressibility of one phase leads to additional difficulties in the proof of the existence result for the system under consideration. Namely, the assumption (A.6) on mass densities in \cite{AJZ2} is not valid for the model under consideration. Although we follow the strategy used in \cite{AJZ2}, that is we first regularize our model and then using the discretization in time, apply the fixed point arguments, still the presence of incompressible phase brings additional difficulties in obtaining a priori estimates and passage to the limit, and makes the proof essentially more involved. Our approach also relies on the compactness result proved in our previous work \cite{AJZ2}. Thus we extend the results of \cite{AJZ2} in the case of water-gas flow in porous media with an incompressible water phase and a compressible gas phase. The rest of the paper is organized as follows. In Section~\ref{Formulation_main} we introduce the notations and formulate the mathematical and physical model under consideration. Then we state the assumptions on the data and present the main result of the paper on the existence of a weak solution of the problem. This result is proved in three steps. In Section~\ref{Regul} we regularize the system, introducing a small regularization parameter $\eta>0$ and state the existence result for a weak solution of the regularized problem. Section~\ref{Discrete} provides a construction of the approximate solutions to the regularized system by replacing the time derivatives with finite differences with a small time step $h>0$ and the existence result for the corresponding system, as well as a maximum principle for the saturation. In Section~\ref{Uniform_h} we establish uniform estimates with respect to $h$ and $\eta$ using suitable test functions. This allows us to pass to the limit when $h\to 0$ which gives the existence of a weak solution for the regularized problem; this is performed in Section~\ref{Limit_h}. Finally, in Section~\ref{Limit_eta} we pass to the limit as $\eta\to 0$ using an adapted compactness result, as in \cite{Simon,KS3,AJZ2}, and prove the existence of the solutions of the problem defined in Section~\ref{Formulation_main}. \section{Problem formulation and the main result}\label{Formulation_main} The model of water-gas flow to be presented herein is formulated in terms of the non-wetting (gas) phase saturation and the global pressure and it is developed in \cite{AJ,AJZ1}. The saturations of the wetting and the non-wetting phases are denoted by $S_w$ and $S:=S_g = 1 - S_w$, respectively, and $\lambda_j = \lambda_j(S)$ stands for the relative mobility of the $j$-phase, $j\in\{w,g\}$. The pressures and the mass densities of the wetting and the non-wetting phases are denoted by $P_w$, $P_g$ and $\rho_w$, $\rho_g$, respectively. The wetting phase (water) is assumed incompressible ($\rho_w =$ const.) and the non-wetting phase (gas) is compressible, $\rho_g=\rho_g(P_g)$. The fully equivalent global pressure formulation of immiscible, compressible two-phase flow introduced in \cite{AJ, AJZ1} is defined in terms of the global pressure $P$ and the saturation potential $\theta$ defined by \begin{align} \theta = \beta(S) = \int_0^{S}\sqrt{\lambda_w(s)\lambda_g(s)}P_c'(s)ds,\label{beta-def} \end{align} where $P_c(S) = P_g - P_w$ is the capillary pressure function. The global pressure $P$ is a pressure-like variable which allows to express the phase pressures $P_w$, $P_g$ as functions of (gas) saturation and the global pressure: $P_w = P_w(S,P)$, $P_g = P_g(S,P)$. The non-wetting phase mass density will be expressed, therefore, as a function of $S$ and $P$, $\rho_g = \rho_g(P_g(S,P)) =: \rho_g(S,P)$. The phase pressures are then given by (see \cite{AJ}): \begin{gather} P_g(S, P) = P + P_c(0) + \int_0^{S}f_w(s, P)P_c'(s)ds,\label{pn_eqn}\\ P_w(S, P) = P_g(S, P) - P_c(S),\label{pw_eqn} \end{gather} where the fractional flow functions are defined by \[ {f}_{w}(S,P) = {\rho_{w}\lambda_{w} (S)}/{\lambda(S,P)}, \quad {f}_{g}(S,P) = {\rho_{g}(S,P)\lambda_{g} (S)}/{\lambda(S,P)}. \] Also, we use the total mobility $\lambda(S,P),$ defined by \[ \lambda(S,P) = \rho_w\lambda_w(S) + \rho_g(S,P)\lambda_g(S). \] %\label{total-mobility} The water-gas flow equations given in a fully equivalent global pressure formulation are described by the following equations (see \cite{AJ}): \begin{gather} \begin{aligned} &-\rho_w\Phi\frac{\partial S}{\partial t} -\operatorname{div}(\Lambda_w(S,P)\mathbb{K}\nabla P)+\operatorname{div}(A(S,P)\mathbb{K}\nabla\theta ) \\ & +\rho_w^2\operatorname{div}(\lambda_w(S) \mathbb{K} \mathbf{g}) % + \rho_w {f_w}(S,P) F_P = F_w , \end{aligned}\label{eqn-w} \\ \begin{aligned} &\Phi \frac{\partial }{\partial t}( \rho_g(S,P) S ) -\operatorname{div}(\Lambda_g(S,P) \mathbb{K}\nabla P) - \operatorname{div}(A(S,P) \mathbb{K}\nabla\theta ) \\ & +\operatorname{div}(\lambda_g(S)\rho_g(S,P)^2 \mathbb{K} \mathbf{g}) = F_g , \end{aligned}\label{eqn-g} \end{gather} where $\Phi(x)$ is the porosity, $\mathbb{K}(x)$ is the absolute permeability tensor of the porous medium, $F_w$, $F_g$ are known source terms and the gravity vector is denoted by $\textbf{g}$. The coefficient $A(S,P)$ is given by \begin{equation} A(S,P) = \rho_w\rho_g(S,P)\frac{\sqrt{\lambda_w(S)\lambda_g(S)}}{\lambda(S,P)} \label{mean-density} \end{equation} and the mobility functions $\Lambda_w, \Lambda_g$ are given by \[ \Lambda_{w}(S,P) = \rho_{w}\lambda_{w}(S) \omega(S,P),\quad \Lambda_{g}(S,P) = \rho_{g}(S,P)\lambda_{g}(S) \omega(S,P), \] where the function $\omega(S, P)$ is defined by (see \cite{AJ,AJZ1}) \[ \omega(S, P) = \frac{\partial P_w(S,P)}{\partial P} = \frac{\partial P_g(S,P)}{\partial P}. \] Let a porous domain $\Omega\subset\mathbb{R}^d$, $ d=1,2,3$, be a bounded, connected, Lipschitz domain. The domain boundary is considered to be decomposed as $\partial\Omega = \Gamma_{D} \cup \Gamma_{N}$. The time interval of interest is $]0, T[$ and $Q = \Omega \times ]0, T[$, $\Gamma_{i}^T = \Gamma_{i}\times ]0,T[$, $i \in \{D, N\}$. We impose the boundary conditions for this system as follows: \begin{gather} \theta = \theta_D, \quad P = P_D \quad\text{on } \Gamma_{D}^T, \label{bc-1}\\ \mathbf{Q}_w\cdot {\bf n} = G_w, \quad \mathbf{Q}_n\cdot {\bf n} = G_g \quad\text{on } \Gamma_{N}^T. \label{bc-2} \end{gather} Here $P_D$, $\theta_D$, $G_w$ and $G_g$ are given functions, ${\bf n}$ is the outward unit normal to $\partial\Omega$ and \begin{gather*} \mathbf{Q}_w = \rho_w\mathbf{q}_w = -\Lambda_w(S,P) \mathbb{K}\nabla P + A(S,P) \mathbb{K} \nabla\theta + \rho_w^2\lambda_w(S) \mathbb{K} \mathbf{g},\\ \mathbf{Q}_g =\rho_g(P_g) \mathbf{q}_g = -\Lambda_g(S,P)\mathbb{K}\nabla P - A(S,P)\mathbb{K} \nabla\theta + \rho_g(S,P)^2\lambda_g(S) \mathbb{K} \mathbf{g} \end{gather*} are the phase mass fluxes with $\mathbf{q}_j$ being the volumetric velocity of the $j$-phase, $j\in\{w,g\}$. The Dirichlet boundary data $P_D$, $\theta_D$ are assumed to be extended to the whole domain $Q$; to express their regularity we introduce the space \[ W = \{\varphi \in L^2(0,T;H^1(\Omega)) \colon \varphi \in L^\infty(0,T;L^1(\Omega)), \partial_t \varphi \in L^1(Q)\} \] with the norm \[ |||\varphi||| = \|\varphi\|_{L^2(0,T;H^1(\Omega))} + \|\varphi\|_{L^\infty(0,T;L^1(\Omega))} + \|\partial_t\varphi\|_{L^1(Q)}. \] We define also $S_D = \mathcal{S}(\theta_D)$, where $\mathcal{S}=\beta^{-1}$, and $P_{wD} = P_w(S_D, P_D)$, $P_{gD} = P_g(S_D, P_D)$. The initial conditions are \begin{align} \theta(x,0) =\theta_0(x), \quad P(x,0)= p_0(x) \quad\text{in } \Omega. \label{ic} \end{align} We are going to prove the existence of weak solutions of the coupled system \eqref{eqn-w}, \eqref{eqn-g} with the boundary and initial conditions \eqref{bc-1}, \eqref{bc-2} and \eqref{ic} under the following assumptions: \begin{itemize} \item[(A1)] The porosity $\Phi$ belongs to $L^{\infty}(\Omega)$, and there exist constants, $0<\phi_m\leq \phi_M <+ \infty$, such that $0<\phi_m \leq \Phi(x) \leq \phi_M$ a.e. in $\Omega$. \item[(A2)] The permeability tensor $\mathbb{K}$ belongs to $\left(L^{\infty}(\Omega)\right)^{d\times d}$, and there exist constants $00$ such that for all $S\in [0,1]$ \[ 0 < \lambda_m \leq\lambda_w(S) + \lambda_g(S) \leq \lambda_M. \] \item[(A4)] There exist constants $p_{c,\rm min}>0$ and $M>0$ such that the capillary pressure function $S\mapsto P_c(S)$, $P_c\in C([0,1[; \mathbb{R}^+)\cap C^1(]0,1[; \mathbb{R}^+)$, for all $S\in ]0,1[$ satisfy \begin{gather} P_c'(S) \geq p_{c,\rm min}>0,\label{H4:pc_min}\\ \int_0^1 P_c(s)\, ds + \sqrt{\lambda_g(S)\lambda_w(S)}P_c'(S)\leq M.\label{H4:M} \end{gather} Moreover, there exist $S^{\#}\in ]0,1[$ and $\gamma > 0$ such that for all $S\in ]0,S^{\#}]$ \begin{gather} S^{2-\gamma}P_c'(S) \leq M, \label{H5:left-new} \\ P_c(S) - P_c(0) \leq M S P_c'(S). \label{H5:right} \end{gather} \item[(A5)] There exist $0<\tau < 1$ and $C>0$ such that for all $S_1, S_2\in [0,1]$ \[ C \Big| \int_{S_1}^{S_2} \sqrt{\lambda_g(s)\lambda_w(s)}\, ds\Big|^\tau \geq |S_1 -S_2|. \] \item[(A6)] $\rho_w > 0$, $\rho_g$ is a $C^1(\mathbb{R})$ increasing function, and there exist $\rho_m, \rho_M >0$ such that for all $p\in\mathbb{R}$ it holds \begin{align*} \rho_m \leq \rho_g(p) \leq \rho_M,\quad 0 < \rho_g'(p) \leq \rho_M. \end{align*} \item[(A7)] $F_w, F_g \in L^2(Q)$; $F_w \geq 0$ a.e.\ in $Q$. \item[(A8)] The boundary and initial data satisfy: \begin{gather*} P_D, P_c(S_D) \in W,\quad 0 \leq S_D \leq 1\text{ a.e. in }Q; \\ G_w, G_g \in L^2(\Gamma_N),\quad G_w \leq 0; \\ p_0, \theta_0 \in L^2(\Omega),\quad 0\leq \theta_0\leq \beta(1) \text{ a.e. in }\Omega. \end{gather*} \end{itemize} \begin{remark}\label{rem-0} \rm Assumptions (A1)-(A3) are standard. Assumptions (A4) and (A5) are used to prove the H\"{o}lder continuity of the functions $\mathcal{S} = \beta^{-1}$ and $(S,P)\mapsto\rho_g(S,P)S$ in the proofs of Lemma~\ref{Lemma-comp-S} and Lemma~\ref{Lemma-comp-rg}. We note that, as a consequence of incompressibility of the wetting phase, the restrictions on the capillary pressure $P_c$ in $(A4)$ are given only at $S=0$, which is less strict compared to the corresponding assumptions in \cite{AJZ2}, where both phases are compressible. The requirements on the sign of the boundary data in (A7) and (A8) are necessary only if the capillary pressure curve is unbounded at $S = S_g = 1$. In that case the restrictions $F_w \geq 0$ and $ G_w \leq 0$ do not allow extraction of the wetting phase from the domain, since otherwise we can not control the growth of the wetting phase pressure to $- \infty$. From the assumptions $P_D, P_c(S_D) \in W$ in (A8) it easily follows that the functions $P_{wD} = P_w(S_D, P_D)$ and $P_{gD} = P_g(S_D, P_D)$ are also in the space $W$. These are the conditions on boundary data that allow us to obtain uniform a priori estimates in Section~\ref{Uniform_h}. Note also that, due to \eqref{H4:pc_min} in $(A4)$, we have also $S_D \in W$. \end{remark} \begin{remark}\label{rem-1} \rm It can be seen \cite{AJ} that $\omega$ is a smooth function for which there is a constant $C$ such that \begin{align} e^{-CS} \leq \omega(S,P) \leq 1 \textrm{ in }[0,1] \times \mathbb{R}. \label{omega-bound} \end{align} It also follows from \eqref{H4:pc_min} and (A5) that $\mathcal{S} = \beta^{-1}$ is H\"older continuous with exponent $\tau$. More precisely, \begin{align} \frac{p_{c,\rm min}^\tau}{C} |S_2 - S_1| \leq |\beta(S_2)-\beta(S_1)|^\tau. \label{beta:inv:hold} \end{align} \end{remark} In order to deal with the Dirichlet boundary condition, we introduce the space \[ V=\{u\in H^1(\Omega),u|_{\Gamma_D}=0\}. \] The existence result for weak solutions of the system \eqref{eqn-w}-\eqref{eqn-g} with the boundary and initial conditions \eqref{bc-1}-\eqref{ic} is given by the following theorem. \begin{theorem}\label{TM1} Let {\rm (A1)--(A8)} hold. Denote $S =\mathcal{S}(\theta)$. Then there exists $(P,\theta)$ such that \begin{gather*} P\in L^2(0,T;V) + P_D,\; \theta\in L^2(0,T; V) + \theta_D,\quad 0\leq \theta\leq \beta(1) \text{ a.e. in } Q,\;\\ \partial_t (\Phi S) \in L^2(0,T; V'),\quad \partial_t( \Phi\rho_g(S,P) S )\in L^2(0,T; V'); \end{gather*} for all $\varphi, \psi\in L^2(0,T; V)$, \begin{equation} \begin{split} &-\rho_w\int_0^T \langle \partial_t (\Phi S),\varphi\rangle\,dt + \int_{Q} [\Lambda_w(S,P) \mathbb{K}\nabla P\cdot\nabla\varphi - A(S,P)\mathbb{K}\nabla\theta \cdot\nabla\varphi] \,dx\,dt \\ & - \int_{Q} \lambda_w(S) \rho_w^2 \mathbb{K} \mathbf{g}\cdot\nabla\varphi \,dx\,dt\\ &= \int_{Q} F_w\varphi \,dx\,dt - \int_{\Gamma_N^T} G_w \varphi \,d\sigma\,dt, \end{split} \label{eqn-weak-w} \end{equation} \begin{equation} \begin{split} &\int_0^T \langle \partial_t(\Phi \rho_g(S,P) S ),\psi\rangle\,dt + \int_{Q} [ \Lambda_g(S,P) \mathbb{K}\nabla P\cdot\nabla\psi + A(S,P)\mathbb{K}\nabla\theta \cdot\nabla\psi ] \,dx\,dt \\ & - \int_{Q} \lambda_g(S)\rho_g(S,P)^2 \mathbb{K} \mathbf{g}\cdot\nabla\psi \,dx\,dt\\ &= \int_{Q} F_g\psi \,dx\,dt - \int_{\Gamma_N^T} G_g \psi \,d\sigma\,dt. \end{split} \label{eqn-weak-g} \end{equation} Furthermore, for all $\psi\in V$ the functions \[ t\mapsto \int_\Omega \Phi S \psi\,dx ,\quad t\mapsto \int_\Omega \Phi \rho_{g}(P_{g}(S,P))S \psi\,dx \] are continuous in $[0,T]$ and the initial conditions are satisfied in the following sense: \begin{gather*} \Big( \int_\Omega \Phi S \psi\,dx\Big)(0) = \int_\Omega \Phi s_0 \psi\,dx,\\ \Big( \int_\Omega \Phi \rho_{g}(P_{g}(S,P))S \psi\,dx\Big)(0) = \int_\Omega \Phi \rho_{g}(P_{g}(s_0,p_0))s_0 \psi\,dx, \end{gather*} where $s_0=\mathcal{S}(\theta_0)$. \end{theorem} The main difficulty in proving Theorem~\ref{TM1} is the degeneracy of the equations caused by vanishing of the coefficient $A(S,P)$ at $S=0$ and $S=1$. Therefore, we will introduce a regularized problem with a strictly positive $A(S,P)$ by adding a small positive constant $\eta$ to it. At the same time we will regularize the unbounded capillary pressure function and prove Theorem~\ref{TM1} by passing to the limit as $\eta\to0$ in the regularized problem. The other difficulty in proving Theorem~\ref{TM1} is vanishing of the time derivative term in equation \eqref{eqn-weak-g} for $S=0$. This will be treated in Section~\ref{Limit_eta} with appropriate compactness theorem. \section{Regularized problem}\label{Regul} The regularized problem that we construct in this section is based on the global pressure $P$ and the non-wetting phase saturation $S$ as primary variables. A priori estimates, uniform with respect to the regularization parameter $\eta$, will be developed in Section~\ref{Uniform_h} by using the phase pressures $P_w(S,P)$ and $P_g(S,P)$ as the test functions in the variational equations of the problem. This use of the phase pressures introduces a new problem since, under hypothesis $(A4)$ on the capillary pressure, the wetting and the non-wetting phase pressure partial derivatives with respect to $S$ can be unbounded at $S = 1$ and $S = 0$ and thus for $P, S \in L^2(0,T;H^1(\Omega))$, $P_w(S,P)$ and $P_g(S,P)$ are not valid test functions. Therefore, following \cite{AJZ2}, we will correct the capillary pressure function by introducing a regularized capillary pressure derivative, a regularized capillary pressure and regularized phase pressures as follows: \begin{gather}\label{cp-der-reg} R_\eta(P_c'(S)) = \begin{cases} 2(1- \frac{S}{\eta})\frac{P_c(\eta)-P_c(0)}{\eta} + (2\frac{S}{\eta} -1)P_c'(\eta) &\text{for } S \leq \eta\\ P_c'(S) & \text{for } \eta \leq S \leq 1-\eta\\ P_c'(1-\eta) & \text{for } 1-\eta\leq S \leq 1 \end{cases}, \\ P_c^\eta(S) = P_c(0) + \int_0^{S} R_\eta(P_c'(s))\, ds,\label{pc-eta} \\ P_g^\eta(S,P) = P + P_c(0) + \int_0^{S} {f}_w(s,P) R_\eta(P_c'(s))\, ds, \label{pg-eta} \\ P_w^\eta(S,P) = P - \int_0^{S} f_g(s,P) R_\eta(P_c'(s))\, ds.\label{pw-eta} \end{gather} It is clear that $P_g^\eta(S,P) - P_w^\eta(S,P) = P_c^\eta(S)$. %\label{pgw-pc-eta} Some properties of a regularized capillary pressure are listed below, for details see \cite{AJZ2}. For any $\eta>0$, $ P_c^\eta(S)$ is a bounded, monotone, $C^1([0,1])$ function, and $ P_c^\eta(S)= P_c(S)$ for $S\in[\eta,1-\eta]$. For sufficiently small $\eta$ it holds \begin{align} \frac{d}{dS} P_c^{\eta}(S) \geq p_{c,\rm min}/2 > 0. \label{pc-min} \end{align} Also, $ | R_\eta(P_c'(S)) | \leq p_{c,\rm max}^\eta < +\infty $ for some constant $p_{c,\rm max}^\eta$ and there is a constant $M\geq 1$ such that \begin{equation} R_\eta(P_c'(S)) \leq M P_c'(S),\quad \text{for } S \in ]0,1[. \label{monoton} \end{equation} Note that in the case $S \geq \eta$, \eqref{monoton} is easy to check and $M=1$. For $S < \eta$ we need to use \eqref{H5:right} in $(A4)$ which is assumed only to obtain \eqref{monoton}. The derivatives of the regularized phase pressures are equal as in the non-regularized case and can be written as \begin{align*} \frac{\partial P_g^\eta}{\partial P} = \frac{\partial P_w^\eta}{\partial P} = \omega^\eta(S,P). \end{align*} It is easily seen that \begin{gather} \nabla P_g^\eta = \omega^\eta(S,P)\nabla P + {f}_w(S,P) R_\eta(P_c'(S))\nabla S ,\label{gradpneta}\\ \nabla P_w^\eta = \omega^\eta(S,P)\nabla P -f_g(S,P) R_\eta(P_c'(S))\nabla S\label{gradpweta}, \end{gather} so that $ P_g^\eta, P_w^\eta\in L^2(0,T; H^1(\Omega))$ for $P,S\in L^2(0,T; H^1(\Omega))$, as intended. We are going to consider the regularized version of the system \eqref{eqn-w}, \eqref{eqn-g} in which we will replace $\rho_{g}(S,P)$ by \begin{equation} \rho_{g}^\eta(S, P) := \rho_{g}( P_{g}^\eta(S, P) ) \label{density-reg} \end{equation} and the function $A(S,P)$ by $A^\eta(S,P)$, for $\eta > 0$, defined by \begin{equation} A^\eta(S,P) = \frac{\rho_w\rho_g(S,P)}{\lambda(S,P)} \lambda_w(S)\lambda_g(S) R_\eta(P_c'(S)) +\eta >0. \label{A-eta} \end{equation} Now we define the regularized system as \begin{equation} \begin{split} &-\rho_w\Phi \frac{\partial S^\eta}{\partial t} -\operatorname{div}(\Lambda_w^\eta(S^\eta,P^\eta) \mathbb{K}\nabla P^\eta) + \operatorname{div}(A^\eta(S^\eta,P^\eta)\mathbb{K}\nabla S^\eta ) \\ & +\rho_w^2\operatorname{div}(\lambda_w(S^\eta) \mathbb{K} \mathbf{g}) = F_w , \end{split}\label{eqn-reg-w} \end{equation} \begin{equation} \begin{split} &\Phi\frac{\partial }{\partial t}(\rho_g^\eta(S^\eta,P^\eta) S^\eta ) -\operatorname{div}(\Lambda_g^\eta(S^\eta,P^\eta) \mathbb{K}\nabla P^\eta) - \operatorname{div}(A^\eta(S^\eta,P^\eta)\mathbb{K}\nabla S^\eta ) \\ & +\operatorname{div}(\lambda_g(S^\eta)\rho_g^\eta(S^\eta,P^\eta)^2 \mathbb{K} \mathbf{g}) = F_g , \end{split}\label{eqn-reg-g} \end{equation} where we define \begin{align} \Lambda_{w}^\eta(S,P) = \rho_{w} \lambda_{w}(S) \omega^\eta(S,P),\quad \Lambda_{g}^\eta(S,P) = \rho_{g}(S,P) \lambda_{g}(S) \omega^\eta(S,P) \label{Lambda_j-reg} \end{align} and introduce the regularized total mobility \begin{align} \Lambda^\eta(S,P) = \Lambda_w^\eta(S,P) + \Lambda_g^\eta(S,P), \label{Lambda-f} \end{align} and the regularized function $\beta$: \begin{equation} \beta^\eta(S) = \int_0^S \sqrt{\lambda_w(s)\lambda_g(s)}R_\eta(P_c'(s))\, ds. \label{beta-eta} \end{equation} We will denote $\mathcal{S}^\eta=(\beta^\eta)^{-1}$. Now we quote some uniform estimates and limits for regularized coefficients, proved in \cite{AJZ2}, Lemma~1. \begin{lemma}\label{Lemma-reg-coef} %\label{lemma-ap1} Assume {\rm (A4)} and {\rm (A6)}. Then there exists a constant $C >0$, independent of $\eta$, such that \begin{gather} | P_g^\eta(S,P) | \leq C ( |P| + 1),\label{pneta-bound}\\ P_w^\eta(S,P) \leq P,\label{pweta-bound}\\ | \lambda_w(S) P_w^\eta(S,P) | \leq C ( |P| + 1),\label{pweta-lambdaw-bound}\\ e^{-CS} \leq \omega^\eta(S,P)\leq 1,\label{omega_eta-bound} \end{gather} and the following sequences converge uniformly in $ [0,1]\times \mathbb{R}$ as $\eta \to 0$: \begin{gather} P_g^\eta(S,P) \to P_g(S,P), \label{conv-ae-pn}\\ \omega^\eta(S,P) \to \omega(S,P), \label{conv-ae-omega}\\ \Lambda_j^\eta(S,P) \to \Lambda_j(S,P),\quad j\in\{w,g\}, \label{conv-ae-Lambda_w}\\ \beta^\eta(S) \to \beta(S)\quad \text{uniformly in } [0,1]. \label{conv-ae-beta} \end{gather} \end{lemma} \begin{remark} \label{rem-3} \rm From the assumption on the boundary data $P_D, P_c(S_D)\in W$ in (A8) it is easy to show, as mentioned in Remark~\ref{rem-0}, that $P_{wD}, P_{gD}, \beta(S_D) \in W$. We define also $P_{wD}^\eta = P_w^\eta(S_D,P_D)$ and $P_{gD}^\eta = P_g^\eta(S_D,P_D)$. Now using the estimate \eqref{monoton} we can also show that the norms $|||P_{wD}^\eta|||$, $|||P_{gD}^\eta|||$ and $|||\beta^\eta(S_{D})|||$ are uniformly bounded with respect to the parameter $\eta$. For example, \begin{align*} \nabla P_{wD}^\eta = \omega^\eta (S_D, P_D) \nabla P_D - f_g(S_D, P_D)R_\eta(P_c'(S_D))\nabla S_D, \end{align*} which, by \eqref{monoton} and \eqref{omega_eta-bound}, gives the estimate \[ |\nabla P_{wD}^\eta| \leq |\nabla P_D| + M P_c'(S_D) |\nabla S_D|, \] leading to \[ \|\nabla P_{wD}^\eta\|_{L^2(0,T;H^1(\Omega))} \leq C (1 + \|\nabla P_D\|_{L^2(0,T;H^1(\Omega))} + \|P_c(S_D)\|_{L^2(0,T;H^1(\Omega))} ). \] Also, due to uniform convergence \eqref{conv-ae-beta} we have \begin{equation} \beta^\eta(S_D) \rightharpoonup \theta_D\text{ weakly in }L^2(0,T;H^1(\Omega))\text{ as }\eta \to 0 \label{beta_eta_S_D_conv}. \end{equation} \end{remark} The variational formulation of the regularized problem and the existence of its weak solution is presented in the following theorem. \begin{theorem}\label{TM2} Assume {\rm (A1)--(A4), (A6)--(A8)} hold and $p_0, s_0 \in H^1(\Omega)$. For all $\eta >0$ sufficiently small there exists $(P^{\eta},S^\eta)$ satisfying \begin{gather*} P^\eta \in L^2(0,T; V) + P_D,\; S^\eta\in L^2(0,T; V) + S_D,\quad 0\leq S^\eta\leq 1 \text{ a.e. in } Q, \\ \partial_t (\Phi S^\eta), \partial_t(\Phi \rho_g^\eta(S^\eta,P^{\eta}) S^\eta ) \in L^2(0,T; V'); \end{gather*} for all $\varphi, \psi\in L^2(0,T; V)$, \begin{equation} \begin{split} &-\rho_w\int_0^T \langle \partial_t (\Phi S^\eta),\varphi\rangle\,dt\\ &+ \int_{Q} [\Lambda_w^\eta(S^\eta,P^{\eta}) \mathbb{K}\nabla P^{\eta}\cdot\nabla\varphi - A^\eta(S^\eta,P^{\eta})\mathbb{K}\nabla S^\eta \cdot\nabla\varphi] \,dx\,dt \\ &- \int_{Q} \lambda_w(S^\eta)\rho_w^2 \mathbb{K} \mathbf{g}\cdot\nabla\varphi \,dx\,dt\\ &= \int_{Q} F_w\varphi \,dx\,dt - \int_{\Gamma_N^T} G_w \varphi \,d\sigma\,dt, \end{split}\label{eqn-reg-weak-w} \end{equation} \begin{equation} \begin{split} &\int_0^T \langle \partial_t( \Phi\rho_g^\eta(S^\eta,P^{\eta}) S^\eta ),\psi\rangle\,dt \\ & + \int_{Q} [ \Lambda_g^\eta(S^\eta,P^{\eta}) \mathbb{K}\nabla P^{\eta}\cdot\nabla\psi + A^\eta(S^\eta,P^{\eta})\mathbb{K}\nabla S^\eta \cdot\nabla\psi ] \,dx\,dt \\ & - \int_{Q} \lambda_g(S^\eta)\rho_g^\eta(S^\eta,P^{\eta})^2 \mathbb{K} \mathbf{g}\cdot\nabla\psi \,dx\,dt\\ &= \int_{Q} F_g\psi \,dx\,dt - \int_{\Gamma_N^T} G_g \psi \,d\sigma\,dt. \end{split}\label{eqn-reg-weak-g} \end{equation} Furthermore, $S^\eta, \rho_g^\eta(S^\eta,P^{\eta}) S^\eta \in C([0,T]; L^2(\Omega))$ and \begin{equation} S^\eta(\cdot, 0) = s_0,\; \rho_g^\eta(S^\eta,P^{\eta}) S^\eta (\cdot, 0) = \rho_g^\eta(s_0,p_0) s_0 \textrm{ a.e. in }\Omega.\label{reg-init} \end{equation} \end{theorem} Theorem~\ref{TM2} will be proved by discretization of the time derivatives, with a small parameter $h>0$, and by passing to the limit as $h\to0$. \section{Time discretization}\label{Discrete} In this Section we deal with the regularized problem \eqref{eqn-reg-weak-w}-\eqref{reg-init} for a fixed $\eta>0$ and for simplicity we skip the dependence of the saturation and the global pressure on the small parameter $\eta$ in writing. In order to discretize the regularized system \eqref{eqn-reg-weak-w}-\eqref{eqn-reg-weak-g} we approximate the time derivative by a backward difference. Namely, for each positive integer $N$ we divide $[0,T]$ into $N$ subintervals, each of length $h = T/N$. We denote $t_n = n h$ and $J_n = ]t_{n-1},t_n]$ for $1\leq n\leq N$, and for any $h>0$ we denote the time difference operator by \[ \partial^h v(t) = \frac{v(t+h)-v(t)}{h}. \] For any Hilbert space $\mathcal{H}$, let \[ l_h(\mathcal{H}) =\{ v\in L^\infty(0,T;\mathcal{H})\colon v \text{ is constant in time on each subinterval } J_n \subset [0,T]\}. \] For any $v^h \in l_h(\mathcal{H})$ we denote $v^n = (v^h)^n = v^h|_{J_n}$ and assign to $v^h$ a piecewise linear in time function \begin{equation} \tilde{v}^h = \sum_{n=1}^{N} \Big( \frac{t_n-t}{h} v^{n-1} + \frac{t-t_{n-1}}{h} v^{n}\Big) \chi_{J_n}(t), \quad \tilde{v}^h(0) = v^h(0) = v^0 \label{tilda-operator} \end{equation} which satisfies \begin{equation} \partial_t \tilde{v}^h(t) = \partial^{-h} v^h(t), \quad \text{for } t\neq nh, n=0,1,\ldots, N. \end{equation} Finally, for any function $f\in L^1(0,T; \mathcal{H})$ we define $f^h \in l_h(\mathcal{H})$ by \[ f^h(t) = \frac{1}{h} \int_{J_n} f(\tau) d\tau, \quad t\in J_n. \] The discrete problem is defined as follows: find $P^{h} \in l_h(V) + P_D^h$ and $S^h\in l_h(V) + S_D^h$ such that for all $\varphi\in l_h(V)$, \begin{equation} \begin{split} & -\rho_w\int_{Q} \Phi \partial^{-h} S^h \varphi \,dx\,dt\\ & + \int_{Q} [\Lambda_w^\eta(S^h,P^{h}) \mathbb{K}\nabla P^{h}\cdot\nabla\varphi - A^\eta(S^h,P^{h})\mathbb{K}\nabla S^h \cdot\nabla\varphi] \,dx\,dt \\ &- \int_{Q} \lambda_w(S^h)\rho_w^2\mathbb{K} \mathbf{g}\cdot\nabla\varphi \,dx\,dt\\ & = \int_{Q} F_w^h\varphi \,dx\,dt - \int_{\Gamma_N^T} G_w^h \varphi \,d\sigma\,dt, \end{split} \label{discrete-w} \end{equation} and for all $\psi\in l_h(V)$, \begin{equation} \begin{split} &\int_{Q} \Phi \partial^{-h}( \rho_g^\eta(S^h,P^{h}) S^h )\psi \,dx\,dt \\ & + \int_{Q} [\Lambda_g^\eta(S^h,P^{h}) \mathbb{K}\nabla P^{h}\cdot\nabla\psi + A^\eta(S^h,P^{h})\mathbb{K}\nabla S^h \cdot\nabla\psi ] \,dx\,dt \\ & - \int_{Q} \lambda_g(S^h)\rho_g^\eta(S^h,P^{h})^2 \mathbb{K} \mathbf{g}\cdot\nabla\psi \,dx\,dt \\ & = \int_{Q} F_g^h\psi \,dx\,dt - \int_{\Gamma_N^T} G_g^h \psi \,d\sigma\,dt, \end{split} \label{discrete-g} \end{equation} and $S^h=s_0$, $P^{h} = p_0$ for $t\leq 0$. The existence for the discrete system \eqref{discrete-w}, \eqref{discrete-g}, \eqref{reg-init} is given by the following Proposition. \begin{proposition}\label{Prop1} Assume {\rm (A1)--(A8)}. Then there exists a solution $P^{h} \in l_h(V) + P_D^h,$ $S^h\in l_h(V) + S_D^h$ of \eqref{discrete-w}, \eqref{discrete-g}; moreover, $0\leq S^h \leq 1$ a.e. in $Q$. \end{proposition} The proof of Proposition~\ref{Prop1} is analogous to that of \cite[Proposition~1]{AJZ2} and therefore we omit it. \section{Uniform estimates}\label{Uniform_h} To the limit as $h\to0$ in \eqref{discrete-w}, \eqref{discrete-g} we need a priori estimates uniform with respect to $h$. We will establish in this section the estimates that are uniform in $h$ and also in $\eta$. \begin{proposition}\label{Prop-apriori-h} Let the assumptions of Proposition~\ref{Prop1} hold. Let $(P^h)_h$ and $(S^h)_h$ be the sequences of solutions to \eqref{discrete-w}-\eqref{discrete-g} and $r_g^k = \rho_g(P_g^\eta(S^{k},P^k))S^k$. Then the following bounds hold, uniform with respect to $h$: \begin{gather} \|P^h\|_{L^2(0,T;H^1(\Omega))} + \|S^h\|_{L^2(0,T;H^1(\Omega))} + \|\beta^\eta(S^h)\|_{L^2(0,T;H^1(\Omega))} \leq C, \label{Ph_Sh_betaetaSh-H1ogr}\\ \|\tilde{S}^h\|_{L^2(0,T;H^1(\Omega))} + \|{r}_{g}^h\|_{L^2(0,T;H^1(\Omega))} + \|\tilde{r}_{g}^h\|_{L^2(0,T;H^1(\Omega))} \leq C, \label{tildeSh_rgh_tildergh-H1ogr}\\ \|\partial_t(\Phi \tilde{S}^h)\|_{L^2(0,T;V')} + \|\partial_t(\Phi \tilde{r}_{g}^h)\|_{L^2(0,T;V')} \leq C.\label{Phi_h_alpha_time_der} \end{gather} \end{proposition} \begin{proof} First, we quote some identities that are going to be used throughout the proof. From the relations \eqref{gradpneta}, \eqref{gradpweta} and the definitions of the functions $A^\eta$ and $\beta^\eta$ we can obtain the following representations of the regularized wetting and non-wetting phase fluxes (without gravity term) \begin{gather} \Lambda_w^\eta(S,P) \mathbb{K}\nabla P - A^\eta(S,P)\mathbb{K}\nabla S = \rho_w \lambda_w(S) \mathbb{K}\nabla P_w^\eta - \eta\mathbb{K}\nabla S,\label{reg-flux-w}\\ \Lambda_g^\eta(S,P) \mathbb{K}\nabla P + A^\eta(S,P) \mathbb{K}\nabla S = \rho_g(S,P) \lambda_g(S) \mathbb{K}\nabla P_g^\eta +\eta \mathbb{K}\nabla S,\label{reg-flux-g} \end{gather} as well as the equality \begin{equation} \begin{aligned} & \rho_w\lambda_w(S) \mathbb{K}\nabla P_w^\eta\cdot \nabla P_w^\eta + \rho_g(S,P)\lambda_g(S) \mathbb{K}\nabla P_g^\eta\cdot \nabla P_g^\eta \\ & = \Lambda^\eta(S,P)\omega^\eta(S,P) \mathbb{K}\nabla P\cdot \nabla P + \frac{\rho_w\rho_g(S,P)}{\lambda(S,P)} \mathbb{K}\nabla \beta^\eta(S)\cdot \nabla \beta^\eta(S) \end{aligned} \label{energy-equality}. \end{equation} In this section, for simplicity, we assume that $P_c(0)=0$. From now on, $C, C_1,\ldots$ denote generic constants that do not depend on $h$ or $\eta$. We consider the discrete problem taken at a time level $k$, that is, the variational equations, \begin{equation} \label{tmp_p1-1} \begin{aligned} & \frac{\rho_w}{h} \int_{\Omega} \Phi (S^{k-1} - S^k)\varphi\,dx\\ &+\int_{\Omega} [\Lambda_w^\eta(S^k,P^{k}) \mathbb{K}\nabla P^{k}\cdot\nabla\varphi - A^\eta(S^k,P^{k})\mathbb{K}\nabla S^k \cdot\nabla\varphi]\,dx \\ &- \int_{\Omega} \lambda_w(S^k)\rho_w^2 \mathbb{K} \mathbf{g}\cdot\nabla\varphi\,dx\\ &= \int_{\Omega} F_w^k \varphi\,dx - \int_{\Gamma_N} G_w^k \varphi\,dx \end{aligned} \end{equation} for all $\varphi\in V$, and \begin{equation} \label{tmp_p1-2} \begin{aligned} &\frac{1}{h} \int_{\Omega} \Phi ( \rho_g^\eta(S^k,P^{k}) S^k - \rho_g^\eta(S^{k-1},P^{k-1}) S^{k-1} )\psi\,dx \\ &+\int_{\Omega} [ \Lambda_g^\eta(S^k,P^{k}) \mathbb{K}\nabla P^{k}\cdot\nabla\psi + A^\eta(S^k,P^{k})\mathbb{K}\nabla S^k \cdot\nabla\psi ]\,dx \\ & - \int_{\Omega} \lambda_g(S^k)\rho_g^\eta(S^k,P^{k})^2 \mathbb{K} \mathbf{g}\cdot\nabla\psi\,dx \\ &= \int_{\Omega} F_g^k\psi\,dx - \int_{\Gamma_N} G_g^k \psi\,dx \end{aligned} \end{equation} for all $\psi\in V$. Similarly as in \cite{GS1,KS3,KS1} we use the following test functions in \eqref{tmp_p1-1}, \eqref{tmp_p1-2}: \[ \varphi = \varphi(P_w^{\eta,k}) = \frac{1}{\rho_w}(P_w^{\eta,k} - P_{wD}^{\eta,k}) \quad \text{ and } \psi = \psi(P_g^{\eta,k}) = \int_{P_{gD}^{\eta,k}}^{P_g^{\eta,k}}\frac{dp}{\rho_g(p)}, \] respectively, where $P_j^{\eta,k}=P_j^\eta(S^k,P^{k})$ and $P_{jD}^{\eta,k}=P_{jD}^\eta(S_D^k,P_D^{k})$, $j\in\{w,g\}$. Taking into account \eqref{reg-flux-w} and \eqref{reg-flux-g}, the sum of the equations \eqref{tmp_p1-1} and \eqref{tmp_p1-2} with the chosen test functions reads \begin{align} &\frac{1}{h} \int_{\Omega} \Phi \big[ (S^{k-1} - S^k) P_w^{\eta,k} + (\rho_g(P_g^{\eta,k})S^k - \rho_g(P_g^{\eta,k-1}) S^{k-1}) \int_{0}^{P_g^{\eta,k}}\frac{dp}{\rho_g(p)}\big]\,dx \nonumber \\ & + \frac{1}{\rho_w}\int_{\Omega} [\lambda_w(S^k)\rho_w\mathbb{K}\nabla P_w^{\eta,k} - \eta\mathbb{K}\nabla S^k] \cdot\nabla P_w^{\eta,k} \,dx \nonumber \\ & + \int_{\Omega}\frac{1}{\rho_g(P_g^{\eta,k})} [\lambda_g(S^k)\rho_g(S^k,P^k)\mathbb{K}\nabla P_g^{\eta,k} + \eta\mathbb{K}\nabla S^k] \cdot\nabla P_g^{\eta,k}\,dx \nonumber\\ & = \frac{1}{h} \int_{\Omega} \Phi \big[(S^{k-1} - S^k) P_{wD}^{\eta,k} + (\rho_g(P_g^{\eta,k})S^k - \rho_g(P_g^{\eta,k-1}) S^{k-1}) \int_{0}^{P_{gD}^{\eta,k}}\frac{dp}{\rho_g(p)} \big]\,dx \nonumber\\ &\quad + \frac{1}{\rho_w}\int_{\Omega} [\lambda_w(S^k)\rho_w\mathbb{K}\nabla P_w^{\eta,k} - \eta\mathbb{K}\nabla S^k] \cdot\nabla P_{wD}^{\eta,k}\,dx \nonumber\\ &\quad + \int_{\Omega}\frac{1}{\rho_g(P_{gD}^{\eta,k})} [\lambda_g(S^k)\rho_g(S^k,P^k)\mathbb{K}\nabla P_g^{\eta,k} + \eta\mathbb{K}\nabla S^k] \cdot\nabla P_{gD}^{\eta,k}\,dx \label{eqn-apr-test}\\ &\quad + \int_{\Omega} [\lambda_w(S^k)\rho_w \mathbb{K} \mathbf{g}\cdot\nabla P_w^{\eta,k} +\lambda_g(S^k)\rho_g(P_g^{\eta,k}) \mathbb{K} \mathbf{g}\cdot\nabla P_g^{\eta,k}]\,dx \nonumber\\ &\quad - \int_{\Omega} [\lambda_w(S^k)\rho_w \mathbb{K} \mathbf{g}\cdot\nabla P_{wD}^{\eta,k} +\lambda_g(S^k)\frac{\rho_g^2(P_g^{\eta,k})}{\rho_g(P_{gD}^{\eta,k})} \mathbb{K} \mathbf{g}\cdot\nabla P_{gD}^{\eta,k}]\,dx \nonumber\\ &\quad + \int_{\Omega} [\frac{1}{\rho_w} F_w^k (P_w^{\eta,k} - P_{wD}^{\eta,k}) + F_g^k \int_{P_{gD}^{\eta,k}}^{P_g^{\eta,k}}\frac{dp}{\rho_g(p)} ]\,dx \nonumber\\ &\quad - \int_{\Gamma_N} [\frac{1}{\rho_w} G_w^k (P_w^{\eta,k} - P_{wD}^{\eta,k}) + G_g^k \int_{P_{gD}^{\eta,k}}^{P_g^{\eta,k}}\frac{dp}{\rho_g(p)} ]\,dx. \nonumber \end{align} Let us denote the integral terms in the expression \eqref{eqn-apr-test} by $Z_1, Z_2,\ldots, Z_{10}$, respectively. Denote the discrete time derivative terms as $$ Z_1 = \frac{1}{h} \int_{\Omega} \Phi X_1^k\,dx,\quad Z_4 = \frac{1}{h} \int_{\Omega} \Phi X_4^k\,dx. $$ We have as in \cite{KS1, AJZ2}, using the monotonicity of the non-wetting phase mass density and the monotonicity of the capillary pressure, \begin{equation} \begin{split} X_1^k &= P_w^{\eta,k} (S^{k-1}-S^k) + \rho_g(P_g^{\eta,k}) S^k \int_{0}^{P_g^{\eta,k}}\frac{dp}{\rho_g(p)} \\ &\quad - \rho_g(P_g^{\eta,k-1}) S^{k-1}\int_{0}^{P_g^{\eta,k-1}}\frac{dp}{\rho_g(p)} + \rho_g(P_g^{\eta,k-1}) S^{k-1} \int_{P_g^{\eta,k}}^{P_g^{\eta,k-1}} \frac{dp}{\rho_g(p)}\\ &\geq {H}^\eta(S^k,P^{k}) - {H}^\eta(S^{k-1},P^{k-1}), \end{split}\label{time-discr-1} \end{equation} where \[ {H}^\eta(S,P) = \Big(\rho_g(P_g^\eta)\int_{0}^{P_g^{\eta}}\frac{dp}{\rho_g(p)} - P_g^\eta \Big) S + \int_0^S P_c^\eta(z) dz. \] The monotonicity of $\rho_g$ implies \begin{equation} {H}^\eta(S,P) \geq 0.\label{H-positivity} \end{equation} The second discrete time derivative term can be transformed as follows: \begin{align*} \sum_{k=1}^N X_4^k & = \sum_{k=1}^N \big[(S^{k-1} - S^k) P_{wD}^{\eta,k} + (\rho_g(P_g^{\eta,k})S^k - \rho_g(P_g^{\eta,k-1}) S^{k-1}) \int_{0}^{P_{gD}^{\eta,k}}\frac{dp}{\rho_g(p)} \big] \\ & = s_0 P_{wD}^{\eta,0} - S^N P_{wD}^{\eta,N} + \sum_{k=1}^N S^{k-1} (P_{wD}^{\eta,k} - P_{wD}^{\eta,k-1}) \\ &\quad - \rho_g(P_g^\eta(s_0,p_0))s_0 \int_{0}^{P_{gD}^{\eta,0}}\frac{dp}{\rho_g(p)} + \rho_g(P_g^\eta(S^N,P^N))S^N \int_{0}^{P_{gD}^{\eta,N}}\frac{dp}{\rho_g(p)} \\ &\quad - \sum_{k=1}^N \rho_g(P_g^{\eta,k-1})S^{k-1} \int_{P_{gD}^{\eta,k-1}}^{P_{gD}^{\eta,k}}\frac{dp}{\rho_g(p)}, \end{align*} and by using (A1) and (A6) \begin{equation} \label{time-discr-2} \begin{aligned} | \sum_{k=1}^N \int_{\Omega} \Phi X_4^k\,dx | & \leq \frac{2 \phi_M \rho_M}{\rho_m} \Big( \sup_{t}\int_{\Omega}|P_{wD}^{\eta}|dx + \sup_{t}\int_{\Omega}|P_{gD}^{\eta}|dx \\ &\quad + \int_Q (|\partial^{-h} P_{wD}^{\eta,h}| + |\partial^{-h} P_{gD}^{\eta,h}| ) \,dx\,dt\Big) \\ & \leq C \Big( \|P_{wD}^{\eta}\|_{L^\infty(0,T;L^1(\Omega))} + \|P_{gD}^{\eta}\|_{L^\infty(0,T;L^1(\Omega))}\\ &\quad+ \|\partial_t P_{wD}^{\eta}\|_{L^1(Q)} + \|\partial_t P_{gD}^{\eta}\|_{L^1(Q)}\Big). \end{aligned} \end{equation} By applying equality \eqref{energy-equality}, relations \eqref{gradpneta} and \eqref{gradpweta}, using (A2), (A3), (A6) and the bounds \eqref{pc-min}, \eqref{omega_eta-bound} it follows that we can find a constant $C_1$ and a constant $\eta_0$, such that for all $0<\eta\leq \eta_0$, \[ Z_2 + Z_3 \geq C_1 \int_{\Omega} (|\nabla P^{k}|^2 + |\nabla \beta^\eta(S^k)|^2 + \eta |\nabla S^k |^2)\,dx . \] Using the relations \eqref{gradpneta}, \eqref{gradpweta}, the definition of $\beta^\eta$ given by \eqref{beta-eta}, and (A2), (A3) and (A6) we obtain \begin{align*} |Z_5 + Z_6| &\leq C_2 \int_{\Omega} (|\nabla P^{k}|^2 + |\nabla \beta^\eta(S^k)|^2 + \eta |\nabla S^k |^2)\,dx \, \\ &\quad + C_3 \int_{\Omega} (|\nabla P_{wD}^{\eta,k} |^2 + |\nabla P_{gD}^{\eta,k} |^2)\,dx, \end{align*} where $C_2>0$ can be chosen arbitrary small. Similarly, using the same arguments it can be seen that for any $C_4>0$, \[ |Z_7 + Z_8| \leq C_4 \int_{\Omega} (|\nabla P^{k}|^2 + |\nabla \beta^\eta(S^k)|^2)\,dx + C_5 \int_{\Omega} (1 + |\nabla P_{wD}^{\eta,k} |^2 + |\nabla P_{gD}^{\eta,k} |^2)\,dx. \] To estimate $Z_9$, we use the uniform bounds \eqref{pneta-bound} and \eqref{pweta-bound}, $(A6)$ and $F_w\geq 0$ in $(A7)$ to obtain for arbitrary $C_6>0$ \[ |Z_9| \leq C_6 \int_{\Omega} |P^{k}|^2\,dx + C_7 \int_{\Omega} (1 + |F_w^k|^2 + |F_g^k|^2 + |P_{wD}^{\eta,k}|^2 + |P_{gD}^{\eta,k}|^2)\,dx. \] Finally, in a same manner, \eqref{pneta-bound}, \eqref{pweta-bound}, (A6) and $G_w\leq0$ in (A8) imply \begin{align*} |Z_{10}| &\leq C_8 \|P^{k}\|_{H^1(\Omega)}^2 + C_9 (1 + \|P_{wD}^{\eta,k}\|_{H^1(\Omega)}^2 + \|P_{gD}^{\eta,k}\|_{H^1(\Omega)}^2 \\ &\quad + \|G_w^{k}\|_{L^2(\Gamma_N)}^2 + \|G_g^{k}\|_{L^2(\Gamma_N)}^2), \end{align*} for any $C_8>0$. Collecting the estimates for $Z_j$, $j=1,\ldots,10$ we get for arbitrary small $C_2>0$, $C_3>0$, \begin{align*} & \frac{1}{h} \int_{\Omega} \Phi (H^\eta(S^k,P^{k}) - H^\eta(S^{k-1},P^{k-1}))\,dx \\ & + C_1 \int_{\Omega} (|\nabla P^{k}|^2 + |\nabla \beta^\eta(S^k)|^2 + \eta |\nabla S^k |^2)\,dx \\ & \leq \frac{1}{h} \int_{\Omega} \Phi X_4^k\,dx + C_2 (\int_{\Omega} |\nabla P^{k}|^2\,dx + \int_{\Omega}|\nabla \beta^\eta(S^k)|^2\,dx + \eta \int_{\Omega}|\nabla S^k |^2\,dx) \\ &\quad + C_3 \int_{\Omega} |P^{k}|^2\,dx + C_4 (1 + \|P_{wD}^{\eta,k}\|_{L^2(0,T;H^1(\Omega))}^2 + \|P_{gD}^{\eta,k}\|_{L^2(0,T;H^1(\Omega))}^2 \\ &\quad + \|F_w^k\|_{L^2(\Omega)}^2 + \|F_g^k\|_{L^2(\Omega)}^2 + \|G_w^{k}\|_{L^2(\Gamma_N)}^2 + \|G_g^{k}\|_{L^2(\Gamma_N)}^2 ). \end{align*} We multiply this inequality by $h$, sum it for $k=1,\ldots,N$, take into account \eqref{time-discr-2} and use Poincar\'{e} inequality for $P$ to find \begin{equation} \label{unif-est} \begin{aligned} & \int_{\Omega} \Phi H^\eta(S^h,P^{h})(T)\,dx + C_1 \int_{Q} (|\nabla P^{h}|^2 + |\nabla \beta^\eta(S^h)|^2 + \eta |\nabla S^h |^2) \,dx\,dt \\ & \leq C_2 (1 + \|F_w\|_{L^2(Q)}^2 + \|F_g\|_{L^2(Q)}^2 +\|P_D\|_{L^2(0,T;H^1(\Omega))}^2 + |||P_{wD}^{\eta}|||^2 \\ &\quad + |||P_{gD}^{\eta}|||^2 + \|G_w\|_{L^2(\Gamma_N^T)}^2 + \|G_g\|_{L^2(\Gamma_N^T)}^2) + \int_{\Omega} \Phi H^\eta(s^{0},p^{0})\,dx. \end{aligned} \end{equation} The last term in \eqref{unif-est} is uniformly bounded with respect to $\eta$ which can be easily seen from the estimate \eqref{pneta-bound}. The other terms on the right-hand side of \eqref{unif-est} are bounded, uniformly in $\eta$, due to $(A7)$, $(A8)$, as seen in Remark~\ref{rem-3}. To obtain \eqref{Ph_Sh_betaetaSh-H1ogr}, we employ the Poincar\'{e} inequality and the fact that $P_D, S_D, \beta^\eta(S_D) \in L^2(0,T;H^1(\Omega))$ (see Remark~\ref{rem-3}). To prove \eqref{tildeSh_rgh_tildergh-H1ogr}, we first note that the functions ${S}^h$, $\tilde{S}^h$, ${r}_{g}^h$ and $\tilde{r}_{g}^h$ are uniformly bounded in $L^\infty(Q)$. Next, using \eqref{gradpneta} and \eqref{gradpweta} we easily obtain \[ | \nabla {r}_g^h | \leq C_\eta (| \nabla P^h| + | \nabla S^h|),\quad | \nabla \tilde{S}^h | \leq C | \nabla S^h|,\quad | \nabla \tilde{r}_g^h | \leq C_\eta (| \nabla P^h| + | \nabla S^h|). \] These estimates with the uniform bound \eqref{Ph_Sh_betaetaSh-H1ogr} yield the estimate \eqref{tildeSh_rgh_tildergh-H1ogr}. The estimates on the time derivatives of $\Phi\tilde{S}^h$ and $\Phi \tilde{r}_{g}^h$ are obtained in a standard way from the variational equations \eqref{discrete-w} and \eqref{discrete-g}, using the estimate \eqref{Ph_Sh_betaetaSh-H1ogr}, the boundedness of the coefficients independently of $h$ and $\eta$ and the density of $\cup_{h>0} l_h(V)$ in $L^2(0,T; V)$. The proof of Proposition~\ref{Prop-apriori-h} is completed. \end{proof} \section{Proof of Theorem~\ref{TM2}}\label{Limit_h} At this point we present an auxiliary result that is used when passing to the limit in this subsection and in Subsection \ref{Limit_eta}. \begin{lemma}\label{Lemma_ident} Let $\eta>0$ be fixed and let $(S^\varepsilon)$, $(P^\varepsilon)$ be sequences satisfying as $\varepsilon \to 0$: \begin{itemize} \item[(i)] $S^\varepsilon \to S$ a.e. in $Q$; $0 \leq S^\varepsilon \leq 1$ a.e. in $Q$; \item[(ii)] $ P^\varepsilon \rightharpoonup P$ in $L^2(Q)$; \item[(iii)] $\rho_g(P_g^\eta(S^\varepsilon, P^\varepsilon))S^\varepsilon \to r_g$ a.e. in $Q$. \end{itemize} Then $r_g = \rho_g(P_g^\eta(S, P))S$. The same is true if $P_g^\eta(S,P)$ is replaced by $P_g(S,P)$. \end{lemma} \begin{proof} Let us denote $Q^+ = \{(x,t)\in Q: S(x,t)>0 \}$ and $Q^0 = \{(x,t)\in Q: S(x,t)=0 \}$. Consider the case $S>0$. From $i)$ and $iii)$ we conclude that, as $\varepsilon \to 0$, \begin{gather*} \rho_g(P_g^\eta(S^\varepsilon, P^\varepsilon)) \to \frac{r_g}{S}\quad \text{a.e. in } Q^+,\\ P_g^\eta(S^\varepsilon, P^\varepsilon) \to (\rho_g)^{-1}(\frac{r_g}{S})\text{ a.e. in } Q^+, \end{gather*} due to the smoothness and monotonicity of $\rho_g$. Using the boundedness of $(\lambda_w P_c')(S)$ in neighborhood of $S=1$ (which follows from (A4)) and (A3), (A6), for any $S_1, S_2$ we obtain \begin{align*} |P_g^\eta(S_1, P) - P_g^\eta(S_2, P)| \leq C(|P_c^\eta(S_1) - P_c^\eta(S_2)| + |S_1 - S_2|) \end{align*} and therefore, \begin{align*} P_g^\eta(S, P^\varepsilon) \to (\rho_g)^{-1}(\frac{r_g}{S})\quad \text{a.e. in } Q^+. \end{align*} Since $P\mapsto P_g^\eta(S,P)$ is invertible (see \eqref{omega_eta-bound}), we have $P^\varepsilon \to X$ a.e. in $Q^+$, for some $X$. From (ii) we have $X=P$ so $r_g=\rho_g(P_g^\eta(S,P))S$ a.e. in $Q^+$. In the case $S=0$, the boundedness of $\rho_g$ and $i)$ imply \[ \rho_g(P_g^\eta(S^\varepsilon, P^\varepsilon)) S^\varepsilon \to r_g = 0 = \rho_g(P_g^\eta(S,P)) S \quad \text{a.e. in } Q^0. \] The same argument holds if $P_g^\eta$ is replaced by $P_g$. \end{proof} Next we obtain the convergence results holding as $h\to 0$. \begin{proposition}\label{Prop-limit-h} If the assumptions of Theorem~\ref{TM2} are satisfied then the following convergence results hold true as $h\to0$, up to a subsequence: \begin{gather} \| {S}^h -\tilde{S}^h\|_{L^2(Q)} + \| {r}_{g}^h -\tilde{r}_{g}^h\|_{L^2(Q)}\to 0, \label{conv_h_Sh_Sh_tilde_rgh_rgh_tilde}\\ {S}^h \to S \in L^2(0,T;V) + S_D \text{ weakly in } L^2(0,T; H^1(\Omega)), \text{ strongly in } L^2(Q), \label{conv_h_Sh}\\ {r}_{g}^h \to r_g = \rho_{g}(P_{g}^\eta(S,P)) S \text{ strongly in } L^2(Q),\label{conv_h_rgh}\\ \beta^\eta(S^h) \rightharpoonup \beta^\eta(S) \in L^2(0,T;V) + \beta^\eta(S_D) \text{ weakly in } L^2(0,T; H^1(\Omega)) \text{ and a.e. in } Q,\label{conv_h_beta}\\ P^{h} \rightharpoonup P \in L^2(0,T;V) + P_D \text{ weakly in } L^2(0,T; H^1(\Omega)). \label{conv_h_Ph} \end{gather} Moreover, $0\leq S\leq 1$ a.e. in $Q$ and \begin{gather} \partial_t (\Phi \tilde{S}^h) \rightharpoonup \partial_t (\Phi S) \label{conv_h_time_S} \quad \text{ weakly in } L^2(0,T; V'), \\ \partial_t (\Phi \tilde{r}_{g}^h) \rightharpoonup \partial_t (\Phi r_g) \label{conv_h_time_rg} \quad \text{ weakly in } L^2(0,T; V'). \end{gather} \end{proposition} \begin{proof} The proof of \eqref{conv_h_Sh_Sh_tilde_rgh_rgh_tilde} is the same as the proof of an analogous claim in \cite[Proposition 4]{AJZ2}. The non-homogenous boundary conditions are eliminated by using a cut-off function $\zeta$ as in Proposition~4 in \cite{AJZ2}. The weak convergence for the sequences $P^h$, $S^h$ and $\beta^\eta(S^h)$ in $L^2(0,T; H^1(\Omega))$ follow from Proposition~\ref{Prop-apriori-h}. Since $P_D^h \rightharpoonup P_D \textrm{ in }L^2(0,T; H^1(\Omega))$, we can conclude $P\in L^2(0,T; V) + P_D$. The analogous conclusion holds for the sequences $S^h$ and $\beta^\eta(S^h)$. Next, Proposition~\ref{Prop-apriori-h} allows us to apply the modification of a classical compactness result (see \cite[Lemma 7]{AJZ2}) to a sequence $\tilde{S}^h$ to obtain the relative compactness of $\tilde{S}^h$ in $L^2(Q)$. Combining \eqref{conv_h_Sh_Sh_tilde_rgh_rgh_tilde} and the weak convergence of $S^h$ to $S$ completes the proof of \eqref{conv_h_Sh}. The fact that $0\leq S^h\leq 1$ a.e. in $Q$ implies that $0\leq S\leq 1$ a.e. in $Q$. The limit of $\beta^\eta(S^h)$ is identified using a.e. convergence of $S^h$. Using the same arguments as for $\tilde{S}^h$, we obtain that $\tilde{r}_{g}^h \to r_{g}$ strongly in $L^2(Q)$ and also ${r}_{g}^h \to r_{g}$ strongly in $L^2(Q)$. Lemma~\ref{Lemma_ident} gives that $r_g = \rho_{g}(P_{g}^\eta(S,P)) S$. The existence of the weak limits \eqref{conv_h_time_S} and \eqref{conv_h_time_rg} is a consequence of the estimate \eqref{Phi_h_alpha_time_der} and the limits are identified in a standard way. This completes the proof of Proposition~\ref{Prop-limit-h}. \end{proof} Our final step in establishing the existence of weak solutions of the regularized system \eqref{eqn-reg-w}-\eqref{eqn-reg-g} is to pass to the limit as $h\to0$ in the discrete system \eqref{discrete-w}-\eqref{discrete-g}. Since we have established the pointwise convergence of the global pressure only on the set where the limit of saturation is strictly positive, we can pass to the limit as $h\to0$ in the nonlinear functions of $S$ and $P$ only if they have a special form which is given in the following lemma, whose proof is elementary. \begin{lemma}\label{Lemma_pass_limit} Let $F\in C([0,1]\times\mathbb{R})$ and let there exist functions $F_1, F_2 \in C([0,1])$ such that $F_1(S) \leq F(S, P) \leq F_2(S)$ and $F_1(0) = F_2(0)$. Then for any two sequences $(S^\varepsilon)$, $(P^\varepsilon)$, such that \begin{itemize} \item[(i)] $S^\varepsilon \to S$ a.e. in $Q$; \item[(ii)] $P^\varepsilon \to P$ a.e. in $Q^+$, \end{itemize} we have \begin{equation} F(S^\varepsilon, P^\varepsilon) \to F(S,P) \text{ a.e. in } Q.\label{limit_general} \end{equation} \end{lemma} \begin{remark}\label{rem-limit-h} \rm It is easy to verify that Lemma~\ref{Lemma_pass_limit} can be applied to all coefficients in \eqref{discrete-w}-\eqref{discrete-g} and that we have the following limits a.e. in $Q$ as $h \to 0$: \begin{gather*} \Lambda_w^\eta(S^h,P^h) \to \Lambda_w^\eta(S,P),\quad A^\eta(S^h,P^h) \to A^\eta(S,P),\quad f_w(S^h,P^h) \to f_w(S,P), \\ \Lambda_g^\eta(S^h,P^h) \to \Lambda_g^\eta(S,P),\quad \lambda_g(S^h)\rho_g^\eta(S^h,P^h)^2 \to \lambda_g(S)\rho_g^\eta(S,P)^2, \\ \rho_g^\eta(S^h,P^h)f_g(S^h,P^h) \to \rho_g^\eta(S,P)f_g(S,P). \end{gather*} \end{remark} Now we employ the convergence results in Proposition~\ref{Prop-limit-h} and Remark~\ref{rem-limit-h} to pass to the limit as $h\to0$ in the discrete system \eqref{discrete-w}-\eqref{discrete-g} and obtain \eqref{eqn-reg-weak-w}-\eqref{eqn-reg-weak-g}. Next, we conclude in a standard way that $\rho_g(P_g^\eta(S,P)) S, S\in C([0,T]; L^2(\Omega))$ and that the initial conditions are satisfied (see \cite{AJZ2}). This completes the proof of Theorem~\ref{TM2}. \section{Proof of Theorem~\ref{TM1}}\label{Limit_eta} In this section we prove the existence of weak solutions for the degenerate problem. From now on we express again the dependence of the solution of the regularized problem on parameter $\eta$. Our final step is passing to the limit as $\eta \to 0$ in the regularized problem \eqref{eqn-reg-weak-w}, \eqref{eqn-reg-weak-g}. In order to be able to apply Theorem~\ref{TM2}, we will replace the initial conditions $s_0$ and $p_0$ with the regularized initial conditions $s_0^\eta$ and $p_0^\eta$ from $H^1(\Omega)$ such that $s_0^\eta\to s_0$ and $p_0^\eta\to p_0$ in $L^2(\Omega)$ and a.e. in $\Omega$ when $\eta$ tends to zero. \begin{proposition}\label{Prop-apr-eta} For sufficiently small $\eta$, let $(P^\eta, S^\eta)_{\eta}$ be the sequence of solutions given by Theorem~\ref{TM2}. Denote $P_{g}^\eta = P_{g}^\eta(S^\eta,P^{\eta})$. The following bounds are valid, uniform with respect to $\eta$: \begin{gather} \|P^{\eta}\|_{L^2(0,T;H^1(\Omega))} \leq C, \label{P_eta-H1ogr}\\ \|\beta^\eta(S^\eta)\|_{L^2(0,T;H^1(\Omega))} \leq C, \label{beta_eta-Seta-H1ogr}\\ \|\sqrt{\eta} \nabla S^\eta\|_{L^2(Q)^d} \leq C,\label{etaSeta-L2ogr}\\ \|\partial_t (\Phi S^\eta)\|_{L^2(0,T;V')} + \|\partial_t(\Phi \rho_g(P_g^\eta) S^\eta)\|_{L^2(0,T;V')} \leq C. \label{time-eta} \end{gather} \end{proposition} \begin{proof} We pass to the limit as $h\to0$ in the estimate \eqref{unif-est} which is uniform in $\eta$, and then make use of the weak lower semicontinuity of the seminorms $f\mapsto\int_Q |\nabla f|^2\,dxdt$ to obtain \begin{align*} & \int_Q (|\nabla P^{\eta}|^2 + |\nabla \beta^\eta(S^\eta)|^2)\, \,dx\,dt + \eta \int_{Q} |\nabla S^\eta |^2 \,dx\,dt \\ & \leq C (1 + \|F_w\|_{L^2(Q)}^2 + \|F_g\|_{L^2(Q)}^2 + \|P_D\|_{L^2(0,T;H^1(\Omega))}^2 + |||P_{wD}^\eta|||^2 + |||P_{gD}^\eta|||^2 \\ &\quad + \|G_w\|_{L^2(\Gamma_N^T)}^2 + \|G_g\|_{L^2(\Gamma_N^T)}^2). \end{align*} Using Remark~\ref{rem-3} and the Poincar\'{e} inequality, \eqref{P_eta-H1ogr}, \eqref{beta_eta-Seta-H1ogr} and \eqref{etaSeta-L2ogr} follow immediately. The uniform estimates for the time derivatives of the functions $(\Phi S^\eta)_\eta$ and $(\Phi\rho_g(P_g^\eta) S^\eta)_\eta$ follow from the estimates \eqref{P_eta-H1ogr}-\eqref{etaSeta-L2ogr} by setting arbitrary $\varphi\in L^2(0,T;V)$ in the weak formulation \eqref{eqn-reg-weak-w}-\eqref{eqn-reg-weak-g}. This completes the proof of Proposition~\ref{Prop-apr-eta}. \end{proof} The compactness results for the families $(S^\eta)_{\eta}$ and $(\rho_g(P_g^\eta(S^\eta,P^\eta))S^\eta)_{\eta}$ will follow from the two following Lemmas, which can be considered as special cases of \cite[Lemma~5]{AJZ2} (see also \cite[Lemma~4.3]{KS1}) and therefore we will not give the proofs. \begin{lemma}\label{Lemma-comp-S} For any $c>0$ and for any $\eta_0>0$, the set \begin{align*} A^{c, \eta_0} = & \big\{ S \colon 0<\eta\leq \eta_0,\; \|\beta^\eta(S) \|_{L^2(0,T; H^1(\Omega))}\leq c,\; \| \partial_t(\Phi S)\|_{ L^2(0,T; V')}\leq c \big\} \end{align*} is relatively compact in $L^2(Q)$. \end{lemma} \begin{lemma}\label{Lemma-comp-rg} For any $c>0$ and for any $\eta_0>0$, the set \begin{align*} B^{c,\eta_0} =\big\{ & \rho_g(P_g^\eta(S,P)) S\colon 0<\eta\leq \eta_0,\; \|P\|_{L^2(0,T; H^1(\Omega))}\leq c, \\ & \|\beta^\eta(S) \|_{L^2(0,T; H^1(\Omega))}\leq c,\; \| \partial_t(\Phi \rho_g(P_g^\eta(S,P)) S)\|_{ L^2(0,T; V')}\leq c \} \end{align*} is relatively compact in $L^2(Q)$. \end{lemma} The limit behavior as $\eta\to0$ of the solutions to the regularized problem given by Theorem~\ref{TM2} is described by the following result. \begin{lemma} \label{Lemma-limit-eta} Let $\theta^\eta = \beta^\eta(S^\eta)$. The following convergence results hold, up to a subsequence, as $\eta\to0$. \begin{gather} P^{\eta}\rightharpoonup P \in L^2(0,T;V) + P_D \;\text{ weakly in } L^2(0,T; H^1(\Omega)),\label{conv-eta-P_eta} \\ \theta^\eta\rightharpoonup \theta \in L^2(0,T;V) + \theta_D \;\text{ weakly in } L^2(0,T; H^1(\Omega)) \text{ and a.e. in } Q, \label{conv-eta-beta}\\ S^\eta\to \mathcal{S}(\theta) \;\text{ strongly in } L^2(Q), \label{conv-eta-S_eta}\\ \partial_t (\Phi S^\eta) \rightharpoonup \partial_t (\Phi\mathcal{S}(\theta)) \;\text{ weakly in } L^2(0,T; V'), \label{conv-eta-time-S}\\ \partial_t(\Phi\rho_g(P_g^\eta(S^\eta,P^{\eta})) S^\eta ) \rightharpoonup \partial_t(\Phi\rho_g(P_g(\mathcal{S}(\theta),P))\mathcal{S}(\theta))\;\text{ weakly in } L^2(0,T; V') . \label{conv-eta-time-rg} \end{gather} Moreover, $ 0\leq \theta\leq \beta(1)$ a.e. in $Q$. \end{lemma} \begin{proof} Weak convergence of $(P^\eta)$, $(\theta^\eta)$, $\partial_t (\Phi S^\eta)$ and $\partial_t(\Phi\rho_g(P_g^\eta(S^\eta,P^{\eta})) S^\eta )$ follow from Proposition~\ref{Prop-apr-eta}. Strong convergence of $S^\eta$ to $S$ in $L^2(Q)$ is a consequence of Lemma~\ref{Lemma-comp-S}; we also get $S=\mathcal{S}(\theta)$. In order to identify the limit in \eqref{conv-eta-time-rg}, we conclude from Lemma~\ref{Lemma-comp-rg} and the uniform convergence \eqref{conv-ae-pn} that \begin{equation} \rho_g(P_g(S^\eta,P^{\eta})) S^\eta \to l_g\quad \text{ in } L^2(Q) \text{ and a.e. in } Q.\label{rgeta_comp} \end{equation} Now we apply Lemma~\ref{Lemma_ident} to obtain $l_g=\rho_g(S,P)S$. This completes the proof of Lemma~\ref{Lemma-limit-eta}. \end{proof} \begin{remark}\label{rem-limit-eta} \rm Using Lemma~\ref{Lemma_pass_limit} we get the following convergence results a.e. in $Q$, as $\eta \to 0$, \begin{gather*} \Lambda_w(S^\eta,P^\eta) \to \Lambda_w(S,P),\quad A(S^\eta,P^\eta) \to A(S,P),\quad f_w(S^\eta,P^\eta) \to f_w(S,P),\\ \Lambda_g(S^\eta,P^\eta) \to \Lambda_g(S,P),\quad \lambda_g(S^\eta)\rho_g^\eta(S^\eta,P^\eta)^2 \to \lambda_g(S)\rho_g(S,P)^2,\\ \rho_g^\eta(S^\eta,P^\eta)f_g(S^\eta,P^\eta) \to \rho_g(S,P)f_g(S,P). \end{gather*} \end{remark} Finally, we insert a test function $\varphi\in C^1([0,T];V)$ such that $\varphi(T)=0$ into \eqref{eqn-reg-w}. After the integration by parts in the time derivative term we obtain \begin{align*} & \rho_w\int_{Q} \Phi S^\eta \partial_t\varphi \,dx\,dt +\int_{Q} [\Lambda_w^\eta(S^\eta,P^{\eta}) \mathbb{K}\nabla P^{\eta}\cdot\nabla\varphi - A(S^\eta,P^{\eta})\mathbb{K}\nabla \theta^\eta \cdot\nabla\varphi] \,dx\,dt \\ &- \int_{Q} \lambda_w(S^\eta)\rho_w^2 \mathbb{K} \mathbf{g}\cdot\nabla\varphi \,dx\,dt - \eta\int_{Q}\mathbb{K}\nabla S^\eta \cdot\nabla\varphi \,dx\,dt \\ & = \int_{Q} F_w\varphi \,dx\,dt - \int_{\Gamma_N^T} G_w \varphi \,d\sigma\,dt - \rho_w\int_\Omega \Phi s_0^\eta \varphi(0)\,dx. \end{align*} Here we take into account the definitions \eqref{A-eta}, \eqref{mean-density} and \eqref{beta-eta} which give \[ A^\eta(S^\eta,P^{\eta})\nabla S^\eta = A(S^\eta,P^{\eta})\nabla \beta^\eta(S^\eta) + \eta\nabla S^\eta. \] We can pass to the limit as $\eta\to 0$ in the nonlinear terms using pointwise convergence in Remark~\ref{rem-limit-eta} and uniform convergence in Lemma~\ref{Lemma-reg-coef}. The penalisation term tends to zero due to \eqref{etaSeta-L2ogr}. Thus we obtain \begin{align*} & \rho_w\int_{Q} \Phi S \partial_t\varphi \,dx\,dt +\int_{Q} [\Lambda_w(S,P) \mathbb{K}\nabla P\cdot\nabla\varphi - A(S,P)\mathbb{K}\nabla \theta \cdot\nabla\varphi] \,dx\,dt \\ &- \int_{Q} \lambda_w(S)\rho_w^2 \mathbb{K} \mathbf{g}\cdot\nabla\varphi \,dx\,dt \\ & = \int_{Q} F_w\varphi \,dx\,dt - \int_{\Gamma_N^T} G_w \varphi \,d\sigma\,dt - \rho_w\int_\Omega \Phi s_0 \varphi(0)\,dx, \end{align*} where $S=\mathcal{S}(\theta)$. In the same way, taking a test function $\psi\in C^1([0,T];V)$ with $\psi(T)=0$ we get by integration by parts from \eqref{eqn-reg-weak-g} \begin{align*} &-\int_{Q} \Phi \rho_g^\eta(S^\eta,P^{\eta}) S^\eta \partial_t\psi \,dx\,dt \\ &+\int_{Q} [\Lambda_g^\eta(S^\eta,P^{\eta}) \mathbb{K}\nabla P^{\eta}\cdot\nabla\psi + A(S^\eta,P^{\eta})\mathbb{K}\nabla \theta^\eta \cdot\nabla\psi] \,dx\,dt \\ &- \int_{Q} \lambda_g(S^\eta)\rho_g^\eta(S^\eta,P^{\eta})^2 \mathbb{K} \mathbf{g}\cdot\nabla\psi \,dx\,dt + \eta\int_{Q}\mathbb{K}\nabla S^\eta \cdot\nabla\psi \,dx\,dt \\ & = \int_{Q} F_g\psi \,dx\,dt - \int_{\Gamma_N^T} G_g \psi \,d\sigma\,dt + \int_\Omega \Phi\rho_g^\eta(s_0^\eta,p_0^\eta)s_0^\eta \psi(0)\,dx, \end{align*} and after passing to the limit as $\eta\to 0$, \begin{align*} &-\int_{Q} \Phi \rho_g(S,P) S \partial_t\psi\,dt + \int_{Q} [ \Lambda_g(S,P) \mathbb{K}\nabla P\cdot\nabla\psi + A(S,P)\mathbb{K}\nabla \theta \cdot\nabla\psi ] \,dx\,dt \\ & - \int_{Q} [ \lambda_g(S)\rho_g(S,P)^2 \mathbb{K} \mathbf{g}\cdot\nabla\psi - \rho_g(S,P) f_g(S,P) F_P\psi ] \,dx\,dt\\ & = \int_{Q} F_g \psi \,dx\,dt - \int_{\Gamma_N^T} G_g \psi \,d\sigma\,dt + \int_\Omega \Phi\rho_g(s_0,p_0) s_0 \psi(0)\,dx. \end{align*} Using the fact that the functions $\Phi\rho_g(S,P) S$ and $\Phi S$ belong to $C([0,T]; V')$ and by an integration by parts, we easily conclude that the initial condition in Theorem~\ref{TM1} is satisfied and the proof of Theorem~\ref{TM1} is completed. \subsection*{Acknowledgments} This research was partially supported by the AUF ``Agence Universitaire de la Francophonie'', Projet de Coop\'eration Scientifique Inter - Universitaire, PHC COGITO 2011-2012 PROJET \# 24775XE, and the GnR MoMaS (PACEN/CNRS, ANDRA, BRGM, CEA, EDF, IRSN) France, their supports are gratefully acknowledged. \begin{thebibliography}{00} \bibitem{AB} H. W. Alt, E. Di Benedetto; \emph{Nonsteady flow of water and oil through inhomogeneous porous media}, Ann. Scu. Norm. Sup. Pisa Cl. Sci. \textbf{12} (1985), 335--392. \bibitem{AJ} B. Amaziane, M. Jurak; \emph{A new formulation of immiscible compressible two-phase flow in porous media}, C. R. Mecanique {\textbf 336} (2008), 600--605. \bibitem{AJZ1} B. Amaziane, M. Jurak, A. \v Zgalji\' c Keko; \emph{ Modeling and numerical simulations of immiscible compressible two-phase flow in porous media by the concept of global pressure}, Transp. Porous Media {\textbf 84} (2010), 133--152. \bibitem{AJZ2} B. Amaziane, M. Jurak, A. \v Zgalji\' c Keko; \emph{ An existence result for a coupled system modeling a fully equivalent global pressure formulation for immiscible compressible two-phase flow in porous media}, J. Differential Equations {\textbf 250} (2011), 1685--1718. \bibitem{AJZ3} B. Amaziane, M. Jurak, A. \v Zgalji\' c Keko; \emph{Numerical simulations of water-gas flow in heterogeneous porous media with discontinuous capillary pressures by the concept of the global pressure} J. Comput. Appl. Math. (2012), in press. \bibitem{AY-HK-ZA-96} Y. Amirat, K. Hamdache, A. Ziani; \emph{Mathematical analysis for compressible miscible displacement models in porous media}, Math. Models Methods Appl. Sci. \textbf{6} (1996), 729--747. \bibitem{AY-MM-95} Y. Amirat, M. Moussaoui; \emph{Analysis of a one-dimensional model for compressible miscible displacement in porous media}, SIAM J. Math. Anal. \textbf{26} (1995), 659--674. \bibitem{AY-SV-07} Y. Amirat, V. Shelukhin; \emph{Global weak solutions to equations of compressible miscible flows in porous media}, SIAM J. Math. Anal. 38 (6) (2007) 1825--1846. \bibitem{ANDRA} ANDRA; \emph{Dossier 2005 Argile, les Recherches de l'Andra sur le Stockage G\'eologique des D\'echets Radioactifs \`a Haute Activit\'e et \`a Vie Longue}, Collection les Rapports, Andra, Ch\^atenay-Malabry, 2005. \bibitem{ant-kaz-mon1990} S.N. Antontsev, A.V. Kazhikhov, V.N. Monakhov; \emph{Boundary Value Problems in Mechanics of Nonhomogeneous Fluids}, North-Holland, Amsterdam, 1990. \bibitem{Arb92} T. Arbogast; \emph{The existence of weak solutions to single-porosity and simple dual-porosity models of two-phase incompressible flow}, Nonlinear Anal. \textbf{19} (1992), 1009--1031. \bibitem{BA-HA-95} A. Bourgeat, A. Hidani; \emph{A result of existence for a model of two-phase flow in a porous medium made of different rock types}, Appl. Anal. \textbf{56} (1995), 381--399. \bibitem{CF-SB-SM} F. Caro, B. Saad, M. Saad; \emph{Two-component two-compressible flow in a porous medium}, Acta Appl. Math. \textbf{117} (2012), 15--46. \bibitem{CG-09} G. Chavent; \emph{A fully equivalent global pressure formulation for three-phases compressible flows}, Appl. Anal. \textbf{88} (2009), 1527--1541. \bibitem{GC-JJ} G. Chavent and J. Jaffr\'e; \emph{Mathematical Models and Finite Elements for Reservoir Simulation}, North--Holland, Amsterdam, 1986. \bibitem{Chen01} Z.~Chen; \emph{Degenerate two-phase incompressible flow I: Existence, uniqueness and regularity of a weak solution}, J. Differential Equations \textbf{171} (2001), 203--232. \bibitem{ZC-GH-YM-06} Z.~Chen, G. Huan, Y. Ma; \emph{Computational Methods for Multiphase Flows in Porous Media}, SIAM, Philadelphia, 2006. \bibitem{CC-08} C. Choquet; \emph{On a fully coupled nonlinear parabolic problem modelling miscible compressible displacement in porous media}, J. Math. Anal. Appl. \textbf{339} (2008), 1112--1133. \bibitem{CroiseEtAl} J. Crois\'e, G. Mayer, J. Talandier, J. Wendling; \emph{ Impact of water consumption and saturation-dependent corrosion rate on hydrogen generation and migration from an intermediate-level radioactive waste repository}, Transp. Porous Media \textbf{90} (2011), 59--75. \bibitem{DFZ-ER-HD-07} 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 the pressure}, J. Math. Anal. Appl. \textbf{326} (2007), 332--351. % \bibitem{FBK-07} B. K. Fadimba; \emph{On existence and uniqueness for a coupled system modeling immiscible flow through a porous medium}, J. Math. Anal. Appl. \textbf{328} (2007), 1034--1056. \bibitem{FX-94} X. Feng; \emph{Strong solutions to a nonlinear parabolic system modeling compressible miscible displacement in porous media}, Nonlinear Anal. \textbf{23} (1994), 1515--1531. \bibitem{GG-MT-96} G. Gagneux, M. Madaune-Tort; \emph{Analyse Math\'ematique de Mod\`eles Non Lin\'eaires de l'Ing\'enierie P\'etroli\`ere}, Springer-Verlag, Berlin, 1996. \bibitem{GS1} C. Galusinski, M. Saad; \emph{ On a degenerate parabolic system for compressible, immiscible, two-phase flows in porous media}, Adv. Differential Equations \textbf{9} (2004), 1235--1278. % \bibitem{GS2} C. Galusinski, M. Saad; \emph{ Water-gas flow in porous media}, Discrete Contin. Dyn. Syst. \textbf{Suppl.} (2005), 307--316. \bibitem{GS3} C. Galusinski, M. Saad; \emph{ A nonlinear degenerate system modelling water-gas flows in porous media}, Discrete Contin. Dyn. Syst. Ser. B, \textbf{9} (2008), 281--308. \bibitem{GS4} C. Galusinski, M. Saad; \emph{ Two compressible immiscible fluids in porous media}, J. Differential Equations, \textbf{244} (2008), 1741--1783. \bibitem{GS5} C. Galusinski, M. Saad; \emph{ Weak solutions for immiscible compressible multifluid flows in porous media}, C. R. Math. Acad. Sci. Paris, \textbf{347} (2009), 249--254. \bibitem{KS1} Z. Khalil, M. Saad; \emph{ Solutions to a model for compressible immiscible two phase flow in porous media}, Electron. J. Diff. Equ. \textbf{2010}, \textbf{122}, (2010), 1--33. \bibitem{KS2} Z. Khalil, M. Saad; \emph{ Degenerate two-phase compressible immiscible flow in porous media: The case where the density of each phase depends on its own pressure}, Math. Comput. Simulation \textbf{81} (2011), 2225--2233. \bibitem{KS3} Z. Khalil, M. Saad; \emph{ On a fully nonlinear degenerate parabolic system modelling immiscible gas-water displacement in porous media}, Nonlinear Analysis: Real World Applications \textbf{12} (2011), 1591--1615. \bibitem{KD-LS-84} D. Kr\"oner, S. Luckhaus; \emph{Flow of oil and water in a porous medium}, J. Differential Equations \textbf{55} (1984), 276--288. \bibitem{LS} M. Lenzinger, B. Schweizer; \emph{ Two-phase flow equations with outflow boundary conditions in the hydrophobic-hydrophilic case}, Nonlinear Analysis \textbf{73} (2010), 840--853. \bibitem{LB-SW} B. Li, W. Sun; \emph{ Global existence of weak solution for nonisothermal multicomponent flow in porous textile media}, SIAM J. Math. Anal. \textbf{42} (2010), 3076--3102. \bibitem{MA-09} A. Mikeli\'c; \emph{An existence result for the equations describing a gas--liquid two--phase flow}, C. R. M\'ecanique \textbf{337} (2009), 226--232. \bibitem{Norris} S. Norris; \emph{Summary of Gas Generation and Migration Current State-of-the-Art}, FORGE D1.2-R, 2009, Available online at http://www.bgs.ac.uk/forge/reports.html. \bibitem{OECD-NEA} OECD/NEA; \emph{Safety of Geological Disposal of High-level and Long-lived Radioactive Waste in France, An International Peer Review of the ``Dossier 2005 Argile'' Concerning Disposal in the Callovo-Oxfordian Formation}, OECD Publishing, 2006, Available online at http://www.nea.fr/html/rwm/reports/2006/nea6178-argile.pdf. \bibitem{SengerEtAl} R. Senger, J. Ewing, K. Zhang, J. Avis, P. Marschall, I. Gauss; \emph{Modeling approaches for investigating gas migration from a deep low/intermediate level waste repository (Switzerland)}, Transp. Porous Media \textbf{90} (2011), 113--133. \bibitem{Simon} J. Simon; \emph{ Compact sets in the space $L^p(O,T; B)$}, Ann. Mat. Pura Appl. \textbf{146} (1987), 65--96. \bibitem{Smai1} F. Sma\"\i; \emph{A model of multiphase flow and transport in porous media applied to gas migration in underground nuclear waste repository}, C. R. Math. Acad. Sci. Paris \textbf{347} (2009), 527--532. \bibitem{Smai2} F. Sma\"\i; \emph{Existence of solutions for a model of multiphase flow in porous media applied to gas migration in underground nuclear waste repository}, Appl. Anal. \textbf{88} (2009), 1609--1616. \bibitem{YLM-02} L.M. Yeh; \emph{On two-phase flow in fractured media}, Math. Models Methods Appl. Sci. 12 (8) (2002) 1075--1107. \bibitem{YLM-06} L.M. Yeh; \emph{H\"older continuity for two-phase flows in porous media}, Math. Methods Appl. Sci. \textbf{29} (2006), 1261--1289. \bibitem{AZK} A. \v{Z}galji\'{c} Keko; \emph{Modelling, Analysis and Numerical Simulations of Immiscible Compressible Two-Phase Fluid Flow in Heterogeneous Porous Media}, PhD thesis, University of Zagreb, 2011. Available online at http://www.zpm.fer.hr/\~\ zgaljic/ \end{thebibliography} \end{document}