# Apriori and aposteriori error estimation of Subgrid multiscale stabilized
finite element method for fully coupled Navier-Stokes
Transport model

B.V. Rathish Kumar, Manisha Chowdhury  Email addresses<EMAIL_ADDRESS>(B.V.R. Kumar); chowdhurymanisha8@gmail.com(M.Chowdhury)

(Indian Institute of Technology Kanpur
Kanpur, Uttar Pradesh, India)

###### Abstract

In this paper a fully coupled system of transient $Navier$-$Stokes$ ($NS$)
fluid flow model and variable coefficient unsteady Advection-Diffusion-
Reaction ($VADR$) transport model has been studied through subgrid multiscale
stabilized finite element method. In particular algebraic approach of
approximating the subscales has been considered to arrive at stabilized
variational formulation of the coupled system. This system is strongly coupled
since viscosity of the fluid depends upon the concentration, whose
transportation is modelled by $VADR$ equation. Fully implicit schemes have
been considered for time discretisation. Further more elaborated derivations
of both $apriori$ and $aposteriori$ error estimates for stabilized finite
element scheme have been carried out. Credibility of the stabilized method is
also established well through various numerical experiments, presented before
concluding.

Keywords Navier-Stokes equation $\cdot$ Advection-Diffusion-Reaction equation
$\cdot$ Subgrid multiscale stabilized method $\cdot$ Apriori error estimation
$\cdot$ Aposteriori error estimation

##  1 Introduction

For more than a decade transport equation coupled with fluid flow model
attracted the attention of researchers due to its significant role in
modelling various real life problems ranging from environmental issues to
physiological importance. For instance contemporary world is tackling with the
challenges of ground water pollution due to diffusion of pollutant transported
through rivers, use of drug-eluting stents to remove stenosis in human
arteries after implanting stents into them etc. Authors of
$\cite[cite]{[\@@bibref{}{RefJ}{}{}]}-\cite[cite]{[\@@bibref{}{RefN}{}{}]}$
have studied various coupled systems involving different fluid flow models and
transport equation. Whereas $Vassilev$ and $Yotov$ in [10] have presented a
mixed finite element analysis of coupled Stokes-Darcy-Transport model, $Hui$
and $Jhang$ in [13] have discussed about a stabilized mixed finite element
method for coupled transient Stokes-Darcy flows with transport.
$Cesmelio\breve{g}lu$ $et$ $al.$ have studied continuous and discontinuous
finite element methods for coupled Navier-Stokes/Darcy and transport problems
in [11] and $Cesmelio\breve{g}lu$ together with $Rivi\grave{e}re$ in [12] have
presented a mathematical analysis of existence and uniqueness of coupled NS-
Darcy- unsteady Transport equation. Recently $Chowdhury$ and $Kumar$ [14] have
considered to study subgrid scale stabilized finite element analysis of
coupled Stokes-Brinkman-Transport problem. Importantly authors in [12]-[14]
have considered the viscosity of fluid flow problem dependent on concentration
of the solute transported into the fluid. In few recent works authors in
[23]-[24] have focused on studying advection-diffusion transport equation
coupled with incompressible Navier-Stokes equation. Whereas $Du$ and $Liu$ in
[23] have worked with lattice Boltzmann model, $Yua$ $et$ $al.$ ([24]) have
studied finite difference method for the coupled model. Both of these studies
have considered constant viscosity coefficient which indicates an one-way or
weak coupling between fluid flow model and transport equation. In this paper
we have presented a stabilized finite element analysis of transient Navier-
Stokes ($NS$) fully-coupled with unsteady Advection-diffusion-reaction
equation with variable coefficients ($VADR$). $Subgrid$ $multiscale$ ($SGS$),
a most general finite element stabilization technique, has been employed to
study this coupled system. We have considered concentration dependent
viscosity in the fluid flow model as well as spatially variable diffusion
coefficients in transport equation. These considerations make this coupling
not only two-sided or strong but also more efficient to model the contemporary
real life challenges accurately. The previous studies on coupled
$NS$-Transport model have neither considered variable viscosity and diffusion
coefficients nor discussed about any error estimation for the method applied
to study the model. In this study we have elaborately carried out both
$apriori$ and $aposteriori$ error estimations for a general finite element
stabilization scheme to study strongly coupled transient $NS$-$VADR$ model.
It is well known fact that lack of stability in standard Galerkin finite
element method has driven researchers to introduce stabilized methods such as
Streamline upwind/Petrov-Galerkin ($SUPG$) formulation, Galerkin/least-squares
($G$ $LS$) method, characteristic-based split ($CBS$) method, Subgrid Scale
($SGS$) method, bubble stabilization etc. Over last four decades huge
developments have been taken place in the study of various stabilization
techniques. Starting with the works of $Brooks$ and $Hughes$ [1] on $SUPG$;
$Hughes$, $Franca$ and $Hulbert$ [2] on $GLS$; $Hughes$ introducing $SGS$ in
[3]; $Russo$ [8] explaining bubble stabilization method for the linearized
incompressible $NS$ equations, the stabilization schemes have been growing
through the studies of $Hannani$ $et$ $al.$ [4] on comparison between $SUPG$
and $GLS$ formulation for steady state incompressible $NS$ equations; $Codina$
and $Zienkiewicz$ [9] on comparison of $CBS$ and $GLS$; $Codina$ $et$ $al.$ on
numerical comparison of $CBS$ and $SGS$ method of the incompressible $NS$
equations in [5]; $Russo$ [7] on comparison of $SUPG$ and residual-free
bubbles($RFB$); $Kirk$ and $Carey$ [6] on development and validation of $SUPG$
for compressible $NS$ equations etc. $Codina$ in [15] has experimentally
established that for solving $ADR$ equation $SGS$ method performs better than
other stabilized methods, such as $SUPG$, $GLS$, $Taylor$-$Galerkin$ etc. and
in fact it is the most general method amongst them.
Generally two approaches of $SGS$ stabilized formulation, namely algebraic
approach, abbreviated as $ASGS$ and orthogonal projection approach, known as
$OSGS$ method, have been studied. Though few studies [25]\- [27] are there
applying only general form of $SGS$ method instead of working with one
specific approach, but authors in [28]-[31] have employed $ASGS$ method,
whereas authors of [32]\- [36] have worked with $OSGS$ method. Again $Badia$
and $Codina$ in [22] have studied both the approaches for unified Stokes-Darcy
fluid flow problem and experimentally established equally well performances of
both the stabilized formulations. In this paper we have considered algebraic
approach of approximating subscales which implies the stabilization parameters
are of algebraic forms. This stabilization method begins with division of weak
solution space into the spaces of the known finite element space and an
unknown subgrid or unresolvable scale space and finally the stabilized
formulation has been reached through expressing the element of subgrid scales
in terms of the element of resolvable finite element space. For time
discretization fully implicit schemes have been chosen. A detailed derivation
of $apriori$ error estimation has been carried out for this stabilized
variational form of the coupled system. Further more residual based
$aposteriori$ error estimate too is derived elaborately. First order
convergence in space has been established with respect to complete norms on
all the variables. This paper also establishes the accuracy of the stabilized
method through various numerical results, which include all possible
combinations of cases containing small and large Reynolds numbers as well as
cases involving concentration dependent viscosity. In every numerical example
$ASGS$ performs consistently well in compared to standard Galerkin finite
element method.
This paper is organised as follows: Section 2 introduces the coupled system
along with important assumptions. In the next section we have introduced weak
formulation, stabilized formulation and finally fully-discrete formulation
after applying time-discretization rule. This section also contains stability
analysis of the fully-discrete stabilized form. Section 4 has elaborately
described the derivations of $apriori$ and $aposteriori$ error estimations for
this stabilized formulation. Before concluding the article section 5 presents
numerical results to verify accuracy of the method.

##  2 Model problem

In this section we introduce the flow problem described through transient
$Navier$-$Stokes$ equations coupled with unsteady transport model over an open
bounded domain $\Omega\subset$ $\mathbb{R}^{d}$, d=2,3 with piece-wise smooth
boundary $\partial\Omega$. Let us first present the $Navier$-$Stokes$ fluid
flow model in the following: Find velocity u: $\Omega$ $\times$ (0,T)
$\rightarrow\mathbb{R}^{d}$ and pressure $p$: $\Omega\times$ (0,T)
$\rightarrow\mathbb{R}$ of the fluid such that,

$\begin{split}\rho\frac{\partial\textbf{u}}{\partial
t}+\rho(\textbf{u}\cdot\bigtriangledown)\textbf{u}-\mu(c)\Delta\textbf{u}+\bigtriangledown
p&=\textbf{f}\hskip 5.69054ptin\hskip 5.69054pt\Omega\times(0,T)\\\
\bigtriangledown\cdot\textbf{u}&=\textbf{0}\hskip 5.69054ptin\hskip
5.69054pt\Omega\times(0,T)\\\ \textbf{u}&=\textbf{0}\hskip 5.69054pton\hskip
5.69054pt\partial\Omega\times(0,T)\\\ \textbf{u}&=\textbf{u}_{0}\hskip
5.69054ptat\hskip 5.69054ptt=0\\\ \end{split}$ (1)

where $\rho$ is the density of the fluid, $\mu(c)$ is the dynamic viscosity of
the fluid depending on concentration $c$ of the dispersing mass of the solute,
f is body force and $\textbf{u}_{0}$ is the initial velocity.
This above flow problem is fully-coupled with the following transient
advection-diffusion-reaction equation with variable coefficients($VADR$),
representing the transportation of solute in $\Omega$.
Find the concentration $c$: $\Omega\times$ (0,T) $\rightarrow\mathbb{R}$ of
the solute such that,

$\begin{split}\frac{\partial c}{\partial
t}-\bigtriangledown\cdot\tilde{\bigtriangledown}c+\textbf{u}\cdot\bigtriangledown
c+\alpha c&=g\hskip 5.69054ptin\hskip 5.69054pt\Omega\times(0,T)\\\ c&=0\hskip
5.69054pton\hskip 5.69054pt\partial\Omega\times(0,T)\\\ c&=c_{0}\hskip
5.69054ptat\hskip 5.69054ptt=0\\\ \end{split}$ (2)

where the notation,
$\tilde{\bigtriangledown}:=\sum_{i=1}^{d}D_{i}\frac{\partial}{\partial
x_{i}}e_{i}$ for $d=2,3$ and $\\{e_{i}\\}_{i=1}^{d}$ is standard basis of
$\mathbb{R}^{d}$. $D_{i}$ are variable diffusion coefficients, $\alpha$ is the
reaction coefficient and $g$ denotes the source of solute mass and $c_{0}$ is
the initial concentration of the solute.
Let us consider a notation U= (u,p,c) and the system of equations can be
written in the following operator form,

$M\partial_{t}\textbf{U}+\mathcal{L}(\textbf{u};\textbf{U})=\textbf{F}$ (3)

where M, a matrix = diag($\rho$,$\rho$,0,1),
$\partial_{t}\textbf{U}=(\frac{\partial\textbf{u}}{\partial t},\frac{\partial
p}{\partial t},\frac{\partial c}{\partial t})^{T}$

$\mathcal{L}(\textbf{u};\textbf{U})=\begin{bmatrix}\rho(\textbf{u}\cdot\bigtriangledown)\textbf{u}-\mu(c)\Delta\textbf{u}+\bigtriangledown
p\\\ \bigtriangledown\cdot\textbf{u}\\\
-\bigtriangledown\cdot\tilde{\bigtriangledown}c+\textbf{u}\cdot\bigtriangledown
c+\alpha c\end{bmatrix}$

and

$\textbf{F}=\begin{bmatrix}\textbf{f}\\\ 0\\\ g\end{bmatrix}$

Let us introduce the adjoint $\mathcal{L}^{*}$ of $\mathcal{L}$ as follows,

$\mathcal{L}^{*}(\textbf{u};\textbf{U})=\begin{bmatrix}-\rho(\textbf{u}\cdot\bigtriangledown)\textbf{u}-\mu(c)\Delta\textbf{u}-\bigtriangledown
p\\\ -\bigtriangledown\cdot\textbf{u}\\\
-\bigtriangledown\cdot\tilde{\bigtriangledown}c-\textbf{u}\cdot\bigtriangledown
c+\alpha c\end{bmatrix}$

Now we assume suitable conditions on the coefficients mentioned above, which
will be useful to conclude the results further.
(i) The fluid viscosity $\mu(c)=\mu\in C^{0}(\mathbb{R}^{+};\mathbb{R}^{+})$,
the space of positive real valued functions defined on positive real numbers
and we will have two positive real numbers $\mu_{l}$ and $\mu_{u}$ such that

$0<\mu_{l}\leq\mu(x)\leq\mu_{u}\hskip 5.69054ptfor\hskip 5.69054ptany\hskip
5.69054ptx\in\mathbb{R}^{+}$ (4)

(ii) $D_{i}=D_{i}(\textbf{x},t)\in
C^{0}(\mathbb{R}^{d}\times(0,T);\mathbb{R})$ (for $i=1,...,d$) where
$C^{0}(\mathbb{R}^{d}\times(0,T);\mathbb{R})$ is the space of real valued
continuous function defined on $\mathbb{R}^{d}$ for fixed $t\in(0,T)$. Both
are bounded quantity that is we can find lower and upper bounds for both of
them.
(iii) $\rho$ and $\alpha$ are positive constants.
(iv) The spaces of continuous solution $(\textbf{u},p,c)$ are assumed as:
$\textbf{u}\in L^{\infty}(0,T;(H^{2}(\Omega))^{d})\bigcap
C^{0}(0,T;H_{0}^{1}(\Omega))$ and
$p\in L^{\infty}(0,T;H^{1}(\Omega))\bigcap C^{0}(0,T;L^{2}_{0}(\Omega))$,
$c\in L^{\infty}(0,T;H^{2}(\Omega))\bigcap C^{0}(0,T;H^{1}_{0}(\Omega))$
(v) Additional assumptions on continuous velocity and concentration solutions
are: $\textbf{u}_{tt},\textbf{u}_{ttt},c_{tt},c_{ttt}$ all are taken to be
bounded functions on $\Omega$ for each $t\in(0,T)$.
Weak formulation: Assuming the body force $\textbf{f}\in
L^{2}(0,T;(L^{2}(\Omega))^{d})$ and the source term $g\in
L^{2}(0,T;L^{2}(\Omega))$ let us consider the spaces suitable to define the
weak form as $V=H^{1}_{0}(\Omega)$ and $Q=L^{2}(\Omega)$ and J= (0,T).
Now denoting the space $V^{d}\times Q\times V$ (for $d=2,3$) by
$\textbf{V}_{F}$ the weak formulation of (3) is to find U= (u,p,c): J
$\rightarrow\textbf{V}_{F}$ such that $\forall$ V=(v,q,d) $\in\textbf{V}_{F}$

$\begin{split}(M\partial_{t}\textbf{U},\textbf{V})+B(\textbf{u};\textbf{U},\textbf{V})=L(\textbf{V})\end{split}$
(5)

where
$B(\textbf{u};\textbf{U},\textbf{V})=c(\textbf{u},\textbf{u},\textbf{v})+a_{NS}(\textbf{u},\textbf{v})-b(\textbf{v},p)+b(\textbf{u},q)+a_{T}(c,d)$
$(M\partial_{t}\textbf{U},\textbf{V})=\rho\int_{\Omega}\frac{\partial\textbf{u}}{\partial
t}\cdot\textbf{v}+\int_{\Omega}\frac{\partial c}{\partial t}d$ and
$L(\textbf{V})=l_{NS}(\textbf{v})+l_{T}(d)$
where the notations are defined in the following:
$c(\textbf{u},\textbf{v},\textbf{w})=\rho\int_{\Omega}((\textbf{u}\cdot\bigtriangledown)\textbf{v})\cdot\textbf{w}+\frac{\rho}{2}\int_{\Omega}(\bigtriangledown\cdot\textbf{u})\textbf{v}\cdot\textbf{w}$
$a_{NS}(\textbf{u},\textbf{v})=\int_{\Omega}\mu(c)\bigtriangledown\textbf{u}:\bigtriangledown\textbf{v}$
$b(\textbf{v},q)=\int_{\Omega}(\bigtriangledown\cdot\textbf{v})q$
$a_{T}(c,d)=\int_{\Omega}\tilde{\bigtriangledown}c\cdot\bigtriangledown
d+\int_{\Omega}d\textbf{u}\cdot\bigtriangledown c+\alpha\int_{\Omega}cd$
$l_{NS}(\textbf{v})=\int_{\Omega}\textbf{f}\cdot\textbf{v}$ and
$l_{T}(d)=\int_{\Omega}gd$
The modified trilinear form $c(\cdot,\cdot,\cdot)$ considered here is
equivalent to the original trilinear form obtained from the non-linear
convective term in (1). By the virtue of this modified form the trilinear term
introduces the following important property.
(a) For any $\textbf{u}\in V^{d}$, $c(\textbf{u},\textbf{v},\textbf{v})=0$
$\forall\textbf{v}\in V^{d}$
Besides the trilinear form $c(\cdot,\cdot,\cdot)$ has the another property
[18] too.
(b) For any $\textbf{u},\textbf{v},\textbf{w}$ $\in V^{d}$

$c(\textbf{u},\textbf{v},\textbf{w})\leq\begin{cases}C\|\textbf{u}\|_{1}\|\textbf{v}\|_{1}\|\textbf{w}\|_{1}&\\\
C\|\textbf{u}\|_{0}\|\textbf{v}\|_{2}\|\textbf{w}\|_{1}&\\\
C\|\textbf{u}\|_{2}\|\textbf{v}\|_{1}\|\textbf{w}\|_{0}&\\\
C\|\textbf{u}\|_{0}\|\textbf{v}\|_{1}\|\textbf{w}\|_{L^{\infty}(\Omega)}\end{cases}$
(6)

where $C$ is a constant and $\|\cdot\|_{j}$ for $j$=0,1,2 denote standard
$L^{2},H^{1},H^{2}$ full norms respectively. In the other forms the bilinear
forms $a_{NS}(\cdot,\cdot)$ is $coercive$ [18] and $a_{T}(\cdot,\cdot)$ is
also $continuous$ and $coercive$ [19]. Again $b(\cdot,\cdot)$ satisfies
$inf$-$sup$ condition [18] too.

##  3 Discrete formulation

###  3.1 Semi-discrete formulation

In this section we introduce the finite element space discretization for the
variational formulation (5) followed by a stabilized finite element
formulation for the same.
Let the domain $\Omega$ be discretized into finite numbers of subdomains
$\Omega_{k}$ for k=1,2,…,$n_{el}$, where $n_{el}$ is the total number element
subdomains. Let $h_{k}$ be the diameter of each subdomain $\Omega_{k}$ and h=
$\underset{k=1,2,...n_{el}}{max}h_{k}$
Let $\tilde{\Omega}=\bigcup_{k=1}^{n_{el}}\Omega_{k}$ be the union of interior
elements.
Let $V^{h}=\\{v\in V:v(\Omega_{k})=\mathcal{P}^{k}(\Omega_{k})\\}$ and
$Q^{h}=\\{q\in Q:q(\Omega_{k})=\mathcal{P}^{k}(\Omega_{k})\\}$
where $V^{h}$ and $Q^{h}$ be finite dimensional subspaces of $V$ and $Q$
respectively and $\mathcal{P}^{k}(\Omega_{k})$ denotes complete polynomial of
order $k$ over each $\Omega_{k}$ for k=1,2,…,$n_{el}$. For regular partitions
the functions belonging to finite dimensional spaces satisfy the following
inverse inequalities:

$\|\Delta v_{h}\|\leq C_{I}h^{-1}\|\bigtriangledown v_{h}\|_{0,k}$ and
$\|\bigtriangledown v_{h}\|_{0,k}\leq C_{I}h^{-1}\|v_{h}\|_{0,k}$

Considering similar notation $\textbf{V}_{F}^{h}$, denoting
$\textbf{V}_{F}^{h}=(V^{h})^{d}\times Q^{h}\times V^{h}$ the finite element
formulation of the variational form (5) in the finite dimensional space
$\textbf{V}_{F}^{h}$ is to find $\textbf{U}_{h}$=
$(\textbf{u}_{h},p_{h},c_{h})$: J $\rightarrow\textbf{V}_{F}^{h}$ such that
$\forall$ $\textbf{V}_{h}=(\textbf{v}_{h},q_{h},d_{h})$
$\in\textbf{V}_{F}^{h}$

$(M\partial_{t}\textbf{U}_{h},\textbf{V}_{h})+B(\textbf{u}_{h};\textbf{U}_{h},\textbf{V}_{h})=L(\textbf{V}_{h})$
(7)

where $(M\partial_{t}\textbf{U}_{h},\textbf{V}_{h})$=
$\rho(\frac{\partial\textbf{u}_{h}}{\partial
t},\textbf{v}_{h})+(\frac{\partial c_{h}}{\partial t},d_{h})$
$B(\textbf{u}_{h};\textbf{U}_{h},\textbf{V}_{h})$ =
$c(\textbf{u}_{h},\textbf{u}_{h},\textbf{v}_{h})+a_{NS}(\textbf{u}_{h},\textbf{v}_{h})-b(\textbf{v}_{h},p_{h})+b(\textbf{u}_{h},q_{h})+a_{T}(c_{h},d_{h})$
and $L(\textbf{V}_{h})=l_{NS}(\textbf{v}_{h})+l_{T}(d_{h})$
In addition let us consider the initial conditions
$(\textbf{u}_{h},\textbf{v}_{h})\mid_{t=0}=(\textbf{u}_{0},\textbf{v}_{h})$
$\forall\textbf{v}_{h}\in(V^{h})^{d}$ and
$(c_{h},d_{h})\mid_{t=0}=(c_{0},d_{h})$ $\forall d_{h}\in V^{h}$.
Now we are going to introduce $subgrid$ $multiscale$ stabilized finite element
method with algebraic approximation of the subscales of (5). It involves
decomposition of the weak solution space $\textbf{V}_{F}$ into the spaces of
resolvable scales and unresolvable or subgrid scales. The finite element space
$\textbf{V}_{F}^{h}$ is chosen to be the space of resolvable scales and in
literature one of the ways of choosing the space of subgrid scales is the
space that completes $\textbf{V}_{F}^{h}$ in $\textbf{V}_{F}$. Then the final
form of subgrid formulation will be arrived while the elements of subgrid
scales will be expressed in the terms of elements of resolvable scales.
The stabilized algebraic subgrid multiscale ($ASGS$) formulation for this
coupled equation to find $\textbf{U}_{h}$= $(\textbf{u}_{h},p_{h},c_{h})$: J
$\rightarrow\textbf{V}_{F}^{h}$ such that $\forall$
$\textbf{V}_{h}=(\textbf{v}_{h},q_{h},d_{h})$ $\in\textbf{V}_{F}^{h}$

