New analysis of Galerkin-mixed FEMs for incompressible miscible flow in porous media

Weiwei Sun and Chengda Wu ¹¹1Department of Mathematics, City University of Hong Kong, Kowloon, Hong Kong. The work of the authors was supported in part by a grant from the Research Grants Council of the Hong Kong Special Administrative Region, China (Project No. CityU 102613) [email protected], [email protected].

Abstract

Analysis of Galerkin-mixed FEMs for incompressible miscible flow in porous media has been investigated extensively in the last several decades. Of particular interest in practical applications is the lowest-order Galerkin-mixed method, in which a linear Lagrange FE approximation is used for the concentration and the lowest-order Raviart-Thomas FE approximation is used for the velocity/pressure. The previous works only showed the first-order accuracy of the method in $L^{2}$ -norm in spatial direction, which however is not optimal and valid only under certain extra restrictions on both time step and spatial mesh. In this paper, we provide new and optimal $L^{2}$ -norm error estimates of Galerkin-mixed FEMs for all three components in a general case. In particular, for the lowest-order Galerkin-mixed FEM, we show unconditionally the second-order accuracy in $L^{2}$ -norm for the concentration. Numerical results for both two and three-dimensional models are presented to confirm our theoretical analysis. More important is that our approach can be extended to the analysis of mixed FEMs for many strongly coupled systems to obtain optimal error estimates for all components.

1 Introduction

In many engineering applications, incompressible miscible flow in porous media can be described by the following miscible displacement system

	$\displaystyle\Phi\frac{\partial c}{\partial t}-\nabla\cdot(D({\bf u})\nabla c)+{\bf u}\cdot\nabla c{+q^{P}c=\hat{c}q^{I}},$		(1.1)
	$\displaystyle-\nabla\cdot\frac{k(x)}{\mu(c)}\nabla p=q^{I}-q^{P},$		(1.2)

with the initial and boundary conditions:

\displaystyle\begin{array}[]{ll}{\bf u}\cdot{\bf n}=0,~{}~{}D({\bf u})\nabla c\cdot{\bf n}=0&\mbox{for}~{}~{}x\in\partial\Omega,~{}~{}t\in[0,T],\\[3.0pt] c(x,0)=c_{0}(x){}{}&\mbox{for}~{}~{}x\in\Omega,\end{array}

(1.5)

where ${\bf u}$ denotes the Darcy velocity of the fluid mixture defined by

\displaystyle{\bf u}=-\frac{k(x)}{\mu(c)}\nabla p,

(1.6)

$p$ is the pressure of the fluid mixture and $c$ is the concentration. Numerical solutions of the above system (1.1)-(1.6) play a key role in these applications. Here, $k(x)$ is the permeability of the medium, $\mu(c)$ is the concentration-dependent viscosity, $\Phi$ is the porosity of the medium, $q^{I}$ and $q^{P}$ are the given injection and production sources, $\hat{c}$ is the concentration in the injection source, and $D({\bf u})=[D_{ij}({\bf u})]_{d\times d}$ is the velocity-dependent diffusion-dispersion tensor, which may be given in different forms (see [8] for details). We assume that the system is defined in a bounded smooth domain $\Omega$ in $\mathbb{R}^{d}$ $(d=2,3)$ and $t\in[0,T]$ .

In the last several decades, numerous effort has been devoted to the development of numerical methods for the system (1.1)-(1.6). Numerical simulations have been made extensively in various engineering areas, such as reservoir simulations and exploration of underground water and oil [19, 26]. Two review articles for numerical methods used in these areas have been written by Ewing and Wang [28] and Scovazzi et al.[43]. The existence of weak solutions of the system was proved by Feng [29] for the 2D model and by Chen and Ewing [13] for the 3D problem, while the existence of semi-classical/classical solutions is unknown so far. More detailed discussion can be found in [30]. The existence and uniqueness of weak solutions for some different models of Darcy flow were studied in [1, 39]. In [1], the model includes mobile and immobile species, with possibly discontinuous reaction rates, and with a variable porosity that depends on the concentration of the immobile species. The model in [39] includes a component for the unsaturated flow (the Richards equation) and another component for reactive transport, with nonlinear reaction terms. Whereas here the reaction terms are linear, the diffusion tensor depends on the fluid velocity and the viscosity depends on the solute concentration. Numerical analysis for the system (1.1)-(1.6) in two-dimensional space was first presented by Ewing and Wheeler [25] for a standard Galerkin-Galerkin approximation $({\cal C}_{h},P_{h})\in(V_{h}^{r},\widehat{V}_{h}^{s})$ to the concentration and pressure in spatial direction, where $V_{h}^{r}$ denotes $C^{0}$ Lagrange finite element space of piecewise polynomials of degree $r$ and $\widehat{V}_{h}^{s}:=V_{h}^{s}/\{constant\}$ . Later, a Galerkin-mixed method was proposed by Douglas et al. [17] for solving the system (1.1)-(1.6). In the Galerkin-mixed method, a standard Lagrange type Galerkin approximation ${\cal C}_{h}\in V_{h}^{r}$ was applied for the concentration equation and a mixed approximation in the Raviart–Thomas finite element space ( $(P_{h},{\bf U}_{h})\in S_{h}^{k}\times{\bf H}_{h}^{k}$ ) was used for the pressure equation. Due to the nature of the continuity of the velocity and the discontinuity of the gradient of the pressure in applications, the Galerkin-mixed method is more popular in many areas, particularly in industries of underground water and oil. The error estimate of the semi-discrete Galerkin-mixed method was first presented in [17] and later, in [18] for a fully discrete semi-implicit Euler scheme. In [18], the error estimate

\displaystyle\|c^{n}-{\cal C}_{h}^{n}\|_{L^{2}}+\|p^{n}-P_{h}^{n}\|_{L^{2}}+\|{\bf u}^{n}-{\bf U}_{h}^{n}\|_{L^{2}}\leq C(\tau+h_{c}^{r+1}+h_{p}^{k+1})\,

(1.7)

was established for $d=2$ under the time step restriction $\tau=o(h)$ and an extra spatial mesh size condition,

\displaystyle h_{c}^{-1}h_{p}^{k+1}=o(1)

(1.8)

where $h_{c}$ and $h_{p}$ denote the mesh sizes of FE discretization for the concentration and pressure equations, respectively. Subsequently, some improvements on time step restriction and spatial mesh condition were presented by several authors [12, 21, 34, 35]. In particular, the error estimate (1.7) was proved in [34] for $d=2,3$ , $h_{p}=h_{c}=h$ and $k\geq 1$ in which no time step restriction was required. Based on superconvergence analysis, further improvement on the spatial mesh condition (1.8) was presented in [15], while the analysis is valid only for regular rectangular meshes.

The most commonly-used Galerkin-mixed method in practical computation is the lowest order one ( $r=1,k=0$ ) [12, 15, 18, 21, 28, 43], a linear approximation to the concentration and the lowest order Raviart–Thomas approximation to the pressure and velocity. The lowest order Galerkin-mixed method has been widely used in a variety of numerical simulations and applications, $e.g.$ , see [18, 26, 43]. In this case, the error estimate (1.7) reduces to

\displaystyle\|c^{n}-{\cal C}_{h}^{n}\|_{L^{2}}+\|p^{n}-P_{h}^{n}\|_{L^{2}}+\|{\bf u}^{n}-{\bf U}_{h}^{n}\|_{L^{2}}\leq C(\tau+h_{p}+h_{c}^{2})

(1.9)

and the spatial mesh condition (1.8) becomes

\displaystyle h_{c}^{-1}h_{p}=o(1)\,.

(1.10)

There are two serious concerns arising from previous analysis for the popular lowest-order Galerkin-mixed FEM. First, it is noted that the error estimate (1.9) is not optimal for the concentration in $L^{2}$ -norm. In previous analysis of the lowest-order Galerkin-mixed method, a linear approximation to the concentration only produces the numerical concentration of the accuracy $O(h)$ in spatial direction, while the optimal accuracy of a linear approximation is $O(h^{2})$ in the traditional sense. Due to the strong coupling of the system, it was assumed that the one-order lower accuracy of the numerical pressure/velocity may pollute the numerical concentration through the diffusion-dispersion tensor $D({\bf u})$ and the viscosity $\mu=\mu(c)$ . Our numerical results show that this assumption is incorrect. Secondly, based on the above spatial mesh condition (1.10), one has to use two types of spatial meshes with a much finer one for the pressure/velocity equation. The Galerkin-mixed method based on the same mesh for both concentration and pressure equations ( $h_{p}=h_{c}$ ) may not satisfy the spatial mesh condition (1.10) although this method is more efficient and most commonly-used in practical computation.

Moreover, mixed finite element methods have been used for solving the system (1.1)-(1.6) by combining with many different schemes in time direction and different approximations to the concentration equation, such as characteristics type mixed method [4, 21, 27, 31, 33, 42, 47], finite volume method [2, 32], ELLAM [45, 46] and SUPG [37]. However, the non-optimality of error estimates for the concentration and a time-step/spatial-mesh condition as mentioned above arise again in the analysis of these methods. Optimal error estimates mentioned above are usually based on strong regularity assumptions on data and solution. A related topic is the convergence of numerical schemes under low regularity assumptions. The convergence of a semidiscrete FEM was proved in [16] for a linear parabolic equation under a weak regularity assumption of the solution in $L^{2}(0,T;H_{0}^{1}(\Omega)\cap H^{1}(0,T;H^{-1}(\Omega))$ . An implicit Euler scheme with a mixed-DG approximation in spatial direction was proposed in [7] for the nonlinear system (1.1)-(1.6). The convergence of the discrete solution to certain weak solution of the system was proved in [7] by applying the Aubin-Lions compactness on a nonconforming space. Since only some weak regularity assumptions on data were made in [39], their analysis implies the existence of weak solution of the nonlinear system under weaker assumptions than those in [13, 29]. Later, a high-order DG scheme in time direction was studied in [41]. However, no optimal convergence rate/error estimate was obtained under these weak assumptions. On the other hand, degenerate cases were analyzed, $e.g.$ , in [4, 14]. Usually, the convergence rate for degenerate system is one-order lower. In particular, in [4], a volume corrected characteristics-mixed method is proposed for a purely transport problem. A lower-order $L^{1}$ -norm error estimate $O(h/\sqrt{\tau}+h+\tau^{m})$ is obtained, where $m$ is related to the accuracy of the characteristic tracing. Different models of coupled Darcy flow and reactive transport are studied, for example, in [1, 6, 39]. These include more complicated, nonlinear reaction terms, but the system is weakly coupled, since the diffusion-dispersion tensor $D({{\bf u}})=I$ and/or the viscosity is constant. An Euler implicit-mixed finite element scheme is analyzed in [39], in which the lowest order mixed FE approximation is used for both the concentration equation and pressure equation. The optimal first-order accuracy is established for the weakly coupled model. Numerical methods and analysis for incompressible and immiscible Darcy flow can be found, $e.g.$ , in [14, 26, 38].

The main purpose of this paper is to establish the optimal error estimate of Galerkin-mixed methods for all three components, concentration, pressure and velocity, without the time-step restriction and spatial mesh condition. In particular, for the lowest order Galerkin-mixed method $(r=1,k=0)$ , we will provide the optimal error estimate

\displaystyle\|c^{n}-{\cal C}^{n}_{h}\|_{L^{2}}+h(\|p^{n}-P^{n}_{h}\|_{L^{2}}+\|{\bf u}^{n}-{\bf U}_{h}^{n})\|_{L^{2}})\leq C(\tau+h^{2})

(1.11)

for $h=h_{p}=h_{c}$ unconditionally. The analysis is based on an elliptic quasi-projection. In terms of the projection and a negative norm estimate of Raviart–Thomas finite element methods for the pressure equation, the low order accuracy of the velocity will not pollute the concentration in our analysis and the lowest order Galerkin-mixed method provides numerical concentration of the accuracy $O(h^{2})$ in $L^{2}$ -norm. Also we extend our analysis to the general approximation ( $(C_{h}^{n},P_{h}^{n},{\bf U}_{h}^{n})\in V_{h}^{r}\times S_{h}^{r-1}\times{\bf H}_{h}^{r-1}$ ) to obtain the optimal error estimate

\displaystyle\|c^{n}-C_{h}^{n}\|_{L^{2}}+h(\|p^{n}-P_{h}^{n}\|_{L^{2}}+\|{\bf u}^{n}-{\bf U}_{h}^{n}\|_{L^{2}})\leq C(\tau+h^{r+1})\,.

(1.12)

With the numerical concentration $C_{h}^{n}$ , a new numerical velocity of the same order accuracy as $C_{h}^{n}$ can be calculated by resolving the (elliptic) pressure equation in a given time level with a higher-order approximation. More important is that such a strong coupling can be found in many other physical systems, $e.g.$ , see [3, 23, 50], where a higher-order approximation was also used for one of the computational components. Our approach can be extended to finite element analysis for these strongly coupled systems to obtain optimal error estimates for all components.

The rest of the paper is organized as follows. In Section 2, we introduce a linearized Euler scheme with Galerkin-mixed approximations in the spatial direction for the system (1.1)-(1.6) and present our main results. In Section 3, we introduce a new elliptic quasi-projection and establish the boundedness of numerical solutions in terms of an error splitting technique proposed in [34]. Then we prove the optimal error estimates of Galerkin-mixed FEMs in $L^{2}$ -norm unconditionally. In Section 4, we establish some basic estimates of the quasi-projection which were used in the proof of the main theorem. Finally, numerical simulations for the system in both two and three dimensional spaces are provided in Section 5. Numerical results confirm our theoretical analysis that the methods provide the optimal accuracy in $L^{2}$ -norm for all three physical components.

