# Sharp-diffuse interfaces model with cross-diffusion in a free aquifer: numerical analysis of a finite volume scheme

Ahmed Ait Hammou Oulhaj, Éloïse Comte, Moussa Mory Moussa Mory Diedhiou, David Maltese

► To cite this version:

Ahmed Ait Hammou Oulhaj, Éloïse Comte, Moussa Mory Moussa Mory Diedhiou, David Maltese. Sharp-diffuse interfaces model with cross-diffusion in a free aquifer: numerical analysis of a finite volume scheme. 2024. hal-04636738

HAL Id: hal-04636738

https://hal.science/hal-04636738v1

Preprint submitted on 5 Jul 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
SHARP-DIFFUSE INTERFACES MODEL WITH
CROSS-DIFFUSION IN A FREE AQUIFER : NUMERICAL
ANALYSIS OF A FINITE VOLUME SCHEME

AHMED AIT HAMMOU OULHAJ, ÊLOÏSE COMTE, MOUSSA MORY DIÉDHIOU, AND DAVID MALTESE

ABSTRACT. Sharp-diffuse interface is a new model for seawater intrusion problem which have the efficiency of the sharp interface approach and the physical realism of the diffuse interface. The introduction of the diffuse interface as a phase field in the new model makes it non-degenerate. Three-dimensional problem is reduced to a two-dimensional model involving a strongly coupled system of PDE of parabolic type describing the evolution of the depths of the interface between salt-and freshwater and the water table. We propose and analyze a finite volume scheme based on two-point flux approximation with upwind mobilities. The scheme preserves at the discrete level the main features of the continuous problem, namely the nonnegativity of the solutions, the decay of the energy and the control of the entropy and its dissipation. Based on these nonlinear stability results, we show that the scheme converges towards a weak solution to the problem. Numerical results are provided to illustrate the behavior of the model and of the scheme.

**Keywords.** Unsaturated porous media flow, seawater intrusion, nonlinear discretization, entropy stability, convergence analysis, cross-diffusion, unconfined aquifer

**AMS subjects classification.** 65M12, 65M08, 76S05

# 1. INTRODUCTION

Seawater intrusion in costal aquifers is a major problem for water supply. If they are densely populated in costal regions, this was because there existe in the ground-water, hydraulic exchanges between fresh groundwater and seawater, moreover the intensive extraction of freshwater leads for instance to local water table depression causing problems of saltwater intrusion in the aquifers. It's motivating the study of efficient and accurate models to simulate the displacement of seawater fronts for the optimal exploitation of freshwater in costal aquifer. Considerable work has been undertaken and for general information about seawater intrusion problems (see [8, 9, 10]). The static model based on Ghyben-Herzberg hypothesis do not enable to follow the evolution of the saltwater front and the existing dynamic models may be classified in four categories : The first approach is to consider that fresh and salt water are two miscible fluids. This modelling approach is physically correct and due to density contrast they tend to separate into two layers with a transition zone characterized by the variation of the salt concentration, but it has the disadvantage that it is impossible to follow explicitly the interface (see [15]). This approach is very heavy from theoretical and numerical points of view. The second consists in assuming that fresh and salt water are two immiscible fluids. That's a certain type of a sharp interface when the salty phase disappears. Nevertheless
the explicit tracking of the interface remains unworkable to implement without further assumptions (see also [15] in unsaturated media). The third approach is to assume that the fluids are immiscible and the domains occupied by each fluid are separated by an interface called sharp interface. Moreover, we consider that no mass transfer occurs between the fresh and the saltwater and the so-called Dupuit approximation. Physically, this approach is often reasonable but it enables to take into account and does not describe the behaviour of the transition zone. We refer to [8, 9, 10, 38, 37, 4] for more details about this approach. Recently, the fourth and new approach is a mixed approach (sharp-diffuse). This approach combining the simplicity of the sharp interface with the realism of diffuse interfaces. It have the advantage of respecting the physics of the problem and that of the computational efficiency. The two key assumptions of this approach are : there is no mass transfer between fresh and salt water and we suppose the existence of the diffuse interface which is modelled using a phase field approach, here an Allen-cahn type model in fluid-fluid context. The last assumption also highlighted a new term in the model which contains the thickness of the transition zone. We refer [16] for the derivation and mathematical analysis of sharp-diffuse interfaces in unconfined aquifer. In [17], there are the numerical simulations and some physical predictions from the model. The comparaison between a mixed approach and a sharp interface approach has been done in [18]. In a confined aquifer, the comparaison of a full 3D model for miscible and compressible displacements in slightly deformable porous media compared the mixed approach with numerical simulations results (see [14]).

In this present paper, we choose this new mixed (sharp-diffuse) interface approach for seawater intrusion in an unconfined aquifer. Thanks to this advantages, the sharp-diffuse model can be rewritten in the form of the formal asymptotic limit as the aspect ratio between the thickness and the horizontal length of the porous medium tends to zero. The new formulation described in [30] is a cross-diffusion system of degenerate parabolic equations for the evolution of the interface between saltwater and freshwater and the interface between the completely saturated part and the dry part in a partially saturated porous media. The mathematical analysis particularly the uniqueness of a cross-diffusion systems was done in [19]. We propose a finite volume scheme for a mixed approach based on a two-point flux approximation with upwind mobilities. We also present some numerical results.

The paper is organized as follows. In Section 2, we detail and present the continue problem, the classical entropy/dissipation property and the decrease energy functional along the time have made. The existence of nonnegative solution is shown in Section 3. Discrete counterparts of the entropy/entropy-dissipation and energy/energy-dissipation relations are established in Section 4. Section 5 is devoted to the convergence proof of the scheme. This proof is based first on the compactness of the sequence of approximate solutions and then on the identification of the limit. We finally present numerical experiments in Section 6, to illustrate the behaviour of the model and of the scheme.

## 2. PRESENTATION OF THE PROBLEM

2.1. **Continuous problem.** In our setting $\xi$ is a nonnegative function expressing the height of the interface between the saltwater and the freshwater while $h \ge \xi$ is
the height of the interface separating the freshwater and the dry soil. We assume that the bottom of the porous medium, which is located at $b \neq 0$, is impermeable, cf. Figure 1. Moreover we assume quasi-horizontal displacements (Dupuit approximation), hence we get a 2D-vertically averaged model. This assumption is reasonable when the thickness of the aquifer is small compared to the horizontal length of the aquifer. This is a sharp interface approach without sources terms.

<figure><img src="image_2.png" /><figcaption>FIGURE 1. Description of an unconfined aquifer</figcaption></figure>

The evolution of $\xi$ and $h$ is given by the following cross-diffusion system of degenerate parabolic equations