$(M\partial_{t}\textbf{U}_{h},\textbf{V}_{h})+B_{ASGS}(\textbf{u}_{h};\textbf{U}_{h},\textbf{V}_{h})=L_{ASGS}(\textbf{V}_{h})$
(8)

where
$B_{ASGS}(\textbf{u}_{h};\textbf{U}_{h},\textbf{V}_{h})=B(\textbf{u}_{h};\textbf{U}_{h},\textbf{V}_{h})+\sum_{k=1}^{n_{el}}(\tau_{k}^{\prime}(M\partial_{t}\textbf{U}_{h}+\mathcal{L}(\textbf{u}_{h};\textbf{U}_{h})-\textbf{d}),-\mathcal{L}^{*}(\textbf{u}_{h};\textbf{V}_{h}))_{\Omega_{k}}-\sum_{k=1}^{n_{el}}((I-\tau_{k}^{-1}\tau_{k}^{\prime})(M\partial_{t}\textbf{U}_{h}+\mathcal{L}(\textbf{u}_{h};\textbf{U}_{h})),\textbf{V}_{h})_{\Omega_{k}}\\\
-\sum_{k=1}^{n_{el}}(\tau_{k}^{-1}\tau_{k}^{\prime}\textbf{d},\textbf{V}_{h})_{\Omega_{k}}$
$L_{ASGS}(\textbf{V}_{h})=L(\textbf{V}_{h})+\sum_{k=1}^{n_{el}}(\tau_{k}^{\prime}\textbf{F},-\mathcal{L}^{*}(\textbf{u}_{h};\textbf{V}_{h}))_{\Omega_{k}}-\sum_{k=1}^{n_{el}}((I-\tau_{k}^{-1}\tau_{k}^{\prime})\textbf{F},\textbf{V}_{h})_{\Omega_{k}}$
where the stabilization parameter $\tau_{k}$ is in matrix form as

$\tau_{k}=diag(\tau_{1k},\tau_{1k},\tau_{2k},\tau_{3k})=\begin{bmatrix}\tau_{1k}I_{d\times
d}&0&0\\\ 0&\tau_{2k}&0\\\ 0&0&\tau_{3k}\end{bmatrix}$

and

$\tau_{k}^{\prime}=(\frac{1}{dt}M+\tau_{k}^{-1})^{-1}=\begin{bmatrix}\frac{\tau_{1k}dt}{dt+\rho\tau_{1k}}I_{d\times
d}&0&0\\\ 0&\tau_{2k}&0\\\
0&0&\frac{\tau_{3k}dt}{dt+\tau_{3k}}\end{bmatrix}\\\
=diag(\tau_{1k}^{\prime},\tau_{1k}^{\prime},\tau_{2k}^{\prime},\tau_{3k}^{\prime})$

$I_{d\times d}$ is an identity matrix for $d=2,3$.
d=
$\sum_{i=1}^{n+1}(\frac{1}{dt}M\tau_{k}^{\prime})^{i}(\textbf{F}-M\partial_{t}\textbf{U}_{h}-\mathcal{L}(\textbf{u}_{h};\textbf{U}_{h}))$
=$[\textbf{d}_{1},d_{2},d_{3}]^{T}$
It can be easily observed that $d_{2}$ is always 0 due to the matrix M.
We have the forms of the stabilization parameters $\tau_{1k},\tau_{2k}$ for
$Navier$-$Stokes$ equation in [21] and $\tau_{3k}$ for $VADR$ equation [20]
and for each k=1,2,…,$n_{el}$ all the coefficients $\tau_{ik}$ coincide with
$\tau_{i}$ for i=1,2,3 and choosing the parameters $c_{1},c_{2},c_{3}$
suitably that $\tau_{i}$’s are as follows:

$\begin{split}\tau_{1k}&=\tau_{1}=(c_{1}\frac{\mu_{\textbf{u}}}{h^{2}}+c_{2}\frac{\rho\|\textbf{u}_{h}\|}{h})^{-1}\\\
\tau_{2k}&=\tau_{2}=\frac{h^{2}}{c_{1}\tau_{1}}\\\
\tau_{3k}&=\tau_{3}=c_{3}(\frac{9D}{4h^{2}}+\frac{3\|\textbf{u}_{h}\|}{2h}+\alpha)^{-1}\end{split}$
(9)

where $\textbf{u}_{h}$ is the computed velocity.

######  Remark 1.

Considering continuity of the solutions at the inter-element boundaries, we
have not encountered with any jump term in the above stabilized formulation.

###  3.2 Fully-discrete formulation

Before introducing time discretization, some notations have been introduced:
for $dt$= $\frac{T}{N}$, where $N$ is a positive integer, $t_{n}=ndt$ and for
given $0\leq\theta\leq 1$,

$\begin{split}f^{n}&=f(\cdot,t_{n})\hskip 11.38109ptfor\hskip 5.69054pt0\leq
n\leq N\\\
f^{n,\theta}&=\frac{1}{2}(1+\theta)f^{(n+1)}+\frac{1}{2}(1-\theta)f^{n}\hskip
11.38109ptfor\hskip 5.69054pt0\leq n\leq N-1\end{split}$ (10)

Later we will see for $\theta=0$ the discretization follows Crank-Nicolson
formula and for $\theta=1$ it is backward Euler discretization rule.
For sufficiently smooth function $f(t)$, using the Taylor series expansion
about t= $t^{n,\theta}$, we will have

$\begin{split}f^{n+1}&=f(t^{n,\theta})+\frac{(1-\theta)dt}{2}\frac{\partial
f}{\partial
t}(t^{n,\theta})+\frac{(1-\theta)^{2}dt^{2}}{8}\frac{\partial^{2}f}{\partial
t^{2}}(t^{n,\theta})+\mathcal{O}(dt^{3})\\\
f^{n}&=f(t^{n,\theta})-\frac{(1+\theta)dt}{2}\frac{\partial f}{\partial
t}(t^{n,\theta})+\frac{(1+\theta)^{2}dt^{2}}{8}\frac{\partial^{2}f}{\partial
t^{2}}(t^{n,\theta})+\mathcal{O}(dt^{3})\end{split}$ (11)

We have considered here $t^{n,\theta}-t^{n}=\frac{(1+\theta)\Delta t}{2}$
Multiplying the above first and second sub-equations in (14) by
$\frac{1+\theta}{2}$ and $\frac{1-\theta}{2}$ respectively and then adding
them we will have the following

$f^{n,\theta}=f(t^{n,\theta})+\frac{1}{8}(1+\theta)(1-\theta)dt^{2}\frac{\partial^{2}f}{\partial
t^{2}}(t^{n,\theta})+\mathcal{O}(dt^{3})$ (12)

Let $\textbf{u}^{n,\theta},p^{n,\theta},c^{n,\theta}$ be approximations of
$\textbf{u}(\textbf{x},t^{n,\theta}),p(\textbf{x},t^{n,\theta}),c(\textbf{x},t^{n,\theta})$
respectively. Now by Taylor series expansion [17],we have

$\begin{split}\frac{\textbf{u}^{n+1}-\textbf{u}^{n}}{dt}&=\frac{\partial\textbf{u}}{\partial
t}(\textbf{x},t^{n,\theta})+\textbf{TE}_{1}\mid_{t=t^{n,\theta}}\hskip
11.38109pt\forall\textbf{x}\in\Omega\\\
\frac{c^{n+1}-c^{n}}{dt}&=\frac{\partial c}{\partial
t}(\textbf{x},t^{n,\theta})+TE_{2}\mid_{t=t^{n,\theta}}\hskip
14.22636pt\forall\textbf{x}\in\Omega\end{split}$ (13)

where the truncation error $TE\mid_{t=t^{n,\theta}}$ $\simeq$ $TE^{n,\theta}$
depends upon time-derivatives of the respective variables and $dt$.

$\begin{split}\|\textbf{TE}_{1}^{n,\theta}\|&\leq\begin{cases}C^{\prime}dt\|\textbf{u}_{tt}^{n,\theta}\|_{L^{\infty}(t^{n},t^{n+1},L^{2})}&if\hskip
2.84526pt\theta=1\\\
C^{\prime\prime}dt^{2}\|\textbf{u}_{ttt}^{n,\theta}\|_{L^{\infty}(t^{n},t^{n+1},L^{2})}&if\hskip
2.84526pt\theta=0\end{cases}\end{split}$ (14)

The above relation holds for $TE_{2}$ in similar manner. Now applying
assumption (v) on $\textbf{u}_{tt}$ and $\textbf{u}_{ttt}$ we will have
another property as follows:

$\begin{split}\|\textbf{TE}_{1}^{n,\theta}\|&\leq\begin{cases}C^{\prime}dt&if\hskip
2.84526pt\theta=1\\\ C^{\prime\prime}dt^{2}&if\hskip
2.84526pt\theta=0\end{cases}\end{split}$ (15)

Similarly

$\begin{split}\|TE_{2}^{n,\theta}\|&\leq\begin{cases}C^{\prime}dt&if\hskip
2.84526pt\theta=1\\\ C^{\prime\prime}dt^{2}&if\hskip
2.84526pt\theta=0\end{cases}\end{split}$ (16)

After introducing all the required definitions finally the fully-discrete
formulation of $subgrid$ form is as follows: For given
$\textbf{U}_{h}^{n}=(\textbf{u}_{h}^{n},p_{h}^{n},c_{h}^{n})\in\textbf{V}_{F}^{h}$
find
$\textbf{U}_{h}^{n+1}=(\textbf{u}_{h}^{n+1},p_{h}^{n+1},c_{h}^{n+1})\in\textbf{V}_{F}^{h}$
such that , $\forall\hskip
2.84526pt\textbf{V}_{h}=(\textbf{v}_{h},q_{h},d_{h})\in\textbf{V}_{F}^{h}$

$(M\frac{(\textbf{U}_{h}^{n+1}-\textbf{U}_{h}^{n})}{dt},\textbf{V}_{h})+B_{ASGS}(\textbf{u}_{h}^{n};\textbf{U}_{h}^{n,\theta},\textbf{V}_{h})=L_{ASGS}(\textbf{V}_{h})+(\textbf{TE}^{n,\theta},\textbf{V}_{h})$
(17)

Again for the exact solution we will have the discrete formulation as follows:
For given $\textbf{U}^{n}=(\textbf{u}^{n},p^{n},c^{n})\in\textbf{V}_{F}$ find
$\textbf{U}^{n+1}=(\textbf{u}^{n+1},p^{n+1},c^{n+1})\in\textbf{V}_{F}$ such
that , $\forall\hskip
2.84526pt\textbf{V}_{h}=(\textbf{v}_{h},q_{h},d_{h})\in\textbf{V}_{F}^{h}$

$(M\frac{(\textbf{U}^{n+1}-\textbf{U}^{n})}{dt},\textbf{V}_{h})+B(\textbf{u}^{n};\textbf{U}^{n,\theta},\textbf{V}_{h})=L(\textbf{V}_{h})+(\textbf{TE}^{n,\theta},\textbf{V}_{h})$
(18)

##  4 Error estimates

We start this section with introducing the projection operator corresponding
to each unknown variable followed by notation of error and it’s component wise
splitting. Later we go to derive $apriori$ and $aposteriori$ error estimates.

###  4.1 Projection operators : Error splitting

Let us introduce the projection operator for each of these error components.
(i)For any $\textbf{u}\in(H^{2}(\Omega))^{d}$ we assume that there exists an
interpolation
$I^{h}_{\textbf{u}}:(H^{2}(\Omega))^{d}\longrightarrow(V^{h})^{d}$ satisfying
$b(\textbf{u}-I^{h}_{\textbf{u}}\textbf{u},q_{h})=0$ $\forall q_{h}\in Q^{h}$
(ii) Let $I^{h}_{p}:H^{1}(\Omega)\longrightarrow Q^{h}$ be the $L^{2}$
orthogonal projection given by
$\int_{\Omega}(p-I^{h}_{p}p)q_{h}=0$ $\forall q_{h}\in Q^{h}$ and for any
$p\in H^{1}(\Omega)$
(iii)Similarly let $I^{h}_{c}:H^{2}(\Omega)\longrightarrow V^{h}$ be the
$L^{2}$ orthogonal projection given by
$\int_{\Omega}(c-I^{h}_{c}c)d_{h}=0$ $\forall d_{h}\in V^{h}$ and for any
$c\in H^{2}(\Omega)$
Let $\textbf{e}=(e_{\textbf{u}},e_{p},e_{c})$ denote the error where the
components are $e_{\textbf{u}}=\textbf{u}-\textbf{u}_{h},e_{p}=(p-p_{h})$ and
$e_{c}=(c-c_{h})$. Now each component of the error can be split into two parts
interpolation part, $E^{I}$ and auxiliary part, $E^{A}$ as follows:
$e_{\textbf{u}}=(\textbf{u}-\textbf{u}_{h})=(\textbf{u}-I^{h}_{\textbf{u}}\textbf{u})+(I^{h}_{\textbf{u}}\textbf{u}-\textbf{u}_{h})=E^{I}_{\textbf{u}}+E^{A}_{\textbf{u}}$
Similarly $e_{p}=E^{I}_{p}+E^{A}_{p}$, and $e_{c}=E^{I}_{c}+E^{A}_{c}$
At this point let us mention the standard interpolation estimation result [17]
in the following: for any exact solution with regularity upto (m+1)

$\|v-I^{h}_{v}v\|_{l}=\|E^{I}_{v}\|_{l}\leq C(p,\Omega)h^{m+1-l}\|v\|_{m+1}$
(19)

where l ($\leq m+1$) is a positive integer and C is a constant depending on m
and the domain. For l=0 and 1 it implies standard $L^{2}(\Omega)$ and
$H^{1}(\Omega)$ norms respectively. For simplicity we will use $\|\cdot\|$
instead of $\|\cdot\|_{0}$ to denote $L^{2}(\Omega)$ norm. Now we put some
results using the properties of projection operators and these results will be
used in error estimations.

######  Result 1.

$(\frac{\partial}{\partial t}E^{I,n},v_{h})=0\hskip 5.69054ptv_{h}\in V^{h}$
(20)

######  Result 2.

For any given auxiliary error $E^{A,n}$ and unknown $E^{A,n+1}$

$(\frac{\partial}{\partial
t}E^{A,n},E^{A,n,\theta})\geq\frac{1}{2dt}(\|E^{A,n+1}\|^{2}-\|E^{A,n}\|^{2})$
(21)

######  Remark 2.

The proof of the results have been discussed in [14] elaborately.

###  4.2 Apriori error estimate

In this section we will find $apriori$ error bound, which depends on the exact
solution. Here we first estimate $auxiliary$ error bound and later using that
we will find $apriori$ error estimate. Before deriving error estimations let
us define norms required for error estimations. Let us consider the space
$\tilde{\textbf{V}}$ := $L^{2}(0,T;V_{s})\bigcap L^{\infty}(0,T;Q_{s})$ and
it’s associated norm is denoted by $\tilde{\textbf{V}}$-norm. For the
functions $g_{1},g_{2},g_{3}$ belonging to the spaces
$L^{2}(0,T;L^{2}(\Omega))$, $L^{2}(0,T;H_{0}^{1}(\Omega))$,
$\tilde{\textbf{V}}$ respectively norms over these spaces, abbreviated as
$L^{2}(L^{2})$, $L^{2}(H^{1})$, $\tilde{\textbf{V}}$ are defined in the
following

$\begin{split}\|g_{1}\|_{L^{2}(L^{2})}^{2}&=\sum_{n=0}^{N-1}\int_{t^{n}}^{t^{n+1}}\int_{\Omega}\mid
g_{1}^{n,\theta}\mid^{2}dt\\\
\|g_{2}\|_{L^{2}(H^{1})}^{2}&=\sum_{n=0}^{N-1}\int_{t^{n}}^{t^{n+1}}(\int_{\Omega}\mid
g_{2}^{n,\theta}\mid^{2}+\int_{\Omega}\mid\frac{\partial g_{2}}{\partial
x}^{n,\theta}\mid^{2}+\int_{\Omega}\mid\frac{\partial g_{2}}{\partial
y}^{n,\theta}\mid^{2})dt\\\
\|g_{3}\|_{\tilde{\textbf{V}}}^{2}&=\underset{0\leq n\leq
N}{max}\|g_{3}^{n}\|^{2}+\|g_{3}\|_{L^{2}(H^{1})}^{2}\\\ \end{split}$ (22)

######  Theorem 1.

(Auxiliary error estimate) For computed velocity $\textbf{u}_{h}$, pressure
$p_{h}$ and concentration $c_{h}$ belonging to $(V^{h})^{d}\times Q^{h}\times
V^{h}$ satisfying (31)-(32), assume $dt$ is sufficiently small and positive,
and sufficient regularity of exact solution in equations (1)-(2). Then there
exists a constant C, depending upon $\textbf{u},p,c$ such that

$\|E^{A}_{\textbf{u}}\|^{2}_{\tilde{\textbf{V}}}+\|E^{A}_{p}\|_{L^{2}(L^{2})}^{2}+\|E^{A}_{c}\|^{2}_{\tilde{\textbf{V}}}\leq
C(h^{2}+dt^{2r})$ (23)

where

$r=\begin{cases}1,&\text{if}\ \theta=1\\\ 2,&\text{if}\ \theta=0\end{cases}$
(24)

###### Proof.

In first part we will find bound for auxiliary error part of velocity u and
concentration c with respect to ${\tilde{\textbf{V}}}$-norm and in the second
part we will estimate auxiliary error for pressure term with respect to $Q$
norm and finally combining them we will arrive at the desired result.
First part Subtracting (17) from (18) and then simplifying the terms, we have
$\forall\hskip 2.84526pt\textbf{V}_{h}\in(V^{h})^{d}\times Q^{h}\times V^{h}$

$(M\frac{(\textbf{U}^{n+1}-\textbf{U}^{n+1}_{h})-(\textbf{U}^{n}-\textbf{U}^{n}_{h})}{dt},\textbf{V}_{h})+B(\textbf{u}^{n};\textbf{U}^{n,\theta},\textbf{V}_{h})-B(\textbf{u}_{h}^{n};\textbf{U}^{n,\theta}_{h},\textbf{V}_{h})\\\
+\sum_{k=1}^{n_{el}}(\tau_{k}^{\prime}(M\partial_{t}(\textbf{U}^{n}-\textbf{U}^{n}_{h})+\mathcal{L}(\textbf{u}^{n};\textbf{U}^{n,\theta})-\mathcal{L}(\textbf{u}^{n}_{h};\textbf{U}^{n,\theta}_{h})),-\mathcal{L}^{*}(\textbf{u}_{h};\textbf{V}_{h}))_{\Omega_{k}}\\\
+\sum_{k=1}^{n_{el}}((I-\tau_{k}^{-1}\tau_{k}^{\prime})(M\partial_{t}(\textbf{U}^{n}-\textbf{U}^{n}_{h})+\mathcal{L}(\textbf{u}^{n};\textbf{U}^{n,\theta})-\mathcal{L}(\textbf{u}^{n}_{h};\textbf{U}^{n,\theta}_{h})),-\textbf{V}_{h})_{\Omega_{k}}\\\
+\sum_{k=1}^{n_{el}}(\tau_{k}^{-1}\tau_{k}^{\prime}\textbf{d},\textbf{V}_{h})_{\Omega_{k}}+\sum_{k=1}^{n_{el}}(\tau_{k}^{\prime}\textbf{d},-\mathcal{L}^{*}(\textbf{u}_{h};\textbf{V}_{h}))_{\Omega_{k}}=(\textbf{TE}^{n,\theta},\textbf{V}_{h})\\\
$ (25)

where d=
$(\sum_{i=1}^{n+1}(\frac{1}{dt}M\tau_{k}^{\prime})^{i})(M\partial_{t}(\textbf{U}^{n}-\textbf{U}^{n}_{h})+\mathcal{L}(\textbf{u}^{n,\theta};\textbf{U}^{n,\theta})-\mathcal{L}(\textbf{u}_{h}^{n,\theta};\textbf{U}_{h}^{n,\theta}))$
Now after applying error splitting for each of the terms and later using the
result obtained in (31) and properties of projection operators we have
rearranged the above equation (36) as follows:
$\forall\textbf{V}_{h}\in(V^{h})^{d}\times Q^{h}\times V^{h}$

$\rho(\frac{E^{A,n+1}_{\textbf{u}}-E^{A,n}_{\textbf{u}}}{dt},\textbf{v}_{h})+(\frac{E^{A,n+1}_{c}-E^{A,n}_{c}}{dt},d_{h})+\int_{\Omega}\mu(c^{n})\bigtriangledown
E^{A,n,\theta}_{\textbf{u}}:\bigtriangledown\textbf{v}_{h}+\\\
\int_{\Omega}\mu(c^{n})E^{A,n,\theta}_{\textbf{u}}\cdot\textbf{v}_{h}+\int_{\Omega}\tilde{\bigtriangledown}E^{A,n,\theta}_{c}\cdot\bigtriangledown\textbf{v}_{h}+\alpha\int_{\Omega}E^{A,n,\theta}_{c}d_{h}=\int_{\Omega}(\bigtriangledown\cdot\textbf{v}_{h})(E^{I,n,\theta}_{p}+E^{A,n,\theta}_{p})\\\
-\int_{\Omega}(\bigtriangledown\cdot
E^{A,n,\theta}_{\textbf{u}})q_{h}-\int_{\Omega}\mu(c^{n})\bigtriangledown
E^{I,n,\theta}_{\textbf{u}}:\bigtriangledown\textbf{v}_{h}-\sigma\int_{\Omega}E^{I,n,\theta}_{\textbf{u}}\cdot\textbf{v}_{h}-\int_{\Omega}\tilde{\bigtriangledown}E^{I,n,\theta}_{c}\cdot\bigtriangledown
d_{h}-\\\
\alpha\int_{\Omega}E^{I,n,\theta}_{c}d_{h}-\int_{\Omega}d_{h}\textbf{u}^{n}\cdot\bigtriangledown
E^{I,n,\theta}_{c}-\int_{\Omega}d_{h}\textbf{u}^{n}\cdot\bigtriangledown
E^{A,n,\theta}_{c}-\int_{\Omega}d_{h}E^{I,n}_{\textbf{u}}\cdot\bigtriangledown
c_{h}^{n,\theta}-\\\
\int_{\Omega}d_{h}E^{A,n}_{\textbf{u}}\cdot\bigtriangledown
c_{h}^{n,\theta}-c(\textbf{u}^{n},E^{I,n,\theta}_{\textbf{u}},\textbf{v}_{h})-c(\textbf{u}^{n},E^{A,n,\theta}_{\textbf{u}},\textbf{v}_{h})-c(E^{I,n}_{\textbf{u}},\textbf{u}^{n,\theta}_{h},\textbf{v}_{h})\\\
-c(E^{A,n}_{\textbf{u}},\textbf{u}^{n,\theta}_{h},\textbf{v}_{h})+\int_{\Omega}\mu(c^{n})E^{A,n,\theta}_{\textbf{u}}\cdot\textbf{v}_{h}-I_{1}-I_{2}-I_{3}-I_{4}-(\textbf{TE}^{n,\theta},\textbf{V}_{h})\hskip
14.22636pt$ (26)