2 The Galerkin-mixed FEM and the main results

2.1 Notations and assumptions

For any integer $m\geq 0$ and $1\leq p\leq\infty$ , let $W^{m,p}$ be the usual Sobolev spaces and $H^{m}:=W^{m,2}$ . In addition, we denote by ${\bf H(div)}$ the space of vector-valued functions ${\vec{f}}\in[L^{2}(\Omega)]^{d}$ such that $\nabla\cdot{\vec{f}}\in L^{2}(\Omega)$ and ${\vec{f}}\cdot\nu=0$ on $\partial\Omega$ . $L^{k}([0;T];W^{m,p}(\Omega))$ denotes the space of time-dependent functions valued in $W^{m,p}(\Omega)$ , which are $L^{k}$ integrable $w.r.t.$ time in the sense of Bochner, while $H^{1}([0,T];W^{m,n}(\Omega))$ denotes the space of time-dependent functions valued in $W^{m,n}(\Omega)$ , which are $H^{1}$ integrable $w.r.t.$ time in the sense of Bochner.

Let $\mathcal{T}_{h}$ be a regular triangular partition of $\Omega$ with $\Omega=\cup_{K}\Omega_{K}$ and the mesh size $h=\max_{\Omega_{K}\in\mathcal{T}_{h}}\{\mathrm{diam}\,\Omega_{K}\}$ . For a given division of $\mathcal{T}_{h}$ , we define the classical Lagrange finite element spaces by

\displaystyle V_{h}^{r}=\{v_{h}\in C^{0}(\Omega):v_{h}|_{K}\in P_{r}(K),\quad\forall K\in\mathcal{T}_{h}\},

and Raviart-Thomas finite element spaces [40, 44] by

	$\displaystyle{\bf H}_{h}^{s}:=\{{\bf v}_{h}\in{\bf H(div)}:{\bf v}_{h}\|_{K}\in[P_{s}(K)]^{d}+{\bf x}P_{s}(K),\quad\forall K\in\mathcal{T}_{h}\}$
	$\displaystyle S_{h}^{s}:=\{v_{h}\in L^{2}:v_{h}\|_{K}\in P_{s}(K),\quad\forall K\in\mathcal{T}_{h}\},\quad\widehat{S}_{h}^{s}:=S_{h}^{s}/\{constants\},$

where $P_{r}(K)$ is the space of polynomials of degree $r\geq 0$ on $K$ . Moreover, we denote by $I_{h}$ the commonly used Lagrange nodal interpolation operator on $V_{h}^{r}$ .

Let $\{t_{n}\}_{n=0}^{n}$ be a uniform partition in the time direction with the step size $\tau=T/N$ and we denote

p^{n}=p(x,t_{n}),\quad{\bf u}^{n}={\bf u}(x,t_{n}),\quad c^{n}=c(x,t_{n})\,.

For any sequence of functions $\{f^{n}\}_{n=0}^{n}$ , we define

D_{t}f^{n+1}=\frac{f^{n+1}-f^{n}}{\tau}\,.

Throughout this article, we make use of the following assumptions:

(A1)

The solution to the initial-boundary value problem (1.1)-(1.6) exists and satisfies

		$\displaystyle\\|p\\|_{L^{\infty}([0,T];W^{2,4}\cap H^{s+1})}+\\|{\bf u}\\|_{L^{\infty}([0,T];W^{1,4}\cap H^{s+1})}+\\|{\bf u}_{t}\\|_{L^{2}([0,T];H^{s+1})}+\\|c\\|_{L^{\infty}([0,T];W^{2,4}\cap H^{r+1})}$
		$\displaystyle+\\|c_{t}\\|_{L^{\infty}([0,T];H^{2})}+\\|c_{t}\\|_{L^{4}([0,T];W^{1,4}\cap H^{r+1})}+\\|c_{tt}\\|_{L^{4}([0,T];L^{4})}\leq K_{1}.$		(2.1)

(A2)

${\|\hat{c}\|_{H^{1}([0,T];H^{1}(\Omega))},\,\|q^{I}\|_{H^{1}([0,T];H^{1}(\Omega))},\,\|q^{P}\|_{H^{1}([0,T];H^{1}(\Omega))}}\leq K_{2}$ .
(A3)

$k\in W^{2,\infty}(\Omega)$ ; $k_{0}^{-1}\leq k(x)\leq k_{0},\forall~{}x\in\Omega.$
(A4)

$\mu\in C^{1}(\mathbb{R})$ ; $\exists~{}\mu_{0}>0,~{}\mu_{0}^{-1}\leq\mu(s)\leq\mu_{0},~{}\forall~{}s\in\mathbb{R}.$
(A5)

Following [8], the diffusion-dispersion tensor is defined by

$\displaystyle D({{\bf u}})=\Phi\left(d_{mt}(|{{\bf u}}|)I+d_{lt}(|{{\bf u}}|){{\bf u}}\otimes{{\bf u}}\right),$ (2.2)

where $d_{mt}(z)>d_{m}>0$ , $d_{lt}(z)>0$ for $z>0$ and ${\bf u}\otimes{\bf u}={\bf u}{\bf u}^{T}$ denotes a $d\times d$ matrix. For simplicity, here we further assume that $d_{mt},d_{lt}\in H^{3}(\mathbb{R})$ and $\partial_{t}\nabla D({\bf v})\in L^{\infty}(0,T;L^{2}(\Omega))$ for all smooth function $\bf v$ as required in [48]. More discussion on the regularity of the diffusion-dispersion tensor was given in [11, 35], where they get optimal $L^{p}$ error estimates under the regularity assumption of the commonly-used Bear–Scheidegger diffusion-dispersion tensor $D({\bf u})\in W^{1,\infty}(\Omega\times(0,T))$ and a similar assumption to the exact solution as in ((A1)).
(A6)

To keep the well-posedness of the initial-boundary value problem (1.1)-(1.6), we require

$\int_{\Omega}q^{I}\,{\rm d}x=\int_{\Omega}q^{P}\,{\rm d}x.$ (2.3)

2.2 Schemes and main results

Before proposing the fully discrete numerical scheme, we introduce the weak formulation of the deeply coupled system (1.1)-(1.6). Find $p\in L^{2}(0,T;L^{2}(\Omega)/\{constants\})$ , ${\bf u}\in L^{2}(0,T;{\bf H(div)})$ and $c\in H^{1}(0,T;H^{1}(\Omega))$ , such that for all ${\bf v}\in{\bf H(div)}$ , $\varphi\in L^{2}(\Omega)$ and $\phi\in H^{1}(\Omega)$ ,

	$\displaystyle\biggl{(}\frac{\mu(c)}{k(x)}{\bf u},\,{\bf v}\biggl{)}=-(p,\,\nabla\cdot{\bf v}),$		(2.4)
	$\displaystyle(\nabla\cdot{\bf u},\,\varphi)=(q^{I}-q^{P},\,\varphi),$		(2.5)
	$\displaystyle(\phi\partial_{t}c,\,\phi)+(D({\bf u})\nabla c,\,\nabla\phi)+({\bf u}\cdot\nabla c,\,\phi){+(q^{P}c,\phi)=(\hat{c}q^{I}},\,\phi)$		(2.6)

for $a.e.$ $t\in(0,T]$ , where the initial concentration is given by $c(x,0)=c_{0}(x)$ .

With the above notations, a fully discrete Galerkin-mixed finite element scheme is to find $P_{h}^{n}\in\widehat{S}_{h}^{s}$ , ${\bf U}_{h}^{n}\in{\bf H}_{h}^{s}$ and ${\cal C}_{h}^{n}\in V_{h}^{r}$ , $n=0,1,\cdots,N$ , such that for all ${\bf v}_{h}\in{\bf H}_{h}^{s}$ , $\varphi_{h}\in S_{h}^{s}$ and $\phi_{h}\in V_{h}^{r}$ ,

	$\displaystyle\biggl{(}\frac{\mu({\cal C}_{h}^{n})}{k(x)}{\bf U}_{h}^{n+1},\,{\bf v}_{h}\biggl{)}=-\Big{(}P_{h}^{n+1},\,\nabla\cdot{\bf v}_{h}\Big{)},$		(2.7)
	$\displaystyle\Big{(}\nabla\cdot{\bf U}_{h}^{n+1},\,\varphi_{h}\Big{)}=\Big{(}q^{I}-q^{P},\,\varphi_{h}\Big{)},$		(2.8)
	$\displaystyle\Big{(}\Phi D_{t}{\cal C}_{h}^{n+1},\,\phi_{h}\Big{)}+\Big{(}D({\bf U}_{h}^{n+1})\nabla{\cal C}_{h}^{n+1},\,\nabla\phi_{h}\Big{)}$
	$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}+\Big{(}{\bf U}_{h}^{n+1}\cdot\nabla{\cal C}_{h}^{n},\,\phi_{h}\Big{)}+{\Big{(}q^{P}{\cal C}_{h}^{n+1},\,\phi_{h}\Big{)}}=\Big{(}\hat{c}q^{I},\,\phi_{h}\Big{)},$		(2.9)

where the initial data ${\cal C}_{h}^{0}=I_{h}c_{0}$ . Some slightly different schemes were investigated by several authors. In [12, 18], a scheme with two different partitions for (2.7) and (2.8)-(2.9) was investigated, while a smaller mesh size was suggested for the pressure/velocity than for the concentration for the lowest-order method. In [39], an extra Lipschitz continuous reaction term was introduced for some applications and the lowest-order mixed FEM was used for both concentration equation and pressure equation. A Crank-Nicolson scheme with the coupled convection term $({\bf U}_{h}^{n+1}\cdot\nabla{\cal C}_{h}^{n+1},\phi_{h})$ was proposed in [10]. Moreover, a fully implicit scheme was studied in [24, 42], where an extra inner iteration was required at each time step for solving a system of nonlinear equations. In this paper, we only focus our attention to the scheme (2.7)-(2.9), while the analysis for schemes mentioned above is similar.

In this paper, we denote by $C$ a generic positive constant and by $\epsilon$ a generic small positive constant, which are independent of $n$ , $h$ and $\tau$ , and may depend upon $K_{1}$ , $K_{2}$ and the physical constants $k_{0}$ and $\mu_{0}$ . The following classical Gagliardo-Nirenburg interpolation inequality will be often used in our proof,

\displaystyle\|\partial^{j}u\|_{L^{p}}\leq C\|\partial^{m}u\|_{L^{k}}^{a}\,\|u\|_{L^{q}}^{1-a}+C\|u\|_{L^{q}},

(2.10)

for $0\leq j<m$ and $\frac{j}{m}\leq a\leq 1$ with $\frac{1}{p}=\frac{j}{d}+a\left(\frac{1}{k}-\frac{m}{d}\right)+(1-a)\frac{1}{q}$ , except $1<k<\infty$ and $m-j-\frac{n}{k}$ is a non-negative integer, in which case the above estimate holds only for $\frac{j}{m}\leq a<1$ .

We present our main results in the following theorem.

Theorem 2.1

Suppose that the initial-boundary value problem (1.1)-(1.6) under the assumptions $(A2)-(A6)$ has a unique solution $(p,{\bf u},c)$ which satisfies ((A1)) with $s=r-1$ . Then there exist positive constants $h_{0}$ and $\tau_{0}$ such that when $h<h_{0}$ and $\tau<\tau_{0}$ , the finite element system (2.7)-(2.9) admits a unique solution $(P_{h}^{n},{\bf U}^{n}_{h},{\cal C}_{h}^{n})\in(\widehat{S}_{h}^{r-1},{\bf H}_{h}^{r-1},V_{h}^{r})$ , $n=1,\cdots,N$ , satisfying

	$\displaystyle\max_{1\leq n\leq N}\\|{\cal C}_{h}^{n}-c^{n}\\|_{L^{2}}\leq C_{0}(\tau+h^{r+1}),$		(2.11)
	$\displaystyle\max_{1\leq n\leq N}\\|P_{h}^{n}-p^{n}\\|_{L^{2}}+\max_{1\leq n\leq N}\\|{\bf U}^{n}_{h}-{\bf u}^{n}\\|_{L^{2}}\leq C_{0}(\tau+h^{r}),\quad\mbox{ for }r\geq 1,$		(2.12)

where $C_{0}$ is a constant, independent of $h$ , $\tau$ and $n$ , and may be dependent on $K_{1}$ , $K_{2}$ , $k_{0}$ and $\mu_{0}$ .

We will present the proof of Theorem 2.1 in the next two sections.

3 Analysis

The key to our analysis is a new elliptic quasi-projection. In this section, we introduce the projection and prove our main results in Theorem 2.1 in terms of an error splitting technique proposed in [34]. Correspondingly to the fully discrete system (2.7)-(2.9), we define the time-discrete solution $(P^{n+1},{\bf U}^{n+1},{\cal C}^{n+1})$ by the following elliptic system:

	$\displaystyle{\bf U}^{n+1}=-\frac{k(x)}{\mu({\cal C}^{n})}\nabla P^{n+1},$		(3.1)
	$\displaystyle\nabla\cdot{\bf U}^{n+1}=q^{I}-q^{P},$		(3.2)
	$\displaystyle\Phi D_{t}{\cal C}^{n+1}-\nabla\cdot(D({\bf U}^{n+1})\nabla{\cal C}^{n+1})+{\bf U}^{n+1}\cdot\nabla{\cal C}^{n}{+q^{P}{\cal C}^{n+1}}=\hat{c}q^{I},$		(3.3)