$$ (1) \qquad \left\{ \begin{array}{ll} \partial_t(h - \xi) - \nabla \cdot \left( (h - \xi)\nabla((1 - \varepsilon_0)h) \right) = 0 & \text{in } \Omega \times (0, T) =: \Omega_T \\ \partial_t\xi - \nabla \cdot \left( (\xi - b)\nabla((1 - \varepsilon_0)h + \varepsilon_0\xi) \right) = 0 & \text{in } \Omega \times (0, T), \end{array} \right. $$

where $\Omega \subset \mathbb{R}^2$ is a polygonal open bounded subset and $T > 0$ is a finite time horizon. The parameter $\varepsilon_0$ is given by

$$ \varepsilon_0 = \frac{\rho_s - \rho_f}{\rho_s}, $$

where $\rho_s$ (resp. $\rho_f$) is the mass density of the fluid saltwater (resp. freshwater) (assumed to be constant with $0 < \rho_f < \rho_s$). We set $f = h - \xi$, $g = \xi - b$ and

$$ (2) \qquad \nu = 1 - \varepsilon_0 = \frac{\rho_f}{\rho_s} \in (0, 1). $$

The system (1) then rewrites :

$$ (3) \qquad \left\{ \begin{array}{ll} \partial_t f - \nabla \cdot (\nu f \nabla(f + g + b)) = 0 & \text{in } \Omega \times (0, T), \\ \partial_t g - \nabla \cdot (g \nabla(\nu f + g + b)) = 0 & \text{in } \Omega \times (0, T). \end{array} \right. $$

When we consider the mixed approach and the bedrok $b = 0$, we have the new terms which appear, thanks to the diffuse interface these terms enable to take into account the thickness ($\delta$) of the two transition zones. We have the following cross-diffusion issuing from a mixed approach. The system is non-degenerate parabolic equations :

$$ (4) \qquad \left\{ \begin{array}{ll} \partial_t f - \nabla \cdot (\nu f \nabla(f + g)) - \nabla \cdot (\delta \nabla f) = 0 & \text{in } \Omega \times (0, T), \\ \partial_t g - \nabla \cdot (g \nabla(\nu f + g)) - \nabla \cdot (\delta \nabla g) = 0 & \text{in } \Omega \times (0, T). \end{array} \right. $$

It is supplemented with no-flux boundary conditions

$$ (5) \qquad \nabla f \cdot \mathbf{n} = \nabla g \cdot \mathbf{n} = 0, \quad \text{on } \partial\Omega \times (0, T), $$
where $\mathbf{n}$ is the unit normal to the boundary $\partial\Omega$. Initial data are

$$ (6) \qquad f|_{t=0} = f_0, \quad g|_{t=0} = g_0, $$

with $f_0, g_0 \in L^\infty(\Omega)$ and

$$ (7) \qquad f_0, g_0 \ge 0, \quad \text{a.e } x \in \Omega. $$

Before describing more precisely our results, let us mention that the problems of kind (3), have been the object of several studies. The authors in [22, 23] studied the classical solutions of the system (3) (with $b = 0$). Moreover, weak solutions are established under different assumptions in [21, 31, 32].

We recall (see [21, 31, 32]) the definition of entropy functional :

$$ \mathcal{H}(f,g) = \int_{\Omega} \left[ \Gamma(g) + \frac{1}{\nu} \Gamma(f) \right] \mathrm{d}\mathbf{x}, \quad \text{where } \Gamma(s) = s \log s - s + 1, $$

and of the energy functional :

$$ \mathcal{E}(f,g) = \int_{\Omega} \left[ \frac{\nu}{2} (f + g + b)^2 + \frac{1-\nu}{2} (g + b)^2 \right] \mathrm{d}\mathbf{x}. $$

Multiplying (formally) the first equation of (3) by $\frac{1}{\nu} \log f$ and the second equation by $\log g$, integrating over $\Omega$ and summing both equations, yields the classical entropy/dissipation property :

$$ (8) \qquad \frac{d}{dt} \mathcal{H}(f,g) + \frac{1-\nu}{2} \int_{\Omega} \left[ (\nabla f)^2 + (\nabla g)^2 \right] \mathrm{d}\mathbf{x} \le \frac{1}{2(\nu+1)} \int_{\Omega} (\nabla b)^2 \mathrm{d}\mathbf{x}. $$

Moreover multiplying (formally) the first equation of (3) by $\nu(f + g + b)$ and the second equation by $\nu f + g + b$, integrating over $\Omega$ and summing both equations, yields that the energy functional decreases along time :

$$ (9) \qquad \frac{d}{dt} \mathcal{E}(f,g) + \int_{\Omega} \left[ \nu^2 f(\nabla(f+g+b))^2 + g(\nabla(\nu f + g+b))^2 \right] \mathrm{d}\mathbf{x} = 0. $$

Let us mention that the cross-diffusion systems are extensively presented in different domain as ecology, biology, chemistry and others. In [28] the author propose and analyze a finite volume scheme for the Patlak-Keller-Segel (PKS) chemotaxis model. In [11] the authors studie the PKS model with additional cross diffusion. We refer to [6] for the analysis of a finite volume method for a cross diffusion model in population dynamics. See [36, 34] for the numerical analysis for a seawater intrusion problem in a unconfined aquifer with finite element method approximation. In [1] the authors propose a finite element method and a finite volume method and compare the results given by these two methods. In [18] the authors address the question of global existence for the sharp interface approach. For an analysis of a finite volume scheme for two-phase immiscible flow in porous media, used in petroleum engineering, we can refer to these papers [35, 27].

In this work, we propose a finite volume scheme for the problem (3). This scheme is based on a two-point flux approximation with upwind mobilities. It is designed in order to preserve at the discrete level the main features of the continuous problem : the nonnegativity of the solutions, the decay of the energy (9), the control of the entropy and its dissipation (8).
2.2. **The numerical scheme.** In this section, we explicit the discretization of the problem (3)-(5) we will study in this paper (see [2]). The time discretization relies on backward Euler scheme, while the space discretization relies on a finite volume approach (see e.g [26]), with two-point flux approximation and upstream mobility.

Let $\Omega \subset \mathbb{R}^2$ be an open, bounded, polygonal subset. An admissible mesh of $\Omega$ is given by a family $\mathcal{T}$ of a control volumes (open and convex polygons), a family $\mathcal{E}$ of edges and a family of points $(x_K)_{K \in \mathcal{T}}$ which satisfy Definition 9.1 in [26]. This definition implies that the straight line between two neighboring centers of cells $(x_K, x_L)$ is orthogonal to the edge $\sigma = K|L$.

We distinguish the interior edges $\sigma \in \mathcal{E}_{\text{int}}$ and the boundary edges $\sigma \in \mathcal{E}_{\text{ext}}$. The set of edges $\mathcal{E}$ equals the union $\mathcal{E}_{\text{int}} \cup \mathcal{E}_{\text{ext}}$. For a control volume $K \in \mathcal{T}$, we denote by $\mathcal{E}_K$ the set of its edges, by $\mathcal{E}_{K,\text{int}}$ the set of its interior edges and by $\mathcal{E}_{K,\text{ext}}$ the set of edges of $K$ included in $\partial\Omega$.

Furthermore, we denote by $d$ the distance in $\mathbb{R}^2$ and by $m$ the Lebesgue measure in $\mathbb{R}^2$ or $\mathbb{R}$. We assume that the family of meshes satisfies the following regularity requirement : there exists $\zeta > 0$ such that for all $K \in \mathcal{T}$ and all $\sigma \in \mathcal{E}_{\text{int},K}$ with $\sigma = K|L$, it holds

$$ (10) \qquad d(x_K, \sigma) \geq \zeta d(x_K, x_L) $$

For all $\sigma \in \mathcal{E}_{\text{int},K}$ with $\sigma = K|L$, we define $d_\sigma = d(x_K, x_L)$ and the transmissibility coefficient

$$ (11) \qquad \tau_\sigma = \frac{m(\sigma)}{d_\sigma}, \quad \sigma \in \mathcal{E}. $$

The size of the mesh is defined by

$$ \delta = \max_{K \in \mathcal{T}}(\text{diam}(K)). $$

In order to avoid heavier notations, we restrict our study to the case of a uniform time discretization of $(0, T)$. However, all the results presented in this paper can be extended to general time discretizations without any technical difficulty. In what follows, we assume that the spatial mesh is fixed and does not change with the time step. Let $T > 0$ be some final time and $M_T$ the number of time steps. Then the time step size and the time points are given by, respectively,

$$ \Delta t = \frac{T}{M_T}, \quad t^n = n\Delta t, \quad 0 \le n \le M_T. $$

We denote by $\mathcal{D}$ an admissible space-time discretization of $\Omega_T = \Omega \times (0, T)$ composed of an admissible mesh $\mathcal{T}$ of $\Omega$ and the values $\Delta t$ and $M_T$. The size of this space-time discretization $\mathcal{D}$ is defined by $\eta = \max(\delta, \Delta t)$.

The initial conditions are discretized by

$$ (12) \qquad f_{\mathcal{T}}^0 = \sum_{K \in \mathcal{T}} f_K^0 \mathbf{1}_K, \quad \text{where } f_K^0 = \frac{1}{m(K)} \int_K f_0(x) \, dx, \quad \forall K \in \mathcal{T}, $$

$$ (13) \qquad g_{\mathcal{T}}^0 = \sum_{K \in \mathcal{T}} g_K^0 \mathbf{1}_K, \quad \text{where } g_K^0 = \frac{1}{m(K)} \int_K g_0(x) \, dx, \quad \forall K \in \mathcal{T}, $$

and $\mathbf{1}_K$ is the characteristic function on $K$. Denoting by $f_K^n$ and $g_K^n$ approximations of the mean value of $f(., t^n)$ and $g(., t^n)$ on $K$, respectively. Taking for $b_K$ the value of $b$ in a fixed point of $K$ (for instance, the center of gravity of $K$), where $b$ is a
regular function and assume that $b_K \ge 0 \ \forall K \in \mathcal{T}$. The discretization of problem (3) is given by the following set of nonlinear equations :

$$ (14) \qquad m(K) \frac{f_K^{n+1} - f_K^n}{\Delta t} + \sum_{\sigma \in \mathcal{E}_{\text{int},K}} \tau_{\sigma} f_{\sigma}^{n+1} \nu \left( (f_K^{n+1} - f_L^{n+1}) + (g_K^{n+1} - g_L^{n+1}) + (b_K - b_L) \right) = 0, $$

and

$$ (15) \qquad m(K) \frac{g_K^{n+1} - g_K^n}{\Delta t} + \sum_{\sigma \in \mathcal{E}_{\text{int},K}} \tau_{\sigma} g_{\sigma}^{n+1} \left( \nu(f_K^{n+1} - f_L^{n+1}) + (g_K^{n+1} - g_L^{n+1}) + (b_K - b_L) \right) = 0, $$

for $K \in \mathcal{T}$ and $0 \le n \le M_T - 1$, where

$$ (16) \qquad f_{\sigma}^{n+1} = \begin{cases} f_K^{n+1} & \text{if } (f_K^{n+1} - f_L^{n+1}) + (g_K^{n+1} - g_L^{n+1}) + (b_K - b_L) \ge 0, \\ f_L^{n+1} & \text{if } (f_K^{n+1} - f_L^{n+1}) + (g_K^{n+1} - g_L^{n+1}) + (b_K - b_L) < 0. \end{cases} $$

and

$$ (17) \qquad g_{\sigma}^{n+1} = \begin{cases} g_K^{n+1} & \text{if } \nu(f_K^{n+1} - f_L^{n+1}) + (g_K^{n+1} - g_L^{n+1}) + (b_K - b_L) \ge 0, \\ g_L^{n+1} & \text{if } \nu(f_K^{n+1} - f_L^{n+1}) + (g_K^{n+1} - g_L^{n+1}) + (b_K - b_L) < 0, \end{cases} $$

where $x^+ = \max(0, x)$. We next define the numerical approximation $(f_D, g_D)$ of $(f, g)$ by

$$ f_D(x, t) = \sum_{\substack{K \in \mathcal{T} \\ 0 \le n \le M_T - 1}} f_K^{n+1} \mathbf{1}_{K \times (t^n, t^{n+1}]}(x, t), \quad \text{and} $$

$$ g_D(x, t) = \sum_{\substack{K \in \mathcal{T} \\ 0 \le n \le M_T - 1}} g_K^{n+1} \mathbf{1}_{K \times (t^n, t^{n+1}]}(x, t). $$

We also define approximations of the gradients $\nabla^D f_D$ and $\nabla^D g_D$ of $f$ and $g$, respectively. To this end, we introduce that : for $K \in \mathcal{T}$

* If $\sigma = K|L \in \mathcal{E}_{\text{int},K}$, $\mathfrak{D}_{K,L}$ is the cell ("diamond") whose vertices are given by $x_K, x_L$ and the end points of the edge $\sigma = K|L$.

* $\mathfrak{D}_{K,\sigma} = \mathfrak{D}_{K,L} \cap K$ is the cell ("triangle") whose vertices are given by $x_K$ and the end points of the edge $\sigma = K|L$.

The approximate gradient $\nabla^D S_D$ (with $S = f$, or $S = g$) is a piecewise constant function, defined in $\Omega_T$ by

$$ (18) \qquad \nabla^D S_D(x, t) = - \frac{m(\sigma)}{m(\mathfrak{D}_{K,L})} (S_K^{n+1} - S_L^{n+1}) \nu_{K,L}, \quad x \in \mathfrak{D}_{K,L}, t \in (t^n, t^{n+1}], $$

where $\nu_{K,L}$ is the unit vector normal to $\sigma$ and outward to $K$.

2.3. Main results and outline of the paper. The scheme (14)-(17) amounts to a nonlinear system to be solved at each time step. The existence of a solution to this system is therefore non trivial. The first result we highlight is thus the existence of a nonnegative solution to the scheme (14)-(17), the stability in terms of the discrete entropy and the decay of the discrete energy.
**Theorem 2.1.** There exists (at least) one solution $(f_K^{n+1}, g_K^{n+1})_{K \in \mathcal{T}, n \in \{0, \dots, M_T - 1\}}$ to the scheme (14)-(17). Moreover, $f_K^n \ge 0$, $g_K^n \ge 0$ for all $K \in \mathcal{T}$ and for all $n \in \{0, \dots, M_T\}$ and there exists $C$ depending only on $\Omega, f_0, g_0, \nu$ and $b$ such that

$$ \sup_{n \in \{0, \dots, M_T - 1\}} \mathfrak{S}^{n+1} + \sum_{n=0}^{M_T-1} \Delta t \sum_{\substack{\sigma \in \mathcal{E}_{\text{int}} \\ \sigma=K|L}} \tau_{\sigma} \left[ (f_K^{n+1} - f_L^{n+1})^2 + (g_K^{n+1} - g_L^{n+1})^2 \right] \le C(1+T), $$

and

$$ \sup_{n \in \{0, \dots, M_T - 1\}} \mathfrak{E}^{n+1} + \sum_{n=0}^{M_T-1} \Delta t \sum_{\substack{\sigma \in \mathcal{E}_{\text{int}} \\ \sigma=K|L}} \tau_{\sigma} f_{\sigma}^{n+1} \nu^2 \left( (f_K^{n+1} - f_L^{n+1}) + (g_K^{n+1} - g_L^{n+1}) + (b_K - b_L) \right)^2 + \sum_{n=0}^{M_T-1} \Delta t \sum_{\substack{\sigma \in \mathcal{E}_{\text{int}} \\ \sigma=K|L}} \tau_{\sigma} g_{\sigma}^{n+1} \left( \nu(f_K^{n+1} - f_L^{n+1}) + (g_K^{n+1} - g_L^{n+1}) + (b_K - b_L) \right)^2 \le C. $$

Our second result concerns the convergence of the scheme to a weak solution of (3)-(5). Let $(\mathcal{D}_m)_{m>0}$ be a family of admissible space-time discretization of $\Omega_T$. We denote by $(\mathcal{T}_m)_{m>0}$ the corresponding meshes of $\Omega$, with size$(\mathcal{T}_m) = \delta_m \to 0$, as $m \to 0$. We define $(f_m, g_m) := (f_{\mathcal{D}_m}, g_{\mathcal{D}_m})$ the sequence of approximate solutions constructed on the discretization $\mathcal{D}_m$. We set $\nabla^m := \nabla^{\mathcal{D}_m}$.

**Theorem 2.2.** Let $(\mathcal{D}_m)_{m>0}$ be a sequence of admissible discretizations satisfying (10) uniformly in $m$ and $\lim_{m \to +\infty} \eta_m = 0$. Let $(f_m, g_m)$ be a sequence of finite volume solutions to (14)-(17). Then there exists $(f, g)$ such that, up a subsequence,

$$ \begin{align*} f_m &\to f \quad \text{in } L^r(\Omega_T), \forall r < 4, \quad \text{and } \nabla^m f_m \to \nabla f \quad \text{weakly in } L^2(\Omega_T)^2, \\ g_m &\to g \quad \text{in } L^r(\Omega_T), \forall r < 4, \quad \text{and } \nabla^m g_m \to \nabla g \quad \text{weakly in } L^2(\Omega_T)^2, \end{align*} $$

and $(f, g) \in L^2(0, T; H^1(\Omega))^2$ is a weak solution to (3)-(5) in the following sense

$$ (19) \quad \int_0^T \int_{\Omega} (f \partial_t \psi - \nu f \nabla (f + g + b) \cdot \nabla \psi) \, dx \, dt + \int_{\Omega} f_0 \psi(\cdot, 0) \, dx = 0, $$

$$ (20) \quad \int_0^T \int_{\Omega} (g \partial_t \psi - g \nabla (\nu f + g + b) \cdot \nabla \psi) \, dx \, dt + \int_{\Omega} g_0 \psi(\cdot, 0) \, dx = 0, $$

for all test functions $\psi \in C_0^\infty(\Omega \times [0, T))$.

### 3. EXISTENCE OF A NONNEGATIVE DISCRETE SOLUTIONS

First all, we prove the positivity of the discrete solutions. This estimate allows to prove the existence of a solution to the nonlinear system (14)-(17).

**Proposition 3.1.** For all $K \in \mathcal{T}, n \ge 0$,

$$ (21) \qquad f_K^n \ge 0, \quad g_K^n \ge 0, $$
hence

$$ (22) \qquad \sum_{K \in \mathcal{T}} m(K) f_K^n = \sum_{K \in \mathcal{T}} m(K) f_K^0 = \|f_0\|_{L^1(\Omega)}, $$

and

$$ (23) \qquad \sum_{K \in \mathcal{T}} m(K) g_K^n = \sum_{K \in \mathcal{T}} m(K) g_K^0 = \|g_0\|_{L^1(\Omega)}. $$

We will prove the existence of a solution, we follow the methodology proposed in [25], using a topological degree argument [20, 33].

**Proposition 3.2.** Let $\mathcal{D}$ be an admissible discretization of $\Omega \times (0, T)$. There exists (at least) one solution to the scheme (14)-(17).

*Proof.* Let $\mu \in [0, 1]$ and define $(f_{K,\mu}^{n+1}, g_{K,\mu}^{n+1})_K$ as the solution of the scheme : $\forall K \in \mathcal{T}$

$$ m(K) \frac{f_{K,\mu}^{n+1} - f_K^n}{\Delta t} + \mu \sum_{\sigma \in \mathcal{E}_{\text{int}}, K} \tau_{\sigma} f_{\sigma,\mu}^{n+1} \nu \left( (f_{K,\mu}^{n+1} - f_{L,\mu}^{n+1}) + (g_{K,\mu}^{n+1} - g_{L,\mu}^{n+1}) + (b_K - b_L) \right) = 0, $$

$$ m(K) \frac{g_{K,\mu}^{n+1} - g_K^n}{\Delta t} + \mu \sum_{\sigma \in \mathcal{E}_{\text{int}}, K} \tau_{\sigma} g_{\sigma,\mu}^{n+1} \left( \nu(f_{K,\mu}^{n+1} - f_{L,\mu}^{n+1}) + (g_{K,\mu}^{n+1} - g_{L,\mu}^{n+1}) + (b_K - b_L) \right) = 0. $$

Reproducing the proof of Proposition 3.1, one can show that

$$ f_{K,\mu}^{n+1} \ge 0, \text{ and } g_{K,\mu}^{n+1} \ge 0 \quad \forall \mu \in [0, 1], $$

hence

$$ \sum_{K \in \mathcal{T}} m(K) f_{K,\mu}^{n+1} = \sum_{K \in \mathcal{T}} m(K) f_K^0 = \|f_0\|_{L^1(\Omega)}, $$

and

$$ \sum_{K \in \mathcal{T}} m(K) g_{K,\mu}^{n+1} = \sum_{K \in \mathcal{T}} m(K) g_K^0 = \|g_0\|_{L^1(\Omega)}. $$

Therefore, for all $K \in \mathcal{T}$, one has

$$ (24) \qquad 0 \le f_{K,\mu}^{n+1} \le \frac{\|f_0\|_{L^1(\Omega)}}{m(K)} \le \frac{\|f_0\|_{L^1(\Omega)}}{\min_{K \in \mathcal{T}} m(K)} := m_f, $$

and

$$ (25) \qquad 0 \le g_{K,\mu}^{n+1} \le \frac{\|g_0\|_{L^1(\Omega)}}{m(K)} \le \frac{\|g_0\|_{L^1(\Omega)}}{\min_{K \in \mathcal{T}} m(K)} := m_g. $$

Define the compact subset $K = [-1, m_f + 1]^{\#\mathcal{T}} \times [-1, m_g + 1]^{\#\mathcal{T}}$ of $\mathbb{R}^{\#\mathcal{T}} \times \mathbb{R}^{\#\mathcal{T}}$ and define the function $\mathcal{H}((f_K, g_K)_K, \mu) : \mathbb{R}^{\#\mathcal{T}} \times \mathbb{R}^{\#\mathcal{T}} \times [0, 1] \to \mathbb{R}^{\#\mathcal{T}} \times \mathbb{R}^{\#\mathcal{T}}$ by :
$\forall K \in \mathcal{T},$

$$
\mathcal{H}((f_K, g_K)_K, \mu) = \left( m(K) \frac{f_{K,\mu}^{n+1} - f_K^n}{\Delta t} + \mu \sum_{\sigma \in \mathcal{E}_{\text{int},K}} \tau_{\sigma} f_{\sigma,\mu}^{n+1} \nu \left( (f_{K,\mu}^{n+1} - f_{L,\mu}^{n+1}) + (g_{K,\mu}^{n+1} - g_{L,\mu}^{n+1}) + (b_K - b_L) \right), \right. \\
\left. m(K) \frac{g_{K,\mu}^{n+1} - g_K^n}{\Delta t} + \mu \sum_{\sigma \in \mathcal{E}_{\text{int},K}} \tau_{\sigma} g_{\sigma,\mu}^{n+1} \left( \nu (f_{K,\mu}^{n+1} - f_{L,\mu}^{n+1}) + (g_{K,\mu}^{n+1} - g_{L,\mu}^{n+1}) + (b_K - b_L) \right) \right).
$$

The function $\mathcal{H}$ is uniformly continuous on $\mathcal{K} \times [0, 1]$ and it follows from (24) that for all $\mu \in [0, 1]$, the nonlinear system

$$
(26) \qquad \mathcal{H}((f_K, g_K)_K, \mu) = (0, 0),
$$

cannot admit any solution on $\partial \mathcal{K}$. Therefore, the corresponding topological degree $\delta(\mathcal{H}, \mathcal{K})(\mu)$ is constant w.r.t $\mu$. For $\mu = 0$ the linear system $\mathcal{H}((f_K, g_K)_K, 0)$ admits a unique solution and the topological degree is equal to 1. Hence, the nonlinear system (26) admits at least one solution for $\mu = 1$, ensuring the existence of a solution to the scheme (14)-(17). $\square$

4. ENTROPY AND ENERGY ESTIMATES

The goal of this section is to establish discrete counterparts to the entropy/entropy-
dissipation estimate (8) and energy/energy-dissipation estimate (9). In what fol-
lows, ($f_K^n, g_K^n$) denotes a solution to the scheme (14)-(17). The proof of Theorem 2.1
is based on suitable estimates, which are shown below. This section also contains
some results that will be useful in the convergence proof of Section 5.

4.1. **Discrete** $L^2(0, T; H^1(\Omega))$ **semi-norm.** We first have to define the space $\mathcal{X}(\mathcal{D})$ the solution belongs to and the discrete $L^2(0, T; H^1(\Omega))$ semi-norm.

**Definition 4.1.** We denote by $\mathcal{X}(\mathcal{D})$ the functional space :

$$
\mathcal{X}(\mathcal{D}) = \left\{ \begin{array}{ll} u \in L^{\infty}(\Omega_T)/u & \text{is constant on } K \times (t^n, t^{n+1}] \\ \forall K \in \mathcal{T}, \forall n \in \{0, \dots, M_T - 1\} & \end{array} \right\}.
$$

**Definition 4.2.** (Discrete $L^2(0, T; H^1(\Omega))$ semi-norm) We define the discrete $L^2(0, T; H^1(\Omega))$ semi-norm on $\mathcal{X}(\mathcal{D})$ by :

$$
|u|_{1,\mathcal{D}} = \left( \sum_{n=0}^{M_T-1} \Delta t \sum_{\substack{\sigma \in \mathcal{E}_{\text{int}} \\ \sigma=K|_L}} \tau_{\sigma} (u_K^{n+1} - u_L^{n+1})^2 \right)^{1/2} .
$$

Remark 4.3. Note that

$$
(27) \qquad \| \nabla^D u \|_{L^2(\Omega_T)} = \sqrt{2} |u|_{1,\mathcal{D}}.
$$

4.2. Entropy estimate. We introduce a discrete version of entropy functional :

$$
\mathfrak{H}^n := \mathfrak{H}(f_K^n, g_K^n) = \sum_{K \in \mathcal{T}} m(K) \left( \frac{1}{\nu} \Gamma(f_K^n) + \Gamma(g_K^n) \right).
$$
**Proposition 4.4.** (Entropy stability) For all $n \in \{0, ..., M_T - 1\}$, one has

$$ (28) \qquad \begin{aligned} & \mathfrak{H}^{n+1} - \mathfrak{H}^n + \frac{1-\nu}{2} \Delta t \sum_{\substack{\sigma \in \mathcal{E}_{\text{int}} \\ \sigma=K|L}} \tau_{\sigma} \left[ (f_K^{n+1} - f_L^{n+1})^2 + (g_K^{n+1} - g_L^{n+1})^2 \right] \\ & \leq \frac{\Delta t}{2(\nu+1)} \sum_{\substack{\sigma \in \mathcal{E}_{\text{int}} \\ \sigma=K|L}} \tau_{\sigma} (b_K - b_L)^2. \end{aligned} $$

*Proof.* We multiply (14) by $\Delta t \frac{\log f_K^{n+1}}{\nu}$ and summing over $K \in \mathcal{T}$ and (15) by $\Delta t \log g_K^{n+1}$ and summing over $K \in \mathcal{T}$, provides that :

$$ A + B + C = 0, $$

where

$$ A = \sum_{K \in \mathcal{T}} m(K) \left[ \frac{1}{\nu} (f_K^{n+1} - f_K^n) \log f_K^{n+1} + (g_K^{n+1} - g_K^n) \log g_K^{n+1} \right], $$

$$ B = \Delta t \sum_{K \in \mathcal{T}} \sum_{\sigma \in \mathcal{E}_{\text{int},K}} \tau_{\sigma} f_{\sigma}^{n+1} \left( (f_K^{n+1} - f_L^{n+1}) + (g_K^{n+1} - g_L^{n+1}) + (b_K - b_L) \right) \log f_K^{n+1}, $$

$$ C = \Delta t \sum_{K \in \mathcal{T}} \sum_{\sigma \in \mathcal{E}_{\text{int},K}} \tau_{\sigma} g_{\sigma}^{n+1} \left( \nu(f_K^{n+1} - f_L^{n+1}) + (g_K^{n+1} - g_L^{n+1}) + (b_K - b_L) \right) \log g_K^{n+1}. $$

By the convexity of $\Gamma$, we find that

$$ \mathfrak{H}^{n+1} - \mathfrak{H}^n = \sum_{K \in \mathcal{T}} m(K) \left[ \frac{1}{\nu} (\Gamma(f_K^{n+1}) - \Gamma(f_K^n)) + \Gamma(g_K^{n+1}) - \Gamma(g_K^n) \right] \leq A. $$

We can rewrite $B$ and $C$ as :

$$ B = \Delta t \sum_{\substack{\sigma \in \mathcal{E}_{\text{int}} \\ \sigma=K|L}} \tau_{\sigma} f_{\sigma}^{n+1} \left( (f_K^{n+1} - f_L^{n+1}) + (g_K^{n+1} - g_L^{n+1}) + (b_K - b_L) \right) \times (\log f_K^{n+1} - \log f_L^{n+1}), $$

$$ C = \Delta t \sum_{\substack{\sigma \in \mathcal{E}_{\text{int}} \\ \sigma=K|L}} \tau_{\sigma} g_{\sigma}^{n+1} \left( \nu(f_K^{n+1} - f_L^{n+1}) + (g_K^{n+1} - g_L^{n+1}) + (b_K - b_L) \right) \times (\log g_K^{n+1} - \log g_L^{n+1}). $$

It follows from the convexity of exp that

$$ a(\log a - \log b) \geq a - b \geq b(\log a - \log b) \quad \forall a, b \in [0, +\infty[, $$

where we have used the convention $\log(0) = -\infty$ and $0 \log(0) = 0$. Hence, in view of the definition (16) and (17) of the upwind mobilities, one has

$$ B \geq \Delta t \sum_{\substack{\sigma \in \mathcal{E}_{\text{int}} \\ \sigma=K|L}} \tau_{\sigma} \left( (f_K^{n+1} - f_L^{n+1})^2 + (g_K^{n+1} - g_L^{n+1})(f_K^{n+1} - f_L^{n+1}) \right. \\ \left. \qquad + (b_K - b_L)(f_K^{n+1} - f_K^{n+1}) \right), $$
$$C \geq \Delta t \sum_{\substack{\sigma \in \mathcal{E}_{\text{int}} \\ \sigma=K|L}} \tau_{\sigma} \left( (g_{K}^{n+1} - g_{L}^{n+1})^2 + \nu(g_{K}^{n+1} - g_{L}^{n+1})(f_{K}^{n+1} - f_{L}^{n+1}) \right. \\ \left. \qquad + (b_{K} - b_{L})(g_{K}^{n+1} - g_{L}^{n+1}) \right).$$

Combining these inequalities, one deduces that

$$\mathfrak{H}^{n+1} - \mathfrak{H}^n + \Delta t \sum_{\substack{\sigma \in \mathcal{E}_{\text{int}} \\ \sigma=K|L}} \tau_{\sigma} (f_K^{n+1} - f_L^{n+1})^2 + \Delta t \sum_{\substack{\sigma \in \mathcal{E}_{\text{int}} \\ \sigma=K|L}} \tau_{\sigma} (g_K^{n+1} - g_L^{n+1})^2 \\ + (\nu+1) \Delta t \sum_{\substack{\sigma \in \mathcal{E}_{\text{int}} \\ \sigma=K|L}} \tau_{\sigma} (f_K^{n+1} - f_L^{n+1})(g_K^{n+1} - g_L^{n+1}) \leq D,$$

where

$$D = -\Delta t \sum_{\substack{\sigma \in \mathcal{E}_{\text{int}} \\ \sigma=K|L}} \tau_{\sigma} (b_{K} - b_{L}) \left[ (f_{K}^{n+1} - f_{L}^{n+1}) + (g_{K}^{n+1} - g_{L}^{n+1}) \right].$$

Using the Young inequality, one has

$$D \leq \frac{1}{2\epsilon} \Delta t \sum_{\substack{\sigma \in \mathcal{E}_{\text{int}} \\ \sigma=K|L}} \tau_{\sigma} (b_{K} - b_{L})^2 + \frac{\epsilon}{2} \Delta t \sum_{\substack{\sigma \in \mathcal{E}_{\text{int}} \\ \sigma=K|L}} \tau_{\sigma} \left[ (f_{K}^{n+1} - f_{L}^{n+1}) + (g_{K}^{n+1} - g_{L}^{n+1}) \right]^2,$$

for all $\epsilon > 0$. We choose $\epsilon = 1 + \nu$, we have

$$D \leq \frac{1}{2(\nu+1)} \Delta t \sum_{\substack{\sigma \in \mathcal{E}_{\text{int}} \\ \sigma=K|L}} \tau_{\sigma} (b_{K} - b_{L})^2 \\ + \frac{\nu+1}{2} \Delta t \sum_{\substack{\sigma \in \mathcal{E}_{\text{int}} \\ \sigma=K|L}} \tau_{\sigma} \left[ (f_{K}^{n+1} - f_{L}^{n+1})^2 + (g_{K}^{n+1} - g_{L}^{n+1})^2 \right] \\ + (\nu+1) \Delta t \sum_{\substack{\sigma \in \mathcal{E}_{\text{int}} \\ \sigma=K|L}} \tau_{\sigma} \left[ (f_{K}^{n+1} - f_{L}^{n+1})(g_{K}^{n+1} - g_{L}^{n+1}) \right].$$

Finally, one has

$$\mathfrak{H}^{n+1} - \mathfrak{H}^n + \frac{1-\nu}{2} \Delta t \sum_{\substack{\sigma \in \mathcal{E}_{\text{int}} \\ \sigma=K|L}} \tau_{\sigma} (f_K^{n+1} - f_L^{n+1})^2 \\ + \frac{1-\nu}{2} \Delta t \sum_{\substack{\sigma \in \mathcal{E}_{\text{int}} \\ \sigma=K|L}} \tau_{\sigma} (g_K^{n+1} - g_L^{n+1})^2 \leq \frac{1}{2(\nu+1)} \Delta t \sum_{\substack{\sigma \in \mathcal{E}_{\text{int}} \\ \sigma=K|L}} \tau_{\sigma} (b_K - b_L)^2,$$

where $1 - \nu > 0$ thanks to (2). $\square$

**Corollary 4.5.** There exists $C_1$ depending only on $T, \Omega, f_0, g_0, \nu$ and $b$ such that

$$(29) \qquad \mathfrak{H}^{M_T} + \frac{1-\nu}{2} \left( |f_D|^2_{1,D} + |g_D|^2_{1,D} \right) \leq C_1.$$
*Proof.* Summing (28) over $n = 0, ..., M_T - 1$ provides

$$ 
\begin{aligned}
\mathfrak{H}^{M_T} + \frac{1-\nu}{2} \sum_{n=0}^{M_T-1} \Delta t \sum_{\substack{\sigma \in \mathcal{E}_{\text{int}} \\ \sigma=K|L}} \tau_{\sigma} \left[ (f_K^{n+1} - f_L^{n+1})^2 + (g_K^{n+1} - g_L^{n+1})^2 \right] & \\
\le \mathfrak{H}^0 + \frac{1}{2(\nu+1)} \sum_{n=0}^{M_T-1} \Delta t \sum_{\substack{\sigma \in \mathcal{E}_{\text{int}} \\ \sigma=K|L}} \tau_{\sigma} (b_K - b_L)^2.
\end{aligned}
$$

As a consequence of Jensen's inequality, one has

$$ \mathfrak{H}^0 = \int_{\Omega} \left[ \frac{1}{\nu} \Gamma(f_K^0) + \Gamma(g_K^0) \right] dx \le \int_{\Omega} \left[ \frac{1}{\nu} \Gamma(f_0) + \Gamma(g_0) \right] dx < +\infty, $$

and

$$ \sum_{n=0}^{M_T-1} \Delta t \sum_{\substack{\sigma \in \mathcal{E}_{\text{int}} \\ \sigma=K|L}} \tau_{\sigma} (b_K - b_L)^2 \le T |\Omega| \|\nabla b\|_{\infty}^2, $$

concluding the proof of Corollary 4.5.

<figure><img src="image_3.png" /></figure>

We obtain immediately thanks to (27), the following discrete $L^2(\Omega_T)$ estimate on the discrete gradients :

$$ (30) \qquad \|\nabla^D f_D\|_{L^2(\Omega_T)}^2 + \|\nabla^D g_D\|_{L^2(\Omega_T)}^2 \le 2C_1. $$

**4.3. Energy estimate.** The current subsection is devoted to the proof of the discrete energy estimate. We introduce a discrete version of energy functional :

$$ \mathfrak{E}^n := \mathfrak{E}(f_K^n, g_K^n) = \sum_{K \in \mathcal{T}} m(K) \left( \frac{\nu}{2} (f_K^n + g_K^n + b_K)^2 + \frac{1-\nu}{2} (g_K^n + b_K)^2 \right). $$

**Proposition 4.6.** For all $n \in \{0, ..., M_T - 1\}$, one has

$$ (31) \qquad \mathfrak{E}^{n+1} + \Delta t \sum_{\substack{\sigma \in \mathcal{E}_{\text{int}} \\ \sigma=K|L}} \tau_{\sigma} f_{\sigma}^{n+1} \nu^2 \left( (f_K^{n+1} - f_L^{n+1}) + (g_K^{n+1} - g_L^{n+1}) + (b_K - b_L) \right)^2 \\ 
+ \Delta t \sum_{\substack{\sigma \in \mathcal{E}_{\text{int}} \\ \sigma=K|L}} \tau_{\sigma} g_{\sigma}^{n+1} \left( \nu(f_K^{n+1} - f_L^{n+1}) + (g_K^{n+1} - g_L^{n+1}) + (b_K - b_L) \right)^2 \le \mathfrak{E}^n. $$

*Proof.* We multiply (14) (resp. (15)) by $\Delta t \nu (f_K^{n+1} + g_K^{n+1} + b_K)$ (resp. $\Delta t (\nu f_K^{n+1} + g_K^{n+1} + b_K)$) and sum over $K \in \mathcal{T}$. Summing both equalities and reorganizing the sums, we get $A + B = 0$, where

$$ A = \sum_{K \in \mathcal{T}} m(K) \left[ \nu(f_K^{n+1} - f_K^n)(f_K^{n+1} + g_K^{n+1} + b_K) \right] \\ 
+ \sum_{K \in \mathcal{T}} m(K) \left[ (g_K^{n+1} - g_K^n)(\nu f_K^{n+1} + g_K^{n+1} + b_K) \right], $$
$$B = \Delta t \sum_{\substack{\sigma \in \mathcal{E}_{\text{int}} \\ \sigma = K|L}} \tau_{\sigma} f_{\sigma}^{n+1} \nu^2 \left( (f_K^{n+1} - f_L^{n+1}) + (g_K^{n+1} - g_L^{n+1}) + (b_K - b_L) \right)^2 \\ + \Delta t \sum_{\substack{\sigma \in \mathcal{E}_{\text{int}} \\ \sigma = K|L}} \tau_{\sigma} g_{\sigma}^{n+1} \left( \nu(f_K^{n+1} - f_L^{n+1}) + (g_K^{n+1} - g_L^{n+1}) + (b_K - b_L) \right)^2.$$

One has

$$A = \sum_{K \in \mathcal{T}} m(K) \left[ \nu \left( (f_K^{n+1} + g_K^{n+1} + b_K) - (f_K^n + g_K^n + b_K) \right) (f_K^{n+1} + g_K^{n+1} + b_K) \right] \\ + \sum_{K \in \mathcal{T}} m(K) \left[ (1 - \nu) \left( (g_K^{n+1} + b_K) - (g_K^n + b_K) \right) (g_K^{n+1} + b_K) \right].$$

We use the following inequality: $(a-b)a \ge \frac{1}{2}(a^2-b^2)$, $\forall a, b \in \mathbb{R}$, to get

$$A \ge \sum_{K \in \mathcal{T}} m(K) \left[ \frac{\nu}{2} \left( (f_K^{n+1} + g_K^{n+1} + b_K)^2 - (f_K^n + g_K^n + b_K)^2 \right) \right] \\ + \sum_{K \in \mathcal{T}} m(K) \left[ \frac{1-\nu}{2} \left( (g_K^{n+1} + b_K)^2 - (g_K^n + b_K)^2 \right) \right] = \mathfrak{E}^{n+1} - \mathfrak{E}^n.$$

□

**Remark 4.7.** In the discrete counterpart (31) of (9), the equality is replaced by an inequality. But as well as in the continuous setting, the function $\mathfrak{E}$ decreases along time.

**Corollary 4.8.** There exists $C_2$ depending only on $f_0, g_0, b, \nu$ and $\Omega$ such that

$$\|f_{\mathcal{D}}\|_{L^{\infty}(0,T;L^2(\Omega))} + \|g_{\mathcal{D}}\|_{L^{\infty}(0,T;L^2(\Omega))} \le C_2.$$

*Proof.* Summing (31) over $n = 0, \dots, M_T - 1$, we obtain immediately thanks to the positivity of $f_{\sigma}^{n+1}$ and $g_{\sigma}^{n+1}$ that

$$\mathfrak{E}^n \le \mathfrak{E}^0.$$

Using the Cauchy-Schwarz inequality, we obtain

$$\mathfrak{E}^0 \le \sum_{K \in \mathcal{T}} m(K) \left( \frac{\nu + \nu^2}{2} (f_K^0)^2 + 2(g_K^0)^2 + 2b_K^2 \right).$$

Moreover, for $s = f$, or $g$ one has

$$\sum_{K \in \mathcal{T}} m(K) (s_K^0)^2 \le \|s_0\|_{L^2(\Omega)}^2, \quad \text{and} \quad \sum_{K \in \mathcal{T}} m(K) b_K^2 \le |\Omega| \|b\|_{\infty}^2.$$

Hence

$$\mathfrak{E}^0 \le \|f_0\|_{L^2(\Omega)}^2 + 2\|g_0\|_{L^2(\Omega)}^2 + 2|\Omega|\|b\|_{\infty}^2 < +\infty.$$

On the other hand since $b_K, f_K$ and $g_K$ are nonnegative for all $K \in \mathcal{T}$, then we have

$$\mathfrak{E}^n \ge \sum_{K \in \mathcal{T}} m(K) \frac{\nu}{2} \left[ (f_K^n)^2 + (g_K^n)^2 \right].$$
We deduce that

$$ \sum_{K \in \mathcal{T}} m(K) \left[ (f_K^n)^2 + (g_K^n)^2 \right] \leq \frac{2}{\nu} \mathfrak{e}^n \leq \frac{2}{\nu} \mathfrak{e}^0, $$

concluding the proof of Corollary 4.8. □

Then we deduce from Proposition 4.6 that

**Corollary 4.9.** There exists $C_3$ depending only on $f_0, g_0, b, \nu$ and $\Omega$ such that

$$ \sum_{n=0}^{M_T-1} \Delta t \sum_{\substack{\sigma \in \mathcal{E}_{\text{int}} \\ \sigma=K|L}} \tau_{\sigma} f_{\sigma}^{n+1} \nu^2 \left( (f_K^{n+1} - f_L^{n+1}) + (g_K^{n+1} - g_L^{n+1}) + (b_K - b_L) \right)^2 + \\ \sum_{n=0}^{M_T-1} \Delta t \sum_{\substack{\sigma \in \mathcal{E}_{\text{int}} \\ \sigma=K|L}} \tau_{\sigma} g_{\sigma}^{n+1} \left( \nu(f_K^{n+1} - f_L^{n+1}) + (g_K^{n+1} - g_L^{n+1}) + (b_K - b_L) \right)^2 \leq C_3. $$

## 5. CONVERGENCE ANALYSIS

This section is devoted to the compactness of the approximate solution. Our goal is to show the strong compactness of the sequences $(f_m)_{m>0}$ and $(g_m)_{m>0}$ in $L^2(\Omega_T)$ and the weak compactness in $L^2(\Omega_T)^2$ of the approximate gradient of $f_m$ and $g_m$ defined by (18). As a first step, we show in §5.1 the appropriate compactness properties on the reconstructed discrete solutions. Then we identify in §5.2 the limit value (whose existence is ensured thanks to the compactness properties) as the weak solution to the problem (3).

### 5.1. Compactness properties of discrete solutions.
As it is classical for unsteady problems, we need to prove some time-compactness for the approximate solutions. We make use of the time-compactness result for degenerate parabolic equations proposed in [7], as an alternative to the classical technique that consists in estimating the time-translates (see [5] in the continuous setting and [26] in the discrete setting).

**Lemma 5.1.** There exists $C_4$ depending only on $\zeta, T, f_0, g_0, \nu$ and $b$ such that

$$ (32) \quad \sum_{n=0}^{M_T-1} \sum_{K \in \mathcal{T}} m(K)(f_K^{n+1} - f_K^n)\varphi(x_K, t_{n+1}) \leq C_4 \|\nabla\varphi\|_{L^\infty(\Omega_T)}, \quad \varphi \in \mathcal{C}_c^\infty(\Omega_T). $$

$$ (33) \quad \sum_{n=0}^{M_T-1} \sum_{K \in \mathcal{T}} m(K)(g_K^{n+1} - g_K^n)\varphi(x_K, t_{n+1}) \leq C_4 \|\nabla\varphi\|_{L^\infty(\Omega_T)}, \quad \varphi \in \mathcal{C}_c^\infty(\Omega_T). $$

*Proof.* For the sake of readability, we denote by $\varphi_K^{n+1} = \varphi(x_K, t_{n+1})$ for all $K \in \mathcal{T}$ and all $n \in \{0, ..., M_T - 1\}$. We multiply the scheme (14) by $\Delta t \varphi_K^{n+1}$ and sum for $K \in \mathcal{T}$, for $n \in \{0, ..., M_T - 1\}$. This yields :

$$ A = B, $$

where

$$ A = \sum_{n=0}^{N} \sum_{K \in \mathcal{T}} m(K)(f_K^{n+1} - f_K^n)\varphi_K^{n+1}, $$
and

$$B = -\nu \sum_{n=0}^{M_T-1} \Delta t \sum_{\substack{\sigma \in \mathcal{E}_{\text{int}} \\ \sigma=K|L}} \tau_{\sigma} f_{\sigma}^{n+1} \left[ (f_K^{n+1} - f_L^{n+1}) + (g_K^{n+1} - g_L^{n+1}) \right. \\ \left. \qquad + (b_K - b_L) \right] (\varphi_K^{n+1} - \varphi_L^{n+1}).$$

Using the Cauchy-Schwarz inequality, we get

$$|B|^2 \le \sum_{n=0}^{M_T-1} \Delta t \sum_{\substack{\sigma \in \mathcal{E}_{\text{int}} \\ \sigma=K|L}} \tau_{\sigma} f_{\sigma}^{n+1} (\varphi_K^{n+1} - \varphi_L^{n+1})^2 \\ \times \sum_{n=0}^{M_T-1} \Delta t \sum_{\substack{\sigma \in \mathcal{E}_{\text{int}} \\ \sigma=K|L}} \tau_{\sigma} f_{\sigma}^{n+1} \nu^2 \left[ (f_K^{n+1} - f_L^{n+1}) + (g_K^{n+1} - g_L^{n+1}) + (b_K - b_L) \right]^2.$$

Moreover $f_{\sigma}^{n+1} \in [\min(f_K^{n+1}, f_L^{n+1}), \max(f_K^{n+1}, f_L^{n+1})]$, hence

$$(34) \qquad 0 \le f_{\sigma}^{n+1} \le f_K^{n+1} + f_L^{n+1}, \quad \forall \sigma \in \mathcal{E}, \forall n \in \{0, \dots, M_T - 1\}.$$

Using Corollary 4.9, we get

$$|B|^2 \le C_3 \sum_{n=0}^{M_T-1} \Delta t \sum_{\substack{\sigma \in \mathcal{E}_{\text{int}} \\ \sigma=K|L}} \tau_{\sigma} (f_K^{n+1} + f_L^{n+1}) d(x_K, x_L)^2 \|\nabla \varphi\|_{L^{\infty}(\Omega_T)}^2.$$

Observe that in two space dimensions,

$$(35) \qquad \sum_{K \in \mathcal{T}} \sum_{\sigma \in \mathcal{E}_{\text{int},K}} d_{\sigma} m(\sigma) \le 2 \sum_{K \in \mathcal{T}} m(K).$$

By (11), (10) one has

$$\begin{align*} \sum_{n=0}^{M_T-1} \Delta t \sum_{\substack{\sigma \in \mathcal{E}_{\text{int}} \\ \sigma=K|L}} \tau_{\sigma} (f_K^{n+1} + f_L^{n+1}) d(x_K, x_L)^2 & \\ & \le \sum_{n=0}^{M_T-1} \Delta t \sum_{\substack{\sigma \in \mathcal{E}_{\text{int}} \\ \sigma=K|L}} m(\sigma) d(x_K, x_L) (f_K^{n+1} + f_L^{n+1}) \\ & \le \sum_{n=0}^{M_T-1} \Delta t \sum_{K \in \mathcal{T}} f_K^{n+1} \sum_{\sigma \in \mathcal{E}_{\text{int},K}} m(\sigma) d(x_K, x_L) \\ & \le \frac{2T}{\zeta} \sum_{K \in \mathcal{T}} m(K) f_K^{n+1} \le \frac{2T}{\zeta} \|f_0\|_{L^1(\Omega)}, \end{align*}$$

thanks to the mass conservation (22). Hence

$$(36) \qquad \sum_{n=0}^{M_T-1} \Delta t \sum_{\substack{\sigma \in \mathcal{E}_{\text{int}} \\ \sigma=K|L}} m(\sigma) d(x_K, x_L) f_{\sigma}^{n+1} \le \frac{2T}{\zeta} \|f_0\|_{L^1(\Omega)}.$$

This concludes the proof of (32). The proof of (33) is similar. □
We can apply the [7, Theorem 3.9], we conclude that

**Proposition 5.2.** Let $(D_m)_{m \ge 1}$ be a family of admissible space-time discretization of $\Omega_T$ such that (10) holds. Let $(f_m)_m$ and $(g_m)_m$ be the corresponding sequence of discrete solution to the scheme (14)-(17), then up to an unlabeled subsequence, there exist $f \in L^2(0, T; H^1(\Omega))$ and $g \in L^2(0, T; H^1(\Omega))$ such that

$$f_m \underset{m \to +\infty}{\longrightarrow} f, \quad \text{a.e in } \Omega_T, \quad \text{and} \quad g_m \underset{m \to +\infty}{\longrightarrow} g, \quad \text{a.e in } \Omega_T.$$

$(f_m)_{m>0}$ and $(g_m)_{m>0}$ are uniformly bounded in $L^\infty(0; T, L^2(\Omega))$ thanks to Corollary 4.8. By Corollary 4.5 and Sobolev embedding, $(f_m)_{m>0}$ and $(g_m)_{m>0}$ are bounded in $L^2(0, T; L^p(\Omega))$, with $p < +\infty$. Then, thanks to Riez-Thorin theorem, $(f_m)_{m>0}$ and $(g_m)_{m>0}$ are bounded in $L^r(\Omega_T)$ with $2 < r < 4$. Hence $(f_m)_{m>0}$ and $(g_m)_{m>0}$ are equi-integrable in $L^r(\Omega_T)$. Applying the Vitali's convergence theorem we deduce that

**Lemma 5.3.** Keeping the assumption and notations of Proposition 5.2, one has

$$f_m \underset{m \to +\infty}{\longrightarrow} f, \quad \text{strongly in } L^r(\Omega_T), \quad \text{for all } r < 4,$$

and

$$g_m \underset{m \to +\infty}{\longrightarrow} g, \quad \text{strongly in } L^r(\Omega_T), \quad \text{for all } r < 4.$$

We show now the weak compactness in $L^2(\Omega_T)^2$ of the approximate gradient of $f_m$ and $g_m$ defined in (18):

**Proposition 5.4.** Keeping the assumption and notations of Proposition 5.2, one has

$$\nabla^m f_m \underset{m \to +\infty}{\longrightarrow} \nabla f \quad \text{weakly in } L^2(\Omega_T)^2, \quad \text{and} \quad \nabla^m g_m \underset{m \to +\infty}{\longrightarrow} \nabla g \quad \text{weakly in } L^2(\Omega_T)^2.$$

*Proof.* Thanks to (30) $\nabla^D f_D$ and $\nabla^D g_D$ are bounded in $L^2(\Omega_T)^2$, then there exists a subsequence of $\nabla^D f_D$ and of $\nabla^D g_D$ (still labeled $\nabla^D f_D$ and $\nabla^D g_D$) and two function $\Theta, \Xi \in L^2(\Omega_T)^2$ such that

$$\nabla^m f_m \underset{m \to +\infty}{\longrightarrow} \Theta, \quad \text{weakly in } L^2(\Omega_T)^2,$$

and

$$\nabla^m g_m \underset{m \to +\infty}{\longrightarrow} \Xi, \quad \text{weakly in } L^2(\Omega_T)^2.$$

We refer to [13, 24] to prove that $\Theta = \nabla f$ and $\Xi = \nabla g$. $\square$

<figure><img src="image_4.png" /></figure>

## 5.2. Identification as a weak solution.

**Proposition 5.5.** Let $(f,g)$ be as in Proposition 5.2, then $f$ and $g$ are the weak solution to (3)-(5) in the sense of (19) and (20).

*Proof.* Let $\psi \in C_0^\infty(\bar{\Omega} \times [0, T))$ be a test function and $\psi_K^{n+1} = \psi(x_K, t_{n+1})$ for all $K \in \mathcal{T}$ and $n \in \{0, \dots, M_T - 1\}$. We first establish (19) from (14) and to obtain (20) from (15) is similar. In order to prove that $f$ is a weak solution, we multiply (14) by $\Delta t_m \psi_K^n$ and sum over $n \in \{0, \dots, M_T - 1\}$ and $K \in \mathcal{T}$, we obtain

$$A_m + B_m = 0,$$
where

$$A_m = \sum_{n=0}^{M_T-1} \sum_{K \in \mathcal{T}} m(K)(f_K^{n+1} - f_L^{n+1})\psi_K^n,$$

$$B_m = \nu \sum_{n=0}^{M_T-1} \Delta t_m \sum_{K \in \mathcal{T}} \psi_K^n \sum_{\sigma \in \mathcal{E}_{\text{int},K}} \tau_{\sigma} f_{\sigma}^{n+1} \left[ (f_K^{n+1} - f_L^{n+1}) + (g_K^{n+1} - g_L^{n+1}) + (b_K - b_L) \right].$$

Note that $\psi_K^{M_T} = 0$ for all $K \in \mathcal{T}$, then a discrete integration part yields

$$A_m = - \sum_{n=0}^{M_T-1} \Delta t_m \sum_{K \in \mathcal{T}} m(K) \frac{\psi_K^{n+1} - \psi_K^n}{\Delta t_m} f_K^{n+1} - \sum_{K \in \mathcal{T}} m(K) f_K^0 \psi_K^0 \\ = - \int_0^T \int_{\Omega} f_m(\delta\psi)_m \, dx \, dt - \int_{\Omega} f_m^0 \psi_m(\cdot, 0) \, dx,$$

where the function $\delta\psi_m(x_K, t) = \frac{\psi_K^{n+1} - \psi_K^n}{\Delta t_m}$, if $(x_K, t) \in K \times (t^n, t^{n+1})$. Thanks to the regularity of $\psi$, the function $\delta\psi_m$ converges uniformly towards $\partial_t\psi$ on $\Omega_T$. Moreover, we have

$$f_m \to f \quad \text{in } L^2(\Omega_T) \text{ as } m \to \infty.$$

Therefore

$$(37) \qquad \int_0^T \int_{\Omega} f_m(\delta\psi)_m \, dx \, dt \to \int_0^T \int_{\Omega} f(x) \partial_t\psi \, dx \, dt \quad \text{as } m \to \infty.$$

Moreover, $f_m^0$ converges strongly in $L^1(\Omega)$ towards the initial data $f_0$ and $\psi_m(\cdot, 0)$ converges uniformly towards $\psi(\cdot, 0)$. Therefore, we get that

$$(38) \qquad \int_{\Omega} f_m^0 \psi_m(\cdot, 0) \, dx \to \int_{\Omega} f_0(x) \psi(\cdot, 0) \, dx \quad m \to \infty.$$

We deduce from (37) and (38) that

$$A_m \to - \int_0^T \int_{\Omega} f(x) \partial_t \psi \, dx \, dt - \int_{\Omega} f_0(x) \psi(\cdot, 0) \, dx \quad m \to \infty.$$

We introduce the term

$$E_m = \int_{\Omega_T} \bar{f}_{D_m} \nabla^m u_m \cdot \nabla \psi \, dx \, dt,$$

where

$$\bar{f}_{D}(x, t) = f_{\sigma}^{n+1} \quad \forall (t, x) \in (t^n, t^{n+1}] \times \mathcal{D}_{K,L}, \quad \text{and} \quad u = f + g + b.$$

We have $f_m \xrightarrow{m \to +\infty} f$, strongly in $L^2(\Omega_T)$. Let us to prove that

$$\bar{f}_{D_m} := \bar{f}_m \xrightarrow{m \to +\infty} f, \quad \text{strongly in } L^2(\Omega_T).$$
Since $\|\bar{f}_m - f\|_{L^2(\Omega_T)} \le \|\bar{f}_m - f_m\|_{L^2(\Omega_T)} + \|f_m - f\|_{L^2(\Omega_T)}$, it is sufficient to prove that $\|\bar{f}_m - f_m\|_{L^2(\Omega_T)} \to 0$, as $\delta_m \underset{m \to +\infty}{\longrightarrow} 0$. One has

$$
\begin{align*}
\|\bar{f}_m - f_m\|_{L^2(\Omega_T)}^2 &= \sum_{n=0}^{M_T-1} \Delta t \sum_{K \in \mathcal{T}} \sum_{\sigma \in \mathcal{E}_{\text{int},K}} m(\mathfrak{D}_{K,\sigma}) (f_K^{n+1} - f_\sigma^{n+1})^2 \\
&\le \sum_{n=0}^{M_T-1} \Delta t \sum_{K \in \mathcal{T}} \sum_{\sigma \in \mathcal{E}_{\text{int},K}} m(\mathfrak{D}_{K,\sigma}) (f_K^{n+1} - f_L^{n+1})^2 \\
&\le \sum_{n=0}^{M_T-1} \Delta t \sum_{\substack{\sigma \in \mathcal{E}_{\text{int}} \\ \sigma=K|L}} m(\mathfrak{D}_{K,L}) (f_K^{n+1} - f_L^{n+1})^2 \\
&\le \frac{1}{2} \sum_{n=0}^{M_T-1} \Delta t \sum_{\substack{\sigma \in \mathcal{E}_{\text{int}} \\ \sigma=K|L}} \tau_\sigma (f_K^{n+1} - f_L^{n+1})^2 d_\sigma^2 \le \frac{C_1}{2} \delta_m^2.
\end{align*}
$$

Since $\nabla^m u_m$ converges weakly in $L^2(\Omega_T)$ to $\nabla u$, since $\bar{f}_m$ converges strongly in $L^2(\Omega_T)$ to $f$, we have

$$E_m \to \int_{\Omega_T} f \nabla u \cdot \nabla \psi \, dx \, dt \quad \text{as } m \to \infty.$$

Let us to prove $E_m - \bar{E}_m$ tends to 0 as $m \to \infty$, where

$$\bar{E}_m = \sum_{n=0}^{M_T-1} \Delta t \sum_{\substack{\sigma \in \mathcal{E}_{\text{int}} \\ \sigma=K|L}} \tau_\sigma f_\sigma^{n+1} (u_K^{n+1} - u_L^{n+1}) (\psi_K^n - \psi_L^n).$$

Using the definition of discrete gradient (18), we have

$$E_m = \sum_{n=0}^{M_T-1} \sum_{\substack{\sigma \in \mathcal{E}_{\text{int}} \\ \sigma=K|L}} f_\sigma^{n+1} \int_{t^n}^{t^{n+1}} \int_{\mathfrak{D}_{K,L}} \frac{m(\sigma)}{m(\mathfrak{D}_{K,L})} (u_K^{n+1} - u_L^{n+1}) \nabla \psi \cdot \nu_{L,K} \, dx \, dt.$$

Therefore by the definition of $\tau_\sigma$,

$$E_m - \bar{E}_m = \sum_{n=0}^{M_T-1} \sum_{\substack{\sigma \in \mathcal{E}_{\text{int}} \\ \sigma=K|L}} m(\sigma) f_\sigma^{n+1} (u_K^{n+1} - u_L^{n+1}) \int_{t^n}^{t^{n+1}} \left( \frac{\psi_K - \psi_L}{d(x_K, x_L)} - \frac{1}{m(\mathfrak{D}_{K,L})} \int_{\mathfrak{D}_{K,L}} \nabla \psi \cdot \nu_{L,K} \, dx \right) dt.$$

On the one hand, since the straight line $(x_K, x_L)$ is orthogonal to $\sigma$, we have $x_K - x_L = d(x_K, x_L)\nu_{L,K}$. It follows from the regularity of $\psi$ that

$$
\begin{align*}
\frac{\psi_K^n - \psi_L^n}{d(x_K, x_L)} &= \nabla \psi(t^n, x_L) \cdot \nu_{L,K} + O(\delta) \\
&= \nabla \psi(t, x) \cdot \nu_{L,K} + O(\eta), \quad \forall (t, x) \in (t^n, t^{n+1}) \times \mathfrak{D}_{K,L}.
\end{align*}
$$
By taking the mean value over $\mathcal{D}_{K,L}$, there exists a constant $C_5 > 0$, depending only on $\psi$, such that

$$ \left| \int_{t^n}^{t^{n+1}} \left( \frac{\psi_K^n - \psi_L^n}{d(x_K, x_L)} - \frac{1}{m(\mathcal{D}_{K,L})} \int_{\mathcal{D}_{K,L}} \nabla \psi \cdot \nu_{L,K} dx \right) dt \right| \le C_5 \Delta t \eta. $$

On the other hand, one has $m(\sigma) = \sqrt{\tau_\sigma} \sqrt{m(\sigma)d(x_K, x_L)}$. Hence by Cauchy-schwarz inequality, we have

$$ \left| \sum_{n=0}^{M_T-1} \Delta t \sum_{\substack{\sigma \in \mathcal{E}_{\text{int}} \\ \sigma=K|L}} m(\sigma) f_{\sigma}^{n+1} (u_K^{n+1} - u_L^{n+1}) \right|^2 \\ \le \sum_{n=0}^{M_T-1} \Delta t \sum_{\substack{\sigma \in \mathcal{E}_{\text{int}} \\ \sigma=K|L}} \tau_{\sigma} f_{\sigma}^{n+1} (u_K^{n+1} - u_L^{n+1})^2 \times \sum_{n=0}^{M_T-1} \Delta t \sum_{\substack{\sigma \in \mathcal{E}_{\text{int}} \\ \sigma=K|L}} m(\sigma) d(x_K, x_L) f_{\sigma}^{n+1}. $$

Using Corollary 4.9 and (36) we conclude that

$$ |E_m - \bar{E}_m|^2 \le \frac{2T}{\zeta} \|f_0\|_{L^1(\Omega)} C_2 C_5^2 \eta^2 \to 0, \text{ as } \eta \to 0. $$

This ensures that $B_m$ converges towards $\nu \int_{\Omega_T} f \nabla u \cdot \nabla \psi dx dt$, as $m \to +\infty$. $\square$

## 6. NUMERICAL RESULTS

Let us provide some illustrations of the behavior of the numerical scheme (14)-(17). The scheme leads to a nonlinear system that we solve thanks to the Newton-Raphson method. In our test case, the domain is the unit square, i.e., $\Omega = (0, 1)^2$. We consider an admissible triangular mesh made of 14336 triangles. An illustration of a mesh type used here is given in Figure 2. The numerical analysis of the scheme was carried out for a uniform time discretization of $(0, T)$ only in order to avoid heavy notations. In order to increase the robustness of the algorithm and to ensure the convergence of the Newton-Raphson iterative procedure, we used an adaptive time step procedure in the practical implementation. More precisely, we associate a maximal time step $\Delta t_{\max} = 0.00004$ for the mesh. If the Newton-Raphson method fails to converge after 30 iterations —we choose that the $\ell^\infty$ norm of the residual has to be smaller than $10^{-10}$ as stopping criterion—, the time step is divided by two. If the Newton-Raphson method converges, the time step is multiplied by two and projected on $[0, \Delta t_{\max}]$. The first time step $\Delta t$ is equal to $\Delta t_{\max}$ in the test case presented below. We perform the numerical experiments with the following data

$$ b(x, y) = \max \left( 0, \frac{1}{2} \left( 1 - 16(x - 1/2)^2 \right) \left( (\cos(\pi y) + 2) \right) \right), \quad \nu = 0.9. $$

As an initial condition we take

$$ f_0(x, y) = \begin{cases} \frac{1}{2} & \text{if } x \le \frac{1}{4}, \\ 0 & \text{elsewhere,} \end{cases} \quad g_0(x, y) = \begin{cases} b\left(\frac{1}{2}, 0\right) - b(x, y) - \left(x - \frac{1}{2}\right) & \text{if } x \le \frac{1}{2}, \\ 0 & \text{elsewhere.} \end{cases} $$

Figure 3 shows the evolution of $\xi(x, t) = b(x) + g(x, t)$ (red) and $h(x, t) = b(x) + g(x, t) + f(x, t)$ (blue) at time $t = 0, t = 0.2, t = 0.72$ and $t = 12$. There is
<figure><img src="image_5.png" /><figcaption>FIGURE 2. Mesh type used in the numerical test</figcaption></figure>

convergence towards an equilibrium state, with horizontal interfaces as expected
(see [21]).

<figure><img src="image_6.png" /><figcaption>FIGURE 3. Behaviour of the model at $t = 0, t = 0.2, t = 0.72, t = 12$</figcaption></figure>

Figure 4 shows the evolution of the energy along time

7. CONCLUSION

We proposed and analyzed a finite volume scheme for solving the seawater in-
trusion model. It preserves at the discrete level the main features of the continuous
problem : the nonnegativity of the solutions, the decay of the energy, the con-
trol of the entropy and its dissipation. Moreover, we were able to carry out a full
convergence analysis based on compactness arguments.
<figure><img src="image_7.png" /><figcaption>FIGURE 4. Evolution of the energy along time</figcaption></figure>

Let us mention that to derive the problem (3), the authors in [30] assume that, for simplification, the porous medium is isotropic. Moreover in our work the mesh satisfies the so-called orthogonality condition (see, e.g., [26, Definition 9.1]) so that the two-point flux approximation is consistent. We know that the finite volume method with two-point flux approximation, used here, does not allow to handle the anisotropic case. Nevertheless in order to treat this case, we can use for instance a Control Volume Finite Element scheme (CVFE) [29, 12, 3].

ACKNOWLEDGEMENTS

The authors thank the team Inria/Rapsodi, the Labex CEMPI (ANR-11-LABX-0007-01) and the project GEOPOR (ANR-13-JS01-0007-01) for their support. The authors also thank Claire Chainais-Hillairet and Clément Cancès for fruitful discussions. Finally, the authors thank the MIA Lab (EA 3165) and the SMAI (Société de Mathématiques Appliquées et Industrielles) for their support.

REFERENCES

[1] A. Abudawia and C. Rosier. Numerical analysis for a seawater intrusion problem in a confined aquifer. *Math. Comput. Simulation*, 118:2-16, 2015.

[2] A. Ait Hammou Oulhaj. A finite volume scheme for a seawater intrusion model with cross-diffusion. In *Finite volumes for complex applications VIII—methods and theoretical aspects*, volume 199 of Springer Proc. Math. Stat., pages 421–429. Springer, Cham, 2017.

[3] A. Ait Hammou Oulhaj, C. Cancès, and C. Chainais-Hillairet. Numerical analysis of a nonlinearly stable and positive control volume finite element scheme for Richards equation with anisotropy. HAL: hal-01372954, 2016.

[4] A. Al Bitar. *Modélisation des écoulements en milieu poreux hétérogènes 2D/3D, avec couplages surface/souterrain et densitaires*. PhD thesis, Institut National Polytechnique de Toulouse, 2007.

[5] H. W. Alt and S. Luckhaus. Quasilinear elliptic-parabolic differential equations. *Math. Z.*, 183(3):311–341, 1983.

[6] B. Andreianov, M. Bendahmane, and R. Ruiz-Baier. Analysis of a finite volume method for a cross-diffusion model in population dynamics. *Math. Models Methods Appl. Sci.*, 21(02):307–344, 2011.

[7] B. Andreianov, C. Cancès, and A. Moussa. A nonlinear time compactness result and applications to discretization of degenerate parabolic-elliptic PDEs. *J. Funct. Anal.*, 273(12):3633–3670, 2017.

[8] J. Bear. *Dynamic of builds in Porous media*. American Elsevier, New York, 1972.

[9] J. Bear. *Hydraulics of groundwater*. McGraw-Hill series in water resources and environmental engineering. McGraw-Hill International Book Co., 1979.
[10] J. Bear, A. H.-D. Cheng, S. Sorek, D. Ouazar, and I. Herrera. *Seawater intrusion in coastal aquifers: concepts, methods and practices*, volume 14. Springer Science & Business Media, 1999.

[11] M. Bessemoulin-Chatard and A. Jüngel. A finite volume scheme for a Keller-Segel model with additional cross-diffusion. *IMA J. Numer. Anal.*, 34(1):96-122, 2014.

[12] C. Cancès and C. Guichard. Convergence of a nonlinear entropy diminishing control volume finite element scheme for solving anisotropic degenerate parabolic equations. *Math. Comp.*, 85(298):549-580, 2016.

[13] C. Chainais-Hillairet and Y.-J. Peng. Finite volume approximation for degenerate drift-diffusion system in several space dimensions. *Math. Models Methods Appl. Sci.*, 14(3):461-481, 2004.

[14] L. Cherfils, C. Choquet, and M. M. Diédhiou. Numerical validation of an upscaled sharp-diffuse interface model for stratified miscible flows. *Mathematics and Computers in Simulation*, 137:246-265, July 2017.

[15] C. Choquet. Parabolic and degenerate parabolic models for pressure-driven transport problems. *Math. Models Methods Appl. Sci.*, 20(4):543-566, 2010.

[16] C. Choquet, M. M. Diédhiou, and C. Rosier. Mathematical analysis of a sharp-diffuse interfaces model for seawater intrusion. *J. Differential Equations*, 259(8):3803-3824, 2015.

[17] C. Choquet, M. M. Diédhiou, and C. Rosier. Derivation of a sharp-diffuse interfaces model for seawater intrusion in a free aquifer. numerical simulations. *SIAM J. Appl. Math.*, 76(1):138-158, 2016.

[18] C. Choquet, J. Li, and C. Rosier. Global existence for seawater intrusion models: comparison between sharp interface and sharp-diffuse interface approaches. *Electron. J. Differential Equations*, pages No. 126, 27, 2015.

[19] C. Choquet, J. Li, and C. Rosier. Uniqueness for cross-diffusion systems issuing from seawater intrusion problems. *Electron. J. Differential Equations*, 017(256):1-22, 2017.

[20] K. Deimling. *Nonlinear functional analysis*. Springer-Verlag, Berlin, 1985.

[21] J. Escher, P. Laurençot, and B.-V. Matioc. Existence and stability of weak solutions for a degenerate parabolic system modelling two-phase flows in porous media. *Ann. Inst. H. Poincaré Anal. Non Linéaire*, 28(4):583-598, 2011.

[22] J. Escher, A.-V. Matioc, and B.-V. Matioc. Modelling and analysis of the Muskat problem for thin fluid layers. *J. Math. Fluid Mech.*, 14(2):267-277, 2012.

[23] J. Escher and B.-V. Matioc. Existence and stability of solutions for a strongly coupled system modelling thin fluid films. *NoDEA Nonlinear Differential Equations Appl.*, 20(3):539-555, 2013.

[24] R. Eymard and T. Gallouët. *H*-convergence and numerical schemes for elliptic problems. *SIAM J. Numer. Anal.*, 41(2):539-562, 2003.

[25] R. Eymard, T. Gallouët, M. Ghilani, and R. Herbin. Error estimates for the approximate solutions of a nonlinear hyperbolic equation given by finite volume schemes. *IMA J. Numer. Anal.*, 18(4):563-594, 1998.

[26] R. Eymard, T. Gallouët, and R. Herbin. Finite volume methods. In *Handbook of numerical analysis, Vol. VII*, Handb. Numer. Anal., VII, pages 713-1020. North-Holland, Amsterdam, 2000.

[27] R. Eymard, R. Herbin, and A. Michel. Mathematical study of a petroleum-engineering scheme. *M2AN Math. Model. Numer. Anal.*, 37(6):937-972, 2003.

[28] F. Filbet. A finite volume scheme for the Patlak-Keller-Segel chemotaxis model. *Numer. Math.*, 104(4):457-488, 2006.

[29] P. A. Forsyth. A control volume finite element approach to NAPL groundwater contamination. *SIAM J. Sci. Statist. Comput.*, 12(5):1029-1057, 1991.

[30] M. Jazar and R. Monneau. Derivation of seawater intrusion models by formal asymptotics. *SIAM J. Appl. Math.*, 74(4):1152-1173, 2014.

[31] P. Laurençot and B.-V. Matioc. A gradient flow approach to a thin film approximation of the Muskat problem. *Calc. Var. Partial Differential Equations*, 47(1-2):319-341, 2013.

[32] P. Laurençot and B.-V. Matioc. A thin film approximation of the Muskat problem with gravity and capillary forces. *J. Math. Soc. Japan*, 66(4):1043-1071, 2014.

[33] J. Leray and J. Schauder. Topologie et équations fonctionnelles. *Ann. Sci. École Norm. Sup. (3)*, 51:45-78, 1934.
[34] P. Marion, K. Najib, and C. Rosier. Numerical simulations for a seawater intrusion problem in a free aquifer. *Appl. Numer. Math.*, 75:48–60, 2014.

[35] A. Michel. A mathematical comparison of two finite volume methods for two phase flow in porous media. In *Finite volumes for complex applications, III (Porquerolles, 2002)*, pages 213–220. Hermes Sci. Publ., Paris, 2002.

[36] K. Najib and C. Rosier. On the global existence for a degenerate elliptic–parabolic seawater intrusion problem. *Math. Comput. Simulation*, 81(10):2282–2295, 2011.

[37] M. E. A. Talibi and M. H. Tber. Existence of solutions for a degenerate seawater intrusion problem. *Electron. J. Differential Equations*, 2005(72):1–14, 2005.

[38] M. H. Tber, M. E. A. Talibi, and D. Ouazar. Parameters identification in a seawater intrusion model using adjoint sensitive method. *Math. Comput. Simulation*, 77(2):301–312, 2008.

Ahmed AIT HAMMOU OULHAJ

aithammououlhaj@gmail.com

Éloïse COMTE

Université Clermont Auvergne, INRAE, UR LISC, Centre de Clermont-Ferrand
F-63178 Aubière

eloise.comte@inrae.fr

Moussa Mory DIÉDHIOU

Université Cheikl Anta Diop de Dakar, LMA, DMI, FST
BP 5005 Dakar-Fann, Sénégal

moussamory.diedhiou@ucad.edu.sm

David MALTESE

Laboratoire d'analyse et de mathématiques appliquées, UMR 8050

Université Gustave Eiffel
5 boulevard Descartes, Bâtiment Copernic, 77420 Champs-sur-Marne

david.maltese@univ-eiffel.fr