Applying various properties of the projection operators we have the final
expression of $I_{1}$ above.

$\begin{split}I_{1}&=\sum_{k=1}^{n_{el}}(\tau_{k}^{\prime}(M\partial_{t}(\textbf{U}^{n}-\textbf{U}^{n}_{h})+\mathcal{L}(\textbf{u}^{n};\textbf{U}^{n,\theta})-\mathcal{L}(\textbf{u}^{n}_{h};\textbf{U}^{n,\theta}_{h})),-\mathcal{L}^{*}(\textbf{u}_{h};\textbf{V}_{h}))_{\Omega_{k}}\\\
&=\sum_{k=1}^{n_{el}}[(\tau_{1k}^{\prime}I_{d\times
d}\\{\rho\partial_{t}E^{I,n}_{\textbf{u}}+\rho
E^{I,n}_{\textbf{u}}\cdot\bigtriangledown\textbf{u}^{n,\theta}+\rho\textbf{u}_{h}^{n}\cdot\bigtriangledown
E^{I,n,\theta}_{\textbf{u}}-\mu(c^{n})\Delta E^{I,n,\theta}_{\textbf{u}}\\\
&\quad+\bigtriangledown
E^{I,n,\theta}_{p}\\},\rho(\textbf{u}_{h}\cdot\bigtriangledown)\textbf{v}_{h}+\mu(c)\Delta\textbf{u}_{h}+\bigtriangledown
p_{h})_{\Omega_{k}}+(\tau_{1k}^{\prime}I_{d\times
d}\\{\rho\partial_{t}E^{A,n}_{\textbf{u}}+\rho\\\ &\quad
E^{A,n}_{\textbf{u}}\cdot\bigtriangledown\textbf{u}^{n,\theta}+\rho\textbf{u}_{h}^{n}\cdot\bigtriangledown
E^{A,n,\theta}_{\textbf{u}}-\mu(c^{n})\Delta
E^{A,n,\theta}_{\textbf{u}}+\bigtriangledown
E^{A,n,\theta}_{p}\\},\rho(\textbf{u}_{h}\cdot\bigtriangledown)\textbf{v}_{h}\\\
&\quad+\mu(c)\Delta\textbf{u}_{h}+\bigtriangledown
p_{h})_{\Omega_{k}}+(\tau_{2k}^{\prime}\bigtriangledown\cdot(E^{I,n,\theta}_{\textbf{u}}+E_{\textbf{u}}^{A,n,\theta}),\bigtriangledown\cdot\textbf{v}_{h})_{\Omega_{k}}+(\tau_{3k}^{\prime}\\{\partial_{t}E^{I,n}_{c}\\\
&\quad-\bigtriangledown\cdot\tilde{\bigtriangledown}E^{I,n,\theta}_{c}+(E_{\textbf{u}}^{I,n}\cdot\bigtriangledown)c^{n,\theta}+(\textbf{u}_{h}^{n}\cdot\bigtriangledown)E^{I,n,\theta}_{c}+\alpha
E_{c}^{I,n,\theta}\\},\bigtriangledown\cdot\tilde{\bigtriangledown}d_{h}\\\
&\quad+\textbf{u}_{h}\cdot\bigtriangledown d_{h}-\alpha
d_{h})_{\Omega_{k}}+(\tau_{3}^{\prime}\\{\partial_{t}E^{A,n}_{c}-\bigtriangledown\cdot\tilde{\bigtriangledown}E^{A,n,\theta}_{c}+(E_{\textbf{u}}^{A,n}\cdot\bigtriangledown)c^{n,\theta}\\\
&\quad+(\textbf{u}_{h}^{n}\cdot\bigtriangledown)E^{A,n,\theta}_{c}+\alpha
E_{c}^{A,n,\theta}\\},\bigtriangledown\cdot\tilde{\bigtriangledown}d_{h}+\textbf{u}_{h}\cdot\bigtriangledown
d_{h}-\alpha d_{h})_{\Omega_{k}}]\\\
&=I_{1}^{1}+I_{1}^{2}+I_{1}^{3}+I_{1}^{4}+I_{1}^{5}\hskip
5.69054pt(say)\end{split}$ (27)

where $I_{1}^{i}$ for $i=1,2,...,5$ are five terms of $I_{1}$ which we will
discuss in the later part of the proof and since
$(1-\tau_{2}^{-1}\tau_{2}^{\prime})=0$ the next term will take the following
form

$\begin{split}I_{2}&=\sum_{k=1}^{n_{el}}((I-\tau_{k}^{-1}\tau_{k}^{\prime})(M\partial_{t}(\textbf{U}^{n}-\textbf{U}^{n}_{h})+\mathcal{L}(\textbf{u}^{n};\textbf{U}^{n,\theta})-\mathcal{L}(\textbf{u}^{n}_{h};\textbf{U}^{n,\theta}_{h})),-\textbf{V}_{h})_{\Omega_{k}}\\\
\end{split}$ (28)
$\begin{split}&=\sum_{k=1}^{n_{el}}[(\frac{\rho\tau_{1k}}{dt+\rho\tau_{1k}}I_{d\times
d}\\{\rho\partial_{t}E^{I,n}_{\textbf{u}}+\rho
E^{I,n}_{\textbf{u}}\cdot\bigtriangledown\textbf{u}^{n,\theta}+\rho\textbf{u}_{h}^{n}\cdot\bigtriangledown
E^{I,n,\theta}_{\textbf{u}}-\mu(c^{n})\\\ &\quad\Delta
E^{I,n,\theta}_{\textbf{u}}+\bigtriangledown
E^{I,n,\theta}_{p}\\},-\textbf{v}_{h})_{\Omega_{k}}+(\frac{\rho\tau_{1k}}{dt+\rho\tau_{1k}}I_{d\times
d}\\{\rho\partial_{t}E^{A,n}_{\textbf{u}}+\rho
E^{A,n}_{\textbf{u}}\cdot\bigtriangledown\textbf{u}^{n,\theta}\\\
&\quad+\rho\textbf{u}_{h}^{n}\cdot\bigtriangledown
E^{A,n,\theta}_{\textbf{u}}-\mu(c^{n})\Delta
E^{A,n,\theta}_{\textbf{u}}+\bigtriangledown
E^{A,n,\theta}_{p}\\},-\textbf{v}_{h})_{\Omega_{k}}+(\frac{\tau_{3k}}{dt+\tau_{3k}}\\{\partial_{t}E^{I,n}_{c}\\\
&\quad-\bigtriangledown\cdot\tilde{\bigtriangledown}E^{I,n,\theta}_{c}+(E_{\textbf{u}}^{I,n}\cdot\bigtriangledown)c^{n,\theta}+(\textbf{u}_{h}^{n}\cdot\bigtriangledown)E^{I,n,\theta}_{c}+\alpha
E^{I,n,\theta}_{c}\\},-d_{h})_{\Omega_{k}}+\\\
&\quad(\frac{\tau_{3k}}{dt+\tau_{3k}}\\{\partial_{t}E^{A,n}_{c}-\bigtriangledown\cdot\tilde{\bigtriangledown}E^{A,n,\theta}_{c}+(E_{\textbf{u}}^{A,n}\cdot\bigtriangledown)c^{n,\theta}+(\textbf{u}_{h}^{n}\cdot\bigtriangledown)E^{A,n,\theta}_{c}+\\\
&\quad\alpha E^{A,n,\theta}_{c}\\},-d_{h})_{\Omega_{k}}]\\\
&=I_{2}^{1}+I_{2}^{2}+I_{2}^{3}+I_{2}^{4}\hskip 5.69054pt(say)\end{split}$
(29)

The next terms of $LHS$ in (36) are as follows:

$\begin{split}I_{3}&=\sum_{k=1}^{n_{el}}(\tau_{k}^{-1}\tau_{k}^{\prime}\textbf{d},\textbf{V}_{h})_{\Omega_{k}}\\\
&=\sum_{k=1}^{n_{el}}(\\{\frac{dt}{dt+\rho\tau_{1k}}I_{d\times
d}\textbf{d},\textbf{v}_{h})_{\Omega_{k}}+(\frac{dt}{dt+\tau_{3k}}d_{4},d_{h})_{\Omega_{k}}\\}\hskip
85.35826pt\\\ \end{split}$ (30)

and

$\begin{split}I_{4}&=\sum_{k=1}^{n_{el}}(\tau_{k}^{\prime}\textbf{d},-\mathcal{L}^{*}(\textbf{u}_{h};\textbf{V}_{h}))_{\Omega_{k}}\\\
&=\sum_{k=1}^{n_{el}}\\{(\tau_{1k}^{\prime}I_{d\times
d}\textbf{d},\rho(\textbf{u}_{h}\cdot\bigtriangledown)\textbf{v}_{h}+\mu(c)\Delta\textbf{v}_{h}+\bigtriangledown
q_{h})_{\Omega_{k}}+(\tau_{3k}^{\prime}d_{3},\bigtriangledown\cdot\tilde{\bigtriangledown}d_{h}+\\\
&\quad\textbf{u}_{h}\cdot\bigtriangledown d_{h}-\alpha
d_{h})_{\Omega_{k}}\\}\\\ \end{split}$ (31)

Now we will treat each term separately to find out the estimate. Our procedure
contains finding two bounds: one is lower bound of $LHS$ and the other one is
upper bound for the terms in $RHS$ and combining those bounds in the equation
(36) we will finally obtain the required estimate. Before further proceeding
let us mention an important consideration: since the above equation holds for
all $\textbf{V}_{h}\in V_{s}^{h}\times V_{s}^{h}\times Q_{s}^{h}\times
V_{s}^{h}$, therefore in each term we replace $v_{1h},v_{2h},q_{h},d_{h}$ by
$E^{A,n,\theta}_{u1},E^{A,n,\theta}_{u2},E^{A,n,\theta}_{p},E^{A,n,\theta}_{c}$
respectively as these auxiliary part of the errors belonging to their
respective finite element spaces. From now onwards we will start derivation of
each expression after considering the replacements directly.
Applying the result obtained in (32) on the first term of $LHS$ and taking out
the $infimum$ of the coefficients of the remaining terms we can easily see
that

$\frac{\rho}{2dt}(\|E^{A,n+1}_{\textbf{u}}\|^{2}-\|E^{A,n}_{\textbf{u}}\|^{2})+\frac{1}{2dt}(\|E^{A,n+1}_{c}\|^{2}-\|E^{A,n}_{c}\|^{2})+\mu_{l}\|E^{A,n,\theta}_{\textbf{u}}\|_{1}^{2}\\\
+D_{l}\mid E^{A,n,\theta}_{c}\mid_{1}^{2}+\alpha\|E^{A,n,\theta}_{c}\|^{2}\leq
LHS=RHS\\\ $ (32)

where $D_{l}$= min
$\\{\underset{\Omega}{inf}D_{1},\underset{\Omega}{inf}D_{2}\\}$.
Now we will find upper bounds of each of the terms in the $RHS$ of the
equation (37). We usually use $Cauchy-Schwarz$ and $Young^{\prime}s$
inequality to reach at the desired bounds. We have already estimated the
bounds of few terms on $RHS$ in [14]. Therefore we only mention the results
here and the estimation of the remaining terms are shown later in details.

$\begin{split}\int_{\Omega}(\bigtriangledown\cdot
E^{A,n,\theta}_{\textbf{u}})E^{I,n,\theta}_{p}&\leq\epsilon_{1}C^{2}h^{2}(\frac{1+\theta}{2}\|p^{n+1}\|_{1}+\frac{1-\theta}{2}\|p^{n}\|_{1})^{2}+\\\
&\quad\frac{1}{2\epsilon_{1}}\mid E^{A,n,\theta}_{\textbf{u}}\mid_{1}^{2}\\\
-\int_{\Omega}\mu(c^{n})\bigtriangledown
E^{I,n,\theta}_{\textbf{u}}:\bigtriangledown
E^{A,n,\theta}_{\textbf{u}}&\leq\epsilon_{2}\mu_{u}C^{2}h^{2}(\frac{1+\theta}{2}\|\textbf{u}^{n+1}\|_{2}+\frac{1-\theta}{2}\|\textbf{u}^{n}\|_{2})^{2}+\\\
&\quad\frac{\mu_{u}}{2\epsilon_{2}}\mid
E^{A,n,\theta}_{\textbf{u}}\mid_{1}^{2}\\\
-\int_{\Omega}\tilde{\bigtriangledown}E^{I,n,\theta}_{c}\cdot\bigtriangledown
E^{A,n,\theta}_{c}&\leq\frac{D_{m}\epsilon_{3}}{2}C^{2}h^{2}(\frac{1+\theta}{2}\|c^{n+1}\|_{2}+\frac{1-\theta}{2}\|c^{n}\|_{2})^{2}+\\\
&\quad\frac{D_{m}}{2\epsilon_{3}}(\|\frac{\partial
E^{A,n,\theta}_{c}}{\partial x}\|^{2}+\|\frac{\partial
E^{A,n,\theta}_{c}}{\partial y}\|^{2})\\\
-\sigma\int_{\Omega}E^{I,n,\theta}_{\textbf{u}}\cdot
E^{A,n,\theta}_{\textbf{u}}&\leq\frac{\epsilon_{4}}{2}\sigma
h^{4}(\frac{1+\theta}{2}\|\textbf{u}^{n+1}\|_{2}+\frac{1-\theta}{2}\|\textbf{u}^{n}\|_{2})^{2}+\\\
&\quad\frac{\sigma}{2\epsilon_{4}}\|E^{A,n,\theta}_{\textbf{u}}\|^{2}\\\
-\alpha\int_{\Omega}E^{I,n,\theta}_{c}E^{A,n,\theta}_{c}&\leq\frac{\epsilon_{5}}{2}\alpha
h^{4}(\frac{1+\theta}{2}\|c^{n+1}\|_{2}+\frac{1-\theta}{2}\|c^{n}\|_{2})^{2}+\\\
&\quad\frac{\alpha}{2\epsilon_{5}}\|E^{A,n,\theta}_{c}\|^{2}\end{split}$ (33)

where $D_{m}$= max
$\\{\underset{\Omega}{sup}D_{1},\underset{\Omega}{sup}D_{2}\\}$
Now applying $Poincare$ inequality in the following we have:

$\begin{split}\int_{\Omega}\mu(c^{n})\\{(E^{A,n,\theta}_{u1})^{2}+(E^{A,n,\theta}_{u2})^{2}\\}&\leq\mu_{u}(\|E^{A,n,\theta}_{u1}\|^{2}+\|E^{A,n,\theta}_{u2}\|^{2})\\\
&\leq\mu_{u}C_{P}(\mid E^{A,n,\theta}_{u1}\mid_{1}^{2}+\mid
E^{A,n,\theta}_{u2}\mid_{1}^{2})\end{split}$ (34)

where $C_{P}$ is the $Poincare$ constant. The next term is estimated following
[14]

$\begin{split}-\int_{\Omega}E^{A,n,\theta}_{c}\textbf{u}^{n}\cdot\bigtriangledown
E^{A,n,\theta}_{c}&\leq 2\bar{C}_{1}^{2}(\|\frac{\partial
E^{A,n,\theta}_{c}}{\partial x}\|^{2}+\|\frac{\partial
E^{A,n,\theta}_{c}}{\partial y}\|^{2})\\\
-\int_{\Omega}E^{A,n,\theta}_{c}\textbf{u}^{n}\cdot\bigtriangledown
E^{I,n,\theta}_{c}&\leq\bar{C}_{1}^{2}C^{2}h^{2}\epsilon_{6}(\frac{1+\theta}{2}\|c^{n+1}\|_{2}+\frac{1-\theta}{2}\|c^{n}\|_{2})^{2}+\\\
&\quad\frac{\bar{C}_{1}^{2}}{\epsilon_{6}}\|E^{A,n,\theta}_{c}\|^{2}\end{split}$
(35)

Now we estimate the trilinear term $c(\cdot,\cdot,\cdot)$ using it’s
properties (a) and (b) given in section 2. Let us start the estimation with
the first trilinear term on $RHS$ of (37) as follows:

$\begin{split}-c(\textbf{u}^{n},E^{I,n,\theta}_{\textbf{u}},E^{A,n,\theta}_{\textbf{u}})&\leq
C\|\textbf{u}^{n}\|_{2}\|E^{I,n,\theta}_{\textbf{u}}\|_{1}\|E^{A,n,\theta}_{\textbf{u}}\|\\\
&\leq\frac{C_{2}\epsilon_{7}}{2}\|E^{I,n,\theta}_{\textbf{u}}\|_{1}^{2}+\frac{C_{2}}{2\epsilon_{7}}\|E^{A,n,\theta}_{\textbf{u}}\|^{2}\\\
&\leq\frac{C_{2}\epsilon_{7}}{2}(\frac{1+\theta}{2}\|E^{I,n+1}_{\textbf{u}}\|_{1}+\frac{1-\theta}{2}\|E^{I,n}_{\textbf{u}}\|_{1})^{2}+\frac{C_{2}}{2\epsilon_{7}}\|E^{A,n,\theta}_{\textbf{u}}\|^{2}\\\
&\leq\frac{\epsilon_{7}}{2}C_{2}h^{2}(\frac{1+\theta}{2}\|\textbf{u}^{n+1}\|_{2}+\frac{1-\theta}{2}\|\textbf{u}^{n}\|_{2})^{2}+\frac{C_{2}}{2\epsilon_{7}}\|E^{A,n,\theta}_{\textbf{u}}\|^{2}\end{split}$
(36)

By the property (a) of trilinear case for both the linear and non-linear
cases:

$-c(\textbf{u}^{n},E^{A,n,\theta}_{\textbf{u}},E^{A,n,\theta}_{\textbf{u}})=0$
(37)

The next term,

$\begin{split}-c(E^{I,n}_{\textbf{u}},\textbf{u}^{n,\theta}_{h},E^{A,n,\theta}_{\textbf{u}})&=c(E^{I,n}_{\textbf{u}},E^{I,n,\theta}_{\textbf{u}},E^{A,n,\theta}_{\textbf{u}})+c(E^{I,n}_{\textbf{u}},E^{A,n,\theta}_{\textbf{u}},E^{A,n,\theta}_{\textbf{u}})-\\\
&\quad
c(E^{I,n}_{\textbf{u}},\textbf{u}^{n,\theta},E^{A,n,\theta}_{\textbf{u}})\\\
&\leq
C\|E^{I,n}_{\textbf{u}}\|\|E^{I,n,\theta}_{\textbf{u}}\|_{1}\|E^{A,n,\theta}_{\textbf{u}}\|_{1}+C\|E^{I,n,\theta}_{\textbf{u}}\|\|\textbf{u}^{n,\theta}\|_{2}\|E^{A,n,\theta}_{\textbf{u}}\|_{1}\\\
&\leq
h^{4}\frac{C_{2}\epsilon_{8}}{2}\|\textbf{u}^{n}\|(\frac{1+\theta}{2}\|\textbf{u}^{n+1}\|_{2}+\frac{1-\theta}{2}\|\textbf{u}^{n}\|_{2})^{2}+\frac{C_{2}}{\epsilon_{8}}\|E^{A,n,\theta}_{\textbf{u}}\|_{1}^{2}\\\
&\quad+h^{4}\frac{C_{2}\epsilon_{8}}{2}(\frac{1+\theta}{2}\|\textbf{u}^{n+1}\|_{2}+\frac{1-\theta}{2}\|\textbf{u}^{n}\|_{2})^{2}\end{split}$
(38)

and

$\begin{split}-c(E^{A,n}_{\textbf{u}},\textbf{u}^{n,\theta}_{h},E^{A,n,\theta}_{\textbf{u}})&=c(E^{A,n}_{\textbf{u}},E^{I,n,\theta}_{\textbf{u}},E^{A,n,\theta}_{\textbf{u}})+c(E^{A,n}_{\textbf{u}},E^{A,n,\theta}_{\textbf{u}},E^{A,n,\theta}_{\textbf{u}})\\\
&\quad-c(E^{A,n}_{\textbf{u}},\textbf{u}^{n,\theta},E^{A,n,\theta}_{\textbf{u}})\\\
&\leq
C\|E^{A,n,\theta}_{\textbf{u}}\|\|E^{A,n,\theta}_{\textbf{u}}\|_{1}\\{\|E^{I,n,\theta}_{\textbf{u}}\|_{2}+\|\textbf{u}^{n,\theta}\|_{2}\\}\\\
&\leq
C\|E^{A,n,\theta}_{\textbf{u}}\|_{1}^{2}(\frac{1+\theta}{2}\|E^{I,n+1}_{\textbf{u}}\|_{2}+\frac{1-\theta}{2}\|E^{I,n}_{\textbf{u}}\|_{2})+\\\
&\quad C_{2}\|E^{A,n,\theta}_{\textbf{u}}\|_{1}^{2}\\\ &\leq
C\|E^{A,n,\theta}_{\textbf{u}}\|_{1}^{2}\\{(\frac{1+\theta}{2}\|\textbf{u}^{n+1}\|_{2}+\frac{1-\theta}{2}\|\textbf{u}^{n}\|_{2})+C_{2}\\}\\\
&\leq C_{2}^{\prime}\|E^{A,n,\theta}_{\textbf{u}}\|_{1}^{2}\end{split}$ (39)

Now we will find bounds for each remaining term of $I_{1}$ separately. Before
going to further calculations let us mention an important observation: By the
virtue of the choices of the finite element spaces $V^{h}$ and $Q^{h}$, we can
clearly say that over each element sub-domain $\Omega_{k}$ every function
belonging to that spaces and their first and second order derivatives all are
bounded functions. We can always find positive finite real numbers to bound
each of the functions over element sub-domain. We will use this fact for
several times further.
Let us start with $I_{1}^{1}$