for $x\in\Omega$ and $t\in[0,T]$ , with the initial and boundary conditions

\displaystyle\begin{array}[]{ll}{\bf U}^{n+1}\cdot{\bf n}=0,~{}~{}D({\bf U}^{n+1})\nabla{\cal C}^{n+1}\cdot{\bf n}=0&\mbox{for}~{}~{}x\in\partial\Omega,~{}~{}t\in[0,T],\\[3.0pt] {\cal C}^{0}(x)=c_{0}(x){}{}&\mbox{for}~{}~{}x\in\Omega,\end{array}

(3.6)

The condition $\int_{\Omega}P^{n+1}dx=0$ is enforced for the uniqueness of solution. The fully discrete FE solution $(P_{h}^{n+1},{\bf U}_{h}^{n+1},{\cal C}_{h}^{n+1})$ can be viewed as a FE solution of the time-discrete system (3.1)-(3.6).

3.1 Preliminary

Before to prove our main results, we present some lemmas in this section. With the solution $(P^{n},{\bf U}^{n},{\cal C}^{n})$ to the time-discrete system, the error functions can be split into

	$\displaystyle\\|{\bf U}^{n}_{h}-{\bf u}^{n}\\|_{L^{2}}\leq\\|{\bf U}^{n}-{\bf U}_{h}^{n}\\|_{L^{2}}+\\|{\bf U}^{n}-{\bf u}^{n}\\|_{L^{2}},$		(3.7)
	$\displaystyle\\|P_{h}^{n}-p^{n}\\|_{L^{2}}\leq\\|P^{n}-P_{h}^{n}\\|_{L^{2}}+\\|P^{n}-p^{n}\\|_{L^{2}},$		(3.8)
	$\displaystyle\\|{\cal C}_{h}^{n}-c^{n}\\|_{L^{2}}\leq\\|{\cal C}^{n}-{\cal C}_{h}^{n}\\|_{L^{2}}+\\|{\cal C}^{n}-c^{n}\\|_{L^{2}}.$		(3.9)

The estimates for the second parts of the above splitting and the regularity of the solution of the time-discrete system (3.1)-(3.6) were given in Theorem 3.1 of [34] under a slightly different assumption. We present these results in the following lemma and the proof is omitted.

Lemma 3.1

Suppose that the initial-boundary value problem (1.1)-(1.6) has a unique solution $(p,{\bf u},c)$ which satisfies ((A1)). Then there exists a positive constant $\tau_{0}^{*}$ such that when $\tau<\tau_{0}^{*}$ , the time-discrete system (3.1)-(3.6) admits a unique solution $(P^{n},{\bf U}^{n},{\cal C}^{n})$ , $n=1,\cdots,N$ , satisfying

	$\displaystyle\\|P^{n}\\|_{W^{1,4}}+\\|{\bf U}^{n}\\|_{W^{1,4}}+\\|{\cal C}^{n}\\|_{W^{2,4}}+\\|D_{t}{\cal C}^{n}\\|_{L^{4}}+\\|\nabla{\cal C}^{n}\\|_{L^{\infty}}$
	$\displaystyle+\biggl{(}\sum_{n=1}^{N}\tau\\|D_{t}U^{n}\\|_{W^{1,4}}^{2}\biggl{)}^{\frac{1}{2}}+\biggl{(}\sum_{n=1}^{N}\tau\\|D_{t}C^{n}\\|_{H^{2}}^{2}\biggl{)}^{\frac{1}{2}}\leq C_{1},$		(3.10)

and

\displaystyle\max_{1\leq n\leq N}\|P^{n}-p^{n}\|_{L^{4}}+\max_{1\leq n\leq N}\|{\bf U}^{n}-{\bf u}^{n}\|_{L^{4}}+\max_{1\leq n\leq N}\|{\cal C}^{n}-c^{n}\|_{L^{4}}\leq C_{1}\tau.

(3.11)

where $C_{1}$ is a constant independent of $h$ , $\tau$ , $n$ and $C_{0}$ in Theorem 2.1 and may be dependent on $K_{1}$ , $K_{2}$ , $k_{0}$ and $\mu_{0}$ .

Now we introduce our elliptic quasi-projection. For any fixed integer $n\geq 1$ , we denote by $(\widetilde{P}^{n}_{h},\widetilde{\bf U}_{h}^{n})$ the mixed projection of $(P^{n},{\bf U}^{n})$ on $\widehat{S}_{h}^{r-1}\times{\bf H}_{h}^{r-1}$ such that

	$\displaystyle\biggl{(}\frac{\mu({\cal C}^{n})}{k(x)}{\widetilde{\bf U}}_{h}^{n+1},\,{\bf v}_{h}\biggl{)}=-\Big{(}{\widetilde{P}}_{h}^{n+1},\,\nabla\cdot{\bf v}_{h}\Big{)},$		(3.12)
	$\displaystyle\Big{(}\nabla\cdot({\widetilde{\bf U}}_{h}^{n+1}-{\bf U}^{n+1}),\,\varphi_{h}\Big{)}=0,\quad\forall(\varphi_{h},{\bf v}_{h})\in S_{h}^{r-1}\times{\bf H}_{h}^{r-1}\,.$		(3.13)

By the classical mixed FE theory [9, 22, 40, 44] and negative norm estimates in [20], we have

	$\displaystyle\\|{\bf U}^{n}-{\widetilde{\bf U}}_{h}^{n}\\|_{L^{p}}+\\|P^{n}-{\widetilde{P}}_{h}^{n}\\|_{L^{p}}\leq Ch,\quad\mbox{for~{}all}~{}~{}2\leq p\leq 4,$		(3.14)
	$\displaystyle\\|{\bf U}^{n}-{\widetilde{\bf U}}_{h}^{n}\\|_{H^{-1}}+\\|P^{n}-{\widetilde{P}}_{h}^{n}\\|_{H^{-1}}\leq Ch^{2},$		(3.15)

For a given ${\bf U}^{n}$ , the quasi-projection $\Pi_{c}^{n}:H^{1}(\Omega)\rightarrow V_{h}^{r}$ is defined by the elliptic problem,

	$\displaystyle\Big{(}D({\bf U}^{n})\nabla(\Pi_{c}^{n}{\cal C}^{n}-{\cal C}^{n}),\,\nabla\phi_{h}\Big{)}+\Big{(}(D({\widetilde{\bf U}}_{h}^{n})-D({\bf U}^{n}))\nabla{\Pi_{c}^{n}\cal C}^{n},\,\nabla\phi_{h}\Big{)}=0,$
	$\displaystyle\quad\mbox{for~{}all}~{}~{}\phi_{h}\in V_{h}^{r},~{}~{}n\geq 1,$		(3.16)

with $\int_{\Omega}(\Pi_{c}^{n}{\cal C}^{n}-{\cal C}^{n})dx=0$ and $\Pi_{c}^{0}:=I_{h}$ . Clearly, $\Pi_{c}^{n}$ is not a projection since $\Pi_{c}^{n}{\cal C}_{h}^{n}\neq{\cal C}_{h}^{n}$ , and reduces to a classical elliptic projection only when ${\bf U}_{p}^{n}={\bf U}^{n}$ . We present some basic estimates of the elliptic quasi-projection $\Pi_{c}^{n}$ in the following lemma and the proof will be given in section 4.

Lemma 3.2

Under the assumptions of Theorem 2.1, there exists $h_{1}>0$ such that for any $h\leq h_{1}$ and $2\leq p\leq 4$

	$\displaystyle\\|{\cal C}^{n}-\Pi_{c}^{n}{\cal C}^{n}\\|_{L^{2}}+h\\|\nabla({\cal C}^{n}-\Pi_{c}^{n}{\cal C}^{n})\\|_{L^{p}}\leq C_{2}h^{2},$		(3.17)
	$\displaystyle\\|\Pi_{c}^{n}{\cal C}^{n}\\|_{W^{1,\infty}}\leq C_{2},$		(3.18)

and

\displaystyle\biggl{(}\sum_{n=0}^{N-1}\tau\|D_{t}({\cal C}^{n}-\Pi_{c}^{n}{\cal C}^{n})\|_{L^{2}}^{2}\biggl{)}^{1/2}\leq C_{2}h^{2}.

(3.19)

where $C_{2}$ is a constant independent of $h$ , $\tau$ , $n$ , $C_{0}$ and may be dependent upon $K_{1}$ , $K_{2}$ , $C_{1}$ , $k_{0}$ and $\mu_{0}$ .

Following the splitting in (3.7)-(3.9), we first present estimates for the first parts of the splitting in the following lemma.

Lemma 3.3

Under the assumptions of Theorem 2.1, there exist positive constants $\widehat{h}_{0}$ and $\widehat{\tau}_{0}$ such that when $h<\widehat{h}_{0}$ and $\tau<\widehat{\tau}_{0}$ , the finite element system (2.7)-(2.9) admits a unique solution $\{P_{h}^{n},{\bf U}_{h}^{n},{\cal C}_{h}^{n}\}_{n=1}^{n}\in(\widehat{S}_{h}^{r-1},{\bf H}_{h}^{r-1},V_{h}^{r})$ , which satisfies

\displaystyle\|{\cal C}_{h}^{n}-{\cal C}^{n}\|_{L^{2}}+h(\|P_{h}^{n}-P^{n}\|_{L^{2}}+\|{\bf U}^{n}_{h}-{\bf U}^{n}\|_{L^{2}})\leq C_{3}h^{2},\quad\mbox{ for }r\geq 1

(3.20)

where $C_{3}$ is a constant independent of $h$ , $\tau$ , $n$ , $C_{0}$ and may be dependent upon $K_{1}$ , $K_{2}$ , $C_{1}$ , $C_{2}$ , $k_{0}$ and $\mu_{0}$ .

Proof Since at each time step (2.7)-(2.8) is a standard saddle point system and the coefficient matrix of the FE system (2.9) is symmetric positive definite, the existence and uniqueness of the numerical solution follow immediately.

Let

\theta_{p}^{n}=P_{h}^{n}-{\widetilde{P}}_{h}^{n},\quad\theta_{\bf u}^{n}={\bf U}^{n}_{h}-{\widetilde{\bf U}}_{h}^{n}\quad\mbox{ and }\quad\theta_{c}^{n}={\cal C}_{h}^{n}-\Pi_{c}^{n}{\cal C}^{n}.

By noting the projection error estimates in (3.14) and (3.17), we only need to prove the following estimate

\displaystyle\|\theta_{p}^{n}\|_{L^{2}}+\|\theta_{\bf u}^{n}\|_{L^{2}}+\|\theta_{c}^{n}\|_{L^{2}}\leq Ch^{2}.

(3.21)

Since the solution of the time-discrete system (3.1)-(3.6) satisfies

	$\displaystyle\biggl{(}\frac{\mu({\cal C}^{n})}{k(x)}{\bf U}^{n+1},\,v_{h}\biggl{)}=-\Big{(}P^{n+1},\,\nabla\cdot{\bf v}_{h}\Big{)},$		(3.22)
	$\displaystyle\Big{(}\nabla\cdot{\bf U}^{n+1},\,\varphi_{h}\Big{)}=\Big{(}q^{I}-q^{P},\,\varphi_{h}\Big{)},$		(3.23)
	$\displaystyle\Big{(}\Phi D_{t}{\cal C}^{n+1},\,\phi_{h}\Big{)}+\Big{(}D({\bf U}^{n+1})\nabla{\cal C}^{n+1},\,\nabla\phi_{h}\Big{)}$
	$\displaystyle~{}~{}~{}+\Big{(}{\bf U}^{n+1}\cdot\nabla{\cal C}^{n},\,\phi_{h}\Big{)}+{\Big{(}q^{P}{\cal C}^{n+1},\,\phi_{h}\Big{)}}=\Big{(}\hat{c}q^{I},\,\phi_{h}\Big{)},$		(3.24)

for any ${\bf v}_{h}\in{\bf H}_{h}^{r-1}$ , $\varphi_{h}\in S_{h}^{r-1}$ and $\phi_{h}\in V_{h}^{r}$ , from the finite element system (2.7)-(2.9), we can see that the error functions $(\theta_{p}^{n},\theta_{\bf u}^{n},\theta_{c}^{n})$ satisfy following system:

	$\displaystyle\biggl{(}\frac{\mu({\cal C}^{n}_{h})}{k(x)}\theta_{\bf u}^{n+1}+(\frac{\mu({\cal C}^{n}_{h})}{k(x)}-\frac{\mu({\cal C}^{n})}{k(x)}){\widetilde{\bf U}}_{h}^{n+1},{\bf v}_{h}\biggl{)}=-\Big{(}\theta_{p}^{n+1},\,\nabla\cdot{\bf v}_{h}\Big{)},$	(3.25)
	$\displaystyle\Big{(}\nabla\cdot\theta_{\bf u}^{n+1},\,\varphi_{h}\Big{)}=0,$	(3.26)
	$\displaystyle\Big{(}\Phi D_{t}\theta_{c}^{n+1},\,\phi_{h}\Big{)}+\Big{(}D({\bf U}_{h}^{n+1})\nabla\theta_{c}^{n+1},\,\nabla\phi_{h}\Big{)}{+\Big{(}q^{P}\theta_{c}^{n+1},\,\phi_{h}\Big{)}}$
