\documentclass[reqno]{amsart}
\usepackage{hyperref}
\usepackage{graphicx}
\usepackage{mathrsfs}

\AtBeginDocument{{\noindent\small
\emph{Electronic Journal of Differential Equations},
Vol. 2012 (2012), No. 204, pp. 1--22.\newline
ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu
\newline ftp ejde.math.txstate.edu}
\thanks{\copyright 2012 Texas State University - San Marcos.}
\vspace{9mm}}

\begin{document}
\title[\hfilneg EJDE-2012/204\hfil Approximation for fluid flows]
{Numerical approximation for a degenerate parabolic-elliptic system modeling
 flows in porous media}

\author[R.-H. Bellout \hfil EJDE-2012/204\hfilneg]
{Rabah-Hacene Bellout}  % in alphabetical order

\address{Rabah-Hacene Bellout \newline
Facult\'e de Math\'ematiques, Universit\'e des Sciences\\
et Technologies Houari Boumediene, Algiers, Algeria}
\email{rbellout@usthb.dz}

\thanks{Submitted August 1, 2012. Published November 24, 2012.}
\subjclass[2000]{76S05, 35K65, 35B65, 65M06}
\keywords{Mixed finite element methods; finite volume methods; porous media}

\begin{abstract}
 We present a numerical scheme for the approximation of the system
 of partial differential equations of the Peaceman model for the
 miscible displacement of one fluid by another in a two dimensional
 porous medium. In this scheme, the velocity-pressure equations are
 treated by a mixed finite element discretization using the
 Raviart-Thomas element, and the concentration equation is
 approximated by a finite volume discretization using the Upstream scheme,
 knowing that the Raviart-Thomas element gives good approximations
 for fluids velocities and that the Upstream scheme is well suited
 for convection dominated equations. We prove a maximum principle
 for our approximate concentration more precisely
 $ 0\le c_h(x,t)\le 1$  a.e. in $\Omega_T $ as long as some
 grid conditions are satisfied - at the difference of
 Chainais and Droniou \cite{CD}who have only observed that their
 approximate concentration remains in $[0;1]$ (and such is the case
 for other proposed numerical methods; e.g., \cite {WEL,SY}).
 Moreover our grid conditions are satisfied even with very large
 time steps and spatial steps. Finally we prove the consistency
 of the proposed scheme and thus are assured of convergence.
 A numerical test is reported.
\end{abstract}

\maketitle
\numberwithin{equation}{section}
\newtheorem{theorem}{Theorem}[section]
\newtheorem{lemma}[theorem]{Lemma}
%\newtheorem{corollary}[theorem]{Corollary}
\newtheorem{definition}[theorem]{Definition}
\newtheorem{proposition}[theorem]{Proposition}
\newtheorem{remark}[theorem]{Remark}
\allowdisplaybreaks

\section{Statement of the problem} \label{intro}

 This article is concerned with the numerical approximation of the
system of partial differential equations modeling the miscible displacement
of one fluid by another in a two dimensional porous medium, when the
molecular diffusion effects are neglected.

Let $\Omega$ be a convex polygonal bounded domain in $\mathbb{R}^2$,
representing the porous reservoir, $(0,T)$ a time interval and
$\Omega_T = \Omega \times (0,T)$. Under appropriate physical assumptions,
the equations describing the displacement of one incompressible fluid by another,
completely miscible with the first, are given by:
\begin{gather}
\operatorname{div}{\mathbf{u}}=q^{+}-q^{-},  \label{intro:eqnpres} \\
{\mathbf{u}}=-\frac{1}{\mu (c)}\operatorname{grad}p, \label{intro:eqnvit} \\
\partial _{t}c+{\mathbf{u}}\cdot \operatorname{grad} c-\operatorname{div}
(D(\mathbf{u})\operatorname{grad}c)+c q^{-}={\hat{c}}
 q^{+}\quad \text{in }\Omega _{T}.  \label{intro:eqncon}
\end{gather}
We refer to Bear~\cite{Be}, Chavent and Jaffr\'{e}~\cite{CJ},
and Scheidegger \cite{Sc}
for a detailed description of the model.
Here, the gravitational terms are omitted for simplicity of exposition, $p$
denotes the pressure in the fluid mixture, $\mathbf{u}$ is the Darcy velocity,
$c$ is the concentration of one of the two components of the fluid mixture,
and $\mu$ is the concentration-dependent viscosity. By definition
\begin{equation*}
0 \le c(x,t) \le 1, \quad (x,t) \in \Omega_T.
\end{equation*}
We assume that the function $\mu$ is such that $\mu$ and $1/\mu$ are
strictly convex, and
\begin{equation}
\mu \in {C}^2([0,1]), \quad 0<\mu_{-}\le\mu(c)\le\mu_{+} \quad \forall c\in
(0,1),  \label{intro:eqnvis}
\end{equation}
where $\mu_{-}$ and $\mu_{+}$ are two fixed real numbers. The following form
is widely used, see {Koval}~\cite{Ko},
\begin{equation}
\mu(c)= \mu(0)\big(1 + (M^{1/4} -1)c\big)^{-4},  \label{intro:koval}
\end{equation}
where $M= \mu(0)/\mu(1)>1$ is the mobility ratio, this function satisfies
\eqref{intro:eqnvis}.
The stability of the flow is characterized by the mobility ratio.
At relatively high flow velocities, corresponding to P\'eclet numbers $Pe\gg 1 $,
the effects of mechanical dispersion are much greater than those of molecular diffusion,
the contribution of molecular diffusion often is negligible, see Bear~\cite{Be}, Pearson and Tardy~\cite{PT}.
In this case, the tensor $D(\mathbf{u})$ has the form:
\begin{equation}
D(\mathbf{u})=\Big(d_l  \frac{ {\mathbf{u}}{\mathbf{u}}^{\top}}{| {\mathbf{u}}|
} + d_t\Big( | {\mathbf{u}}| \mathbf{I}- \frac{ {\mathbf{u}}{\mathbf{u}}
^{\top}}{| {\mathbf{u}}|}\Big)\Big),  \label{intro:eqnbearn}
\end{equation}
where ${I}$ is the identity matrix, and $d_l\ge d_t> 0$ are the longitudinal and
transverse dispersion coefficients, respectively.

The functions $q^{+}$ and $q^{-}$ are the injection and production source
terms, respectively, and ${\hat c}$ is specified at the sources and is equal
to the resident concentration at the sinks. System~\eqref{intro:eqnpres}--
\eqref{intro:eqncon}, which is a partial differential system is of
degenerate elliptic-parabolic type, is completed by the boundary conditions
\begin{gather}
 {\mathbf{u}}\cdot \nu =0 \quad\text{ on } \Gamma_T,  \label{intro:eqncdtp} \\
 D(\mathbf{u})\operatorname{grad}c\cdot \nu =0 \quad\text{on }
\Gamma_T,  \label{intro:eqncdtc}
\end{gather}
and by the initial condition
\begin{equation}
c(x,0)=c_0(x) \quad \text{for } x\in \Omega.  \label{intro:eqncdtic}
\end{equation}
Here $\Gamma_T= \Gamma\times (0,T)$, $\Gamma$ denoting the boundary of
$\Omega$, and $\nu$ is the unit normal pointing outward $\Omega$.

Since the pressure is only determined up to a constant we additionally
require that the pressure is normalized; i.e.,
\begin{equation}
\int_{\Omega} p(x,t) dx=0, \quad t\in (0,T).  \label{intro:eqnnormaliz}
\end{equation}
Equation \eqref{intro:eqnpres} (with \eqref{intro:eqnvit}) is the pressure
equation derived from the conservation of the total mass and
\eqref{intro:eqncon} is the concentration equation derived from the
conservation of mass for one of the two components of the mixture. The
pressure equation is of elliptic nature while the concentration equation is
an advection-diffusion equation, advection being the dominant phenomenon.
Notice that, because the molecular diffusion effects are neglected, Equation
\eqref{intro:eqncon} is of degenerate parabolic type; the tensor $D(\mathbf{u}) $
is semi-positive, satisfying
\begin{equation*}
D(\mathbf{u})\boldsymbol{\xi}\cdot\boldsymbol{\xi}\ge {d_t}
|{\mathbf{u}} |\,|  \boldsymbol{\xi}|^2, 
\quad | D(\mathbf{u})\boldsymbol{\xi}| \le {d_l} |{
\mathbf{u}} | |\boldsymbol{\xi}| \quad\text{for }
\boldsymbol{\xi}\in \mathbb{R}^2. \label{intro:tensor2}
\end{equation*}
The diffusion operator may be zero pointwise, it can be small or zero in
regions of the solution space, and fairly large for other values of the
solution. Consequently, solutions of the concentration equation will, in
general, possess minimal smoothness. We do not know \textit{a priori} the
regions where $u$ has zero, small or large values. Taking into account the molecular
diffusion effects makes the tensor $D(\mathbf{u}) $ more regular and our next
analysis will be still valid. We did not also include the permeability in our model
since we were mainly concerned by the degenerate type of the concentration equation \eqref{intro:eqncon}.
We assume that $q^{+}$ and $q^{-}$ are non-negative functions satisfying $q^{+}, \; q^{-} \in L^\infty(0,T; L^2(\Omega))$  with
\begin{equation}
\int_\Omega ( q^{+}- q^{-})(x,t)dx =0 \quad\text{a.e. for } t\in (0,T)
\label{intro:eqnsourcep}
\end{equation}
while $c_{0}$ and ${\hat c}$ satisfy
\begin{equation}
{\hat c},  c_0 \in L^\infty (\Omega),\quad
0\le {\hat c}(x)\le 1, \quad
0\le c_0(x) \le 1\text{ a.e. in } \Omega.  \label{intro:eqnsourcec}
\end{equation}
System \eqref{intro:eqnpres}--\eqref{intro:eqncon}, \eqref{intro:eqncdtp}--
\eqref{intro:eqncdtic} is referred in the following as Problem
${\mathscr P}$. This problem is a nonlinear system coupling an elliptic
equation and a degenerate parabolic equation. A weak solution of Problem
${\mathscr P}$ is defined by the following sense.

\begin{definition}\label{intro:weaksol} \rm
A pair $(p,c)$ is said to be a weak solution of Problem ${\mathscr P}$ if:
\begin{itemize}
\item[(i)] $p\in L^\infty(0,T; H^1(\Omega))$ and $p$ is a solution of the
elliptic problem \eqref{intro:eqnpres}, \eqref{intro:eqnvit},
\eqref{intro:eqncdtp}, \eqref{intro:eqnnormaliz};

\item[(ii)] $c\in L^\infty(\Omega_T)$ with $0\le c(x,t)\le 1$ for a.
e. $(x,t)\in \Omega_T$, $| {\mathbf{u}}|^{1/2}\operatorname{grad}
 c \in (L^2(\Omega_T))^2$, and $c$ satisfies \eqref{intro:eqncon},