$\begin{split}I_{1}^{1}&=\sum_{k=1}^{n_{el}}(\rho\tau_{1k}^{\prime}I_{d\times
d}\partial_{t}E^{I,n}_{\textbf{u}},\rho(\textbf{u}_{h}\cdot\bigtriangledown)E^{A,n,\theta}_{\textbf{u}}+\mu(c)\Delta
E^{A,n,\theta}_{\textbf{u}}+\bigtriangledown
E^{A,n,\theta}_{p})_{\Omega_{k}}+\\\
&\quad\sum_{k=1}^{n_{el}}(\tau_{1k}^{\prime}I_{d\times d}\\{\rho
E^{I,n,\epsilon\theta}_{\textbf{u}}\cdot\bigtriangledown\textbf{u}^{n,\theta}+\rho\textbf{u}_{h}^{n,\theta}\cdot\bigtriangledown
E^{I,n,\theta}_{\textbf{u}}-\mu(c^{n})\Delta
E^{I,n,\theta}_{\textbf{u}}+\bigtriangledown E^{I,n,\theta}_{p}\\},\\\
&\quad\rho(\textbf{u}_{h}\cdot\bigtriangledown)E^{A,n,\theta}_{\textbf{u}}+\mu(c)\Delta
E^{A,n,\theta}_{\textbf{u}}+\bigtriangledown
E^{A,n,\theta}_{p})_{\Omega_{k}}\\\ &=I_{1}^{11}+I_{1}^{12}\hskip
5.69054pt(say)\end{split}$ (40)

Now we present the estimations of these two terms separately in details.
According to the above observation we can find bounds on each of the terms
$\textbf{u}_{h},E^{A,n,\theta}_{\textbf{u}},E^{A,n,\theta}_{p}$ and their
first and second order derivatives over each sub-domain $\Omega_{k}$. Applying
these bounds in the following we will have

$\begin{split}I_{1}^{11}&=\int_{\Omega^{\prime}}\rho\tau_{1k}^{\prime}I_{d\times
d}\frac{E^{I,n+1}_{\textbf{u}}-E^{I,n}_{\textbf{u}}}{dt}\\{\rho(\textbf{u}_{h}\cdot\bigtriangledown)E^{A,n,\theta}_{\textbf{u}}+\mu(c)\Delta
E^{A,n,\theta}_{\textbf{u}}+\bigtriangledown E^{A,n,\theta}_{p}\\}\\\
&\leq\rho\frac{\mid\tau_{1}^{\prime}\mid}{dt}\\{\sum_{k=1}^{n_{el}}D_{B_{1k}}\\}(\|E^{I,n+1}_{\textbf{u}}\|+\|E^{I,n}_{\textbf{u}}\|)\\\
&\leq\rho Ch^{2}\frac{\mid\tau_{1}\mid}{\mid
dt+\rho\tau_{1}\mid}\\{\sum_{k=1}^{n_{el}}D_{B_{1k}}\\}(\|\textbf{u}^{n+1}\|_{2}+\|\textbf{u}^{n}\|_{2})\\\
&\leq h^{2}\frac{\rho CC_{\tau_{1}}}{(T_{0}-\rho
C_{\tau_{1}})}\\{\sum_{k=1}^{n_{el}}D_{B_{1k}}\\}(\|\textbf{u}^{n+1}\|_{2}+\|\textbf{u}^{n}\|_{2})\\\
\end{split}$ (41)

where the constant $D_{B_{1k}}$ is obtained after imposing bounds on the above
bracketed terms over each sub-domain $\Omega_{k}$. $C_{\tau_{1}}$and $T$ are
upper bounds on respectively $\tau_{1}$ and $dt$. Since $dt$ is a non-zero
positive real number, let $T_{0}$ is lower bound on $dt$. In order to make
$(T_{0}-\rho C_{\tau_{1}})$ positive we have to take $h$ very small.
Now the estimation of the second term is as follows:

$\begin{split}I_{1}^{12}&=\sum_{k=1}^{n_{el}}(\tau_{1k}^{\prime}I_{d\times
d}\\{\rho
E^{I,n}_{\textbf{u}}\cdot\bigtriangledown\textbf{u}^{n,\theta}+\rho\textbf{u}^{n}\cdot\bigtriangledown
E^{I,n,\theta}_{\textbf{u}}-\rho E^{I,n}_{\textbf{u}}\cdot\bigtriangledown
E^{I,n,\theta}_{\textbf{u}}-\\\ &\quad\rho
E^{A,n}_{\textbf{u}}\cdot\bigtriangledown
E^{I,n,\theta}_{\textbf{u}}-\mu(c^{n})\Delta
E^{I,n,\theta}_{\textbf{u}}+\bigtriangledown
E^{I,n,\theta}_{p}\\},\rho(\textbf{u}_{h}\cdot\bigtriangledown)E^{A,n,\theta}_{\textbf{u}}+\\\
&\quad\mu(c)\Delta E^{A,n,\theta}_{\textbf{u}}+\bigtriangledown
E^{A,n,\theta}_{p})_{\Omega_{k}}\\\
&\leq\mid\tau_{1k}^{\prime}\mid[\sum_{k=1}^{n_{el}}D_{B_{1k}}\\{\sum_{i=1}^{d}(\rho\|E^{I,n}_{ui}\|\|\frac{\partial
u_{i}^{n,\theta}}{\partial x_{i}}\|+\rho\|u_{i}^{n}\|\|\frac{\partial
E^{I,n,\theta}_{ui}}{\partial x_{i}}\|+\|E^{I,n}_{ui}\|\\\
&\quad\rho\|\frac{\partial E^{I,n,\theta}_{ui}}{\partial
x_{i}}\|+\rho\|E^{A,n}_{ui}\|_{k}\|\frac{\partial
E^{I,n,\theta}_{ui}}{\partial x_{i}}\|+\mu(c)\|\Delta
E^{I,n,\theta}_{ui}\|+\|\frac{\partial E^{I,n,\theta}_{p}}{\partial
x_{i}}\|)\\}]\\\ &\leq\frac{\mid\tau_{1}\mid T}{(T_{0}-\rho
C_{\tau_{1}})}C(\sum_{k=1}^{n_{el}}D_{B_{1k}})[\sum_{i=1}^{d}\\{\rho
h^{2}\|\frac{\partial u_{i}^{n,\theta}}{\partial x_{i}}\|+\rho
h\|u_{i}^{n,\epsilon\theta}\|+\rho hB^{i}_{1k}+\mu_{u}\\\
&\quad+h^{3}(\frac{1+\theta}{2}\|u_{i}^{n+1}\|_{2}+\frac{1-\theta}{2}\|u_{i}^{n}\|_{2})\\}(\frac{1+\theta}{2}\|u_{i}^{n+1}\|_{2}+\frac{1-\theta}{2}\|u_{i}^{n}\|_{2})+\\\
&\quad
h(\frac{1+\theta}{2}\|p^{n+1}\|_{1}+\frac{1-\theta}{2}\|p^{n}\|_{1})]\end{split}$
(42)

where $B^{i}_{1k}$ are bounds on $E^{A,n,\theta}_{ui}$ for $i=1,...,d$. Now
the next term $I_{1}^{2}$ can be estimated by dividing it into two terms
$I_{1}^{21}$ and $I_{1}^{22}$ as above. Therefore we directly start the
estimation here with the first term of $I_{1}^{2}$ denoting that by
$I_{1}^{21}$.

$\begin{split}I_{1}^{21}&=\int_{\Omega^{\prime}}\rho\tau_{1k}^{\prime}I_{d\times
d}\frac{E^{A,n+1}_{\textbf{u}}-E^{A,n}_{\textbf{u}}}{dt}\\{\rho(\textbf{u}_{h}\cdot\bigtriangledown)E^{A,n,\theta}_{\textbf{u}}+\mu(c)\Delta
E^{A,n,\theta}_{\textbf{u}}+\bigtriangledown E^{A,n,\theta}_{p}\\}\\\
&\leq\frac{\rho TC_{\tau_{1}}}{dt(T_{0}-\rho
C_{\tau_{1}})}(\sum_{k=1}^{n_{el}}D_{B_{1k}})\\{\|E^{A,n+1}_{\textbf{u}}\|^{2}-\|E^{A,n}_{\textbf{u}}\|^{2}\\}\end{split}$
(43)

and

$\begin{split}I_{1}^{22}&=\sum_{k=1}^{n_{el}}(\tau_{1k}^{\prime}I_{d\times
d}\\{\rho
E^{A,n}_{\textbf{u}}\cdot\bigtriangledown\textbf{u}^{n,\theta}+\rho\textbf{u}^{n}\cdot\bigtriangledown
E^{A,n}_{\textbf{u}}-\rho E^{I,n}_{\textbf{u}}\cdot\bigtriangledown
E^{A,n,\theta}_{\textbf{u}}-\rho E^{A,n}_{\textbf{u}}\cdot\bigtriangledown\\\
&\quad E^{A,n,\theta}_{\textbf{u}}-\mu(c^{n})\Delta
E^{A,n,\theta}_{\textbf{u}}+\bigtriangledown
E^{A,n,\theta}_{p}\\},\rho(\textbf{u}_{h}\cdot\bigtriangledown)E^{A,n,\theta}_{\textbf{u}}+\mu(c)\Delta
E^{A,n,\theta}_{\textbf{u}}+\\\ &\quad\bigtriangledown
E^{A,n,\theta}_{p})_{\Omega_{k}}\\\
&\leq\mid\tau_{1}^{\prime}\mid[\sum_{k=1}^{n_{el}}D_{B_{1k}}\\{\sum_{i=1}^{d}(\rho\|E^{A,n}_{ui}\|_{k}\|\frac{\partial
u_{i}^{n}}{\partial x_{i}}\|+\rho\|u_{i}^{n,\theta}\|\|\frac{\partial
E^{A,n,\theta}_{ui}}{\partial x_{i}}\|_{k}+\rho\|E^{I,n}_{ui}\|\\\
&\quad\|\frac{\partial E^{A,n,\theta}_{ui}}{\partial
x_{i}}\|_{k}+\|E^{A,n}_{ui}\|_{k}\|\frac{\partial
E^{A,n,\theta}_{ui}}{\partial x_{i}}\|_{k}+\mu_{u}\|\Delta
E^{A,n,\theta}_{ui}\|_{k}+\|\frac{\partial E^{A,n,\theta}_{p}}{\partial
x_{i}}\|_{k})\\}]\end{split}$ (44)

Applying bounds on the functions belonging to $V^{h}$ and $Q^{h}$ on the above
equation and denoting that bound by $\bar{D}_{B_{1k}}$ we have

$I_{1}^{22}\leq\frac{\mid\tau_{1}\mid T}{(T_{0}-\rho
C_{\tau_{1}})}(\sum_{k=1}^{n_{el}}D_{B_{1k}}\bar{D}_{B_{1k}})$ (45)

Now expanding out the next term of $I_{1}$ we can proceed to estimate that in
the following way:

$\begin{split}I_{1}^{3}&=\sum_{k=1}^{n_{el}}\int_{\Omega_{k}}\tau_{2}\\{(\bigtriangledown\cdot
E^{I,n,\theta}_{\textbf{u}})\cdot(\bigtriangledown\cdot
E^{A,n,\theta}_{\textbf{u}})+(\bigtriangledown\cdot
E^{A,n,\theta}_{\textbf{u}})^{2}\\}\\\ &\leq
C_{\tau_{2}}\\{\sum_{i=1}^{d}(\|\frac{\partial E^{I,n,\theta}_{ui}}{\partial
x_{i}}\|^{2}+C_{1}\|\frac{\partial E^{A,n,\theta}_{ui}}{\partial
x_{i}}\|^{2})\\}\\\ &\leq
C_{\tau_{2}}\\{h^{2}\sum_{i=1}^{d}(\frac{1+\theta}{2}\|u_{i}^{n+1}\|_{2}+\frac{1-\theta}{2}\|u_{i}^{n}\|_{2})^{2}+\mid
E^{A,n,\theta}_{\textbf{u}}\mid_{1}^{2}\\}\end{split}$ (46)

where $C_{\tau_{2}}$ is the maximum numerical value of $\tau_{2}$ over
$\Omega$. Now the remaining terms of $I_{1}$ associated with the variable $c$
representing concentration are estimated as follows:

$\begin{split}I_{1}^{4}&=\sum_{k=1}^{n_{el}}\tau_{3k}^{\prime}(\partial_{t}E^{I,n}_{c},\bigtriangledown\cdot\tilde{\bigtriangledown}E^{A,n,\theta}_{c}+\textbf{u}_{h}\cdot\bigtriangledown
E^{A,n,\theta}_{c}-\alpha E^{A,n,\theta}_{c})_{\Omega_{k}}+\\\
&\quad\sum_{k=1}^{n_{el}}\tau_{3k}^{\prime}(-\bigtriangledown\cdot\tilde{\bigtriangledown}E^{I,n,\theta}_{c}+(E^{I,n}_{\textbf{u}}\cdot\bigtriangledown)c^{n,\theta}+(\textbf{u}_{h}^{n}\cdot\bigtriangledown)E^{I,n,\theta}_{c}+\alpha
E^{I,n,\theta}_{c},\\\
&\quad\bigtriangledown\cdot\tilde{\bigtriangledown}E^{A,n,\theta}_{c}+\textbf{u}_{h}\cdot\bigtriangledown
E^{A,n,\theta}_{c}-\alpha E^{A,n,\theta}_{c})_{\Omega_{k}}\\\ &\leq
Ch^{2}\mid\tau_{3}^{\prime}\mid\\{\sum_{k=1}^{n_{el}}D_{B_{3k}}\\}(\|c^{n+1}\|_{2}+\|c^{n}\|_{2})+\mid\tau_{3}^{\prime}\mid[\sum_{k=1}^{n_{el}}D_{B_{3k}}\\{\sum_{i=1}^{d}(\|E^{I,n}_{ui}\|\\\
&\quad\|\frac{\partial c^{n,\theta}}{\partial
x_{i}}\|+D_{im}\|\frac{\partial^{2}E^{I,n,\theta}_{c}}{\partial
x_{i}^{2}}\|+\bar{D}_{im}\|\frac{\partial E^{I,n,\theta}_{c}}{\partial
x_{i}}\|+(\|u_{i}^{n}\|+\|E^{I,n}_{ui}\|+\|E^{A,n}_{ui}\|)\\\
&\quad\|\frac{\partial E^{I,n,\theta}_{c}}{\partial
x_{i}}\|+\mid\alpha\mid\|E^{I,n,\theta}_{c}\|)\\}]\\\ &\leq
h^{2}\frac{CC_{\tau_{3}}}{(T_{0}-C_{\tau_{3}})}\\{\sum_{k=1}^{n_{el}}D_{B_{3k}}\\}(\|c^{n+1}\|_{2}+\|c^{n}\|_{2})+\frac{C\mid\tau_{3}\mid}{(T_{0}-C_{\tau_{3}})}[\sum_{k=1}^{n_{el}}D_{B_{3k}}\sum_{i=1}^{d}\\{h^{2}\\\
&\quad(\frac{1+\theta}{2}\|u_{i}^{n+1}\|_{2}+\frac{1-\theta}{2}\|u_{i}^{n}\|_{2})\|\frac{\partial
c^{n,\theta}}{\partial x_{i}}\|+(D_{im}+h\bar{D}_{im}+\mid\alpha\mid
h^{2}+h\|u_{i}^{n,\theta}\|\\\
&\quad+hB_{1k}^{i}+h^{3}(\frac{1+\theta}{2}\|u_{i}^{n+1}\|_{2}+\frac{1-\theta}{2}\|u_{i}^{n}\|_{2}))(\frac{1+\theta}{2}\|c^{n+1}\|_{2}+\frac{1-\theta}{2}\|c^{n}\|_{2})\\}]\end{split}$
(47)

Let $D_{B_{3k}}$ be the summation of the bounds imposed on the elements
$\frac{\partial^{2}E^{A,n,\theta}_{c}}{\partial x_{i}^{2}}$, $\frac{\partial
E^{A,n,\theta}_{c}}{\partial x_{i}}$, $E^{A,n,\theta}_{c}$ belonging to finite
element space $V^{h}$ and $D_{im},\bar{D}_{im}$ be the $supremum$ of $D_{i}$
and $\frac{\partial D_{i}}{\partial x_{i}}$ respectively over each sub-domain
$\Omega_{k}$ for $i=1,...,d$. Now the estimation of the last term $I_{1}^{5}$
follows the same way as above and considering $\bar{D}_{B_{3k}}$ as an
expression to denote the estimated result briefly the derivation of the bound
of $I_{1}^{5}$ is in the following:

$\begin{split}I_{1}^{5}&=\sum_{k=1}^{n_{el}}\tau_{3k}^{\prime}(\partial_{t}E^{A,n}_{c},\bigtriangledown\cdot\tilde{\bigtriangledown}E^{A,n,\theta}_{c}+\textbf{u}_{h}\cdot\bigtriangledown
E^{A,n,\theta}_{c}-\alpha E^{A,n,\theta}_{c})_{\Omega_{k}}+\\\
&\quad\sum_{k=1}^{n_{el}}\tau_{3k}^{\prime}(-\bigtriangledown\cdot\tilde{\bigtriangledown}E^{A,n,\theta}_{c}+(E^{A,n}_{\textbf{u}}\cdot\bigtriangledown)c^{n,\theta}+(\textbf{u}_{h}^{n}\cdot\bigtriangledown)E^{A,n,\theta}_{c}+\alpha
E^{A,n,\theta}_{c},\\\
&\quad\bigtriangledown\cdot\tilde{\bigtriangledown}E^{A,n,\theta}_{c}+\textbf{u}_{h}\cdot\bigtriangledown
E^{A,n,\theta}_{c}-\alpha E^{A,n,\theta}_{c})_{\Omega_{k}}\\\
&\leq\frac{TC_{\tau_{3}}}{dt(T_{0}-C_{\tau_{3}})}\\{\sum_{k=1}^{n_{el}}D_{B_{3k}}\\}(\|E^{A,n+1}_{c}\|^{2}-\|E^{A,n}_{c}\|^{2})+\\\
&\quad\frac{C\mid\tau_{3}\mid}{(T_{0}-C_{\tau_{3}})}\\{\sum_{k=1}^{n_{el}}\bar{D}_{B_{3k}}D_{B_{3k}}\\}\end{split}$
(48)

This completes estimation of $I_{1}$ finally. Now we see that the terms
$I_{2}^{i}$ (for $i=1,...,4$) are same as that of $I_{1}$. Therefore we only
mention the results here for each of them as follows:

$\begin{split}I_{2}^{1}&=\sum_{k=1}^{n_{el}}(\frac{\rho\tau_{1k}}{dt+\rho\tau_{1k}}I_{d\times
d}\partial_{t}E^{I,n}_{\textbf{u}},-E^{A,n,\theta}_{\textbf{u}})_{\Omega_{k}}+\sum_{k=1}^{n_{el}}(\frac{\rho\tau_{1k}}{dt+\rho\tau_{1k}}I_{d\times
d}\\{\rho E^{I,n}_{\textbf{u}}\cdot\bigtriangledown\textbf{u}^{n,\theta}\\\
&\quad+\rho\textbf{u}^{n}\cdot\bigtriangledown
E^{I,n,\theta}_{\textbf{u}}-\rho E^{I,n}_{\textbf{u}}\cdot\bigtriangledown
E^{I,n,\theta}_{\textbf{u}}-\rho E^{A,n}_{\textbf{u}}\cdot\bigtriangledown
E^{I,n,\theta}_{\textbf{u}}-\mu(c^{n})\Delta E^{I,n,\theta}_{\textbf{u}}\\\
&\quad+\bigtriangledown
E^{I,n,\theta}_{p}\\},-E^{A,n,\theta}_{\textbf{u}})_{\Omega_{k}}\\\ &\leq
h^{2}\frac{\rho CC_{\tau_{1}}}{(T_{0}-\rho
C_{\tau_{1}})}[\sum_{i=1}^{d}\\{\sum_{k=1}^{n_{el}}B_{1k}^{i}\\}(\|u_{i}^{n+1}\|_{2}+\|u_{i}^{n}\|_{2})]+\frac{\mid\tau_{1}\mid}{(T_{0}-\rho
C_{\tau_{1}})}\sum_{i=1}^{d}(\sum_{k=1}^{n_{el}}B_{1k}^{i})\\\ &\quad[\\{\rho
h^{2}\|\frac{\partial u_{i}^{n,\theta}}{\partial x_{i}}\|+\rho
h\|u_{i}^{n}\|+\rho
hB^{i}_{1k}+\mu_{u}+h^{3}(\frac{1+\theta}{2}\|u_{i}^{n+1}\|_{2}+\frac{1-\theta}{2}\|u_{i}^{n}\|_{2})\\}\\\
&\quad(\frac{1+\theta}{2}\|u_{i}^{n+1}\|_{2}+\frac{1-\theta}{2}\|u_{i}^{n}\|_{2})+h(\frac{1+\theta}{2}\|p^{n+1}\|_{1}+\frac{1-\theta}{2}\|p^{n}\|_{1})]\end{split}$
(49)

and

$\begin{split}I_{2}^{2}&=\sum_{k=1}^{n_{el}}(\frac{\rho\tau_{1k}}{dt+\rho\tau_{1k}}I_{d\times
d}\partial_{t}E^{A,n}_{\textbf{u}},-E^{A,n,\theta}_{\textbf{u}})_{\Omega_{k}}+\sum_{k=1}^{n_{el}}(\frac{\rho\tau_{1k}}{dt+\rho\tau_{1k}}I_{d\times
d}\\{\rho E^{A,n}_{\textbf{u}}\cdot\bigtriangledown\textbf{u}^{n,\theta}\\\
&\quad+\rho\textbf{u}^{n}\cdot\bigtriangledown
E^{A,n,\theta}_{\textbf{u}}-\rho E^{I,n}_{\textbf{u}}\cdot\bigtriangledown
E^{A,n,\theta}_{\textbf{u}}-\rho E^{A,n}_{\textbf{u}}\cdot\bigtriangledown
E^{A,n,\theta}_{\textbf{u}}-\mu(c^{n})\Delta E^{A,n,\theta}_{\textbf{u}}\\\
&\quad+\bigtriangledown
E^{A,n,\theta}_{p}\\},-E^{A,n,\theta}_{\textbf{u}})_{\Omega_{k}}\\\
&\leq\frac{\rho TC_{\tau_{1}}}{dt(T_{0}-\rho
C_{\tau_{1}})}[\sum_{i=1}^{d}(\sum_{k=1}^{n_{el}}B^{i}_{1k})]\\{\|E^{A,n+1}_{\textbf{u}}\|^{2}-\|E^{A,n}_{\textbf{u}}\|^{2}\\}+\\\
&\quad\frac{\mid\tau_{1}\mid T}{(T_{0}-\rho
C_{\tau_{1}})}[\sum_{i=1}^{d}\sum_{k=1}^{n_{el}}\bar{D}_{B_{1k}}B^{i}_{1k}]\end{split}$
(50)