$\displaystyle=$	$\displaystyle\Big{(}\Phi D_{t}({\cal C}^{n+1}-\Pi_{c}^{n+1}{\cal C}^{n+1}),\,\phi_{h}\Big{)}-\Big{(}{\bf U}^{n+1}\cdot\nabla({\cal C}_{h}^{n}-{\cal C}^{n}),\,\phi_{h}\Big{)}$
	$\displaystyle-\Big{(}({\bf U}_{h}^{n+1}-{\bf U}^{n+1})\cdot\nabla{\cal C}_{h}^{n},\phi_{h}\Big{)}-\Big{(}({\Pi_{c}^{n+1}{\cal C}^{n+1}}-{\cal C}^{n+1})q^{P},\,\phi_{h}\Big{)}$
	$\displaystyle+\Big{(}(D({\widetilde{\bf U}}_{h}^{n+1})-D({\bf U}_{h}^{n+1}))\nabla\Pi_{c}^{n+1}{\cal C}^{n+1},\,\nabla\phi_{h}\Big{)}$
$\displaystyle:=$	$\displaystyle J_{1}(\phi_{h})+J_{2}(\phi_{h})+J_{3}(\phi_{h})+J_{4}(\phi_{h})+J_{5}(\phi_{h}).$	(3.27)

Taking ${\bf v}_{h}=\theta_{\bf u}^{n+1}$ in (3.25) leads to

\displaystyle\|\theta_{\bf u}^{n+1}\|_{L^{2}}\leq C\|\theta^{n}_{c}\|_{L^{2}}+Ch^{2}.

(3.28)

Taking $\phi_{h}=\theta_{c}^{n+1}$ in (3.27) and by Lemma 3.2, we get

	$\displaystyle\|J_{1}(\theta_{c}^{n+1})\|$	$\displaystyle\leq C(\\|\theta_{c}^{n+1}\\|_{L^{2}}^{2}+\\|D_{t}({\cal C}^{n+1}-\Pi_{c}^{n+1}{\cal C}^{n+1})\\|_{L^{2}}^{2}),$
	$\displaystyle\|J_{4}(\theta_{c}^{n+1})\|$	$\displaystyle\leq C\\|q^{P}\\|_{L^{3}}{\\|{\cal C}^{n+1}-\Pi_{c}^{n+1}{\cal C}^{n+1}\\|_{L^{2}}}\\|\theta_{c}^{n+1}\\|_{L^{6}}$
		$\displaystyle\leq{\epsilon\\|\nabla\theta_{c}^{n+1}\\|_{H^{1}}^{2}+C\epsilon^{-1}h^{4}},$

and by (3.28), we have

	$\displaystyle\|J_{5}(\theta_{c}^{n+1})\|$	$\displaystyle\leq C\\|\nabla\Pi_{c}^{n+1}{\cal C}^{n+1}\\|_{L^{\infty}}\\|\theta_{u}^{n+1}\\|_{L^{2}}\\|\nabla\theta_{c}^{n+1}\\|_{L^{2}}$
		$\displaystyle\leq\epsilon\\|\nabla\theta_{c}^{n+1}\\|_{L^{2}}^{2}+\epsilon^{-1}C\\|\theta_{c}^{n}\\|_{L^{2}}^{2}.$

Moreover, using integration by part and noting the fact that $\nabla\cdot{\bf U}^{n+1}=q^{I}-q^{P}$ and ${\bf U}^{n+1}\cdot{\bf n}=0$ on the boundary,

	$\displaystyle\|J_{2}(\theta_{c}^{n+1})\|$	$\displaystyle=\|({\bf U}^{n+1}\cdot\nabla(\theta_{c}^{n}+\Pi_{c}^{n}{\cal C}^{n}-{\cal C}^{n}),\,\theta_{c}^{n+1})\|$
		$\displaystyle=\|((q^{I}-q^{P})(\theta_{c}^{n}+\Pi_{c}^{n}{\cal C}^{n}-{\cal C}^{n}),\,\theta_{c}^{n+1})$
		$\displaystyle~{}~{}~{}+(\theta_{c}^{n}+\Pi_{c}^{n}{\cal C}^{n}-{\cal C}^{n},\,{\bf U}^{n+1}\cdot\nabla\theta_{c}^{n+1})\|$
		$\displaystyle\leq\epsilon\\|\theta_{c}^{n+1}\\|_{H^{1}}^{2}+C\epsilon^{-1}\\|\theta_{c}^{n}\\|_{L^{2}}^{2}+C\epsilon^{-1}h^{4}\,.$

Finally, we estimate

\displaystyle|J_{3}(\theta_{c}^{n+1})|=

\displaystyle|\big{(}({\bf U}_{h}^{n+1}-{\bf U}^{n+1})\cdot\nabla({\cal C}_{h}^{n}-{\cal C}^{n}),\,\theta_{c}^{n+1}\big{)}+\big{(}({\bf U}_{h}^{n+1}-{\bf U}^{n+1})\cdot\nabla{\cal C}^{n},\,\theta_{c}^{n+1}\big{)}\big{|}.

By (3.14)-(3.15) and (3.28),

		$\displaystyle\|\big{(}({\bf U}^{n+1}-{\bf U}_{h}^{n+1})\cdot\nabla({\cal C}_{h}^{n}-{\cal C}^{n}),\,\theta_{c}^{n+1}\big{)}\|$
	$\displaystyle\leq$	$\displaystyle C(\\|\theta_{u}^{n+1}\\|_{L^{2}}+\\|{\widetilde{\bf U}}_{h}^{n+1}-{\bf U}_{h}^{n+1}\\|_{L^{2}})(\\|\nabla\theta_{c}^{n}\\|_{L^{3}}+\\|\nabla(\Pi_{c}^{n}{\cal C}^{n}-{\cal C}^{n})\\|_{L^{3}})\\|\theta_{c}^{n+1}\\|_{L^{6}}$
	$\displaystyle\leq$	$\displaystyle C(\\|\theta^{n}_{c}\\|_{L^{2}}+h)(h^{-d/6}\\|\nabla\theta_{c}^{n}\\|_{L^{2}}+h)\\|\theta_{c}^{n+1}\\|_{L^{6}}$
	$\displaystyle\leq$	$\displaystyle C(h+h^{-d/6}\\|\nabla\theta_{c}^{n}\\|_{L^{2}})(\\|\theta_{c}^{n}\\|_{L^{2}}\\|\theta_{c}^{n+1}\\|_{H^{1}}+h\\|\theta_{c}^{n+1}\\|_{H^{1}})\,$
	$\displaystyle\leq$	$\displaystyle(Ch^{-d/6}\\|\theta^{n}_{c}\\|_{L^{2}}+Ch+\epsilon)(\\|\theta_{c}^{n}\\|_{H^{1}}^{2}+\\|\theta_{c}^{n+1}\\|_{H^{1}}^{2})+C\epsilon^{-1}h^{4},$

and

		$\displaystyle\big{\|}\big{(}({\bf U}_{h}^{n+1}-{\bf U}^{n+1})\cdot\nabla{\cal C}^{n},\,\theta_{c}^{n+1}\big{)}\big{\|}$
	$\displaystyle\leq$	$\displaystyle\|\big{(}\theta_{u}^{n+1}\cdot\nabla{\cal C}^{n},\,\theta_{c}^{n+1}\big{)}\big{\|}+\|\big{(}({\widetilde{\bf U}}_{h}^{n+1}-{\bf U}^{n+1})\cdot\nabla{\cal C}^{n},\,\theta_{c}^{n+1}\big{)}\big{\|}$
	$\displaystyle\leq$	$\displaystyle C\\|\theta_{u}^{n+1}\\|_{L^{2}}\\|\theta_{c}^{n+1}\\|_{L^{2}}+C\\|{\widetilde{\bf U}}_{h}^{n+1}-{\bf U}^{n+1}\\|_{H^{-1}}\\|\nabla{\cal C}^{n}\theta_{c}^{n+1}\\|_{H^{1}}$
	$\displaystyle\leq$	$\displaystyle\epsilon\\|\theta_{c}^{n+1}\\|_{H^{1}}^{2}+C\\|\theta^{n}_{c}\\|_{L^{2}}^{2}+C\\|\theta_{c}^{n+1}\\|_{L^{2}}^{2}+C\epsilon^{-1}h^{4}.$

Then we further have the estimate

	$\displaystyle\|J_{3}(\theta_{c}^{n+1})\|$	$\displaystyle\leq(Ch^{-d/6}\\|\theta^{n}_{c}\\|_{L^{2}}+Ch+\epsilon)(\\|\theta_{c}^{n}\\|_{H^{1}}^{2}+\\|\theta_{c}^{n+1}\\|_{H^{1}}^{2})$
		$\displaystyle+C\\|\theta^{n}_{c}\\|_{L^{2}}^{2}+C\\|\theta_{c}^{n+1}\\|_{L^{2}}^{2}+C\epsilon^{-1}h^{4}.$

It follows that

	$\displaystyle\frac{1}{2}D_{t}\\|\theta_{c}^{n+1}\\|^{2}_{L^{2}}+\big{\\|}\sqrt{D({\bf U}_{h}^{n+1})}\nabla\theta_{c}^{n+1}\big{\\|}^{2}_{L^{2}}$
$\displaystyle\leq$	$\displaystyle C(\epsilon+h^{-d/6}\\|\theta^{n}_{c}\\|_{L^{2}}+h)(\\|\nabla\theta_{c}^{n}\\|_{L^{2}}^{2}+\\|\nabla\theta_{c}^{n+1}\\|_{L^{2}}^{2})+C\epsilon^{-1}(\\|\theta_{c}^{n+1}\\|_{L^{2}}^{2}+\\|\theta^{n}_{c}\\|_{L^{2}}^{2})+C\epsilon^{-1}h^{4}$
	$\displaystyle~{}~{}~{}+C\big{\\|}D_{t}({\cal C}^{n+1}-\Pi_{c}^{n+1}{\cal C}^{n+1})\big{\\|}_{L^{2}}^{2}\,.$	(3.29)

Now we prove the $\tau$ -independent estimate

\displaystyle\|\theta_{c}^{n}\|_{L^{2}}\leq h

(3.30)

from (3.29) by mathematical induction. It is easy to see that $\|\theta_{c}^{0}\|_{L^{2}}\leq h$ , when $h<h_{2}$ for some $h_{2}>0$ . We assume that the inequality (3.30) holds for $0\leq n\leq k$ . Then there exists a positive constant $h_{3}$ such that when $h<h_{3}$ , (3.29) reduces to

		$\displaystyle D_{t}\\|\theta_{c}^{n+1}\\|^{2}_{L^{2}}+\big{\\|}\sqrt{D({\bf U}^{n+1}_{c})}\nabla\theta_{c}^{n+1})\big{\\|}^{2}_{L^{2}}$
	$\displaystyle\leq$	$\displaystyle Ch^{4}+C(\\|\theta_{c}^{n+1}\\|_{L^{2}}^{2}+\\|\theta_{c}^{n}\\|_{L^{2}}^{2})+C\big{\\|}D_{t}({\cal C}^{n+1}-\Pi_{c}^{n+1}{\cal C}^{n+1})\big{\\|}_{L^{2}}^{2}$

for $0\leq n\leq k$ . By Gronwall’s inequality and (3.19), there exists $\tau_{1}>0$ such that,

\displaystyle\|\theta_{c}^{n+1}\|_{L^{2}}\leq Ch^{2}

(3.31)

for $0\leq n\leq k$ and $\tau<\tau_{1}$ , which further implies that

\displaystyle\|\theta_{c}^{k+1}\|_{L^{2}}<h.

Taking $\widehat{\tau}_{0}\leq\min\{\tau_{0}^{*},\tau_{1}\}$ and $\widehat{h}_{0}\leq\min\{h_{1},h_{3}\}$ , the mathematical induction is closed and (3.30) holds for $1\leq n\leq N$ . Moreover, inequalities (3.28) and (3.31) hold for all $0\leq n\leq N-1$ .

It remains to estimate $\theta_{p}$ . In a traditional way, we consider the equation

\displaystyle-\nabla\cdot\biggl{(}\frac{k(x)}{\mu({\cal C}^{n})}\nabla g\biggl{)}=\theta_{p}^{n+1}

with the boundary condition $\frac{k(x)}{\mu({\cal C}^{n})}\nabla g\cdot{\bf n}=0$ on $\partial\Omega$ . It is easy to see that

\|g\|_{H^{2}}\leq C\|\theta_{p}^{n+1}\|_{L^{2}}.

Let

{\bf v}_{h}=Q_{h}\biggl{(}\frac{k(x)}{\mu({\cal C}^{n})}\nabla g\biggl{)}

where $Q_{h}:{\bf H(div)}\rightarrow{\bf H}_{h}^{r-1}$ is a projection such that [44] for ${\bf w}\in{\bf H(div)}$ ,

\displaystyle\Big{(}\nabla\cdot({\bf w}-Q_{h}{\bf w})\,,\chi_{h}\Big{)}=0,\quad\mbox{for~{}all}~{}~{}\chi_{h}\in S_{h}^{r-1}.

(3.32)

Then

(\varphi_{h},\nabla\cdot{\bf v}_{h})=-(\varphi_{h},\theta_{p}^{n+1}),~{}~{}\mbox{ for }\varphi_{h}\in S_{h}^{r-1}

and from (3.25) and the classical result $\|Q_{h}w\|_{L^{2}}\leq C\|w\|_{H^{1}}$ , we obtain

	$\displaystyle\\|\theta_{p}^{n+1}\\|_{L^{2}}^{2}$	$\displaystyle=\biggl{(}\frac{\mu({\cal C}^{n}_{h})}{k(x)}{\bf U}_{h}^{n+1}-\frac{\mu(\Pi_{c}^{n}{\cal C}^{n})}{k(x)}{\widetilde{\bf U}}_{h}^{n+1},Q_{h}\biggl{(}\frac{k(x)}{\mu({\cal C}^{n})}\nabla g\biggl{)}\biggl{)}$
		$\displaystyle\leq C(\\|\theta^{n}_{c}\\|_{L^{2}}+\\|\theta_{u}^{n+1}\\|_{L^{2}})\biggl{\\|}\frac{k(x)}{\mu({\cal C}^{n})}\nabla g\biggl{\\|}_{H^{1}}$
		$\displaystyle\leq Ch^{2}\\|\theta_{p}^{n+1}\\|_{L^{2}},$

