\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 $$\mu \in {C}^2([0,1]), \quad 0<\mu_{-}\le\mu(c)\le\mu_{+} \quad \forall c\in (0,1), \label{intro:eqnvis}$$ where $\mu_{-}$ and $\mu_{+}$ are two fixed real numbers. The following form is widely used, see {Koval}~\cite{Ko}, $$\mu(c)= \mu(0)\big(1 + (M^{1/4} -1)c\big)^{-4}, \label{intro:koval}$$ 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: $$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}$$ 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 $$c(x,0)=c_0(x) \quad \text{for } x\in \Omega. \label{intro:eqncdtic}$$ 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., $$\int_{\Omega} p(x,t) dx=0, \quad t\in (0,T). \label{intro:eqnnormaliz}$$ 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 $$\int_\Omega ( q^{+}- q^{-})(x,t)dx =0 \quad\text{a.e. for } t\in (0,T) \label{intro:eqnsourcep}$$ while $c_{0}$ and ${\hat c}$ satisfy $${\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}$$ 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 \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} 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 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 $$\sum_{k\in S_h} m_k p^{n+1}_k = 0 \label{sec2:eqnpmn}$$ We now turn to the the discretization of the concentration equation \eqref{intro:eqncon}. First, we write \eqref{intro:eqncon} in the conservative form $$\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}$$ 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: $$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}$$ 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 $${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}$$ 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{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}$$ The semi-implicit scheme \eqref{sec2:eqc} writes also: $$\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}$$ where $$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}$$ 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 $$\theta_{j,k}^n= 1 - 2 {\frac{\Delta t} { m_k}} m_{jk} ( u^{n+1}_{j,k})^{+}, \label{sec2:hyptjk1}$$ 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 $${\Delta t}\leq C h^2,\quad \text{where C is an arbitrary positive constant,} \label{sec2:cfl}$$ 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 $$\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}$$ 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 $$\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}$$ We deduce from \eqref{sec3:n1} and \eqref{sec3:n2} that $$\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}$$ 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 $$\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}$$ and $$\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}$$ 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 $${\mathscr A}^n=\{ (j,k), \; j\in S_h,\; k\in N_j \,; \; u^n_{j,k}>0 \} . \label{sec3:n7}$$ 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 $$\operatorname{TV}(c_h) \le C h^{-1/2}. \label{sec3:n9}$$ \end{lemma} \begin{proof} Let us write \eqref{sec2:eqc} in the form $$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}$$ with $$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}$$ Using \eqref{sec2:eq3}, we transform \eqref{sec3:n12} into $$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}$$ 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 $$\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}$$ 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{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}$$ 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{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}$$ 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: $$\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}$$ 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 $$E^h \equiv E^h_1 + E^h_{2} + E^h_3 + E^h_{4}=0, \label{sec3:n25}$$ 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, $$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}$$ 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, $$| R^h_{2}| \le C h \operatorname{TV}(c_h). \label{sec3:n29}$$ Now we discuss the term $E^h_3$, proceeding as previously. We introduce $$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}$$ Note that $F_3^h$, which has a meaning, represents formally the quantity $$\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.$$ 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) $$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}$$ 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 $$\|\mathbf{v}\|_{L^2(AC)}^2\leq\gamma|\int_{AC}\mathbf{v}\cdot\mathbf{v} \nu_{x}d\sigma|.\label{2}$$ Combining \eqref{1} and \eqref{2} and using Cauchy - Schwartz we obtain $$\|\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}$$ Now we simply use the fact that $$\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}$$ 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}. $$| \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}$$ 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}