From these results it is clear that estimations of the remaining terms of
$I_{2}$ follow the same path as we have done for $I_{1}^{4}$ and $I_{1}^{5}$.
Hence considering $B_{2k}$ as bound for $E^{A,n,\theta}_{c}$ over each sub-
domain $\Omega_{k}$ we are skipping the repetition in mentioning the similar
kind of results, though they will be added up in the final stage of combining
all the results.
Now the job is to estimate next part denoted by $I_{3}$ and $I_{4}$ of the
equation (43) which contain the matrix d. Earlier we have mentioned that
$d_{2}$ is $zero$. Let us look at the other three terms explicitly.

$\begin{split}\textbf{d}_{1}&=\\{\sum_{i=1}^{n+1}(\frac{\rho}{dt}\tau_{1}^{\prime})^{i}\\}I_{d\times
d}[\rho\partial_{t}(\textbf{u}^{n}-\textbf{u}_{h}^{n})+\rho((\textbf{u}^{n}-\textbf{u}_{h}^{n})\cdot\bigtriangledown)\textbf{u}^{n,\theta}+\rho(\textbf{u}_{h}^{n}\cdot\bigtriangledown)\\\
&\quad(\textbf{u}^{n,\theta}-\textbf{u}_{h}^{n,\theta})-\mu(c^{n})\Delta(\textbf{u}^{n,\theta}-\textbf{u}_{h}^{n,\theta})+\bigtriangledown(p^{n,\theta}-p_{h}^{n,\theta})]\\\
&\leq\\{\sum_{i=1}^{\infty}(\frac{\rho}{dt}\tau_{1}^{\prime})^{i}\\}I_{d\times
d}[\rho\partial_{t}(E^{I,n}_{\textbf{u}}+E^{A,n}_{\textbf{u}})+\rho((E^{I,n}_{\textbf{u}}+E^{A,n}_{\textbf{u}})\cdot\bigtriangledown)\textbf{u}^{n,\theta}+\\\
&\quad\rho(\textbf{u}_{h}^{n}\cdot\bigtriangledown)(E^{I,n,\theta}_{\textbf{u}}+E^{A,n,\theta}_{\textbf{u}})-\mu(c^{n})\Delta(E^{I,n,\theta}_{\textbf{u}}+E^{A,n,\theta}_{\textbf{u}})+\bigtriangledown
E^{I,n,\theta}_{p}\\\ &\quad+\bigtriangledown E^{A,n,\theta}_{p}]\\\
&=\frac{\rho\tau_{1}^{\prime}}{dt-\rho\tau_{1}^{\prime}}I_{d\times
d}[\\{\rho\partial_{t}E^{I,n}_{\textbf{u}}+\rho(E^{I,n}_{\textbf{u}}\bigtriangledown)\textbf{u}^{n,\theta}+\rho(\textbf{u}_{h}^{n}\cdot\bigtriangledown)E^{I,n,\theta}_{\textbf{u}}-\mu(c^{n})\\\
&\quad\Delta E^{I,n,\theta}_{\textbf{u}}+\bigtriangledown
E^{I,n,\theta}_{p}\\}+\\{\rho\partial_{t}E^{A,n}_{\textbf{u}}+\rho(E^{A,n}_{\textbf{u}}\bigtriangledown)\textbf{u}^{n,\theta}+\rho(\textbf{u}_{h}^{n}\cdot\bigtriangledown)E^{A,n,\theta}_{\textbf{u}}\\\
&\quad-\mu(c^{n})\Delta E^{A,n,\theta}_{\textbf{u}}+\bigtriangledown
E^{A,n,\theta}_{p}\\}]\end{split}$ (51)

Since $\frac{\rho\tau_{1}}{dt+\rho\tau_{1}}<1$, which implies
$\frac{\rho\tau_{1}^{\prime}}{dt}<1$ and therefore the series
$\sum_{i=1}^{\infty}(\frac{\rho\tau_{1}^{\prime}}{dt})^{i}$ converges to
$\frac{\rho\tau_{1}^{\prime}}{dt-\rho\tau_{1}^{\prime}}$.
Similar to $\textbf{d}_{1}$, the other component $d_{3}$ is as follows:

$\begin{split}d_{3}&\leq\frac{\tau_{3}^{\prime}}{dt-\tau_{3}^{\prime}}[\partial_{t}E^{I,n}_{c}-\bigtriangledown\cdot\tilde{\bigtriangledown}E^{I,n,\theta}_{c}+(\textbf{u}_{h}^{n}\cdot\bigtriangledown)E^{I,n,\theta}_{c}+(E^{I,n}_{\textbf{u}}\cdot\bigtriangledown)c^{n,\theta}+\alpha
E^{I,n,\theta}_{c}\\\
&\quad+\partial_{t}E^{A,n}_{c}-\bigtriangledown\cdot\tilde{\bigtriangledown}E^{A,n,\theta}_{c}+(\textbf{u}_{h}^{n}\cdot\bigtriangledown)E^{A,n,\theta}_{c}+(E^{A,n}_{\textbf{u}}\cdot\bigtriangledown)c^{n,\theta}+\alpha
E^{A,n,\theta}_{c}]\end{split}$ (52)

It is clearly seen in the expansion of $\textbf{d}_{1}$ and $d_{3}$ that the
terms in $I_{3}$ and $I_{4}$ exactly match with the terms in $I_{2}$ and
$I_{1}$ respectively. Hence their estimations also follow the same way as we
have done earlier. Therefore skipping the repetition of presenting same
results, here we have mentioned the estimated results only for one term from
each of $I_{3}$ and $I_{4}$ in the following. Denoting first term of $I_{3}$
by the notation $I_{3}^{1}$ we have the estimated result as follows:

$\begin{split}I_{3}^{1}&=\sum_{k=1}^{n_{el}}(\frac{\rho\tau_{1k}}{dt+\rho\tau_{1k}}I_{d\times
d}\\{\rho\partial_{t}E^{I,n}_{\textbf{u}}+\rho(E^{I,n}_{\textbf{u}}\bigtriangledown)\textbf{u}^{n,\theta}+\rho(\textbf{u}_{h}^{n}\cdot\bigtriangledown)E^{I,n,\theta}_{\textbf{u}}\\\
&\quad-\mu(c^{n})\Delta E^{I,n,\theta}_{\textbf{u}}+\bigtriangledown
E^{I,n,\theta}_{p}\\},E^{A,n,\theta}_{\textbf{u}})_{\Omega_{k}}\\\ &\leq
h^{2}\frac{\rho^{2}CC_{\tau_{1}}}{(T_{0}-\rho
C_{\tau_{1}})}[\sum_{i=1}^{d}\\{\sum_{k=1}^{n_{el}}B_{1k}^{i}\\}(\|u_{i}^{n+1}\|_{2}+\|u_{i}^{n}\|_{2})]+\frac{\rho\mid\tau_{1}\mid}{(T_{0}-\rho
C_{\tau_{1}})}\sum_{i=1}^{d}(\sum_{k=1}^{n_{el}}B_{1k}^{i})\\\ &\quad[\\{\rho
h^{2}\|\frac{\partial u_{i}^{n,\theta}}{\partial x_{i}}\|+\rho
h\|u_{i}^{n}\|+\rho
hB^{i}_{1k}+\mu_{u}+h^{3}(\frac{1+\theta}{2}\|u_{i}^{n+1}\|_{2}+\frac{1-\theta}{2}\|u_{i}^{n}\|_{2})\\}\\\
&\quad(\frac{1+\theta}{2}\|u_{i}^{n+1}\|_{2}+\frac{1-\theta}{2}\|u_{i}^{n}\|_{2})+h(\frac{1+\theta}{2}\|p^{n+1}\|_{1}+\frac{1-\theta}{2}\|p^{n}\|_{1})]\end{split}$
(53)

and denoting first term of $I_{4}$ by $I_{4}^{1}$ we have the estimated result
in the following

$\begin{split}I_{4}^{1}&=\sum_{k=1}^{n_{el}}(\frac{\rho\tau_{1k}^{2}}{dt+\rho\tau_{1k}}I_{d\times
d}\\{\rho\partial_{t}E^{I,n}_{\textbf{u}}+\rho(E^{I,n}_{\textbf{u}}\cdot\bigtriangledown)\textbf{u}^{n,\theta}+\rho(\textbf{u}_{h}^{n}\cdot\bigtriangledown)E^{I,n,\theta}_{\textbf{u}}-\mu(c^{n})\\\
&\quad\Delta E^{I,n,\theta}_{\textbf{u}}+\bigtriangledown
E^{I,n,\theta}_{p}\\},\rho(\textbf{u}_{h}\cdot\bigtriangledown)E^{A,n,\theta}_{\textbf{u}}+\mu(c)\Delta
E^{A,n,\theta}_{\textbf{u}}+\bigtriangledown
E^{A,n,\theta}_{p})_{\Omega_{k}}\\\ &\leq
h^{2}\frac{\rho^{2}CC_{\tau_{1}}^{2}}{(T_{0}-\rho
C_{\tau_{1}})}(\sum_{k=1}^{n_{el}}D_{B_{1k}})\\{\|\textbf{u}^{n+1}\|_{2}+\|\textbf{u}^{n}\|_{2}\\}+\frac{\rho\mid\tau_{1}\mid^{2}}{(T_{0}-\rho
C_{\tau_{1}})}(\sum_{k=1}^{n_{el}}D_{B_{1k}})[\sum_{i=1}^{d}\\\ &\quad\\{\rho
h^{2}\|\frac{\partial u_{i}^{n,\theta}}{\partial x_{i}}\|+\rho
h\|u_{i}^{n}\|+\rho
hB^{i}_{1k}+\mu_{u}+h^{3}(\frac{1+\theta}{2}\|u_{i}^{n+1}\|_{2}+\frac{1-\theta}{2}\|u_{i}^{n}\|_{2})\\}\\\
&\quad(\frac{1+\theta}{2}\|u_{i}^{n+1}\|_{2}+\frac{1-\theta}{2}\|u_{i}^{n}\|_{2})+h(\frac{1+\theta}{2}\|p^{n+1}\|_{1}+\frac{1-\theta}{2}\|p^{n}\|_{1})]\end{split}$
(54)

Now the estimations of the remaining terms are quite obvious. Therefore we
directly add those results while combining them into (43) at last. Finally the
last term containing truncation error can be estimated as follows:

$\begin{split}(\textbf{TE}^{n,\theta},E^{A,n,\theta}_{\textbf{U}})&=(\textbf{TE}_{1}^{n,\theta},E^{A,n,\theta}_{\textbf{u}})+(TE_{2}^{n,\theta},E^{A,n,\theta}_{c})\\\
&\leq\frac{\epsilon_{9}}{2}(\|\textbf{TE}^{n,\theta}_{1}\|^{2}+\|TE^{n,\theta}_{2}\|^{2})+\frac{1}{2\epsilon_{9}}(\|E^{A,n,\theta}_{\textbf{u}}\|^{2}+\|E^{A,n,\theta}_{c}\|^{2})\end{split}$
(55)

This completes estimation of all the terms in the $RHS$ of (43). Now we start
with putting all the bounds, obtained for each of the terms in the right hand
side of (43). Then we take out few common terms in the left hand side and
consequently we have left with those terms multiplied by
$h^{2},\mid\tau_{1}\mid$ and $\mid\tau_{3}\mid$. Now we multiply both sides by
2 and taking integration over $(t^{n},t^{n+1})$ for $n$=0,1,…,$(N-1)$ to both
the sides. Finally we have (43) as follows:

$\\{1-\frac{2TC_{\tau_{3}}(1-C_{\tau_{3}})}{T_{0}+\rho
C_{\tau_{3}}}\sum_{k=1}^{n_{el}}D_{B_{3k}}-\frac{4TC_{\tau_{3}}}{T_{0}-\rho
C_{\tau_{3}}}\sum_{k=1}^{n_{el}}B_{2k}\\}\sum_{n=0}^{N-1}(\|E^{A,n+1}_{c}\|^{2}-\|E^{A,n}_{c}\|^{2})\\\
+\rho\\{1-\frac{2(1+\rho C_{\tau_{1}})TC_{\tau_{1}}}{T_{0}-\rho
C_{\tau_{1}}}\sum_{k=1}^{n_{el}}D_{B_{1k}}-\frac{2(1+\rho)TC_{\tau_{1}}}{T_{0}-\rho
C_{\tau_{1}}}\sum_{i=1}^{d}\sum_{k=1}^{n_{el}}B_{1k}^{i}\\}\sum_{n=0}^{N-1}(\|E^{A,n+1}_{\textbf{u}}\|^{2}\\\
-\|E^{A,n}_{\textbf{u}}\|^{2})+\\{2\mu_{l}-\frac{1}{\epsilon_{1}}-\frac{\mu_{u}}{\epsilon_{2}}-\frac{4C_{2}}{\epsilon_{8}}-4C_{2}^{\prime}-2C_{\tau_{2}}\\}\sum_{n=0}^{N-1}\int_{t^{n}}^{t^{n+1}}\mid
E^{A,n,\theta}_{\textbf{u}}\mid_{1}^{2}dt\\\
+\\{2\sigma-\frac{\sigma}{\epsilon_{4}}-\frac{2C_{2}}{\epsilon_{7}}-\frac{4C_{2}}{\epsilon_{8}}-4C_{2}^{\prime}-\frac{1}{\epsilon_{9}}\\}\sum_{n=0}^{N-1}\int_{t^{n}}^{t^{n+1}}\|E^{A,n,\theta}_{\textbf{u}}\|^{2}dt\\\
+\\{2D_{l}-\frac{D_{m}}{\epsilon_{2}}-4\bar{C}_{1}^{2}\\}\sum_{n=0}^{N-1}\int_{t^{n}}^{t^{n+1}}\mid
E^{A,n,\theta}_{c}\mid_{1}^{2}dt+\\\ $ (56)
$\\{2\alpha-\frac{\alpha}{\epsilon_{5}}-\frac{4\bar{C}_{1}^{2}}{\epsilon_{6}}-\frac{1}{\epsilon_{9}}\\}\sum_{n=0}^{N-1}\int_{t^{n}}^{t^{n+1}}\|E^{A,n,\theta}_{c}\|^{2}dt\\\
\leq
h^{2}\sum_{n=0}^{N-1}\int_{t^{n}}^{t^{n+1}}[2C^{2}\epsilon_{1}(\frac{1+\theta}{2}\|p^{n+1}\|_{2}+\frac{1-\theta}{2}\|p^{n}\|_{2})^{2}+\\{2C^{2}\mu_{u}\epsilon_{2}+h^{2}\sigma\epsilon_{4}+2C_{2}\epsilon_{7}\\\
+3C_{2}h^{2}\epsilon_{8}+C_{\tau_{2}}\\}(\frac{1+\theta}{2}\|\textbf{u}^{n+1}\|_{2}+\frac{1-\theta}{2}\|\textbf{u}^{n}\|_{2})^{2}+\\{C^{2}D_{m}\epsilon_{3}+h^{2}\alpha\epsilon_{5}+\bar{C}_{1}^{2}C^{2}\epsilon_{6}\\}\\\
(\frac{1+\theta}{2}\|c^{n+1}\|_{2}+\frac{1-\theta}{2}\|c^{n}\|_{2})^{2}+2\frac{\rho
CC_{\tau_{1}}}{T_{0}-\rho C_{\tau_{1}}}\\{(1+\rho
C_{\tau_{1}})(\sum_{k=1}^{n_{el}}D_{B_{1k}})+(1+\rho)\\\
(\sum_{i=1}^{d}\sum_{k=1}^{n_{el}}B_{1k}^{i})\\}(\|\textbf{u}^{n+1}\|_{2}+\|\textbf{u}^{n}\|_{2})+\frac{2CC_{\tau_{3}}}{T_{0}-C_{\tau_{3}}}\\{\sum_{k=1}^{n_{el}}D_{B_{3k}}+(1+C_{\tau_{3}})\sum_{k=1}^{n_{el}}B_{2k}\\}\\\
(\|c^{n+1}\|_{2}+\|c^{n}\|_{2})]dt+2\mid\tau_{1}\mid\sum_{n=0}^{N-1}\int_{t^{n}}^{t^{n+1}}[\frac{TC+\rho
C_{\tau_{1}}}{T_{0}-\rho
C_{\tau_{1}}}\sum_{k=1}^{n_{el}}D_{B_{1k}}+\frac{1+\rho}{T_{0}-\rho
C_{\tau_{1}}}\sum_{i=1}^{d}\\\ \sum_{k=1}^{n_{el}}B_{1k}^{i}][\\{\rho
h^{2}\|\frac{\partial u_{i}^{n,\theta}}{\partial x_{i}}\|+\rho
h\|u_{i}^{n,\epsilon\theta}\|+\rho
hB^{i}_{1k}+\mu_{u}+h^{3}(\frac{1+\theta}{2}\|u_{i}^{n+1}\|_{2}+\frac{1-\theta}{2}\|u_{i}^{n}\|_{2})\\}\\\
(\frac{1+\theta}{2}\|u_{i}^{n+1}\|_{2}+\frac{1-\theta}{2}\|u_{i}^{n}\|_{2})+h(\frac{1+\theta}{2}\|p^{n+1}\|_{1}+\frac{1-\theta}{2}\|p^{n}\|_{1})]dt+\\\
\mid\tau_{1}\mid\sum_{n=0}^{N-1}\int_{t^{n}}^{t^{n+1}}\frac{2T}{T_{0}-\rho
C_{\tau_{1}}}[(1+\rho)\sum_{k=1}^{n_{el}}D_{B_{1k}}\bar{D}_{B_{1k}}+(1+\rho
C_{\tau_{1}})\sum_{i=1}^{d}\sum_{k=1}^{n_{el}}\bar{D}_{B_{1k}}B_{1k}^{i}]dt\\\
+2\mid\tau_{3}\mid\sum_{n=0}^{N-1}\int_{t^{n}}^{t^{n+1}}[\frac{2C}{T_{0}-C_{\tau_{3}}}\sum_{k=1}^{n_{el}}D_{B_{3k}}+\frac{1+C_{\tau_{3}}}{T_{0}-C_{\tau_{3}}}\sum_{k=1}^{n_{el}}B_{2k}][\sum_{i=1}^{d}\\{h^{2}\|\frac{\partial
c^{n,\theta}}{\partial x_{i}}\|\\\
(\frac{1+\theta}{2}\|u_{i}^{n+1}\|_{2}+\frac{1-\theta}{2}\|u_{i}^{n}\|_{2})+(D_{im}+h\bar{D}_{im}+\mid\alpha\mid
h^{2}+h\|u_{i}^{n,\theta}\|+hB_{1k}^{i}\\\
+h^{3}(\frac{1+\theta}{2}\|u_{i}^{n+1}\|_{2}+\frac{1-\theta}{2}\|u_{i}^{n}\|_{2}))(\frac{1+\theta}{2}\|c^{n+1}\|_{2}+\frac{1-\theta}{2}\|c^{n}\|_{2})\\}]dt\\\
+\mid\tau_{3}\mid\sum_{n=0}^{N-1}\int_{t^{n}}^{t^{n+1}}\frac{2C}{T_{0}-C_{\tau_{3}}}[\sum_{k=1}^{n_{el}}\bar{D}_{B_{3k}}D_{B_{3k}}+(1+C_{\tau_{3}})\sum_{k=1}^{n_{el}}\bar{D}_{B_{3k}}B_{2k}]dt\\\
+\epsilon_{9}\sum_{n=0}^{N-1}\int_{t^{n}}^{t^{n+1}}(\|\textbf{TE}^{n,\theta}_{1}\|^{2}+\|TE^{n,\theta}_{2}\|^{2})dt\hskip
119.50157pt$ (57)

We can choose the values of the arbitrary parameters in such a manner that we
can make all the coefficients in the left hand side positive. In order to
satisfy such condition it is inevitable to choose the characteristic lengths
small. Now after taking minimum of all the coefficients in left hand side, let
us divide both the sides with that minimum, which turns out to be a positive
real number. Applying assumption (iv) it can be seen that
$\|\textbf{u}^{n}\|_{2}$, $\|p^{n}\|_{1}$ and $\|c^{n}\|_{2}$ are bounded for
$n=0,1,2,...,N$. The initial conditions considered in section 3.1, imply
$\|E^{A,0}_{\textbf{u}}\|=0$ and $\|E^{A,0}_{c}\|=0$.
After performing all these intermediate steps and applying the properties
(15)-(16) on truncation errors we finally arrive at the following expression
since $\tau_{1}$ and $\tau_{3}$ are of order $h^{2}$:

$\|E^{A,N}_{\textbf{u}}\|^{2}+\|E^{A,N}_{c}\|^{2}+\sum_{n=0}^{N-1}\int_{t^{n}}^{t^{n+1}}\|E^{A,n,\theta}_{\textbf{u}}\|_{1}^{2}dt+\sum_{n=0}^{N-1}\int_{t^{n}}^{t^{n+1}}\|E^{A,n,\theta}_{c}\|^{2}_{1}dt\\\
\leq C(T,\textbf{u},p,c)(h^{2}+dt^{2r})\hskip 85.35826pt$ (58)

This implies

$\boxed{\|E^{A}_{\textbf{u}}\|_{\tilde{\textbf{V}}}^{2}+\|E^{A}_{c}\|_{\tilde{\textbf{V}}}^{2}\leq
C(T,\textbf{u},p,c)(h^{2}+dt^{2r})}$ (59)

where

$r=\begin{cases}1,&\text{if}\ \theta=1\\\ 2,&\text{if}\ \theta=0\end{cases}$
(60)

We have used the fact that $\sum_{n=0}^{N-1}\int_{t^{n}}^{t^{n+1}}Mdt\leq MT$.
This completes the first part of the proof.
Second part Using this above result we are going to estimate auxiliary error
part of pressure. We will use inf-sup condition to find estimate for
$E_{p}^{A}$. Applying Galerkin orthogonality only for variational form of
Navier-Stokes flow problem we have obtained