which implies that

\displaystyle\|\theta_{p}^{n+1}\|_{L^{2}}\leq Ch^{2}.

(3.33)

(3.21) follows (3.28), (3.31) and (3.33). The proof of Lemma 3.3 is completed.

3.2 Proof of Theorem 2.1

For $r=1$ , Theorem 2.1 can be proved by combining Lemma 3.1 and Lemma 3.3. In this section, we only prove the case $r\geq 2$ .

From Lemma 3.1, Lemma 3.2, (3.14), (3.21) and the Gagliardo-Nirenburg inequality (2.10), we can see the boundedness of numerical solution

\displaystyle\|P_{h}^{n}\|_{L^{\infty}}+\|{\bf U}_{h}^{n}\|_{L^{\infty}}+\|{\cal C}_{h}^{n}\|_{W^{1,6}}\leq C.

(3.34)

Similarly, for any fixed integer $n\geq 1$ we denote by $({\widetilde{p}}_{h}^{n},{\widetilde{\bf u}}_{h}^{n})$ the classical mixed projection of $(p^{n},{\bf u}^{n})$ on $(\widehat{S}_{h}^{r-1},{\bf H}_{h}^{r-1})$ such that

	$\displaystyle\biggl{(}\frac{\mu(c^{n+1})}{k(x)}{\widetilde{\bf u}}_{h}^{n+1},\,{\bf v}_{h}\biggl{)}=-\Big{(}{\widetilde{p}}_{h}^{n+1},\,\nabla\cdot{\bf v}_{h}\Big{)},$		(3.35)
	$\displaystyle\Big{(}\nabla\cdot({\widetilde{\bf u}}_{h}^{n+1}-{\bf u}^{n+1}),\,\varphi_{h}\Big{)}=0,\quad\mbox{ for all }(\varphi_{h},v_{h})\in S_{h}^{r-1}\times{\bf H}_{h}^{r-1}.$		(3.36)

By the classical mixed method theory [9, 22, 44] and negative norm estimates in [20], we have

	$\displaystyle\\|{\bf u}^{n}-{\widetilde{\bf u}}_{h}^{n}\\|_{L^{2}}+\\|p^{n}-{\widetilde{p}}_{h}^{n}\\|_{L^{2}}\leq Ch^{r},$		(3.37)
	$\displaystyle\\|{\bf u}^{n}-{\widetilde{\bf u}}_{h}^{n}\\|_{H^{-1}}+\\|p^{n}-{\widetilde{p}}_{h}^{n}\\|_{H^{-1}}\leq Ch^{r+1},$		(3.38)

and by inverse inequalities and noting $r\geq 2$ ,

\displaystyle\|{\widetilde{\bf u}}_{h}^{n}\|_{L^{\infty}}\leq C.

(3.39)

For a given ${\bf u}^{n}$ , we define an elliptic quasi-projection ${\widetilde{\Pi}}_{c}^{n}:H^{1}(\Omega)\rightarrow V_{h}^{r}$ , slightly different from one in section 3.1, by

	$\displaystyle\Big{(}D({\bf u}^{n})\nabla({\widetilde{\Pi}}_{c}^{n}c^{n}-c^{n}),\,\nabla\phi_{h}\Big{)}+\Big{(}(D({\widetilde{\bf u}}_{h}^{n})-D({\bf u}^{n}))\nabla{\widetilde{\Pi}}_{c}^{n}c^{n},\,\nabla\phi_{h}\Big{)}=0,$
	$\displaystyle\quad\mbox{for~{}all}~{}~{}\phi_{h}\in V_{h}^{r},~{}~{}n\geq 1,$		(3.40)

with $\int_{\Omega}({\widetilde{\Pi}}_{c}^{n}c^{n}-c^{n})dx=0$ and ${\widetilde{\Pi}}_{c}^{0}:=I_{h}$ . By a proof similar to Lemma 3.2, we can get basic estimates of the elliptic quasi-projection ${\widetilde{\Pi}}_{c}^{n}$ as follows.

Lemma 3.4

Under the assumptions of Theorem 2.1, there exists $\widehat{h}_{1}>0$ such that for any $h\leq\widehat{h}_{1}$ ,

	$\displaystyle\\|c^{n}-{\widetilde{\Pi}}_{c}^{n}c^{n}\\|_{L^{2}}+h\\|\nabla(c^{n}-{\widetilde{\Pi}}_{c}^{n}c^{n})\\|_{L^{2}}\leq Ch^{r+1},$		(3.41)
	$\displaystyle\\|{\widetilde{\Pi}}_{c}^{n}c^{n}\\|_{W^{1,\infty}}\leq C,$		(3.42)

and

\displaystyle\biggl{(}\sum_{n=0}^{N-1}\tau\|D_{t}(c^{n}-{\widetilde{\Pi}}_{c}^{n}c^{n})\|_{L^{2}}^{2}\biggl{)}^{1/2}\leq Ch^{r+1}.

(3.43)

Now we start to prove Theorem 2.1. Let

{\widetilde{\theta}}_{p}^{n}=P_{h}^{n}-{\widetilde{p}}_{h}^{n},\quad{\widetilde{\theta}}_{u}^{n}={\bf U}_{h}^{n}-{\widetilde{\bf u}}_{h}^{n}\quad\mbox{ and }\quad{\widetilde{\theta}}_{c}^{n}={\cal C}_{h}^{n}-{\widetilde{\Pi}}_{c}^{n}c^{n}.

We prove below the estimate

\displaystyle\|{\widetilde{\theta}}_{p}^{n}\|_{L^{2}}+\|{\widetilde{\theta}}_{\bf u}^{n}\|_{L^{2}}+\|{\widetilde{\theta}}_{c}^{n}\|_{L^{2}}\leq C(\tau+h^{r+1}).

(3.44)

From (1.1)-(1.5) and the finite element system (2.7)-(2.9), the error function $({\widetilde{\theta}}_{p}^{n}$ , ${\widetilde{\theta}}_{u}^{n}$ , ${\widetilde{\theta}}_{c}^{n})$ satisfies

	$\displaystyle\biggl{(}\frac{\mu({\cal C}_{h}^{n})}{k(x)}{\widetilde{\theta}}_{\bf u}^{n+1}+(\frac{\mu({\cal C}_{h}^{n})}{k(x)}-\frac{\mu(c^{n+1})}{k(x)}){\widetilde{\bf u}}_{h}^{n+1},{\bf v}_{h}\biggl{)}=-\Big{(}{\widetilde{\theta}}_{p}^{n+1},\,\nabla\cdot{\bf v}_{h}\Big{)},$	(3.45)
	$\displaystyle\Big{(}\nabla\cdot{\widetilde{\theta}}_{\bf u}^{n+1},\,\varphi_{h}\Big{)}=0,$	(3.46)
	$\displaystyle\Big{(}\Phi D_{t}{\widetilde{\theta}}_{c}^{n+1},\,\phi_{h}\Big{)}+\Big{(}D({\bf U}_{h}^{n+1})\nabla{\widetilde{\theta}}_{c}^{n+1},\,\nabla\phi_{h}\Big{)}{+\Big{(}q^{P}{\widetilde{\theta}}_{c}^{n+1},\,\phi_{h}\Big{)}}$
$\displaystyle=$	$\displaystyle(T^{n+1}_{c},\phi_{h})+\Big{(}\Phi D_{t}(c^{n+1}-{\widetilde{\Pi}}_{c}^{n+1}c^{n+1}),\,\phi_{h}\Big{)}-\Big{(}{\bf u}^{n+1}\cdot\nabla({\cal C}_{h}^{n}-c^{n}),\,\phi_{h}\Big{)}$
	$\displaystyle-\Big{(}({\bf U}_{h}^{n+1}-{\bf u}^{n+1})\cdot\nabla{\cal C}_{h}^{n},\,\phi_{h}\Big{)}-\Big{(}{({\widetilde{\Pi}}_{c}^{n+1}c^{n+1}-c^{n+1})}q^{P},\,\phi_{h}\Big{)}$
	$\displaystyle+\Big{(}(D({\widetilde{\bf u}}_{h}^{n+1})-D({\bf U}_{h}^{n+1}))\nabla{\widetilde{\Pi}}_{c}^{n+1}c^{n+1},\,\nabla\phi_{h}\Big{)}$
$\displaystyle:=$	$\displaystyle(T^{n+1}_{c},\phi_{h})+{\widetilde{J}}_{1}(\phi_{h})+{\widetilde{J}}_{2}(\phi_{h})+{\widetilde{J}}_{3}(\phi_{h})+{\widetilde{J}}_{4}(\phi_{h})+{\widetilde{J}}_{5}(\phi_{h}),$	(3.47)

where $T^{n+1}_{c}$ denotes the truncation error. By the regularity assumption ((A1)), we have

\displaystyle\sum_{k=1}^{n}\tau\|T^{k}_{c}\|_{L^{2}}^{2}\leq C\tau^{2}.

(3.48)

Letting $v_{h}={\widetilde{\theta}}_{\bf u}^{n+1}$ in (3.45) and by (3.39), we further have

\displaystyle\|{\widetilde{\theta}}_{\bf u}^{n+1}\|_{L^{2}}\leq C\|{\widetilde{\theta}}^{n}_{c}\|_{L^{2}}+Ch^{r+1}+C\tau.

(3.49)

Taking $\phi_{h}={\widetilde{\theta}}_{c}^{n+1}$ in (3.47) and using Lemma 3.4 gives

	$\displaystyle\|{\widetilde{J}}_{1}({\widetilde{\theta}}_{c}^{n+1})\|$	$\displaystyle\leq C(\\|{\widetilde{\theta}}_{c}^{n+1}\\|_{L^{2}}^{2}+\\|D_{t}(c^{n+1}-{\widetilde{\Pi}}_{c}^{n+1}c^{n+1})\\|_{L^{2}}^{2}),$
	$\displaystyle\|{\widetilde{J}}_{4}({\widetilde{\theta}}_{c}^{n+1})\|$	$\displaystyle\leq C\\|q^{P}\\|_{L^{3}}{\\|c^{n+1}-{\widetilde{\Pi}}_{c}^{n+1}c^{n+1}\\|_{L^{2}}}\\|{\widetilde{\theta}}_{c}^{n+1}\\|_{L^{6}}$
		$\displaystyle\leq{\epsilon\\|{\widetilde{\theta}}_{c}^{n+1}\\|_{H^{1}}^{2}+C\epsilon^{-1}h^{2r+2}},$

and by (3.49), we get

	$\displaystyle\|{\widetilde{J}}_{5}({\widetilde{\theta}}_{c}^{n+1})\|$	$\displaystyle\leq C\\|\nabla{\widetilde{\Pi}}_{c}^{n+1}c^{n+1}\\|_{L^{\infty}}\\|{\widetilde{\theta}}_{\bf u}^{n+1}\\|_{L^{2}}\\|\nabla{\widetilde{\theta}}_{c}^{n+1}\\|_{L^{2}}$
		$\displaystyle\leq C\epsilon\\|\nabla{\widetilde{\theta}}_{c}^{n+1}\\|_{L^{2}}^{2}+\epsilon^{-1}\\|{\widetilde{\theta}}_{c}^{n}\\|_{L^{2}}^{2}+C\epsilon^{-1}(h^{2r+2}+\tau^{2}),$

Moreover, using integration by part and noting the fact that $\nabla\cdot{\bf u}^{n+1}=q^{I}-q^{P}$ and ${\bf u}^{n+1}\cdot{\bf n}=0$ on the boundary,

	$\displaystyle\|{\widetilde{J}}_{2}({\widetilde{\theta}}_{c}^{n+1})\|$	$\displaystyle=\|({\bf u}^{n+1}\cdot\nabla({\widetilde{\theta}}_{c}^{n}+{\widetilde{\Pi}}_{c}^{n}c^{n}-c^{n}),\,{\widetilde{\theta}}_{c}^{n+1})\|$
		$\displaystyle=\|((q^{I}-q^{P})({\widetilde{\theta}}_{c}^{n}+{\widetilde{\Pi}}_{c}^{n}c^{n}-c^{n}),\,{\widetilde{\theta}}_{c}^{n+1})$
		$\displaystyle~{}~{}~{}+({\widetilde{\theta}}_{c}^{n}+{\widetilde{\Pi}}_{c}^{n}c^{n}-c^{n},\,{\bf u}^{n+1}\cdot\nabla{\widetilde{\theta}}_{c}^{n+1})\|$
		$\displaystyle\leq\epsilon\\|{\widetilde{\theta}}_{c}^{n+1}\\|_{H^{1}}^{2}+C\epsilon^{-1}\\|{\widetilde{\theta}}_{c}^{n}\\|_{L^{2}}^{2}+C\epsilon^{-1}h^{2r+2}\,.$

Finally, we rewrite $\widetilde{J}_{3}$ by

\displaystyle|{\widetilde{J}}_{3}({\widetilde{\theta}}_{c}^{n+1})|=

