\documentclass[reqno]{amsart}
\usepackage{hyperref}

\AtBeginDocument{{\noindent\small
\emph{Electronic Journal of Differential Equations},
Vol. 2015 (2015), No. 126, pp. 1--27.\newline
ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu
\newline ftp ejde.math.txstate.edu}
\thanks{\copyright 2015 Texas State University - San Marcos.}
\vspace{9mm}}

\begin{document}
\title[\hfilneg EJDE-2015/126\hfil Seawater intrusion models]
{Global existence for seawater intrusion models:
comparison between sharp interface and sharp-diffuse interface approaches}

\author[C. Choquet, J. Li, C. Rosier \hfil EJDE-2015/126\hfilneg]
{Catherine Choquet, Ji Li, Carole Rosier}

\address{Catherine Choquet \newline
Univ. La Rochelle, MIA, Avenue M. Cr\'epeau, F-17042,
La Rochelle, France. \newline
CNRS EA 3165, France}
\email{cchoquet@univ-lr.fr}

\address{Ji Li\newline
Univ. Lille Nord de France, F-59000, Lille, France.\newline
ULCO, LMPA J. Liouville, BP 699, F-62 228 Calais, France.\newline
CNRS FR 2956, France}
\email{ji.li@lmpa.univ-littoral.fr}

\address{Carole Rosier \newline
Univ. Lille Nord de France, F-59000, Lille, France.\newline
ULCO, LMPA J. Liouville, BP 699, F-62 228 Calais, France.\newline
CNRS FR 2956, France}
\email{rosier@lmpa.univ-littoral.fr}

\thanks{Submitted April 2, 2015. Published May 6, 2015.}
\subjclass[2010]{35R35, 35K20, 35J60, 76S05, 76T05}
\keywords{Seawater intrusion problem; sharp-diffuse interface; existence;
\hfill\break\indent strongly coupled system;  elliptic - degenerate parabolic equations}

\begin{abstract}
 We study seawater intrusion problems in confined and unconfined aquifers.
 We compare from a mathematical point of view the sharp interface approach
 with the sharp-diffuse interface approach. We demonstrate that,
 if the diffuse interface allows to establish a more efficient and logical
 maximum principle in the unconfined case, this advantage fails in the
 confined case.
 Problems can be formulated as strongly coupled systems of partial
 differential equations which include elliptic and parabolic equations
 (that can be degenerate), the degeneracy appearing only in the sharp interface
 case. Global in time existence results of weak solutions are established
 under realistic assumptions on the data.
\end{abstract}

\maketitle
\numberwithin{equation}{section}
\newtheorem{theorem}{Theorem}[section]
\newtheorem{lemma}[theorem]{Lemma}
\newtheorem{remark}[theorem]{Remark}
\allowdisplaybreaks



%\def\N{\mathbbN}
%\def\A{\mathcal A}
%\def\C{\mathcal C}

\section{Introduction} \label{s1}

In coastal zones, which are densely populated areas, the intensive
extraction of freshwater yields to local water table depression causing sea
intrusion problems.
In order to get an optimal exploitation of fresh water and also to control
seawater intrusion in coastal aquifers, we need to develop efficient and
accurate models for simulating the transport of salt water front in coastal aquifer.
We refer to the textbooks \cite{BO,Brezis,Bear}
for general information about seawater intrusion problems.

We distinguish two important cases: the case of free aquifer and
the one of confined aquifer.
In each case, the aquifer is bounded by two layers with  lower layer always
supposed to be impermeable. The upper surface is assumed to be
impermeable in the confined case and  permeable in the unconfined case
(the interface between the saturated and unsaturated zones is thus free).

The basis of the modeling is the mass conservation law for each
species (fresh and salt water) combined with
the classical Darcy law for porous media.
In the present work we essentially have chosen to adopt the simplicity of a
sharp interface approach.
This approach is based on the assumption that the two fluids are immiscible.
We assume that each fluid is confined to a well defined portion of the flow
domain with a smooth interface separating them called sharp interface.
No mass transfer occurs between the fresh and the salt area and capillary
pressure's type effects are neglected.
This  approximation is often reasonable (see  e.g. \cite{BO} and below).
Of course, this type of model does not describe the behavior of the real
transition zone but gives information concerning the movement of the saltwater
front.

Following \cite{CDR1}, we can mix this abrupt interface approach with a phase
field approach (here an Allen--Cahn type model  in fluid-fluid context
see  e.g. \cite{Alfaro2014,Alfaro2005,CahHil58,Dube})
for re-including the existence of a diffuse interface  between fresh and salt
water where mass exchanges occur.
We thus combine the advantage of respecting the physics of the problem and
that of the computational efficiency.
The same process is applied to model the transition between the saturated
and unsaturated zones in unconfined aquifers.

From a theoretical point of view, in the unconfined case, two advantages resulting
from the addition of diffuse areas compared to the sharp interface approximation
are stated in \cite{CDR2}:

$\bullet$ If diffuse interfaces are both present,
the system has a parabolic structure, so it is not necessary to introduce
viscous terms in a preliminary fixed point for treating degeneracy as in the
case of sharp interface approach.

$\bullet$ The main advantage is that we can now demonstrate a more efficient
 maximum principle and logical from the point of view of physics, which can
not be established in the case of sharp interface approximation.
(see for instance \cite{Essaid1990,NR,Tber05}).

However the latter is no longer valid in the confined case. Indeed, we need
to assume a freshwater thickness strictly positive in the interior of the
 aquifer to ensure uniform estimation in the $L^2$ space of
 the gradient of the freshwater hydraulic head. This artificial condition
is always necessary in the case of diffuse interface.
The maximum principle is then identical in both cases
(sharp interface and sharp-diffuse interface).

The outline of this article is as follows.
Section \ref{s2} is devoted to models and their derivation: we model the
evolution of the depth $h$ of the interface between freshwater and saltwater and
of the freshwater hydraulic head (in the confined case) and of depths $h$
and $h_1$, the interface between the saturated and unsaturated zone
(in the unconfined case).
The resulting models consist in a system of strongly and nonlinearly coupled PDEs
of parabolic type in the case of free aquifer
and a system of strongly and nonlinearly coupled PDEs of elliptic-parabolic
type in the case of confined aquifer.
 In section \ref{s3} all mathematical notations are stated and global in time
existence results are established in the two following cases:
the confined case with sharp-diffuse interface approach and the unconfined
case with sharp interface approach.
The section \ref{s4} is devoted to the proof of the existence results:
we apply a Schauder fixed point strategy to a regularized and truncated
system  then we establish uniform estimates  allowing us to turn back to
the original problem.

\section{ Modeling}\label{s2}
 Introducing specific index for the fresh ($f$) and salt ($s$) waters,
we write the mass conservation law for each species (fresh and salt water)
combined with the classical Darcy law for porous media.
Hydraulic heads $\Phi_i$, $i=f,s$ are defined at elevation $z$ by
\[
\Phi_i= \frac{P_i}{\rho_i g}+z,
\]
where $P_i$ denotes the pressure.
 The Darcy law relating together the effective
velocity $q_i$ of the flow and the hydraulic head $\Phi_i$ reads:
\begin{equation} \label{Darcy}
q_i = - K_i  \nabla  (\Phi_i), \quad K_i=\frac{\kappa \rho_i g}{\mu_i}.
\end{equation}
Characteristics  $\rho_i$ and $\mu_i$ are respectively the density and
the viscosity of the fluid, $\kappa$ is the permeability of the soil and
$g$ the gravitational acceleration constant.
The matrix $K_i$ is the hydraulic conductivity.
It expresses the ability of the ground to conduct water, $K_i$ is
 proportional to $\kappa$ the permeability of
the ground which only depends on the characteristics of the porous medium
and not on the fluid.

At this point, using \eqref{Darcy}, we  derive from the mass conservation
law for each species (fresh and salt water)  the following model:
\[
S_{ i} \partial_t \Phi_i+\nabla \cdot q_i=Q_i , \quad q_i=-K_i \nabla\Phi_i,
 \quad K_i=kg\rho_i/\mu_i.
\]
The coefficient of water storage  $S_i$ ($i=f,s$) characterizes the workable
water volume. It accounts for the rock and fluid compressibility.
In general, this coefficient is extremely small because of the weak
compressibility of the fluid and of the rock. In the present work, we
choose to neglect it but we emphasize that, in the case of free aquifer,
$S_{ f} \partial_t \Phi_f$ is of order of $\phi \partial_t \Phi_f$,
with $\phi$ the porosity of the medium.\\
Let us now exploit Dupuit approximation which legitimates the upscaling
of the $3D$ problem to a $2D$ model by vertical averaging.
We integrate the mass conservation law
between the interfaces depths $h$ and $h_1$ in the fresh layer
and between $h$ and the lower topography $h_2$ , in the salty zone.
The averaged mass conservation laws for the fresh and salt water thus read
\begin{gather} \label{fresh}
  S_{ f} B_f \partial_t \tilde{\Phi}_f
 = \nabla ' \cdot (B_f  \tilde{K}_f  \nabla ' \tilde{\Phi}_f)- q_f\big|_{z=h_1}
\cdot \nabla (z-h_1) + q_f\big|_{z=h} \cdot \nabla (z-h) + B_f \tilde Q_f, \\
\label{salt}
 S_{ s} B_s \partial_t \tilde{\Phi}_s
 = \nabla ' \cdot (B_s  \tilde{K}_s  \nabla '  \tilde{\Phi}_s)
+ q_s\big|_{z=h_2}\cdot \nabla (z-h_2) - q_s\big|_{z=h} \cdot \nabla (z-h)
+  B_s\tilde Q_s,
\end{gather}
where $ \nabla ' =( \partial_{x_1},\partial_{x_2}$).
The coefficients $B_f=h_1-h$ and $B_s=h-h_2 $ denote  the thickness of the fresh
and salt water zones and $\tilde \Phi_i,  \, i=f,s$, the vertically averaged
 hydraulic heads
$$
\tilde \Phi_f = \frac{1}{B_f} \int_{h}^{h_1} \Phi_f dz \quad\text{and} \quad
\tilde \Phi_s = \frac{1}{B_s} \int_{h_2}^{h} \Phi_sdz .
$$
The source terms  $\tilde{Q}_i$,  $i=f,s $ represent distributed surface
supplies of fresh and salt water into the aquifer.
 Besides sharp interface assumption implies the continuity of the pressure
at the interface between salt and fresh water, it follows that
    \begin{equation}
    (1+\alpha)\, \tilde \Phi_s= \tilde \Phi_f +\alpha h , \quad
\alpha=\frac{\rho_s}{\rho_f} -1.
     \label{reduc3}
\end{equation}
Here the parameter $\alpha$ characterizes the densities contrast.
 Equation \eqref{reduc3} allows us to avoid $\tilde \Phi_s$ in the final system.

Our aim is now to include in the model the continuity properties across
interfaces in view of expressing the four flux terms in \eqref{fresh}-\eqref{salt}.
First, since the lower layer is impermeable, there is no flux across the
boundary $z=h_2$:
\begin{equation}
   q_{s|z=h_2} \cdot \nabla (z-h_2) =0.\label{flux1}
\end{equation}
In the same way, in the case of confined aquifer, the upper layer is impermeable, thus
\begin{equation}
   q_{f|z=h_1} \cdot \nabla(z-h_1)=0.
\label{flux2}
\end{equation}
At the interface between fresh and salt water, we present the two following
 approaches:

\noindent$\bullet$ Sharp interface approach.
 With the traditional sharp interface  characterization,
there is no mass transfer across the interface  between fresh and salt water,
 i.e.  the normal component of the effective velocity $\vec{v}$ is continue
at the interface  $z=h$,
\[
 \Big(\frac{q_{f|z=h}}{\phi}-\vec{v}\Big)\cdot \vec{n}
=\Big(\frac{q_{s|z=h}}{\phi}-\vec{v}\Big)\cdot \vec{n}=0,
 %\label{interf2}
\]
 where $\vec{n}$ denotes the normal unit vector to the interface.
Thus we obtain
 \begin{equation}
   q_{f|z=h}\cdot \nabla(z-h)=q_{s|z=h}\cdot \nabla(z-h)=\phi \partial_t h
   \label{flux3}
 \end{equation}

\noindent$\bullet$  Sharp-diffuse interface between fresh and salt water.
 This approach includes now existence of miscible zone, taking the form of
diffuse interface of characteristic thickness $\delta$ between fresh and salt water.
Upscaling the 3D-dynamics of the diffuse interface assumed ruled by a phase
field model, we obtain the following continuity equation instead of \eqref{flux3}
(see \cite{CDR1} for more details about the derivation of this equation):
 \begin{equation}
   q_{f|z=h}\cdot \nabla(z-h)=q_{s|z=h}\cdot \nabla(z-h)