$(\frac{\partial(\textbf{u}-\textbf{u}_{h})}{\partial
t},\textbf{v}_{h})+c(\textbf{u},\textbf{u},\textbf{v}_{h})-c(\textbf{u}_{h},\textbf{u}_{h},\textbf{v}_{h})+a_{NS}(\textbf{u}-\textbf{u}_{h},\textbf{v}_{h})-b(\textbf{v}_{h},p-p_{h})=0\\\
$ (61)

Splitting of the errors implies the following

$b(\textbf{v}_{h},p-I_{h}p)+b(\textbf{v}_{h},I_{h}p-p_{h})=(\partial_{t}E^{A}_{\textbf{u}},\textbf{v}_{h})+c(E^{I}_{\textbf{u}},\textbf{u},\textbf{v}_{h})+c(E^{A}_{\textbf{u}},\textbf{u},\textbf{v}_{h})+\\\
c(\textbf{u}_{h},E^{I}_{\textbf{u}},\textbf{v}_{h})+c(\textbf{u}_{h},E^{A}_{\textbf{u}},\textbf{v}_{h})+a_{NS}(E^{I}_{\textbf{u}},\textbf{v}_{h})+a_{NS}(E^{A}_{\textbf{u}},\textbf{v}_{h})$
(62)

Without loss of generality considering the inclusion $\bigtriangledown\cdot
V^{h}\subset Q^{h}$ and the property of the $L^{2}$ orthogonal projection of
$I^{h}_{p}$ we have

$b(\textbf{v}_{h},p-I_{h}p)=\int_{\Omega}(p-I_{h}p)(\bigtriangledown\cdot\textbf{v}_{h})=0$
(63)

Now according to inf-sup condition we will have the following expression

$\begin{split}\|I_{h}p-p_{h}\|_{L^{2}(L^{2})}^{2}&=\|E_{p}^{A}\|_{L^{2}(L^{2})}^{2}\\\
&=\sum_{n=0}^{N-1}\int_{t^{n}}^{t^{n+1}}\|E_{p}^{A,n,\theta}\|^{2}dt\\\
&\leq\sum_{n=0}^{N-1}\int_{t^{n}}^{t^{n+1}}\underset{\textbf{v}_{h}}{sup}\frac{b(\textbf{v}_{h},E_{p}^{A,n,\theta})}{\|\textbf{v}_{h}\|_{1}}dt\end{split}$
(64)

Using (74) on (73) we will have

$\begin{split}\sum_{n=0}^{N-1}\int_{t^{n}}^{t^{n+1}}b(\textbf{v}_{h},E_{p}^{A,n,\theta})dt&=\sum_{n=0}^{N-1}\int_{t^{n}}^{t^{n+1}}\\{(\frac{E^{A,n+1}_{\textbf{u}}-E^{A,n}_{\textbf{u}}}{dt},\textbf{v}_{h})+a_{NS}(E^{I,n,\theta}_{\textbf{u}},\textbf{v}_{h})\\\
&\quad+c(E^{I,n}_{\textbf{u}},\textbf{u}^{n,\theta},\textbf{v}_{h})+c(E^{A,n}_{\textbf{u}},\textbf{u}^{n,\theta},\textbf{v}_{h})+\\\
&\quad
c(\textbf{u}_{h}^{n},E^{I,n,\theta}_{\textbf{u}},\textbf{v}_{h})+c(\textbf{u}_{h}^{n},E^{A,n,\theta}_{\textbf{u}},\textbf{v}_{h})+\\\
&\quad a_{NS}(E^{A,n,\theta}_{\textbf{u}},\textbf{v}_{h})\\}dt\end{split}$
(65)

Now applying the results obtained in the previous part we will have,

$\sum_{n=0}^{N-1}\int_{t^{n}}^{t^{n+1}}b(\textbf{v}_{h},E_{p}^{A,n,\theta})dt\leq
C(T,\textbf{u},p,c)(h^{2}+dt^{2r})\|\textbf{v}_{h}\|_{1}$ (66)

Using this above result into (75), we will have the estimate for the pressure
term

$\|I_{h}p-p_{h}\|_{L^{2}(L^{2})}^{2}\leq C(T,\textbf{u},p,c)(h^{2}+dt^{2r})$
(67)

Now combining the results obtained in the first and second part we have
finally arrived at the following auxiliary error estimate:

$\|E^{A}_{\textbf{u}}\|^{2}_{\tilde{\textbf{V}}}+\|E^{A}_{p}\|_{L^{2}(L^{2})}^{2}+\|E^{A}_{c}\|^{2}_{\tilde{\textbf{V}}}\leq
C(T,\textbf{u},p,c)(h^{2}+dt^{2r})$ (68)

where

$r=\begin{cases}1,&\text{if}\ \theta=1\\\ 2,&\text{if}\ \theta=0\end{cases}$
(69)

This completes the proof. ∎

######  Theorem 2.

(Apriori error estimate) Assuming the same condition as in the previous
theorem,

$\|\textbf{u}-\textbf{u}_{h}\|_{\tilde{\textbf{V}}}^{2}+\|p-p_{h}\|_{L^{2}(L^{2})}^{2}+\|c-c_{h}\|^{2}_{\tilde{\textbf{V}}}\leq
C^{\prime}(h^{2}+dt^{2r})$ (70)

where $C^{\prime}$ depends on T, u,p,c and

$r=\begin{cases}1,&\text{if}\ \theta=1\\\ 2,&\text{if}\ \theta=0\end{cases}$
(71)

###### Proof.

By applying triangle inequality, the interpolation inequalities and the result
of the previous theorem we will have,

$\|\textbf{u}-\textbf{u}_{h}\|_{\tilde{\textbf{V}}}^{2}+\|p-p_{h}\|_{L^{2}(L^{2})}^{2}+\|c-c_{h}\|^{2}_{\tilde{\textbf{V}}}\\\
=\|E^{I}_{\textbf{u}}+E^{A}_{\textbf{u}}\|_{\tilde{\textbf{V}}}^{2}+\|E^{I}_{p}+E^{A}_{p}\|_{L^{2}(L^{2})}^{2}+\|E^{I}_{c}+E^{A}_{c}\|_{\tilde{\textbf{V}}}^{2}\hskip
88.20354pt\\\
\leq\bar{C}(\|E^{I}_{\textbf{u}}\|_{\tilde{\textbf{V}}}^{2}+\|E^{I}_{p}\|_{L^{2}(L^{2})}^{2}+\|E^{I}_{c}\|_{\tilde{\textbf{V}}}^{2}+\|E^{A}_{\textbf{u}}\|_{\tilde{\textbf{V}}}^{2}+\|E^{A}_{p}\|_{L^{2}(L^{2})}^{2}+\|E^{A}_{c}\|_{\tilde{\textbf{V}}}^{2})\\\
\leq C^{\prime}(T,\textbf{u},p,c)(h^{2}+dt^{2r})\hskip 184.9429pt$ (72)

This completes $apriori$ error estimation. ∎

###  4.3 Aposteriori error estimation

In this section we are going to derive residual based aposteriori error
estimation. This estimation is also comprised of two parts similar to the
auxiliary apriori error estimate derived in the earlier section.

######  Theorem 3.

For computed velocity $\textbf{u}_{h}$, pressure $p_{h}$ and concentration
$c_{h}$ belonging to $(V^{h})^{d}\times Q^{h}\times V^{h}$ satisfying
(15)-(16), assume $dt$ is sufficiently small and positive, and sufficient
regularity of exact solution in equations (1)-(2). Then there exists a
constant $\bar{C}$, independent of $\textbf{u},p,c$ and depending on the
residual such that

$\|\textbf{u}-\textbf{u}_{h}\|^{2}_{\tilde{\textbf{V}}}+\|p-p_{h}\|_{L^{2}(L^{2})}^{2}+\|c-c_{h}\|^{2}_{\tilde{\textbf{V}}}\leq\bar{C}(\textbf{R})(h^{2}+dt^{2r})$
(73)

where R is the residual vector and

$r=\begin{cases}1,&\text{if}\ \theta=1\\\ 2,&\text{if}\ \theta=0\end{cases}$
(74)

###### Proof.

We estimate $aposteriori$ error by dividing the procedure into two parts. In
the first part we find error bound corresponding to $velocity$ and
$concentration$ followed by the second part estimating error associated with
the $pressure$ term. Let us first introduce the residual vector corresponding
to each equations

$\textbf{R}=\begin{bmatrix}\textbf{f}-\\{\rho\frac{\partial\textbf{u}_{h}}{\partial
t}+\rho(\textbf{u}_{h}\cdot\bigtriangledown)\textbf{u}_{h}-\mu(c)\Delta\textbf{u}_{h}+\bigtriangledown
p_{h}\\}\\\ -\bigtriangledown\cdot\textbf{u}_{h}\\\ g-(\frac{\partial
c_{h}}{\partial
t}-\bigtriangledown\cdot\tilde{\bigtriangledown}c_{h}+\textbf{u}\cdot\bigtriangledown
c_{h}+\alpha c_{h})\end{bmatrix}=\begin{bmatrix}\textbf{R}_{1}\\\ R_{2}\\\
R_{3}\end{bmatrix}$

First part: We have $\forall\textbf{V}\in\textbf{V}_{F}$

$\mu_{l}\mid\textbf{v}\mid_{1}^{2}+D_{\alpha}\|d\|_{1}^{2}\leq
B(\textbf{u};\textbf{V},\textbf{V})=a_{NS}(\textbf{v},\textbf{v})+a_{T}(d,d)$
(75)

Since $\textbf{e}\in\textbf{V}_{F}$ we substitute the errors
$e_{\textbf{u}},e_{c}$ into the above relation:

$\mu_{l}\|e_{\textbf{u}}\|_{1}^{2}+D_{\alpha}\|e_{c}\|_{1}^{2}\leq
a_{NS}(e_{\textbf{u}},e_{\textbf{u}})+a_{T}(e_{c},e_{c})+\mu_{l}\|e_{\textbf{u}}\|^{2}\\\
$ (76)

By adding few terms in both sides the above equation becomes

$(\frac{\partial e_{\textbf{u}}}{\partial t},e_{\textbf{u}})+(\frac{\partial
e_{c}}{\partial
t},e_{c})+\mu_{l}\|e_{\textbf{u}}\|_{1}^{2}+D_{\alpha}\|e_{c}\|_{1}^{2}\\\
\leq(\frac{\partial e_{\textbf{u}}}{\partial
t},e_{\textbf{u}})+(\frac{\partial e_{c}}{\partial
t},e_{c})+c(\textbf{u},e_{\textbf{u}},e_{\textbf{u}})+a_{NS}(e_{\textbf{u}},e_{\textbf{u}})+a_{T}(e_{c},e_{c})\hskip
28.45274pt\\\
+b(e_{\textbf{u}},e_{p})-b(e_{\textbf{u}},e_{p})+\mu_{l}\|e_{\textbf{u}}\|^{2}$
(77)

Now first we will find a lower bound of $LHS$ and then upper bound for $RHS$
and finally combining them we will get $aposteriori$ error estimate. To find
the lower bound the $LHS$ can be written as

$LHS=(\frac{e^{n+1}_{\textbf{u}}-e^{n}_{\textbf{u}}}{dt},e_{\textbf{u}}^{n,\theta})+(\frac{e_{c}^{n+1}-e_{c}^{n}}{dt},e_{c}^{n,\theta})+\mu_{l}\|e_{\textbf{u}}^{n,\theta}\|_{1}^{2}+D_{\alpha}\|e_{c}^{n,\theta}\|^{2}_{1}$
(78)

Applying (32) on first two terms of $LHS$ we have the following relations

$\begin{split}(\frac{e_{\textbf{u}}^{n+1}-e_{\textbf{u}}^{n}}{dt},e_{\textbf{u}}^{n,\theta})&\geq\frac{1}{2dt}(\|e_{\textbf{u}}^{n+1}\|^{2}-\|e_{\textbf{u}}^{n}\|^{2})\\\
\end{split}$ (79)

and

$\begin{split}(\frac{e_{c}^{n+1}-e_{c}^{n}}{dt},e_{c}^{n,\theta})&\geq\frac{1}{2dt}(\|e_{c}^{n+1}\|^{2}-\|e_{c}^{n}\|^{2})\\\
\end{split}$ (80)

Hence

$\frac{1}{2dt}(\|e_{\textbf{u}}^{n+1}\|^{2}-\|e_{\textbf{u}}^{n}\|^{2})+\frac{1}{2dt}(\|e_{c}^{n+1}\|^{2}-\|e_{c}^{n}\|^{2})+\mu_{l}\|e_{\textbf{u}}^{n,\theta}\|_{1}^{2}+D_{\alpha}\|e_{c}^{n,\theta}\|^{2}_{1}\\\
\leq LHS\leq RHS$ (81)

Now our job is to find upper bound for $RHS$ and to reach at the desired
estimates let us divide it into two broad parts by splitting errors in each of
the terms in the following way:

$\begin{split}RHS&=[(\frac{e_{\textbf{u}}^{n+1}-e_{\textbf{u}}^{n}}{dt},E_{\textbf{u}}^{I,n,\theta})+(\frac{e_{c}^{n+1}-e_{c}^{n}}{dt},E^{I,n,\theta}_{c})+c(\textbf{u}^{n},e_{\textbf{u}}^{n,\theta},E^{I,n,\theta}_{\textbf{u}})+\\\
&\quad
a_{NS}(e^{n,\theta}_{\textbf{u}},E^{I,n,\theta}_{\textbf{u}})+b(e^{n,\theta}_{\textbf{u}},E^{I,n,\theta}_{p})-b(E^{I,n,\theta}_{\textbf{u}},e^{n,\theta}_{p})+a_{T}(e^{n,\theta}_{c},E^{I,n,\theta}_{c})]+\\\
&\quad[(\frac{e_{\textbf{u}}^{n+1}-e_{\textbf{u}}^{n}}{dt},E_{\textbf{u}}^{A,n,\theta})+(\frac{e^{n+1}_{c}-e_{c}^{n}}{dt},E^{A,n,\theta}_{c})+c(\textbf{u}^{n},e_{\textbf{u}}^{n,\theta},E^{A,n,\theta}_{\textbf{u}})+\\\
&\quad
a_{NS}(e^{n,\theta}_{\textbf{u}},E^{A,n,\theta}_{\textbf{u}})+b(e^{n,\theta}_{\textbf{u}},E^{A,n,\theta}_{p})-b(E^{A,n,\theta}_{\textbf{u}},e^{n,\theta}_{p})+a_{T}(e^{n,\theta}_{c},E^{A,n,\theta}_{c})]\\\
&\quad+\mu_{l}\|e_{\textbf{u}}^{n,\theta}\|^{2}\\\
&=RHS^{I}+RHS^{A}+\mu_{l}\|e_{\textbf{u}}^{n,\theta}\|^{2}\\\ \end{split}$
(82)

Our aim is to bring residual into context and for this purpose $RHS^{I}$
involving interpolation error terms can be written as follows:

$\begin{split}RHS^{I}&=[\rho(\frac{e_{\textbf{u}}^{n+1}-e_{\textbf{u}}^{n}}{dt},E^{I,n,\theta}_{\textbf{u}})+(\frac{e_{c}^{n+1}-e_{c}^{n}}{dt},E^{I,n,\theta}_{c})+c(\textbf{u}^{n},\textbf{u}^{n,\theta},E^{I,n,\theta}_{\textbf{u}})-\\\
&\quad
c(\textbf{u}^{n}_{h},\textbf{u}^{n,\theta}_{h},E^{I,n,\theta}_{\textbf{u}})+a_{NS}(e^{n,\theta}_{\textbf{u}},E^{I,n,\theta}_{\textbf{u}})-b(E^{I,n,\theta}_{\textbf{u}},e^{n,\theta}_{p})+b(e^{n,\theta}_{\textbf{u}},E^{I,n,\theta}_{p})\\\
&\quad+a_{T}(e_{c}^{n,\theta},E^{I,n,\theta}_{c})]-c(e_{\textbf{u}}^{n},\textbf{u}_{h}^{n,\theta},E^{I,n,\theta}_{\textbf{u}})\\\
&=RHS^{I}_{1}-c(e_{\textbf{u}}^{n},\textbf{u}_{h}^{n,\theta},E^{I,n,\theta}_{\textbf{u}})\end{split}$
(83)

The bracketed term in the above equation is denoted by $RHS^{I}_{1}$.
$RHS^{A}$ involving auxiliary part of error can also be decomposed in the
similar manner as above and let us denote the alike term corresponding to
$RHS^{A}$ by $RHS^{A}_{1}$. Therefore

$RHS^{A}=RHS^{A}_{1}-c(e_{\textbf{u}}^{n},\textbf{u}_{h}^{n,\theta},E^{A,n,\theta}_{\textbf{u}})\\\
$ (84)

Hence combining these above two results (93) becomes

$\begin{split}RHS&=RHS^{I}_{1}+RHS^{A}_{1}-c(e_{\textbf{u}}^{n},\textbf{u}_{h}^{n,\theta},e_{\textbf{u}}^{n,\theta})+\mu_{l}\|e_{\textbf{u}}^{n,\theta}\|^{2}\\\
&=RHS^{I}_{1}+RHS^{A}_{1}-c(e_{\textbf{u}}^{n},\textbf{u}^{n,\theta},e_{\textbf{u}}^{n,\theta})+c(e_{\textbf{u}}^{n},e_{\textbf{u}}^{n,\theta},e_{\textbf{u}}^{n,\theta})+\mu_{l}\|e_{\textbf{u}}^{n,\theta}\|^{2}\\\
&=RHS^{I}_{1}+RHS^{A}_{1}-c(e_{\textbf{u}}^{n},\textbf{u}^{n,\theta},e_{\textbf{u}}^{n,\theta})+\mu_{l}\|e_{\textbf{u}}^{n,\theta}\|^{2}\end{split}$
(85)

Property (a) of trilinear term $c(\cdot,\cdot,\cdot)$ implies
$c(e_{\textbf{u}}^{n,\epsilon\theta},e_{\textbf{u}}^{n,\theta},e_{\textbf{u}}^{n,\theta})=0$.
In the most general way for $Navier$-$Stokes$ flow problem we have for all
$\textbf{v}\in(V)^{d}$

$\rho(\frac{e_{\textbf{u}}^{n+1}-e_{\textbf{u}}^{n}}{dt},\textbf{v})+c(\textbf{u}^{n},\textbf{u}^{n,\theta},\textbf{v})-c(\textbf{u}^{n}_{h},\textbf{u}^{n,\theta}_{h},\textbf{v})+a_{NS}(e^{n,\theta}_{\textbf{u}},\textbf{v})-b(\textbf{v},e_{p}^{n,\theta})\\\
=\int_{\Omega}\textbf{R}_{1}^{n,\theta}\cdot\textbf{v}\hskip 256.0748pt\\\
Similarly\hskip 5.69054pt\int_{\Omega}(\bigtriangledown\cdot
e_{\textbf{u}}^{n,\theta})q=\int_{\Omega}R_{2}^{n,\theta}q\hskip
5.69054pt\forall q\in Q\hskip 133.72795pt\\\
\int_{\Omega}(\frac{e_{c}^{n+1}-e_{c}^{n}}{dt}d+\tilde{\bigtriangledown}e_{c}^{n,\theta}\cdot\bigtriangledown
d+d\textbf{u}^{n}\cdot\bigtriangledown e_{c}^{n,\theta}+\alpha
e_{c}^{n,\theta}d)=\int_{\Omega}R_{3}^{n,\theta}d\hskip 8.53581pt\forall d\in
V$ (86)

Now substituting $\textbf{v},q,d$ in the above expressions by
$E^{I,n,\theta}_{\textbf{u}},E^{I,n,\theta}_{p},E^{I,n,\theta}_{c}$
respectively, we have the $RHS^{I}_{1}$ as,

$\begin{split}RHS^{I}_{1}&=\int_{\Omega}(\textbf{R}_{1}^{n,\theta}\cdot
E^{I,n,\theta}_{\textbf{u}}+R_{2}^{n,\theta}E^{I,n,\theta}_{p}+R_{3}^{n,\theta}E^{I,n,\theta}_{c})\\\
&\leq
h^{2}\\{\|\textbf{R}_{1}^{n,\theta}\|(\frac{1+\theta}{2}\|\textbf{u}^{n+1}\|_{2}+\frac{1-\theta}{2}\|\textbf{u}^{n}\|_{2})+C_{2}\|R_{2}^{n,\theta}\|(\frac{1+\theta}{2}\|p^{n+1}\|_{1}\\\
&\quad+\frac{1-\theta}{2}\|p^{n}\|_{1})+\|R_{3}^{n,\theta}\|(\frac{1+\theta}{2}\|c_{n+1}\|_{2}+\frac{1-\theta}{2}\|c^{n}\|_{2})\\}\\\
&\leq
h^{2}(\bar{C}_{1}\|\textbf{R}_{1}^{n,\theta}\|+\bar{C}_{2}\|R_{2}^{n,\theta}\|+\bar{C}_{3}\|R_{3}^{n,\theta}\|)\\\
\end{split}$ (87)

The parameters $\bar{C}_{i}$, for i=1,2,3,4, are coming from imposing
assumption (iv). Now we are going to estimate of $RHS^{A}_{1}$. For that we
employ $subgrid$ formulation (8). Subtracting (8) from the variational finite
element formulation satisfied by the exact solution we have
$\forall\textbf{V}_{h}\in\textbf{V}_{F}^{h}$