\displaystyle|\big{(}({\bf U}_{h}^{n+1}-{\bf u}^{n+1})\cdot\nabla({\cal C}_{h}^{n}-c^{n}),\,{\widetilde{\theta}}_{c}^{n+1}\big{)}+\big{(}({\bf U}_{h}^{n+1}-{\bf u}^{n+1})\cdot\nabla c^{n},\,{\widetilde{\theta}}_{c}^{n+1}\big{)}\big{|}.

By (3.37)-(3.38) and (3.49), we have

		$\displaystyle\|\big{(}({\bf u}^{n+1}-{\bf U}_{h}^{n+1})\cdot\nabla({\cal C}_{h}^{n}-c^{n}),\,{\widetilde{\theta}}_{c}^{n+1}\big{)}\|$
	$\displaystyle\leq$	$\displaystyle C(\\|\theta_{\bf u}^{n+1}\\|_{L^{3}}+\\|{\widetilde{\bf U}}_{h}^{n+1}-{\bf U}^{n+1}\\|_{L^{3}}+\\|{\bf U}^{n+1}-{\bf u}^{n+1}\\|_{L^{3}})$
		$\displaystyle(\\|\nabla{\widetilde{\theta}}_{c}^{n}\\|_{L^{2}}+\\|\nabla({\widetilde{\Pi}}_{c}^{n}c^{n}-c^{n})\\|_{L^{2}})\\|{\widetilde{\theta}}_{c}^{n+1}\\|_{L^{6}}$
	$\displaystyle\leq$	$\displaystyle C(\tau+h)(\\|\nabla{\widetilde{\theta}}_{c}^{n}\\|_{L^{2}}+h^{r})\\|{\widetilde{\theta}}_{c}^{n+1}\\|_{L^{6}}$
	$\displaystyle\leq$	$\displaystyle(C\tau+Ch+\epsilon)(\\|{\widetilde{\theta}}_{c}^{n}\\|_{H^{1}}^{2}+\\|{\widetilde{\theta}}_{c}^{n+1}\\|_{H^{1}}^{2})+C\epsilon^{-1}(h^{2r+2}+\tau^{2}),$

and

		$\displaystyle\big{\|}\big{(}({\bf U}_{h}^{n+1}-{\bf u}^{n+1})\cdot\nabla c^{n},\,{\widetilde{\theta}}_{c}^{n+1}\big{)}\big{\|}$
	$\displaystyle\leq$	$\displaystyle\|\big{(}{\widetilde{\theta}}_{\bf u}^{n+1}\cdot\nabla c^{n},\,{\widetilde{\theta}}_{c}^{n+1}\big{)}\big{\|}+\|\big{(}({\widetilde{\bf u}}_{h}^{n+1}-{\bf u}^{n+1})\cdot\nabla c^{n},\,{\widetilde{\theta}}_{c}^{n+1}\big{)}\big{\|}$
	$\displaystyle\leq$	$\displaystyle C\\|{\widetilde{\theta}}_{u}^{n+1}\\|_{L^{2}}\\|{\widetilde{\theta}}_{c}^{n+1}\\|_{L^{2}}+C\\|{\widetilde{\bf u}}_{h}^{n+1}-{\bf u}^{n+1}\\|_{H^{-1}}\\|\nabla c^{n}{\widetilde{\theta}}_{c}^{n+1}\\|_{H^{1}}$
	$\displaystyle\leq$	$\displaystyle\epsilon\\|{\widetilde{\theta}}_{c}^{n+1}\\|_{H^{1}}^{2}+C\\|{\widetilde{\theta}}^{n}_{c}\\|_{L^{2}}^{2}+C\\|{\widetilde{\theta}}_{c}^{n+1}\\|_{L^{2}}^{2}+C\epsilon^{-1}(h^{2r+2}+\tau^{2}).$

It follows that

\displaystyle|J_{3}({\widetilde{\theta}}_{c}^{n+1})|\leq(C\tau+Ch+\epsilon)(\|{\widetilde{\theta}}_{c}^{n}\|_{H^{1}}^{2}+\|{\widetilde{\theta}}_{c}^{n+1}\|_{H^{1}}^{2})+C\|{\widetilde{\theta}}^{n}_{c}\|_{L^{2}}^{2}+C\|{\widetilde{\theta}}_{c}^{n+1}\|_{L^{2}}^{2}+C\epsilon^{-1}(h^{2r+2}+\tau^{2}).

Therefore,

	$\displaystyle\frac{1}{2}D_{t}\\|{\widetilde{\theta}}_{c}^{n+1}\\|^{2}_{L^{2}}+\big{\\|}\sqrt{D({\bf U}_{h}^{n+1})}\nabla{\widetilde{\theta}}_{c}^{n+1}\big{\\|}^{2}_{L^{2}}$
$\displaystyle\leq$	$\displaystyle C(\epsilon+\tau+h)(\\|\nabla{\widetilde{\theta}}_{c}^{n}\\|_{L^{2}}^{2}+\\|\nabla{\widetilde{\theta}}_{c}^{n+1}\\|_{L^{2}}^{2})+C\epsilon^{-1}(\\|{\widetilde{\theta}}_{c}^{n+1}\\|_{L^{2}}^{2}+\\|{\widetilde{\theta}}^{n}_{c}\\|_{L^{2}}^{2})+C\epsilon^{-1}(h^{2r+2}+\tau^{2})$
	$\displaystyle~{}~{}~{}+C\big{\\|}D_{t}(c^{n+1}-{\widetilde{\Pi}}_{c}^{n+1}c^{n+1})\big{\\|}_{L^{2}}^{2}+C\\|T_{c}^{n+1}\\|_{L^{2}}^{2}\,.$	(3.50)

By Gronwall’s inequality and (3.43), there exists $\tau_{0}>0$ and $h_{0}>0$ such that,

\displaystyle\|{\widetilde{\theta}}_{c}^{n+1}\|_{L^{2}}\leq C(\tau+h^{r+1})

(3.51)

for $0\leq n\leq k$ , when $\tau<\tau_{0}\leq\widehat{\tau}_{0}$ and $h<h_{0}\leq\widehat{h}_{0}$ .

Similarly to the estimate for $\theta_{p}^{n}$ in section 3.1, we can get the estimate to ${\widetilde{\theta}}_{p}^{n}$

\displaystyle\|{\widetilde{\theta}}_{p}^{n+1}\|_{L^{2}}\leq C(\tau+h^{r+1}).

(3.52)

(3.44) follows (3.49), (3.51) and (3.52). Theorem 2.1 is proved by combining (3.37)-(3.38), (3.44) and the basic projection error estimates in Lemma 3.4. The proof is complete.

4 Proof to Lemma 3.2

To prove Lemma 3.2, we define an extra classical elliptic projection $R_{h}^{n}$ : $H^{1}(\Omega)\rightarrow V_{h}^{r}$ by

\displaystyle\Big{(}D({\bf U}^{n})\nabla(R_{h}^{n}\psi-\psi),\nabla\phi_{h}\Big{)}=0,\qquad\qquad\mbox{for~{}all}~{}~{}\phi_{h}\in V_{h}^{r},

(4.1)

with $\int_{\Omega}(\psi-R_{h}^{n}\psi)dx=0$ . Similarly, by classical FE theory [9, 44],

\displaystyle\|R_{h}^{n}\psi\|_{W^{1,p}}\leq C\|\psi\|_{W^{1,p}},\quad 1<p\leq 4.

(4.2)

For $2\leq p\leq 4$ and $1/q+1/p=1$ , by (4.1) we have

		$\displaystyle\|(D({\bf U}^{n})\nabla(\Pi_{c}^{n}{\cal C}^{n}-{\cal C}^{n}),\nabla v)\|$
	$\displaystyle\leq$	$\displaystyle\|(D({\bf U}^{n})\nabla(\Pi_{c}^{n}{\cal C}^{n}-I_{h}{\cal C}^{n}),\nabla v)\|+\|(D({\bf U}^{n})\nabla(I_{h}{\cal C}^{n}-{\cal C}^{n}),\nabla v)\|$
	$\displaystyle\leq$	$\displaystyle\|(D({\bf U}^{n})\nabla(\Pi_{c}^{n}{\cal C}^{n}-I_{h}{\cal C}^{n}),\nabla R_{h}^{n}v)\|+C\\|\nabla(I_{h}{\cal C}^{n}-{\cal C}^{n})\\|_{L^{p}}\\|\nabla v\\|_{L^{q}}$
	$\displaystyle\leq$	$\displaystyle\|(D({\bf U}^{n})\nabla(\Pi_{c}^{n}{\cal C}^{n}-{\cal C}^{n}),\nabla R_{h}^{n}v)\|+Ch(\\|\nabla v\\|_{L^{q}}+\\|\nabla R_{h}^{n}v\\|_{L^{q}}),$

where we have noted the fact that $D({\bf U}^{n})$ and $D({\widetilde{\bf U}}_{h}^{n})$ are symmetric. With (3.1) and (4.2) we further have

	$\displaystyle\|(D({\bf U}^{n})\nabla(\Pi_{c}^{n}{\cal C}^{n}-{\cal C}^{n}),\nabla v)\|$
$\displaystyle\leq$	$\displaystyle\|((D({\bf U}^{n})-D({\widetilde{\bf U}}_{h}^{n}))\nabla\Pi_{c}^{n+1}{\cal C}^{n+1},\nabla R_{h}^{n}v)\|+Ch\\|\nabla v\\|_{L^{q}}$
$\displaystyle\leq$	$\displaystyle C\left(\\|{\bf U}^{n}-{\widetilde{\bf U}}_{h}^{n}\\|_{L^{p}}\\|\nabla{\cal C}^{n}\\|_{L^{\infty}}+\\|{\bf U}^{n}-{\widetilde{\bf U}}_{h}^{n}\\|_{L^{\infty}}\\|\nabla(\Pi_{c}^{n}{\cal C}^{n}-{\cal C}^{n})\\|_{L^{p}}+h\right)\\|\nabla v\\|_{L^{q}}$
$\displaystyle\leq$	$\displaystyle C\left(h+h^{1/4}\\|\nabla(\Pi_{c}^{n+1}{\cal C}^{n+1}-{\cal C}^{n+1})\\|_{L^{p}}\right)\\|\nabla v\\|_{L^{q}},$	(4.3)

where we have also used some basic estimates to $R_{h}^{n}$ . It follows that

\displaystyle\|\nabla(\Pi_{c}^{n}{\cal C}^{n}-{\cal C}^{n})\|_{L^{p}}\leq C\|D({\bf U}^{n})\nabla(\Pi_{c}^{n}{\cal C}^{n}-{\cal C}^{n})\|_{L^{p}}\leq Ch,

(4.4)

when $h\leq h_{1}$ for some $h_{1}>0$ . By using an inverse inequality, we see that

	$\displaystyle\\|\Pi_{c}^{n}{\cal C}^{n}\\|_{W^{1,\infty}}\leq$	$\displaystyle\\|I_{h}{\cal C}^{n}\\|_{W^{1,\infty}}+\\|\Pi_{c}^{n}{\cal C}^{n}-I_{h}{\cal C}^{n}\\|_{W^{1,\infty}}$
	$\displaystyle\leq$	$\displaystyle Ch^{-d/4}\\|\Pi_{c}^{n}{\cal C}^{n}-I_{h}{\cal C}^{n}\\|_{W^{1,4}}+C$
	$\displaystyle\leq$	$\displaystyle Ch^{1-d/4}+C.$

We have proved (3.18) and the second part of (3.17).

To get the $L^{2}$ -norm estimate in (3.17), we use the duality argument and consider the equation

\displaystyle\left.\begin{aligned} &-\nabla\cdot\Big{(}D({\bf U}^{n})\nabla w\Big{)}={\cal C}^{n}-\Pi_{c}^{n}{\cal C}^{n},\qquad&&\mbox{in }\Omega,\\ &\hskip 72.26999ptD({\bf U}^{n})\nabla w\cdot{\bf n}=0,\qquad&&\mbox{on }\partial\Omega,\end{aligned}\right.

(4.5)

with $\int_{\Omega}wdx=0$ . Its solution satisfies

\displaystyle\|w\|_{H^{2}}\leq C\|{\cal C}^{n}-\Pi_{c}^{n}{\cal C}^{n}\|_{L^{2}}\,.

(4.6)

By (3.1),

	$\displaystyle\\|{\cal C}^{n}-\Pi_{c}^{n}{\cal C}^{n}\\|_{L^{2}}^{2}$
$\displaystyle=$	$\displaystyle\Big{(}D({\bf U}^{n})\nabla(w-I_{h}w),\nabla({\cal C}^{n}-\Pi_{c}^{n}{\cal C}^{n})\Big{)}$
	$\displaystyle+\Big{(}\nabla(I_{h}w),D({\bf U}^{n})\nabla({\cal C}^{n}-\Pi_{c}^{n}{\cal C}^{n})\Big{)}$
$\displaystyle=$	$\displaystyle\Big{(}D({\bf U}^{n})\nabla(w-I_{h}w),\nabla({\cal C}^{n}-\Pi_{c}^{n}{\cal C}^{n})\Big{)}$
	$\displaystyle+\Big{(}(D({\widetilde{\bf U}}_{h}^{n})-D({\bf U}^{n}))\nabla\Pi_{c}^{n}{\cal C}^{n},\nabla I_{h}w\Big{)}$