\eqref{intro:eqncdtc}, \eqref{intro:eqncdtic} in the sense
\begin{equation} \label{intro:eqnconfaible}
\begin{aligned}
& \int_{\Omega_{T}}\{c \partial_{t}\varphi+{c {\mathbf{u}}}\cdot
\operatorname{grad}\varphi-D(\mathbf{u})\operatorname{grad}
c\cdot\operatorname{grad}\varphi-c q^{-}\varphi
\}\,dx\,dt\\
& =-\int_{\Omega_{T}}{\hat{c}} q^{+}\varphi\,dx\,dt-\int_{\Omega}
c_{0}(x)\varphi(x,0)dx,
\end{aligned}
\end{equation}
for any $\varphi$ in $C^2(\overline{\Omega}_T)$ with compact support in
$\overline{\Omega} \times [0,T[$.
\end{itemize}
\end{definition}

The existence of a weak solution of Problem ${\mathscr P}$ is proved by
Amirat and Ziani in~\cite{AZ}.
In this paper we are concerned with the numerical approximation of a weak
solution of Problem ${\mathscr P}$. Several numerical procedures for
hyperbolic or advection-diffusion equations arising from flow problems
in porous media have been proposed and analyzed. We refer to
Douglas and Hornung \cite{DH}, Wheeler \cite{Wh} and references therein.

The following works are closely related to our subject. Jaffr\'e and Roberts~
\cite{JR} considered the non degenerate elliptic-parabolic system modeling a
miscible flow, including molecular diffusion effects. They derived error
estimates for a semi-discretized problem.
Eymard and Gallou\"et~\cite{EG} considered a system of an elliptic equation
and a hyperbolic equation arising in a flow problem in porous media. They
proved the convergence of a finite element - weighted finite volume scheme.
Chainais and Droniou \cite{CD} proved the convergence of a finite volume scheme
for an elliptic-parabolic system. Ameziane and El Ossmani \cite{AO} proposed a mixed finite
element - finite volume method based on a better Control Volume than we did but
using many more assumptions, among them their assumption $ (A10) $ which they say is determined by the
geometry of the grid and by the "character" of the anisotropy for the diffusion-dispersion tensor $D$.
We also mention the papers dealing with finite volume schemes for flow in porous media
by Eymard et al \cite{EHV}, Sun and Yuan \cite{SY}, Marpeau and Saad \cite {MS}.
In this last reference, a particular emphasis is made on the control of CPU time and we
are instead mainly concerned with the degenerate nature of the concentration equation.
More recently, Choquet and Zimmermann \cite {CZ}, proved a maximum principle for the concentration
after a lengthy proof based on projection theorems at the difference of our simple direct proof.

Here, we present a numerical scheme in which the velocity-pressure equation
is treated by a space mixed finite element discretization and the
concentration equation is approximated by a space finite volume
discretization with a time semi-implicit Euler scheme. Our main contribution
concern the proof of uniform estimates of our approximate concentration, more precisely $ 0\le c_h(x,t)\le 1 \; \text{ a.e. in } \Omega_T $ under suitable CFL conditions \eqref{sec2:eqhydt} and \eqref{sec2:hyptjk2}.

The outline of the paper is as follows. In Sect.~\ref{sec2} we present the
discrete equations and state the main theorem concerning uniform estimates
of the approximate solutions and the convergence of the proposed scheme.
Sect.~\ref{sec3} is devoted to the proof of the theorem. In
Sect. \ref{sec4} we report some numerical results we obtained in the case of
a benchmark incompressible flow within a horizontal reservoir over a period
of ten years with injection and production wells at the corners studied by
R.E.\ Ewing and al \cite{WEL}.

\section{The discrete problem }
\label{sec2}

\subsection{Notation and Definitions}
\label{sec2:1} Let $h>0$ be a small parameter and let ${\mathscr  T}_h$ be
an admissible mesh of $\Omega$, in the sense of
\cite[Definition 3.5]{EGH}, by triangles with diameter bounded by $h$.
Let $S_h$ be the set
of triangles of ${\mathscr T}_h$. For each triangle $T_j\in S_h$,
 $m_j$ denotes the area of $T_j$ and $x_j$ stands for the orthocenter of
$T_j$.
For two adjacent triangles $T_j$ and $T_k$, we denote
 $\sigma_{jk}= T_j\cap T_k$,
$\tau_{jk}=m_{jk}/ d_{jk}$ where $m_{jk}=m(\sigma_{jk})$ is the length of
$\sigma_{jk}$ and $d_{jk}$ is the distance between $x_j$ and $x_k$, and
 $\nu_{jk}$ is the outward unit normal on $\sigma_{jk}$ in the direction $T_k$.

We denote by ${\mathscr E}_h^+$ the set of all edges of the triangles of ${
\mathscr T}_h$, ${\mathscr E}_h^0$ the set of the edges $\sigma$ located
on the boundary $\Gamma$, and
${\mathscr E}_h ={\mathscr E}_h^+\setminus{\mathscr E}_h^0$.
Any edge $\sigma$ in ${\mathscr E}_h$ will be denoted $\sigma_{jk}$,
with the convention $j < k$, since it is a common edge to two
triangles $T_j$ and $T_k$. For convenience we denote by $j$
($j\in S_h$) the triangle $T_j$ of ${\mathscr  T}_h$.
For $j\in S_h$, $N_j$ is the set of triangles $T_k$ of ${\mathscr T}_h$
that share an edge $\sigma_{jk}$, $\partial T_j \cap \partial T_k \not= \emptyset$.
We will use the following property concerning the triangulation:
$ \sum_{(j,k)\in {\mathscr E}_h} m_{jk} \le C  h^{-1},$
where $C$ is a number which depends only on the diameter of $\Omega$ and on
the regularity of the triangulation.

Let the spaces $V$ and $W$ be defined as
\begin{equation*}
W=L^2(\Omega), \quad V=\{ \mathbf{v}\in L^2(\Omega)^2; \; \operatorname{div}
 \mathbf{v}\in L^2(\Omega), \; \mathbf{v}\cdot \nu \mid_{\Gamma}=0 \},
\end{equation*}
The space $V$ being equipped with the norm
$ \| \mathbf{v}\|_{V}^2 = \| \mathbf{v}\|_{L^2(\Omega)}^2
+ \| \operatorname{div} \mathbf{v}\|_{L^2(\Omega)}^2$.
We introduce two discrete subspaces of $V$ and $W$ due to Raviart and
 Thomas~\cite{RT}:
\begin{gather*}
V_h = \{ \mathbf{v}_h; \; \mathbf{v}_h\in V,\; \mathbf{v}_h\mid_{T_j}=
a_j+ c_j   x   \text{ for any } j \in S_h \}, \\
W_h = \{ w_h \in W\; ; \; w_h\mid_{T_j} \text{ is constant for any
} j \in S_h \},
\end{gather*}
where $a_j=({a_j}_1, {a_j}_{2})\in \mathbb{R}^2$,
 $c_j\in \mathbb{R}$, and $x=(x_1, x_{2})$ is the space variable.
We also define the spaces $ W^{0}= W/{\mathbb{R}} $ ,
and  $ \quad W_h^{0}= W_h/{\mathbb{R}} $.

The space $W_h$ is generated by the basis functions
$(\psi_j)_{j\in S_h}$ that are piecewise constant,
$\psi_j(x)= 1 $ for $  x\in T_j $  and zero elsewhere.
For the flux function $ \mathbf{v}_h$ in  $ V_h$, we have
$ \mathbf{v}_h=\sum_{\sigma\in {\mathscr E}_h}
v_{\sigma} \boldsymbol{\phi}_{\sigma}, \quad v_\sigma \in \mathbb{R}$.
Note that  ${\sigma} \in {\mathscr E}_h$   stands for an edge
$\sigma_{jk}$ of triangles $T_j$ and $ T_k$ associated with
the normal unit vector ${\nu}_{jk}$. The basis function
$ \boldsymbol{\phi}_{\sigma}$ is non zero only on $T_j\cup T_k$, the
restriction of  $\boldsymbol{\phi}_{\sigma}$ on $T_j$ or $T_k$
writes $(a_1 + c x_1, a_{2} + c x_2)$, and $ \boldsymbol{\phi}_{\sigma}\cdot \nu_{\sigma'} = \delta_{\sigma,
\sigma'}$, where $\delta_{\sigma, \sigma'}$ is the
Kronecker delta symbol. These basis functions can be found
in Raviart and Thomas~\cite{RT} or in Ern and Guermond~\cite[p.41-42]{EGu}.

Let $N$ be a positive integer. We set as time step $\Delta t=T/N$ and
$t_n=n\Delta t$ for $n=0, N$. The discrete unknowns are the values of the
pressure and the concentration $(p_j^n,c_j^n )$, and that of the flux
$u_\sigma^n$, for $1 \le n \le N$,  $j \in S_h$,
$\sigma \in {\mathscr E}_h$, given by the numerical procedure described below,
see Subsect.~\ref{sec2:2}. We use the following approximation for the initial
 data $c_0$:
\begin{equation*}
 c_j^0= {\frac{1} {m_j}} \int_{T_j} c_0(x)dx\quad
  \text{for any $j \in S_h$}.
\end{equation*}
For $0 \le n \le N$, we define $p^n_h$ and $c^n_h$ in $W_h$ and
${\mathbf{u}}^n_h \in V_h$ by
\begin{equation*}
p^n_h= \sum_{j\in S_h} p^n_j  \psi_j,\quad
c^n_h= \sum_{j\in S_h} c^n_j  \psi_j,\quad
{\mathbf{u}}^n_h=\sum_{\sigma\in {\mathscr E }_h} u^n_{\sigma}
 \boldsymbol{\phi}_{\sigma}.
\end{equation*}
We also define, for $x\in T_j$, $j\in S_h$, $t_{n} \le t <t_{n+1}$, and
$0 \le n \le N-1$,
\begin{gather}
 p_h(x,t)= p^n_h(x) + {\frac{t-t_{n}}{{\Delta t}}}
  (p^{n+1}_h(x)- p^n_h(x)), \quad
   {\mathbf{u}}_h(x,t)= {\mathbf{u}}^n_h(x), \label{defph}\\
 c_h(x,t)= c^n_h(x) + {\frac{t-t_{n}}{{\Delta t}}}
(c^{n+1}_h(x)-c^n_h(x))\label{defch}.
\end{gather}
 For $f= q^{+}$, $q^{-}$, or ${\hat c}$, we denote
\[
 f_{n,j} ={\frac{1} {m_j  {\Delta t}}} 
\int_{t_n}^{t_{n+1}} \int_{T_j} f(x,s)dx\,ds,\;  f_{n,h}(x)
= f_{n,j}, \; f_h(x,t) = f_{n,h} (x),
\]
for $x\in T_j$, $j\in S_h$, $t_{n} \le t < t_{n+1}$, and $0 \le n \le
N-1$.

\subsection{The discrete equations}
\label{sec2:2} To discretize the pressure equation, we use a mixed finite
element scheme for the spatial discretization. Suppose that the approximate
solution $(p_h^n,{\mathbf{u}}_h^n,c_h^n )$ of $(p,{\mathbf{u}},c)$ at time $
t_n $ is given, $(p,c)$ denoting a weak solution of Problem $(P)$. Then $
(p_h^{n+1},{\mathbf{u}}_h^{n+1}) $ is determined by
\begin{gather}
p_h^{n+1}\in W_h^{0}, \quad {\mathbf{u}}^{n+1}_h\in V_h,  \\
\int_{\Omega} \operatorname{div} {\mathbf{u}}^{n+1}_h(x)
\psi_h(x) \,dx = \int_{\Omega}( q^{+}_{n,h} - q^{-}_{n,h})(x)
 \psi_h(x) \,dx \quad \forall \psi_h\in W_h \label{sec2:eq1}\\
\int_{\Omega} \mu(c^n_h(x)){\mathbf{u}}^{n+1}_h(x)\mathbf{v}
_h(x)dx-\hspace*{-0.2cm}\int_{\Omega} p^{n+1}_h(x)\operatorname{div}
( \mathbf{v}_h(x))dx =0\quad \forall\mathbf{v}_h\in V_h \label{sec2:eq2}
\end{gather}
The existence and uniqueness of a solution of Problem \eqref{sec2:eq1},
\eqref{sec2:eq2} is well-known, see Raviart and Thomas~\cite{RT}, and Brezzi
and Fortin~\cite{BF}, since the space $( W_h, V_h)$ satisfy the inf-sup
condition; there is $\gamma>0$ such that
\begin{equation*}
 \inf_{\psi\in W_h}  \sup_{\mathbf{v}\in V_h} \frac{(
\operatorname{div} \mathbf{v}, \psi)}{\| \mathbf{v}\|_{H(div)} 
\| \psi\|_{L^2}} \ge \gamma.
\end{equation*}
The pair $(p_h^{n+1},{\mathbf{u}}^{n+1}_h) \in W_h^{0}\times V_h$
is then uniquely defined.
Denoting
\begin{equation*}
{u}^{n+1}_{j,k} ={\frac{1} {m_{jk}}}  \int_{\sigma_{jk}} {\mathbf{u}}
^{n+1}_h(x)\cdot \nu_{jk}\, d\gamma(x), \quad j\in S_h, \; k \in N_j,
\end{equation*}
System \eqref{sec2:eq1}, \eqref{sec2:eq2} can be written in the form:
\begin{gather}
\sum_{k\in N_j} m_{jk} u^{n+1}_{j,k} = m_j( q^{+}_{n,j} -
q^{-}_{n,j}) \label{sec2:eq3} \\
 \sum_{\sigma'\in {\mathscr E}_h} {u}
^{n+1}_{\sigma'}\int_{\Omega}\mu(c^n_h(x)) \boldsymbol{\phi}
_{\sigma'}(x) \boldsymbol{\phi}_{\sigma}(x)dx-\hspace*{-0.2cm}
\sum_{k\in S_h}p^{n+1}_k\sum_{i\in N_k}\int_{\sigma_{ki}}
\boldsymbol{\phi}_{\sigma}(x)\nu_{ki}d\gamma=0 \label{sec2:eq4}
\end{gather}
for any $j\in S_h$ and $\sigma\in {\mathscr E}_h$, with the constraint
\begin{equation}
\sum_{k\in S_h} m_k  p^{n+1}_k = 0  \label{sec2:eqnpmn}
\end{equation}
We now turn to the the discretization of the concentration
equation \eqref{intro:eqncon}. First, we write \eqref{intro:eqncon} in the
conservative form
\begin{equation}
\partial_t c + \operatorname{div} (c {\mathbf{u}} - D(\mathbf{u})
\operatorname{grad} c) + {c} q^{-} ={\hat c} q^{+ } \text{ for }
(x,t)\in \Omega_T.  \label{sec2:eqnnc}
\end{equation}
We use a semi-implicit Euler scheme in time and an upwind finite volume
scheme in space. We write formally, for each $j\in S_h$,
\begin{align*}
&\int_{T_j}(c(x,t_{n+1}) - c(x,t_{n})) \,dx + {\Delta t} \int_{T_j}
\operatorname{div} (c(x,t_{n}){\mathbf{u}}_h^{n+1}(x))dx \\
&- {\Delta t} \int_{T_j} \operatorname{div}(D({\mathbf{u}}
_h^{n+1}(x))\operatorname{grad}c(x,t_{n+1}))) \,dx + {\Delta t}
 \int_{T_j} q_{n,h}^{-}(x)c(x,t_{n})dx \\
& = {\Delta t} \int_{T_j }  {\hat c}_{n,h}(x)   q_{n,h}^{+}(x) \,dx.
\end{align*}
We approximate the convective term by
\begin{align*}
\int_{T_j}\operatorname{div} (c(x,t_{n}){\mathbf{u}}
_h^{n+1}(x))dx
&= \int_{\partial T_j} c(x,t_{n}){\mathbf{u}}
_h^{n+1}(x)\cdot\nu_j \,d\gamma(x) \\
&\equiv   \sum_{k\in N_j} Q_{jk}( {\mathbf{u}}^{n+1}_h,
c_j^n, c_k^n).
\end{align*}
Here $Q_{jk}$ is the numerical flux constructed with the upwind method:
\begin{equation}
Q_{jk}( {\mathbf{u}}^{n+1}_h, c_j^n, c_k^n) = m_{jk}
(u^{n+1}_{j,k})^{+} c_j^n + m_{jk} (u^{n+1}_{j,k})^{-} c_k^n,
\label{sec2:eq5flu}
\end{equation}
where
\begin{equation*}
(u^{n+1}_{j,k})^{+} = \begin{cases}
u^{n+1}_{j,k} &\text{if } u^{n+1}_{j,k} >0, \\
0 & \text{otherwise,}
\end{cases}
 \quad
  (u^{n+1}_{j,k})^{-} = \begin{cases}
u^{n+1}_{j,k} & \text{if } u^{n+1}_{j,k} < 0, \\
0 & \text{otherwise. }
\end{cases}
\end{equation*}
For the dispersive term, we write
\[
\int_{T_j}\operatorname{div}(D({\mathbf{u}}_h^{n+1}(x))
\operatorname{grad}c(x,t_{n+1})) \,dx
= \int_{\partial T_j}  D({\mathbf{u}}
_h^{n+1}(x))\operatorname{grad}c(x,t_{n+1}) \cdot\nu_j
\,d\gamma(x)
\]
Taking into account the form \eqref{intro:eqnbearn} of the tensor $D(\mathbf{u})$,
we have to handle with two types of terms (which give the weighted flux
through the edges $\sigma_{jk}$)
\begin{equation*}
 {\mathbf{u}}^{n+1}_h(x) \cdot\operatorname{grad} c(x,t_{n+1}) \quad
 \text{and}\quad {| {\mathbf{u}}^{n+1}_h(x)|} \operatorname{grad}
  c(x,t_{n+1})\cdot\nu_j .
\end{equation*}
 We have for the first one,
\begin{align*}
&\int_{\partial T_j} \frac{ {\mathbf{u}}^{n+1}_h(x) ({\mathbf{u}}
^{n+1}_h(x))^{\top}}{| {\mathbf{u}}^{n+1}_h(x)|} \operatorname{grad}
 c(x,t_{n+1}) \cdot\nu_j \,d\gamma(x) \\
&= \int_{\partial T_j}  \frac {{\mathbf{u}}^{n+1}_h(x)}{| {\mathbf{u}}
^{n+1}_h(x)|} \cdot\operatorname{grad} c(x,t_{n+1})\; {\mathbf{u}}
^{n+1}_h(x)\cdot\nu_j\,d\gamma(x) \\
&\equiv  \sum_{k\in N_j}   \frac{(c_k^{n+1}-c_j^{n+1})
}{d_{jk}}  \int_{\sigma_{jk}} \frac{|{\mathbf{u}}
^{n+1}_h(x)\cdot\nu_{jk}|^2} {| {\mathbf{u}}^{n+1}_h(x)|}
\,d\gamma(x).
\end{align*}
We have neglected the tangential component of the vector
$ \operatorname{grad} c(x,t_{n+1}) $
since we are integrating over the control volume $T_j$ and approximating
 $ c(x,t_{n+1})$ by  $c^{n+1}_h$ in $W_h$. The integration over a better
 control volume will be more precise certainly but once more we are mainly
concerned by the degenerate type of the concentration equation \eqref{intro:eqncon}.
The second term is given by
\begin{align*}
&\int_{\partial T_j}| {\mathbf{u}}^{n+1}_h(x)|\operatorname{grad}
 c(x,t_{n+1})\cdot\nu_j d\gamma(x)\\
&\equiv \sum_{k\in N_j} \frac{(c_k^{n+1}-c_j^{n+1})}{d_{jk}}  
\int_{\sigma_{jk}}  | {\mathbf{u}}^{n+1}_h(x)|\, d\gamma(x).
\end{align*}
Then
\[
\int_{T_j}\operatorname{div}(D({\mathbf{u}}_h^{n+1}(x))
\operatorname{grad}c(x,t_{n+1})) \,dx
\equiv \sum_{k\in N_j} {D }_{jk}( {\mathbf{u}}^{n+1}_h, c_j^{n+1}, c_k^{n+1}),
\]
 where ${D}_{jk}$ is given by
\begin{equation}
 {D}_{jk}({\mathbf{u}}^{n+1}_h, c_j^{n+1}, c_k^{n+1}) = f_{jk}^{n+1}\;
\frac{(c_k^{n+1}-c_j^{n+1})}{d_{jk}}  \label{sec2:eq6disp}
\end{equation}
with
\begin{align*}
f_{jk}^{n+1}
& = (d_l-d_t)    \int_{\sigma_{jk}} \frac{|{
\mathbf{u}}^{n+1}_h(x)\cdot\nu_{jk}|^2} {| {\mathbf{u}}^{n+1}_h(x)|}
\,d\gamma(x) + d_t  \int_{\sigma_{jk}}  | {\mathbf{u}}^{n+1}_h(x)|\,
d\gamma(x).   \\
&= \int_{\sigma_{jk}} D({\mathbf{u}}_h^{n+1}(x))\nu_{jk}\cdot\nu_{jk}
\,d\gamma(x)
\end{align*}
Then we define the approximate solution $c^{n+1}_h$, corresponding to the
discretization of \eqref{sec2:eqnnc}, by the scheme:
\begin{equation}
\begin{split}
& m_j (c_j^{n+1} - c_j^n) + {\Delta t}\sum_{k\in N_j} Q_{jk}(
{\mathbf{u}}^{n+1}_h, c_j^n, c_k^n)   \\
& - {\Delta t}\sum_{k\in N_j} {D}_{jk}({\mathbf{u}}^{n+1}_h,
c_j^{n+1}, c_k^{n+1}) + {\Delta t}  m_j  q_{n,j}^{-}  c_j^n
= {\Delta t}  m_j  {\hat c}_{n,j}  q_{n,j}^{+}
\end{split}\label{sec2:eqc}
\end{equation}
The semi-implicit scheme \eqref{sec2:eqc} writes also:
\begin{equation}
\begin{split}
\Big(1 + {\frac{\Delta t} {m_j}}  B^{n+1}_j\Big) c_j^{n+1}
& = c_j^n\Big(1 - {\frac{\Delta t} {m_j}}  A^n_j\Big)
- {\frac{\Delta t} {m_j}} \sum_{k\in N_j}m_{jk} (u^{n+1}_{j,k})^{-}  c_k^n   \\
&\quad + {\frac{\Delta t} {m_j}} \sum_{k\in N_j}{\frac{1} {d_{jk}}}
 f_{jk}^{n+1}  c_k^{n+1} + {\Delta t}   {\hat c}_{n,j}  q_{n,j}^{+},
\end{split}\label{sec2:e6}
\end{equation}
where
\begin{equation}
A^n_j= \sum_{k\in N_j} m_{jk} (u^{n+1}_{j,k})^{+} +m_j  q_{n,j}^{-},\quad
B^{n+1}_j= \sum_{k\in N_j}  \frac{1} {d_{jk}}  f_{jk}^{n+1},  \label{sec2:cnij}
\end{equation}
for $j\in S_h$, $0\le n\le N-1$. Note that $f_{jk}^{n+1}\ge 0$, so that
$A^n_j\ge 0$ and $B^{n+1}_j\ge 0$.

%\label{sec2:3}
We make the following assumptions on the time step $\Delta t$. Let
\begin{equation}
\theta_{j,k}^n= 1 - 2 {\frac{\Delta t} { m_k}} m_{jk} ( u^{n+1}_{j,k})^{+},  \label{sec2:hyptjk1}
\end{equation}
for $j\in S_h,\; k\in N_j, \; 0\le n\le N-1$. Let $\theta_{0}$ be a
fixed number, $0<\theta_{0}\le 1$. We assume:
\begin{gather}
 {\frac{\Delta t} {m_j}}  A^n_j\le 1 \quad \text{for }
j\in S_h,\; 0\le n\le N-1,  \label{sec2:eqhydt} \\
 \theta_{j,k}^n \ge \theta_{0} \quad \text{for } j\in S_h,\; k\in
N_j,\; 0\le n\le N-1.  \label{sec2:hyptjk2}
\end{gather}
Condition \eqref{sec2:eqhydt} is imposed to ensure $0\leq c_h(x,t)\leq 1$
a.e. in $\Omega _{T}$, and \eqref{sec2:hyptjk2} allows to show an estimate
on the total variation of the function $c_h$. By uniform estimates on $
p_h$ and ${\mathbf{u}}_h$ that will be established in the next section, we
will see that \eqref{sec2:eqhydt} and \eqref{sec2:hyptjk2} are available if
the time step restriction
\begin{equation}
{\Delta t}\leq C h^2,\quad \text{where $C$ is an arbitrary positive
constant,}  \label{sec2:cfl}
\end{equation}
is imposed
 ($h^2$ is due to the fact that the velocity field $\mathbf{u}$
is not in $L^{\infty }$).
 However this restriction is only used for the proof of the stability
and is easily satisfied in practise - see the numerical experiments
in Sect.~\ref{sec4}.

\section{Stability and consistency of the approximation}

\label{sec3} Let $(p_h)$, $({\mathbf{u}}_h)$ and $(c_h)$ be defined
by \eqref{defph}--\eqref{sec2:eq2}, \eqref{sec2:eqc}.
We prove uniform estimates of $(p_h)$ and
$({\mathbf{u}}_h)$ in $L^2(\Omega _{T})$ and $(L^2(\Omega _{T}))^2$,
respectively. Using \eqref{sec2:eqhydt} we show that
$0\leq c_h(x,t)\leq 1$  a.e. in $\Omega _{T}$. Using \eqref{sec2:hyptjk2} we
establish a uniform estimate on the total variation of $c_h$. Finally we
prove that the numerical scheme \eqref{sec2:eq3}, \eqref{sec2:eq4},
\eqref{sec2:eqc} is consistent with the concentration equation \eqref{sec2:eqnnc}.
In the following, we will use often $C$ to represent a generic positive
constant depending only on fixed data.


\subsection{Estimates of the pressure and the flux}
\label{sec3:1}
We have the following result.

\begin{proposition}
\label{lempu} The sequence $(p_h)$ is bounded in
$L^\infty(0,T;L^2(\Omega)) $ and $({\mathbf{u}}_h)$ is bounded
in $L^\infty(0,T;V)$. Moreover, the discrete $H^1$-norm of $(p_h)$ is
uniformly bounded; i.e.,
\begin{equation*}
 \max_{1\le n \le N} \sum_{(j,k)\in {\mathscr E}_h} {\frac{
m_{jk}}{d_{jk} }}  | p^n_j- p^n_k |^2 \le C.
%\label{sec3:n0}
\end{equation*}
\end{proposition}

\begin{proof}
From \eqref{sec2:eq4} we have
\begin{equation*}
\int_{T_j\cup T_k}\mu (c_h^n(x)){\mathbf{u}}
_h^{n+1}(x)\cdot \boldsymbol{\phi}_{\sigma
_{jk}}(x)dx=m_{jk}(p_j^{n+1}-p_k^{n+1}).
\end{equation*}
Since
\begin{equation*}
\big|  \int_{T_j\cup T_k}\mu (c_h^n(x)){\mathbf{u}}
_h^{n+1}(x)\cdot \boldsymbol{\phi}_{\sigma _{jk}}(x)dx\big| ^2\leq
C h^2 \int_{T_j\cup T_k}|{\mathbf{u}}_h^{n+1}(x)|^2\,dx
\end{equation*}
and $m_{jk}\,d_{jk}\geq C h^2$, we have
\begin{equation*}
{\frac{m_{jk}}{d_{jk}}} |p_j^{n+1}-p_k^{n+1}|^2
\leq C \int_{T_j\cup T_k}|{\mathbf{u}}_h^{n+1}(x)|^2\,dx,
\end{equation*}
then
\begin{equation}  \label{sec3:n1}
\begin{split}
\max_{1\leq n\leq N}\sum_{(j,k)\in {\mathscr E}_h}{\frac{
m_{jk}}{d_{jk}}} |p_j^n-p_k^n|^2
&\leq C\max_{1\leq n\leq N}\sum_{j\in S_h}\sum_{k\in N_j}\int_{T_j\cup T_k}|
{\mathbf{u}}_h^n(x)|^2\,dx  \\
&\leq C \| {\mathbf{u}}_h\| _{L^{\infty }(0,T;L^2(\Omega ))}\,.
\end{split}
\end{equation}
Taking $\mathbf{v}_h={\mathbf{u}}_h^{n+1}$ in \eqref{sec2:eq2} and
$\psi_h=p_h^{n+1}$ in \eqref{sec2:eq1} yields
\begin{equation*}
\int_{\Omega }\mu (c_h^n){\mathbf{u}}_h^{n+1}(x)\cdot {\mathbf{u}}
_h^{n+1}(x)dx=\int_{\Omega
}(q_{n,h}^{+}-q_{n,h}^{-})(x)p_h^{n+1}(x)dx.
\end{equation*}
Using \eqref{intro:eqnvis}, the Cauchy-Schwarz inequality and the discrete
Poincar\'{e} inequality for the function $p_h^{n+1}\in W_h^{0}$, we
obtain
\begin{equation} \label{sec3:n2}
\begin{split}
\| {\mathbf{u}}_h^{n+1}\| _{(L^2(\Omega ))^2}^2
 &\leq C \| q_{n,h}^{+}-q_{n,h}^{-}\| _{L^2(\Omega )}\Big( \sum_{(j,k)\in {
\mathscr E}_h}{\frac{m_{jk}}{d_{jk}}} |p_j^{n+1}-p_k^{n+1}|^2
\Big) ^{1/2}   \\
&\leq  C \Big( \sum_{(j,k)\in {\mathscr E}_h}{\frac{m_{jk}}{d_{jk}}}
|p_j^{n+1}-p_k^{n+1}|^2\Big) ^{1/2}.
\end{split}
\end{equation}
We deduce from \eqref{sec3:n1} and \eqref{sec3:n2} that
\begin{equation}
\max_{1\leq n\leq N}\sum_{(j,k)\in {\mathscr E}_h}{\frac{
m_{jk}}{d_{jk}}} |p_j^n-p_k^n|^2\leq C,\quad
\| {\mathbf{u}} _h\| _{L^{\infty }(0,T;L^2(\Omega ))}\leq C .  \label{sec3:n3}
\end{equation}
Taking $\psi _h=\operatorname{div}{\mathbf{u}}_h^{n+1}$ in
\eqref{sec2:eq1} yields
\begin{equation*}
\int_{\Omega }|\operatorname{div}{\mathbf{u}}_h^{n+1}(x)|^2\,dx=
\int_{\Omega }(q_{n,h}^{+}-q_{n,h}^{-})(x)\operatorname{div}{\mathbf{u}}
_h^{n+1}(x)dx,
\end{equation*}
which implies, by the Cauchy-Schwarz inequality,
\begin{equation*}
\| \operatorname{div}{\mathbf{u}}_h^{n+1}\| _{L^2(\Omega )}\leq
\| q_{n,h}^{+}-q_{n,h}^{-}\| _{L^2(\Omega )}.
\end{equation*}
We conclude that $(\operatorname{div}{\mathbf{u}}_h)$ is bounded in
$L^{\infty }(0,T;L^2(\Omega ))$ and then $({\mathbf{u}}_h)$ is bounded in
$L^2(0,T;V)$.
\end{proof}

Since the scheme \eqref{sec2:eqc} involves numerical fluxes, we need the
following estimates.

\begin{lemma}\label{lemflux}
There exists a positive number $C$ such that
\begin{equation}
\sum_{k\in N_j} \int_{\sigma_{jk}}| {\mathbf{u}}
^{n+1}_h(x)\cdot \nu_{jk}|^2 \, d\gamma(x) \le C  h^{-1}
 \| {\mathbf{u}}^{n+1}_h \|_{L^2(T_j)}^2,  \label{sec3:n4}
\end{equation}
and
\begin{equation}
\max_{1\le n\le N} \Big( \sum_{(j,k)\in {\mathscr E}_h }
\int_{\sigma_{jk}}| {\mathbf{u}}^n_h(x)\cdot \nu_{jk}| \, d\gamma(x)
\Big) \le C h^{-1}.  \label{sec3:n5}
\end{equation}
uniformly with respect to $h$, ${\Delta t}$ and $j \in S_h$.
\end{lemma}

\begin{proof}
We use the property of ${\mathbf{u}}_h$, that is the restriction of
 ${\mathbf{u}}_h$ to any triangle $T_j$ in $S_h$ belongs to
$L^\infty(0,T; H^{1}(T_j))$.
We will make use of the following well-known result(a proof of which
is attached as an appendix at the end for the convenience of the reader.),
 there is a constant $C$ such that, for $\varphi\in H^{1}(T_j)$,
the local inverse estimate holds:
\begin{equation*}
\| \varphi \|^2_{L^2(\partial T_j)} \le C  \|
\varphi\|_{L^2(T_j)}  \left( \| \operatorname{grad}
\varphi\|_{L^2(T_j)} + h^{-1} \| \varphi\|_{L^2(T_j)}\right).
\end{equation*}
We have
\begin{align*}
&\sum_{k\in N_j} \int_{\sigma_{jk}}| {\mathbf{u}}^{n+1}_h(x)\cdot
\nu_{jk}|^2 \, d\gamma(x) \\
&\le C  \| {\mathbf{u}}^{n+1}_h\|_{L^2(T_j)} \left( | {\mathbf{u}}
^{n+1}_h|_{1,T_j} + h^{-1} \| {\mathbf{u}}^{n+1}_h
\|_{L^2(T_j)}\right),
\end{align*}
where $| {\mathbf{u}}^{n+1}_h|_{1,T_j}$ is the
semi-norm involving the $L^2$-norms of the spatial derivatives of order 1 of
${\mathbf{u}}^{n+1}_h$. Recall that
\begin{equation*}
{\mathbf{u}}^{n+1}_h\big|_{T_j} =  \sum_{k\in N_j} {u}
^{n+1}_{j,k}  \boldsymbol{\phi}_{\sigma_{jk}}.
\end{equation*}
Using properties of the basis functions $\boldsymbol{\phi}_{\sigma_{jk}}$, we have
\begin{equation*}
 | {\mathbf{u}}^{n+1}_h|_{1,T_j}^2
\le C \sum_{k\in N_j} | {u}^{n+1}_{j,k} |^2.
\end{equation*}
Since
\begin{equation*}
| {u}^{n+1}_{j,k} |^2 \le C h^{-1}  \int_{\sigma_{jk}}
| {\mathbf{u}}^{n+1}_h(x)\cdot \nu_{jk}|^2 \, d\gamma(x),
\end{equation*}
we obtain
\begin{equation*}
 | {\mathbf{u}}^{n+1}_h|_{1,T_j} \le C  \Big(
h^{-1}\sum_{k\in N_j} \int_{\sigma_{jk}}| {\mathbf{u}}^{n+1}_h(x)\cdot
\nu_{jk}|^2 \, d\gamma(x) \Big)^{1/2}.
\end{equation*}
Thus
\begin{align*}
&\sum_{k\in N_j} \int_{\sigma_{jk}}| {\mathbf{u}}^{n+1}_h(x)\cdot
\nu_{jk}|^2 \, d\gamma(x) \\
&\le  C  \| {\mathbf{u}}^{n+1}_h \|_{L^2(T_j)}
 \Big( \Big(h^{-1}\sum_{k\in N_j} \int_{\sigma_{jk}}| {
\mathbf{u}}^{n+1}_h(x)\cdot \nu_{jk}|^2 \, d\gamma(x) \Big)^{1/2} \\
&\quad + h^{-1} \| {\mathbf{u}}^{n+1}_h \|_{L^2(T_j)}\Big),
\end{align*}
and using the inequality $2 a b\le a^2+b^2$, we get the first estimate
\eqref{sec3:n4}. Summing \eqref{sec3:n4} over all triangles $T_j$, we
obtain
\[
\sum_{j\in S_h} \sum_{k\in N_j} \int_{\sigma_{jk}}| {\mathbf{u}}
^{n+1}_h(x)\cdot \nu_{jk}|^2 \, d\gamma(x)
\le C  h^{-1}  \| {\mathbf{u}}^{n+1}_h \|_{L^2(\Omega)}^2.
\]
Using the inequality
\begin{align*}
& \Big( \sum_{(j,k)\in {\mathscr E}_h } \int_{\sigma_{jk}}| {\mathbf{u}}
^{n+1}_h(x)\cdot \nu_{jk}| \, d\gamma(x) \Big)^2 \\
&\le \Big( \sum_{(j,k)\in {\mathscr E}_h } m_{jk} \Big)
\Big(\sum_{(j,k)\in {\mathscr E}_h } \int_{\sigma_{jk}} | {\mathbf{u}}
^{n+1}_h(x)\cdot \nu_{jk}|^2 \, d\gamma(x) \Big),
\end{align*}
and $ \sum_{(j,k)\in {\mathscr E}_h} m_{jk} \le C  h^{-1}$,
it follows that
\begin{equation*}
\max_{1\le n\le N} \Big( \sum_{(j,k)\in {\mathscr E}_h }
\int_{\sigma_{jk}}| {\mathbf{u}}^n_h(x)\cdot \nu_{jk}| \, d\gamma(x)
\Big) \le C h^{-1}  \max_{1\le n\le N} \| {\mathbf{u}}^n_h\|_{L^2(\Omega)},
\end{equation*}
from which, according to the boundedness of $(u_h)$ in
$(L^\infty(0,T;L^2(\Omega)))^2$, we deduce the estimate on the flux \eqref{sec3:n5}.
The proof of Lemma~\ref{lemflux} is complete.
\end{proof}

\subsection{Estimates for the concentration}
\label{sec3:2} Clearly, $A_j^n$ defined in \eqref{sec2:cnij} is bounded
by a constant $C_h$ which may depend on $h$ but does not depend on $j\in
S_h$ and on $0\leq n\leq N-1$. By the Cauchy-Schwarz inequality we have
\begin{align*}
\sum_{k\in N_j}m_{jk}({\mathbf{u}}_{j,k}^{n+1})^{+}
&\leq \sum_{k\in N_j}\int_{\sigma _{jk}}|{\mathbf{u}}_h^{n+1}(x)\cdot
\nu _{jk}|\,d\gamma (x) \\
&\leq C h^{1/2}\Big( \sum_{k\in N_j}\int_{\sigma _{jk}}|{
\mathbf{u}}_h^{n+1}(x)\cdot \nu _{jk}|^2\,d\gamma (x)\Big) ^{1/2}.
\end{align*}
Using \eqref{sec3:n4} we deduce that
\[
\sum_{k\in N_j}m_{jk}(u_{j,k}^{n+1})^{+}
\leq C \| {\mathbf{u}} _h^{n+1}\| _{L^2(T_j)}
\]
for all $j\in S_h$ and  $n$, $0\leq n\leq N-1$. Since $({\mathbf{u}}_h)$
is bounded in $(L^{\infty }(0,T;L^2(\Omega )))^2$,
 $\| {\mathbf{u}} _h^{n+1}\| _{L^2(T_j)}$ tends to 0 as $h$ and ${\Delta t}$ tend to
0. We also have, for any $j\in S_h$ and $0\leq n\leq N-1$,
\begin{equation*}
m_j q_{n,j}^{-}\leq C {\frac{h}{{\Delta t}^{1/2}}}\Big(
\int_{t_{n}}^{t_{n+1}}\int_{T_j}| q_{n,h}^{-}(x)|^2\,dx\,dt\Big) ^{1/2},
\end{equation*}
and $\int_{t_{n}}^{t_{n+1}}\int_{T_j}|q_{n,h}^{-}(x)| ^2\,dx\,dt$ tends to 0
as $h$ and ${\Delta t}$ tend to 0.
For the coefficient $B_j^{n+1}$ involving the term
$\sum_{k\in N_j}{\frac{1}{d_{jk}}} f_{jk}^{n+1}$, we note as in the proof
of Lemma~\ref{lemflux} that
\begin{equation*}
\sum_{k\in N_j}{\frac{1}{d_{jk}}} f_{jk}^{n+1}\leq
C h^{-1} \| {\mathbf{u}}_h^{n+1}\| _{L^2(T_j)},  \label{sec3:ed5p}
\end{equation*}
thus
\begin{equation*}
\max_{1\leq n\leq N}\sum_{j\in S_h}\sum_{k\in N_j}{\frac{1}{d_{jk}}}
 |f_{jk}^n|\leq C h^{-1} \max_{1\leq n\leq N}\| u_h^n\|
_{L^2(\Omega )}.  \label{sec3:ed5}
\end{equation*}
Consequently, if we suppose
\begin{equation*}
{\Delta t} \le C  h^2, \text{ where $C$ is an a positive
constant,}  \label{sec3:n6}
\end{equation*}
then  relations \eqref{sec2:eqhydt} and \eqref{sec2:hyptjk2} hold.

We remark that
if  $q^{+}$ and $q^{-}$ belong to $L^\infty(\Omega_T)$,
 then \eqref{sec2:eqhydt} and \eqref{sec2:hyptjk2} are
fulfilled under the condition
$\Delta t \le C  h$ for a suitable positive constant $C$.
The following proposition provides the stability of the scheme \eqref{sec2:eqc} in
the $L^{\infty}(\Omega_T)$ norm.

\begin{proposition}\label{lemcon}
The sequence $(c_h)$ constructed by the numerical scheme
\eqref{sec2:eqc} satisfies
\begin{equation*}
0\le c_h(x,t)\le 1 \quad \text{a.e. in } \Omega_T,
\end{equation*}
provided the stability condition \eqref{sec2:eqhydt} is valid.
\end{proposition}

\begin{proof}
Starting from $(c_j^n)$ with $0\le c_j^n\le 1$, and according to
\eqref{intro:eqnsourcep}, \eqref{intro:eqnsourcec} and \eqref{sec2:eqhydt},
we deduce from the semi-implicit scheme \eqref{sec2:e6} that
$c_j^{n+1}\ge 0$ for all $j\in S_h$. Using \eqref{sec2:eq3},
we transform \eqref{sec2:e6} to obtain
\begin{align*}
\Big(1 + {\frac{\Delta t} {m_j}}  B^{n+1}_j\Big) (1- c_j^{n+1})
&= (1- c_j^n)\Big(1 - {\frac{\Delta t} {m_j}} A^n_j\Big) -
{\frac{\Delta t} {m_j}} \sum_{k\in N_j} m_{jk} (u^{n+1}_{j,k})^{-}
(1- c_k^n) \\
&\quad + {\frac{\Delta t} {m_j}} \sum_{k\in N_j}{\frac{1} {d_{jk}}}
 f_{jk}^{n+1} (1- c_k^{n+1}) + {\Delta t} (1- {\hat c}_{n,j})
q_{n,j}^{+},
\end{align*}
then as above we infer that
$1- c_j^{n+1}\ge 0$ for all $j\in S_h$. We conclude that
\begin{equation*}
0 \le c_j^n\le 1 \quad \text{for all $j\in S_h$  and $0\le n\le N-1$}.
\end{equation*}
This completes the proof.
\end{proof}

At time level $t_{n}$, for any edge in $\mathscr {E}_h$, the
flux $u^n_{j,k}$ associated with the edge $\sigma_{jk}$ is either greater than $ 0$ or
less than or equal to $0$. We then select the edges such that
\begin{equation}
{\mathscr A}^n=\{ (j,k), \; j\in S_h,\; k\in N_j \,; \;
u^n_{j,k}>0 \} .  \label{sec3:n7}
\end{equation}
We introduce the total variation of
the function $c_h$ by
\begin{equation*}
\operatorname{TV}(c_h) = \sum_{n=0}^{N-1} {\Delta t}\sum_{(j,k)\in {\mathscr A}
^{n+1}} m_{jk} u^{n+1}_{j,k}   | c_k^n- c_j^n|.
\label{sec3:n8}
\end{equation*}
We have the following estimate.

\begin{lemma}\label{lemvbc}
There is a number $C$, independent of $h$ and ${\Delta t}$,such that
\begin{equation}
\operatorname{TV}(c_h) \le C h^{-1/2}.  \label{sec3:n9}
\end{equation}
\end{lemma}

\begin{proof}
Let us write \eqref{sec2:eqc} in the form
\begin{equation}
c_j^{n+1} - r_j^n = \frac{{\Delta t}}{m_j} \sum_{k\in N_j}
D_{jk}({\mathbf{u}}^{n+1}_h, c_j^{n+1}, c_k^{n+1}) \label{sec3:n11}
\end{equation}
with
\begin{equation}
r_j^n= c_j^n - \frac{{\Delta t}}{m_j} \sum_{k\in
N_j} Q_{jk}( {\mathbf{u}}^{n+1}_h, c_j^n, c_k^n) - {\Delta t}
q_{n,j}^{-} c_j^n + {\Delta t}  q_{n,j}^{+}  {\hat c}_{n,j}.
\label{sec3:n12}
\end{equation}
Using \eqref{sec2:eq3}, we transform \eqref{sec3:n12} into
\begin{equation}
r_j^n= c_j^n - \frac{{\Delta t}}{m_j} \sum_{k\in
N_j} \left( Q_{jk}( {\mathbf{u}}^{n+1}_h, c_j^n, c_k^n) - m_{jk}
u^{n+1}_{j,k} c_j^n\right) + {\Delta t} g_{n,j} \label{sec3:n13}
\end{equation}
where $g_{n,j}= q_{n,j}^{+} ({\hat c}_{n,j} - c_j^n)$.
Multiplying \eqref{sec3:n11} by $c_j^{n+1}$ and using the identity
\[
(c^{n+1}_j - r^n_j) c^{n+1}_j =  {\frac{1}{2}}
(c^{n+1}_j - r^n_j)^2 +  {\frac{1}{2}} (c^{n+1}_j)^2-
 {\frac{1}{2}} (r^n_j)^2,
\]
we obtain
\[
(c^{n+1}_j)^2 +(c^{n+1}_j - r^n_j)^2 - (r^n_j)^2
=2 \frac{{\Delta t}}{m_j} \sum_{k\in N_j} D_{jk}({\mathbf{u}}^{n+1}_h,
c_j^{n+1}, c_k^{n+1}) c^{n+1}_j.  %\label{sec3:n14}
\]
We observe that
\[
2\sum_{j\in S_h} \sum_{k\in N_j} D_{jk}({\mathbf{u}}^{n+1}_h,
c_j^{n+1}, c_k^{n+1}) c^{n+1}_j = -  \sum_{j\in S_h}
 \sum_{k\in N_j} \frac{f_{jk}^{n+1}}{d_{jk}}  (c^{n+1}_k
- c^{n+1}_j)^2.  \label{sec3:n15}
\]
According to \eqref{sec3:n13}, we have
\begin{align*}
(r^n_j)^2&= (c^n_j)^2 + \Big(\frac{{\Delta t}}{m_j}
\sum_{k\in N_j} \left( Q_{jk}( {\mathbf{u}}^{n+1}_h, c_j^n,
c_k^n) - m_{jk} u^{n+1}_{j,k} c_j^n\right) \Big)^2 + {\Delta t}
^2(g_{n,j})^2 \\
&\quad - 2 \frac{{\Delta t}}{m_j} \sum_{k\in N_j} \left(
Q_{jk}( {\mathbf{u}}^{n+1}_h, c_j^n, c_k^n) - m_{jk} u^{n+1}_{j,k} c_j^n\right)c_j^n \\
&\quad + 2{\Delta t} c_j^n g_{n,j} - 2{\Delta t}   g_{n,j} 
\frac{{\Delta t}}{m_j} \sum_{k\in N_j}
\left(Q_{jk}({\mathbf{u}}^{n+1}_h, c_j^n, c_k^n) - m_{jk} u^{n+1}_{j,k} c_j^n\right).
\end{align*}
Simple computations give the  inequality
\begin{align*}
(r^n_j)^2
&\le  (c^n_j)^2 + 2 \Big(
\frac{{\Delta t}}{m_j} \sum_{k\in N_j} \left( Q_{jk}( {\mathbf{u}}^{n+1}_h, c_j^n,
c_k^n) - m_{jk} u^{n+1}_{j,k} c_j^n\right) \Big)^2 +2{\Delta
t}^2(g_{n,j})^2   \\
&\quad - 2\frac{{\Delta t}}{m_j} \sum_{k\in N_j} \left(
Q_{jk}( {\mathbf{u}}^{n+1}_h, c_j^n, c_k^n) - m_{jk}
u^{n+1}_{j,k} c_j^n\right)c_j^n + 2{\Delta t}|
g_{n,j}|. % \label{sec3:n166}
\end{align*}
By definition of $Q_{jk} $ we have
\begin{equation}
\sum_{k\in N_j} Q_{jk}({\mathbf{u}}^{n+1}_h,c_j^n,c_k^n)
c^n_j=\sum_{k\in N_j} m_{jk}(u^{n+1}_{j,k})^{+} ( c_j^n)^2
+\sum_{k\in N_j } m_{jk}(u^{n+1}_{j,k})^{-} c_k^nc_j^n  \label{sec3:n17}
\end{equation}
According to the identity
\begin{equation*}
c^n_j c^n_k =  {\frac{1}{2}} (c^n_j)^2 +
 {\frac{1}{2}} (c^n_k)^2 - {\frac{1}{2}}(c^n_j - c^n_k)^2,
\end{equation*}
we write \eqref{sec3:n17} in the form
\begin{equation}
\begin{split}
&\sum_{k\in N_j} Q_{jk}({\mathbf{u}}^{n+1}_h, c_j^n, c_k^n)  c_j^n \\
&=  -  {\frac{1}{2}} \sum_{k\in N_j} m_{jk} (u^{n+1}_{j,k})^{-}  (c^n_j - c^n_k)^2
+  {\frac{1}{2}} \sum_{k\in N_j } m_{jk} (u^{n+1}_{j,k})^{-} ( c_k^n)^2   \\
&\quad +  {\frac{1}{2}} \sum_{k\in N_j} m_{jk} (u^{n+1}_{j,k})^{-}  (c_j^n)^2
 + \sum_{k\in N_j} m_{jk} (u^{n+1}_{j,k})^{+}  (c_j^n)^2 .
\end{split} \label{sec3:n18}
\end{equation}
Since
\begin{equation*}
 \sum_{j\in S_h} \sum_{k\in N_j} m_{jk} (u^{n+1}_{j,k})^{-}  (c^n_k)^2
= - \sum_{(j,k)\in {\mathscr A}^{n+1}} m_{jk} u^{n+1}_{j,k}  (c_j^n )^2.
\end{equation*}
 From \eqref{sec3:n18} we deduce that
\begin{align*}
&\sum_{j\in S_h} \sum_{k\in N_j} Q_{jk}({\mathbf{u}}^{n+1}_h, c_j^n,
c_k^n)  c_j^n\\
&= -  {\frac{1}{2}} \sum_{j\in S_h}
\sum_{k\in N_j } m_{jk} (u^{n+1}_{j,k})^{-}  (c^n_j - c^n_k)^2
+  {\frac{1}{2}} \sum_{j\in S_h} \sum_{k\in N_j}m_{jk} u^{n+1}_{j,k} ( c_j^n)^2 .
\end{align*}
Using \eqref{sec2:eq3} and the identity
\begin{equation*}
\sum_{j\in S_h}\sum_{k\in N_j} m_{jk} (u^{n+1}_{j,k})^{-} (c^n_j
- c^n_k)^2 = - \sum_{(j,k)\in {\mathscr A}^{n+1}} m_{jk} u^{n+1}_{j,k}
 (c^n_j - c^n_k)^2,
\end{equation*}
we obtain
\begin{align*}
&\sum_{j\in S_h}  \sum_{k\in N_j} Q_{jk}(\mathbf{u}^{n+1}_h, c_j^n, c_k^n)  c_j^n\\
&=  {\frac{1}{2}} \sum_{(j,k)\in {\mathscr A}^{n+1}} m_{jk} u^{n+1}_{j,k}
  (c^n_j - c^n_k)^2
 +  {\frac{1}{2}} \sum_{ j\in S_h} m_j ( q^{+}_{n,j}- q^{-}_{n,j}) (c_j^n )^2.
\end{align*}
Thus
\[
 2\sum_{j\in S_h} \sum_{k\in N_j} \left( Q_{jk}( {\mathbf{u}}^{n+1}_h,
c_j^n, c_k^n) - m_{jk} u^{n+1}_{j,k} c_j^n\right)c_j^n
= \sum_{(j,k)\in {\mathscr A}^{n+1}} m_{jk} u^{n+1}_{j,k}  (c^n_j -
c^n_k)^2 .
\]
We also have the relations
\begin{equation*}
 \sum_{k\in N_j} \left(Q_{jk}({\mathbf{u}}^{n+1}_h, c_j^n,
c_k^n)- m_{jk} u^{n+1}_{j,k} c_j^n\right)
 = \sum_{k\in N_j} m_{jk} (u^{n+1}_{j,k})^{-} ( c_k^n- c_j^n),
\end{equation*}
and
\begin{align*}
&\sum_{j\in S_h} {\frac{1} {m_j}}
\Big(\sum_{k\in N_j} Q_{jk}({\mathbf{u}}^{n+1}_h, c_j^n, c_k^n)-
m_{jk} u^{n+1}_{j,k} c_j^n \Big)^2\\
&=  \sum_{j\in S_h} {\frac{1} {m_j}} \Big(
\sum_{k\in N_j} m_{jk} (u^{n+1}_{j,k})^{-} ( c_k^n- c_j^n)
\Big)^2
 =  \sum_{(j,k)\in {\mathscr A}^{n+1}} {\frac{1} {m_k}}
( m_{jk} ( c_k^n- c_j^n)u^{n+1}_{j,k})^2.
\end{align*}
Now, multiplying \eqref{sec3:n13} by $m_j$, summing over
 $j\in S_h $, and using the above results we obtain
\begin{align*}
& \sum_{j\in S_h} m_j (c^{n+1}_j)^2 + \sum_{j\in S_h}
m_j(c^{n+1}_j - r^n_j)^2 +{\Delta t}\sum_{j\in
S_h}  \sum_{k\in N_j} \frac{f_{jk}^{n+1}}{d_{jk}}  (c^{n+1}_k - c^{n+1}_j)^2   \\
& +{\Delta t}\sum_{(j,k)\in {\mathscr A}^{n+1}}
\Big(1 -2 \frac{\Delta t}{m_k}  m_{jk} u^{n+1}_{j,k} \Big)
 m_{jk} u^{n+1}_{j,k}  ( c_k^n- c_j^n)^2   \\
& \le \sum_{j\in S_h} m_j (c^n_j)^2 +2 ({\Delta t})^2\sum_{j\in S_h} m_j(g_{n,j})^2
 + {\Delta t}\sum_{j\in S_h} m_j | g_{n,j}|.  %\label{sec3:n19}
\end{align*}
Summing from $n=0$ to $N-1$, we obtain
\begin{equation}
\begin{split}
&\sum_{j\in S_h} m_j (c^{N}_j)^2 +
\sum_{n=0}^{N-1}\sum_{j\in S_h} m_j(c^{n+1}_j - r^n_j)^2
 + \sum_{n=0}^{N-1} {\Delta t}\sum_{j\in
S_h}  \sum_{k\in N_j} \frac{f_{jk}^{n+1}}{d_{jk}}  (c^{n+1}_k - c^{n+1}_j)^2   \\
&+ \sum_{n=0}^{N-1}{\Delta t}\sum_{(j,k)\in {
\mathscr A}^{n+1}} \Big(1 - 2 \frac{\Delta t}{m_k}  m_{jk} u^{n+1}_{j,k} \Big)
 m_{jk} u^{n+1}_{j,k}  ( c_k^n- c_j^n)^2
 \\
&\le \sum_{j\in S_h} m_j(c^{0}_j)^2 +2 {\Delta t}
 \sum_{n=0}^{N-1}{\Delta t}  \sum_{j\in S_h}
m_j(g_{n,j})^2 +  \sum_{n=0}^{N-1}{\Delta t}
 \sum_{j\in S_h} m_j| g_{n,j}|
\end{split} \label{sec3:n20}
\end{equation}
The right-hand side of \eqref{sec3:n20} is bounded by a constant $K$ which
depends only on the data.
Bearing in mind the definition \eqref{sec2:hyptjk1} of $\theta_{j,k}^n$
which is greater than some $\theta_{0}>0$ we derive the following estimates
\begin{gather}
\sum_{n=0}^{N-1} {\Delta t}\sum_{(j,k)\in {\mathscr A}^{n+1}} m_{jk} u^{n+1}_{j,k}  ( c_k^n- c_j^n)^2 \le  {\ \frac{K} {
\theta_{0}}},  \label{sec3:n21}
\\
\sum_{n=0}^{N-1} {\Delta t}\sum_{(j,k)\in {\mathscr E}_h} f_{jk}^{n+1} 
\frac{(c_k^{n+1}-c_j^{n+1})^2}{d_{jk}} \le K.  \label{sec3:n22}
\end{gather}
Note that the latter estimate corresponds to the discrete approximation
of the function $D^{1/2}({\mathbf{u}}_h)\operatorname{grad} c_h$ in the space
 $(L^2(\Omega_T))^2 $.
By the Cauchy-Schwarz inequality, we infer that
\begin{align*}
&\operatorname{TV}(c_h)\\
& \le  \Big( \sum_{n=0}^{N-1} {\Delta t}\sum_{(j,k)
 \in {\mathscr A}^{n+1}}m_{jk} u^{n+1}_{j,k}  (c_k^n- c_j^n)^2\Big)^{1/2}
\Big(  \sum_{n=0}^{N-1} {\Delta t}\sum_{(j,k)\in {\mathscr A}^{n+1}}
  m_{jk} u^{n+1}_{j,k} \Big)^{1/2} \\
& \le  {\frac{K} {\theta_{0}}}  \Big( \sum_{n=0}^{N-1} {
\Delta t}\sum_{(j,k)\in {\mathscr A}^{n+1}} m_{jk} u^{n+1}_{j,k}
\Big)^{1/2}.
\end{align*}
To estimate the right-hand side of the latter inequality we note that
\[
\max_{1\le n\le N} \sum_{(j,k)\in {\mathscr A}^n} m_{jk} u^n_{j,k} \le
C  \max_{1\le n\le N} \sum_{(j,k)\in {\mathscr E}_h }
\int_{\sigma_{jk}}| {\mathbf{u}}^n_h(x)\cdot \nu_{jk}| \, d\gamma(x).
\]
Using the estimate on the flux \eqref{sec3:n5} we obtain the desired
estimate. The proof of Lemma~\ref{lemvbc} is complete.
\end{proof}

\subsection{Consistency of the numerical scheme \eqref{sec2:eq3},
\eqref{sec2:eq4}, \eqref{sec2:eqc} with the concentration equation \eqref{sec2:eqnnc}}
\label{sec3:4}

Let $c_h(x,t)$ be the approximate concentration obtained by the numerical scheme
\eqref{sec2:eq3},\eqref{sec2:eq4}, \eqref{sec2:eqc}; we define the
consistency error $\epsilon(h)\ $ to be:
\begin{equation}
\partial_{t}c_h+\operatorname{div}(c_h {\mathbf{u}}_h-D({\mathbf{u}}
_h)\operatorname{grad}c_h)+{c}_h q_h^{-}-{\hat{c}}
_h q_h^{+}=\epsilon(h),\label{sec3:eqnconservh}
\end{equation}
and we will show that $\lim_{h\to0}\varepsilon(h)=0$.
The exact concentration $c(x,t)$ satisfy the equation \eqref{sec2:eqnnc}
whose right hand side is equal to zero. Equation \eqref{sec3:eqnconservh}
will be understood in distributional sense since $c(x,t)$ satisfy also
the equation \eqref{intro:eqnconfaible}.
For this let $\varphi \in C^2(\overline{\Omega}_T) $ with compact
support contained in $\overline{\Omega}\times [0,T[$ and let
\[
{\tilde \varphi}_n(x) = {\frac{1}{\Delta t}}
 \int_{t_n}^{t_{n+1}} \varphi(x,t) \,dt,
\quad
 \varphi_{n,j}= {\frac{1} {m_j}} \int_{T_j} {\tilde
\varphi}_n(x) \,dx.
\]
 We multiply the discretized equation \eqref{sec2:eqc}
by ${\tilde \varphi}_n(x)/ m_j$, integrate on $T_j$ and take the sum
over $j\in S_h$ and over $n$, from 0 to $N-1$. This yields
\begin{equation}
E^h \equiv E^h_1 + E^h_{2} + E^h_3 + E^h_{4}=0,
\label{sec3:n25}
\end{equation}
with
\begin{gather*}
E^h_1 = \sum_{n=0}^{N-1} \sum_{j\in S_h} m_j (c_j^{n+1}
 - c_j^n) \varphi_{n,j}; \\
 E^h_{2}=  \sum_{n=0}^{N-1} {\Delta t}\sum_{j\in S_h}
\sum_{k\in N_j} Q_{jk}({\mathbf{u}}^{n+1}_h, c_j^n, c_k^n)  \varphi_{n,j}; \\
 E^h_3= -  \sum_{n=0}^{N-1} {\Delta t}\sum_{j\in S_h}
\sum_{k\in N_j} {D}_{jk}({\mathbf{u}}^{n+1}_h, c_j^{n+1}, c_k^{n+1})  \varphi_{n,j}; \\
 E^h_{4}=  \sum_{n=0}^{N-1} {\Delta t}\sum_{j\in S_h}
m_j q_{n,j}^{-}  c_j^n \varphi_{n,j} -
\sum_{n=0}^{N-1} {\Delta t}\sum_{j\in S_h} m_j {\hat c}_{n,j}
q_{n,j}^{+}  \varphi_{n,j}.
\end{gather*}
Let us first consider $E_1^h$. We have
\begin{equation*}
\sum_{n=0}^{N-1} \sum_{j\in S_h} m_j (c_j^{n+1} - c_j^n)  \varphi_{n,j}
= - \sum_{n=1}^{N-1} \sum_{j\in S_h} m_j (\varphi_{n,j}
- \varphi_{n-1,j})c_j^n - \sum_{j\in S_h}
m_j\varphi_{0,j}  c_j^{0},
\end{equation*}
then
\begin{equation*}
 E^h_1= - \sum_{n=1}^{N-1} \int_{\Omega} c_h^n(x,t)
(\varphi_{n,h}(x) - \varphi_{n-1,h}(x)) \,dx - \int_{\Omega}
c^{0}_h(x) \varphi(x,0)dx.  \label{sec3:n26}
\end{equation*}
Clearly,
\begin{equation}
E^h_{4} = \sum_{n=0}^{N-1} {\Delta t}  \int_{\Omega} c_{n,h}(x,t)
q^{-}_{n,h}(x,t)   {\tilde \varphi}_n(x) \, dx - \sum_{n=0}^{N-1} {\Delta t}
  \int_{\Omega} {\hat c}_h  q^{+}_{n,h}(x,t)  {\tilde \varphi}_n(x) \,
dx.  \label{sec3:n27}
\end{equation}
Consider now $E^h_{2}$. Multiplying \eqref{sec2:eq3} by $c_j^n$ and using
the identity
\begin{equation*}
 u^{n+1}_{j,k} = (u^{n+1}_{j,k})^{+} + (u^{n+1}_{j,k})^{-},
\end{equation*}
we have
\[
 {\Delta t} \sum_{k\in N_j} m_{jk} c_j^n ( u^{n+1}_{j,k})^{+}
= - {\Delta t} \sum_{k\in N_j} m_{jk}  c_j^n (
u^{n+1}_{j,k})^{-} + {\Delta t}   m_j c_j^n  ( q^{+}_{n,j} -
q^{-}_{n,j}).
\]
Using \eqref{sec2:eq5flu}, one can write $E^h_{2}$ in the form
$E^h_{2}= E^h_{21} + E^h_{22}$ with
\begin{align*}
E^h_{21} &=   \sum_{n=0}^{N-1} {\Delta t} \sum_{j\in
S_h} \sum_{k\in N_j} m_{jk} (c_k^n- c_j^n)(
u^{n+1}_{j,k})^{-} \varphi_{n,j} \\
& = -  \sum_{n=0}^{N-1} {\Delta t} \sum_{(j,k)\in {\mathscr
A}^{n+1}} m_{jk} (c_j^n- c_k^n)u^{n+1}_{j,k} \varphi_{n,k},
\end{align*}
and
\begin{equation*}
E^h_{22} =   \sum_{n=0}^{N-1} {\Delta t} \sum_{j\in
S_h} c_j^n   \int_{T_j} {\tilde \varphi}_n(x) \operatorname{div}
 {\mathbf{u}}_h^{n+1}(x) \,dx.
\end{equation*}
Let us define
\begin{equation*}
F^h_{2} = - \sum_{n=0}^{N-1} {\Delta t}  \int_{\Omega} c_h^n(x){
\mathbf{u}}_h^{n+1}(x) \cdot \operatorname{grad} {\tilde \varphi}_n(x)
\,dx.
\end{equation*}
One can write $F^h_{2} = F^h_{21} + F^h_{22}$ with
\begin{align*}
F^h_{21} &=  - \sum_{n=0}^{N-1} {\Delta t} \sum_{j\in S_h}
c_j^n\sum_{k\in N_j} \int_{\sigma_{jk}} {\tilde \varphi}_n(x) {
\mathbf{u}}_h^{n+1} \cdot \nu_{jk}\, d\gamma(x) \\
&= -\sum_{n=0}^{N-1} {\Delta t} \sum_{(j,k)\in {\mathscr A}^{n+1}}
(c_j^n- c_k^n)u^{n+1}_{j,k}   \int_{\sigma_{jk}}
{\tilde \varphi}_n(x) d\gamma(x),
\end{align*}
and
\begin{equation*}
F^h_{22} = \sum_{n=0}^{N-1} {\Delta t} \sum_{j\in S_h} c_j^n 
\int_{T_j} {\tilde \varphi}_n(x) \operatorname{div} {\mathbf{u}}
_h^{n+1}(x) \,dx.
\end{equation*}
To compare $E^h_{2}$ and $F^h_{2}$, we introduce
$R^h_{2}\equiv E^h_{2} -F^h_{2}$; i.e.,
\begin{equation*}
R^h_{2} =\sum_{n=0}^{N-1} {\Delta t} \sum_{(j,k)\in {\mathscr A}^{n+1}}
m_{jk} (c_j^n- c_k^n)u^{n+1}_{j,k}  \Big( \frac{1}{ m_{jk}}
\int_{\sigma_{jk}} {\tilde \varphi}_n(x) d\gamma(x) - \varphi_{n,k}\Big).
\end{equation*}
According to the regularity assumption on $\varphi$, we have
$$
\big| \varphi_{n,k}- \frac{1}{ m_{jk}}\int_{\sigma_{jk}} {\tilde \varphi
}_n(x) d\gamma(x)\big| \le C h;
$$
therefore,
\begin{equation}
| R^h_{2}| \le C h \operatorname{TV}(c_h).  \label{sec3:n29}
\end{equation}

Now we discuss the term $E^h_3$, proceeding as previously. We introduce
\begin{equation}
F_3^h = -  \sum_{n=0}^{N-1} {\Delta t}  \int_{\Omega}
c_h^{n+1}  \operatorname{div} ( D({\mathbf{u}}^{n+1}_h)
\operatorname{grad} {\tilde \varphi}_n(x))dx.  \label{sec3:n30}
\end{equation}
Note that $F_3^h$, which has a meaning, represents formally the quantity
\begin{equation}
 \sum_{n=0}^{N-1} {\Delta t}  \int_{\Omega} D({\mathbf{u}}
^{n+1}_h) \operatorname{grad} c_h^{n+1}\cdot \operatorname{grad}
 {\tilde \varphi}_n(x)dx.
\end{equation}
The term $F_3^h$ reads
\begin{align*}
F_3^h &= -  \sum_{n=0}^{N-1} {\Delta t}\sum_{j\in
S_h} c_j^{n+1} \int_{T_j} \operatorname{div} (D({\mathbf{u}}
^{n+1}_h) \operatorname{grad} {\tilde \varphi}_n(x))dx
\\
&= -  \sum_{n=0}^{N-1} {\Delta t}\sum_{j\in S_h}
c_j^{n+1}\sum_{k\in N_j} \int_{\sigma_{jk}} D({\mathbf{u}}
^{n+1}_h) \operatorname{grad} {\tilde \varphi}_n(x) \cdot\nu_{jk}
\,d\gamma(x)   \\
&= -  \sum_{n=0}^{N-1} {\Delta t}\sum_{(j,k)
\in {\mathscr E}_h} (c_j^{n+1} -c_k^{n+1}) 
\int_{\sigma_{jk}} D({\mathbf{u}}^{n+1}_h) \operatorname{grad} {
\tilde \varphi}_n(x) \cdot\nu_{jk} \,d\gamma(x).  %\label{sec3:n31}
\end{align*}
Following \eqref{sec2:eq6disp}, we put the term $E^h_3$ in the form
\begin{align*}
E^h_3&= -  \sum_{n=0}^{N-1}{\Delta t}\sum_{j\in S_h}
\varphi_{n,j} \sum_{k\in N_j} f_{jk}^{n+1} \frac{
(c_k^{n+1}-c_j^{n+1})}{d_{jk}}   \\
&= -  \sum_{n=0}^{N-1}{\Delta t}\sum_{(j,k)\in {\mathscr E}_h} f_{jk}^{n+1} \frac{
(c_k^{n+1}-c_j^{n+1})}{d_{jk}} \;(\varphi_{n,j} -\varphi_{n,k}).
%\label{sec3:n32}
\end{align*}
We introduce $R^h_3 \equiv E^h_3 -F_3^h$; i.e.,
\[
R^h_3= \sum_{n=0}^{N-1}{\Delta t}\sum_{(j,k)
\in {\mathscr E}_h} r_{jk}^{n+1} (c_j^{n+1}-c_k^{n+1})
%\label{sec3:n33}
\]
with
\begin{align*}
r_{jk}^{n+1}
&=  \frac{f_{jk}^{n+1}}{d_{jk}} (\varphi_{n,j} -\varphi_{n,k})
- \int_{\sigma_{jk}} D({\mathbf{u}}^{n+1}_h) \operatorname{grad} {
\tilde \varphi}_n(x) \cdot\nu_{jk} \,d\gamma(x) \\
&=   \frac{\varphi_{n,j} -\varphi_{n,k}}{d_{jk}} 
\int_{\sigma_{jk}} D({\mathbf{u}}^{n+1}_h)\nu_{jk}\cdot\nu_{jk} \,d\gamma(x)\\
&\quad -\int_{\sigma_{jk}} D({\mathbf{u}}^{n+1}_h) \operatorname{grad} {
\tilde \varphi}_n(x) \cdot\nu_{jk} \,d\gamma(x)   \\
&=   \int_{\sigma_{jk}} D({\mathbf{u}}^{n+1}_h)\Big( \frac{
\varphi_{n,j} -\varphi_{n,k}}{d_{jk}} \nu_{jk} - \operatorname{grad}
{\tilde \varphi}_n(x) \Big)\cdot\nu_{jk} \,d\gamma(x).
\end{align*}
Due to the regularity of $\varphi$, $\varphi\in C^2({\overline \Omega}_T)$,
and since $ \frac{\varphi_{n,j} -\varphi_{n,k}}{d_{jk}} $
is an approximation of $\operatorname{grad} {\tilde \varphi}_n(x).\nu_{jk}$,
and since we are integrating along $\nu_{jk}$, we have
\begin{equation*}
 | r_{jk}^{n+1}| \le C  h  \int_{\sigma_{jk}} | D(
{\mathbf{u}}^{n+1}_h)\nu_{jk} |\,d\gamma(x).
\end{equation*}
We also have
\begin{align*}
&\int_{\sigma_{jk}} | D({\mathbf{u}}^{n+1}_h)\nu_{jk}
|\,d\gamma(x)\\
& = \int_{\sigma_{jk}} \frac{| D({\mathbf{u}}
^{n+1}_h)\nu_{jk}|} {(D({\mathbf{u}}^{n+1}_h)\nu_{jk}\cdot
\nu_{jk})^{1/2}}  {(D({\mathbf{u}}^{n+1}_h)\nu_{jk}\cdot\nu_{jk})^{1/2}}
\,d\gamma(x) \\
&\le  \Big(\int_{\sigma_{jk}} D({\mathbf{u}}
^{n+1}_h)\nu_{jk}\cdot\nu_{jk} \,d\gamma(x)\Big)^{1/2}
 \Big(\int_{\sigma_{jk}} \frac{| D({\mathbf{u}}
^{n+1}_h)\nu_{jk}|^2} {D({\mathbf{u}}^{n+1}_h)\nu_{jk}\cdot\nu_{jk}}
\,d\gamma(x)\Big)^{1/2} \\
& \le  \Big(\int_{\sigma_{jk}} D({\mathbf{u}}
^{n+1}_h)\nu_{jk}\cdot\nu_{jk} \,d\gamma(x)\Big)^{1/2}
 \Big(\frac{d_{l}}{d_{t}}\int_{\sigma_{jk}} | {\mathbf{u}}
^{n+1}_h | \,d\gamma(x)\Big)^{1/2}.
\end{align*}
These implies the  estimate
\begin{align*}
| r_{jk}^{n+1} (c_k^{n+1}-c_j^{n+1})|
&\le  C h \frac{ | c_k^{n+1}-c_j^{n+1}|}{d_{jk}^{1/2}} 
\Big(\int_{\sigma_{jk}} D({\mathbf{u}}^{n+1}_h)\nu_{jk}\cdot\nu_{jk}
\,d\gamma(x)\Big)^{1/2}\\
&\quad\times d_{jk}^{1/2}  \Big(\int_{\sigma_{jk}} | {\mathbf{u}
}^{n+1}_h | \,d\gamma(x)\Big)^{1/2}.
\end{align*}
By the Cauchy-Schwarz inequality, we obtain
\begin{align*}
&\sum_{(j,k)\in {\mathscr E}_h} | r_{jk}^{n+1}
(c_k^{n+1}-c_j^{n+1})| \\
&\le  C h\Big(\sum_{(j,k)\in {\mathscr
E}_h} f_{jk}^{n+1} \frac{| c_k^{n+1}-c_j^{n+1}|^2}{d_{jk}}
\Big)^{1/2}
 \Big(\sum_{(j,k)\in {\mathscr E}_h} {d_{jk}} 
\int_{\sigma_{jk}} | {\mathbf{u}}^{n+1}_h | \,d\gamma(x) \Big)^{1/2}.
\end{align*}
Similarly to \eqref{sec3:n5}, we have
\begin{equation*}
\max_{0\le n\le N-1} \sum_{(j,k)\in {\mathscr E}_h}
\int_{\sigma_{jk}} | {\mathbf{u}}^{n+1}_h | \,d\gamma(x) \le C h^{-1}
\end{equation*}
so that
\begin{align*}
| R^h_3| &\le & C  h\sum_{n=0}^{N-1}{\Delta t}
\Big(\sum_{(j,k)\in {\mathscr E}_h} f_{jk}^{n+1} \frac{|
c_k^{n+1}-c_j^{n+1}|^2}{d_{jk}}\Big)^{1/2}.
\end{align*}
According to \eqref{sec3:n22}, we obtain
\begin{equation*}
|R_3^h|\leq C h.  \label{sec3:n34}
\end{equation*}
From \eqref{sec3:n9}, \eqref{sec3:n25}--\eqref{sec3:n27},
\eqref{sec3:n29}, and \eqref{sec3:n34} we deduce that
\begin{align*}
E^h &= -\int_{\Omega _{T}}c_h(x,t)\partial _{t}\varphi
(x,t)\,dx\,dt-\int_{\Omega _{T}}c_h(x,t){\mathbf{u}}_h(x,t)\cdot
\operatorname{grad}\varphi (x,t)\,dx\,dt   \\
&\quad +\int_{\Omega _{T}}c_h(x,t)q_h^{-}(x,t)\varphi
(x,t)\,dx\,dt-\int_{\Omega _{T}}{\hat{c}}_h(x)q_h^{+}(x,t)\varphi
(x,t)\,dx\,dt   \\
&\quad +\int_{\Omega _{T}}D({\mathbf{u}}_h(x,t))\operatorname{grad}
c_h(x,t)\cdot \operatorname{grad}\varphi (x,t)\,dx\,dt-\int_{\Omega
}c_h^{0}(x)\varphi (x,0)dx  
 +\varepsilon (h)  %\label{sec3:eh}
\end{align*}
with $\lim_{h\to 0}\varepsilon (h)=0$. This order of approximation
is superior or equal to one since all the previous estimations have been done
with at least a constant times $ h$.
We summarize the result in the following statement.

\begin{lemma}\label{lemconsist}
The scheme \eqref{sec2:eq3}, \eqref{sec2:eq4},
\eqref{sec2:eqc} is consistent with the concentration equation \eqref{sec2:eqnnc}.
\end{lemma}

From the stability  proved in proposition \ref{lemcon} and this lemma,
we conclude that the numerical scheme \eqref{sec2:eq3}, \eqref{sec2:eq4},
\eqref{sec2:eqc} converge to the solution of the concentration
equation \eqref{sec2:eqnnc} and the convergence is at least of order $1$.


\section{Numerical experiments} \label{sec4}

We have applied the numerical scheme presented in Sect.~\ref{sec2} for the
simulation of a miscible flow within a horizontal reservoir of one unit
thickness over ten years period with injection and production wells at the
corners reported by  Ewing et al \cite{WEL}. In our tests we have taken
for spatial domain $\Omega =[0;400]\times [0;400]ft^2$ and $T=3600$
days. The injection well is located at the upper-right of the domain with
a volumetric injection rate of $Q=30\,$ft$^2$/day. The production
well is located at the lower-left corner with a production rate of
$Q=-30\,$ft$^2$/day.The mass flow rate $q$ in equations
\eqref{intro:eqnpres}-\eqref{intro:eqncon}
is equal to the product of the mass density $\varrho $ and the quantity $Q$ per unit
volume; i.e., $q=Q/m_j=0.024$/day for the uniform coarse spatial grid of
$ \Delta x=\Delta y=50$ft and where $m_j=\Delta x\times \Delta
y/2=1250$ft$^2$ is the uniform element $T_j$ area.
The viscosity $\mu =\mu (c)$ is taken according to  \eqref{intro:koval}
with a large adverse mobility ratio $M=\mu (0)/\mu (1)=41$ and a viscosity
of the residential fluid $\mu (0)=1.0$cp. The longitudinal and transverse
dispersion coefficients, are taken respectively, as $d_{l}=5$ft,
$d_{t}=0.5$ft. The initial concentration is $c_{0}(x,y)=0.0$ and the
injection concentration is $\hat{c}=1.0$.

The discrete equations \eqref{sec2:eq3}-\eqref{sec2:eq4} with the constrain
\eqref{sec2:eqnpmn} corresponding to the mixed finite element method
approximation of the velocity-pressure equation
\eqref{intro:eqnpres}-\eqref{intro:eqnvit},\eqref{intro:eqnnormaliz}
were solved with an augmented lagrangian method - see \cite{BF},\cite{BS}-
 which converged in very few
iterations -two or three iterations. The discrete equations
\eqref{sec2:eqc}-\eqref{sec2:cnij} corresponding to the finite volume
approximation of the concentration equation \eqref{sec2:eqnnc} were solved with a conjugate
gradient method which converged also in few iterations - less than ten
iterations.

The contour plots for the concentration of the invading fluid at time $t=3$
years, $5$ years, $7$ years and $10$ years for the data above are presented
in figures $1-4.$ They indicate the uniform fluid front move from the
injection source towards the production well and the reservoir invasion with
time. We have also observed that this reservoir invasion is closely related
to the injection and production rates. We did not take into our model the permeability and the
porosity since we are mainly concerned with the degenerate parabolic nature
of the concentration equation \eqref{sec2:eqnnc}.
That is why we did not treat the other examples reported in \cite {CD,WEL,SY},
but at the difference of Chainais and Droniou \cite{CD} who have only
observed that their
approximate concentration remains in [0;1]( and such is the case for other
proposed numerical methods; e.g., \cite {WEL,SY}), we have rigourously
proved the boundedness of our approximate concentration, more precisely
 $ 0\le c_h(x,t)\le 1$  a.e. in $\Omega_T $
under the grid conditions \eqref{sec2:eqhydt} and \eqref{sec2:hyptjk2}.
It is possible to include the permeability in our model and still having
the uniform estimates for the approximate concentration since the permeability tensor
is only used in the determination of the velocity fluxes ${u}^{n+1}_{j,k}$
on which we did not make any restrictive assumption.
The time step restriction\eqref{sec2:cfl}is only used in the proof of the grid
conditions \eqref{sec2:eqhydt} and \eqref{sec2:hyptjk2}and is easily satisfied
with the large steps $h$ and $\Delta t$ typical in petroleum engineering -
we used $h=50$ft, $\Delta t=3.6$day
in our reported numerical test.

Our grid conditions \eqref{sec2:eqhydt} and \eqref{sec2:hyptjk2} depends
of the injections rates but are easily satisfied even for very large time
steps and spatial steps as long as the injections rates are not very high.

We have also observed in the tests corresponding to different data that our
numerical scheme \eqref{sec2:eq3}-\eqref{sec2:eq4}, \eqref{sec2:eqc}
generates stable numerical approximations as long as the grid conditions
\eqref{sec2:eqhydt} and \eqref{sec2:hyptjk2} were satisfied.
This suggests that these conditions  are in fact necessary and sufficient conditions.


\begin{figure}[ht]
\begin{center}
\includegraphics[width=0.41\textwidth]{fig1a}
\includegraphics[width=0.50\textwidth]{fig1b}\\
 3 years \hfil 5 years  \\
\includegraphics[width=0.44\textwidth]{fig1c}
\includegraphics[width=0.50\textwidth]{fig1d}\\
 7 years \hfil 10 years
\end{center}
\caption{Concentration contour plots of invading fluid at 3, 5, 7 and 10
years}
\label{fig1}
\end{figure}


\section{Appendix}

Let $T$ be the triangle with vertices $A=(0,0)$; $B=(b,0)$; $C=(c_1,c_{2})$, with
$c_1>0$, $c_{2}>0$. Let $\mathbf{v}\in V_h$ and such that $\mathbf{v}$
restricted to $T$ is a polynomial of degree 1. Now assume that
$\mathbf{v}(B)=\mathbf{v}(C)=0$. Now we have by integration by parts
\[
\int_{T}\mathbf{v}\cdot\frac{\partial\mathbf{v}}{\partial x}\,dx\,dy
=-\int_{T}\frac{\partial\mathbf{v}}{\partial x}\mathbf{v}\,dx\,dy
+\int_{\partial T}\mathbf{v}\cdot\mathbf{v}\nu_{x}d\sigma
\]
where $\nu=(\nu_{x},\nu_{y})$ is the outward unit normal to $\partial T$ and
$d\sigma$ is the usual measure on the boundary.

It follows from the integral above that (by transfering to the right hand side)
\begin{equation}
2\int_{T}\mathbf{v}\cdot\frac{\partial\mathbf{v}}{\partial x}\,dx\,dy
=\int_{\partial T}\mathbf{v}\cdot\mathbf{v}\nu_{x}d\sigma\label{1}
\end{equation}
Now given that $\mathbf{v}(B)=\mathbf{v}(C)=0$ it follows that
$\int_{BC}\mathbf{v}^2\nu_{x}d\sigma=0$ similarly, since $\nu_{x}=0$ on $AB$,
 we have also $\int_{AB}\mathbf{v}^2\nu_{x}d\sigma=0$.

It is easy to see that there exist a positive constant $\gamma$ such that
\begin{equation}
\|\mathbf{v}\|_{L^2(AC)}^2\leq\gamma|\int_{AC}\mathbf{v}\cdot\mathbf{v}
\nu_{x}d\sigma|.\label{2}
\end{equation}
Combining \eqref{1} and \eqref{2} and using Cauchy - Schwartz we obtain
\begin{equation}
\|\mathbf{v}\|_{L^2(AC)}^2\leq2\gamma|\int_{T}\mathbf{v}\cdot
\frac{\partial\mathbf{v}}{\partial x}\,dx\,dy|
\leq2\gamma(\int_{T}\mathbf{v} ^2\,dx\,dy)^{1/2}
\cdot(\int_{T}(\frac{\partial\mathbf{v}}{\partial x})^2\,dx\,dy)^{1/2}\label{3}
\end{equation}
Now we simply use the fact that
\begin{equation}
\Big(\int_{T}(\frac{\partial\mathbf{v}}{\partial x})^2\,dx\,dy\Big)^{1/2}
\leq \|\ |\nabla\mathbf{v}|^2\|_{L^2(T)}\label{4}
\end{equation}
The proof obviously can be easily extended to the other sides of the triangle.
To get the estimate we want we will also use the following inverse inequality
that can be found in Ciarlet \cite[formula (3.2.33) p. 141]{CI}.
\begin{equation}
| \mathbf{v}_h| _{m,q,K}\leq C\frac{(h^n)^{\frac{1}
{q}-\frac{1}{r}}}{h^{m-l}}| \mathbf{v}_h| _{l,r,K}
\label{5}
\end{equation}
with $m=1$, $q=r=2$, $l=0 $.


\subsection*{Acknowledgement}
The author wishes to thank Professor Y. Amirat for useful conversations
and suggestions, and the anonymous referees for their
insightful comments and suggestions.

\begin{thebibliography}{00}

\bibitem{AO} B. Amaziane, M. El Ossmani;
\emph{Convergence Analysis of an Aproximation to Miscible Fluid Flows
in Porous Media by Combining Mixed Finite Element and Finite Volume Methods},
Numer. Methods Partial Differential Eq.   Vol. 24(3)(2008), pp. 799-832.

\bibitem{AZ} Y. Amirat, A. Ziani;
\emph{Asymptotic behavior of the solutions of an elliptic-parabolic system
arising in flow in porous media}, Z. Anal. Anwendungen, 23 (2004), pp. 335-351.

\bibitem{Be} J. Bear;
\emph{Dynamics of Fluids in Porous Media}, American Elsevier, New York, 1972.

\bibitem{BS} S. C. Brenner  L. R. Scott;
\emph{The Mathematical Theory of Finite Element Method, Springer Series
in Applied Mathematics}, Vol. 1, Springer, Berlin, (2002).

\bibitem{BF} F. Brezzi, M. Fortin;
\emph{ Mixed and hybrid finite element methods}, Springer Series in Computational
Mathematics Vol. 1, Springer, Berlin, (1991).

\bibitem{CD} C. Chainais-Hillairet, J. Droniou;
\emph{Convergence analysis of a mixed finite volume scheme for
 an elliptic-parabolic system modeling miscible fluid flows in porous media}.
SIAM J. Numer Anal., Vol. 45(5)(2007) 2228-2258.

\bibitem{CJ} G. Chavent, J. Jaffr\'e;
\emph{Mathematical models and
finite elements for reservoir simulation. Single phase, multiphase and
multicomponent flows through porous media}, in: Studies in Mathematics and
its Applications Vol. 17, North-Holland, (1986).

\bibitem{CZ} C. Choquet, S. Zimmermann;
\emph{Analysis of a finite-volume-finite-element scheme for a nuclear transport model},
 IMA Journal of Numerical Analysis, 31 (1), (2011), pp. 86-115.

\bibitem{DH} J. Douglas, U. Hornung (eds.);
\emph{Flow in porous media}, International Series of Numerical Mathematics,
Vol. 114, Birkhaeuser, (1993).

\bibitem{CI}  P.G. Ciarlet;
\emph{The finite element method for elliptic problems}, Studies in Mathematics and
its Applications Vol. 4, 2nd printing, North-Holland, (1979).

\bibitem{EGu}  A. Ern, J. L. Guermond;
\emph{El\'ments finis: Th\'eorie, Applications}, mise en oeuvre,
Math\'e\-matiques \& Applications (36), Springer, (2002).

\bibitem{EG} R. Eymard, T. Gallou\"et;
\emph{Convergence of a finite element-finite volume type scheme for
a system formed by an elliptic equation and a hyperbolic equation},
RAIRO, Modelisation Math. Anal. Numer. 27(7) (1993) 843-861.

\bibitem{EGH} R. Eymard, T. Gallou\"et,  R. Herbin;
\emph{Finite volume methods},  in Ph. Ciarlet and J. L. Lions, (Eds.),
Handbook for Numerical Analysis,
 Vol.7, North Holland, 2000, pp. 713-1020.

\bibitem{EHV} R. Eymard, D. Hilhorst, M. Vohralik;
\emph{A combined finite volume nonconforming/mixed-hybrid finite
element scheme for degenerate parabolic problems},
Numer. Math. 105(1) (2006).

\bibitem{JR} J. Jaffr\'e, J.-E. Roberts;
\emph{Upstream weighting and mixed
finite elements in the simulation of miscible displacements}, RAIRO,
Mod\'elisation Math. Anal. Num\'er. 19 (1985) 443-460.

\bibitem{Ko} E. J. Koval;
\emph{A Method for Predicting the Performance of
Unstable Miscible Displacements in Heterogeneous Media},  SPEJ
Trans AIME, 228 (1963) 145-154.

\bibitem {MS} F. Marpeau, M. Saad;
\emph{3D simulation of radionuclide transport in porous media},
International Journal forNumerical Methods in Fluids, 84(1) (2010) pp. 44-70.

\bibitem{PT} J. R. A. Pearson, P. M. J. Tardy;
\emph{Models for flow of non-Newtonian and complex fluids through porous media},
J. Non-Newtonian Fluid Mech. 102 (2002) 447-473.

\bibitem{RT} P.-A. Raviart, J.-M. Thomas;
\emph{A mixed finite element method for second order elliptic problems},
Math. Aspects of Finite Elem. Meth. Lect. Notes Math. 606 (1977) 292-315.

\bibitem{Sc}  A. E. Scheidegger;
\emph{The Physics of Flow through Porous Media}, Univ. Toronto Press, 1974.

\bibitem{SY}  T. Sun, Y. Yuan;
\emph{An approximation of incompressible miscible displacement
in porous media by mixed finite element method and characteristics-mixed
finite element method}, J. Computational and Applied Mathematics, 228 (2009) 391-411.

\bibitem{WEL} H. Wang, D. Liang, R. E. Ewing, S. L. Lyons, G. Quin;
\emph{An ELLAM-MFEM solution Technique for Compressible Fluid Flows in Porous Media
with Point Sources and Sinks}, Journal of Computational Physics, 159(2000) 344-376.

\bibitem{Wh} M.-F. Wheeler (Ed.);
\emph{Numerical simulation in oil recovery}, The IMA Volumes in Mathematics
and Its Applications Vol.11, Springer-Verlag, 1988.

\end{thebibliography}

\section*{Addendum posted by the editor on September 16, 2013} 

An anonymous reader informed us that the results in this article are incorrect.
Here is an extract of the reader's message and of the author's response.
\medskip 

(1) READER: 
This scheme is unfortunately based on a completely false mathematical 
approximation of the fluxes for the dispersion operator. The author
approximates $\mathbf{u}\cdot\operatorname{grad} c$, with $\mathbf{u}$ 
the Darcy velocity given by the elliptic equation and $c$ the concentration, 
along an edge of the mesh by neglecting the tangential component
of $\operatorname{grad} c$ 
along this edge. 
The author writes
$$
\mathbf{u}\cdot \operatorname{grad} c \approx (\mathbf{u}\cdot \nu)
(\operatorname{grad} c \cdot \nu), 
$$
instead of
$$
\mathbf{u}\cdot \operatorname{grad} c \approx (\mathbf{u}\cdot \nu)(\operatorname{grad} c \cdot \nu)
+ (\mathbf{u}\cdot \mathbf{t}) (\operatorname{grad} c \cdot \mathbf{t})
$$
where $\nu$ is the unit normal vector, and $\mathbf{t}$ is a tangential unit 
vector to the edge.
In general, the  term $\operatorname{grad} c \cdot \mathbf{t}$
is not negligible in front of 
$(\mathbf{u}\cdot \nu)(\operatorname{grad} c \cdot \nu)$.
\medskip 

\noindent AUTHOR:
We approximate the concentration $c(x; t_n )$ by
$c_h^n (x)$ which is piecewise constant on each triangle and therefore 
the $\operatorname{grad} c_h^n (x)$ has no tangential component,
and  the approximation
$$
\operatorname{grad} c_h\cdot \nu_{jk} =\frac{c_k-c_j}{d_{jk}}+O(h)
$$
holds. 
This approximation of the gradient may not be  precise  enough - it is
only of order one - but it is correct in our case.
\medskip 

\noindent READER:
It is unacceptable to justify this approximation by saying that 
``$c$ is a piecewise constant function on each triangle
and therefore grad $c$ has no tangential component''. 
If we follow this reasoning, then grad $c$ has no normal component either 
(since $c$ is constant in the triangle)
and it can be safely approximated by 0, thus removing any convection or diffusion
from the equation in practice.

A classical Mixed Finite Element scheme (for example the P0-P1 method)
approximates $c$ by piecewise constant functions on the triangles and
 grad $c$ by piecewise linear functions on the triangles. 
All components of the gradients must be approximated if we are 
to obtain a consistent scheme in the case of anisotropic diffusion.
It is a very active research to find proper approximations of the
gradient while not increasing the degrees of freedom (on c) too much. 
But nobody expects, because $c$ is piecewise constant,
grad $c$ to be $0$ or only along some particular (normal) directions.
It has also been shown, for a long time ago, that neglecting one component of the
gradient while approximating anisotropic diffusion leads to schemes that do
not converge to the proper solution. 
\bigskip

(2) READER:
The numerical results are completely inconsistent with
the ones found in the literature (compare with [8]). In the numerical tests 
performed by the author, the fluid simply does not invade the domain, 
it remains around the injection well. That is an unphysical behaviour which 
should have warned the authors that their scheme simply does not converge 
to an appropriate solution, because it is in fact inconsistent. 
Simply put, it is not possible to approximate a dispersive flux, 
involving a full matrix $D(u)$, by a 2-points flux.
Compare the invasion speed with test 2 in the article: 

Hong Wang, Dong Liang, Richard E. Ewings Stephen L. Lyons, Guan Qin;
An approximation to miscible fluid flows in porous
media with point sources and sinks by an
Eulerian–Lagrangian localized adjoint method and
mixed finite element methods
SIAM J. Sci. Compt. Vol. 22, No. 2, pp. 561–581
\medskip 

\noindent AUTHOR:
I did not take into account in our model the permeability
and the porosity of the domain, that is what makes the difference between our
numerical results and those mentioned. Certainly a more complete model would
be more realistic but will require a complete reprogramming. We were more
interested in proving that under  non restrictive grid conditions, our 
scheme is stable.
\medskip

End of addendum.

\end{document} 


\end{document}