$\rho(\frac{e_{\textbf{u}}^{n+1}-e_{\textbf{u}}^{n}}{dt},\textbf{v}_{h})+(\frac{e_{c}^{n+1}-e_{c}^{n}}{dt},d_{h})+c(\textbf{u}^{n},\textbf{u}^{n,\theta},\textbf{v}_{h})-c(\textbf{u}^{n}_{h},\textbf{u}^{n,\theta}_{h},\textbf{v}_{h})\\\
+a_{NS}(e_{\textbf{u}}^{n,\theta},\textbf{v}_{h})-b(\textbf{v}_{h},e^{n,\theta}_{p})+b(e_{\textbf{u}}^{n,\theta},q_{h})+a_{T}(e_{c}^{n,\theta},d_{h})\\\
=\sum_{k=1}^{n_{el}}\\{(\tau_{k}^{\prime}(\textbf{R}^{n,\theta}+\textbf{d}),-\mathcal{L}^{*}\textbf{V}_{h})_{\Omega_{k}}-((I-\tau_{k}^{-1}\tau_{k})\textbf{R}^{n,\theta},\textbf{V}_{h})_{\Omega_{k}}+(\tau_{k}^{-1}\tau_{k}\textbf{d},\textbf{V}_{h})_{\Omega_{k}}\\}\\\
+(\textbf{TE}^{n,\theta}_{1},\textbf{v}_{h})+(TE^{n,\theta}_{2},d_{h})\\\ $
(88) $=\sum_{k=1}^{n_{el}}[(\tau_{1}^{\prime}I_{d\times
d}\\{\textbf{R}_{1}^{n,\theta}+\textbf{d}_{1}\\},\rho(\textbf{u}_{h}\cdot\bigtriangledown)\textbf{v}_{h}+\mu(c)\Delta\textbf{v}_{h}+\bigtriangledown
q_{h})_{\Omega_{k}}+\tau_{2}^{\prime}(R_{2}^{n,\theta},\bigtriangledown\cdot\textbf{v}_{h})_{\Omega_{k}}\\\
\quad+\tau_{3}^{\prime}(R_{3}^{n,\theta}+d_{3},\bigtriangledown\cdot\tilde{\bigtriangledown}d_{h}+\textbf{u}_{h}\cdot\bigtriangledown
d_{h}-\alpha d_{h})_{\Omega_{k}}+((1-\tau_{1}^{-1}\tau_{1}^{\prime})I_{d\times
d}\textbf{R}_{1}^{n,\theta},\textbf{v}_{h})_{\Omega_{k}}\\\
+((1-\tau_{3}^{-1}\tau_{3}^{\prime})R_{3}^{n,\theta},d_{h})_{\Omega_{k}}+(\tau_{1}^{-1}\tau_{1}^{\prime}I_{d\times
d}\textbf{d}_{1},\textbf{v}_{h})_{\Omega_{k}}+\tau_{3}^{-1}\tau_{3}^{\prime}(d_{3},d_{h})_{\Omega_{k}}]\\\
+(\textbf{TE}_{1}^{n,\theta},\textbf{v}_{h})+(TE_{2}^{n,\theta},d_{h})\hskip
156.49014pt$ (89)

Now substituting $\textbf{V}_{h}$ by
$(E^{A,n,\theta}_{\textbf{u}},E^{A,n,\theta}_{p},E^{A,n,\theta}_{c})$ in the
above equation we have $RHS^{A}_{1}$ as follows

$\begin{split}RHS^{A}_{1}&=\sum_{k=1}^{n_{el}}[(\tau_{1}^{\prime}I_{d\times
d}\\{\textbf{R}_{1}^{n,\theta}+\textbf{d}_{1}\\},\rho(\textbf{u}_{h}\cdot\bigtriangledown)E^{A,n,\theta}_{\textbf{u}}+\mu(c)\Delta
E^{A,n,\theta}_{\textbf{u}}+\bigtriangledown
E^{A,n,\theta}_{p})_{\Omega_{k}}\\\
&\quad+\tau_{2}^{\prime}(R_{2}^{n,\theta},\bigtriangledown\cdot
E^{A,n,\theta}_{\textbf{u}})_{\Omega_{k}}+\tau_{3}^{\prime}(R_{3}^{n,\theta}+d_{4},\bigtriangledown\cdot\tilde{\bigtriangledown}E^{A,n,\theta}_{c}+\textbf{u}_{h}\cdot\bigtriangledown
E^{A,n,\theta}_{c}\\\ &\quad-\alpha
E^{A,n,\theta}_{c})_{\Omega_{k}}+((1-\tau_{1}^{-1}\tau_{1}^{\prime})I_{d\times
d}\textbf{R}_{1}^{n,\theta},E^{A,n,\theta}_{\textbf{u}})_{\Omega_{k}}+(1-\tau_{3}^{-1}\tau_{3}^{\prime})(R_{3}^{n,\theta},\\\
&\quad
E^{A,n,\theta}_{c})_{\Omega_{k}}+(\tau_{1}^{-1}\tau_{1}^{\prime}I_{d\times
d}\textbf{d}_{1},E^{A,n,\theta}_{\textbf{u}})_{\Omega_{k}}+\tau_{3}^{-1}\tau_{3}^{\prime}(d_{3},E^{A,n,\theta}_{c})_{\Omega_{k}}]+\\\
&\quad(\textbf{TE}_{1}^{n,\theta},E^{A,n,\theta}_{\textbf{u}})+(TE_{2}^{n,\theta},E^{A,n,\theta}_{c})\end{split}$
(90)

Now we estimate each term separately. We use the results mentioned earlier
during $apriori$ error estimation.

$\begin{split}\sum_{k=1}^{n_{el}}(\tau_{1}^{\prime}I_{d\times
d}\textbf{R}_{1}^{n,\theta},\rho(\textbf{u}_{h}\cdot\bigtriangledown)E^{A,n,\theta}_{\textbf{u}}+\mu(c)\Delta
E^{A,n,\theta}_{\textbf{u}}+\bigtriangledown
E^{A,n,\theta}_{p})_{\Omega_{k}}\\\ \leq\frac{\mid\tau_{1}\mid T}{T_{0}-\rho
C_{\tau_{1}}}(\sum_{k=1}^{n_{el}}D_{B_{1k}})\|\textbf{R}_{1}^{n,\theta}\|\hskip
139.4185pt\\\
\sum_{k=1}^{n_{el}}\tau_{3}^{\prime}(R_{3}^{n,\theta},\bigtriangledown\cdot\tilde{\bigtriangledown}E^{A,n,\theta}_{c}+\textbf{u}_{h}\cdot\bigtriangledown
E^{A,n,\theta}_{c}-\alpha E^{A,n,\theta}_{c})_{k}\hskip 39.83385pt\\\
\leq\frac{\mid\tau_{3}\mid
T}{T_{0}-C_{\tau_{3}}}(\sum_{k=1}^{n_{el}}D_{B_{3k}})\|R_{3}^{n,\theta}\|\hskip
142.26378pt\\\ \end{split}$ (91)

and the other set of terms can be estimated as follows:

$\begin{split}\sum_{k=1}^{n_{el}}((1-\tau_{1}^{-1}\tau_{1}^{\prime})I_{d\times
d}\textbf{R}_{1}^{n,\theta},E^{A,n,\theta}_{\textbf{u}})_{\Omega_{k}}&\leq\frac{\rho\mid\tau_{1}\mid}{T_{0}-\rho
C_{\tau_{1}}}(\sum_{i=1}^{d}\sum_{k=1}^{n_{el}}B_{1k}^{i})\|\textbf{R}_{1}^{n,\theta}\|\\\
\sum_{k=1}^{n_{el}}(1-\tau_{3}^{-1}\tau_{3}^{\prime})(R_{3}^{n,\theta},E^{A,n,\theta}_{c})_{\Omega_{k}}&\leq\frac{\mid\tau_{3}\mid}{T_{0}-C_{\tau_{3}}}(\sum_{k=1}^{n_{el}}B_{2k})\|R_{3}^{n,\theta}\|\\\
\end{split}$ (92)

let us look into the form of the column vector d which has components
$\textbf{d}_{1},d_{2},d_{3}$.
d=
$\sum_{i=1}^{n+1}(\frac{1}{dt}M\tau_{k}^{\prime})^{i}(\textbf{F}-M\partial_{t}\textbf{U}_{h}-\mathcal{L}\textbf{U}_{h})=\sum_{i=1}^{n+1}(\frac{1}{dt}M\tau_{k}^{\prime})^{i}\textbf{R}$
Hence we have the components
$\textbf{d}_{1}=(\sum_{i=1}^{n+1}(\frac{1}{dt}\tau_{1}^{\prime})^{i})I_{d\times
d}\textbf{R}_{1}^{n,\theta}$, $d_{2}=0$ and
$d_{3}=(\sum_{i=1}^{n+1}(\frac{1}{dt}\tau_{3}^{\prime})^{i})R_{3}^{n,\theta}$
Now the terms containing the components of d can be estimated in the following
way:

$\sum_{k=1}^{n_{el}}(\tau_{1}^{\prime}I_{d\times
d}\textbf{d}_{1},\rho(\textbf{u}_{h}\cdot\bigtriangledown)E^{A,n,\theta}_{\textbf{u}}+\mu(c)\Delta
E^{A,n,\theta}_{\textbf{u}}+\bigtriangledown
E^{A,n,\theta}_{p})_{\Omega_{k}}\\\
=\sum_{k=1}^{n_{el}}(\tau_{1}^{\prime}\\{\sum_{i=1}^{n+1}(\frac{1}{dt}\tau_{1}^{\prime})^{i}\\}I_{d\times
d}\textbf{R}_{1}^{n,\theta},\rho(\textbf{u}_{h}\cdot\bigtriangledown)E^{A,n,\theta}_{\textbf{u}}+\mu(c)\Delta
E^{A,n,\theta}_{\textbf{u}}+\bigtriangledown
E^{A,n,\theta}_{p})_{\Omega_{k}}\\\
\leq\sum_{k=1}^{n_{el}}(\tau_{1}^{\prime}\\{\sum_{i=1}^{\infty}(\frac{1}{dt}\tau_{1}^{\prime})^{i}\\}I_{d\times
d}\textbf{R}_{1}^{n,\theta},\rho(\textbf{u}_{h}\cdot\bigtriangledown)E^{A,n,\theta}_{\textbf{u}}+\mu(c)\Delta
E^{A,n,\theta}_{\textbf{u}}+\bigtriangledown
E^{A,n,\theta}_{p})_{\Omega_{k}}\\\
=\sum_{k=1}^{n_{el}}(\frac{\rho\tau_{1}^{2}}{dt+\rho\tau_{1}}I_{d\times
d}\textbf{R}_{1}^{n,\theta},\rho(\textbf{u}_{h}\cdot\bigtriangledown)E^{A,n,\theta}_{\textbf{u}}+\mu(c)\Delta
E^{A,n,\theta}_{\textbf{u}}+\bigtriangledown
E^{A,n,\theta}_{p})_{\Omega_{k}}\hskip 22.76219pt\\\
\leq\frac{\mid\tau_{1}\mid}{T_{0}-\rho C_{\tau_{1}}}\rho
C_{\tau_{1}}(\sum_{k=1}^{n_{el}}D_{B_{1k}})\|\textbf{R}_{1}^{n,\theta}\|\hskip
150.79959pt$ (93)

Similarly the next few terms will follow the same way as above.

$\sum_{k=1}^{n_{el}}\tau_{3}^{\prime}(d_{3},\bigtriangledown\cdot\tilde{\bigtriangledown}E^{A,n,\theta}_{c}+\textbf{u}_{h}\cdot\bigtriangledown
E^{A,n,\theta}_{c}-\alpha E^{A,n,\theta}_{c})_{k}\hskip 99.58464pt\\\
\leq\frac{\mid\tau_{3}\mid}{T_{0}-C_{\tau_{3}}}C_{\tau_{3}}(\sum_{k=1}^{n_{el}}D_{B_{3k}})\|R_{3}^{n,\theta}\|\hskip
179.25235pt\\\ \sum_{k=1}^{n_{el}}(\tau_{1}^{-1}\tau_{1}^{\prime}I_{d\times
d}\textbf{d}_{1},E^{A,n,\theta}_{\textbf{u}})_{\Omega_{k}}\leq\frac{\rho\mid\tau_{1}\mid}{T_{0}-\rho
C_{\tau_{1}}}(\sum_{i=1}^{d}\sum_{k=1}^{n_{el}}B_{1k}^{i})\|\textbf{R}_{1}^{n,\theta}\|\hskip
56.9055pt\\\
\sum_{k=1}^{n_{el}}\tau_{3}^{-1}\tau_{3}^{\prime}(d_{3},E^{A,n,\theta}_{c})_{\Omega_{k}}\leq\frac{\mid\tau_{3}\mid}{T_{0}-C_{\tau_{3}}}(\sum_{k=1}^{n_{el}}B_{2k})\|R_{3}^{n,\theta}\|\hskip
54.06023pt$ (94)

The terms containing truncation errors already have been estimated earlier
during apriori error estimation. Now we estimate the remaining terms as
follows:

$\begin{split}&\sum_{k=1}^{n_{el}}\tau_{2}^{\prime}(R_{3}^{h,n,\theta},\bigtriangledown\cdot
E^{A,n,\theta}_{\textbf{u}})_{\Omega_{k}}\\\
&=\sum_{k=1}^{n_{el}}\tau_{2}^{\prime}(R_{3}^{h,n,\theta},\bigtriangledown\cdot
e^{n,\theta}_{\textbf{u}})_{\Omega_{k}}-\sum_{k=1}^{n_{el}}\tau_{2}^{\prime}(R_{3}^{h,n,\theta},\bigtriangledown\cdot
E^{I,n,\theta}_{\textbf{u}})_{\Omega_{k}}\\\
&=\sum_{k=1}^{n_{el}}\tau_{2}^{\prime}(\bigtriangledown\cdot
e^{n,\theta}_{\textbf{u}},\bigtriangledown\cdot
e^{n,\theta}_{\textbf{u}})_{\Omega_{k}}-\sum_{k=1}^{n_{el}}\tau_{2}^{\prime}(\bigtriangledown\cdot
e^{n,\theta}_{\textbf{u}},\bigtriangledown\cdot
E^{I,n,\theta}_{\textbf{u}})_{\Omega_{k}}\\\ \end{split}$ (95)
$\begin{split}&\leq\sum_{k=1}^{n_{el}}\tau_{2}\int_{\Omega_{k}}(\bigtriangledown\cdot
e^{n,\theta}_{\textbf{u}})^{2}+\sum_{k=1}^{n_{el}}\tau_{2}\int_{\Omega_{k}}\mid(\bigtriangledown\cdot
e^{n,\theta}_{\textbf{u}})(\bigtriangledown\cdot
E^{I,n,\theta}_{\textbf{u}})\mid\\\ &\leq
C_{\tau_{2}}(\sum_{i=1}^{d}\|\frac{\partial e^{n,\theta}_{ui}}{\partial
x_{i}}\|)^{2}+C_{\tau_{2}}(\sum_{i=1}^{d}\|\frac{\partial
e^{n,\theta}_{ui}}{\partial x_{i}}\|)(\sum_{i=1}^{d}\|\frac{\partial
E^{I,n,\theta}_{ui}}{\partial x_{i}}\|)\\\ &\leq
2C_{\tau_{2}}\sum_{i=1}^{d}\|\frac{\partial e^{n,\theta}_{ui}}{\partial
x_{i}}\|^{2}+\epsilon_{1}^{\prime}C_{\tau_{2}}\sum_{i=1}^{d}\|\frac{\partial
e^{n,\theta}_{ui}}{\partial
x_{i}}\|^{2}+\frac{C_{\tau_{2}}}{\epsilon_{1}^{\prime}}\sum_{i=1}^{d}\|\frac{\partial
E^{I,n,\theta}_{ui}}{\partial x_{i}}\|^{2}\\\ &\leq
C_{\tau_{2}}[(2+\epsilon_{1}^{\prime})\sum_{i=1}^{d}\|e^{n,\theta}_{ui}\|_{1}^{2}+\frac{h^{2}}{\epsilon_{1}^{\prime}}\sum_{i=1}^{d}(\frac{1+\theta}{2}\|u_{i}^{n+1}\|_{2}+\frac{1-\theta}{2}\|u_{i}^{n}\|_{2})^{2}]\\\
&\leq(2+\epsilon_{1}^{\prime})C_{\tau_{2}}\|e^{n,\theta}_{\textbf{u}}\|_{1}^{2}+h^{2}\frac{C_{\tau_{2}}}{\epsilon_{1}^{\prime}}\bar{C}_{5}\end{split}$
(96)

where the parameter $\bar{C}_{5}$ comes for applying assumption (iv). Now the
terms involving trancation error can be estimated in slightly different way as
we have done in the previous section. Let us present here a detailed
derivation of one term only and the other follows the same way.

$\begin{split}(TE_{2}^{n,\theta},E^{A,n,\theta}_{c})&=(TE_{2}^{n,\theta},e^{n,\theta}_{c})-(TE_{2}^{n,\theta},E^{I,n,\theta}_{c})\\\
&\leq\mid(TE_{2}^{n,\theta},e^{n,\theta}_{c})\mid+\mid(TE_{2}^{n,\theta},E^{I,n,\theta}_{c})\mid\\\
&\leq\|TE_{2}^{n,\theta}\|(\|e^{n,\theta}_{c}\|+\|E^{I,n,\theta}_{c}\|)\\\
&\leq\frac{1}{\epsilon_{2}^{\prime}}\|TE_{2}^{n,\theta}\|^{2}+\frac{\epsilon_{2}^{\prime}}{2}(\|e^{n,\theta}_{c}\|^{2}+\|E^{I,n,\theta}_{c}\|^{2})\\\
&\leq\frac{1}{\epsilon_{2}^{\prime}}\|TE_{2}^{n,\theta}\|^{2}+\frac{\epsilon_{2}^{\prime}}{2}\\{\|e^{n,\theta}_{c}\|^{2}+h^{4}(\frac{1+\theta}{2}\|c^{n+1}\|_{2}+\frac{1-\theta}{2}\|c^{n}\|_{2})^{2}\\}\\\
&\leq\frac{1}{\epsilon_{2}^{\prime}}\|TE_{2}^{n,\theta}\|^{2}+\frac{\epsilon_{2}^{\prime}}{2}\|e^{n,\theta}_{c}\|^{2}_{1}+h^{4}\frac{\epsilon_{2}^{\prime}}{2}\bar{C}_{6}\end{split}$
(97)

Similarly the estimated result for the remaining term is

$(\textbf{TE}_{1}^{n,\theta},E^{A,n,\theta}_{\textbf{u}})\leq\frac{1}{\epsilon_{2}^{\prime}}\|\textbf{TE}_{1}^{n,\theta}\|^{2}+\frac{\epsilon_{2}^{\prime}}{2}\|e^{n,\theta}_{\textbf{u}}\|^{2}_{1}+h^{4}\frac{\epsilon_{2}^{\prime}}{2}\bar{C}_{7}$
(98)

and this completes estimating all the terms of $RHS^{A}_{1}$. Now the
$trilinear$ term in $RHS$ in (96) can be estimated as follows using property
(b) of the $trilinear$ term as follows:

$\begin{split}c(e_{\textbf{u}}^{n},\textbf{u}^{n,\theta},e_{\textbf{u}}^{n,\theta})&\leq
C\|e_{\textbf{u}}^{n}\|\|\textbf{u}^{n,\theta}\|_{2}\|e_{\textbf{u}}^{n,\theta}\|_{1}\\\
&\leq\bar{C}_{8}\|e_{\textbf{u}}^{n}\|\|e_{\textbf{u}}^{n,\theta}\|_{1}\\\
&\leq\bar{C}_{8}\|e_{\textbf{u}}^{n,\theta}\|^{2}\end{split}$ (99)

The term $\|\textbf{u}^{n,\theta}\|_{2}$ is bounded by the virtue of
assumption (iv) and applying $Poincare$ inequality on last term in $RHS$ in
(96) we have

$\mu_{l}\|e_{\textbf{u}}^{n,\theta}\|^{2}\leq\mu_{l}C_{P}\mid
e_{\textbf{u}}^{n,\theta}\mid_{1}^{2}\leq\mu_{l}C_{P}\|e_{\textbf{u}}^{n,\theta}\|^{2}$
(100)

Now this completes finding bounds for each term in the $RHS$ of (93). Putting
common terms all together in the left hand side and multiplying them by $2$
and then integrating both sides over $(t^{n},t^{n+1})$ for $n=0,...,(N-1)$ ,
we will finally have

$\sum_{n=0}^{N-1}(\|e^{n+1}_{\textbf{u}}\|^{2}-\|e^{n}_{\textbf{u}}\|^{2})+\sum_{n=0}^{N-1}(\|e^{n+1}_{c}\|^{2}-\|e^{n}_{c}\|^{2})+\\{2\mu_{l}-2(2+\epsilon_{1}^{\prime})C_{\tau_{2}}\\\
\quad-2C_{8}-\mu_{l}C_{P}-\epsilon_{2}^{\prime}\\}\sum_{n=0}^{N-1}\int_{t^{n}}^{t^{n+1}}\|e^{n,\theta}_{\textbf{u}}\|_{1}^{2}dt+(2D_{\alpha}-\epsilon_{2}^{\prime})\sum_{n=0}^{N-1}\int_{t^{n}}^{t^{n+1}}\|e^{n,\theta}_{c}\|^{2}_{1}dt\\\
\leq
h^{2}\sum_{n=0}^{N-1}\int_{t^{n}}^{t^{n+1}}\\{\bar{C}_{1}\|\textbf{R}_{1}^{n,\theta}\|+\bar{C}_{2}\|R_{2}^{n,\theta}\|+\bar{C}_{3}\|R_{3}^{n,\theta}\|+\frac{2C_{\tau_{2}}}{\epsilon_{1}^{\prime}}\bar{C}_{5}+h^{2}\epsilon_{2}^{\prime}(\bar{C}_{6}+\bar{C}_{7})\\}dt\\\
\quad+2\frac{\mid\tau_{1}\mid}{T_{0}-\rho C_{\tau_{1}}}[(T+\rho
C_{\tau_{1}})(\sum_{k=1}^{n_{el}}D_{B_{1k}})+2\rho\sum_{i=1}^{d}\sum_{k=1}^{n_{el}}B_{1k}^{i}]\sum_{n=0}^{N-1}\int_{t^{n}}^{t^{n+1}}\|\textbf{R}_{1}^{n,\theta}\|^{2}dt\\\
+2\frac{\mid\tau_{3}\mid}{T_{0}-\rho
C_{\tau_{3}}}[(T+C_{\tau_{3}})\sum_{k=1}^{n_{el}}D_{B_{3k}}+2\sum_{k=1}^{n_{el}}B_{2k}]\sum_{n=0}^{N-1}\int_{t^{n}}^{t^{n+1}}\|R_{3}^{n,\theta}\|^{2}dt+\\\
\frac{2}{\epsilon_{2}^{\prime}}\sum_{n=0}^{N-1}\int_{t^{n}}^{t^{n+1}}(\|\textbf{TE}^{n,\theta}_{1}\|^{2}+\|TE^{n,\theta}_{2}\|^{2})dt\hskip
56.9055pt$ (101)

Choose the arbitrary parameters including $C_{\tau_{2}}$ and the $Poincare$
constant $C_{P}$ in such a way that all the coefficients in the left hand side
can be made positive. Then taking minimum over the coefficients in the left
hand side let us divide both sides by them. Using properties (15)-(16)
associated with both implicit time discretisation scheme and the fact that
$\tau_{1},\tau_{3}$ are of order $h^{2}$, we have arrived at the following
relation:

$\boxed{\|\textbf{u}-\textbf{u}_{h}\|_{\tilde{\textbf{V}}}^{2}+\|c-c_{h}\|_{\tilde{\textbf{V}}}^{2}\leq
C^{\prime}(\textbf{R})(h^{2}+dt^{2r})}$ (102)