$\displaystyle=$	$\displaystyle\Big{(}D({\bf U}^{n})\nabla(w-I_{h}w),\nabla({\cal C}^{n}-\Pi_{c}^{n}{\cal C}^{n})\Big{)}$
	$\displaystyle+\Big{(}(D({\widetilde{\bf U}}_{h}^{n})-D({\bf U}^{n}))\nabla{\cal C}^{n},\nabla I_{h}w\Big{)}$
	$\displaystyle+\Big{(}(D({\widetilde{\bf U}}_{h}^{n})-D({\bf U}^{n}))(\nabla\Pi_{c}^{n}{\cal C}^{n}-\nabla{\cal C}^{n}),\nabla I_{h}w\Big{)}$	(4.7)

For the second term in the right hand side of the last equation, we have the following estimate

	$\displaystyle\|((D({\widetilde{\bf U}}_{h}^{n})-D({\bf U}^{n}))\nabla{\cal C}^{n},\nabla I_{h}w)\|$
$\displaystyle\leq$	$\displaystyle\|((D({\widetilde{\bf U}}_{h}^{n})-D({\bf U}^{n}))\nabla{\cal C}^{n},\nabla w)\|+\|((D({\widetilde{\bf U}}_{h}^{n})-D({\bf U}^{n}))\nabla{\cal C}^{n},\nabla(I_{h}w-w))\|$
$\displaystyle\leq$	$\displaystyle C\\|{\widetilde{\bf U}}_{h}^{n}-{\bf U}^{n}\\|_{H^{-1}}\\|\nabla{\cal C}^{n}\nabla w\\|_{H^{1}}+C\\|{\widetilde{\bf U}}_{h}^{n}-{\bf U}^{n}\\|_{L^{2}}\\|\nabla(I_{h}w-w)\\|_{L^{2}}\\|\nabla{\cal C}^{n}\\|_{L^{\infty}}$
$\displaystyle\leq$	$\displaystyle Ch^{2}\\|w\\|_{H^{2}}.$	(4.8)

where we have used (3.14)-(3.15).

	$\displaystyle\\|{\cal C}^{n}-\Pi_{c}^{n}{\cal C}^{n}\\|_{L^{2}}^{2}\leq$	$\displaystyle Ch\\|w\\|_{H^{2}}\\|\nabla({\cal C}^{n}-\Pi_{c}^{n}{\cal C}^{n})\\|_{L^{2}}+Ch^{2}\\|{\cal C}^{n}-\Pi_{c}^{n}{\cal C}^{n}\\|_{L^{2}}$
		$\displaystyle+\\|{\widetilde{\bf U}}_{h}^{n}-{\bf U}^{n}\\|_{L^{3}}\\|\nabla(\Pi_{c}^{n}{\cal C}^{n}-{\cal C}^{n})\\|_{L^{2}}\\|\nabla I_{h}w\\|_{L^{6}}$
	$\displaystyle\leq$	$\displaystyle Ch^{2}\\|{\cal C}^{n}-\Pi_{c}^{n}{\cal C}^{n}\\|_{L^{2}}\,,$

and (3.17) follows immediately.

It remains to prove (3.19). From (4.1), we see that

	$\displaystyle\Big{(}D({\bf U}^{n})\nabla D_{\tau}(\Pi_{c}^{n}{\cal C}^{n}-{\cal C}^{n}),\nabla\phi_{h}\Big{)}+\Big{(}D_{\tau}(D({\bf U}^{n}))\nabla(\Pi_{c}^{n-1}{\cal C}^{n-1}-{\cal C}^{n-1}),\nabla\phi_{h}\Big{)}$
	$\displaystyle+\Big{(}(D({\widetilde{\bf U}}_{h}^{n})-D({\bf U}^{n}))\nabla(D_{\tau}\Pi_{c}^{n}{\cal C}^{n}),\nabla\phi_{h}\Big{)}+\Big{(}D_{\tau}(D({\widetilde{\bf U}}_{h}^{n})-D({\bf U}^{n}))\Pi_{c}^{n-1}\nabla{\cal C}^{n-1},\nabla\phi_{h}\Big{)}=0.$

Using similar proof for (3.17), we can get the desired result. The proof is complete.

5 Numerical examples

In this section, we present numerical results for incompressible miscible flows in both two and three-dimensional porous media to confirm our theoretical analysis and show the efficiency of linearized Galerkin FEMs. Here we always assume that the solution of the system is smooth. The problem with non-smooth solutions was considered in [7], where the convergence of a discontinuous-mixed FEM was proved. All computations in this section are performed by using the software FEniCS [36].

We rewrite the system (1.1)-(1.6) by

	$\displaystyle\frac{\partial c}{\partial t}-\nabla\cdot(D({\bf u})\nabla c)+{\bf u}\cdot\nabla c=g,$		(5.1)
	$\displaystyle\nabla\cdot{\bf u}=f,$		(5.2)
	$\displaystyle{\bf u}=-\frac{1}{\mu(c)}\nabla p,$		(5.3)

where $D({\bf u})=1+|{\bf u}|^{2}/(1+|{\bf u}|^{2})+{\bf u}\otimes{\bf u}$ and $\mu(c)=1+c^{2}$ .

Example 5.1. First, we consider a two-dimensional model where $\Omega=[0,1]\times[0,1]$ . We set the terminal time $T=1.0$ . The functions $f$ and $g$ are chosen correspondingly to the exact solution

	$\displaystyle p=1+1000x^{2}(1-x)^{3}y^{2}(1-y)^{3}t^{2}e^{-t},$		(5.4)
	$\displaystyle c=0.2+50x^{2}(1-x)^{2}y^{2}(1-y)^{2}te^{t},$		(5.5)

which satisfies the boundary condition (1.5).

Refer to caption — Figure 1: A uniform triangular mesh on the unit square domain with M = 8

Table 1:

L^{2}

-norm errors of linearized Galerkin-mixed FEM (2.7)-(2.9) in 2D (Example 5.1).

$\tau=\frac{8}{M^{2}}(r=1)$	$Err_{P}$	$Err_{\bf U}$	$Err_{\cal C}$
M = 8	2.63e-02	1.99e-01	5.09e-02
M = 16	1.29e-02	1.01e-01	1.20e-02
M = 32	6.38e-03	5.07e-02	2.93e-03
M = 64	3.18e-03	2.54e-02	7.29e-04
M =128	1.59e-03	1.27e-02	1.82e-04
order	1.01	0.99	2.03
$\tau=\frac{16}{M^{3}}(r=2)$	$Err_{P}$	$Err_{\bf U}$	$Err_{\cal C}$
M = 8	3.48e-03	2.81e-02	4.66e-03
M = 16	8.53e-04	7.23e-03	5.64e-04
M = 32	2.12e-04	1.82e-03	6.94e-05
M = 64	5.30e-05	4.56e-04	8.62e-06
M =128	1.33e-05	1.14e-04	1.08e-06
order	2.01	1.99	3.02

We use a uniform triangular mesh with M+1 vertices in each direction, where $h=\frac{\sqrt{2}}{M}$ (see Figure 1 for the illustration with $M=8$ ). We solve the system (5.1)-(5.3) by the scheme (2.7)-(2.9) with $r=1$ and $r=2$ , respectively.

As the expected optimal convergence rate in $L^{2}$ -norm is $O(\tau+h^{r+1})$ , we set $\tau=\frac{8}{M^{r+1}}$ in our computation. The $L^{2}$ -norm errors of $P_{h}^{N}$ , ${\bf U}_{h}^{N}$ and ${\cal C}_{h}^{N}$ are presented in Table 1 for $r=1,2$ , where $Err_{v}=\|v_{h}^{N}-v(\cdot,t_{N})\|_{L^{2}}$ . From Table 1, we can see clearly that the convergence rates for $P_{h}^{N}$ and ${\bf U}_{h}^{N}$ are optimal with the order $O(h^{r})$ , while the rate for ${\cal C}_{h}^{n}$ is optimal with the order $O(h^{r+1})$ for both $r=1,2$ . The one-order lower approximation to $(p,{\bf u})$ does not affect the accuracy of the numerical concentration and the mesh-size restriction in (1.10) is not necessary.

Example 5.2. Secondly, we study the system (5.1)-(5.3) in a three-dimensional cube $[0,1]\times[0,1]\times[0,1]$ ., where the functions $f$ and $g$ are chosen correspondingly to the smooth exact solution

	$\displaystyle p=1.0+1000x^{2}(1-x)^{3}y^{2}(1-y)^{3}z^{2}(1-z)^{3}t^{2}e^{-t},$		(5.6)
	$\displaystyle c=0.2+50x^{2}(1-x)^{2}y^{2}(1-y)^{2}z^{2}(1-z)^{2}te^{t}.$		(5.7)

We also set the terminal time T = 1.0 in this example.

A uniform triangular mesh with M+1 vertices in each direction of the cube is used in our computation, where $h=\frac{\sqrt{3}}{M}$ . We solve the system by the linearized Galerkin FEMs in (2.7)-(2.9) by the lowest-order Galerkin-mixed FEM ( $r=1$ ). Such a Galerkin-mixed FEM is most popular in practical computation, particular for problems with discontinuous media. To show the accuracy in spatial direction, we set $\tau=\frac{8}{M^{2}}$ in our computation. We present in Table 2 the $L^{2}$ -norm errors of concentration, velocity and pressure. Again numerical results confirm our theoretical analysis and the accuracy of the numerical concentration is in the order $O(h^{2})$ , while all previous analyses only showed the first-order accuracy of numerical concentration.

Table 2:

L^{2}

-norm errors of linearized Galerkin FEM (2.7)-(2.9) in 3D (Example 5.2)

$\tau=\frac{8}{M^{2}}$	$Err_{P}$	$Err_{\bf U}$	$Err_{\cal C}$
M = 8	5.70e-04	5.36e-03	9.05e-04
M = 16	2.82e-04	2.72e-03	2.40e-04
M = 32	1.40e-04	1.36e-03	6.10e-05
M = 64	7.55e-05	7.13e-04	1.38e-05
order	0.98	0.97	2.01

6 Conclusions

We have presented unconditionally optimal error analysis of commonly-used Galerkin-mixed FEMs for a nonlinear and strongly coupled parabolic system from incompressible miscible flow in porous media. In particular, for the most popular lowest-order Galerkin-mixed method, we show unconditionally the second-order accuracy $O(h^{2})$ for the numerical concentration ${\cal C}_{h}^{n}$ in spatial direction. In all previous works, only the first-order accuracy was obtained under certain time-step restriction and the mesh-size condition. Moreover, our analysis is based on a quasi-projection and the approach presented in this paper can be extended to many other coupled nonlinear parabolic PDEs to obtain optimal error estimates for all components. With the numerical concentration $C_{h}^{n}\in V_{h}^{r}$ obtained by the lowest-order Galerkin-mixed FEM, one can get new numerical velocity and pressure of the accuracy $O(h^{2})$ (same as the concentration) by resolving the elliptic pressure equation (1.5) at a given time with a higher-order Galerkin FEM or a higher-order mixed method.

Optimal error estimates presented in this paper is based on strong regularity assumptions of the solution of the system and physical parameters as usual. Since the present paper focuses on the new analysis of classical Galerkin-mixed FEMs, the regularity assumptions in $(A1)-(A6)$ may not be optimal. Some related works under weaker assumptions were done by several authors [11, 35, 38]. Also the existence and uniqueness of the strong solution of the system (1.1)-(1.2) in a general setting remains open. A systematic numerical simulation on the above scheme with many other approximations will be presented in the coming article [49] and the problems with non-smooth domains and discontinuous coefficients will be also studied there.

Acknowledgements The authors would like to thank the anonymous referees for their valuable suggestions and comments.

References