=\phi (\partial_t h -\delta\Delta ' h )
   \label{flux4}
 \end{equation}

The same approach for the capillary fringe  in the unconfined case yields
\begin{equation}
   q_{f|z=h_1}\cdot \nabla(z-h_1)=\phi (\partial_t h_1 -\delta \Delta ' h_1 )
   \label{flux5}
 \end{equation}
Finally,  the following assumptions are introduced for sake of simplicity
in the notation. The medium is assumed to be isotrope and the viscosity
the same for the salt and fresh water, then
 \begin{equation}
\tilde K_s=(1+  \alpha ) \tilde K_f.
\end{equation}
We re-write models with some notational simplifications.
The `primes' are suppressed in the differentiation operators in
$\mathbb{R}^2$ and source terms are denoted without `tildes'.
We also reverse the vertical axis thus changing  $h_1$ into $-h_1$, $h$
into $-h$, $h_2$ into $-h_2$, $z$ into $-z$
(bearing in mind that now $B_s=h_2 -h$, $B_f=h-h_1$).

 In the case of confined aquifer, the well adapted unknowns are the interface
depth $h$ and the freshwater hydraulic head $\Phi_f $. We set
$\alpha \tilde K_f = K$ and $\tilde \Phi_f=\alpha f$.
The final  model then reads
 \begin{gather*}
-\nabla\cdot(K(h_2 - h_1)\nabla f )+ \nabla\cdot(K (h_2 - h)\nabla h)
=B_f Q_f +B_s Q_s,\\
\phi\frac{\partial h}{\partial t}+\nabla\cdot( K (h_2 - h)\nabla f)
- \nabla\cdot( K (h_2 - h)\nabla h )-\beta \delta_h \nabla\cdot(\phi\nabla{h})
=-B_s  Q_s.
\end{gather*}
The coefficient $\beta$ is equal to $0$ in the case of sharp interface and
to $1$ in the case of sharp-diffuse interface.

 In the case of an unconfined aquifer, the unknowns are the interfaces
 depths $h$ and $h_1$.
Since quantities $h$ and $h_1$ are only meaningful inside the aquifer,
we introduce in the final model $h ^+  = \sup (0,h)$ and
$h_1^+  = \sup (0,h_1)$.
Neglecting the storage coefficient $S_f$ and introducing the characteristic
function $\mathcal{X}_0$ on the interval $(0, + \infty)$,
the sharp-diffuse interface model reads
\begin{gather*}
\begin{aligned}
&\phi\mathcal{X}_0(h_1)\partial_t h_1-\nabla\cdot( \tilde{K}_f
 \mathcal{X}_0(h_1)((h-h_1)+(h_2-h))\nabla h_1)\\
&-\beta \nabla\cdot(\delta \phi  \tilde{K}_f \mathcal{X}_0(h_1) \nabla h_1)
 -\nabla\cdot ( \tilde{K}_f \alpha(h_2-h) \mathcal{X}_0(h) \nabla h)
=- B_f Q_f - B_s Q_s,
\end{aligned} \\
\begin{aligned}
&\phi\mathcal{X}_0(h)\partial_t h-\nabla\cdot(\alpha  \tilde{K}_f(h_2-h)
\mathcal{X}_0(h_1) \nabla h)-\beta\nabla\cdot(\delta \phi
\mathcal{X}_0(h)\nabla h)\\
&-\nabla\cdot( \tilde{K}_f \mathcal{X}_0(h_1) (h_2-h)\nabla h_1)=-B_s Q_s.
\end{aligned}
\end{gather*}
Again the coefficient $\beta$ is equal to $0$ in the case of sharp interfaces
 and to $1$ in the case of sharp-diffuse interfaces.

In the previous two systems, the first equation models the conservation
of total mass of water, while the second is modeling the mass conservation
of fresh water. This is a 2D model, the third dimension being preserved by
the  upscaling process  via the depth information $h$ and $h_1$.

\section{Mathematical setting and main results}
\label{s3}

We consider a bounded and open domain $\Omega$ of $\mathbb{R}^2$ describing
the projection of the aquifer on the horizontal plane.
The boundary of $\Omega$, assumed  $\mathcal{C}^1$, is denoted by  $\Gamma$.
The time interval of interest is $(0,T)$, $T$ being any nonnegative real number,
and we set $\Omega _T = (0,T) \times \Omega$.


\subsection{Some auxiliary results}

For any $n\in N ^*$ and  any
$p\in (1,+\infty )$, let $W^{n,p}(\Omega )$ be the usual Sobolev space,
 with the norm $\| \phi \| _{W^{n,p}(\Omega )}=
\sum_{\alpha \in  \mathbb{N}^2, \alpha \le n}
\| \partial ^{\alpha }\phi \| _{L^p(\Omega )} $.
For the sake of brevity we shall write $H^1(\Omega )=W^{1,2}(\Omega )$ and
\[
V=H_0^1(\Omega ),\quad  E=H_0^1(\Omega) \cap L^\infty(\Omega),\quad
 H=L^2(\Omega ).
\]
The embeddings  $V\subset H=H'\subset V'$ are dense and compact.
For any $T>0$, let $W(0,T)$ denote the space
$$
W(0,T):=\bigl\{ \omega \in L^2(0,T;V ),\ \partial_t \omega \in L^2(0,T;V ')\bigr\}
$$
endowed with the Hilbertian norm
$ \| \omega \| _{W(0,T)}=
\bigl( \| \omega \| ^2_{L^2(0,T;V )}+
\| \partial_t \omega \| ^2_{ L^2(0,T;V ')} \bigr) ^{1/2}$.
The following embeddings are continuous
\cite[prop. 2.1 and thm 3.1, chapter 1]{LM}
\[
W(0,T)\subset \mathcal{C}([0,T];[V ,V ']_{1/2})=\mathcal{C}([0,T];H )
\]
while the embedding
\begin{equation} \label{eq:10}
W(0,T)\subset L^2(0,T;H  )
\end{equation}
is compact (Aubin's Lemma, see \cite{Simon}).
The following result by Mignot \cite{GM} is used  in the sequel.

\begin{lemma} \label{lem1}
Let $f:\mathbb{R} \to \mathbb{R}$ be a continuous and nondecreasing
function such that
$\limsup_{|\lambda | \to +\infty}
| f(\lambda )/\lambda | <+\infty$.
Let $\omega \in L^2(0,T;H)$ be such that
$\partial_t \omega \in L^2(0,T;V')$ and
$f(\omega )\in L^2(0,T;V)$.
Then
\[
\langle \partial_t\omega ,f(\omega) \rangle_{V',V}
=\frac{d}{dt}\int_{\Omega}\Bigl( \int_0^{\omega (\cdot ,y)}f(r)\, dr\Bigr)\, dy
\hbox{ in } \mathcal{D}'(0,T).
\]
Hence for all $0\le t_1<t_2\le T$,
\[
\int_{t_1}^{t_2}
<\partial_t\omega ,f(\omega )\rangle_{V',V}\, dt
=\int_{\Omega}\Bigl( \int_{\omega (t_1,y)}^{\omega (t_2,y)}f(r)\, dr\Bigr)\, dy.
\]
\end{lemma}

\subsection{Main results}

\subsubsection{Case of confined aquifer}
We focus here on models in confined case.
We aim giving existence results of physically admissible weak solutions
for these models completed by initial and boundary conditions.
We consider that the confined aquifer is bounded by two layers,
 the lower surface corresponds to $z=h_1$ and the upper surface
 $z=h_2$.
Quantity $h_2- h_1$ is the thickness of the groundwater,  we assume that
depths $h_1, h_2$ are constant,
such that $h_2>\delta_1>0$ and without lost of generality we can set $h_1 = 0$.
We introduce functions $T_s$ and $T_f$ defined by
 $$
T_s(u)=h_2-u \quad \forall  u \in (\delta_1,h_2) \quad \text{and} \quad
 T_f(u)=\begin{cases}
u &  u \in (\delta_1,h_2)\\
0 & u\leq  \delta_1
\end{cases}
$$
Functions $T_s $ and $T_f$ are extended continuously and constantly outside
$(\delta_1 ,h_2)$ for $T_s $  and for
  $  u\geq h_2$ for  $T_f $.
  $T_s(h)$ represents the thickness of the salt water zone, the previous
extension of $T_s$ for $h\leq \delta_1$ enables us to ensure a thickness
of freshwater zone always $\geq \delta_1$ in the aquifer.
 We also emphasize that the function $T_f$ only acts on the
source term $Q_f$ for avoiding the pumping when the thickness of freshwater
zone is smaller than $\delta_1$.
Then we consider the following set of equations in $\Omega_T$:
\begin{gather}
 \phi \partial_t h -\nabla \cdot \Big(KT_s(h)
\nabla h\Big)-\nabla \cdot \Big(\beta \delta\phi  \nabla h\Big)
+ \nabla \cdot \Big(KT_s(h)\nabla f\Big) = - Q_sT_s(h) ,
\label{hconfine} \\
 -\nabla \cdot \Big( h_2 K \nabla f\Big)
    + \nabla \cdot \Big(KT_s(h) \nabla h\Big)
= Q_f T_f(h)  + Q_s T_s(h).
\label{fconfine}
\end{gather}
This system is complemented with the  boundary and initial conditions:
\begin{gather}
 h=h_D, \quad f =f_{D}\quad \text{in }\Gamma \times (0,T) ,
\label{boundary}\\
 h(0,x)=h_0(x), \quad   \text{in }\Omega ,
\label{initialc}
\end{gather}
with the compatibility conditions
  \[
h_{0}(x)=h_D(0,x), \quad x \in \Gamma .
\]

Let us now detail the mathematical assumptions.
We begin with the characteristics of the porous structure.
We assume the existence of two positive real numbers $K_{-}$ and $K_{+}$ such
that  the hydraulic conductivity tensor is a bounded elliptic and uniformly
positive definite tensor:
\[
0<K_{-}|\xi|^2\leq \sum_{i,j=1,2}K_{i,j}(x)\xi_i\xi_j\leq K_{+}|\xi|^2
<\infty\quad x\in\Omega,\; \xi\in\mathbb{R}^2, \; \xi\neq 0.
\]
We assume that porosity is constant in the aquifer.
Indeed, in the field envisaged here, the effects due to variations in
$\phi$ are negligible compared with those due to  density contrasts.
From a mathematical point of view, these assumptions do not change the
complexity of the analysis but rather avoid  complicated computations.

The source terms $Q_f$ and $Q_s$ are given  functions in $L^2(0,T;H)$
such that $Q_s\leq 0$.
Notice for instance that pumping of freshwater corresponds to
assumption $Q_f \le 0$ a.e. in $\Omega \times (0,T)$.
Functions $h_D$ and $f_{D}$ belong to the space
$\big(L^2(0,T;H^1(\Omega)) \cap H^1(0,T;(H^1(\Omega))')\big)
\times L^2(0,T;H^1(\Omega))$ while function $h_0$ is in  $H^1(\Omega)$.
Finally, we assume that the boundary and initial data satisfy conditions on
the hierarchy of interfaces depths:
\[0
< \delta _1 \leq h_D\leq h_2 \quad  \text{a.e. in } \Gamma \times (0,T),
\quad  0 <  \delta _1 \leq h_0\leq h_2 \quad \text{a.e. in } \Omega.
\]
We  state and prove  the following existence result.

\begin{theorem} \label{Theo1}
        Assume  a low spatial heterogeneity  for the hydraulic conductivity tensor:
\[
K_-\leq K_+\leq \frac{3}{2}K_-.
\]
Then for any $ T > 0 $,  problem \eqref{hconfine}-\eqref{initialc}
admits a weak solution  $(h,f) $ satisfying
$$
(h-h_D,f -f_{D}) \in W(0,T)\times L^2(0,T;H_0^1(\Omega)).
$$
Furthermore the following maximum principle holds
$$
0 < \delta _1 \leq h(t,x)\leq h_2 \quad \text{for a.e. } x \in \Omega
\text{ and for any } t \in (0,T) .
$$
  \end{theorem}

Theorem \ref{Theo1} is proven in  \cite{NR} in the degenerated case
 $\beta =0$. The main difficulty is the handling of the
degeneracy since the classical Aubin's Lemma can not be applied.
Furthermore, we need to assume the thickness of freshwater zone
 $  \geq \delta _1 > 0$ inside the aquifer to
ensure an uniform estimate in $L^2$ space of the gradient of fresh water
hydraulic head $f$.

With the additional diffuse interface (corresponding to the case $\beta =1$),
the system has a parabolic structure, it is thus no longer necessary to
introduce viscous terms  in a preliminary fixed point step for avoiding degeneracy .
But we still need  to impose a freshwater thickness strictly positive inside
the aquifer to prove an uniform estimate of the gradient of $f$ since the
 presence of the diffuse interface does not allow us to
get this estimate. We can then establish the same maximum principle
for the sharp interface approximation than for that
of  the diffuse interface.

\subsubsection{Case of unconfined aquifer}

We focus now on the unconfined case.
$\tilde K_f$ is now denoted by $K$ and we set $\alpha =1$.
We assume that depth $h_2$ is constant, $h_2>0$.
We distinguish the two approaches as follows :
\smallskip

\noindent$\bullet$   $\beta =0$.
We define functions $T_s$ and $T_f$ by
\[
T_s(u)=\begin{cases}
h_2-u &\  u\in (0,h_2)\\
0& u\leq 0\,.
\end{cases}
\quad T_f(u)=u,\quad \forall u\in(\delta_1 ,h_2).
\]
Function $T_s$ is extended continuously and constantly for $u\geq h_2$
 and $T_f$ is extended continuously and constantly outside  $(\delta_1 ,h_2)$.
This condition on $T_f$ imposes a thickness of freshwater always $\geq\delta_1$
inside the aquifer.
 \smallskip

\noindent$\bullet$   $\beta =1$.
We define functions $T_s$  and $T_f$ by
  \[
T_s(u)=h_2-u, \quad   T_f(u)=u,
  \quad \text{for } u\in (0,h_2)
\]
and  $T_s$ and $T_f$ are extended continuously and constant outside  $(0,h_2)$.

Then we consider the following set of equations in $\Omega_T$,
\begin{gather}
\begin{aligned}
&\phi \partial_t h -\nabla \cdot \Big(KT_s(h) \nabla h\Big)
 -\nabla \cdot \Big(\beta\delta\phi  \nabla h\Big)
 -\nabla \cdot \Big(KT_s(h)\mathcal{X}_0(h_1)\nabla h_1\Big)\\
&=-Q_sT_s(h) ,
\end{aligned} \label{probleml} \\
\begin{aligned}
& \phi \partial_t h_1 -\nabla \cdot \Big(K\Big(T_f(h-h_1)+T_s(h)\Big)\nabla h_1\Big)
   -\nabla \cdot \Big(\beta \delta\phi \nabla h_1\Big)\\
&-\nabla \cdot \Big(KT_s(h)\mathcal{X}_0(h_1)\nabla h\Big)   \\
&=-\mathcal{X}_0(h_1) \Big(Q_fT_f(h-h_1)+Q_sT_s(h)\Big).
\end{aligned}\label{problem}
\end{gather}
Notice that we do not use $h^+=\sup(0,h)$ and $h^+_1=\sup(0,h_1)$ in
functions $T_s$ and $T_f$ because a maximum principle will ensure that
these supremums are useless.
Likewise, we have canceled the terms  $\mathcal{X}_0(h)$
(resp. $\mathcal{X}_0(h_1)$) in front of $\partial_t h$ and $\nabla h$
 (resp.  $\partial_t h_1$).  System \eqref{problem} is completed by the
following boundary and initial conditions:
\begin{gather}
 h=h_D, \quad h_1=h_{1,D}\quad \text{in }\Gamma \times (0,T) ,
\label{boundary2}\\
 h(0,x)=h_0(x), \quad h_1(0,x)=h_{1,0}(x)\quad  \text{in }\Omega ,
\label{initial}
\end{gather}
with the compatibility conditions
  \[
h_{0}(x)=h_D(0,x),\quad h_{1,0}(x)=h_{1,D}(0,x), \quad x \in \Gamma .
\]

We make the same mathematical assumptions than above for  the porosity and the
 hydraulic conductivity tensor $K$ but we do not make any assumptions on the
sign of the source terms $Q_f$ and $Q_s$.
Functions $h_D$ and $h_{1,D}$ belong to the space
$L^2(0,T;H^1(\Omega)) \cap H^1(0,T;(H^1(\Omega))')$ while functions
 $h_0$ and $h_{1,0}$ are in  $H^1(\Omega)$.
Finally, we assume that the boundary and initial data satisfy physically
realistic conditions on the hierarchy of interfaces depths:
  \[
0 \leq h_{1,D}\leq h_D\leq h_2 \quad  \text{a.e. in } \Gamma \times (0,T), \quad
 0 \leq h_{1,0}\leq h_0\leq h_2 \quad  \text{a.e. in } \Omega.
\]
Now we state and prove the following existence result.

\begin{theorem} \label{Theo2}
Assume  a spatial heterogeneity  for the hydraulic conductivity tensor:
     \[
  K_+ \leq   2\,\sqrt{\gamma} K_-, \quad       0 < \gamma < \frac{8}{9}.
\]
 Then for any $ T > 0 $,  problem \eqref{probleml}-\eqref{initial}
admits a weak solution  $(h,h_1) $ satisfying
$$
(h-h_D,h_1-h_{1,D})
\in \bigl(L^2(0,T;H_0^1(\Omega))\times L^2(0,T;H_0^1(\Omega))
\cap H^1(0,T;(H_0^1(\Omega))')^{2}$$
Furthermore the following maximum principle holds,
 \begin{itemize}
\item If   $\beta =0$,
$0\leq h_1(t,x)$ and $0\leq h(t,x)\leq h_2$  a.e. $x \in \Omega$ and for any
$t \in (0,T) $.

\item  If  $\beta =1$, $ 0\leq h_1(t,x)\leq h(t,x)\leq h_2$  a.e.
$x \in \Omega$ and for any $t \in (0,T)$.
 \end{itemize}
  \end{theorem}

Theorem \ref{Theo2} is proven in \cite{CDR2} in the non degenerated case
$\beta =1$, with  condition $K_-\leq K_+\leq \frac{3}{2}K_-$
on the spatial heterogeneity for the hydraulic conductivity.
We aim to give an existence result of weak solutions for this model when
$\beta =0$. We introduce a viscous term depending on a parameter $ \epsilon$
in the preliminary fixed point step for avoiding degeneracy.
We again suppose the thickness of freshwater zone $ \geq \delta _1 > 0$
inside the aquifer  to ensure an uniform estimate in $L^2$ of the gradient of $h_1$.
But, since $ \epsilon$ is expected to tend to zero, we only can establish
a weaker maximum principle without hierarchy between $h_1$ and $h$ .

\begin{remark} \rm
We can prove Theorem \ref{Theo1}  without any restrictions on the sign of the
source terms $Q_f$ and $Q_s$, but in this case, we have to impose assumptions
on additional leakage terms $q_{Lf}$ and $q_{Ls}$ like in  \cite{CDR2}.

 Depths $h_1$ and $h_2$ are assumed to be constant  for sake of  simplicity
 but the proof extends directly to $h_i \in L^\infty(\Omega)$, $i=1,2$.
\end{remark}

Next section is devoted to proofs of Theorem \ref{Theo1} for $ \beta =1$
and of Theorem \ref{Theo2} for $ \beta = 0$.
Let us sketch our strategy.
First step consists in using a Schauder fixed point theorem for proving
an existence result for an auxiliary regularized and truncated problem.
More precisely, in the unconfined case,  we regularize the equations
by adding a viscous term and we also regularize the step function
$\mathcal{X}_0$ with a parameter $\epsilon >0$. Furthermore we introduce
a weight based on the velocity of the fresh front in the two equations.
We show that the regularized solution satisfies the maximum principles
announced in Theorem \ref{Theo1} and in Theorem \ref{Theo2}.
We then prove that we have sufficient control on the velocity of the fresh
front to ignore the latter weight.
We finally show sufficient uniform estimates to let the regularization
$\epsilon$ tends to zero.


\section{Proofs}\label{s4}

\subsection{Proof of Theorem \ref{Theo1}} \quad

\noindent \textbf{Step 1:} Existence for the truncated system.
Let $M $ be a positive constant to be determine later.
For $x \in \mathbb{R}_+^*$, we set
\[
L_M(x)=\min\Big(1,\frac{M}{x}\Big) .
\]
Such a truncation $L_M$ was originally introduced in \cite{RR}.
It allows to  use the following point in the estimates hereafter.

 For any $(g,g_1)\in (L^2(0,T;H^1(\Omega)))^2$, setting
$$
d(g,g_1)=-T_s(g)L_M\big(\|\nabla g_1\|_{L^2(\Omega_T)^2}\big)\nabla g_1,
$$
we have
$$
\| d(g,g_1)\|_{L^2(0,T;H)}=
 \| T_s(g)L_M\big(  \|\nabla g_1  \|_{L^2(\Omega_T)^2}\big)\nabla g_1
  \|_{{L^2(\Omega_T)^2}}\leq M   h_2  .
$$
Now, we denote $L_M\big(  \|\nabla g_1  \|_{L^2(\Omega_T)^2}\big) $ by
$L_M\big(  \|\nabla g_1  \|_{L^2}\big)$.
The variational formulation of the problem under consideration involves
the two following integral equations:
\begin{gather}
\begin{aligned}
&\int_0^T \phi \langle\partial_t h,w\rangle_{V,V'}
 +\int_{\Omega_T} \delta { \phi}\nabla h\cdot\nabla w\\
&+\int_{\Omega_T}{T_s({ h})\big(}
K \nabla h\cdot\nabla w-L_M(\|\nabla { f}\|_{{L^2(\Omega_T)^2}})
K\nabla {f }\cdot\nabla w)\,dx\,dt\\
&+\int_{\Omega_T}{Q}_s T_s({h}\big)w\,dx\,dt=0\,,
\end{aligned}\label{tconh}
\\
\begin{aligned}
&\int_{\Omega_T}h_2K \nabla f\cdot\nabla w\,dx\,dt
-\int_{\Omega_T}T_s({h}) K\nabla h\cdot\nabla w\,dx\,dt\\
&-\int_{\Omega_T}({Q}_s T_s({ h})+{Q}_f T_f({ h}))w\,dx\,dt=0.
\end{aligned} \label{tconf}
\end{gather}

For the  fixed point strategy, we define the application
\begin{gather*}
\mathcal{F}: L^2(0,T;H^1(\Omega))\times L^2(0,T;H^1(\Omega))\to
 L^2(0,T;H^1(\Omega))\times L^2(0,T;H^1(\Omega))\\
({\bar h},{\bar f})\to \mathcal{F}({\bar h},{\bar f})
=(\mathcal{F}_1({\bar h},  {\bar f})=h,\mathcal{F}_2({\bar h},{\bar f})=f),
\end{gather*}
where the pair (h,f) is a solution of next variational problem:
for all $w\in V$,
\begin{gather}
\begin{aligned}
&\int_0^T \phi \langle\partial_t h,w\rangle_{V,V'}
+\int_{\Omega_T} \delta { \phi}\nabla h\cdot\nabla w
+\int_{\Omega_T} T_s({\bar h}) \Big(K \nabla h\cdot\nabla w\\
&-L_M(\|\nabla {\bar f}\|_{ L^2})K\nabla {\bar f }\cdot\nabla w\Big)\,dx\,dt
+\int_{\Omega_T}{Q}_s T_s({\bar h}\big)w\,dx\,dt=0\,,
\end{aligned}\label{confh1}
\\
\begin{aligned}
&\int_{\Omega_T}h_2K \nabla f\cdot\nabla w\,dx\,dt
-\int_{\Omega_T}T_s({\bar h}) K\nabla h\cdot\nabla w\,dx\,dt\\
&-\int_{\Omega_T}({Q}_s T_s({\bar h})+{Q}_f T_f({\bar h}))w\,dx\,dt=0.
\end{aligned} \label{conff1}
\end{gather}
Indeed we know from  classical  parabolic theory  (see {\it e.g.}  \cite{LM})
that the linear variational system \eqref{confh1}-\eqref{conff1} admits an
unique solution.
The end of the present subsection is devoted to the proof a fixed point
property for application $\mathcal{F}$.
\smallskip

\noindent Continuity of $\mathcal{F}_1$.
Let $(\bar {h^n},\bar{ f^n})$ be a sequence of functions of
$L^2(0,T;H^1(\Omega))\times L^2(0,T;H^1(\Omega))$
 and $(\bar h,\bar f)$ be a function of
$L^2(0,T;H^1(\Omega))\times L^2(0,T;H^1(\Omega))$ such that
\[
(\bar {h^n},\bar{ f^n})\to  (\bar h,\bar f)\quad \text{in }
L^2(0,T;H^1(\Omega))\times L^2(0,T;H^1(\Omega)).
\]
We set $ h_n=\mathcal{F}_1(\bar{ h^n},\bar{ f^n})$ and
 $ h=\mathcal{F}_1(\bar h,\bar f)$.
We aim showing that $ h_n \to  h$ in $L^2(0,T;H^1(\Omega))$.

For all $n\in \mathbb{N}$, $ h_n$ satisfies \eqref{confh1}.
Choosing  $w=h_n-h_D$ in the $n$-dependent counterpart of \eqref{confh1} yields
\begin{gather*}
\begin{aligned}
&\int_0^T\phi \langle\partial_t(h_n-h_D),h_n-h_D\rangle_{V',V} dt
+\int_{\Omega_T}(\delta  {\phi+K T_s({\bar h^n}))}\nabla h_n\cdot\nabla h_n\,dx\,dt \\
&=\int_{\Omega_T}\big(T_s({\bar h^{n}})\,
L_M(\|\nabla {\bar f^{n}}\|_{L^2})K\nabla {\bar f^n}
\cdot\nabla(h_n-h_D)\big) \,dx\,dt\,,
\end{aligned}\\
\begin{aligned}
&\int_{\Omega_T} -{Q}_s T_s({\bar h^n})(h_n-h_D) \,dx\,dt
-\int_0^T \langle\partial_t h_D,h_n-h_D\rangle_{V',V}dt\\
&+\int_{\Omega_T}(\delta{ \phi+K \,
T_s({\bar h^n}))}\nabla h_n\cdot\nabla h_D\,dx\,dt
\end{aligned}
\end{gather*}
Function $h_n-h_D$ belongs to $L^2(0,T;V)\cap H^1(0,T;V')$ and
then to $\mathcal{C}(0,T;L^2(\Omega))$.
Thus, thanks moreover  to Lemma \ref{lem1}, we write
\begin{equation*}
\int_0^T\phi\langle\partial_t (h_n-h_D),(h_n-h_D)\rangle_{V',V}dt
=\frac{\phi}{2}\|h_n(\cdot, T)-h_D\|_H^2-
\frac{\phi}{2}\|h_0-h_{D|t=0}\|_H^2 .
\end{equation*}
Also
\begin{equation*}
\int_{\Omega_T}\big( \delta\phi+KT_s(\bar h^n)\big) \nabla h_n
\cdot\nabla h_n \,dx\,dt\geq \delta\phi \| \nabla h_n\|_{L^2(0,T;H)}^2.
\end{equation*}
Then applying Cauchy-Schwarz and Young inequalities, for all $ \epsilon_1 >0$
we obtain
 \begin{align*}
&\Big|\int_{\Omega_T} \big(\delta\phi
+ KT_s(\bar h^n)\big)\nabla h_n\cdot\nabla h_D \,dx\,dt\Big| \\
&\leq  (\delta\phi  + K_+ h_2)\|\nabla h_n\|_{L^2(0,T;H)}\|\nabla h_D\|_{L^2(0,T;H)}
 \\
&\leq  \frac{\varepsilon_1}{2}\|\nabla h_n\|_{L^2(0,T;H)}^2+\frac{ (\delta\phi  + K_+h_2 )^2}{2\varepsilon_1}
\|\nabla h_D\|_{L^2(0,T;H)}^2,
\end{align*}
\begin{align*}
&\Big|-\int_{\Omega_T} KT_s(\bar h^n)L_M
\big(\|\nabla  {\bar { f^n}}\|_{L^2}\big)
 \nabla  {\bar { f^n}}\cdot\nabla h_n\,dx\,dt\Big|\\
&\leq  K_+\|d(\bar h^n,\bar f^n)\|_{L^2(0,T;H)}\|\nabla h_n\|_{L^2(0,T;H)}
\\
&\leq MK_+h_2 \|\nabla h_n\|_{L^2(0,T;H)}\\
&\leq \frac{ K_+^2 M^2 }{2\varepsilon_1} h_2^2
+\frac{\varepsilon_1}{2}\|\nabla h_n\|_{L^2(0,T;H)}^2.
 \end{align*}
Since it  depends on $h_D$,  the next term is simply estimated by
\begin{align*}
&\Big|\int_{\Omega_T} KT_s(\bar h^n)L_M\big(\|\nabla \bar f^n\|_{L^2}\big)
\nabla {\bar  f^n} \cdot\nabla h_D\,dx\,dt\Big|\\
&\leq K_+\|d(\bar h^n,\bar f^n)\|_{L^2(0,T;H)}\|h_D\|_{L^2(0,T;H^1)}\\
&\leq  MK_+h_2 \|h_D\|_{L^2(0,T;H^1)}.
\end{align*}
Finally we have
  \begin{align*}
&\Big| - \int_0^T\phi\langle\partial_t h_D,(h_n-h_D)\rangle_{V',V}dt\Big|\\
&\leq  \frac{\phi}{\delta}
 \| \partial_t h_D \|^2_{L^2(0,T;(H^1(\Omega))')}
+\frac{\delta\phi}{2} \|h_n\|_{L^2(0,T; {V})}^2
+\frac{\delta\, \phi}{2} \|h_D\|_{L^2(0,T;{V})}^2  ,
\end{align*}
and
\[
\Big|-\int_{\Omega_T} Q_sT_s(\bar h^n)(h_n-h_D)\,dx\,dt\Big|
\leq  \frac{ \| Q_s\|_{L^2(0,T;H)}^2}{2\phi}h_2^2
+\frac{\phi}{2}\|h_n-h_D\|_{L^2(0,T;H)}^2.
\]
Summing  up all these estimates,  after simplifications, we obtain
\begin{equation}
\begin{aligned}
& \frac{\phi}{2}\|h_n(\cdot, T)-h_D\|_H^2
 +(\frac{\delta\phi}{2}-{\varepsilon_1}) \|\nabla h_n\|_{L^2(0,T;H)}^2\\
&\leq \frac{\phi}{2}\|h_0-h_{D|t=0}\|_H^2+\frac{{\phi}}{2}
\int_0^T\|h_n-h_D\|_{H}^2\, dt +\frac{\phi\delta}{2} \|h_D\|_{L^2(0,T;{V})}^2\\
&\quad +\Big(\frac{ \,\| Q_s\|_{{L^2(\Omega_T)}}^2}{2\phi}
 +\frac{K_+^2 M^2 }{2  \varepsilon_1}\Big) h_2^2
 +\frac{(\delta\phi  + K_+h_2)^2}{2 \varepsilon_1}\|h_D\|_{L^2(0,T;{V})}^2\\
&\quad+ \frac{\phi}{{\delta}}
 \| \partial_t h_D \|^2_{L^2(0,T;(H^1(\Omega))')}
+MK_+h_2 \|h_D\|_{L^2(0,T;H^1)}.
\end{aligned} \label{confh2}
\end{equation}
We choose $\varepsilon_1$ such that
$ \delta \phi/2-\varepsilon_1 \geq \epsilon_0 >0$ for some $\epsilon_0>0$.
Relation  \eqref{confh2} with Gronwall lemma enables to conclude
that there exists real numbers
$A_M= A_M(\phi, \delta, K, h_0, h_D, h_2, Q_s, M, T) $ and
$B_M=B_M(\phi, \delta, K, h_0, h_D, h_2, Q_s, M, T)$ depending only on
the data of the problem such that
  \begin{equation}
\|h_n\|_{L^\infty(0,T;H)} \leq  A_M , \quad
\|h_n\|_{L^2(0,T;H^1)}\leq  B_M .
\label{Chap2_15}
\end{equation}
Hence  sequence $(h_n)_n$ is uniformly bounded in
$L^2(0,T;H^1(\Omega))\cap L^\infty(0,T;H)$.
Notice that the estimate in $L^\infty(0,T;H)$ is justified by the fact that
we could make the same computations replacing  $T$ by any $\tau \leq T$
in the time integration. In the sequel, we set
$$
C_M = \max(A_M,B_M).
$$
Now we  prove that $(\partial_t( h_n-h_D))_n$ is bounded in $L^2(0,T;V')$.
  \begin{align*}
&\|\partial_t (h_n-h_D)\|_{L^2(0,T;V')}\\
&= \sup_{\|w\|_{L^2(0,T;V)}\leq1}\Big|\int_0^T\langle\partial_t (h_n-h_D),
 w\rangle_{V',V}dt\Big|\\
&= \sup_{\|w\|_{L^2(0,T;V)}\leq 1}\Big|\int_0^T
 -\langle\partial_t h_D,w\rangle_{V',V}dt
 - \frac{1}{\phi} \Big(\int_{\Omega_T} \big(\delta\phi
  + KT_s(\bar h^n)\big)\nabla h_n \cdot\nabla w \,dx\,dt\\
&\quad +\int_{\Omega_T} KT_s(\bar h^n)
L_M\big(\|\nabla \bar f^n\|_{L^2}\big)\nabla
\bar f^n\cdot\nabla w \,dx\,dt
 -\int_{\Omega_T} Q_sT_s(\bar h^n)w  \,dx\,dt   \Big)\Big|.
\end{align*}
Since
\[\Big| \int_{\Omega_T} \Big(\delta\phi  + KT_s(\bar h^n)\Big)\nabla h_n.\nabla w
 \,dx\,dt\Big|\leq \big(\delta\phi+K_+h_2\big)
\|h_n\|_{L^2(0,T;H^1(\Omega))}\|w\|_{L^2(0,T;V)},
\]
and since $h_n$ is uniformly bounded in $L^2(0,T;H^1(\Omega))$, we write
\begin{equation}
\Big|\int_{\Omega_T} \Big(\delta\phi  + KT_s(\bar h^n)\Big)\nabla h_n\cdot\nabla w
 \,dx\,dt\Big|
\leq \big(\delta\phi+K_+h_2\big)C_M\|w\|_{L^2(0,T;V)}.\label{Chap2_16}
\end{equation}
Furthermore,
\begin{gather}
\Big|\int_{\Omega_T} T_s(\bar h^n)L_M\big(\|\nabla \bar
f^n\|_{L^2}\big)\nabla \bar f^n\cdot\nabla w \,dx\,dt\Big|
\leq M h_2 \|w\|_{L^2(0,T;V)}, \label{Chap2_17}  \\
\Big|\int_{\Omega_T} Q_sT_s(\bar h^n)w \,dx\,dt\Big|
\leq  \| Q_s\|_{L^2(\Omega_T)} h_2 \|w\|_{L^2(0,T;V)}.
\label{Chap2_18}
\end{gather}
Summing up \eqref{Chap2_16}--\eqref{Chap2_18}, we conclude that
\begin{equation}
\begin{aligned}
\|\partial_t h_n \|_{L^2(0,T;V')}
&\leq \frac{1}{\phi}\Big(
 \| \partial_t h_D \|^2_{L^2(0,T;(H^1(\Omega))')}  + \delta\phi C_M \\
&\quad +h_2 ( K_+C_M+M + \| Q_s\|_{L^2(\Omega_T)} ) \Big)
:= D_M.
\end{aligned}\label{Chap2_19}
\end{equation}
We have proved that $\big( h_n \big)_n$ is uniformly bounded in
$L^2(0,T;H^1(\Omega)) \cap H^1(0,T;V')$.
 Using Aubin's lemma, we  extract a subsequence, not relabeled for convenience,
$(h_n)_n$, converging strongly in $L^2(\Omega_T)$ and weakly in the space
$L^2(0,T;H^1(\Omega)) \cap H^1(0,T;V')$ to some  limit denoted by $\ell$.
Using in particular the strong convergence in $L^2(\Omega _T)$ and thus
the convergence a.e. in $\Omega_T$, we check that $\ell$ is a solution of
\eqref{confh1}.
The solution of \eqref{confh1} being unique,  we actually have $\ell=h$.

It remains to prove that $(h_n)_n$ actually tends to $h$ strongly in
$L^2(0,T;H^1(\Omega))$.
Subtracting the weak formulation \eqref{confh1} to its $n$-dependent
counterpart for the test function $w=h_n-h$, we obtain
\begin{equation}
\begin{aligned}
&\int_0^T\phi\langle\partial_t (h_n -h),h_n-h\rangle_{V',V}dt\\
&+\int_{\Omega_T} \bigl( \delta\phi + KT_s(\bar h^n) \bigr)
 \nabla (h_n-h) \cdot\nabla (h_n-h) \,dx\,dt  \\
&-\int_{\Omega_T} K \bigl( T_s(\bar h^n) -T_s(\bar h) \bigr)
 \nabla (h_n-h) \cdot \nabla h \,dx\,dt  \\
& + \int_{\Omega_T} K \bigl( T_s(\bar h^n)L_M\big(\|\nabla
 \bar f^n\|_{L^2}\big)\nabla \bar f^n - T_s(\bar h)
L_M\big(\|\nabla \bar f\|_{L^2}\big) \nabla \bar f \bigr)
\cdot \nabla (h_n-h) \,dx\,dt\\
&+\int_{\Omega_T} Q_s \bigl( T_s(\bar h^n) - T_s(\bar h)\bigr)
(h_n-h) \,dx\,dt=0.
\end{aligned} \label{strong}
\end{equation}
Using assumption $(\bar h^n, \bar f^n) \to (\bar h, \bar f)$ in
$L^2(0,T;H^1(\Omega))\times L^2(0,T;H^1(\Omega)) $ and the above results
of convergence for $h_n$, the limit as $n \to \infty$ in \eqref{strong} reduces to
$$
\lim_{n \to \infty} \Bigl(
\int_{\Omega_T}  \bigl( \delta\phi + KT_s(\bar h^n) \bigr)
\nabla (h_n-h) \cdot\nabla (h_n-h) \,dx\,dt\Bigr) = 0 .
$$
Due to the positiveness of $K$, we infer from the latter relation that
$$
\lim_{n \to \infty} \Bigl(\int_{\Omega_T}   \delta \phi | \nabla (h_n-h)
 |^2 \,dx\,dt
+\int_{\Omega_T}    K_-T_s(\bar h^n) | \nabla (h_n-h)  |^2 \,dx\,dt
\Bigr) \leq 0 .
$$
Hence $\nabla h_n \to \nabla h$ strongly in $L^2(0,T;H)$.
Continuity of $\mathcal{F}_1$ for the strong topology of
$L^2(0,T;H^1(\Omega))$ is proved.
\smallskip


\noindent Continuity of $\mathcal{F}_2$.
 Likewise, we prove the continuity of   $\mathcal{F}_2$ by setting
$f_{n}=\mathcal{F}_2(\bar h^n,\bar f^n)$ and $f=\mathcal{F}_2(\bar h,\bar f)$
 and showing  that $f_{n} \to f$ in $L^2(0,T;H^1(\Omega))$.
 The key estimates are obtained using the same type of arguments than in the
proof of the continuity of   $\mathcal{F}_1$.
 We thus do not detail the computations.
 Let us only emphasize that we can now use the estimate  \eqref{Chap2_15}
previously derived for $h^n$, thus the dependence with regard to $C_M$
in the  estimate
\begin{equation}
\|f_{n}\|_{L^2(0,T;H^1)}\leq E_M=F_M\big(\phi, \delta, K,{ f_{D}},
 h_2, Q_s,  Q_f, M, C_M, T\big) .
\label{Chap2_30}
\end{equation}
\smallskip

\noindent Conclusion.
$\mathcal{F}$ is continuous in $(L^2(0,T;H^1(\Omega)) )^2$ because
its two components $\mathcal{F}_1$ and $\mathcal{F}_2$ are.
 Furthermore, let $A\in \mathbb{R}_+^*$ be the real number defined by
\[
A=\max(C_M,D_M,E_{M}),
\]
and  $W$ be the nonempty (strongly) closed convex bounded set
in $(L^2(0,T;H^1(\Omega)) )^2$ defined by
\begin{align*}
W&=\Big\{(g,g_1)\in \big(L^2(0,T;H^1(\Omega))\cap H^1(0,T;V')\big)^2;
 (g(0),g_1(0))=(h_0,f_{0}),\\
&\quad (g|_{\Gamma},g_1|_{|\Gamma})=(h_D,f_{D});
\|(g,g_1)\|_{(L^2(0,T;H^1(\Omega))\cap H^1(0,T;V'))
\times L^2(0,T;H^1(\Omega))}\leq A\Big\} .
\end{align*}
We have shown that $\mathcal{F}(W)\subset W$.
It follows from Schauder theorem \cite[cor. 9.7]{Z} that there exists
$(h,f)\in W$ such that $\mathcal{F}(h,f)=(h,f)$.
This fixed point for $\mathcal{F}$  is a weak solution of truncated
problem \eqref{tconh}-\eqref{tconf}.
\smallskip

\noindent\textbf{Step 2:} Maximum Principles.
We are going to prove that for almost every $ x \in \Omega$ and for
all $t \in (0,T)$,
 \[
\delta_1 \leq h(t,x)\leq h_2.
\]
First we show that  $ h(t,x)\leq h_2 $ a.e.  $x\in  \Omega$ and for all
$t\in (0,T)$.
 We set
 \[
h_m=\big(h-h_2\big)^{+} = \sup (0,h-h_2) \in L^2(0,T;V).
\]
It satisfies $\nabla h_m=\chi _{\{h>h_2\}}\nabla h$ and $h_m(t,x)\neq 0$
if and only if $h(t,x) > h_2$,
 where $\chi$ denotes the characteristic function.
 Let $\tau\in (0,T)$.
 Taking $w(t,x)=h_m(t,x) \chi_{(0,\tau)}(t)$  in \eqref{tconh} yields
\begin{equation}
\begin{aligned}
&\int_0^{\tau}\phi\langle\partial_th,h_{m} \chi_{(0,\tau)}\rangle_{V',V}dt
+\int_0^\tau\int_\Omega \delta \phi \nabla h\cdot \nabla h_{m}
+\int_0^\tau \int_\Omega KT_s(h)  \nabla h\cdot\nabla h_{m}  \,dx\,dt\\
&+\int_0^\tau\int_\Omega K T_s(h)L_M\big(\|\nabla f\|_{L^2}\big)
 \nabla f\cdot\nabla h_{m}  \,dx\,dt
+\int_0^\tau\int_\Omega Q_sT_s(h) h_{m} \,dx\,dt=0;
\end{aligned}\label{Chap2_42}
\end{equation}
that is,
\begin{equation}
\begin{aligned}
&\int_0^{\tau}\phi\langle\partial_th,h_m\rangle_{V',V}dt
+\int_0^{\tau}\int_\Omega \delta\phi\chi_{\{h>h_2\}}|\nabla h|^2  \,dx\,dt\\
&+\int_0^{\tau}\int_\Omega KT_s(h) \chi_{\{h>h_2\}}|\nabla h|^2  \,dx\,dt\\
&+\int_0^\tau\int_\Omega KT_s(h)L_M\big(\|\nabla f\|_{L^2}\big)
\nabla f\cdot\nabla h_{m}(x,t)  \,dx\,dt\\
&+\int_0^{\tau}\int_\Omega  Q_sT_s(h) h_{m}(x,t) \,dx\,dt=0.
\end{aligned}\label{Chap2_43}
\end{equation}
To evaluate the first term in the left-hand side of \eqref{Chap2_43},
we apply  Lemma 1 with  function $f$ defined by $f(\lambda)=\lambda-h_2$,
$\lambda \in \mathbb{R}$.
We write
\[
\int_0^{\tau}\phi\langle\partial_t h,h_m\rangle_{V',V}dt
= \frac{\phi}{2}\int_\Omega\Big(h_m^2(\tau,x)-h_m^2(0,x)\Big)\, dx
=\frac{\phi}{2}\int_\Omega h_m^2(\tau,x)\,dx,
\]
since $h_m(0,\cdot)=\big(h_0(\cdot)-h_2(\cdot)\big)^{+}=0$.
Moreover $T_s(h)\chi_{\{h>h_2\}}=0$ by definition of $T_s$,
the three last terms  in the left-hand side of \eqref{Chap2_43}  are null.
Hence \eqref{Chap2_43} becomes
\[
\frac{\phi}{2}\int_\Omega h_m^2(\tau,x)\, dx \leq
- \int_0^{\tau}\int_\Omega \delta\phi\chi_{\{h>h_2\}}|\nabla h|^2  \,dx\,dt \leq 0.
\]
 Then $h_m=0$ a.e. in  $\Omega_T$.

 Now we claim that  $  \delta_1 \leq h(t,x)$ a.e.  $x\in  \Omega$ and for all
$t\in (0,T)$. We  set
\[
h_m=\big(h - \delta_1 \big)^{-}\in L^2(0,T;V).
\]
  Let $\tau\in (0,T)$.
  We recall that $h_m(0,\cdot) =0$ a.e. in $\Omega$ thanks to the maximum
 principle satisfied by the initial data $h_0$.
  Moreover,
$\nabla (h -\delta_1)\cdot \nabla h_m=\chi_{\{\delta_1-h> 0\}}
| \nabla (h - \delta_1) |^2$.
  Thus, taking $w(t,x)=h_{m}(x,t) \chi_{(0,\tau)}(t)$ in \eqref{tconh}
  and
\[
w(t,x)=\frac{h_2-\delta_1}{h_2}L_M(\|\nabla f\|_{L^2})
 h_{m}(x,t) \chi_{(0,\tau)}(t)
\] in \eqref{tconf}
 and adding the two equations gives
\begin{align*}
&\int_0^{\tau}  \phi \langle\partial_t h,h_{m}(x,t)
\rangle_{V',V}dt
 +\int_{\Omega_{\tau}}(\delta \phi+K T_s(h))\nabla h\cdot\nabla h_m
 \,dx\,dt\\
& -\int_{\Omega_{\tau} K}T_s(h)L_M(\|\nabla f\|_{L^2})\nabla
  f\cdot\nabla h_m  {dx\,dt}\\
&+\int_{\Omega_T}(h_2-\delta_1)L_M(\|\nabla f\|_{L^2})
 {K}\nabla f\cdot\nabla h_m\,dx\,dt\\
& -\int_{\Omega_{\tau}}T_s(h)\frac{h_2-\delta_1}{h_2}
 L_M(\|\nabla f\|_{L^2})K \nabla h\cdot
 \nabla h_m \,dx\,dt\\
&+\int_{\Omega_{\tau}}\Big({Q}_s T_s(h)\big(1-\frac{h_2-\delta_1}
 {h_2}L_M(\|\nabla f\|_{L^2})\big)h_m\\
&-{Q}_f T_f(h)\frac{h_2-\delta_1}{h_2}L_M(\|\nabla f\|_{L^2}) h_m\Big) \,dx\,dt=0.
\end{align*}
By definition of $T_s(h)$,  $T_s(h)\chi_{\{h<\delta_1\}}=h_2-\delta_1$,
we can simplify the above equation as follows
\begin{align*}
&\frac{\phi}{2}\int_\Omega h_m^2(\tau,x)dx
 +\int_{\Omega_{\tau}}\chi_{\{h<\delta_1\}} \delta
 \phi \nabla h \cdot \nabla h  \,dx\,dt\\
& +\int_{\Omega_{\tau}} {Q}_f T_f(h)\chi_{\{h<\delta_1\}}
  (\delta_1 -h) \frac{h_2-\delta_1}{h_2}L_M(\|\nabla f\|_{L^2}) \,dx\,dt\\
&+\int_{\Omega_{\tau}}(h_2-\delta_1)\big(1-\frac{h_2-\delta_1}{h_2}
 L_M(\|\nabla f\|_{L^2})\big)\chi_{\{h<\delta_1\}} { K \nabla h
 \cdot \nabla h}  \,dx\,dt\\
&+\int_{\Omega_{\tau}}\chi_{\{h<\delta_1\}}(h-\delta_1){Q}_s(h_2-\delta_1)
 \big(1-\frac{h_2-\delta_1}{h_2}L_M(\|\nabla f\|_{L^2})\big)\,dx\,dt=0.
\end{align*}
We first note that
$$
\big(1-\frac{h_2-\delta_1}{h_2}L_M(\|\nabla f\|_{L^2})\big) \geq 0.
$$
Moreover, since $T_f(h)\chi_{\{h<\delta_1\}}=0$ by definition of $T_f$ and
 ${Q}_s \leq 0$, the previous equation leads to
 $$
\frac{1}{2}\int_\Omega h_m^2(\tau,x)dx\leq 0
$$
and then $h_m=0$ a.e. in $\Omega_T$.
\smallskip

\noindent\textbf{Step 3:}  Elimination of the auxiliary function $L_M$.
We now claim that there exist a real number $B>0$, not depending on $M$,
such that any weak solution $(h,f)\in W$  of problem \eqref{tconh}-\eqref{tconf}
satisfies
\begin{equation}
\|\nabla h\|_{L^2(0,T;H)}\leq B\quad \text{and} \quad
 \|\nabla f\|_{L^2(0,T;H)}\leq B.\label{Chap2_35}
\end{equation}
Taking $w=h-h_D$ (resp.  $w=f -f_{D}$) in \eqref{tconh}
(resp. \eqref{tconf}) leads to
\begin{align*}
&\int_0^T\phi\langle\partial_t h,h-h_D\rangle_{V',V} \,dt
+\int_{\Omega_T}{(\delta\,  \phi+T_s(h)K)}
\nabla h\cdot\nabla(h-h_D)\,dx\,dt\\
&=\int_{\Omega_T}
\Big(T_s(h)L_M(\|\nabla f\|_{L^2})))K \nabla f\cdot\nabla(h-h_D)
- {Q}_s T_s(h)(h-h_D)\Big) \,dx\,dt
\end{align*}
and
\begin{align*}
&\int_{\Omega_T}h_2 K \nabla f\cdot\nabla(f-f_D) \,dx dt
-\int_{\Omega_T} T_s(h)K \nabla h\cdot\nabla(f-f_D)\,dx\,dt\\
&=\int_{\Omega_T} ({Q}_s T_s(h)(f-f_D)+{Q}_f T_f(h)(f-f_D))\,dx\,dt.
\end{align*}
Summing up the previous equations yields
\begin{align*}
&\frac{\phi }{2}\int_\Omega [(h-h_D)^2(T,x)-(h-h_D)^2(0,x)]dx
 +\int_{\Omega_T}\delta\, {\phi}|\nabla h|^2 \,dx\,dt\\
&+\int_{\Omega_T} K_-T_s(h)|\nabla(h-f)|^2 \,dx\,dt
+\int_{\Omega_T}K_- (h_2-T_s(h))|\nabla f|^2 \,dx\,dt\\
& +\int_{\Omega_T} K_-T_s(h)(1-L_M(\|\nabla f\|_{L^2}))|\nabla h|^2 \,dx\,dt
\\
&\leq -\int_0^T \phi \langle\partial_t h_D,h-h_D\rangle_{V',V}dt\\
&\quad +\int_{\Omega_T} T_s(h)(1-L_M(\|\nabla f\|_{L^2})) K \nabla h\cdot\nabla(h-f)
 \,dx\,dt\\
&\quad+\int_{\Omega_T}\delta  \phi \nabla h\cdot\nabla h_D \,dx\,dt
 +\int_{\Omega_T} T_s(h) K \nabla (h-f)\cdot\nabla h_D \,dx\,dt\\
&\quad +\int_{\Omega_T} T_s(h)(1-L_M(\|\nabla f\|_{L^2}))
 K \nabla f\cdot\nabla h_D \,dx\,dt\\
&\quad -\int_{\Omega_T}(T_s(h)-h_2) K \nabla f\cdot\nabla f_D \,dx\,dt \\
&\quad -\int_{\Omega_T} T_s(h) K \nabla(h-f)\cdot\nabla f_D \,dx\,dt
  + \int_{\Omega_T} T_s(h) K \nabla h\cdot\nabla f_D \,dx\,dt\\
&\quad +\int_{\Omega_T} {Q}_s T_s(h)[(f-h)+(h_D-f_D)]\,dx\,dt
 +\int_{\Omega_T}{Q}_f h(f-f_D)\,dx\,dt\\
&:= I_1+I_2+I_3+I_4+I_5+I_6+I_7+I_8+I_9+I_{10}.
\end{align*}
Using the Cauchy-Schwarz and Young inequalities, we obtain
\begin{align*}
|I_1|
&\leq  {\phi}\|\partial_t h_D\|_{L^2(0,T,V')}\|h-h_D\|_{L^2(0,T,V)}\\
&\leq  \phi \|\partial_t h_D\|_{L^2(0,T,V')}
 (\|\nabla h\|_{L^2(0,T;L^2(\Omega))}+\|\nabla h_D\|_{L^2(0,T;L^2(\Omega))})\\
&\leq \frac{\delta \phi }{6}\|\nabla h\|^2_{L^2(0,T;L^2(\Omega))}
 +\frac{3  \phi }{2\delta}\|\partial_t h_D\|^2_{L^2(0,T,V')}
+\frac{ {\phi }}{2}\|\partial_t h_D\|^2_{L^2(0,T,V')}\\
&\quad  + \frac{ \phi}{2} \|\nabla h_D\|^2_{L^2(0,T,V')},
\end{align*}
\begin{align*}
|I_2|&\leq K_+ \Big(\int_{\Omega_T} T_s(h)(1-L_M(\|\nabla f\|_{L^2}))|\nabla h|^2
 \,dx\,dt\Big)^{1/2}\\
&\quad\times \Big(\int_{\Omega_T} T_s(h)|\nabla (h-f)|^2 \,dx\,dt\Big)^{1/2}\\
&\leq \frac{K_+}{2}\Big(\int_{\Omega_T} T_s(h)(1-L_M(\|\nabla
f\|_{L^2}))|\nabla h|^2 \,dx\,dt\Big)\\
&\quad +\frac{K_+}{2}\Big(\int_{\Omega_T} T_s(h)|\nabla (h-f)|^2 \,dx\,dt\Big),
\end{align*}
\begin{align*}
|I_3|&\leq  {\phi \,}\delta\Big(\int_{\Omega_T} |\nabla h|^2 \,dx\,dt\Big)^{1/2}
\Big(\int_{\Omega_T}|\nabla h_D|^2 \,dx\,dt\Big)^{1/2}\\
&\leq \frac{\delta  \phi}{6}\int_{\Omega_T}|\nabla h|^2 \,dx\,dt
+\frac{3 \delta{\phi}}{ 2 }\int_{\Omega_T}|\nabla h_D|^2 \,dx\,dt,
\end{align*}
\begin{align*}
|I_4|&\leq K_+ \Big(\int_{\Omega_T} T_s(h)|\nabla (h-f)|^2 \,dx\,dt\Big)^{1/2}
\Big(\int_{\Omega_T}h_2|\nabla h_D|^2 \,dx\,dt\Big)^{1/2}\\
&\leq\frac{K_-}{12}\int_{\Omega_T}T_s(h)|\nabla (h-f)|^2 \,dx\,dt
+\frac{6K_+^2 h_2}{2K_-}\int_{\Omega_T}|\nabla h_D|^2 \,dx\,dt,
\end{align*}
\begin{align*}
|I_5|&\leq K_+ \sqrt{h_2}\Big(\Big(\int_{\Omega_T}T_s(h) |\nabla (h-f)|^2
\,dx\,dt\Big)^{1/2}\\
&\quad +\Big(\int_{\Omega_T}T_s(h)(1-L_M(\|\nabla f\|_{L^2}))
 |\nabla h|^2 \,dx\,dt\Big)^{1/2}{\Big)}(\int_{\Omega_T}|\nabla h_D|^2 \,dx\,dt)^{1/2}\\
&\leq \frac{K_-}{12}\int_{\Omega_T}T_s(h) |\nabla (h-f)|^2 \,dx\,dt\\
&\quad+\frac{K_-}{12}\int_{\Omega_T}T_s(h)(1-L_M(\|\nabla f\|_{L^2})) |\nabla h|^2 \,dx\,dt
 +\frac{6 h_2 K_+^2}{K_-}\int_{\Omega_T}|\nabla h_D|^2 \,dx\,dt,
\end{align*}
\begin{align*}
|I_6|&\leq K_+ \sqrt{h_2}\Big(\int_{\Omega_T} (h_2-T_s(h))|\nabla f|^2 \,dx\,dt
 \Big)^{1/2}\Big(\int_{\Omega_T}|\nabla f_D|^2 \,dx\,dt\Big)^{1/2}\\
&\leq \frac{K_-}{6}\int_{\Omega_T}(h_2-T_s(h))|\nabla f|^2 \,dx\,dt
+\frac{3h_2K_+^2}{2K_-}\int_{\Omega_T}|\nabla f_D|^2 \,dx\,dt,
\end{align*}
\begin{align*}
|I_7|&\leq K_+\sqrt{h_2}\Big(\int_{\Omega_T} T_s(h)|\nabla(h- f)|^2 \,dx\,dt
\Big)^{1/2}\Big(\int_{\Omega_T}|\nabla f_D|^2 \,dx\,dt\Big)^{1/2}\\
&\leq \frac{K_-}{12}\int_{\Omega_T}T_s(h)|\nabla (h-f)|^2 \,dx\,dt
+\frac{6h_2K_+^2}{2K_-}\int_{\Omega_T}|\nabla f_D|^2 \,dx\,dt,
\end{align*}
\begin{align*}
|I_8|&\leq K_+ h_2(\int_{\Omega_T} |\nabla h|^2 \,dx\,dt)^{1/2}
\Big(\int_{\Omega_T}|\nabla f_D|^2 \,dx\,dt\Big)^{1/2}\\
&\leq \frac{\delta\,  \phi}{6}\int_{\Omega_T}|\nabla h|^2 \,dx\,dt
+\frac{3h_2^2K_+^2}{2\delta\,  \phi}\int_{\Omega_T}|\nabla f_D|^2 \,dx\,dt,
\end{align*}
\begin{align*}
|I_9|&\leq h_2\|{Q}_s\|_{L^2(0,T;H)}(\|h-h_D\|_{L^2(0,T;H)}+\|f-f_D\|_{L^2(0,T;H)})\\
&\leq h_2\|{Q}_s\|_{L^2(0,T;H)}(\|h-h_D\|_{L^2(0,T;H)}
 +C_p\|\nabla(f-f_D)\|_{L^2(0,T;H)})\\
&\leq \frac{\phi}{2}\int_{\Omega_T}(h-h_D)^2 \,dx\,dt
 +\frac{h_2^2}{2\phi}\|{Q}_s\|^2_{L^2(0,T;H)}\\
&\quad +\frac{\delta_1K_-}{6}\int_{\Omega_T}|\nabla f|^2 \,dx\,dt
 +\frac{3C_p^2h_2^2}{2\delta_1K_-}\|{Q}_s\|_{L^2(0,T;H)}\\
&\quad +h_2\|{Q}_s\|_{L^2(0,T;H)}C_p\|\nabla f_D\|_{L^2(0,T;H)},
\end{align*}
\begin{align*}
|I_{10}|&\leq h_2\|{Q}_s\|_{L^2(0,T;H)}\|f-f_D\|_{L^2(0,T;H)}\\
&\leq \frac{\delta_1K_-}{6}\int_{\Omega_T}|\nabla f|^2 \,dx\,dt
 +\frac{3C_p^2h_2^2}{2\delta_1K_-}\|{Q}_s\|_{L^2(0,T;H)}\\
&\quad +h_2\|{Q}_s\|_{L^2(0,T;H)}C_p\|\nabla f_D\|_{L^2(0,T;H)}.
\end{align*}
Since $ \delta_1 \leq  h \leq  h_2$ a.e. in $\Omega_T$, it follows that
 $ T_s(h)-h_2=h \geq   \delta_1$ and then
\begin{align*}
&\frac{\phi}{2}\int_{\Omega_T}(h-h_D)^2(\tau,x)dx
 +\frac{\delta \phi}{2}\int_{\Omega_T}|\nabla h|^2 \,dx\,dt
 +\frac{K_-\delta_1}{2}\int_{\Omega_T}|\nabla f|^2 \,dx\,dt\\
&+\int_{\Omega_T} \frac{(3K_- - 2K_+)}{4}T_s(h)|\nabla(h-f)|^2 \,dx\,dt\\
& +\int_{\Omega_T} \frac{11K_- - 6 K_+}{12}T_s(h)
 (1-L_M(\|\nabla f\|_{L^2}))|\nabla h|^2 \,dx\,dt\\
&\leq \frac{\phi}{2}\int_{\Omega_T}(h-h_D)^2 \,dx\,dt
 +  C(h_2,\delta,\delta_1,h_0,h_D,f_D ,Q_f,Q_s)
\end{align*}
Now we apply the Gronwall lemma and we deduce that  there exists a real number
$B$, that does not depend on  $M$, such that
 \[
\|h\|_{L^\infty(0,T;H)\cap L^2(0,T;H^1(\Omega))}\leq B\quad \text{and} \quad
\|f\|_{L^2(0,T;H^1(\Omega))}\leq B .
\]
In particular, $\|\nabla f\|_{L^2(\Omega_T)^2}\leq B$ and this estimate
does not depend on the choice of the real number $M$ that defines  function $L_M$.
Hence if we choose $M=B$,  any weak solution of the system
\begin{gather*}
\begin{aligned}
&\phi\partial_t h -\nabla \cdot (\delta\, {\phi}\nabla h)-\nabla \cdot
\big( KT_s(h) \nabla h\big)+\nabla \cdot \big(KT_s(h)
L_B(\|\nabla f\|_{L^2})\nabla f\big)\\
&= - Q_sT_s(h), 
\end{aligned}\\
 -\nabla \cdot \big(Kh_2\nabla f\big)  +\nabla \cdot \big(KT_s(h)\nabla h\big)
=Q_fT_f(h-h_1) + Q_sT_s(h)
\end{gather*}
in $\Omega_T$, with the initial and boundary conditions
$h=h_D$ and $f=f_{D}$ on $\Gamma$,
$h(0,x)=h_0$ a.e. in $\Omega$,
satisfies $L_B\big(\|\nabla f\|_{L^2}\big)=1$.
Then the term $L_B\big(\|\nabla f\|_{L^2}\big)=1$ may be dropped.
The proof of Theorem  \ref{Theo1} is complete.

\subsection{Proof of Theorem \ref{Theo2}}

Let $\epsilon >0$  and let $M$ be a positive constant to be determined later.
For $x \in \mathbb{R}_+^*$, we set
\[
L_M(x)=\min\Big(1,\frac{M}{x}\Big) .
\]
Such a truncation $L_M$ allows again to  use the following point in the
estimates hereafter. For any $(g,g_1)\in (L^2(0,T;H^1(\Omega)))^2$, setting
$$
d(g,g_1)=-T_s(g)L_M\big(\|\nabla g_1\|_{L^2(\Omega_T)}\big)\nabla g_1,
$$
we have
$$
\| d(g,g_1)\|_{L^2(0,T;H)}=
 \| T_s(g)L_M\big(  \|\nabla g_1  \|_{L^2(\Omega _T)}\big)\nabla g_1  \|_{H}\leq
 M   h_2  .
$$
We also define a regularized step function for $\mathcal{X}_0$ by
 \[
\mathcal{X}_0(h_1)=\begin{cases}
0& \text{if }  h_1\leq0\\
1& \text{if }   h_1>0\,,
\end{cases}\quad
\mathcal{X}_0^\epsilon(h_1)=\begin{cases}
0& \text{if }   h_1\leq0\\
 h_1/(h_1^2+\epsilon)^{1/2} &\text{if }  h_1>0.
\end{cases}
 \]
Then $ 0\leq\mathcal{X}_0^\epsilon \leq 1 $ and
$\mathcal{X}_0^\epsilon \to \mathcal{X}_0$ as $\epsilon \to 0$.
Introducing the regularization  $ \mathcal{X}_0^\epsilon$ of $\mathcal{X}_0$,
we replace system \eqref{problem} by the system
\begin{gather*}
\begin{aligned}
&\phi\partial_t h^\epsilon -\nabla \cdot (\epsilon \nabla h^\epsilon)
-\nabla \cdot \big( KT_s(h^\epsilon) \nabla h^\epsilon\big)\\
&-\nabla \cdot \big(KT_s(h^\epsilon)\mathcal{X}_0^\epsilon(h_1^\epsilon)
L_M\big(\|\nabla h_1^\epsilon\|_{L^2}\big)\nabla h_1^\epsilon\big)\\
&=-Q_sT_s(h^\epsilon) ,
\end{aligned} \\
\begin{aligned}
&\phi \partial_th_1^\epsilon  -\nabla \cdot (\epsilon \nabla h_1^\epsilon)
-\nabla \cdot \Big(K\big(T_f(h^\epsilon-h_1^\epsilon)
+T_s(h^\epsilon)\big)\nabla h_1^\epsilon\Big)\\
& -\nabla \cdot \big(KT_s(h^\epsilon)\mathcal{X}_0^\epsilon(h_1^\epsilon)
\nabla h^\epsilon\big)
 \\
&=-Q_fT_f(h^\epsilon-h_1^\epsilon)\mathcal{X}_0^\epsilon(h_1^\epsilon)
-Q_sT_s(h^\epsilon)\mathcal{X}_0^\epsilon(h_1^\epsilon).
\end{aligned}
\end{gather*}
The proof is outlined as follows: In the first step, using Schauder theorem,
we prove that for every $T>0$
and every $\epsilon>0$, the above regularized system completed by the initial
and boundary conditions
\begin{gather*}
 h^\epsilon=h_D, \quad   h_1^\epsilon=h_{1,D} \quad \text{in } \Gamma \times (0,T),\\
h^\epsilon(0,x)=h_0(x),\quad
h_1^\epsilon(0,x)=h_{1,0}(x) \quad \text{a.e. in } \Omega ,
\end{gather*}
has a solution $(h^{\epsilon}-h_D, h_1^{\epsilon}-h_{1,D})\in W(0,T)\times W(0,T)$.
We observe that the sequence ($h^\epsilon-h_D,\, h_1^{\epsilon}-h_{1,D}) $
is uniformly bounded in $(L^2(0,T;V ))^2$ and we show the maximum principle
$ 0\leq  h_1^{\epsilon}(t,x)$ and $0 \leq  h^{\epsilon}(t,x)\leq h_2$ a.e.
in $\Omega _T$ for every $\epsilon > 0$
Finally we prove that any (weak) limit  $(h,h_1)$ in
$(L^2(0,T;H^1(\Omega)) \cap H^1(0,T;V'))^2$ of the sequence
$(h^\epsilon,\, h_1^{\epsilon})$ is a solution of the original problem.
\smallskip


\noindent\textbf{Step 1:} Existence for the regularized system.
We now omit $\epsilon$ for the sake of simplicity in the notation.
Then the weak formulation of the latter problem reads: for any $w\in V$,
\begin{gather}
\begin{aligned}
&\int_0^T\phi\langle\partial_t h,w\rangle_{V',V}dt
 +\int_{\Omega_T} \epsilon \nabla h \cdot \nabla w \,dx\,dt\\
&+\int_{\Omega_T} KT_s(h)\nabla h \cdot\nabla w \,dx\,dt 
+\int_{\Omega_T} KT_s(h)\mathcal{X}_0^\epsilon(h_1)
 L_M\big(\|\nabla h_1\|_{L^2}\big)\nabla h_1 \cdot\nabla w \,dx\,dt\\
& +\int_{\Omega_T} Q_sT_s(h)w \,dx\,dt=0,
\end{aligned}\label{Chap2_1}
\\
\begin{aligned}
& \int_0^T\phi\langle\partial_t h_1,w\rangle_{V',V}dt
 +\int_{\Omega_T} \epsilon\nabla h_1 \cdot\nabla w \,dx\,dt\\
& + \int_{\Omega_T} K\Big(T_s(h)+T_f(h-h_1)\Big)\nabla h_1\cdot\nabla w \, dx\, dt
 +\int_{\Omega_T} K\mathcal{X}_0^\epsilon(h_1)T_s(h)\nabla h\cdot
\nabla w\,dx\,dt\\
& +\int_{\Omega_T}\Big(Q_fT_f(h-h_1)+Q_sT_s(h)\Big)
\mathcal{X}_0^\epsilon(h_1^\epsilon)w \,dx\,dt =0.
\end{aligned} \label{Chap2_2}
\end{gather}
For the  fixed point strategy, we define the application $\mathcal{F}$ by
\begin{gather*}
\mathcal{F}:L^2(0,T;H^1(\Omega))\times L^2(0,T;H^1(\Omega))
\to L^2(0,T;H^1(\Omega))\times L^2(0,T;H^1(\Omega))\\
 (\bar h,\bar h_1)\longmapsto\mathcal{F}(\bar h,\bar h_1)
=\Big(\mathcal{F}_1(\bar h,\bar h_1)=h,\mathcal{F}_2(\bar h,\bar h_1)=h_1\Big) ,
\end{gather*}
where $( h,h_1)$ is the solution of
\begin{gather}
\begin{aligned}
&\int_0^T\phi\langle\partial_t h,w\rangle_{V',V}dt
+\int_{\Omega_T} \epsilon \nabla h \cdot\nabla w \,dx\,dt
+\int_{\Omega_T} KT_s(\bar h)\nabla h \cdot\nabla w \,dx\,dt \\
&+\int_{\Omega_T} KT_s(\bar h)L_M\big(\|\nabla \bar h_1\|_{L^2}\big)
 \mathcal{X}_0^\epsilon(\bar h_1)\nabla \bar h_1\cdot \nabla w \,dx\,dt
+\int_{\Omega_T} Q_sT_s(\bar h)w \,dx\,dt\\
&=0 \quad \forall w\in V,
\end{aligned}\label{Chap2_3}
\\
\begin{aligned}
&\int_0^T\phi\langle\partial_t h_1,w\rangle_{V',V}dt
+\int_{\Omega_T} \epsilon \nabla h_1\cdot\nabla w \,dx\,dt\\
&+ \int_{\Omega_T} K\Big(T_s(\bar h)+T_f(\bar h-\bar h_1)\Big)
\nabla h_1\cdot\nabla w \, dx\, dt  
+\int_{\Omega_T} KT_s(\bar h)\mathcal{X}_0^\epsilon(\bar h_1)
 \nabla h\cdot\nabla w\,dx\,dt\\
&+\int_{\Omega_T}\Big(Q_fT_f(\bar h-\bar h_1)
+Q_sT_s(\bar h)\Big)\mathcal{X}_0^\epsilon(\bar h_1)w \,dx\,dt=0\quad
 \forall w\in V.
\end{aligned} \label{Chap2_4}
\end{gather}
Indeed we know from  classical  parabolic theory  (see {\it e.g.}
\cite{LM}) that the linear variational system \eqref{Chap2_3}-\eqref{Chap2_4}
admits an unique solution.
It remains to prove a fixed point property for application $\mathcal{F}$.

Since the proofs of the continuity of $\mathcal{F}_1$ and $\mathcal{F}_2$
are very similar to the previous ones, we do not reproduce here the computations.
We also emphasize that this step is proven in \cite{CDR2} by considering that
 the parameter $ \epsilon$ plays the same role as $\delta$, the thickness
of the diffuse interface, of course the estimates depend on $ \epsilon$,
but it is sufficient for this step.
We can thus conclude that $\mathcal{F}$ is continuous in
 $(L^2(0,T;H^1(\Omega)) )^2$ because its two components $\mathcal{F}_1$
and $\mathcal{F}_2$ are. Furthermore, there exist a real number
$A\in \mathbb{R}_+^*$ (depending on the data and on the
parameter $ \epsilon$) and a non empty (strongly) closed convex bounded
set $W$ in $(L^2(0,T;H^1(\Omega)) )^2$ defined by
\begin{align*}
W=\Big\{&(g,g_1)\in \big(L^2(0,T;H^1(\Omega))\cap H^1(0,T;V')\big)^2;
(g(0),g_1(0))=(h_0,h_{1,0}), \\
& (g|_{\Gamma},g_1|_{\Gamma})=(h_D,h_{1,D});
\|(g,g_1)\|_{(L^2(0,T;H^1(\Omega))\cap H^1(0,T;V')))^2}\leq A\Big\} .
\end{align*}
such that $\mathcal{F}(W)\subset W$.
It follows from Schauder theorem \cite[cor. 9.7]{Z} that there exists
$(h,h_1)\in W$ such that $\mathcal{F}(h,h_1)=(h,h_1)$.
This fixed point for $\mathcal{F}$  is a weak solution
of problem \eqref{Chap2_1}-\eqref{Chap2_2}.
\smallskip


\noindent\textbf{Step 2:} Maximum Principles.
We aim at proving that for almost all $x\in \Omega$ and for all $t\in(0,T)$,
$0\leq h_1(t,x)$ and $0\leq h(t,x)\leq h_2$.

  First we show that  $ h(t,x) \leq h_2 $ a.e.  $x\in  \Omega$ and  all $t\in (0,T)$.
 We set
 \[
h_m=\big(h-h_2\big)^{+} = \sup (0,h-h_2) \in L^2(0,T;V).
\]
It satisfies $\nabla h_m=\chi _{\{h>h_2\}}\nabla h$ and
$h_m(t,x)\neq 0$ if and only if $h(t,x) > h_2$,
 where $\chi$ denotes the characteristic function.
 Let $\tau\in (0,T)$.
 Taking $w(t,x)=h_m(t,x) \chi_{(0,\tau)}(t)$  in \eqref{Chap2_1} yields
\begin{equation}
\begin{aligned}
&\int_0^{\tau}\phi\langle\partial_t h,h_m\rangle_{V',V}dt
+\int_0^{\tau}\int_\Omega  \epsilon \chi_{\{h>h_2\}}|\nabla h|^2  \,dx\,dt\\
&+\int_0^{\tau}\int_\Omega K_- T_s(h) \chi_{\{h>h_2\}}|\nabla h|^2\,dx\,dt \\
&+\int_0^\tau\int_\Omega \mathcal{X}_0^\epsilon(h_1^\epsilon)
\chi_{\{h>h_2\}}KT_s(h)L_M\big(\|\nabla {  h_1}\|_{L^2}\big) \nabla
{ h_1} \cdot\nabla h(x,t)  \,dx\,dt\\
&+\int_0^{\tau}\int_\Omega  Q_sT_s(h) h_{m}(x,t) \,dx\,dt  \leq 0.
\end{aligned} \label{ppe1}
\end{equation}
Lemma 1 applied to \eqref{ppe1},  gives
\[
\int_0^{\tau}\phi\langle\partial_t h,h_m\rangle_{V',V}dt=
\frac{\phi}{2}\int_\Omega\Big(h_m^2(\tau,x)-h_m^2(0,x)\Big)\, dx
=\frac{\phi}{2}\int_\Omega h_m^2(\tau,x)\,dx,
\]
taking into account $h_m(0,\cdot)=\big(h_0(\cdot)-h_2(\cdot)\big)^{+}=0$.
Since $T_s(h)\chi_{\{h>h_2\}}=0$ by definition of $T_s$,
the three last terms of \eqref{ppe1} are null.
Hence \eqref{ppe1} becomes
\[
\frac{\phi}{2}\int_\Omega h_m^2(\tau,x)\, dx
\leq - \int_0^{\tau}\int_\Omega \epsilon \chi_{\{h>h_2\}}|\nabla h|^2  \,dx\,dt
\leq 0.
\]
 Then $h_m=0$ a.e. in  $\Omega_T$.

Now we claim that  $  0 \leq h(t,x)$ a.e.  $x\in  \Omega$ and all $t\in (0,T)$.
We set
\[
h_m=\big(-h  \big)^{+}\in L^2(0,T;V) \quad \text{since }
 h_D(\cdot,\cdot) \geq 0.
\]
Let $\tau\in (0,T)$.  We recall that $h_m(0,\cdot) =0$ a.e. in $\Omega$
thanks to the maximum principle satisfied by the initial data $h_0$.
Moreover, $ \nabla h_m=-\chi_{\{h <0\}} \nabla h $.
Thus, taking $w(t,x)=h_{m}(x,t) \chi_{(0,\tau)}(t)$ in \eqref{Chap2_1}
yields
\begin{equation}
\begin{aligned}
&\int_0^{\tau}\phi\langle\partial_t h,h_m\rangle_{V',V}dt
-\int_0^{\tau}\int_\Omega  \epsilon \chi_{\{h<0\}}|\nabla h|^2  \,dx\,dt\\
&=\int_0^{\tau}\int_\Omega  T_s(h) \chi_{\{h<0\}} { K \nabla h \cdot \nabla h}
\,dx\,dt \\
&\quad + \int_0^\tau\int_\Omega \mathcal{X}_0^\epsilon(h_1)
 \chi_{\{h <0\}} KT_s(h)L_M\big(\|\nabla  h_1\|_{L^2}\big)
 \nabla h_1\cdot\nabla h(x,t)  \,dx\,dt\\
&\quad -\int_0^{\tau} \int_\Omega  Q_sT_s(h) h_{m}(x,t) \,dx\,dt.
\end{aligned} \label{ppe2}
\end{equation}
Applying Lemma 1 to the first term of \eqref{ppe2} and taking into
account $h_m(0,\cdot)=\big(-h_0)^{+}=0$,  gives
\[
\int_0^{\tau}\phi\langle\partial_t h,h_m\rangle_{V',V}dt
= \frac{-  \phi}{2}\int_\Omega\Big(h_m^2(\tau,x)-h_m^2(0,x)\Big)\, dx
=-\frac{\phi}{2}\int_\Omega h_m^2(\tau,x)\,dx.
\]
Since $T_s(h)\chi_{\{h < 0\}}=0$ by definition of $T_s$,  the three
last terms  in the left-hand side of \eqref{ppe2}  are null.
Hence \eqref{ppe2} becomes
\[
\frac{\phi}{2}\int_\Omega h_m^2(\tau,x)\, dx \leq
- \int_0^{\tau}\int_\Omega \epsilon \chi_{\{h>h_2\}}|\nabla h|^2  \,dx\,dt \leq 0.
\]
 Then $h_m=0$ a.e. in  $\Omega_T$.


 We finally prove that  $  0 \leq h_1(t,x)$ a.e.  $x\in  \Omega$ and all $t\in (0,T)$.
We set
\[
h_m=\big(-h_1 \big)^{+}\in L^2(0,T;V) \quad \text{since }
 h_{1,D}(\cdot,\cdot) \geq 0.
\]
 Let $\tau\in (0,T)$. We recall that $h_m(0,\cdot) =0$ a.e. in $\Omega$
thanks to the maximum principle satisfied by the initial data $h_{1,0}$.
  Moreover, $ \nabla h_m=-\chi_{\{h_1 <0\}} \nabla  h_1 $.
  Thus, taking $w(t,x)=h_{m}(x,t) \chi_{(0,\tau)}(t)$ in \eqref{Chap2_2} yields
\begin{equation}
\begin{aligned}
&\int_0^{\tau}\phi\langle\partial_t h_1,h_m\rangle_{V',V}dt
-\int_{\Omega_{\tau}}\chi_{\{h_1<0\}}(\epsilon+T_s(h)+T_f(h-h_1))
\nabla h_1 \cdot \nabla h_1 \,dx\,dt \\
&-\int_{\Omega_{\tau}} \mathcal{X}_0^\epsilon(h_1)T_s(h)\chi_{\{h_1<0\}}
 \nabla h\cdot\nabla  h_1\,dx\,dt\\
&-\int_{\Omega_{\tau}}\big({Q}_s T_s(h)+{Q}_f T_f(h-h_1)\big)
 \mathcal{X}_0^\epsilon(h_1)\chi_{\{h_1<0\}}h_1 \,dx\,dt=0.
\end{aligned} \label{ppe3}
\end{equation}
Applying Lemma 1 to \eqref{ppe3} and taking into account
$h_m(0,\cdot)=\big(-h_0)^{+}=0$,  gives
\[
\int_0^{\tau}\phi\langle\partial_t h_1,h_m\rangle_{V',V}dt
= -\frac{\phi}{2}\int_\Omega\Big(h_m^2(\tau,x)-h_m^2(0,x)\Big)\, dx
=-\frac{\phi}{2}\int_\Omega h_m^2(\tau,x)\,dx.
\]
Since $\mathcal{X}_0^\epsilon(h_1)\chi_{\{h_1 < 0\}}=0$ by definition
of $\mathcal{X}_0^\epsilon(\cdot)$,  the two last terms of \eqref{ppe3}  are null.
Hence \eqref{ppe3} becomes (since $T_s$ and  $T_f$ are nonnegative functions)
\[
\frac{\phi}{2}\int_\Omega h_m^2(\tau,x)\, dx
\leq - \int_0^{\tau}\int_\Omega \epsilon \chi_{\{h>h_2\}}|\nabla h|^2  \,dx\,dt
\leq 0.
\]
Then $h_m=0$ a.e. in  $\Omega_T$.
\smallskip

\textbf{Step 3:} Elimination of the auxiliary function $L_M$.
 We now claim that there exist a real number $B>0$, not depending on
$\epsilon$ neither on $M$,   such that any weak solution $(h,h_1)\in W$
of problem \eqref{Chap2_1}-\eqref{Chap2_2} satisfies
\begin{equation}
\|\nabla h\|_{L^2(0,T;H)}\leq B\quad \text{and} \quad
\|\nabla h_{1}\|_{L^2(0,T;H)}\leq B.\label{Chap2_35b}
\end{equation}
Taking $w=h-h_D$ (resp.  $w=h_{1}-h_{1,D}$) in \eqref{Chap2_1}
 (resp. \eqref{Chap2_2}) leads to
   \begin{equation}
\begin{aligned}
&\int_0^T\phi\langle\partial_th,h-h_D\rangle_{V',V}dt
+\int_{\Omega_T} \epsilon\nabla h \cdot \nabla (h-h_D )\,dx\,dt\\
&+\int_{\Omega_T} K T_s(h)\nabla h\cdot\nabla (h-h_D )\,dx\,dt\\
&=-\int_{\Omega_T} KT_s(h) \mathcal{X}_0^\epsilon(h_{1}
 )L_M\big(\|\nabla h_{1}\|_{L^2}\big)
\nabla h_{1}\cdot\nabla (h-h_D )\,dx\,dt\\
&-\int_{\Omega_T} Q_s T_s(h)(h-h_D )\,dx\,dt
\end{aligned}\label{Chap2_36}
\end{equation}
and
\begin{equation}
\begin{aligned}
&\int_0^T\phi \langle\partial_th_{1},h_{1}-h_{1,D}\rangle_{V',V}dt
+ \int_{\Omega_T} \epsilon\nabla h_{1}\cdot\nabla(h_{1}-h_{1,D})\,dx\,dt \\
&+\int_{\Omega_T} K\big(T_s(h)+T_f(h-h_{1})\big)
 \nabla h_{1}\cdot\nabla (h_{1}-h_{1,D}) \, dx dt\\
&=-\int_{\Omega_T} KT_s(h) {\mathcal{X}_0^\epsilon(h_{1})}
\nabla h\cdot\nabla (h_{1}-h_{1,D})\,dx\,dt\\
&\quad -\int_{\Omega_T}\mathcal{X}_0^\epsilon(h_{1})\big(Q_fT_f(h-h_{1})
+Q_sT_s(h)\big)(h_{1}-h_{1,D}) \,dx\,dt.
\end{aligned} \label{Chap2_37}
\end{equation}
 Summing up relations \eqref{Chap2_36} and \eqref{Chap2_37}, and using
the decomposition
\begin{equation}
\begin{aligned}
&K\nabla h\cdot\nabla h+   { K\mathcal{X}_0^\epsilon(h_{1})
\big(L_M\big(\|\nabla h_{1}\|_{L^2}\big)+1\big)}\nabla h_1\cdot\nabla h
+K\nabla h_1\cdot\nabla h_1\\
&=K\nabla(h+h_1)\cdot\nabla(h+h_1)
 +K\Big(1-\mathcal{X}_0^\epsilon(h_{1})L_M\big(\|\nabla h_{1}\|_{L^2}\big)\Big)
\nabla h_1\cdot\nabla h_1 \\
&\quad - K \Big(1-\mathcal{X}_0^\epsilon(h_{1})
L_M\big(\|\nabla h_{1}\|_{L^2}\big)\Big)\nabla h_1\cdot\nabla(h+h_1),
\end{aligned} \label{Chap2_38}
\end{equation}
we write
\begin{align*}
&\underbrace{\int_0^T\phi \Bigl( \langle\partial_t(h-h_D),h-h_D\rangle_{V',V}dt
+ \langle\partial_t(h_{1}-h_{1,D}),h_{1}-h_{1,D}\rangle_{V',V}\Bigr) \, dt}_{J_1}\\
&+\underbrace{\int_{\Omega_T} \epsilon\big(\nabla h\cdot\nabla h
 +\nabla h_{1}\cdot\nabla h_1\big)\,dx\,dt}_{J_2}
+\underbrace{\int_{\Omega_T} K T_s(h)\nabla (h+h_1)\cdot\nabla (h+h_1)\,dx\,dt}_{J_3}
\\
&+\underbrace{\int_{\Omega_T} K \Big(\big(1-\mathcal{X}_0^\epsilon(h_{1})
 L_M(\|\nabla h_{1}\|_{L^2})\big)T_s(h) +T_f(h-h_{1})\Big)
 \nabla h_1\cdot\nabla h_1\,dx\,dt}_{J_4}\\
&=\underbrace{\int_{\Omega_T}  { K }\big(1-\mathcal{X}_0^\epsilon(h_{1})
 L_M(\|\nabla h_{1}\|_{L^2})\big)T_s(h)\nabla h_1\cdot\nabla (h+h_1)\,dx\,dt}_{J_5}
\\
&\quad +\underbrace{\int_{\Omega_T}\big(\epsilon+K T_s(h)\big)\nabla h
 \cdot\nabla h_D\,dx\,dt}_{J_6}
\\
&\quad +\underbrace{\int_{\Omega_T}\big(\epsilon+KT_s(h)
 +KT_f(h-h_{1})\big)\nabla h_{1}\cdot\nabla h_{1,D}\,dx\,dt}_{J_7}\\
&\quad +\underbrace{\int_{\Omega_T} KT_s(h)L_M(\|\nabla h_{1}\|_{L^2})
\mathcal{X}_0^\epsilon(h_{1})\nabla h_{1}\cdot\nabla h_D\,dx\,dt}_{J_8}
\\
&\quad +\underbrace{\int_{\Omega_T} KT_s(h) {\mathcal{X}_0^\epsilon(h_{1})}
\nabla h\cdot\nabla h_{1,D}\,dx\,dt}_{J_{9}}\\
&\quad -\underbrace{\int_{\Omega_T} \bigl( Q_s T_s(h)(h-h_D )
+ {\mathcal{X}_0^\epsilon(h_{1})}
\big(Q_fT_f(h-h_{1})+Q_sT_s(h)\big)(h_{1}-h_{1,D})\bigr) \,dx\,dt}_{J_{10}}
 \\
&\quad -\underbrace{\int_0^T\phi\bigl(\langle\partial_th_D,h-h_D\rangle_{V',V}
+\langle\partial_th_{1,D},h_{1}-h_{1,D}\rangle_{V',V}\bigr)\, dt}_{J_{11}} .
\end{align*} %\label{Chap2_40}
We now estimate all the terms in the above expression.
We recall that
\[
|\mathcal{X}_0^\epsilon(h_{1}| \leq 1, \quad
L_M\big(\|\nabla h_{1}\|_{L^2}\big) \leq 1, \quad
0 \leq T_s(h) \leq h_2 , \quad
\delta_1 \leq T_f(h-h_1) \leq h_2.
\]
Then, we note that
\begin{gather*}
\begin{aligned}
|J_1|&=\frac{\phi}{2}\int_\Omega\big((h-h_{D})^2(T,x)-(h_0-h_{0,D})^2(x)\big)\, dx\\
&\quad + \frac{\phi}{2}\int_\Omega\big((h_1-h_{1,D})^2(T,x)-(h_{1,0}-h_{0,D})^2(x)
 \big)\, dx,
\end{aligned}\\
|J_2|= \int_{\Omega_T} \epsilon |\nabla h|^2 \,dx\,dt
 +\int_{\Omega_T} \epsilon |\nabla h_1|^2 \,dx\,dt ,\\
|J_3|\geq \int_{\Omega_T} K_- T_s(h)|\nabla(h+h_{1})|^2 \,dx\,dt , \\
|J_4|\geq \int_{\Omega_T}  { K_-} \Big(\big(1-\mathcal{X}_0^\epsilon(h_{1})
L_M(\|\nabla h_{1}\|_{L^2})\big)T_s(h)+ \delta_1\Big) |\nabla h_{1}|^2\,dx\,dt .
\end{gather*}
Next, applying  the Cauchy-Schwarz and Young inequalities,
 for some real $\gamma > 0$, we obtain
\begin{gather*}
 \begin{aligned}
|J_5|& \leq \int_{\Omega_T}T_s(h)\, \Bigl( {\frac{1}{4 \gamma}}
\Big(1-\mathcal{X}_0^\epsilon(h_{1})L_M\big(\|\nabla h_{1}\|_{L^2}\big)\Big)^2\\
&\quad\times  \frac{K_+^2}{K_-} |\nabla h_{1}|^2 +
 {\gamma K_-}|\nabla (h+h_{1})|^2  \Bigr)\,dx\,dt ,\\
\end{aligned} \\
\begin{aligned}
|J_6|& \leq \frac{\epsilon}{2}\int_{\Omega_T}   |\nabla h|^2\,dx\,dt
+\frac{\epsilon}{2} \int_{\Omega_T}   |\nabla h_D|^2\,dx\,dt
+  { \frac{\gamma K_- }{16}}\int_{\Omega_T}   T_s(h)|\nabla(h+ h_1)|^2\,dx\,dt\\
&\quad +  \frac{{ 4}\,h_2   K_+^2}{ \gamma K_-}\int_{\Omega_T} |\nabla h_D|^2\,dx\,dt
+ \frac{\delta_1 K_-}{12}\int_{\Omega_T} |\nabla h_1|^2\,dx\,dt \\
&\quad +  \frac{3 \,h_2^2  K_+^2}{ K_-\delta_1}\int_{\Omega_T} |\nabla h_D|^2\,dx\,dt,
\end{aligned}\\
\begin{aligned}
|J_7|&\leq  \frac{\epsilon}{2}\int_{\Omega_T}  |\nabla h_1|^2\,dx\,dt
+ \frac{ \delta_1 K_-}{6}\int_{\Omega_T} |\nabla h_{1}|^2\, dx\,dt \\
&\quad + \int_{\Omega_T}\Big(\frac{\epsilon}{2}+\frac{3 \,h_2^2 K_+^2}{\delta_1 K_- }\Big)
|\nabla h_{1,D}|^2\,dx\,dt
\end{aligned} \\
|J_8| \leq \frac{\delta_1 K_-}{12}\int_{\Omega_T} |\nabla h_{1}|^2\,dx\,dt
 +\frac{6 K_+^2h_2^2}{2\,\delta_1 K_-}\int_{\Omega_T} |\nabla h_D|^2 \,dx\,dt,
\\
\begin{aligned}
| J_{9}|&\leq \frac{\delta_1 K_-}{12}\int_{\Omega_T}|\nabla h_1|^2 \,dx\,dt
+ { \frac{\gamma K_- }{16}}\int_{\Omega_T}T_s(h)|\nabla(h+ h_1)|^2 \,dx\,dt\\
&\quad +\frac{K_+^2}{K_-}\big(\frac{3\,h_2^2}{\delta_1}+ { \frac{4{h_2}}{\gamma}}\big)
\int_{\Omega_T}  |\nabla h_{1,D}|^2\,dx\,dt,
\end{aligned} \\
\begin{aligned}
| J_{10}| 
&\leq \int_{\Omega_T}   T_s(h)| Q_s(h-h_D)|\,dx\,dt
+\int_{\Omega_T} T_f(h-h_{1})|Q_f (h_1-h_{1,D})|\,dx\,dt\\
&\quad +\int_{\Omega_T} T_s(h)| Q_s(h_1-h_{1,D})|\,dx\,dt \\
&\leq \frac{3 \| Q_s\|_{L^2(0,T;H)} ^2+2 \| Q_f\|_{L^2(0,T;H)} ^2}{2  \phi} h_2^2
+\frac{\phi}{2}\int_{\Omega_T}|h-h_{D}|^2\,dx\,dt\\
&\quad +\frac{\phi}{2}\int_{\Omega_T}  |h_1-h_{1,D}|^2\,dx\,dt,
\end{aligned}\\
\begin{aligned}
|J_{11}|&\leq \frac{\phi}{2}\int_{\Omega_T} |h-h_D|^2\,dx\,dt
+\frac{\delta_1 K_-}{24}\int_{\Omega_T}|\nabla (h_1-h_{1,D})|^2\,dx\,dt\\
&\quad + \frac{\phi}{2}\| \partial_t h_D\|^2_{L^2(0,T;H)}
+\frac{12\phi^2}{\delta_1 K_-}\|\partial_t h_{1,D}\|^2_{L^2(0,T;V')}
\\
&\leq \frac{\phi}{2}\int_{\Omega_T}  |h-h_{D}|^2\,dx\,dt
+ \frac{\delta_1 K_-}{12}\int_{\Omega_T} |\nabla h_1|^2  \,dx\,dt\\
&\quad +  \frac{\delta_1 K_-}{12}\int_{\Omega_T}\! |\nabla h_{1,D}|^2\,dx\,dt
+\frac{\phi}{2}\| \partial_t h_D\|^2_{L^2(0,T;H)}
+\frac{12\phi^2}{\delta_1 K_-}\|\partial_t h_{1,D}\|^2_{L^2(0,T;V')} .
\end{aligned}
\end{gather*}
Summing up all these estimates, we obtain
\begin{equation}
\begin{aligned}
&\phi \int_\Omega (h-h_{D})^2(T,x)\, dx+\phi \int_\Omega (h_1-h_{1,D})^2(T,x)\, dx
 +\int_{\Omega_T} {\epsilon} \bigl( |\nabla h|^2  + |\nabla h_1|^2\bigr) \,dx\,dt
  \\
&+ 2\int_{\Omega_T}  \Big(  \big(K_- - \frac{1 K_+^2}{4 \gamma K_-}\big)
\big(1- \mathcal{X}_0^\epsilon(h_{1})
 L_M(\|\nabla h_{1}\|_{L^2})\big)T_s(h)\Big)%_{(1)}
|\nabla h_{1}|^2 \,dx\,dt\\
&+2\int_{\Omega_T}   \frac{\delta_1 K_-}{2} |\nabla h_{1}|^2 \,dx\,dt
 + 2 \int_{\Omega_T}
 K_- (1-\frac{9 \gamma}{8})T_s(h)|\nabla (h+h_{1})|^2 \,dx\,dt\\
&\leq  {\phi}\int_{\Omega_T} |h-h_D|^2 \,dx\,dt
 +{\phi}\int_{\Omega_T} |h_1-h_{1,D}|^2\,dx\,dt+C,
\end{aligned} \label{Chap2_41}
\end{equation}
where $C=C(u_0,v_0,h_D,h_{1,D},h_2,Q_s,Q_f)$.
We now aim at applying the Gronwall lemma in \eqref{Chap2_41}.
By the hypotheses on $K_-$ and $K_+$, 
$\big(K_- - \frac{1 K_+^2}{4 \gamma K_-}\big) \geq0 $
and taking
$\gamma$ such that $0<\gamma < 8/9$, we obtain that
\[
\big(K_- - \frac{1 K_+^2}{4 \gamma K_-}\big)\geq 0  
\]
is always true because
 $ 0 \leq 1-\mathcal{X}_0^\epsilon(h_{1})L_M(x) \leq 1$  and
$ K_+ \leq   2 \sqrt{\gamma} K_- $.


Now we apply the Gronwall lemma and we deduce that  there exists a real number
$B$, that does not depend on $\epsilon$ nor on  $M$, such that
 \begin{equation}
\begin{gathered}
\|h\|_{L^\infty(0,T;H)\cap L^2(0,T;(H^1(\Omega))')}\leq B,\quad
\|\Psi(h)\|_{ L^2(0,T;H^1(\Omega))}\leq B, \\
\|h_1\|_{L^\infty(0,T;H)\cap L^2(0,T;H^1(\Omega))}\leq B,
\end{gathered} \label{unifest}
 \end{equation}
 where we set
\[
\Psi(x)=\int_0^x(h_2-t)^{1/2}dt
=\frac{2}{3}\big(h_2^{3/2}-(h_2-x)^{3/2}\big).
\]
In particular, $\|\nabla h_1\|_{L^2(0,T;H)}\leq B$ and this estimate does not
depend on the choice of the real number $M$ that defines  function $L_M$.
Then the term $L_B\big(\|\nabla h_1\|_{L^2}\big)=1$ may be dropped.
\smallskip

\noindent\textbf{Remark on the Maximum Principle.}
 Even if we establish step 3, we can not prove that
$  h_1(t,x)+ \delta_1 \leq h(t,x)$ a.e.  $x\in  \Omega$ and $\forall t\in (0,T)$.
We set
\[
h_m=\big(\delta_1  + h_1 -h \big)^{+}\in L^2(0,T;V) \quad \text{since }
  h_{1,D}(\cdot,\cdot)+ \delta_1 \leq h_{D}(\cdot,\cdot).
\]
Likewise,  we recall that $h_m(0,\cdot) =0$ a.e. in $\Omega$ thanks to the
maximum principle satisfied by the initial data
   $h_0$ and $h_{1,0}$.
  Moreover, $ \nabla h_m=\chi_{\{h < \delta_1 + h_1\}} \nabla(h_1 - h) $.
Let $\tau\in (0,T)$, thus, taking $w(t,x)=h_{m}(x,t) \chi_{(0,\tau)}(t)$
in \eqref{Chap2_2}--\eqref{Chap2_1} yields
\begin{align*}
&\int_0^{\tau}\phi\langle\partial_t (h_1 - h),h_m\rangle_{V',V}dt
+\int_{\Omega_{\tau}} \chi_{\{h < \delta_1 + h_1\}}
 (\epsilon+T_s(h))|\nabla (h_1-h)|^2\, dx\,dt \\
&+\int_{\Omega_{\tau}} T_f(h-h_1)) \nabla h\cdot\nabla h_m \,dx\,dt\\
&+\int_{\Omega_{\tau}} \mathcal{X}_0^\epsilon(h_1)T_s(h)
 \big(\nabla h\cdot\nabla h_m - L_M\big(\|\nabla h_1\|_{L^2})
\nabla h_1\cdot\nabla h_m\big) \,dx\,dt\\
&-\int_{\Omega_{\tau}}{Q}_f T_f(h-h_1)\mathcal{X}_0^\epsilon(h_1)h_m \,dx\,dt=0.
\end{align*}
Applying Lemma 1 to \eqref{ppe3} and taking into account $h_m(0,\cdot)=0$,  gives
\[
\int_0^{\tau}\phi\langle\partial_t h_1,h_m\rangle_{V',V}dt
= \frac{\phi}{2}\int_\Omega\Big(h_m^2(\tau,x)-h_m^2(0,x)\Big)\, dx
= \frac{\phi}{2}\int_\Omega h_m^2(\tau,x)\,dx.\]
Moreover, $ L_M\big(\|\nabla h_1\|_{L^2}) =1$, then the above equation becomes
\begin{equation}
\begin{aligned}
&\frac{\phi}{2}\int_\Omega h_m^2(\tau,x)\,dx
+\int_{\Omega_{\tau}} \chi_{\{h < \delta_1 + h_1\}}
 \Big(\epsilon+T_s(h)(1-\mathcal{X}_0^\epsilon(h_1))\Big)|\nabla (h_1-h)|^2\, dx\,dt\\
&+\int_{\Omega_{\tau}} T_f(h-h_1)) \nabla h\cdot\nabla h_m \,dx\,dt
-\int_{\Omega_{\tau}}{Q}_f T_f(h-h_1)\mathcal{X}_0^\epsilon(h_1)h_m \,dx\,dt=0.
\end{aligned} \label{ppe4}
\end{equation}
Since $T_f(h-h_1) \chi_{\{h < \delta_1 + h_1\}}=\delta_1>0$ by definition of
 $T_f$,  if we assume $Q_f \leq  0$
hence \eqref{ppe3} becomes
\begin{align*}
&\frac{\phi}{2}\int_\Omega h_m^2(\tau,x)\, dx
+\int_{\Omega_{\tau}} \delta_1 \nabla h\cdot\nabla h_m \,dx\,dt\\
&\leq - \int_0^{\tau}\int_\Omega \chi_{\{h < \delta_1 + h_1\}}
 \Big(\epsilon+T_s(h)(1-\mathcal{X}_0^\epsilon(h_1))\Big)|\nabla (h_1-h)|^2
 \,dx\,dt \leq 0.
\end{align*}
and because of the term
$\int_{\Omega_{\tau}} \delta_1 \nabla h\cdot\nabla h_m \,dx\,dt $,
we can no more conclude that $h_m=0$ a.e. on  $\Omega_T$.
 \smallskip

\noindent\textbf{Step 4: Existence for the initial  system.}
We now proceed to the last step in the proof of Theorem \ref{Theo2},
namely we let  $\epsilon \to 0$.
 We infer from the above estimates that  $(h_1^\epsilon-h_{1,D})_\epsilon$
is uniformly bounded in $W(0,T)$.
 We deduce thanks to the compactness result of Aubin that
$(h_1^\epsilon-h_{1,D})_\epsilon$ is sequentially compact in $L^2(0,T;H)$.

Concerning the sequence $\{h^\epsilon-h_D\}_\epsilon$, we proceed as in
the case of the confined aquifer when $\beta =0$.
We first  observe that $\Psi$ is a strictly decreasing function of class
$\mathcal{C}^1$ on $[0,h_2]$ and $\Psi^{-1}$ is continue on
$(0, \frac{2}{3}h^{3/2}_2)$.
Using the time translate estimates (see e. g. \cite{Alt}), we deduce that
 $\{\Psi(h^\epsilon)\}_\epsilon$ is sequentially compact in $L^2(\Omega_T)$.
Namely, \eqref{unifest} yields the following time translate estimate in $L_2$
 norm of sequence $\{h^\epsilon\}_\epsilon$
\begin{align*}
&\int_{\xi}^T \langle(h^\epsilon(t,.) - h^\epsilon(t-\xi,.),
\Psi(h^\epsilon(t,x)) - \Psi(h^\epsilon(t-\xi,.)) \rangle_{V',V} \, dt
\\
& \leq  \Big( \int_{\xi}^T \| h^\epsilon(t,\cdot) - h^\epsilon(t-\xi,\cdot)
\|^{2}_{V'} \, dt \Big)^{1/2}
  \Big( \int_{\xi}^T \| \Psi(h^\epsilon(t,\cdot))
- \Psi(h^\epsilon(t-\xi,\cdot) \|^{2}_{V} \, dt \Big)^{1/2} \\
&\leq   C  \Big( \int_{\xi}^T \| h^\epsilon(t,\cdot)
- h^\epsilon(t-\xi,\cdot) \|^{2}_{V'} \, dt \Big)^{1/2}
\quad   \text{thanks to } \eqref{unifest}.
\end{align*}
However, we know that for all $\xi \in (0,T)$ (see \cite{Brezis})
\[
\frac{1}{\xi ^2} \int_{\xi}^T \| h^\epsilon(t,\cdot)
 - h^\epsilon(t-\xi,\cdot) \|^{2}_{V'}  \, dt
 \leq  \int_{\xi}^T \|  \partial _t h^\epsilon \|^2_{(H^1(\Omega))'} \, dt
 \leq  C,
 \]
 Finally we obtain
\[
\int_{\xi}^T \langle(h^\epsilon(t,.) - h^\epsilon(t-\xi,.) ,
\Psi(h^\epsilon(t,x)) - \Psi(h^\epsilon(t-\xi,.)) \rangle_{V',V} \, dt
 \leq  C  \xi, \quad \forall \xi \in (0,T).
\]
Thanks to the regularity of $\Psi$ , we deduce
\[
 \int_{\xi}^T \Big(\Psi(h^\epsilon(t,\cdot))
- \Psi(h^\epsilon(t-\xi,\cdot)) ,\Psi(h^\epsilon(t,\cdot))
- \Psi(h^\epsilon(t-\xi,\cdot)) \Big)_{L^2(\Omega)} \, dt
\leq C  \xi, 
\]
for all $\xi \in (0,T)$.
Therefore, as in [\cite{Chen}, Lemma 2.6], we deduce that
$\{\Psi(h^\epsilon)\}_\epsilon$  converges strongly in $L^2(\Omega_T)$.

 Up to the extraction of a subsequence, not relabeled for convenience,
we claim that there exists functions $h$ and $h_1$ such that
 $(h-h_D,h_1-h_{1,D})\in W(0,T)^2$ and
 \begin{gather*}
h^\epsilon\to h \quad\text{in $L^2(0,T;H)$ and a.e. in $\Omega\times (0,T)$},\\
\Psi(h^\epsilon) \rightharpoonup \Psi(h) \quad
\text{weakly  in $L^2(0,T;H^1(\Omega))$},\\
\partial_t h^\epsilon\rightharpoonup \partial_t h  \quad
\text{weakly   in $L^2(0,T;V')$},\\
h_1^\epsilon\to h_1\quad\text{in $L^2(0,T;H)$ and a.e. in $\Omega\times (0,T)$},\\
h_1^\epsilon\rightharpoonup h_1 \quad\text{weakly  in $L^2(0,T;H^1(\Omega))$},\\
\partial_t h_1^\epsilon\rightharpoonup \partial_th_1
\quad\text{weakly   in $L^2(0,T;V')$}.
\end{gather*}
Letting $\epsilon \to 0$ in  the weak formulation of the regularized problem
and using Lebesgue Theorem (thanks to the uniform estimates \ref{unifest}),
we obtain at once  \eqref{probleml}-\eqref{problem}.
The boundary and initial condition \eqref{boundary}-\eqref{initial}
holds true since the map $i \in W(0,T)  \mapsto i(0) \in H$ is continuous.
 Furthermore $(h,h_1)$ satisfies a maximum principle
which is not consistent with physical reality because we lost the information
 between $h_1$ and $h$:
\[
0\leq h_1(x,t) \quad \text{and}\quad
0 \leq h(x,t)\leq h_2, \quad  \forall t\in (0,T),  \text{ a.e. } x \in\Omega.
\]
The proof of Theorem \ref{Theo2} is complete.

\begin{thebibliography}{00}

\bibitem{Alfaro2014} M. Alfaro, P. Alifrangis;
\emph{Convergence of a mass conserving Allen-Cahn equation whose Lagrange
 multiplier is nonlocal and local}, Interfaces Free Bound., to appear.

\bibitem{Alfaro2005} M. Alfaro, D. Hilhorst, M. Hiroshi;
\emph{Optimal interface width for the Allen-Cahn equation},
RIMS Kokyuroku, 1416, 148--160, 2005.

\bibitem{Alt} H. W. Alt, S. Luckhaus;
\emph{Quasilinear elliptic-parabolic differential equations},
Math. Z., Vol. 1,  311--341, 1983.

\bibitem{BO} J. Bear, A. H. D. Cheng, S. Sorek, D. Ouazar, I. Herrera;
\emph{Seawater intrusion in coastal aquifers: Concepts, Methods and Practices},
Kluwer Academic Pub., 1999.

\bibitem{Brezis} H. Brezis;
\emph{Op\'erateurs maximaux monotones et semi-groupes de constructions
dans les espaces de Hilbert}, North-Holland, Mathematics studies, 1973.

\bibitem{BV} J. Bear, A. Verruijt;
\emph{Modeling groundwater flow and pollution, D. Reidel Publishing Company},
 Dordecht, Holland, 1987.

\bibitem{Bear} J. Bear;
\emph{Dynamics of Fluids in Porous Media}, Elsevier, 1972.

\bibitem{CahHil58} J. W. Cahn, J. E. Hilliard;
\emph{Free energy of non-uniform systems. I. Interfacial free energy},
 J. Chem. Phys., 28, 258--267, 1958.
\bibitem{Chen} Z. Chen, R. Ewing;
\emph{Mathematical analysis for reservoir models}, SIAM J. Math Anal., 30(2),
pp 431-453, 1999.

\bibitem{CDR1} C. Choquet,, M. M. Di\'edhiou, C. Rosier;
\emph{Derivation of a sharp-diffuse interface model in a free aquifer}.
Numerical simulations, submitted, 2015.

\bibitem{CDR2} C. Choquet, M. M. Di\'edhiou, C. Rosier;
\emph{Mathematical analysis of a sharp-diffuse interfaces model for seawater
intrusion}, to appear in J. Diff. Equations, 2015.

\bibitem{Dube} M. Dub\'e, M. Rost, K. R. Elder, M. Alava, S. Majaniemi,  T. Ala-Nissila;
\emph{Liquid Conservation and Nonlocal Interface Dynamics in Imbibition},
Phys. Rev. Lett., Vol. 83 (8), 1628--1631, 1999.

\bibitem{Essaid1990} H. L. Essaid;
\emph{A Multilayered Sharp Interface Model of Coupled Freshwater and Saltwater
 Flow in Coastal Systems: Model Developpement and Application.},
 Water Res. Res., Vol. 26, 1431--1454, 1990.

\bibitem{GM} G. Gagneux, M. Madaune-Tort;
\emph{Analyse math\'ematique de mod\`eles non lin\'eaires de l'ing\'enierie
 p\'etroli\`ere}. Math\'ematiques \& Applications, 22, Springer, 1996.

\bibitem{LM} J. L. Lions, E. Magenes;
\emph{Probl\`emes aux limites non homog\`enes}, Vol. 1, Dunod, 1968.

\bibitem{NR} K. Najib, C. Rosier;
 \emph{On the global existence for a degenerate elliptic-parabolic seawater
intrusion problem}, Math. Comput. Simulation, Vol. 81 Issue 1, 2282--2295, 2011.

\bibitem{RR} C. Rosier, L. Rosier;
\emph{Well-posedness of a degenerate parabolic equation issuing from
two-dimensional perfect fluid dynamics}, Applicable Anal.,
Vol. 75 (3-4), pp 441--465, 2000.

\bibitem{Simon} J. Simon;
\emph{Compact sets in the space $L^p(0,T,B)$}, Ann. Mat. Pura  Appl.,
vol. 146 (4), 65--96, 1987.

\bibitem{Tber05} M. E. Talibi, M. H. Tber;
\emph{Existence of solutions for a degenerate seawater intrusion problem},
Electron. J. differential Equ. vol. 2005, no. 72, pp 1-14, 2005.

\bibitem{Z} E. Zeidler;
\emph{Nonlinear functional analysis and its applications, Part 1,}
 Springer Verlag, 1986.

\end{thebibliography}

\end{document}