where

$r=\begin{cases}1,&\text{if}\ \theta=1\hskip 2.84526ptfor\hskip
2.84526ptbackward\hskip 2.84526ptEuler\hskip 2.84526ptrule\\\ 2,&\text{if}\
\theta=0\hskip 2.84526ptfor\hskip 2.84526ptCrank-Nicolson\hskip
2.84526ptscheme\end{cases}$ (103)

This only completes one part of $aposteriori$ estimation and in the next part
we combine the corresponding pressure part.
Second part: Using the result (74) we can rewrite (73) in the following form:

$b(\textbf{v}_{h},I_{h}p-p_{h})=(\frac{\partial e_{\textbf{u}}}{\partial
t},\textbf{v}_{h})+c(e_{\textbf{u}},\textbf{u},\textbf{v}_{h})+c(\textbf{u}_{h},e_{\textbf{u}},\textbf{v}_{h})+a_{NS}(e_{\textbf{u}},\textbf{v}_{h})$
(104)

Integrating both sides with respect time

$\begin{split}\sum_{n=0}^{N-1}\int_{t^{n}}^{t^{n+1}}b(\textbf{v}_{h},E_{p}^{A,n,\theta})dt&=\sum_{n=0}^{N-1}\int_{t^{n}}^{t^{n+1}}\\{(\frac{e_{\textbf{u}}^{n+1}-e_{\textbf{u}}^{n}}{dt},\textbf{v}_{h})+c(e_{\textbf{u}}^{n},\textbf{u}^{n,\theta},\textbf{v}_{h})+\\\
&\quad
c(\textbf{u}_{h}^{n},e_{\textbf{u}}^{n,\theta},\textbf{v}_{h})+a_{NS}(e_{\textbf{u}}^{n,\theta},\textbf{v}_{h})\\}dt\\\
&=\sum_{n=0}^{N-1}\int_{t^{n}}^{t^{n+1}}\\{(\frac{e_{\textbf{u}}^{n+1}-e_{\textbf{u}}^{n}}{dt},\textbf{v}_{h})+c(e_{\textbf{u}}^{n},\textbf{u}^{n,\theta},\textbf{v}_{h})+\\\
&\quad
c(\textbf{u}^{n},\textbf{u}^{n,\theta},\textbf{v}_{h})-c(e_{\textbf{u}}^{n},e_{\textbf{u}}^{n,\theta},\textbf{v}_{h})+a_{NS}(e_{\textbf{u}}^{n,\theta},\textbf{v}_{h})\\}dt\end{split}$
(105)

Now applying $Cauchy-Schwarz$’s inequality, $Young$’s inequality, property (b)
of the $trilinear$ form $c(\cdot,\cdot,\cdot)$ and the above result (113) on
(116) we have

$\sum_{n=0}^{N-1}\int_{t^{n}}^{t^{n+1}}b(\textbf{v}_{h},E_{p}^{A,n,\theta})dt\leq\bar{C}^{\prime}(\textbf{R})(h^{2}+dt^{2r})\|\textbf{v}_{h}\|_{1}$
(106)

Applying this result on (75) we have

$\|I_{h}p-p_{h}\|^{2}_{L^{2}(L^{2})}\leq\bar{C}^{\prime\prime}(\textbf{R})(h^{2}+dt^{2r})$
(107)

Now combining the results obtained in the first and second part and applying
interpolation estimate on pressure interpolation term $E^{I}_{p}$, we finally
arrive at the following

$\boxed{\|\textbf{u}-\textbf{u}_{h}\|^{2}_{\tilde{\textbf{V}}}+\|p-p_{h}\|_{L^{2}(L^{2})}^{2}+\|c-c_{h}\|^{2}_{\tilde{\textbf{V}}}\leq\bar{C}(\textbf{R})(h^{2}+dt^{2r})}$
(108)

Now this finally completes derivation of $aposteriori$ error estimation. ∎

######  Remark 3.

These estimations clearly imply that the scheme is $first$ order convergent in
space with respect to total norm, whereas in time it is $first$ order
convergent for backward Euler time discretization scheme and $second$ order
convergent for Crank-Nicolson method.

##  5 Numerical Experiment

In this section we verify the credibility of $ASGS$ method for this coupled
transient $Navier$-$Stokes$-$VADR$ model through several numerical examples.
Here we present a comparative study between standard Galerkin and $ASGS$
finite element method. We have considered two broad cases based on one way
coupling and two-way or strong coupling. First case is further divided into
three sub-cases consisting of different values Reynolds number and in the
later one the viscosity of the fluid is taken to be dependent upon
concentration of the solute and variable diffusion coefficients have been
considered. This case too consists of two sub-cases involving different
viscosity coefficients.
Let us take $\Omega$ to be a square bounded domain (0,1) $\times$ (0,1).
Piecewise continuous linear finite element(P1) space is considered for
approximating velocity, pressure and concentration. Now renaming the error in
the following way we have examined the performances of both Galerkin and
$ASGS$ methods.

$Total$ $error$ =
$\\{\|\textbf{u}-\textbf{u}_{h}\|^{2}_{\tilde{\textbf{V}}}+\|p-p_{h}\|_{L^{2}(L^{2})}^{2}+\|c-c_{h}\|^{2}_{\tilde{\textbf{V}}}\\}^{\frac{1}{2}}$

The exact solutions for all the cases are taken as follows:
$\textbf{u}=(e^{-t}x^{2}(x-1)^{2}y(y-1)(2y-1),-e^{-t}x(x-1)(2x-1)y^{2}(y-1)^{2})$,
$p=e^{-t}(3x^{2}+3y^{2}-2)$ and $c=e^{-t}xy(x-1)(y-1)$

Time | Grid | Galerkin method | ASGS method
---|---|---|---
step | size | Total error | RoC | Total error | RoC
0.1 | 10 $\times$ 10 | 0.158556 |  | 0.158435 |
0.05 | 20 $\times$ 20 | 0.0833 | 0.928605 | 0.0833011 | 0.927481
0.025 | 40 $\times$ 40 | 0.0430609 | 0.95194 | 0.0430864 | 0.951103
0.0125 | 80 $\times$ 80 | 0.0219347 | 0.973161 | 0.0219556 | 0.972645
0.00625 | 160 $\times$ 160 | 0.0110526 | 0.98883 | 0.011068 | 0.988194
Table 1: Total error and Rate of convergence(RoC) under both Galerkin and $ASGS$ method for small Reynolds number(Re=50) at $T=1$ Time | Grid | Galerkin method | ASGS method
---|---|---|---
step | size | Total error | RoC | Total error | RoC
0.1 | 10 $\times$ 10 | 0.170253 |  | 0.158437 |
0.05 | 20 $\times$ 20 | 0.0871451 | 0.966187 | 0.0833212 | 0.92715
0.025 | 40 $\times$ 40 | 0.043821 | 0.991797 | 0.0431014 | 0.950949
0.0125 | 80 $\times$ 80 | 0.022057 | 0.990389 | 0.0219237 | 0.975243
0.00625 | 160 $\times$ 160 | 0.011189 | 0.979155 | 0.011076 | 0.985054
Table 2: Total error and Rate of convergence(RoC) under both Galerkin and $ASGS$ method for small Reynolds number(Re=500) at $T=1$ Time | Grid | Galerkin method | ASGS method
---|---|---|---
step | size | Total error | RoC | Total error | RoC
0.1 | 10 $\times$ 10 | 0.226209 |  | 0.158438 |
0.05 | 20 $\times$ 20 | 0.164603 | 0.458663 | 0.0833293 | 0.927026
0.025 | 40 $\times$ 40 | 0.0822173 | 1.00148 | 0.0431091 | 0.950832
0.0125 | 80 $\times$ 80 | 0.0310324 | 1.37098 | 0.0219345 | 0.974791
0.00625 | 160 $\times$ 160 | 0.022146 | 0.486729 | 0.0110245 | 0.992489
Table 3: Total error and Rate of convergence(RoC) under both Galerkin and $ASGS$ method for small Reynolds number(Re=5000) at $T=1$ Time | Grid | Galerkin method | ASGS method
---|---|---|---
step | size | Total error | RoC | Total error | RoC
0.1 | 10 $\times$ 10 | 0.159204 |  | 0.158826 |
0.05 | 20 $\times$ 20 | 0.0834547 | 0.931812 | 0.0834583 | 0.928321
0.025 | 40 $\times$ 40 | 0.0430917 | 0.953584 | 0.043141 | 0.951997
0.0125 | 80 $\times$ 80 | 0.021942 | 0.973715 | 0.0219748 | 0.973209
0.00625 | 160 $\times$ 160 | 0.011122 | 0.980278 | 0.011022 | 0.993703
Table 4: Total error and Rate of convergence(RoC) under both Galerkin and $ASGS$ method for variable viscosity and diffusion coefficients (first sub-case) at $T=1$ Time | Grid | Galerkin method | ASGS method
---|---|---|---
step | size | Total error | RoC | Total error | RoC
0.1 | 10 $\times$ 10 | 0.236613 |  | 0.161085 |
0.05 | 20 $\times$ 20 | 0.201906 | 0.22884 | 0.0855817 | 0.912444
0.025 | 40 $\times$ 40 | 0.128248 | 0.654755 | 0.0445594 | 0.941572
0.0125 | 80 $\times$ 80 | 0.0495898 | 1.37082 | 0.0227193 | 0.971811
0.00625 | 160 $\times$ 160 | 0.041146 | 0.269291 | 0.011148 | 1.027133
Table 5: Total error and Rate of convergence(RoC) under both Galerkin and
$ASGS$ method for variable viscosity and diffusion coefficients (second sub-
case) at $T=1$

(I)First case: Here we have considered constant viscosity coefficient and
therefore the coupled system becomes an one-way coupling. The importance
behind considering this case is here that we want to verify the performance of
$ASGS$ method for different Reynolds number. Here diffusion coefficients are
also taken constant.
(a)Small Reynolds number The exact solutions remain same. The values of
Reynolds number $Re$=50, diffusion coefficient $D$=2 and reaction coefficient
$\beta$=0.01.
Table 1 presents total errors and rates of convergence (RoC) of the coupled
system for this case under Galerkin and $ASGS$ methods for different time
steps $dt$ and grid sizes. It is clearly seen that both Galerkin and $ASGS$
method performs equally well for small Reynolds number. We can conclude the
order of convergence for each of the methods is 1.
(b) Medium Reynolds number: For this case the values of coefficients are taken
as $Re$=500, $D$=2 and $\beta$=0.01. Similar to the previous case table 2
represents the total errors and rates of convergence of the coupled system for
this case under Galerkin and $ASGS$ methods for different time steps $dt$ and
grid sizes. In this case though both Galerkin and $ASGS$ method perform
equally well and retain the desired first order convergence, but total error
obtained in $ASGS$ method is less compared to that of Galerkin method.
(c)Large Reynolds number The values of the coefficients are considered as
$Re$=5000, $D$=2 and $\beta$=0.01. Table 3 presents the total errors and rates
of convergence of the coupled system for this case under Galerkin and $ASGS$
methods. It can be observed that Galerkin method behaves in somewhat
oscillatory manner and it is not possible to conclude a definite order of
convergence for this case, whereas $ASGS$ method performs consistently well at
every time steps and grid sizes and rate of convergence in this case again
turns out to be 1.
(II)Second case: Here we consider the viscosity to be dependent upon
concentration and hence $Navier$-$Stokes$ and Transport equations are coupled
in two-way manner. The proposed expression of concentration dependent
viscosity is taken from [16] and depending upon different viscosity
coefficients we have divided this case into two sub-cases. In this case we
have considered variable diffusion coefficients as follows:
$D_{1}$= $e^{-t}y^{2}(y-1)^{2}(2y-1)^{2}x^{4}(x-1)^{4}$ and $D_{2}$=
$e^{-t}x^{2}(x-1)^{2}(2x-1)^{2}y^{4}(y-1)^{4}$
and the reaction coefficient $\beta$=0.01.
(a) First sub-case: The viscosity coefficient is $\mu(c)=0.00954e^{27.93\times
0.028c}$. Table 4 presents the total errors and rates of convergence of the
coupled system for this case under Galerkin and $ASGS$ methods for different
time steps $dt$ and grid sizes. Both Galerkin and $ASGS$ method performs
equally well and order of convergence for both the methods is 1.
(b) Second sub-case: Here we have considered slightly small viscosity
coefficient $\mu(c)=0.0000954e^{27.93\times 0.028c}$. Table 5 presents the
total errors and rates of convergence of the coupled system for this case
under Galerkin and $ASGS$ methods. This table shows that Galerkin method
performs poorly, whereas the $ASGS$ method performs far better and obtains the
desired first order convergence.

##  6 Conclusion

This paper presents algebraic $subgrid$ $multiscale$ stabilized finite element
analysis of transient $Navier$-$Stokes$ fluid flow equation strongly coupled
with unsteady $VADR$ transport problem. Consideration of concentration
dependent viscosity makes this time dependent coupling more accurate to model
real life based contemporary problems. To ensure the efficiency of the
stabilized finite element method for this model, both $apriori$ and
$aposteriori$ error estimates have been derived in detail. It is essential to
mention that the norm employed for error estimation consists of the full norms
corresponding to each variable belonging to their respective spaces. Therefore
it provides a wholesome information about convergence of the method.
Theoretically the rate of convergence for both $apriori$ and $aposteriori$
error estimations is $O(h)$ in space and first and second order convergences
have come out for two implicit time discretization schemes viz. backward Euler
and Crank-Nicolson methods respectively. The accuracy of the stabilized method
has been numerically tested through considering two different kind of examples
and various possible combinations among them. Numerical results both in
tabular and figure representations show better performance of the stabilized
$ASGS$ method than standard $Galerkin$ finite element method and verify
theoretically established results too.
Acknowledgement
This work has been supported by grant from Innovation in Science Pursuit for
Inspired Research (INSPIRE) programme sponsored and managed by the Department
of Science and Technology(DST), Ministry of Science and Technology, Govt.of
India.

## References

  * [1] A.N. Brooks, T.J.R. Hughes, Streamline upwind/Petrov–Galerkin formulations for convection dominated flows with particular emphasis on the incompressible Navier–Stokes equations, Computer Methods in Applied Mechanics and Engineering, 32: 199–259(1982).
  * [2] T.J.R. Hughes, L.P. Franca, G.M. Hulbert, A new finite element formulation for fluid dynamics: VIII. The Galerkin/least-squares method for advective–diffusive equations, Computer Methods in Applied Mechanics and Engineering, 73: 173–189(1989).
  * [3] T.J.R. Hughes, Multiscale phenomena: Green’s functions, the Dirichlet-to-Neumann formulation, subgrid scale models, bubbles and the origins of stabilized methods, Computer Methods in Applied Mechanics and Engineering, 127: 387–401(1995).
  * [4] S.K. Hannani, M. Stanislas, P. Dupont, Incompressible Navier-Stokes computations with SUPG and GLS formulations — A comparison study, Computer Methods in Applied Mechanics and Engineering, 124: 153-170(1995).
  * [5] R. Codina, H.Coppola-Owen, P. Nithiarasu, C.B. Liu, Numerical comparison of CBS and SGS as stabilization techniques for the incompressible Navier–Stokes equations, International Journal for Numerical Methods in Engineering, 66:1672–1689 (2006).
  * [6] B.S. Krik, G.F. Carey, Development and validation of a SUPG finite element scheme for the compressible Navier–Stokes equations using a modified inviscid flux discretization, International Journal for Numerical Methods in Fluids, 57:265-293(2008).
  * [7] A. R$\acute{u}$sso, Streamline-upwind Petrov/Galerkin method (SUPG) vs residual-free bubbles (RFB), Comput. Methods Appl. Mech. Engrg. 195: 1608–1620 (2006).
  * [8] A. R$\acute{u}$sso, Bubble stabilization of the finite element methods for the linearized incompressible Navier–Stokes equations. Computer Methods in Applied Mechanics and Engineering, 132: 335–343(1996).
  * [9] R. Codina, O.C. Zienkiewicz, CBS versus GLS stabilization of the incompressible Navier-Stokes equations and the role of the time step as stabilization parameter, Communications in Numerical Methods in Engineering, 18: 99-112(2002).
  * [10] D. Vassilev, I. Yotov, Coupling Stokes-Darcy flow with transport, SIAM J. Sci. Comput. 3661-3684(2009).
  * [11] A. Cesmelio$\breve{g}$lu, P. Chidyagwai, B. Rivi$\grave{e}$re, Continuous and discontinuous finite element methods for coupled surface-subsurface flow and transport problems.
  * [12] A. Cesmelio$\breve{g}$lu, B. Rivi$\grave{e}$re, Existence of a weak solution for the fully coupled Navier-Stokes/Darcy-transport problem, J. Differential Equations, 252, 4138-4175(2012).
  * [13] H. Rui, J. Zhang, A stabilized mixed finite element method for coupled Stokes and Darcy flows with transport, Comput. Methods Appl. Mech. Engrg. 315, 169-189(2017).
  * [14] M. Chowdhury, B.V. Rathish Kumar, Apriori and aposteriori error estimation of Subgrid multiscale stabilized finite element method for coupled unified Stokes-Brinkman/Transport model, arxiv (pre-print), math.AP, 2004.01782.
  * [15] R. Codina, Comparison of some finite element methods for solving the diffusion-convection-reaction equation, Comput. Methods Appl. Mech. Engrg. 156: 185-210 (1998).
  * [16] J. Chirife, M. P. Buera, A Simple Model for Predicting the Viscosity of Sugar and Oligosaccharide Solutions, Journal of Food Engineering 33, 221-236(1997).
  * [17] B. Rivi$\grave{e}$re, M. F. Wheeler, A Discontinuous Galerkin Method Applied to Nonlinear Parabolic Equations, Discontinuous Galerkin Methods, Springer, pp. 231-244(2000).
  * [18] J. Blasco, R. Codina, Space and time error estimates for a first order pressure stabilized finite element method for the incompressible Navier-Stokes equations, Applied Numerical Mathematics 28, 475-497(2001).
  * [19] M. Chowdhury, B.V.R. Kumar, A priori and a posteriori error estimation for finite element approximation of advection-diffusion-reaction equation with spatially variable coefficients, arxiv (pre-print), math.AP, 1811.05283.
  * [20] M. Chowdhury, B.V.R. Kumar, On subgrid multiscale stabilized finite element method for advection-diffusion-reaction equation with variable coefficients, Applied Numerical Mathematics 150, 576-586(2020).
  * [21] G.R. Barrenechea, E. Castillo, R. Codina, Time-dependent semi-discrete analysis of the viscoelastic fluid flow problem using a variational multiscale stabilised formulation, IMA Journal of Numerical Analysis, 1-25(2018).
  * [22] S. Badia, R. Codina, Unified stabilized finite element formulations for the Stokes and Darcy problems, SIAM J. Numer. Anal., 47(3), 1971–2000(2009).
  * [23] R. Du, Z. Liu, A lattice Boltzmann model for the fractional advection–diffusion equation coupled with incompressible Navier–Stokes equation, Applied Mathematics Letters, 101(2020).
  * [24] X. Yua, K. Regenauer-Lieba, F. Tian, A hybrid immersed boundary-lattice Boltzmann/finite difference method for coupled dynamics of fluid flow, advection, diffusion and adsorption in fractured and porous media, Computers and Geosciences, 128, 70-78(2019).
  * [25] G. Hauke, A simple subgrid scale stabilized method for the advection–diffusion-reaction equation, Comput. Methods Appl. Mech. Engrg., 191, 2925–2947(2002).
  * [26] R. Codina, A stabilized finite element method for generalized stationery incompressible flows, Comput. Methods Appl. Mech. Engrg., 190, 2681-2706(2001).
  * [27] R. Codina, J. Blasco, Stabilized finite element method for the transient Navier-Stokes equations based on a pressure gradient projection, Comput. Methods Appl. Mech. Engrg., 182, 277-300(2000).
  * [28] A. Rasam, S. Wallin, G. Brethouwer , A. V. Johansson, Large eddy simulation of channel flow with and without periodic constrictions using the explicit algebraic subgrid-scale model, Journal of Turbulence, 15(11), 752–775(2014).
  * [29] R. Codina, On stabilized finite element methods for linear systems of convection-diffusion-reaction equations, Comput. Methods Appl. Mech. Engrg., 188, 61-82 (2000).
  * [30] M. Bischoff, K. Bletzinger, Improving stability and accuracy of Reissner–Mindlin plate finite elements via algebraic subgrid scale stabilization, Comput. Methods Appl. Mech. Engrg., 193, 1517–1528 (2004).
  * [31] O. Guasch, R. Codina, An algebraic subgrid scale finite element method for the convected Helmholtz equation in two dimensions with applications in aeroacoustics, Comput. Methods Appl. Mech. Engrg., 196, 4672–4689(2007).
  * [32] E. Castillo, R. Codina, Dynamic term-by-term stabilized finite element formulation using orthogonal subgrid-scales for the incompressible Navier-Stokes problem, Comput. Methods Appl. Mech. Engrg., 349, 701-721(2019).
  * [33] R. Codina, J. Blasco, Analysis of a stabilized finite element approximation of the transient convection-diffusion-reaction equation using orthogonal subscales, Comput Visual Sci, 4, 167–174 (2002).
  * [34] J. Baigesa, R. Codina, Variational Multiscale error estimators for solid mechanics adaptive simulations: An Orthogonal Subgrid Scale approach, Comput. Methods Appl. Mech. Engrg., 325, 37-55(2017).
  * [35] R. Codina, Analysis of a stabilized finite element approximation of the Oseen equations using orthogonal subscales, Applied Numerical Mathematics, 58, 264-283 (2008).
  * [36] C. Bayona, J. Baiges, R. Codina, Variational multiscale approximation of the one-dimensional forced Burgers equation: The role of orthogonal subgrid scales in turbulence modeling, Int J Numer Meth Fluids.,86, 313–328(2018).
  * [37] V. Girault, P.A. Raviart, Finite Element Methods for Navier-Stokes Equations: Theory and Algorithms, vol. 5 of Springer series in computational mathematics. Springer, Berlin (1986).
  * [38] E. Burman, M. A. Fernández, Continuous interior penalty finite element method for the time-dependent Navier–Stokes equations: space discretization and convergence, Numer. Math., 107,39–77(2007).
  * [39] J. A. Wheeler, M. F. Wheeler, I. Yotovc, Enhanced velocity mixed finite element methods for flow in multiblock domains, Computational Geosciences 6: 315–332(2002).