[1] A. Agosti, L. Formaggia and A. Scotti, Analysis of a model for precipitation and dissolution coupled with a Darcy flux, J. Math. Anal. Appl., 431(2015), 752–781.
[2] B. Amaziane and M. El Ossmani, Convergence analysis of an approximation to miscible fluid flows in porous media by combining mixed finite element and finite volume methods, Numer. Methods Partial Diff. Eq., 24(2008), 799–832.
[3] R. An and J. Su, Optimal error estimates of semi-implicit Galerkin method for time-dependent Nematic liquid crystal flows, J. Scientific Computing, 74(2018), 979–1008.
[4] T. Arbogast and W. Wang, Stability, monotonicity, maximum and minimum principles, and implementation of the volume corrected characteristic method, SIAM J. Sci. Comput., 33(2011), 1549–1573.
[5] T. Arbogast, M.F. Wheeler and N. Zhang, A nonlinear mixed finite element method for a degenerate parabolic equation arising in flow in porous media, SIAM J. Numer. Anal., 33(1996), 1669–1687.
[6] J.W. Barrett and P. Knabner, Finite element approximation of the transport of reactive solutes in porous media. Part 1: error estimates for nonequilibrium adsorption processes, SIAM J. Numer. Anal., 34(1997), 201–227.
[7] S. Bartels, M. Jensen and R. Muller, Discontinuous Galerkin finite element convergence for incompressible miscible displacement problems of low regularity, SIAM J. Numer. Anal., 47(2009), 3720–3742.
[8] J. Bear and Y. Bachmat, Introduction to Modeling of Transport Phenomena in Porous Media, Springer-Verlag, New York, 1990.
[9] S. Brenner and L. Scott, The Mathematical Theory of Finite Element Methods, Springer, New York, 2002.
[10] W. Cai, Max Gunzburger and J. Wang, Convergence analysis of Crank-Nicolson Galerkin-Galerkin FEMs for miscible displacement in porous media, to appear.
[11] W. Cai, B. Li, Y. Lin and W. Sun, Analysis of fully discrete FEM for miscible displacement in porous media with Bear–Scheidegger diffusion tensor, Numer. Math., 141(2019), 1009–1042.
[12] F. Chen, H. Chen and H. Wang, An optimal-order error estimate for a Galerkin-mixed finite-element time-stepping procedure for porous media flows, Numer. Methods Partial Differ. Equations, 28.2(2012), 707–719.
[13] Z. Chen and R. Ewing, Mathematical analysis for reservoir models, SIAM J. Math. Anal., 30 (1999), 43–453.
[14] Z. Chen, G. Huan and Y. Ma, Computational Methods for Multiphase Flows in Porous Media, Computational science and engineering, SIAM, PA, 2006.
[15] A. Cheng, K. Wang and H. Wang, Superconvergence for a time-discretization procedure for the mixed finite element approximation of miscible displacement in porous media, Numer. Methods Partial Differ. Equations, 28(2012), 1382–1398.
[16] K. Chrysafinos and L.S. Hou, Error estimates for semidiscrete finite element approximations of linear and semilinear parabolic equations under minimal regularity assumptions, SIAM J. Numer. Anal., 40(2002), 282–306.
[17] J. Douglas, JR., R.E. Ewing and M.F. Wheeler, The approximation of the pressure by a mixed method in the simulation of miscible displacement, RAIRO Analyse numérique, 17(1983), 17–33.
[18] J. Douglas, JR., R. Ewing and M.F. Wheeler, A time-discretization procedure for a mixed finite element approximation of miscible displacement in porous media, RAIRO Anal. Numer., 17(1983), 249–265.
[19] J. Douglas, JR., F. Furtada and F. Pereira, On the numerical simulation of waterflooding of heterogeneous petroleum reservoirs, Comput. Geosciences, 1(1997), 155–190.
[20] J. Douglas, JR. and J.E. Roberts, Global estimates for mixed methods for second order elliptic equations, Math. Comput., 44(1985), 39–52.
[21] R.G. Durán, On the approximation of miscible displacement in porous media by a method of characteristics combined with a mixed method, SIAM J. Numer. Anal., 25(1988), 989–1001.
[22] R.G. Durán, Error analysis in $L^{p}$ , $1\leq p\leq\infty$ , for mixed finite element methods for linear and quasi-linear elliptic problems, RAIRO Mod. Math. Anal. Numér., 22(1988), 371–387.
[23] V.J. Ervin and W.W. Miles, Approximation of time-dependent viscoelastic fluid flow: SUPG approximation, SIAM J. Numer. Anal., 41(2003), 457–486.
[24] R.E. Ewing and T.F. Russell, Efficient time-stepping methods for miscible displacement problems in porous media, SIAM J. Numer. Anal., 19(1982), 1–67.
[25] R.E. Ewing and M.F. Wheeler, Galerkin methods for miscible displacement problems in porous media, SIAM J. Numer. Anal., 17(1980), 351–365.
[26] R.E. Ewing, ed, The mathematics of Reservoir Simulation, Frontiers in Applied Mathematics, SIAM, Philadelphia, PA, 1983.
[27] R.E. Ewing, T.F. Russell and M.F. Wheeler, Convergence analysis of an approximation of miscible displacement in porous media by mixed finite elements and a modified method of characteristics, Comput. Methods Appl. Mech. Engrg., 47(1984), 73–92.
[28] R.E. Ewing and H. Wang, A summary of numerical methods for time-dependent advection-dominated partial differential equations, J. Comput. Appl. Math., 128(2001), 423–445.
[29] X. Feng, On existence and uniqueness results for a coupled system modeling miscible displacement in porous media, J. Math. Anal. Appl., 194(1995), 883–910.
[30] X. Feng, Recent developments on modeling and analysis of flow of miscible fluids in porous media, Contemp. Math., 295(2002), 229–240.
[31] X. Feng and M. Neilan, A modified characteristic finite element method for a fully nonlinear formulation of the semigeostrophic flow equations, SIAM J. Numer. Anal., 47(2009), 2952–2981.
[32] S. Kumar, N. Nataraj and A.K. Pani, Finite volume element method for the incompressible miscible displacement problems in porous media, Proc. Appl. Math. Mech., 7(2007), 2020015–2020016.
[33] S. Kumar and S. Yadav, Modified method of characteristics combined with finite volume element methods for incompressible miscible displacement problems in porous media. Int. J. PDEs, 2014(2014).
[34] B. Li and W. Sun, Unconditional convergence and optimal error estimates of a Galerkin-mixed FEM for incompressible miscible flow in porous media, SIAM J. Numer. Anal., 51(2013), 1959–1977.
[35] B. Li and W. Sun, Regularity of the diffusion-dispersion tensor and error analysis of FEMs for a porous media flow, SIAM J. Numer. Anal., 53(2015), 1418–1437.
[36] A. Logg, K. Mardal and G. Wells (Eds.), Automated Solution of Differential Equations by the Finite Element Method, Springer, Berlin, 2012.
[37] S.M.C. Malta, A.F.D. Loula and E.L.M Garcia, Numerical analysis of a stabilized finite element method for tracer injection simulations, Comput. Methods Appl. Mech. Engrg., 187(2000), 119–136.
[38] F.A. Radu, K. Kumar, J.M. Nordbotten and I.S. Pop, A robust, mass conservative scheme for two-phase flow in porous media including Hölder continuous nonlinearities, IMA J. Numer. Anal., 38(2018), 884–920.
[39] F.A. Radu, I.S. Pop and S. Attinger, Analysis of an Euler implicit-mixed finite element scheme for reactive solute transport in porous media, Numer. Methods Partial Differential Equations, 26(2010), 320–344.
[40] P.A. Raviart and J.M. Thomas, A mixed finite element method for 2nd order elliptic problems, Mathematical Aspects of Finite Element Methods, Lecture Notes in Math., vol. 606, Springer-Verlag, 1977, 292–315.
[41] B.M. Rivière and N.J. Walkington, Convergence of a discontinuous Galerkin method for the miscible displacement equation under low regularity, SIAM J. Numer. Anal., 49(2011), 1085–1110.
[42] T. F. Russell, Time stepping along characteristics with incomplete iteration for a Galerkin approximation of miscible displacement in porous media, SIAM J. Numer. Anal., 22(1985), 970–1013.
[43] G. Scovazzi, M.F. Wheeler, A. Mikelić and S. Lee, Analytical and variational numerical methods for unstable miscible displacement flows in porous media, J. Comput. Phys., 335(2017), 444–496.
[44] V. Thomée, Galerkin finite element methods for parabolic problems, Springer-Verkag Berkub Geudekberg 1997.
[45] H. Wang, An optimal-order error estimate for a family of ELLAM-MFEM approximations to porous medium flow, SIAM J. Numer. Anal., 46(2008), 2133–2152.
[46] H. Wang, D. Liang, R.E. Ewing, S.L. Lyons and G. 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. Comput., 22(2000), 561–581.
[47] J. Wang, Z. Si and W. Sun, A new error analysis of characteristics-mixed FEMs for miscible displacement in porous media, SIAM J. Numer. Anal., 52(2014), 3300–3020.
[48] M.F. Wheeler, A priori $L^{2}$ error estimates for Galerkin approximations to parabolic partial differential equations, SIAM J. Numer. Anal., 10(1973), 723–759.
[49] C. Wu and W. Sun, Efficient fully-discrete Galerkin-mixed FEMs for incompressible miscible flow in porous media, in preparation.
[50] H. Zheng, J. Yu and L. Shan, Unconditional error estimates for time dependent viscoelastic fluid flow, Appl. Numer. Math., 119(2017), 1–17.

	$\displaystyle\|J_{1}(\theta_{c}^{n+1})\|$	$\displaystyle\leq C(\\|\theta_{c}^{n+1}\\|_{L^{2}}^{2}+\\|D_{t}({\cal C}^{n+1}-\Pi_{c}^{n+1}{\cal C}^{n+1})\\|_{L^{2}}^{2}),$
	$\displaystyle\|J_{4}(\theta_{c}^{n+1})\|$	$\displaystyle\leq C\\|q^{P}\\|_{L^{3}}{\\|{\cal C}^{n+1}-\Pi_{c}^{n+1}{\cal C}^{n+1}\\|_{L^{2}}}\\|\theta_{c}^{n+1}\\|_{L^{6}}$
		$\displaystyle\leq{\epsilon\\|\nabla\theta_{c}^{n+1}\\|_{H^{1}}^{2}+C\epsilon^{-1}h^{4}},$

	$\displaystyle\|J_{2}(\theta_{c}^{n+1})\|$	$\displaystyle=\|({\bf U}^{n+1}\cdot\nabla(\theta_{c}^{n}+\Pi_{c}^{n}{\cal C}^{n}-{\cal C}^{n}),\,\theta_{c}^{n+1})\|$
		$\displaystyle=\|((q^{I}-q^{P})(\theta_{c}^{n}+\Pi_{c}^{n}{\cal C}^{n}-{\cal C}^{n}),\,\theta_{c}^{n+1})$
		$\displaystyle~{}~{}~{}+(\theta_{c}^{n}+\Pi_{c}^{n}{\cal C}^{n}-{\cal C}^{n},\,{\bf U}^{n+1}\cdot\nabla\theta_{c}^{n+1})\|$
		$\displaystyle\leq\epsilon\\|\theta_{c}^{n+1}\\|_{H^{1}}^{2}+C\epsilon^{-1}\\|\theta_{c}^{n}\\|_{L^{2}}^{2}+C\epsilon^{-1}h^{4}\,.$

		$\displaystyle\|\big{(}({\bf U}^{n+1}-{\bf U}_{h}^{n+1})\cdot\nabla({\cal C}_{h}^{n}-{\cal C}^{n}),\,\theta_{c}^{n+1}\big{)}\|$
	$\displaystyle\leq$	$\displaystyle C(\\|\theta_{u}^{n+1}\\|_{L^{2}}+\\|{\widetilde{\bf U}}_{h}^{n+1}-{\bf U}_{h}^{n+1}\\|_{L^{2}})(\\|\nabla\theta_{c}^{n}\\|_{L^{3}}+\\|\nabla(\Pi_{c}^{n}{\cal C}^{n}-{\cal C}^{n})\\|_{L^{3}})\\|\theta_{c}^{n+1}\\|_{L^{6}}$
	$\displaystyle\leq$	$\displaystyle C(\\|\theta^{n}_{c}\\|_{L^{2}}+h)(h^{-d/6}\\|\nabla\theta_{c}^{n}\\|_{L^{2}}+h)\\|\theta_{c}^{n+1}\\|_{L^{6}}$
	$\displaystyle\leq$	$\displaystyle C(h+h^{-d/6}\\|\nabla\theta_{c}^{n}\\|_{L^{2}})(\\|\theta_{c}^{n}\\|_{L^{2}}\\|\theta_{c}^{n+1}\\|_{H^{1}}+h\\|\theta_{c}^{n+1}\\|_{H^{1}})\,$
	$\displaystyle\leq$	$\displaystyle(Ch^{-d/6}\\|\theta^{n}_{c}\\|_{L^{2}}+Ch+\epsilon)(\\|\theta_{c}^{n}\\|_{H^{1}}^{2}+\\|\theta_{c}^{n+1}\\|_{H^{1}}^{2})+C\epsilon^{-1}h^{4},$

		$\displaystyle\big{\|}\big{(}({\bf U}_{h}^{n+1}-{\bf U}^{n+1})\cdot\nabla{\cal C}^{n},\,\theta_{c}^{n+1}\big{)}\big{\|}$
	$\displaystyle\leq$	$\displaystyle\|\big{(}\theta_{u}^{n+1}\cdot\nabla{\cal C}^{n},\,\theta_{c}^{n+1}\big{)}\big{\|}+\|\big{(}({\widetilde{\bf U}}_{h}^{n+1}-{\bf U}^{n+1})\cdot\nabla{\cal C}^{n},\,\theta_{c}^{n+1}\big{)}\big{\|}$
	$\displaystyle\leq$	$\displaystyle C\\|\theta_{u}^{n+1}\\|_{L^{2}}\\|\theta_{c}^{n+1}\\|_{L^{2}}+C\\|{\widetilde{\bf U}}_{h}^{n+1}-{\bf U}^{n+1}\\|_{H^{-1}}\\|\nabla{\cal C}^{n}\theta_{c}^{n+1}\\|_{H^{1}}$
	$\displaystyle\leq$	$\displaystyle\epsilon\\|\theta_{c}^{n+1}\\|_{H^{1}}^{2}+C\\|\theta^{n}_{c}\\|_{L^{2}}^{2}+C\\|\theta_{c}^{n+1}\\|_{L^{2}}^{2}+C\epsilon^{-1}h^{4}.$

	$\displaystyle\frac{1}{2}D_{t}\\|\theta_{c}^{n+1}\\|^{2}_{L^{2}}+\big{\\|}\sqrt{D({\bf U}_{h}^{n+1})}\nabla\theta_{c}^{n+1}\big{\\|}^{2}_{L^{2}}$
$\displaystyle\leq$	$\displaystyle C(\epsilon+h^{-d/6}\\|\theta^{n}_{c}\\|_{L^{2}}+h)(\\|\nabla\theta_{c}^{n}\\|_{L^{2}}^{2}+\\|\nabla\theta_{c}^{n+1}\\|_{L^{2}}^{2})+C\epsilon^{-1}(\\|\theta_{c}^{n+1}\\|_{L^{2}}^{2}+\\|\theta^{n}_{c}\\|_{L^{2}}^{2})+C\epsilon^{-1}h^{4}$
	$\displaystyle~{}~{}~{}+C\big{\\|}D_{t}({\cal C}^{n+1}-\Pi_{c}^{n+1}{\cal C}^{n+1})\big{\\|}_{L^{2}}^{2}\,.$	(3.29)