# Apriori and aposteriori error estimation of Subgrid multiscale stabilized
finite element method for coupled unified Stokes-Brinkman/Transport model

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

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

###### Abstract

In this study, we present a stabilized finite element analysis for completely
unified Stokes-Brinkman problems fully coupled with variable coefficient
transient Advection-Diffusion-Reaction equation($VADR$). As well we have
carried out the stabilized finite element analysis for Stokes-Brinkman model
with interface conditions fully coupled with $VADR$. The viscosity of the
fluid, involved in flow problem, depends on the concentration of the solute,
whose transport is described by $VADR$ equation. The algebraic subgrid
multiscale approach has been employed to arrive at the stabilized coupled
variational formulation. For the time discretization the fully implicit Euler
scheme has been used. A detailed derivation of both the apriori and
aposteriori estimates for the stabilized subgrid multiscale finite element
scheme have been presented. Few numerical experiments have been carried out to
verify the credibility of the method.

Keywords Stokes-Brinkman equation $\cdot$ Advection-diffusion-reaction
equation $\cdot$ Subgrid multiscale stabilized method $\cdot$ A priori error
estimation $\cdot$ A posteriori error estimation

##  1 Introduction

Solute transport model coupled with fluid flow equations in porous media
region plays a crucial role in modelling physical phenomena of physiological
and environmental importance. For several decades many studies have focused on
transport model and fluid flow equations separately. For instance [22]-[28]
have studied various numerical methods for solving advection-diffusion-
reaction equation where as [14]-[18] have focused on studying Stokes-Darcy
flow equations and [28]-[33] and [7]-[13] have respectively worked with
Navier-Stokes fluid flow model and Brinkman model as well. There are many
other researches done in this regard. Here we have mentioned only few of them.
Contemporary complex problems, such as surface water and ground water
pollution due to various human activities, stenosis in human arteries etc,
have made it inevitable to model the problems in more effective way and hence
study of coupled fluid flow-transport equations is attracting great attention
in today’s research. [1] presents finite element analysis of Stokes-Darcy flow
coupled with transport equation in which the fluid flow equations are solved
through mixed finite element method and transport equation for solute
concentration is solved using local discontinuous Galerkin ($LDG$) method. As
well another study [3] discusses continuous and discontinuous finite element
methods for coupled Navier-Stokes/Darcy-transport problem. Though these are
one way or weak coupling in the sense that velocity field obtained after
solving flow problem becomes an input data in transport model. Cesmelioglu and
Rivière in [2] has introduced a two-way strong coupling of Navier-Stokes/Darcy
with transport equation through considering the viscosity of the fluid
depending upon concentration of the solute. The strong coupling is proven to
be more accurate in modelling problems. In [2] authors prove the existence and
uniqueness of the weak solution of the variational form of Navier-
Stokes/Darcy-transport model. Further under constrained viscosity
consideration in [6] authors derive apriori error estimates for a stabilized
mixed finite element scheme for the strongly coupled Stokes-Darcy-ADR model
though they present numerical results to a one-way coupling problem.
In this paper we present Stokes-Brinkman in a completely unified manner
strongly coupled with transport equation. The Brinkman equation is an
extension to Darcy’s law when boundary layer regions have to be taken into
account. The Stokes-Brinkman system has an important role in modelling highly
heterogeneous porous media problem. As per our knowledge adequate attention is
not paid to study numerical method for solving this model. Only few works
[36]-[39] are available in the literature dealing with Stokes-Brinkman model.
Here we have derived algebraic subgrid scale ($ASGS$) stabilized finite
element method for the coupled problem using continuous velocities and
pressure spaces across the inter-element boundaries. In [6] Rui and Zhang have
studied a mixed stabilized finite element method for coupled Stokes-
Darcy/transport problem, but $ASGS$ approach to study coupled Stokes-
Brinkman/$VADR$ model is very new. Hughes in [34] has introduced the concept
of stabilized multiscale subgrid method for Helmholtz equation and further
developments are going on afterwards. In [35] Codina presents a study on
comparison of stabilized finite element methods viz. $SUPG$, $GLS$, $SGS$,
$Taylor-Galerkin$ etc. for solving diffusion-convection-reaction equation and
experimentally shows that $SGS$ performs well in compared to other stabilized
method.
In particular the $ASGS$ approach consists of algebraic approximation of the
subscales that arise from the decomposition of the exact solution field into
resolvable scale and unresolvable scale, have been used for finite element
scheme stabilization. Stabilization parameters are derived following the
approach in [5], [7] for $ASGS$ method. Apriori error estimates for the
current stabilized ASGS finite element method for the unified strongly coupled
Stokes-Brinkman/$VADR$ have been derived. Further the aposteriori error
estimates following the residual approach have also been carried out. Further
more this paper has considered a porous media flow governed by Stokes-Brinkman
with interface condition strongly coupled with transport equation and
presented the corresponding stabilized formulation of the problem. The
interface conditions are taken care of by the standard continuity conditions
of normal velocities, normal stresses and concentration at the interface and
the Beavers-Joseph-Saffman condition at the interface allows Stokes fluid to
slip in contact with porous medium. Numerical studies have shown the
realization of theoretical order of convergence and the robustness of current
stabilized $ASGS$ finite element method for Stokes-Brinkman-$VADR$ tightly
coupled system.
Organisation of the paper is as follows: Section 2 starts from introducing the
model and finishes at Subgrid formulation going through two more subsections
presenting weak formulation and semi-discrete formulation. Next section has
elaborately described the derivation of apriori and aposteriori error
estimations for this subgrid formulation. Section 4 presents stabilization
formulation of coupled Stokes-Brinkman/transport equations through interfaces.
At last section 5 contains numerical results to verify the numerical
performance of the method.

##  2 Model problem

Let $\Omega\subset R^{d}$, d=2,3 be an open bounded domain with piecewise
smooth boundary $\partial\Omega$. For the sake of simplicity in further
calculations, we have considered two dimensional model, but it can be easily
extended for three dimensional model. Let us first mention the Stokes-Brinkman
flow problem for a fluid as follows:
Find u: $\Omega$ $\times$ (0,T) $\rightarrow R^{2}$ and $p$: $\Omega\times$
(0,T) $\rightarrow R$ such that,

$\begin{split}-\mu(c)\Delta\textbf{u}+\sigma\textbf{u}+\bigtriangledown
p&=\textbf{f}_{1}\hskip 5.69054ptin\hskip 5.69054pt\Omega\times(0,T)\\\
\bigtriangledown\cdot\textbf{u}&=f_{2}\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 u= ($u_{1},u_{2}$) is the velocity of the fluid or solvent, p is the
pressure, $\mu(c)$ is the viscosity of the fluid depending on concentration
$c$ of the dispersing mass of the solute, $\sigma$ is the inverse of
permeability, $\textbf{f}_{1}$ is the body force, $f_{2}$ is source term and
$\textbf{u}_{0}$ is the initial velocity. When $\sigma=0$ the flow problem is
fully Stokes. For Stokes flow $f_{2}=0$ too.
This Stokes-Brinkman flow problem is fully-coupled with the following ADR
equation with variable coefficients($VADR$), which represents the
transportation of solute in the same domain $\Omega$.
Find $c$: $\Omega\times$ (0,T) $\rightarrow R$ such that,

$\begin{split}\phi\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)\\\
\tilde{\bigtriangledown}c\cdot\textbf{n}&=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}:=(D_{1}\frac{\partial}{\partial
x},D_{2}\frac{\partial}{\partial y})$
$\phi$ is the porosity,$D_{1},D_{2}$ are variable diffusion coefficients,
$\alpha$ is the reaction coefficient and $g$ denotes the source of solute
mass, n is the outward normal to $\partial\Omega$ and $c_{0}$ is the initial
concentration of the solute. For purely Stokes flow problem $\phi$ takes value
1.
Letting U= (u,p,c) the equations all together can be written in the following
operator form,

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

where M, a matrix = diag(0,0,0,$\phi$),
$\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}=\begin{bmatrix}-\mu(c)\Delta\textbf{u}+\sigma\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}_{1}\\\ f_{2}\\\ g\end{bmatrix}$

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

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

Now we impose suitable assumptions, that are necessary to conclude the results
further, on the coefficients mentioned above.

(i) The fluid viscosity $\mu(c)=\mu\in C^{0}(R^{+};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 R^{+}$ (4)

(ii) $D_{1}=D_{1}((x,y),t)\in C^{0}(R^{2}\times(0,T);R)$ and
$D_{2}=D_{2}((x,y),t)\in C^{0}(R^{2}\times(0,T);R)$ where
$C^{0}(R^{2}\times(0,T);R)$ is the space of real valued continuous function
defined on $R^{2}$ for fixed $t\in(0,T)$. Both are bounded quantity that is we
can find lower and upper bounds for both of them.

(iii) $\sigma$ and $\alpha$ are positive constants.

(iv) The spaces of continuous solution $(\textbf{u},p,c)$ are assumed as:
$u_{1},u_{2}\in L^{\infty}(0,T;H^{2}(\Omega))\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) One additional assumption on continuous velocity solution is that $u_{1}$
and $u_{2}$ are taken to be bounded functions on $\Omega$.

###  2.1 Weak formulation

Let us first introduce the spaces as follows,
$H^{1}(\Omega)=\\{v\in L^{2}(\Omega):\bigtriangledown v\in L^{2}(\Omega)\\}$
Let $V_{s}=H^{1}_{0}(\Omega)=\\{v\in H^{1}(\Omega):v=0\hskip 2.84526pton\hskip
2.84526pt\partial\Omega\\}$ and $Q_{s}=L^{2}(\Omega)$ and J= (0,T)
Let $\tilde{\textbf{V}}$ := $L^{2}(0,T;V_{s})\bigcap L^{\infty}(0,T;Q_{s})$
Let us introduce another notation $\textbf{V}_{F}=V_{s}\times V_{s}\times
Q_{s}\times V_{s}$
The weak formulation of (1)-(2) is to find U= (u,p,c): J
$\rightarrow\textbf{V}_{F}$ such that $\forall$ V=(v,q,d) $\in\textbf{V}_{F}$

$(\frac{\partial c}{\partial
t},d)+a_{S}(\textbf{u},\textbf{v})-b(\textbf{v},p)+b(\textbf{u},q)+a_{T}(c,d)=l_{S}^{1}(\textbf{v})+l_{S}^{2}(q)+l_{T}(d)$
(5)

where
$a_{S}(\textbf{u},\textbf{v})=\int_{\Omega}\mu(c)\bigtriangledown\textbf{u}:\bigtriangledown\textbf{v}+\sigma\int_{\Omega}\textbf{u}\cdot\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_{S}^{1}(\textbf{v})=\int_{\Omega}\textbf{f}_{1}\cdot\textbf{v}$,
$l_{S}^{2}(q)=\int_{\Omega}f_{2}q$ and $l_{T}(d)=\int_{\Omega}gd$
Again the above formulation can be written as,
Find $\textbf{U}\in\textbf{V}_{F}$ such that

$(M\partial_{t}\textbf{U},\textbf{V})+B(\textbf{U},\textbf{V})=L(\textbf{V})\hskip
5.69054pt\forall\textbf{V}\in\textbf{V}_{F}$ (6)

where
$B(\textbf{U},\textbf{V})=a_{S}(\textbf{u},\textbf{v})-b(\textbf{v},p)+b(\textbf{u},q)+a_{T}(c,d)$
$L(\textbf{V})=l_{S}^{1}(\textbf{v})+l_{S}^{2}(q)+l_{T}(d)$

######  Remark 1.

[7] discusses about well posedness of unified Stokes-Darcy equation for
positive viscosity coefficient.

######  Remark 2.

The existence of the weak solution of the variational form for coupled Stokes-
Darcy/transport equation has been discussed in [2]. Under the assumptions
[(i)-(iv)] the existence of unique weak solution of the variational form (6)
can be established easily following the approach presented in [2], as this
model contains only linear terms.

###  2.2 Semi-discrete formulation

In this section we will introduce the standard Galerkin finite element space
discretization for the above variational form (6).
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_{s}^{h}$ and $Q_{s}^{h}$ be finite dimensional subspaces of $V_{s}$ and
$Q_{s}$ respectively. They are taken as follows,
$V_{s}^{h}=\\{v\in V_{s}:v(\Omega_{k})=\mathcal{P}^{2}(\Omega_{k})\\}$
$Q_{s}^{h}=\\{q\in Q_{s}:q(\Omega_{k})=\mathcal{P}^{1}(\Omega_{k})\\}$
where $\mathcal{P}^{1}(\Omega_{k})$ and $\mathcal{P}^{2}(\Omega_{k})$ denote
complete polynomial of order 1 and 2 respectively over each $\Omega_{k}$ for
k=1,2,…,$n_{el}$.
Let us consider similar notation $\textbf{V}_{F}^{h}$ for corresponding finite
dimensional subspace of $\textbf{V}_{F}$ where
$\textbf{V}_{F}^{h}=V_{s}^{h}\times V_{s}^{h}\times Q_{s}^{h}\times V_{s}^{h}$
Now the Galerkin formulation of the variational form (6) will be as follows:
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{V}_{h})=L(\textbf{V}_{h})$
(7)

where $(M\partial_{t}\textbf{U}_{h},\textbf{V}_{h})$= $(\frac{\partial
c_{h}}{\partial t},d_{h})$
$B(\textbf{U}_{h},\textbf{V}_{h})=a_{S}(\textbf{u}_{h},\textbf{v}_{h})-b(\textbf{v}_{h},p_{h})+b(\textbf{u}_{h},q_{h})+a_{T}(c_{h},d_{h})$
$L(\textbf{V}_{h})=l_{S}^{1}(\textbf{v}_{h})+l_{S}^{2}(q_{h})+l_{T}(d_{h})$

###  2.3 Subgrid multiscale formulation

This stabilization method has been introduced to correct the lack of stability
that the Galerkin method suffers due to small diffusion coefficient. It
involves decomposition of the solution space $\textbf{V}_{F}$ into the spaces
of resolved scales and unresolved scales. The finite element space
$\textbf{V}_{F}^{h}$ is considered as the space of resolved scales. Then the
final form of subgrid formulation will be arrived while the elements of
unresolved scales will be expressed in the terms of elements of resolved
scales.
Following the procedure described in [4] the variational subgrid scale model
for this coupled equation will be written as follows,
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{V}_{h})=L_{ASGS}(\textbf{V}_{h})$
(8)

where
$B_{ASGS}(\textbf{U}_{h},\textbf{V}_{h})=B(\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{d}),-\mathcal{L}^{*}\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{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{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&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}\tau_{1k}I&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})\hskip
5.69054pt(say)$

I is an identity matrix.
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})$
considering $d_{i}$ for i=1,2,3,4 are components of the matrix d and it can be
easily observed that $d_{1},d_{2},d_{3}$ are always 0 because of the matrix M.
We have the forms of the stabilization parameters $\tau_{1k},\tau_{2k}$ for
unified Stokes-Darcy problem in [7] and $\tau_{3k}$ for ADR equation with
variable coefficients in [5] and for each k all the coefficients $\tau_{ik}$
coincide with $\tau_{i}$ for i=1,2,3 that is, for each k=1,2,…,$n_{el}$

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

where $c_{1}^{\textbf{u}},c_{2}^{\textbf{u}},c_{1}^{p}$ are the suitable
parameters and $h$ is the mesh size.

######  Remark 3.

Since we are working with continuous velocities and pressure at the inter-
element boundaries, therefore we will not have any jump term in the above
stabilized formulation.

##  3 Error estimates

We start this section with the introduction of the notion of error terms,
followed by splitting of those error terms through introducing the projection
operator corresponding to each unknown variable. Later we have introduced
fully-discrete formulation and then conducted apriori and aposteriori error
estimates.

###  3.1 Projection operators : Error splitting

Let $\textbf{e}=(e_{\textbf{u}},e_{p},e_{c})$ denote the error where the
components are
$e_{\textbf{u}}=(e_{u1},e_{u2})=(u_{1}-u_{1h},u_{2}-u_{2h}),e_{p}=(p-p_{h})$
and $e_{c}=(c-c_{h})$. Here all the remaining notations carry their respective
meanings.
Let us introduce the projection operator for each of these error components.
(i)For any $\textbf{u}\in H^{2}(\Omega)\times H^{2}(\Omega)$ we assume that
there exists an interpolation $I^{h}_{\textbf{u}}:H^{2}(\Omega)\times
H^{2}(\Omega)\longrightarrow V_{s}^{h}\times V_{s}^{h}$ satisfying
(a) $b(\textbf{u}-I^{h}_{\textbf{u}}\textbf{u},q_{h})=0$ $\forall q_{h}\in
Q_{s}^{h}$ and
each component of the projection map that is
$I^{h}_{u_{1}}:H^{2}(\Omega)\longrightarrow V_{s}^{h}$ and
$I^{h}_{u_{2}}:H^{2}(\Omega)\longrightarrow V_{s}^{h}$ are $L^{2}$ orthogonal
projections, satisfying
(b) for any $u_{1}\in H^{2}(\Omega)$ $(u_{1}-I^{h}_{u_{1}}u_{1},v_{1h})=0$
$\forall v_{1h}\in V_{s}^{h}$ and
(c)for any $u_{2}\in H^{2}(\Omega)$ $(u_{2}-I^{h}_{u_{2}}u_{2},v_{2h})=0$
$\forall v_{2h}\in V_{s}^{h}$
(ii) Let $I^{h}_{p}:H^{1}(\Omega)\longrightarrow Q_{s}^{h}$ be the $L^{2}$
orthogonal projection given by
$\int_{\Omega}(p-I^{h}_{p}p)q_{h}=0$ $\forall q_{h}\in Q_{s}^{h}$ and for any
$p\in H^{1}(\Omega)$
(iii)Let $I^{h}_{c}:H^{2}(\Omega)\longrightarrow V_{s}^{h}$ be the $L^{2}$
orthogonal projection given by
$\int_{\Omega}(c-I^{h}_{c}c)d_{h}=0$ $\forall d_{h}\in V_{s}^{h}$ and for any
$c\in H^{2}(\Omega)$
Now each component of the error can be split into two parts interpolation
part, $E^{I}$ and auxiliary part, $E^{A}$ as follows:
$e_{u1}=(u_{1}-u_{1h})=(u_{1}-I^{h}_{u1}u_{1})+(I^{h}_{u1}u_{1}-u_{1h})=E^{I}_{u1}+E^{A}_{u1}$
Similarly, $e_{u2}=E^{I}_{u2}+E^{A}_{u2}$, $e_{p}=E^{I}_{p}+E^{A}_{p}$, and
$e_{c}=E^{I}_{c}+E^{A}_{c}$
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}_{c},d_{h})=0\hskip 5.69054ptd_{h}\in
V_{s}^{h}$ (10)

Proof: We have $(c-I^{h}_{c}c,d_{h})=0=(E^{I}_{c},d_{h})$ $\forall d_{h}\in
V_{s}^{h}$
Therefore

$\begin{split}\frac{d}{dt}(E^{I}_{c},d_{h})&=0\hskip 2.84526pt\forall d_{h}\in
V_{s}^{h}\\\ (\frac{\partial}{\partial
t}E^{I}_{c},d_{h})+(E^{I}_{c},\frac{\partial}{\partial t}d_{h})&=0\hskip
2.84526pt\forall d_{h}\in V_{s}^{h}\\\ (\frac{\partial}{\partial
t}E^{I}_{c},d_{h})&=0\hskip 2.84526pt\forall d_{h}\in V_{s}^{h}\\\
\end{split}$ (11)

Since $\frac{\partial}{\partial t}d_{h}\in V_{s}^{h}$, the second term in
second equation $(E^{I}_{c},\frac{\partial}{\partial t}d_{h})=0$
Useful interpolation estimation results [19] are as follows: 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}$
(12)

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.

###  3.2 Fully-discrete form

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}$ (13)

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}$ (14)

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})$ (15)

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 [40],we have

$\begin{split}\frac{c^{n+1}-c^{n}}{dt}&=c_{t}(\textbf{x},t^{n,\theta})+TE\mid_{t=t^{n,\theta}}\hskip
14.22636pt\forall\textbf{x}\in\Omega\end{split}$ (16)

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}\|TE^{n,\theta}\|&\leq\begin{cases}C^{\prime}dt\|c_{tt}^{n,\theta}\|_{L^{\infty}(t^{n},t^{n+1},L^{2})}&if\hskip
2.84526pt\theta=1\\\
C^{\prime\prime}dt^{2}\|c_{ttt}^{n,\theta}\|_{L^{\infty}(t^{n},t^{n+1},L^{2})}&if\hskip
2.84526pt\theta=0\end{cases}\end{split}$ (17)

Now for backward Euler scheme ($\theta$=1) applying assumption (iv) we will
have another property as follows:

$\begin{split}\|TE^{n,\theta}\|&\leq
C^{\prime}dt\|c_{tt}^{n+1}\|_{L^{\infty}(t^{n},t^{n+1},L^{2})}\\\
&\leq\tilde{C}dt\end{split}$ (18)

After introducing all the required definitions finally the fully-discrete
formulation of sub-grid 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,\theta},\textbf{V}_{h})=L_{ASGS}(\textbf{V}_{h})+(TE^{n,\theta},d_{h})$
(19)

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,\theta},\textbf{V}_{h})=L(\textbf{V}_{h})+(TE^{n,\theta},d_{h})$
(20)

###  3.3 Apriori error estimation

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
mention few definitions of norm in which we are going to estimate the errors:

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

######  Theorem 1.

(Auxiliary error estimate) For velocity $\textbf{u}_{h}=(u_{1h},u_{2h})$,
pressure $p_{h}$ and concentration $c_{h}$ belonging to $V_{s}^{h}\times
V_{s}^{h}\times Q_{s}^{h}\times V_{s}^{h}$ satisfying (10), 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 u,p,c ,
such that

$\|E^{A}_{u1}\|^{2}_{L^{2}(H^{1})}+\|E^{A}_{u2}\|^{2}_{L^{2}(H^{1})}+\|E^{A}_{p}\|_{L^{2}(L^{2})}^{2}+\|E^{A}_{c}\|^{2}_{\tilde{\textbf{V}}}\leq
C(h^{2}+h+dt^{2r})$ (22)

where

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

###### 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 (19) from (20) and then simplifying the terms, we have
$\forall\hskip 2.84526pt\textbf{V}_{h}\in V_{s}^{h}\times V_{s}^{h}\times
Q_{s}^{h}\times V_{s}^{h}$

$(M\frac{(\textbf{U}^{n+1}-\textbf{U}^{n+1}_{h})-(\textbf{U}^{n}-\textbf{U}^{n}_{h})}{dt},V_{h})+B(\textbf{U}^{n,\theta}-\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,\theta}-\textbf{U}_{h}^{n,\theta}),-\mathcal{L}^{*}\textbf{V}_{h})_{\Omega_{k}}-\sum_{k=1}^{n_{el}}\tau_{k}^{\prime}(\textbf{d},-\mathcal{L}^{*}\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,\theta}-\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}}\\\
=(TE^{n,\theta},d_{h})$ (24)

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}_{h}^{n,\theta}))$
Let us divide the big expressions into small parts, then using error splitting
in each of them and simplifying further, we will have them as follows:
Let

$\begin{split}I_{1}&=(M\frac{(\textbf{U}^{n+1}-\textbf{U}^{n+1}_{h})-(\textbf{U}^{n}-\textbf{U}^{n}_{h})}{dt},V_{h})\\\
&=(\frac{(c^{n+1}-c^{n+1}_{h})-(c^{n}-c^{n}_{h})}{dt},d_{h})\\\
&=(\frac{(E^{I,n+1}_{c}+E^{A,n+1}_{c})-(E^{I,n}_{c}+E^{A,n}_{c})}{dt},d_{h})\\\
&=(\frac{E^{A,n+1}_{c}-E^{A,n}_{c}}{dt},d_{h})\end{split}$ (25)

We arrive at the last line after using result 1, deduced in the previous
section.

$\begin{split}I_{2}&=B(\textbf{U}^{n,\theta}-\textbf{U}^{n,\theta}_{h},\textbf{V}_{h})\\\
&=\int_{\Omega}\mu(c^{n})\bigtriangledown(\textbf{u}^{n,\theta}-\textbf{u}^{n,\theta}_{h}):\bigtriangledown\textbf{v}_{h}+\int_{\Omega}\sigma(u_{1}^{n,\theta}-u^{n,\theta}_{1h})v_{1h}\\\
&\quad+\int_{\Omega}\sigma(u_{1}^{n,\theta}-u^{n,\theta}_{1h})v_{2h}-\int_{\Omega}(\bigtriangledown\cdot\textbf{v}_{h})(p^{n,\theta}-p_{h}^{n,\theta})+\int_{\Omega}(\bigtriangledown\cdot\textbf{u}^{n,\theta}-\textbf{u}_{h}^{n,\theta})q_{h}\\\
&\quad+\int_{\Omega}\tilde{\bigtriangledown}(c^{n,\theta}-c^{n,\theta}_{h})\cdot\bigtriangledown
d_{h}+\int_{\Omega}d_{h}\textbf{u}^{n}\cdot\bigtriangledown(c^{n,\theta}-c^{n,\theta}_{h})+\int_{\Omega}\alpha(c^{n,\theta}-c^{n,\theta}_{h})d_{h}\\\
&=\int_{\Omega}\mu(c^{n})\bigtriangledown
E^{I,n,\theta}_{\textbf{u}}:\bigtriangledown\textbf{v}_{h}+\int_{\Omega}\mu(c^{n})\bigtriangledown
E^{A,n,\theta}_{\textbf{u}}:\bigtriangledown\textbf{v}_{h}\\\
&\quad+\int_{\Omega}\sigma E^{A,n,\theta}_{u1}v_{1h}+\int_{\Omega}\sigma
E^{A,n,\theta}_{u2}v_{2h}-\int_{\Omega}(\bigtriangledown\cdot\textbf{v}_{h})(E^{I,n,\theta}_{p}+E^{A,n,\theta}_{p})\\\
&\quad+\int_{\Omega}(\bigtriangledown\cdot
E^{A,n,\theta}_{\textbf{u}})q_{h}+\int_{\Omega}\tilde{\bigtriangledown}E^{I,n,\theta}_{c}\cdot\bigtriangledown
d_{h}+\int_{\Omega}\tilde{\bigtriangledown}E^{A,n,\theta}_{c}\cdot\bigtriangledown
d_{h}\\\ &\quad+\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}\alpha E^{A,n,\theta}_{c}d_{h}\end{split}$
(26)

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

$\begin{split}I_{3}&=\sum_{k=1}^{n_{el}}(\tau_{k}^{\prime}(M\partial_{t}(\textbf{U}^{n}-\textbf{U}^{n}_{h})+\mathcal{L}(\textbf{U}^{n,\theta}-\textbf{U}_{h}^{n,\theta})),-\mathcal{L}^{*}\textbf{V}_{h})_{\Omega_{k}}\\\
&=\sum_{k=1}^{n_{el}}\\{(\tau_{1}^{\prime}(-\mu(c)\Delta(E^{I,n,\theta}_{u1}+E_{u1}^{A,n,\theta})+\sigma(E^{I,n,\theta}_{u1}+E_{u1}^{A,n,\theta})+\frac{\partial(E^{I,n,\theta}_{p})}{\partial
x}+\frac{\partial(E^{A,n,\theta}_{p})}{\partial x}),\\\ &\quad(\mu(c)\Delta
v_{1h}-\sigma v_{1h}+\frac{\partial q_{h}}{\partial
x}))_{\Omega_{k}}+(\tau_{1}^{\prime}(-\mu(c)\Delta(E^{I,n,\theta}_{u2}+E_{u2}^{A,n,\theta})+\sigma(E^{I,n,\theta}_{u2}+E_{u2}^{A,n,\theta})+\\\
&\quad\frac{\partial(E^{I,n,\theta}_{p}+E_{p}^{A,n,\theta})}{\partial
y}),(\mu(c)\Delta v_{2h}-\sigma v_{2h}+\frac{\partial q_{h}}{\partial
y}))_{\Omega_{k}}+(\tau_{2}^{\prime}\bigtriangledown\cdot(E^{I,n,\theta}_{\textbf{u}}+E_{\textbf{u}}^{A,n,\theta}),\bigtriangledown\cdot\textbf{v}_{h})_{\Omega_{k}}\\\
&\quad+(\tau_{3}^{\prime}(\partial_{t}(E^{I,n}_{c}+E_{c}^{A,n})-\bigtriangledown\cdot\tilde{\bigtriangledown}(E^{I,n,\theta}_{c}+E_{c}^{A,n,\theta})+\textbf{u}^{n}\cdot\bigtriangledown(E^{I,n,\theta}_{c}+E_{c}^{A,n,\theta})+\\\
&\quad\alpha(E^{I,n,\theta}_{c}+E_{c}^{A,n,\theta})-d_{4}),\bigtriangledown\cdot\tilde{\bigtriangledown}d_{h}+\textbf{u}^{n}\cdot\bigtriangledown
d_{h}-\alpha d_{h})_{\Omega_{k}}\\}\\\
&=I_{3}^{1}+I_{3}^{2}+I_{3}^{3}+I_{3}^{4}\end{split}$ (27)

where $I_{3}^{1},I_{3}^{2},I_{3}^{3},I_{3}^{4}$ are four terms of $I_{3}$
which we will discuss in the later part of the proof.

$\begin{split}I_{4}&=\sum_{k=1}^{n_{el}}\tau_{k}^{\prime}(-\textbf{d},-\mathcal{L}^{*}\textbf{V}_{h})_{\Omega_{k}}=\sum_{k=1}^{n_{el}}\\{\tau_{3}^{\prime}(d_{4},\bigtriangledown\cdot\tilde{\bigtriangledown}d_{h}+\textbf{u}^{n}\cdot\bigtriangledown
d_{h}-\alpha d_{h})_{\Omega_{k}}\\}\end{split}$ (28)

and since
$(1-\tau_{1}^{-1}\tau_{1}^{\prime})=0=(1-\tau_{2}^{-1}\tau_{2}^{\prime})$ the
next term will take the following form

$\begin{split}I_{5}&=\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,\theta}-\textbf{U}^{n,\theta}_{h})),-\textbf{V}_{h})_{\Omega_{k}}\\\
&=\sum_{k=1}^{n_{el}}((1-\tau_{3}^{-1}\tau_{3}^{\prime})(\partial_{t}E^{I,n}_{c}+\partial_{t}E^{A,n}_{c}-\bigtriangledown\cdot\tilde{\bigtriangledown}(E^{I,n,\theta}_{c}+E_{c}^{A,n,\theta})+\textbf{u}\cdot\bigtriangledown(E^{I,n,\theta}_{c}+E_{c}^{A,n,\theta})+\\\
&\quad\alpha(E^{I,n,\theta}_{c}+E_{c}^{A,n,\theta})),-d_{h})_{\Omega_{k}}\end{split}$
(29)

and the last term,

$\begin{split}I_{6}&=\sum_{k=1}^{n_{el}}(\tau_{k}^{-1}\tau_{k}^{\prime}\textbf{d},\textbf{V}_{h})_{\Omega_{k}}=\sum_{k=1}^{n_{el}}(\tau_{3}^{-1}\tau_{3}^{\prime}d_{4},d_{h})_{\Omega_{k}}\end{split}$
(30)

Now taking all these terms together, (23) becomes

$I_{1}+I_{2}+I_{3}+I_{4}+I_{5}+I_{6}=(TE^{n,\theta},d_{h}),\hskip
2.84526pt\forall\hskip 2.84526pt\textbf{V}_{h}\in V_{s}^{h}\times
V_{s}^{h}\times Q_{s}^{h}\times V_{s}^{h}\\\ $ (31)

This implies

$I_{1}+\int_{\Omega}\mu(c^{n})\bigtriangledown
E^{A,n,\theta}_{\textbf{u}}:\bigtriangledown\textbf{v}_{h}+\int_{\Omega}\tilde{\bigtriangledown}E^{A,n,\theta}_{c}\cdot\bigtriangledown
d_{h}+\int_{\Omega}\sigma E^{A,n,\theta}_{u1}v_{1h}+\int_{\Omega}\sigma
E^{A,n,\theta}_{u2}v_{2h}\\\ +\int_{\Omega}\alpha
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}\tilde{\bigtriangledown}E^{I,n,\theta}_{c}\cdot\bigtriangledown
d_{h}-\int_{\Omega}d_{h}\textbf{u}^{n}\cdot\bigtriangledown(E^{I,n,\theta}_{c}+E^{A,n,\theta}_{c})\\\
-\int_{\Omega}\mu(c^{n})\bigtriangledown
E^{I,n,\theta}_{\textbf{u}}:\bigtriangledown\textbf{v}_{h}-I_{3}-I_{4}-I_{5}-I_{6}+(TE^{n,\theta},d_{h})\\\
\hskip 2.84526pt\forall\hskip 2.84526pt\textbf{V}_{h}\in V_{s}^{h}\times
V_{s}^{h}\times Q_{s}^{h}\times V_{s}^{h}$ (32)

Now we will treat each term separately to find out the 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.
Let us start with $I_{1}$ as follows:

$\begin{split}(\frac{E^{A,n+1}_{c}-E^{A,n}_{c}}{dt},E^{A,n,\theta}_{c})&=(\frac{E^{A,n+1}_{c}-E^{A,n}_{c}}{dt},\frac{1+\theta}{2}E^{A,n+1}_{c}+\frac{1-\theta}{2}E^{A,n}_{c})\\\
&=\frac{1+\theta}{2dt}\|E^{A,n+1}_{c}\|^{2}-\frac{1-\theta}{2dt}\|E^{A,n}_{c}\|^{2}-\frac{\theta}{dt}(E^{A,n}_{c},E^{A,n+1}_{c})\\\
&=\frac{1}{2dt}(\|E^{A,n+1}_{c}\|^{2}-\|E^{A,n}_{c}\|^{2})+\frac{\theta}{2dt}(\|E^{A,n+1}_{c}\|^{2}-\|E^{A,n}_{c}\|^{2})^{2}\\\
&\geq\frac{1}{2dt}(\|E^{A,n+1}_{c}\|^{2}-\|E^{A,n}_{c}\|^{2})\end{split}$ (33)

Let us estimate the remaining terms of $LHS$ as follows:

$\begin{split}\int_{\Omega}\mu(c^{n})\bigtriangledown
E^{A,n,\theta}_{\textbf{u}}:\bigtriangledown
E^{A,n,\theta}_{\textbf{u}}&=\int_{\Omega}\mu(c^{n})\\{\sum_{i=1}^{2}(\frac{\partial
E^{A,n,\theta}_{ui}}{\partial x})^{2}+\sum_{i=1}^{2}(\frac{\partial
E^{A,n,\theta}_{ui}}{\partial y})^{2}\\}\\\ &\geq\mu_{l}\\{\|\frac{\partial
E^{A,n,\theta}_{u1}}{\partial x}\|^{2}+\|\frac{\partial
E^{A,n,\theta}_{u1}}{\partial y}\|^{2}+\|\frac{\partial
E^{A,n,\theta}_{u2}}{\partial x}\|^{2}+\|\frac{\partial
E^{A,n,\theta}_{u2}}{\partial y}\|^{2}\\}\end{split}$ (34)

and

$\begin{split}\int_{\Omega}\tilde{\bigtriangledown}E^{A,n,\theta}_{c}\cdot\bigtriangledown
E^{A,n,\theta}_{c}&=\int_{\Omega}D_{1}(\frac{\partial
E^{A,n,\theta}_{c}}{\partial x})^{2}+\int_{\Omega}D_{2}(\frac{\partial
E^{A,n,\theta}_{c}}{\partial y})^{2}\\\ &\geq D_{l}\\{\|\frac{\partial
E^{A,n,\theta}_{c}}{\partial x}\|^{2}+\|\frac{\partial
E^{A,n,\theta}_{c}}{\partial y}\|^{2}\\}\end{split}$ (35)

where $D_{l}$= min
$\\{\underset{\Omega}{inf}D_{1},\underset{\Omega}{inf}D_{2}\\}$.
Another few terms of $I_{2}$ can be easily simplified as,

$\begin{split}\int_{\Omega}\sigma
E^{A,n,\theta}_{ui}E^{A,n,\theta}_{ui}&=\sigma\|E^{A,n,\theta}_{ui}\|^{2}\hskip
5.69054ptfor\hskip 2.84526pti=1,2\\\ \int_{\Omega}\alpha
E^{A,n,\theta}_{c}E^{A,n,\theta}_{c}&=\alpha\|E^{A,n,\theta}_{c}\|^{2}\\\
\int_{\Omega}(\bigtriangledown\cdot
E^{A,n,\theta}_{\textbf{u}})(E^{I,n,\theta}_{p}+E^{A,n,\theta}_{p})-\int_{\Omega}(\bigtriangledown\cdot
E^{A,n,\theta}_{\textbf{u}})E^{A,n,\theta}_{p}&=\int_{\Omega}(\bigtriangledown\cdot
E^{A,n,\theta}_{\textbf{u}})E^{I,n,\theta}_{p}\end{split}$ (36)

Combining all these inequalities (32) becomes,

$\frac{1}{2dt}(\|E^{A,n+1}_{c}\|^{2}-\|E^{A,n}_{c}\|^{2})+\mu_{l}\\{\|\frac{\partial
E^{A,n,\theta}_{u1}}{\partial x}\|^{2}+\|\frac{\partial
E^{A,n,\theta}_{u1}}{\partial y}\|^{2}+\|\frac{\partial
E^{A,n,\theta}_{u2}}{\partial x}\|^{2}\\\ +\|\frac{\partial
E^{A,n,\theta}_{u2}}{\partial y}\|^{2}\\}+D_{l}\\{\|\frac{\partial
E^{A,n,\theta}_{c}}{\partial x}\|^{2}+\|\frac{\partial
E^{A,n,\theta}_{c}}{\partial y}\|^{2}\\}+\sigma\|E^{A,n,\theta}_{u1}\|^{2}+\\\
\sigma\|E^{A,n,\theta}_{u2}\|^{2}+\alpha\|E^{A,n,\theta}_{c}\|^{2}\\\
\leq\int_{\Omega}(\bigtriangledown\cdot
E^{A,n,\theta}_{\textbf{u}})E^{I,n,\theta}_{p}-\int_{\Omega}\tilde{\bigtriangledown}E^{I,n,\theta}_{c}\cdot\bigtriangledown
E^{A,n,\theta}_{c}-\int_{\Omega}E^{A,n,\theta}_{c}\textbf{u}\cdot\bigtriangledown
E^{A,n,\theta}_{c}\\\
-\int_{\Omega}E^{A,n,\theta}_{c}\textbf{u}\cdot\bigtriangledown
E^{I,n,\theta}_{c}-\int_{\Omega}\mu(c^{n})\bigtriangledown
E^{I,n,\theta}_{\textbf{u}}:\bigtriangledown
E^{A,n,\theta}_{\textbf{u}}-I_{3}-I_{4}-I_{5}-I_{6}\\\
+(TE^{n,\theta},E^{A,n,\theta}_{c})\hskip 227.62204pt$ (37)

Now we will find upper bounds of the terms in the $RHS$ of the above
equation.We will use $Cauchy-Schwarz$ and $Young^{\prime}s$ inequality to
reach at the desired bounds.Let us start with the first term as follows:

$\begin{split}\int_{\Omega}(\bigtriangledown\cdot
E^{A,n,\theta}_{\textbf{u}})E^{I,n,\theta}_{p}&=\int_{\Omega}(\frac{\partial
E^{A,n,\theta}_{u1}}{\partial x}+\frac{\partial E^{A,n,\theta}_{u2}}{\partial
x})E^{I,n,\theta}_{p}\\\ &\quad(applying\hskip 2.84526ptCauchy-Schwarz\hskip
2.84526ptinequality)\\\ &\leq(\|\frac{\partial E^{A,n,\theta}_{u1}}{\partial
x}\|+\|\frac{\partial E^{A,n,\theta}_{u2}}{\partial
y}\|)\|E^{I,n,\theta}_{p}\|\\\ &\quad(applying\hskip
2.84526ptYoung^{\prime}s\hskip 2.84526ptinequality\hskip 2.84526ptfor\hskip
2.84526pteach\hskip 2.84526ptof\hskip 2.84526ptthe\hskip 2.84526pttwo\hskip
2.84526ptterms)\\\ &\leq\frac{1}{2\epsilon_{1}}(\|\frac{\partial
E^{A,n,\theta}_{u1}}{\partial x}\|^{2}+\|\frac{\partial
E^{A,n,\theta}_{u2}}{\partial
y}\|^{2})+\epsilon_{1}\|E^{I,n,\theta}_{p}\|^{2}\\\
&\leq\frac{1}{2\epsilon_{1}}(\|\frac{\partial E^{A,n,\theta}_{u1}}{\partial
x}\|^{2}+\|\frac{\partial E^{A,n,\theta}_{u2}}{\partial
y}\|^{2})+\epsilon_{1}(\frac{1+\theta}{2}\|E^{I,n+1}_{p}\|+\frac{1-\theta}{2}\|E^{I,n}_{p}\|)^{2}\\\
&\leq\frac{1}{2\epsilon_{1}}(\|\frac{\partial E^{A,n,\theta}_{u1}}{\partial
x}\|^{2}+\|\frac{\partial E^{A,n,\theta}_{u2}}{\partial
y}\|^{2})+\epsilon_{1}C^{2}h^{2}(\frac{1+\theta}{2}\|p^{n+1}\|_{1}+\frac{1-\theta}{2}\|p^{n}\|_{1})^{2}\\\
\end{split}$ (38)

Similarly for each term we will use $Cauchy-Schwarz$ inequality and
$Young^{\prime}s$ inequality wherever it will be needed, but without
mentioning about them now onwards.
Proceeding in the same way the second term becomes

$\begin{split}-\int_{\Omega}\tilde{\bigtriangledown}E^{I,n,\theta}_{c}\cdot\bigtriangledown
E^{A,n,\theta}_{c}&=-\int_{\Omega}(D_{1}\frac{\partial
E^{I,n,\theta}_{c}}{\partial x}\frac{\partial E^{A,n,\theta}_{c}}{\partial
x}+D_{2}\frac{\partial E^{I,n,\theta}_{c}}{\partial y}\frac{\partial
E^{A,n,\theta}_{c}}{\partial y})\\\
&\leq\frac{D_{m}}{2\epsilon_{2}}(\|\frac{\partial E^{A,n,\theta}_{c}}{\partial
x}\|^{2}+\|\frac{\partial E^{A,n,\theta}_{c}}{\partial
y}\|^{2})+\frac{D_{m}\epsilon_{2}}{2}\mid E^{I,n,\theta}_{c}\mid_{1}^{2}\\\
&\leq\frac{D_{m}}{2\epsilon_{2}}(\|\frac{\partial E^{A,n,\theta}_{c}}{\partial
x}\|^{2}+\|\frac{\partial E^{A,n,\theta}_{c}}{\partial y}\|^{2})+\\\
&\quad\frac{D_{m}\epsilon_{2}}{2}C^{2}h^{2}(\frac{1+\theta}{2}\|c^{n+1}\|_{2}+\frac{1-\theta}{2}\|c^{n}\|_{2})^{2}\end{split}$
(39)

where $D_{m}$= max
$\\{\underset{\Omega}{sup}D_{1},\underset{\Omega}{sup}D_{2}\\}$
Next term,

$\begin{split}-\int_{\Omega}E^{A,n,\theta}_{c}\textbf{u}^{n}\cdot\bigtriangledown
E^{A,n,\theta}_{c}&=-\int_{\Omega}(u_{1}^{n}E^{A,n,\theta}_{c}\frac{\partial
E^{A,n,\theta}_{c}}{\partial x}+u_{2}^{n}E^{A,n,\theta}_{c}\frac{\partial
E^{A,n,\theta}_{c}}{\partial y})\\\
&\leq\frac{1}{2\epsilon_{3}}(C_{1}^{n}\|\frac{\partial
E^{A,n,\theta}_{c}}{\partial x}\|^{2}+C_{2}^{n}\|\frac{\partial
E^{A,n,\theta}_{c}}{\partial
y}\|^{2})+\frac{\epsilon_{3}}{2}(C_{1}^{n}+C_{2}^{n})\\\
&\quad\|E^{A,n,\theta}_{c}\|^{2}\end{split}$ (40)

where $C_{1}^{n}$= $\underset{\Omega}{sup}$ $\mid u_{1}^{n}\mid$ and
$C_{2}^{n}$= $\underset{\Omega}{sup}$ $\mid u_{2}^{n}\mid$ (applying
assumption (v))
Similarly the next term

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

The next term,

$\begin{split}-\int_{\Omega}\mu(c^{n})\bigtriangledown
E^{I,n,\theta}_{\textbf{u}}:\bigtriangledown
E^{A,n,\theta}_{\textbf{u}}&\leq\mu_{u}(\|\frac{\partial
E^{I,n,\theta}_{u1}}{\partial x}\|\|\frac{\partial
E^{A,n,\theta}_{u1}}{\partial x}\|+\|\frac{\partial
E^{I,n,\theta}_{u1}}{\partial y}\|\|\frac{\partial
E^{A,n,\theta}_{u1}}{\partial y}\|\\\ &\quad+\|\frac{\partial
E^{I,n,\theta}_{u2}}{\partial x}\|\|\frac{\partial
E^{A,n,\theta}_{u2}}{\partial x}\|+\|\frac{\partial
E^{I,n,\theta}_{u2}}{\partial y}\|\|\frac{\partial
E^{A,n,\theta}_{u2}}{\partial y}\|)\\\ \end{split}$ (42)
$\begin{split}&\leq\frac{\epsilon_{4}\mu_{u}}{2}\sum_{i=1}^{2}(\|\frac{\partial
E^{I,n,\theta}_{ui}}{\partial x}\|^{2}+\|\frac{\partial
E^{I,n,\theta}_{ui}}{\partial y}\|^{2})+\\\
&\quad\frac{\mu_{u}}{2\epsilon_{4}}\sum_{i=1}^{2}(\|\frac{\partial
E^{A,n,\theta}_{ui}}{\partial x}\|^{2}+\|\frac{\partial
E^{A,n,\theta}_{ui}}{\partial y}\|^{2})\\\
&\leq\epsilon_{4}\mu_{u}\sum_{i=1}^{2}C^{2}h^{2}(\frac{1+\theta}{2}\|u_{i}^{n+1}\|_{2}+\frac{1-\theta}{2}\|u_{i}^{n}\|_{2})^{2}+\\\
&\quad\frac{\mu_{u}}{2\epsilon_{4}}\sum_{i=1}^{2}(\|\frac{\partial
E^{A,n,\theta}_{ui}}{\partial x}\|^{2}+\|\frac{\partial
E^{A,n,\theta}_{ui}}{\partial y}\|^{2})\end{split}$ (43)

Now we will find bounds for each remaining term of $I_{3}$. Before going to
further calculations let us mention an important observation:

######  Observation 1.

According to the choice of the finite element spaces $V_{s}^{h}$ and
$Q_{s}^{h}$, we can clearly say that over each element sub-domain 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 take the first term of $(-I_{3})$ along with earlier mentioned
replacements. $I_{3}$ has four terms and we will find bounds for each of them
separately. We have already denoted them by the notations
$I_{3}^{1},I_{3}^{2},I_{3}^{3},I_{3}^{4}$. Here we start with $I_{3}^{1}$,

$-I_{3}^{1}=-\sum_{k=1}^{n_{el}}(\tau_{1}^{\prime}(-\mu(c^{n})\Delta(E^{I,n,\theta}_{u1}+E_{u1}^{A,n,\theta})+\sigma(E^{I,n,\theta}_{u1}+E_{u1}^{A,n,\theta})+\frac{\partial(E^{I,n,\theta}_{p})}{\partial
x}+\\\ \quad\frac{\partial(E^{A,n,\theta}_{p})}{\partial x}),\mu(c^{n})\Delta
E_{u1}^{A,n,\theta}-\sigma E_{u1}^{A,n,\theta}+\frac{\partial
E_{p}^{A,n,\theta}}{\partial x})_{\Omega_{k}}\\\
=\tau_{1}^{\prime}(\mu(c^{n})\Delta E^{I,n,\theta}_{u1}-\sigma
E^{I,n,\theta}_{u1}-\frac{\partial(E^{I,n,\theta}_{p})}{\partial
x},\mu(c^{n})\Delta E_{u1}^{A,n,\theta}-\sigma
E_{u1}^{A,n,\theta}+\frac{\partial E_{p}^{A,n,\theta}}{\partial
x})_{\tilde{\Omega}}\\\ \quad+\tau_{1}^{\prime}(\mu(c^{n})\Delta
E^{A,n,\theta}_{u1}-\sigma
E^{A,n,\theta}_{u1}-\frac{\partial(E^{A,n,\theta}_{p})}{\partial
x},\mu(c^{n})\Delta E_{u1}^{A,n,\theta}-\sigma
E_{u1}^{A,n,\theta}+\frac{\partial E_{p}^{A,n,\theta}}{\partial
x})_{\tilde{\Omega}}\\\ $ (44)

We calculate the bounds for the above two terms separately. Applying $Cauchy-
Schwarz$ inequality on each term the first part is as follows:

$\tau_{1}^{\prime}(\mu(c^{n})\Delta E^{I,n,\theta}_{u1}-\sigma
E^{I,n,\theta}_{u1}-\frac{\partial(E^{I,n,\theta}_{p})}{\partial
x},\mu(c^{n})\Delta E_{u1}^{A,n,\theta}-\sigma
E_{u1}^{A,n,\theta}+\frac{\partial E_{p}^{A,n,\theta}}{\partial
x})_{\tilde{\Omega}}\\\ $ (45)
$\leq\sum_{k=1}^{n_{el}}\mid\tau_{1}\mid(\mu_{u}^{2}\|\Delta
E^{I,n,\theta}_{u1}\|_{k}\|\Delta
E^{A,n,\theta}_{u1}\|_{k}+\sigma\mu_{u}\|E^{I,n,\theta}_{u1}\|_{k}\|\Delta
E^{A,n,\theta}_{u1}\|_{k}+\mu_{u}\\\ \|\frac{\partial
E_{p}^{I,n,\theta}}{\partial x}\|_{k}\|\Delta
E_{u1}^{A,n,\theta}\|_{k}+\sigma\mu_{u}\|E^{A,n,\theta}_{u1}\|_{k}\|\Delta
E_{u1}^{I,n,\theta}\|_{k}+\sigma^{2}\|E^{A,n,\theta}_{u1}\|_{k}\\\
\|E_{u1}^{I,n,\theta}\|_{k}+\sigma\|E^{A,n,\theta}_{u1}\|_{k}\|\frac{\partial
E_{p}^{I,n,\theta}}{\partial x}\|_{k}+\mu_{u}\|\Delta
E^{I,n,\theta}_{u1}\|_{k}\|\frac{\partial E_{p}^{A,n,\theta}}{\partial
x}\|_{k}+\\\ \sigma\|E^{I,n,\theta}_{u1}\|_{k}\|\frac{\partial
E_{p}^{A,n,\theta}}{\partial x}\|_{k}+\|\frac{\partial
E_{p}^{A,n,\theta}}{\partial x}\|_{k}\|\frac{\partial
E_{p}^{I,n,\theta}}{\partial x}\|_{k})\\\ $ (46)

Let $B_{1k},B_{2k},B_{3k}$ be the bounds on $E^{A,n,\theta}_{u1},\Delta
E^{A,n,\theta}_{u1},\frac{\partial E^{A,n,\theta}_{p}}{\partial x}$
respectively on each element sub domain under the above observation 1.

$\leq\mid\tau_{1}\mid\sum_{k=1}^{n_{el}}(\mu_{u}^{2}B_{2k}\|\Delta
E^{I,n,\theta}_{u1}\|_{k}+\sigma\mu_{u}B_{2k}\|E^{I,n,\theta}_{u1}\|_{k}+\mu_{u}B_{2k}\mid
E_{p}^{I,n,\theta}\mid_{1,k}+\\\ \sigma\mu_{u}B_{1k}\|\Delta
E_{u1}^{I,n,\theta}\|_{k}+\sigma^{2}B_{1k}\|E_{u1}^{I,n,\theta}\|_{k}+\sigma
B_{1k}\mid E_{p}^{I,n,\theta}\mid_{1,k}+\\\ \mu_{u}B_{3k}\|\Delta
E^{I,n,\theta}_{u1}\|_{k}+\sigma B_{3k}\|E^{I,n,\theta}_{u1}\|_{k}+B_{3k}\mid
E_{p}^{I,n,\theta}\mid_{1,k})\\\
\leq\mid\tau_{1}\mid\\{(\sum_{k=1}^{n_{el}}(\mu_{u}^{2}B_{2k}+\sigma\mu_{u}B_{1k}+\mu_{u}B_{3k}))C(\frac{1+\theta}{2}\|u_{1}^{n+1}\|_{2}+\frac{1-\theta}{2}\|u_{1}^{n}\|_{2})+\\\
(\sum_{k=1}^{n_{el}}\sigma(\mu_{u}B_{2k}+\sigma
B_{1k}+B_{3k}))Ch^{2}(\frac{1+\theta}{2}\|u_{1}^{n+1}\|_{2}+\frac{1-\theta}{2}\|u_{1}^{n}\|_{2})+\\\
(\sum_{k=1}^{n_{el}}(\mu_{u}B_{2k}+\sigma
B_{1k}+B_{3k}))Ch(\frac{1+\theta}{2}\|p^{n+1}\|_{1}+\frac{1-\theta}{2}\|p^{n}\|_{1})\\}$
(47)

This completes the first part. Now we see that the second part has alike
expression with auxiliary error terms in the place of interpolation error
terms. Hence proceeding in the same way as above and applying bounds for
elements belonging to $V_{s}^{h}$ and $Q_{s}^{h}$ spaces we will bound the
second part as follows:

$\tau_{1}^{\prime}(\mu(c)\Delta E^{A,n,\theta}_{u1}-\sigma
E^{A,n,\theta}_{u1}-\frac{\partial(E^{A,n,\theta}_{p})}{\partial
x},\mu(c)\Delta E_{u1}^{A,n,\theta}-\sigma E_{u1}^{A,n,\theta}+\frac{\partial
E_{p}^{A,n,\theta}}{\partial x})_{\tilde{\Omega}}\\\
\leq\mid\tau_{1}\mid\sum_{k=1}^{n_{el}}(\mu_{u}^{2}B_{2k}^{2}+2\sigma\mu_{u}B_{1k}B_{2k}+\mu_{u}B_{2k}B_{3k}+\sigma^{2}B_{1k}^{2}+\sigma
B_{1k}B_{3k}+\\\ \mu_{u}B_{3k}B_{2k}+\sigma B_{1k}B_{3k}+B_{3k}^{2})\\\
\leq\mid\tau_{1}\mid\sum_{k=1}^{n_{el}}M_{1k}\hskip 227.62204pt$ (48)

where
$M_{1k}=(\mu_{u}^{2}B_{2k}^{2}+2\sigma\mu_{u}B_{1k}B_{2k}+\mu_{u}B_{2k}B_{3k}+\sigma^{2}B_{1k}^{2}+\sigma
B_{1k}B_{3k}+\mu_{u}B_{3k}B_{2k}+\sigma B_{1k}B_{3k}+B_{3k}^{2}$)
Combining all these results and putting into (42) we will have

$\begin{split}-I_{3}^{1}&\leq\mid\tau_{1}\mid\\{(\sum_{k=1}^{n_{el}}(\mu_{u}^{2}B_{2k}+\sigma\mu_{u}B_{1k}+\mu_{u}B_{3k}))C(\frac{1+\theta}{2}\|u_{1}^{n+1}\|_{2}+\\\
&\quad\frac{1-\theta}{2}\|u_{1}^{n}\|_{2})+(\sum_{k=1}^{n_{el}}\sigma(\mu_{u}B_{2k}+\sigma
B_{1k}+B_{3k}))Ch^{2}(\frac{1+\theta}{2}\|u_{1}^{n+1}\|_{2}+\\\
&\quad\frac{1-\theta}{2}\|u_{1}^{n}\|_{2})+(\sum_{k=1}^{n_{el}}(\mu_{u}B_{2k}+\sigma
B_{1k}+B_{3k}))Ch(\frac{1+\theta}{2}\|p^{n+1}\|_{1}+\\\
&\quad\frac{1-\theta}{2}\|p^{n}\|_{1})\\}+\mid\tau_{1}\mid\sum_{k=1}^{n_{el}}M_{1k}\end{split}$
(49)

This completes the derivation of bound on the first term of $(-I_{3})$. Now we
see that the second term of $I_{3}$ in (27) is exactly similar to its first
term, only the subscripts are different that is $u_{2}$ replaces $u_{1}$ in
subscript. Therefore considering the constants
$B_{1k}^{\prime},B_{2k}^{\prime},B_{3k}^{\prime}$ as the bounds for
$E^{A,n,\theta}_{u2},\Delta E^{A,n,\theta}_{u2},\frac{\partial
E_{p}^{A,n,\theta}}{\partial y}$ respectively on each element sub domain, we
can bound the term as follows:

$\begin{split}-I_{3}^{2}&\leq\mid\tau_{1}\mid\mid\\{(\sum_{k=1}^{n_{el}}(\mu_{u}^{2}B_{2k}^{\prime}+\sigma\mu_{u}B_{1k}^{\prime}+\mu_{u}B_{3k}^{\prime}))C(\frac{1+\theta}{2}\|u_{2}^{n+1}\|_{2}+\\\
&\quad\frac{1-\theta}{2}\|u_{2}^{n}\|_{2})+(\sum_{k=1}^{n_{el}}\sigma(\mu_{u}B_{2k}^{\prime}+\sigma
B_{1k}^{\prime}+B_{3k}^{\prime}))Ch^{2}(\frac{1+\theta}{2}\|u_{2}^{n+1}\|_{2}+\\\
&\quad\frac{1-\theta}{2}\|u_{2}^{n}\|_{2})+(\sum_{k=1}^{n_{el}}(\mu_{u}B_{2k}^{\prime}+\sigma^{2}B_{1k}^{\prime}+B_{3k}^{\prime}))Ch(\frac{1+\theta}{2}\|p^{n+1}\|_{1}+\\\
&\quad\frac{1-\theta}{2}\|p^{n}\|_{1})\\}+\mid\tau_{1}\mid\sum_{k=1}^{n_{el}}M_{2k}\end{split}$
(50)

where
$M_{2k}=(\mu_{u}^{2}B_{2k}^{{}^{\prime}2}+2\sigma\mu_{u}B^{\prime}_{1k}B^{\prime}_{2k}+\mu_{u}B^{\prime}_{2k}B^{\prime}_{3k}+\sigma^{2}B_{1k}^{{}^{\prime}2}+\sigma
B^{\prime}_{1k}B^{\prime}_{3k}+\mu_{u}B^{\prime}_{3k}B^{\prime}_{2k}+\sigma
B^{\prime}_{1k}B^{\prime}_{3k}+B_{3k}^{{}^{\prime}2})$
Now we are going to derive bounds for the third term of $I_{3}$ as follows:

$\begin{split}-I_{3}^{3}&=-\tau_{2}^{\prime}\sum_{k=1}^{n_{el}}(\bigtriangledown\cdot
E^{I,n,\theta}_{\textbf{u}},\bigtriangledown\cdot
E_{\textbf{u}}^{A,n,\theta})_{\Omega_{k}}-\tau_{2}^{\prime}\sum_{k=1}^{n_{el}}(\bigtriangledown\cdot
E^{A,n,\theta}_{\textbf{u}},\bigtriangledown\cdot
E_{\textbf{u}}^{A,n,\theta})_{\Omega_{k}}\\\
&\leq\mid\tau_{2}\mid\sum_{k=1}^{n_{el}}(\|\frac{\partial
E^{I,n,\theta}_{u1}}{\partial x}\|_{k}\|\frac{\partial
E^{A,n,\theta}_{u1}}{\partial x}\|_{k}+\|\frac{\partial
E^{I,n,\theta}_{u2}}{\partial y}\|_{k}\|\frac{\partial
E^{A,n,\theta}_{u2}}{\partial y}\|_{k}+\|\frac{\partial
E^{A,n,\theta}_{u1}}{\partial x}\|_{k}^{2}+\\\ &\quad\|\frac{\partial
E^{I,n,\theta}_{u1}}{\partial x}\|_{k}\|\frac{\partial
E^{A,n,\theta}_{u2}}{\partial y}\|_{k}+\|\frac{\partial
E^{A,n,\theta}_{u1}}{\partial x}\|_{k}\|\frac{\partial
E^{I,n,\theta}_{u2}}{\partial y}\|_{k}+2\|\frac{\partial
E^{A,n,\theta}_{u1}}{\partial x}\|_{k}\|\frac{\partial
E^{A,n,\theta}_{u2}}{\partial y}\|_{k}+\\\ \end{split}$ (51)
$\begin{split}&\quad\|\frac{\partial E^{A,n,\theta}_{u2}}{\partial
y}\|_{k}^{2})\\\
&\leq\mid\tau_{2}\mid\sum_{k=1}^{n_{el}}((B_{4k}+B_{5k}^{\prime})(\|\frac{\partial
E^{I,n,\theta}_{u1}}{\partial x}\|_{k}+\|\frac{\partial
E^{I,n,\theta}_{u2}}{\partial y}\|_{k})+\mid\tau_{2}\mid
C_{1}\\{(1+\epsilon_{5})\|\frac{\partial E^{A,n,\theta}_{u1}}{\partial
x}\|^{2}\\\ &\quad+(1+\epsilon_{5})\|\frac{\partial
E^{A,n,\theta}_{u2}}{\partial y}\|^{2}\\}\\\
&\leq\mid\tau_{2}\mid(\sum_{k=1}^{n_{el}}(B_{4k}+B_{5k}^{\prime}))(\|E^{I,n,\theta}_{u1}\|_{1}+\|E^{I,n,\theta}_{u2}\|_{1})+C_{\tau_{2}}C_{1}\\{(1+\epsilon_{5})\|\frac{\partial
E^{A,n,\theta}_{u1}}{\partial x}\|^{2}+\\\
&\quad(1+\epsilon_{5})\|\frac{\partial E^{A,n,\theta}_{u2}}{\partial
y}\|^{2}\\}\\\
&\leq\mid\tau_{2}\mid(\sum_{k=1}^{n_{el}}(B_{4k}+B_{5k}^{\prime}))Ch\\{(\frac{1+\theta}{2}\|u_{1}^{n+1}\|_{2}+\frac{1-\theta}{2}\|u_{1}^{n}\|_{2})+(\frac{1+\theta}{2}\|u_{2}^{n+1}\|_{2}+\\\
&\quad\frac{1-\theta}{2}\|u_{2}^{n}\|_{2})\\}+C_{\tau_{2}}C_{1}\\{(1+\epsilon_{5})\|\frac{\partial
E^{A,n,\theta}_{u1}}{\partial x}\|^{2}+(1+\epsilon_{5})\|\frac{\partial
E^{A,n,\theta}_{u2}}{\partial y}\|^{2}\\}\end{split}$ (52)

where the constants $B_{4k},B_{5k},B_{4k}^{\prime}$ and $B_{5k}^{\prime}$ are
bounds on $\frac{\partial E^{A,n,\theta}_{u1}}{\partial x},\frac{\partial
E^{A,n,\theta}_{u1}}{\partial y},\frac{\partial E^{A,n,\theta}_{u2}}{\partial
x}$ and $\frac{\partial E^{A,n,\theta}_{u2}}{\partial y}$ respectively on each
element sub domain and $C_{\tau_{2}}$ is the maximum numerical value for
$\tau_{2}$ over $\Omega$. Now we will focus on the fourth term of $I_{3}$. We
will divide $I_{3}^{4}$ into three parts $P_{1},P_{2}$ and $P_{3}$ and then
calculate bounds for each of them separately.

$\begin{split}-I_{3}^{4}&=\sum_{k=1}^{n_{el}}\tau_{3}^{\prime}(\partial_{t}(E^{I,n}_{c}+E_{c}^{A,n}),-\bigtriangledown\cdot\tilde{\bigtriangledown}E_{c}^{A,n,\theta}-\textbf{u}\cdot\bigtriangledown
E_{c}^{A,n,\theta}+\alpha
E_{c}^{A,n,\theta})_{\Omega_{k}}+\sum_{k=1}^{n_{el}}\tau_{3}^{\prime}\\\
&\quad(\bigtriangledown\cdot\tilde{\bigtriangledown}E^{I,n,\theta}_{c}-\textbf{u}\cdot\bigtriangledown
E^{I,n,\theta}_{c}-\alpha
E^{I,n,\theta}_{c},\bigtriangledown\cdot\tilde{\bigtriangledown}E_{c}^{A,n,\theta}+\textbf{u}\cdot\bigtriangledown
E_{c}^{A,n,\theta}-\alpha E_{c}^{A,n,\theta})_{\Omega_{k}}\\\
&\quad+\sum_{k=1}^{n_{el}}\tau_{3}^{\prime}(\bigtriangledown\cdot\tilde{\bigtriangledown}E^{A,n,\theta}_{c}-\textbf{u}\cdot\bigtriangledown
E^{A,n,\theta}_{c}-\alpha
E^{A,n,\theta}_{c},\bigtriangledown\cdot\tilde{\bigtriangledown}E_{c}^{A,n,\theta}+\textbf{u}\cdot\bigtriangledown
E_{c}^{A,n,\theta}-\\\ &\quad\alpha E_{c}^{A,n,\theta})_{\Omega_{k}}\\\
&=P_{1}+P_{2}+P_{3}\end{split}$ (53)

Let us start with $P_{1}$

$\begin{split}P_{1}&=\sum_{k=1}^{n_{el}}\tau_{3}^{\prime}\alpha(\partial_{t}E^{I,n}_{c}+\partial_{t}E_{c}^{A,n},E_{c}^{A,n,\theta})_{\Omega_{k}}-\sum_{k=1}^{n_{el}}\tau_{3}^{\prime}(\partial_{t}E^{I,n}_{c}+\partial_{t}E_{c}^{A,n},\\\
&\quad\bigtriangledown\cdot\tilde{\bigtriangledown}E_{c}^{A,n,\theta}+\textbf{u}\cdot\bigtriangledown
E_{c}^{A,n,\theta})_{\Omega_{k}}\\\
&=\alpha\tau_{3}^{\prime}\sum_{k=1}^{n_{el}}\int_{\Omega_{k}}\frac{E^{A,n+1}_{c}-E^{A,n}_{c}}{dt}E_{c}^{A,n,\theta}-\tau_{3}^{\prime}\sum_{k=1}^{n_{el}}\int_{\Omega_{k}}\frac{E^{I,n+1}_{c}-E^{I,n}_{c}}{dt}(D_{1}\frac{\partial^{2}E_{c}^{A,n,\theta}}{\partial
x^{2}}+\\\ \end{split}$ (54) $\begin{split}&\quad
D_{2}\frac{\partial^{2}E_{c}^{A,n,\theta}}{\partial
y^{2}}+(u_{1}+\frac{\partial D_{1}}{\partial x})\frac{\partial
E_{c}^{A,n,\theta}}{\partial x}+(u_{2}+\frac{\partial D_{2}}{\partial
y})\frac{\partial E_{c}^{A,n,\theta}}{\partial y})-\\\
&\quad\tau_{3}^{\prime}\sum_{k=1}^{n_{el}}\int_{\Omega_{k}}\frac{E^{A,n+1}_{u1}-E^{A,n}_{u1}}{dt}(D_{1}\frac{\partial^{2}E_{c}^{A,n,\theta}}{\partial
x^{2}}+D_{2}\frac{\partial^{2}E_{c}^{A,n,\theta}}{\partial
y^{2}}+(u_{1}+\frac{\partial D_{1}}{\partial x})\frac{\partial
E_{c}^{A,n,\theta}}{\partial x}+\\\ &\quad(u_{2}+\frac{\partial
D_{2}}{\partial y})\frac{\partial E_{c}^{A,n,\theta}}{\partial y})\\\
&\leq\frac{\alpha\mid\tau_{3}^{\prime}\mid}{dt}(\sum_{k=1}^{n_{el}}B_{6k})(\|E^{A,n+1}_{c}\|^{2}-\|E^{A,n}_{c}\|^{2})+\frac{\mid\tau_{3}^{\prime}\mid}{dt}\\{\sum_{k=1}^{n_{el}}(D_{1m}B_{7k}+D_{2m}B_{7k}^{\prime}+D_{u1}\\\
&\quad
B_{8k}+D_{u2}B_{8k}^{\prime})\\}(\|E^{I,n+1}_{c}\|+\|E^{I,n}_{c}\|)+\frac{\mid\tau_{3}^{\prime}\mid}{dt}\\{\sum_{k=1}^{n_{el}}(D_{1m}B_{7k}+D_{2m}B_{7k}^{\prime}+D_{u1}B_{8k}+\\\
&\quad D_{u2}B_{8k}^{\prime})\\}(\|E^{A,n+1}_{c}\|^{2}-\|E^{A,n}_{c}\|^{2})\\\
&\leq\frac{C_{\tau_{3}}T}{T_{0}(T_{0}-C_{\tau_{3}})}\\{\sum_{k=1}^{n_{el}}(\alpha
B_{6k}+D_{B_{1k}})\\}(\|E^{A,n+1}_{c}\|^{2}-\|E^{A,n}_{c}\|^{2})+\frac{C_{\tau_{3}}Ch^{2}}{(T_{0}-C_{\tau_{3}})}\\{\sum_{k=1}^{n_{el}}D_{B_{1k}}\\}\\\
&\quad(\|c^{n+1}\|_{2}+\|c^{n}\|_{2})\end{split}$ (55)

where the constants $B_{6k},B_{7k},B_{7k}^{\prime},B_{8k}$ and
$B_{8k}^{\prime}$ are upper bounds on
$E^{A,n,\theta}_{c},\frac{\partial^{2}E^{A,n,\theta}_{c}}{\partial
x^{2}},\frac{\partial^{2}E^{A,n,\theta}_{c}}{\partial y^{2}},\\\
\frac{\partial E_{c}^{A,n,\theta}}{\partial x}$ and $\frac{\partial
E_{c}^{A,n,\theta}}{\partial y}$ respectively on each element sub domain and
$D_{1m},D_{2m},{D}_{u1},{D}_{u2}$ are maximum of the functions
$D_{1},D_{2},(\frac{\partial D_{1}}{\partial x}+u_{1}),(\frac{\partial
D_{2}}{\partial y}+u_{2})$ respectively over $\Omega$. $C_{\tau_{3}}$ and
$T_{0}$ are maximum bound of $\tau_{3}$ and minimum bound of time step $dt$
respectively. At the last line new notation $D_{B_{1k}}$ represents the big
sum.

$\begin{split}P_{2}&=\sum_{k=1}^{n_{el}}\tau_{3}^{\prime}(D_{1}\frac{\partial^{2}E^{I,n,\theta}_{c}}{\partial
x^{2}}+D_{2}\frac{\partial^{2}E^{I,n,\theta}_{c}}{\partial
y^{2}}+(\frac{\partial D_{1}}{\partial x}-u_{1})\frac{\partial
E^{I,n,\theta}_{c}}{\partial x}+(\frac{\partial D_{2}}{\partial y}-u_{2})\\\
&\quad\frac{\partial E^{I,n,\theta}_{c}}{\partial y}-\alpha
E^{I,n,\theta}_{c},D_{1}\frac{\partial^{2}E^{A,n,\theta}_{c}}{\partial
x^{2}}+D_{2}\frac{\partial^{2}E^{A,n,\theta}_{c}}{\partial
y^{2}}+(\frac{\partial D_{1}}{\partial x}+u_{1})\frac{\partial
E^{A,n,\theta}_{c}}{\partial x}+\\\ &\quad(\frac{\partial D_{2}}{\partial
y}+u_{2})\frac{\partial E^{A,n,\theta}_{c}}{\partial y}-\alpha
E^{A,n,\theta}_{c})_{\Omega_{k}}\\\ \end{split}$ (56)

Further simplifying and applying the bounds on auxiliary error terms over each
sub-domain, $P_{2}$ becomes

$\begin{split}&\leq\sum_{k=1}^{n_{el}}\mid\tau_{3}^{\prime}\mid(D_{1m}^{2}B_{7k}\|\frac{\partial^{2}E^{I,n,\theta}_{c}}{\partial
x^{2}}\|_{k}+D_{1m}D_{2m}B_{7k}^{\prime}\|\frac{\partial^{2}E^{I,n,\theta}_{c}}{\partial
y^{2}}\|_{k}+D_{1m}\bar{D}_{u1}B_{7k}\|\frac{\partial
E^{I,n,\theta}_{c}}{\partial x}\|_{k}\\\ &\quad
D_{1m}\bar{D}_{u2}B_{7k}\|\frac{\partial E^{I,n,\theta}_{c}}{\partial
x}\|_{k}+\alpha
D_{1m}B_{7k}\|E^{I,n,\theta}_{c}\|_{k}+D_{1m}D_{2m}B_{7k}^{\prime}\|\frac{\partial^{2}E^{I,n,\theta}_{c}}{\partial
x^{2}}\|_{k}+D_{2m}^{2}B_{7k}^{\prime}\\\ \end{split}$ (57)
$\begin{split}&\quad\|\frac{\partial^{2}E^{I,n,\theta}_{c}}{\partial
y^{2}}\|_{k}+D_{2m}\bar{D}_{u1}B_{7k}^{\prime}\|\frac{\partial^{2}E^{I,n,\theta}_{c}}{\partial
y^{2}}\|_{k}+D_{2m}\bar{D}_{u2}B_{7k}^{\prime}\|\frac{\partial
E^{I,n,\theta}_{c}}{\partial y}\|_{k}+\alpha
D_{2m}B_{7k}^{\prime}\|E^{I,n,\theta}_{c}\|_{k}+\\\ &\quad
D_{1m}D_{u1}B_{8k}\|\frac{\partial^{2}E^{I,n,\theta}_{c}}{\partial
x^{2}}\|_{k}+D_{2m}D_{u1}B_{8k}\|\frac{\partial^{2}E^{I,n,\theta}_{c}}{\partial
y^{2}}\|_{k}+\bar{D}_{u1}\bar{D}_{u2}B_{8k}\|\frac{\partial
E^{I,n,\theta}_{c}}{\partial x}\|_{k}+D_{u1}\bar{D}_{u2}\\\ &\quad
B_{8k}\|\frac{\partial E^{I,n,\theta}_{c}}{\partial y}\|_{k}+\alpha
D_{u1}B_{8k}\|E^{I,n,\theta}_{c}\|_{k}+D_{1m}D_{u2}B_{8k}^{\prime}\|\frac{\partial^{2}E^{I,n,\theta}_{c}}{\partial
x^{2}}\|_{k}+D_{2m}D_{u2}B_{8k}^{\prime}\|\frac{\partial^{2}E^{I,n,\theta}_{c}}{\partial
y^{2}}\|_{k}\\\ &\quad+\bar{D}_{u1}D_{u2}B_{8k}^{\prime}\|\frac{\partial
E^{I,n,\theta}_{c}}{\partial
x}\|_{k}+D_{u2}\bar{D}_{u2}B_{8k}^{\prime}\|\frac{\partial
E^{I,n,\theta}_{c}}{\partial y}\|_{k}+\alpha
D_{1m}B_{6k}\|\frac{\partial^{2}E^{I,n,\theta}_{c}}{\partial
x^{2}}\|_{k}+\alpha D_{2m}B_{6k}\\\
&\quad\|\frac{\partial^{2}E^{I,n,\theta}_{c}}{\partial
y^{2}}\|_{k}+\alpha\bar{D}_{u1}B_{6k}\|\frac{\partial
E^{I,n,\theta}_{c}}{\partial x}\|_{k}+\alpha\bar{D}_{u2}B_{6k}\|\frac{\partial
E^{I,n,\theta}_{c}}{\partial
y}\|_{k}+\alpha^{2}B_{6k}\|E^{I,n,\theta}_{c}\|_{k}+\alpha
D_{u2}B_{8k}^{\prime}\\\ &\quad\|E^{I,n,\theta}_{c}\|_{k})\\\
&\leq\mid\tau_{3}^{\prime}\mid\\{\sum_{k=1}^{n_{el}}(D_{1m}^{2}B_{7k}+2D_{1m}D_{2m}B_{7k}^{\prime}+D_{2l}^{2}B_{7k}^{\prime}+D_{2m}\bar{D}_{u1}B_{7k}^{\prime}+D_{1m}D_{u1}B_{8k}+\\\
&\quad
D_{2m}D_{u1}B_{8k}+D_{1m}D_{u2}B_{8k}^{\prime}+D_{2m}D_{u2}B_{8k}^{\prime}+\alpha
D_{1m}B_{6k}+\alpha D_{2m}B_{6k})\\}C(\frac{1+\theta}{2}\|c^{n+1}\|_{2}+\\\
&\quad\frac{1-\theta}{2}\|c^{n}\|_{2})+\mid\tau_{3}^{\prime}\mid\\{\sum_{k=1}^{n_{el}}(D_{1m}\bar{D}_{u1}B_{7k}+D_{1m}\bar{D}_{u2}B_{7k}+D_{2m}\bar{D}_{u2}B_{7k}^{\prime}+\bar{D}_{u1}\bar{D}_{u2}B_{8k}+\\\
&\quad
D_{u1}\bar{D}_{u2}B_{8k}+D_{u2}\bar{D}_{u1}B_{8k}^{\prime}+D_{u2}\bar{D}_{u2}B_{8k}^{\prime}+\alpha\bar{D}_{u1}B_{6k}+\alpha\bar{D}_{u2}B_{6k})\\}Ch(\frac{1+\theta}{2}\|c^{n+1}\|_{1}+\\\
&\quad\frac{1-\theta}{2}\|c^{n}\|_{1})+\mid\tau_{3}^{\prime}\mid\\{\sum_{k=1}^{n_{el}}(\alpha
D_{1m}B_{7k}+\alpha D_{2m}B_{7k}^{\prime}+\alpha D_{u1}B_{8k}+\alpha
D_{u2}B_{8k}^{\prime}+\alpha^{2}B_{6k})\\}Ch^{2}\\\
&\quad(\frac{1+\theta}{2}\|c^{n+1}\|+\frac{1-\theta}{2}\|c^{n}\|)\\\
&\leq\frac{\mid\tau_{3}\mid
T}{(T_{0}-C_{\tau_{3}})}C\\{(\sum_{k=1}^{n_{el}}D_{B_{2k}})(\frac{1+\theta}{2}\|c^{n+1}\|_{2}+\frac{1-\theta}{2}\|c^{n}\|_{2})+h(\sum_{k=1}^{n_{el}}D_{B_{3k}})(\frac{1+\theta}{2}\|c^{n+1}\|_{1}\\\
&\quad+\frac{1-\theta}{2}\|c^{n}\|_{1})+h^{2}(\sum_{k=1}^{n_{el}}D_{B_{4k}})(\frac{1+\theta}{2}\|c^{n+1}\|+\frac{1-\theta}{2}\|c^{n}\|)\\}\end{split}$
(58)

where $D_{B_{2k}},D_{B_{3k}}$ and $D_{B_{4k}}$ are denoting respectively the
summations in which the notations $\bar{D}_{u1},\bar{D}_{u2}$ are the maximum
of the functions $(\frac{\partial D_{1}}{\partial x}-u_{1}),(\frac{\partial
D_{2}}{\partial y}-u_{2})$ respectively over $\Omega$. The next term is
similar to the previous one. Therefore the simplification will be same as
above. Hence skipping the calculations we directly put the result as follows:

$\begin{split}P_{3}&=\sum_{k=1}^{n_{el}}\tau_{3}^{\prime}(\bigtriangledown\cdot\tilde{\bigtriangledown}E^{A,n,\theta}_{c}-\textbf{u}\cdot\bigtriangledown
E^{A,n,\theta}_{c}-\alpha
E^{A,n,\theta}_{c},\bigtriangledown\cdot\tilde{\bigtriangledown}E_{c}^{A,n,\theta}+\textbf{u}\cdot\bigtriangledown
E_{c}^{A,n,\theta}\\\ &\quad-\alpha E_{c}^{A,n,\theta})_{\Omega_{k}}\\\
\end{split}$ (59)
$\begin{split}&\leq\mid\tau_{3}^{\prime}\mid\\{\sum_{k=1}^{n_{el}}(D_{1m}^{2}B_{7k}^{2}+2D_{1m}D_{2m}B_{7k}B_{7k}^{\prime}+D_{1m}\bar{D}_{u1}B_{8k}B_{7k}+D_{1m}\bar{D}_{u2}B_{8k}^{\prime}B_{7k}+\\\
&\quad\alpha D_{1m}B_{7k}B_{6k}+D_{2m}^{2}B_{7k}^{\prime
2}+D_{2m}\bar{D}_{u1}B_{8k}B_{7k}^{\prime}+D_{2m}\bar{D}_{u2}B_{8k}^{\prime}B_{7k}^{\prime}+\alpha
D_{2m}B_{6k}B_{7k}^{\prime}+\\\ &\quad
D_{1m}D_{u1}B_{8k}B_{7k}+D_{2m}D_{u1}B_{8k}B_{7k}^{\prime}+D_{u1}\bar{D}_{u1}B_{7k}^{2}+D_{u1}\bar{D}_{u2}B_{8k}B_{8k}^{\prime}+\alpha
D_{u1}B_{6k}\\\ &\quad
B_{8k}+D_{1m}D_{u2}B_{8k}^{\prime}B_{7k}+D_{2m}D_{u2}B_{8k}^{\prime}B_{7k}^{\prime}+D_{u2}\bar{D}_{u1}B_{8k}B_{8k}^{\prime}+D_{u2}\bar{D}_{u2}B_{8k}^{\prime
2}+\alpha\\\
&\quad(D_{u2}B_{8k}^{\prime}B_{6k}+D_{1m}B_{6k}B_{7k}+D_{2m}B_{6k}B_{7k}^{\prime}+\bar{D}_{u1}B_{8k}B_{6k}+\bar{D}_{u2}B_{8k}^{\prime}B_{6k}+\alpha
B_{6k}^{2}))\\}\\\ &\leq\frac{\mid\tau_{3}\mid
T}{(T_{0}-C_{\tau_{3}})}\sum_{k=1}^{n_{el}}D_{B_{5k}}\end{split}$ (60)

where $D_{B_{5k}}$ is a notation denoting the big sum of the constants.
Now combining all the bounds obtained for $P_{1},P_{2},P_{3}$ and putting them
into the expression of $I^{4}_{3}$ we will have

$\begin{split}-I^{4}_{3}&\leq\frac{C_{\tau_{3}}T}{dt(T_{0}-C_{\tau_{3}})}\\{\sum_{k=1}^{n_{el}}(\alpha
B_{6k}+D_{B_{1k}})\\}(\|E^{A,n+1}_{c}\|^{2}-\|E^{A,n}_{c}\|^{2})+\frac{C_{\tau_{3}}Ch^{2}}{(T_{0}-C_{\tau_{3}})}\\{\sum_{k=1}^{n_{el}}D_{B_{1k}}\\}\\\
&\quad(\|c^{n+1}\|_{2}+\|c^{n}\|_{2})+\frac{\mid\tau_{3}\mid
TC}{(T_{0}-C_{\tau_{3}})}\\{(\sum_{k=1}^{n_{el}}D_{B_{2k}})(\frac{1+\theta}{2}\|c^{n+1}\|_{2}+\frac{1-\theta}{2}\|c^{n}\|_{2})+h\\\
&\quad(\sum_{k=1}^{n_{el}}D_{B_{3k}})(\frac{1+\theta}{2}\|c^{n+1}\|_{1}+\frac{1-\theta}{2}\mid
c^{n}\mid_{1})+h^{2}(\sum_{k=1}^{n_{el}}D_{B_{4k}})(\frac{1+\theta}{2}\|c^{n+1}\|+\frac{1-\theta}{2}\|c^{n}\|)\\}\\\
&\quad+\frac{\mid\tau_{3}\mid
T}{(T_{0}-C_{\tau_{3}})}\sum_{k=1}^{n_{el}}D_{B_{5k}}\end{split}$ (61)

Finally here the process of finding bound for each term of $I_{3}$ is
completed. Now we will focus on finding bounds for the terms of $I_{4}$.
Before going to derivation let us see the term $d_{4}$ explicitly.

$\begin{split}d_{4}&=\sum_{i=1}^{n+1}(\frac{1}{dt}\tau_{3}^{\prime})^{i}(\partial_{t}(c^{n}-c_{h}^{n})-\bigtriangledown\cdot\tilde{\bigtriangledown}(c^{n,\theta}-c_{h}^{n,\theta})+\textbf{u}^{n}\cdot\bigtriangledown(c^{n,\theta}-c_{h}^{n,\theta})+\alpha(c^{n,\theta}-c_{h}^{n,\theta}))\\\
&\leq\sum_{i=1}^{\infty}(\frac{1}{dt}\tau_{3}^{\prime})^{i}(\partial_{t}(E^{I,n}_{c}+E_{c}^{A,n})-\bigtriangledown\cdot\tilde{\bigtriangledown}(E^{I,n,\theta}_{c}+E_{c}^{A,n,\theta})+\textbf{u}^{n}\cdot\bigtriangledown(E^{I,n,\theta}_{c}+E_{c}^{A,n,\theta})\\\
&\quad+\alpha(E^{I,n,\theta}_{c}+E_{c}^{A,n,\theta}))\\\
&=\frac{\tau_{3}^{\prime}}{(dt-\tau_{3}^{\prime})}(\partial_{t}E^{I,n}_{c}+\partial_{t}E^{A,n}_{c})-\frac{\tau_{3}^{\prime}}{(dt-\tau_{3}^{\prime})}(\bigtriangledown\cdot\tilde{\bigtriangledown}E^{I,n,\theta}_{c}-\textbf{u}^{n}\cdot\bigtriangledown
E^{I,n,\theta}_{c}-\alpha E^{I,n,\theta}_{c})\\\
&\quad-\frac{\tau_{3}^{\prime}}{(dt-\tau_{3}^{\prime})}(\bigtriangledown\cdot\tilde{\bigtriangledown}E^{A,n,\theta}_{c}-\textbf{u}^{n}\cdot\bigtriangledown
E^{A,n,\theta}_{c}-\alpha E^{A,n,\theta}_{c})\end{split}$ (62)

Since $\frac{\tau_{3}}{dt+\tau_{3}}<1$, which implies
$\frac{\tau_{3}^{\prime}}{dt}<1$ and therefore the series
$\sum_{i=1}^{\infty}(\frac{1}{dt}\tau_{3}^{\prime})^{i}$ converges to
$\frac{\tau_{3}^{\prime}}{(dt-\tau_{3}^{\prime})}$

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

This completes finding the bounds for $I_{4}$.
Now we will find bounds for $I_{5}$ and $I_{6}$ in similar manner as many
terms of $I_{5},I_{6}$ coincide with the terms of $I_{3}$ and $I_{4}$.

$\begin{split}-I_{5}&=(1-\tau_{3}^{-1}\tau_{3}^{\prime})\sum_{k=1}^{n_{el}}(\partial_{t}E^{A,n}_{c},E_{c}^{A,n,\theta})_{\Omega_{k}}+(1-\tau_{3}^{-1}\tau_{3}^{\prime})\sum_{k=1}^{n_{el}}(\textbf{u}^{n}\cdot\bigtriangledown
E^{I,n,\theta}_{c}-\bigtriangledown\cdot\tilde{\bigtriangledown}E^{I,n,\theta}_{c},\\\
&\quad
E_{c}^{A,n,\theta})_{\Omega_{k}}+(1-\tau_{3}^{-1}\tau_{3}^{\prime})\sum_{k=1}^{n_{el}}(\textbf{u}^{n}\cdot\bigtriangledown
E^{A,n,\theta}_{c}-\bigtriangledown\cdot\tilde{\bigtriangledown}E^{A,n,\theta}_{c}+\alpha
E^{A,n,\theta}_{c},E^{A,n,\theta}_{c})_{\Omega_{k}}\\\
&=Q_{1}+Q_{2}+Q_{3}\end{split}$ (64)

where

$\begin{split}Q_{1}&=(1-\tau_{3}^{-1}\tau_{3}^{\prime})\sum_{k=1}^{n_{el}}(\partial_{t}E_{c}^{A,n},E_{c}^{A,n,\theta})_{\Omega_{k}}\\\
&=\frac{\tau_{3}}{dt+\tau_{3}}\sum_{k=1}^{n_{el}}(\partial_{t}E_{c}^{A,n},E_{c}^{A,n,\theta})_{\Omega_{k}}\\\
&\leq\frac{C_{\tau_{3}}}{dt(T_{0}-C_{\tau_{3}})}(\sum_{k=1}^{n_{el}}B_{6k})(\|E_{c}^{A,n+1}\|^{2}-\|E_{c}^{A,n}\|^{2})\end{split}$
(65)
$\begin{split}Q_{2}&=(1-\tau_{3}^{-1}\tau_{3}^{\prime})\sum_{k=1}^{n_{el}}(\textbf{u}^{n}\cdot\bigtriangledown
E^{I,n,\theta}_{c}-\bigtriangledown\cdot\tilde{\bigtriangledown}E^{I,n,\theta}_{c},E_{c}^{A,n,\theta})_{\Omega_{k}}\\\
&\leq\frac{C_{\tau_{3}}}{(T_{0}-C_{\tau_{3}})}\\{\sum_{k=1}^{n_{el}}(\bar{D}_{u1}+\bar{D}_{u2})B_{6k}\\}\|E^{I,n,\theta}_{c}\|_{1}\\\
&\leq\frac{\mid\tau_{3}\mid
Ch}{(T_{0}-C_{\tau_{3}})}\\{\sum_{k=1}^{n_{el}}(\bar{D}_{u1}+\bar{D}_{u2})B_{6k}\\}(\frac{1+\theta}{2}\|c^{n+1}\|_{1}+\frac{1-\theta}{2}\|c^{n}\|_{1})\end{split}$
(66)

and

$\begin{split}Q_{3}&=(1-\tau_{3}^{-1}\tau_{3}^{\prime})\sum_{k=1}^{n_{el}}(\textbf{u}^{n}\cdot\bigtriangledown
E^{A,n,\theta}_{c}-\bigtriangledown\cdot\tilde{\bigtriangledown}E^{A,n,\theta}_{c}+\alpha
E^{A,n,\theta}_{c},E^{A,n,\theta}_{c})_{\Omega_{k}}\\\
&\leq\frac{C_{\tau_{3}}}{(T_{0}-C_{\tau_{3}})}\sum_{k=1}^{n_{el}}(D_{1m}B_{7k}+D_{2m}B_{7k}^{\prime}+\bar{D}_{u1}B_{8k}+\bar{D}_{u2}B_{8k}^{\prime}+\alpha
B_{6k})B_{6k}\end{split}$ (67)

Combining all these results we will have

$\begin{split}-I_{5}&\leq\frac{C_{\tau_{3}}}{(T_{0}-C_{\tau_{3}})}(\sum_{k=1}^{n_{el}}\frac{B_{6k}}{dt})(\|E_{c}^{A,n+1}\|^{2}-\|E_{c}^{A,n}\|^{2})+\\\
&\quad\frac{\mid\tau_{3}\mid}{(T_{0}-C_{\tau_{3}})}\\{(\sum_{k=1}^{n_{el}}(\bar{D}_{u1}+\bar{D}_{u2})B_{6k})Ch(\frac{1+\theta}{2}\|c^{n+1}\|_{1}+\frac{1-\theta}{2}\|c^{n}\|_{1})\\\
&\quad+\sum_{k=1}^{n_{el}}(D_{1m}B_{7k}+D_{2m}B_{7k}^{\prime}+\bar{D}_{u1}B_{8k}+\bar{D}_{u2}B_{8k}^{\prime}+\alpha
B_{6k})B_{6k}\\}\end{split}$ (68)

This completes finding the bound for each term of $I_{5}$. Now we focus on
deriving bounds of $I_{6}$

$\begin{split}-I_{6}&=-\sum_{k=1}^{n_{el}}\tau_{3}^{-1}\tau_{3}^{\prime}(d_{4},E^{A,n,\theta}_{c})_{\Omega_{k}}\\\
&\leq\frac{\mid\tau_{3}^{-1}\mid\tau_{3}^{\prime
2}}{(dt-\mid\tau_{3}^{\prime}\mid)}\\{(\partial_{t}E^{I,n}_{c}+\partial_{t}E^{A,n}_{c},E^{A,n,\theta}_{c})_{\Omega_{k}}-(\bigtriangledown\cdot\tilde{\bigtriangledown}E^{I,n,\theta}_{c}-\textbf{u}\cdot\bigtriangledown
E^{I,n,\theta}_{c}-\\\ &\quad\alpha
E^{I,n,\theta}_{c},E^{A,n,\theta}_{c})_{\Omega_{k}}-(\bigtriangledown\cdot\tilde{\bigtriangledown}E^{A,n,\theta}_{c}-\textbf{u}\cdot\bigtriangledown
E^{A,n,\theta}_{c}-\alpha
E^{A,n,\theta}_{c},E^{A,n,\theta}_{c})_{\Omega_{k}}\\}\\\
&\leq\frac{\tau_{3}}{(dt+\tau_{3})}\\{(\sum_{k=1}^{n_{el}}\frac{B_{6k}}{dt})(\|E^{A,n+1}_{c}\|^{2}-\|E^{A,n}_{c}\|^{2})+(\sum_{k=1}^{n_{el}}(D_{1m}+D_{2m})B_{6k})\\\
&\quad\mid
E^{I,n,\theta}_{c}\mid_{2}+(\sum_{k=1}^{n_{el}}(\bar{D}_{u1}+\bar{D}_{u2})B_{6k})\mid
E^{I,n,\theta}_{c}\mid_{1}+\sum_{k=1}^{n_{el}}(D_{1m}B_{7k}+D_{2m}B_{7k}^{\prime}+\\\
&\quad\bar{D}_{u1}B_{8k}+\bar{D}_{u2}B_{8k}^{\prime}+\alpha
B_{6k})B_{6k}\\}\\\
&\leq\frac{C_{\tau_{3}}}{(T_{0}-C_{\tau_{3}})}\\{(\sum_{k=1}^{n_{el}}\frac{B_{6k}}{dt})(\|E^{A,n+1}_{c}\|^{2}-\|E^{A,n}_{c}\|^{2})+\frac{\mid\tau_{3}\mid}{(T_{0}-C_{\tau_{3}})}(\sum_{k=1}^{n_{el}}C(D_{1m}+\\\
&\quad
D_{2m})B_{6k})(\frac{1+\theta}{2}\|c^{n+1}\|_{2}+\frac{1-\theta}{2}\|c^{n}\|_{2})+Ch(\sum_{k=1}^{n_{el}}(\bar{D}_{u1}+\bar{D}_{u2})B_{6k})(\frac{1+\theta}{2}\|c^{n+1}\|_{1}\\\
&\quad+\frac{1-\theta}{2}\|c^{n}\|_{1})+\sum_{k=1}^{n_{el}}(D_{1m}B_{7k}+D_{2m}B_{7k}^{\prime}+\bar{D}_{u1}B_{8k}+\bar{D}_{u2}B_{8k}^{\prime}+\alpha
B_{6k})B_{6k}\\}\end{split}$ (69)

Again

$(TE^{n,\theta},E^{A,n,\theta}_{c})\leq\frac{\epsilon_{6}}{2}\|TE^{n,\theta}\|^{2}+\frac{1}{2\epsilon_{6}}\|E^{A,n,\theta}_{c}\|$
(70)

Finally we have completed finding bounds for each of the terms in the right
hand side of (37). Now we explain the further proceeding in language as
follows:
First we put all the bounds, obtained for each of the terms in the right hand
side of (37). Then we take out few common terms in the left hand side and
consequently we have left 3 types of terms in the right hand side. One type
will be few constant terms multiplied by $h^{2}$, other type will be another
few constant terms multiplied by h and the remaining constant terms will be
free of h. 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 (37)
as follows:

$\\{1-\frac{2C_{\tau_{3}}(T+C_{\tau_{3}})}{T_{0}-C_{\tau_{3}}}\sum_{k=1}^{n_{el}}(\alpha
B_{6k}+D_{B_{1k}})-\frac{4C_{\tau_{3}}}{T_{0}-C_{\tau_{3}}}(\sum_{k=1}^{n_{el}}B_{6k})\\}\sum_{n=0}^{N-1}\int_{t^{n}}^{t^{n+1}}(\|E^{A,n+1}_{c}\|^{2}-\|E^{A,n}_{c}\|^{2})\\\
\quad+(\mu_{l}-\frac{1}{\epsilon_{1}}-\frac{\mu_{u}}{\epsilon_{4}})\sum_{n=0}^{N-1}\int_{t^{n}}^{t^{n+1}}(\|\frac{\partial
E^{A,n,\theta}_{u1}}{\partial x}\|^{2}+\|\frac{\partial
E^{A,n,\theta}_{u2}}{\partial
y}\|^{2})dt+(\mu_{l}-\frac{\mu_{u}}{\epsilon_{4}})\sum_{n=0}^{N-1}\int_{t^{n}}^{t^{n+1}}(\|\frac{\partial
E^{A,n,\theta}_{u1}}{\partial y}\|^{2}+\\\ \quad\|\frac{\partial
E^{A,n,\theta}_{u2}}{\partial
x}\|^{2})dt+(D_{l}-\frac{D_{m}}{\epsilon_{2}}-\frac{C_{1}^{n}}{\epsilon_{3}})\sum_{n=0}^{N-1}\int_{t^{n}}^{t^{n+1}}\|\frac{\partial
E^{A,n,\theta}_{c}}{\partial
x}\|^{2}dt+(D_{l}-\frac{D_{m}}{\epsilon_{2}}-\frac{C_{2}^{n}}{\epsilon_{3}})\sum_{n=0}^{N-1}\int_{t^{n}}^{t^{n+1}}\|\frac{\partial
E^{A,n,\theta}_{c}}{\partial y}\|^{2}dt\\\
\quad+\sigma\sum_{n=0}^{N-1}\int_{t^{n}}^{t^{n+1}}(\|E^{A,n,\theta}_{u1}\|^{2}+\|E^{A,n,\theta}_{u2}\|^{2})dt+(\alpha-2\epsilon_{3}(C_{1}^{n}+C_{2}^{n})-\frac{1}{\epsilon_{6}})\sum_{n=0}^{N-1}\int_{t^{n}}^{t^{n+1}}\|E^{A,n,\theta}_{c}\|^{2}dt\\\
\leq
2Ch^{2}\sum_{n=0}^{N-1}\int_{t^{n}}^{t^{n+1}}[C\mu_{u}\epsilon_{4}\sum_{i=1}^{2}(\frac{1+\theta}{2}\|u_{i}^{n+1}\|_{2}+\frac{1-\theta}{2}\|u_{i}^{n}\|_{2})^{2}+C\epsilon_{1}(\frac{1+\theta}{2}\|p^{n+1}\|_{1}+\frac{1-\theta}{2}\|p^{n}\|_{1})^{2}+\\\
\quad
C(\frac{D_{m}\epsilon_{2}}{2}+\frac{C_{1}^{n}+C_{2}^{n}}{2\epsilon_{3}})(\frac{1+\theta}{2}\|c^{n+1}\|_{2}+\frac{1-\theta}{2}\|c^{n}\|_{2})^{2}+\frac{\mid\tau_{3}\mid(T_{0}+\mid\tau_{3}\mid)}{T_{0}(T_{0}-C_{\tau_{3}})}(\sum_{k=1}^{n_{el}}D_{B_{1k}})(\|c^{n+1}\|_{2}+\\\
\quad\|c^{n}\|_{2})+\mid\tau_{1}\mid(\sum_{k=1}^{n_{el}}\sigma(\mu_{u}B_{2k}+\sigma
B_{1k}+B_{3k}))(\frac{1+\theta}{2}\|u_{1}^{n+1}\|_{2}+\frac{1-\theta}{2}\|u_{1}^{n}\|_{2})+\mid\tau_{1}\mid(\sum_{k=1}^{n_{el}}\sigma(\mu_{u}B_{2k}^{\prime}+\\\
\quad\sigma
B_{1k}^{\prime}+B_{3k}^{\prime}))(\frac{1+\theta}{2}\|u_{2}^{n+1}\|_{2}+\frac{1-\theta}{2}\|u_{2}^{n}\|_{2})+(\sum_{k=1}^{n_{el}}D_{B_{4k}})\frac{\mid\tau_{3}\mid}{T_{0}}(\frac{1+\theta}{2}\|c^{n+1}\|+\frac{1-\theta}{2}\|c^{n}\|)]dt\\\
\quad+2Ch\sum_{n=0}^{N-1}\int_{t^{n}}^{t^{n+1}}[\mid\tau_{1}\mid(\sum_{k=1}^{n_{el}}(\mu_{u}(B_{2k}+B_{2k}^{\prime})+\sigma(B_{1k}+B_{1k}^{\prime})+(B_{3k}+B_{3k}^{\prime}))(\frac{1+\theta}{2}\|p^{n+1}\|_{1}+\frac{1-\theta}{2}\|p^{n}\|_{1})\\\
\quad+\mid\tau_{2}\mid(\sum_{k=1}^{n_{el}}B_{4k})(\frac{1+\theta}{2}\|u_{1}^{n+1}\|_{1}+\frac{1-\theta}{2}\|u_{1}^{n}\|_{1})+\mid\tau_{2}\mid(\sum_{k=1}^{n_{el}}B_{5k}^{\prime})(\frac{1+\theta}{2}\|u_{2}^{n+1}\|_{1}+\frac{1-\theta}{2}\|u_{2}^{n}\|_{1})\\\
\quad+\\{(\frac{\mid\tau_{3}\mid
T}{(T_{0}-C_{\tau_{3}})}+\frac{\mid\tau_{3}\mid}{T_{0}})(\sum_{k=1}^{n_{el}}D_{B_{3k}})+\frac{4\mid\tau_{3}\mid}{T_{0}-C_{\tau_{3}}}(\sum_{k=1}^{n_{el}}(\bar{D}_{u1}+\bar{D}_{u2})B_{6k})\\}(\frac{1+\theta}{2}\|c^{n+1}\|_{1}+\frac{1-\theta}{2}\\\
\quad\|c^{n}\|_{1})]dt+2C\sum_{n=0}^{N-1}\int_{t^{n}}^{t^{n+1}}[\mid\tau_{1}\mid(\sum_{k=1}^{n_{el}}(\mu_{u}^{2}B_{2k}+\sigma\mu_{u}B_{1k}+\mu_{u}B_{3k}))(\frac{1+\theta}{2}\|u_{1}^{n+1}\|_{2}+\frac{1-\theta}{2}\|u_{1}^{n}\|_{2})+\\\
\quad\mid\tau_{1}\mid(\sum_{k=1}^{n_{el}}(\mu_{u}^{2}B_{2k}^{\prime}+\sigma\mu_{u}B_{1k}^{\prime}+\mu_{u}B_{3k}^{\prime}))(\frac{1+\theta}{2}\|u_{2}^{n+1}\|_{2}+\frac{1-\theta}{2}\|u_{2}^{n}\|_{2})+(\frac{\mid\tau_{3}\mid
T}{(T_{0}-C_{\tau_{3}})}+\frac{\mid\tau_{3}\mid}{T_{0}})\\\
\quad(\sum_{k=1}^{n_{el}}D_{B_{2k}})(\frac{1+\theta}{2}\|c^{n+1}\|_{1}+\frac{1-\theta}{2}\|c^{n}\|_{1})+\mid\tau_{1}\mid(\sum_{k=1}^{n_{el}}(M_{1k}+M_{2k}))+(\frac{\mid\tau_{3}\mid
T}{(T_{0}-C_{\tau_{3}})}+\frac{\mid\tau_{3}\mid}{T_{0}})\\\
\quad\\{(\sum_{k=1}^{n_{el}}D_{B_{5k}})+(\sum_{k=1}^{n_{el}}(D_{1m}+D_{2m})B_{6k})(\frac{1+\theta}{2}\|c^{n+1}\|_{2}+\frac{1-\theta}{2}\|c^{n}\|_{2})+2(\sum_{k=1}^{n_{el}}(D_{1m}B_{7k}+D_{2m}B_{7k}^{\prime}\\\
\quad+\bar{D}_{u1}B_{8k}+\bar{D}_{u2}B_{8k}^{\prime}+\alpha
B_{6k})B_{6k})\\}]dt+\epsilon_{6}\sum_{n=0}^{N-1}\int_{t^{n}}^{t^{n+1}}\|TE^{n,\theta}\|^{2}dt\hskip
28.45274pt$ (71)

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 $h$ 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 $\|u_{i}^{n}\|_{2}$ for $i=1,2$,
$\|p^{n}\|_{1}$ and $\|c^{n}\|_{2}$ are bounded for $n=0,1,2,...,N$. Now by
applying initial condition on $c$ we will have $\|E^{A,0}_{c}\|=0$.
After performing all these intermediate steps we will finally arrive at the
following expression since $\tau_{1}$ and $\tau_{3}$ are of order $h^{2}$:

$\|E^{A,N}_{c}\|^{2}+\sum_{n=0}^{N-1}\int_{t^{n}}^{t^{n+1}}(\|\frac{\partial
E^{A,n,\theta}_{u1}}{\partial x}\|^{2}+\|\frac{\partial
E^{A,n,\theta}_{u1}}{\partial y}\|^{2}+\|E^{A,n,\theta}_{u1}\|^{2})dt+\\\
\quad\sum_{n=0}^{N-1}\int_{t^{n}}^{t^{n+1}}(\|\frac{\partial
E^{A,n,\theta}_{u2}}{\partial x}\|^{2}+\|\frac{\partial
E^{A,n,\theta}_{u2}}{\partial y}\|^{2}+\|E^{A,n,\theta}_{u2}\|^{2})dt+\\\
\quad\sum_{n=0}^{N-1}\int_{t^{n}}^{t^{n+1}}(\|\frac{\partial
E^{A,n,\theta}_{c}}{\partial x}\|^{2}+\|\frac{\partial
E^{A,n,\theta}_{c}}{\partial y}\|^{2}+\|E^{A,n,\theta}_{c}\|^{2})dt\leq
C(T,\textbf{u},p,c)(h^{2}+h+dt^{2r})$ (72)

This implies

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

where

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

We have used the fact that $\sum_{n=0}^{N-1}\int_{t^{n}}^{t^{n+1}}Mdt\leq MT$
and the property of $TE$ given in (17). 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
Stokes-Darcy flow problem we have obtained

$\begin{split}a_{S}(\textbf{u}-\textbf{u}_{h},\textbf{v}_{h})-b(\textbf{v}_{h},p-p_{h})&=0\\\
b(\textbf{v}_{h},p-I_{h}p)+b(\textbf{v}_{h},I_{h}p-p_{h})&=a_{S}(E^{I}_{\textbf{u}},\textbf{v}_{h})+a_{S}(E^{A}_{\textbf{u}},\textbf{v}_{h})\end{split}$
(75)

Assuming the inclusion $\bigtriangledown\cdot V_{s}^{h}\subset Q_{s}^{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$
(76)

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}$
(77)

Now from (75)

$\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}}\\{a_{S}(E^{I,n,\theta}_{\textbf{u}},\textbf{v}_{h})+a_{S}(E^{A,n,\theta}_{\textbf{u}},\textbf{v}_{h})\\}dt\\\
&=\sum_{n=0}^{N-1}\int_{t^{n}}^{t^{n+1}}\\{\int_{\Omega}\mu(c^{n})\bigtriangledown
E^{I,n,\theta}_{u1}\cdot\bigtriangledown
v_{1h}+\int_{\Omega}\mu(c^{n})\bigtriangledown
E^{I,n,\theta}_{u2}\cdot\bigtriangledown v_{2h}+\\\
&\quad\sigma\int_{\Omega}(E^{I,n,\theta}_{u1}v_{1h}+E^{I,n,\theta}_{u2}v_{2h})+\int_{\Omega}\mu(c^{n})\bigtriangledown
E^{A,n,\theta}_{u1}\cdot\bigtriangledown v_{1h}+\\\
&\quad\int_{\Omega}\mu(c^{n})\bigtriangledown
E^{A,n,\theta}_{u2}\cdot\bigtriangledown
v_{2h}+\sigma\int_{\Omega}(E^{A,n,\theta}_{u1}v_{1h}+E^{A,n,\theta}_{u2}v_{2h})\\}dt\\\
&\leq(\mu_{u}+\sigma)\sum_{n=0}^{N-1}\int_{t^{n}}^{t^{n+1}}\\{(\|E^{I,n,\theta}_{u1}\|_{1}+\|E^{A,n,\theta}_{u1}\|_{1})\|v_{1h}\|_{1}+(\|E^{I,n,\theta}_{u2}\|_{1}+\\\
&\quad\|E^{A,n,\theta}_{u2}\|_{1})\|v_{2h}\|_{1}\\}dt\\\
&\leq(\mu_{u}+\sigma)\sum_{n=0}^{N-1}\int_{t^{n}}^{t^{n+1}}\\{(\|E^{I,n,\theta}_{u1}\|_{1}+\|E^{I,n,\theta}_{u2}\|_{1}+\|E^{A,n,\theta}_{u1}\|_{1}+\|E^{A,n,\theta}_{u2}\|_{1})\\\
&\quad(\|v_{1h}\|_{1}+\|v_{2h}\|_{1})\\}dt\\\
&\leq(\mu_{u}+\sigma)(\|E^{A}_{u1}\|_{L^{2}(H^{1})}+\|E^{A}_{u2}\|_{L^{2}(H^{1})})(\|v_{1h}\|_{1}+\|v_{2h}\|_{1})+\\\
&\quad(\mu_{u}+\sigma)\sum_{n=0}^{N-1}\int_{t^{n}}^{t^{n+1}}(\|E^{I,n,\theta}_{u1}\|_{1}+\|E^{I,n,\theta}_{u2}\|_{1})(\|v_{1h}\|_{1}+\|v_{2h}\|_{1})dt\\\
&\leq(\mu_{u}+\sigma)(\|E^{A}_{u1}\|_{L^{2}(H^{1})}^{2}+\|E^{A}_{u2}\|_{L^{2}(H^{1})}^{2}+\|E^{A}_{c}\|_{\tilde{\textbf{V}}}^{2})\|\textbf{v}_{h}\|_{1}+(\mu_{u}+\sigma)Ch\\\
&\quad\sum_{n=0}^{N-1}\int_{t^{n}}^{t^{n+1}}\\{\sum_{i=1}^{2}(\frac{1+\theta}{2}\mid
u_{i}^{n+1}\mid_{2}+\frac{1-\theta}{2}\mid
u_{i}^{n}\mid_{2})dt\\}\|\textbf{v}_{h}\|_{1}\\\ &\leq
C^{\prime}(T,\textbf{u},p,c)(h^{2}+h+dt^{2r})\|\textbf{v}_{h}\|_{1}\end{split}$
(78)

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

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

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

$\|E^{A}_{u1}\|^{2}_{L^{2}(H^{1})}+\|E^{A}_{u2}\|^{2}_{L^{2}(H^{1})}+\|E^{A}_{p}\|_{L^{2}(L^{2})}^{2}+\|E^{A}_{c}\|^{2}_{\tilde{\textbf{V}}}\leq\bar{C}(T,\textbf{u},p,c)(h^{2}+h+dt^{2r})$
(80)

where

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

This completes the proof. ∎

######  Theorem 2.

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

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

where C’ depends on T, u,p,c and

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

###### Proof.

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

$\|u_{1}-u_{1h}\|_{L^{2}(H^{1})}^{2}+\|u_{2}-u_{2h}\|_{L^{2}(H^{1})}^{2}+\|p-p_{h}\|_{L^{2}(L^{2})}^{2}+\|c-c_{h}\|^{2}_{\tilde{\textbf{V}}}\\\
=\|E^{I}_{u1}+E^{A}_{u1}\|_{L^{2}(H^{1})}^{2}+\|E^{I}_{u2}+E^{A}_{u2}\|_{L^{2}(H^{1})}^{2}+\|E^{I}_{p}+E^{A}_{p}\|_{L^{2}(L^{2})}^{2}+\|E^{I}_{c}+E^{A}_{c}\|_{\tilde{\textbf{V}}}^{2}\\\
\leq\bar{C}(\|E^{I}_{u1}\|_{L^{2}(H^{1})}^{2}+\|E^{I}_{u2}\|_{L^{2}(H^{1})}^{2}+\|E^{I}_{p}\|_{L^{2}(L^{2})}^{2}+\|E^{I}_{c}\|_{\tilde{\textbf{V}}}^{2}+\|E^{A}_{u1}\|_{L^{2}(H^{1})}^{2}+\|E^{A}_{u2}\|_{L^{2}(H^{1})}^{2}\\\
\quad+\|E^{A}_{p}\|_{L^{2}(L^{2})}^{2}+\|E^{A}_{c}\|_{\tilde{\textbf{V}}}^{2})\hskip
227.62204pt\\\ \leq C^{\prime\prime}(T,\textbf{u},p,c)(h^{2}+h+dt^{2r})\hskip
199.16928pt$ (84)

This completes apriori error estimation. ∎

###  3.4 Aposteriori error estimation

In this section we are going to derive residual based aposteriori error
estimation.
We have
$B(\textbf{V},\textbf{V})=a_{S}(\textbf{v},\textbf{v})+a_{T}(d,d)\geq\mu_{l}(\|v_{1}\|_{1}^{2}+\|v_{2}\|_{1}^{2})+D_{\alpha}\|d\|_{1}^{2}$
$\forall\textbf{V}\in\textbf{V}_{F}$
Now we substitute the errors $e_{u1},e_{u2},e_{c}$ into the relation we will
similarly have

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

By adding few terms in both sides the above equation becomes

$\underbrace{(\frac{\partial e_{c}}{\partial
t},e_{c})+\mu_{l}(\|e_{u1}\|_{1}^{2}+\|e_{u2}\|_{1}^{2})+\sigma\|e_{p}\|^{2}+D_{\alpha}\|e_{c}\|_{1}^{2}}_{\textit{LHS}}\\\
\leq\underbrace{(\frac{\partial e_{c}}{\partial
t},e_{c})+a_{S}(e_{\textbf{u}},e_{\textbf{u}})+a_{T}(e_{c},e_{c})+(e_{p},e_{p})+b(e_{\textbf{u}},e_{p})-b(e_{\textbf{u}},e_{p})}_{\textit{RHS}}$
(86)

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_{c}^{n+1}-e_{c}^{n}}{dt},e_{c}^{n,\theta})+\mu_{l}\sum_{i=1}^{2}(\|e_{ui}^{n,\theta}\|^{2}+\|\frac{\partial
e_{ui}^{n,\theta}}{\partial x}\|^{2}+\|\frac{\partial
e_{ui}^{n,\theta}}{\partial y}\|^{2})+\\\
\sigma\|e_{p}^{n,\theta}\|^{2}+D_{\alpha}(\|e_{c}^{n,\theta}\|^{2}+\|\frac{\partial
e_{c}^{n,\theta}}{\partial x}\|^{2}+\|\frac{\partial
e_{c}^{n,\theta}}{\partial y}\|^{2})\hskip 56.9055pt$ (87)

Using the same argument done in (33) we have

$\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}$ (88)

Hence

$\frac{1}{2dt}(\|e_{c}^{n+1}\|^{2}-\|e_{c}^{n}\|^{2})+\mu_{l}\sum_{i=1}^{2}(\|e_{ui}^{n,\theta}\|^{2}+\|\frac{\partial
e_{ui}^{n,\theta}}{\partial x}\|^{2}+\|\frac{\partial
e_{ui}^{n,\theta}}{\partial y}\|^{2})+\\\
\sigma\|e_{p}^{n,\theta}\|^{2}+D_{\alpha}(\|e_{c}^{n,\theta}\|^{2}+\|\frac{\partial
e_{c}^{n,\theta}}{\partial x}\|^{2}+\|\frac{\partial
e_{c}^{n,\theta}}{\partial y}\|^{2})\leq LHS\leq RHS$ (89)

Now our aim is to find upper bound for $RHS$ through dividing it into two
broad parts by splitting errors in each of the terms as follows:

$RHS=\\{(\frac{e_{c}^{n+1}-e_{c}^{n}}{dt},E^{I,n,\theta}_{c})+a_{S}(e^{n,\theta}_{\textbf{u}},E^{I,n,\theta}_{\textbf{u}})+a_{T}(e^{n,\theta}_{c},E^{I,n,\theta}_{c})\\\
-b(E^{I,n,\theta}_{\textbf{u}},e^{n,\theta}_{p})+b(e^{n,\theta}_{\textbf{u}},E^{I,n,\theta}_{p})\\}+\\{(\frac{e^{n+1}_{c}-e_{c}^{n}}{dt},E^{A,n,\theta}_{c})\\\
+a_{S}(e^{n,\theta}_{\textbf{u}},E^{A,n,\theta}_{\textbf{u}})+a_{T}(e^{n,\theta}_{c},E^{A,n,\theta}_{c})-b(E^{A,n,\theta}_{\textbf{u}},e^{n,\theta}_{p})+b(e^{n,\theta}_{\textbf{u}},E^{A,n,\theta}_{p})\\}$
(90)

In the expression of $RHS$ the first under brace part is first part and second
one is second part. Before proceeding further let us introduce the residuals
corresponding to each equations

$\textbf{R}^{h}=\begin{bmatrix}\textbf{f}_{1}-(-\mu(c)\Delta\textbf{u}_{h}+\sigma\textbf{u}_{h}+\bigtriangledown
p_{h})\\\ f_{2}-\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}$

This column vector $\textbf{R}^{h}$ has four components
$R_{1}^{h},R_{2}^{h},R_{3}^{h}$ and $R_{4}^{h}$ denoting four rows
respectively. Let us start finding bound for the first part as follows: for
all $\textbf{v}=(v_{1},v_{2})\in V_{s}\times V_{s}$

$\int_{\Omega}\mu(c^{n})\bigtriangledown
e_{\textbf{u}}^{n,\theta}:\bigtriangledown\textbf{v}+\int_{\Omega}\sigma
e_{\textbf{u}}^{n,\theta}\cdot\textbf{v}-\int_{\Omega}(\bigtriangledown\cdot\textbf{v})e_{p}^{n,\theta}\\\
=\\{\int_{\Omega}\mu(c^{n})\bigtriangledown\textbf{u}^{n,\theta}:\bigtriangledown\textbf{v}+\int_{\Omega}\sigma\textbf{u}^{n,\theta}\cdot\textbf{v}-\int_{\Omega}(\bigtriangledown\cdot\textbf{v})p^{n,\theta}\\}-\\{\int_{\Omega}\mu(c^{n})\bigtriangledown\textbf{u}_{h}^{n,\theta}:\bigtriangledown\textbf{v}+\hskip
5.69054pt\\\
\int_{\Omega}\sigma\textbf{u}_{h}^{n,\theta}\cdot\textbf{v}-\int_{\Omega}(\bigtriangledown\cdot\textbf{v})p_{h}^{n,\theta}\\}\hskip
142.26378pt\\\
=\int_{\Omega}(-\mu(c^{n})\Delta\textbf{u}^{n,\theta}+\sigma\textbf{u}^{n,\theta}+\bigtriangledown
p^{n,\theta})\cdot\textbf{v}-\int_{\Omega}(-\mu(c^{n})\Delta\textbf{u}_{h}^{n,\theta}+\sigma\textbf{u}_{h}^{n,\theta}+\bigtriangledown
p_{h}^{n,\theta})\cdot\textbf{v}\hskip 28.45274pt\\\
=(R_{1}^{h,n,\theta},v_{1})+(R_{2}^{h,n,\theta},v_{2})\hskip 227.62204pt$ (91)

Similarly $\int_{\Omega}(\bigtriangledown\cdot
e_{\textbf{u}}^{n,\theta})q=\int_{\Omega}R_{3}^{h,n,\theta}q$ $\forall q\in
Q_{s}$

$\int_{\Omega}(\frac{e_{c}^{n+1}-e_{c}^{n}}{dt}d+\tilde{\bigtriangledown}e_{c}^{n,\theta}\cdot\bigtriangledown
d+d\textbf{u}\cdot\bigtriangledown e_{c}^{n,\theta}+\alpha
e_{c}^{n,\theta}d)=\int_{\Omega}R_{4}^{h,n,\theta}d\hskip 8.53581pt\forall
d\in V_{s}\hskip 85.35826pt$ (92)

Now substituting $v_{1},v_{2},q,d$ in the above expressions by
$E^{I,n,\theta}_{u1},E^{I,n,\theta}_{u2},E^{I,n,\theta}_{p},E^{I,n,\theta}_{c}$
respectively, we will have the first part of the $RHS$ as,

$\begin{split}First\hskip 2.84526ptpart\hskip 2.84526ptof\hskip
2.84526ptRHS&=\int_{\Omega}\\{R_{1}^{h,n,\theta}E^{I,n,\theta}_{u1}+R_{2}^{h,n,\theta}E^{I,n,\theta}_{u2}+R_{3}^{h,n,\theta}E^{I,n,\theta}_{p}+R_{4}^{h,n,\theta}E^{I,n,\theta}_{c}\\}\\\
&\leq\|R_{1}^{h,n,\theta}\|\|E^{I,n,\theta}_{u1}\|+\|R_{2}^{h,n,\theta}\|\|E^{I,n,\theta}_{u2}\|+\|R_{3}^{h,n,\theta}\|\|E^{I,n,\theta}_{p}\|+\\\
&\quad\|R_{4}^{h,n,\theta}\|\|E^{I,n,\theta}_{c}\|\hskip 28.45274pt(by\hskip
2.84526ptCauchy-Schwarz\hskip 2.84526ptinequality)\\\
&\leq\|R_{1}^{h,n,\theta}\|h\|u_{1}^{n,\theta}\|_{1}+\|R_{2}^{h,n,\theta}\|h\|u_{2}^{n,\theta}\|_{1}+\|R_{3}^{h,n,\theta}\|\|p^{n,\theta}\|+\\\
&\quad\|R_{4}^{h,n,\theta}\|h\|c^{n,\theta}\|_{1}\\\ &\leq
C_{2}(\|R_{1}^{h,n,\theta}\|h\|e_{u1}^{n,\theta}\|_{1}+\|R_{2}^{h,n,\theta}\|h\|e_{u2}^{n,\theta}\|_{1}+\\\
&\quad\|R_{3}^{h,n,\theta}\|\|e_{p}^{n,\theta}\|+\|R_{4}^{h,n,\theta}\|h\|e_{c}^{n,\theta}\|_{1})\\\
&\leq
C_{2}\frac{h^{2}}{2\epsilon_{1}}(\|R_{1}^{h,n,\theta}\|^{2}+\|R_{2}^{h,n,\theta}\|^{2}+\|R_{3}^{h,n,\theta}\|^{2}+\|R_{4}^{h,n,\theta}\|^{2})+\\\
&\quad\frac{\epsilon_{1}}{2}(\|e_{u1}^{n,\theta}\|_{1}^{2}+\|e_{u2}^{n,\theta}\|_{1}^{2}+\|e_{p}^{n,\theta}\|^{2}+\|e_{c}^{n,\theta}\|_{1}^{2})\hskip
8.53581pt(by\hskip 2.84526ptYoung^{\prime}s\hskip 2.84526ptinequality)\\\
\end{split}$ (93)

This completes finding bound for first part of $RHS$. Now we are going to
estimate remaining second part of $RHS$. For that we will use 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}$

$\int_{\Omega}\frac{e_{c}^{n+1}-e_{c}^{n}}{dt}d_{h}+\int_{\Omega}\mu(c^{n})\bigtriangledown
e_{\textbf{u}}^{n,\theta}:\bigtriangledown\textbf{v}_{h}+\int_{\Omega}\sigma
e_{\textbf{u}}^{n,\theta}\cdot\textbf{v}_{h}-\int_{\Omega}(\bigtriangledown\cdot\textbf{v}_{h})e_{p}^{n,\theta}+\int_{\Omega}(\bigtriangledown\cdot
e_{\textbf{u}}^{n,\theta})q_{h}\\\
+\int_{\Omega}\tilde{\bigtriangledown}e_{c}^{n,\theta}\cdot\bigtriangledown
d_{h}+\int_{\Omega}d_{h}\textbf{u}\cdot\bigtriangledown
e_{c}^{n,\theta}+\int_{\Omega}\alpha e_{c}^{n,\theta}d_{h}\\\
=\sum_{k=1}^{n_{el}}\\{(\tau_{k}^{\prime}(\textbf{R}^{h,n,\theta}+\textbf{d}),-\mathcal{L}^{*}\textbf{V}_{h})_{\Omega_{k}}-((I-\tau_{k}^{-1}\tau_{k})\textbf{R}^{h,n,\theta},\textbf{V}_{h})_{\Omega_{k}}+(\tau_{k}^{-1}\tau_{k}\textbf{d},\textbf{V}_{h})_{\Omega_{k}}\\}\\\
+(TE^{n,\theta},d_{h})\hskip 227.62204pt\\\
=\sum_{k=1}^{n_{el}}\\{\tau_{1}^{\prime}(R_{1}^{h,n,\theta},\mu(c)\Delta
v_{1h}-\sigma v_{1h}+\frac{\partial q_{h}}{\partial
x})_{k}+\tau_{1}^{\prime}(R_{2}^{h,n,\theta},\mu(c)\Delta v_{2h}-\sigma
v_{2h}+\frac{\partial q_{h}}{\partial y})_{k}\\\
\quad+\tau_{2}^{\prime}(R_{3}^{h,n,\theta},\bigtriangledown\cdot\textbf{v}_{h})_{k}+\tau_{3}^{\prime}(R_{4}^{h,n,\theta}+d_{4},\bigtriangledown\cdot\tilde{\bigtriangledown}d_{h}+\textbf{u}\cdot\bigtriangledown
d_{h}-\alpha d_{h})_{k}\\\
\quad+(1-\tau_{3}^{-1}\tau_{3}^{\prime})(R_{4}^{h,n,\theta},d_{h})_{k}+\tau_{3}^{-1}\tau_{3}^{\prime}(d_{4},d_{h})_{k}\\}+(TE^{n,\theta},d_{h})\hskip
14.22636pt$ (94)

Here $(\cdot,\cdot)_{k}$ in simple form denotes $(\cdot,\cdot)_{\Omega_{k}}$
Now substituting $\textbf{V}_{h}$ by
$(E^{A,n,\theta}_{u1},E^{A,n,\theta}_{u2},E^{A,n,\theta}_{p},E^{A,n,\theta}_{c})$
in the above equation we will get the second part of $RHS$ as follows

$(\frac{e_{c}^{n+1}-e_{c}^{n}}{dt},E^{A,n,\theta}_{c})+a_{S}(e_{\textbf{u}}^{n,\theta},E^{A,n,\theta}_{\textbf{u}})+a_{T}(e_{c}^{n,\theta},E^{A,n,\theta}_{c}e_{c})-b(E^{A,n,\theta}_{\textbf{u}},e_{p}^{n,\theta})+b(e_{\textbf{u}}^{n,\theta},E^{A,n,\theta}_{p})\\\
=\sum_{k=1}^{n_{el}}\\{\tau_{1}^{\prime}(R_{1}^{h,n,\theta},\mu(c)\Delta
E^{A,n,\theta}_{u1}-\sigma E^{A,n,\theta}_{u1}+\frac{\partial
E^{A,n,\theta}_{p}}{\partial
x})_{k}+\tau_{2}^{\prime}(R_{3}^{h,n,\theta},\bigtriangledown\cdot
E^{A,n,\theta}_{\textbf{u}})_{k}+\\\
\quad\tau_{1}^{\prime}(R_{2}^{h,n,\theta},\mu(c)\Delta
E^{A,n,\theta}_{u2}-\sigma E^{A,n,\theta}_{u2}+\frac{\partial
E^{A,n,\theta}_{p}}{\partial
y})_{k}+(1-\tau_{3}^{-1}\tau_{3}^{\prime})(R_{4}^{h,n,\theta},E^{A,n,\theta}_{c})_{k}+\\\
\quad\tau_{3}^{\prime}(R_{4}^{h,n,\theta}+d_{4},\bigtriangledown\cdot\tilde{\bigtriangledown}E^{A,n,\theta}_{c}+\textbf{u}\cdot\bigtriangledown
E^{A,n,\theta}_{c}-\alpha
E^{A,n,\theta}_{c})_{k}+\tau_{3}^{-1}\tau_{3}^{\prime}(d_{4},E^{A,n,\theta}_{c})_{k}\\}+(TE^{n,\theta},E^{A,n,\theta}_{c})$
(95)

Now we will bound each of the term starting with 4th term of the right hand
side of the above equation.

$\begin{split}\sum_{k=1}^{n_{el}}(1-\tau_{3}^{-1}\tau_{3}^{\prime})(R_{4}^{h,n,\theta},E^{A,n,\theta}_{c})_{k}&\leq\frac{\mid\tau_{3}\mid}{T_{0}-C_{\tau_{3}}}(\sum_{k=1}^{n_{el}}B_{6k})\|R_{4}^{h,n,\theta}\|\end{split}$
(96)

We have obtained this using Cauchy-Schwarz inequality and then imposing bound
on auxiliary error corresponding to $u_{1}$ over each sub-domain $\Omega_{k}$.
Before proceeding further let us look into the form of the column vector d
which has components $d_{1},d_{2},d_{3}$ and $d_{4}$
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}^{h}$
Hence clearly $d_{1}=0$, $d_{2}=0$, $d_{3}=0$ and
$d_{4}=(\sum_{i=1}^{n+1}(\frac{1}{dt}\tau_{3}^{\prime})^{i})R_{4}^{h,n,\theta}$
Now we can bound the last term as follows

$\begin{split}\sum_{k=1}^{n_{el}}\tau_{3}^{-1}\tau_{3}^{\prime}(d_{4},E^{A,n,\theta}_{c})_{k}&=\tau_{3}^{-1}\tau_{3}^{\prime}(\sum_{i=1}^{n+1}(\frac{1}{dt}\tau_{3}^{\prime})^{i})\sum_{k=1}^{n_{el}}(R_{4}^{h,n,\theta},E^{A,n,\theta}_{c})_{k}\\\
&\leq\tau_{3}^{-1}\tau_{3}^{\prime}(\sum_{i=1}^{\infty}(\frac{1}{dt}\tau_{3}^{\prime})^{i})\sum_{k=1}^{n_{el}}(R_{4}^{h,n,\theta},E^{A,n,\theta}_{c})_{k}\\\
&\leq\frac{\mid\tau_{3}\mid}{(T_{0}-C_{\tau_{3}}})(\sum_{k=1}^{n_{el}}B_{6k})\|R_{4}^{h,n,\theta}\|\end{split}$
(97)

Now it will be easy enough to bound the remaining terms of the right hand
side.

$\begin{split}\sum_{k=1}^{n_{el}}\tau_{1}^{\prime}(R_{1}^{h,n,\theta},\mu(c)\Delta
E^{A,n,\theta}_{u1}-\sigma E^{A,n,\theta}_{u1}+\frac{\partial
E^{A,n,\theta}_{p}}{\partial x})_{k}&\\\
\leq\mid\tau_{1}^{\prime}\mid(\sum_{k=1}^{n_{el}}(\mu_{u}B_{2k}+\mid\sigma\mid
B_{1k}+B_{3k}))\|R_{1}^{h,n,\theta}\|\hskip 19.91692pt\\\
\leq\mid\tau_{1}\mid(\sum_{k=1}^{n_{el}}\bar{B}_{1k})\|R_{1}^{h,n,\theta}\|\hskip
113.81102pt\\\
\sum_{k=1}^{n_{el}}\tau_{1}^{\prime}(R_{2}^{h,n,\theta},\mu(c)\Delta
E^{A,n,\theta}_{u2}-\sigma E^{A,n,\theta}_{u2}+\frac{\partial
E^{A,n,\theta}_{p}}{\partial y})_{k}&\\\
\leq\mid\tau_{1}^{\prime}\mid(\sum_{k=1}^{n_{el}}(\mu_{u}B_{2k}^{\prime}+\mid\sigma\mid
B_{1k}^{\prime}+B_{3k}^{\prime}))\|R_{2}^{h,n,\theta}\|\hskip 19.91692pt\\\
\leq\mid\tau_{1}\mid(\sum_{k=1}^{n_{el}}\bar{B}_{2k})\|R_{2}^{h,n,\theta}\|\hskip
5.69054pt(say)\hskip 85.35826pt\\\ \end{split}$ (98)

where $\bar{B}_{1k}$= $(\mu_{u}B_{2k}+\mid\sigma\mid B_{1k}+B_{3k})$ and
$\bar{B}_{2k}$= $(\mu_{u}B_{2k}^{\prime}+\mid\sigma\mid
B_{1k}^{\prime}+B_{3k}^{\prime})$ and

$\sum_{k=1}^{n_{el}}\tau_{3}^{\prime}(R_{4}^{h,n,\theta}+d_{4},\bigtriangledown\cdot\tilde{\bigtriangledown}E^{A,n,\theta}_{c}+\textbf{u}\cdot\bigtriangledown
E^{A,n,\theta}_{c}-\alpha E^{A,n,\theta}_{c})_{k}\\\
=\sum_{k=1}^{n_{el}}\tau_{3}^{\prime}(R_{4}^{h,n,\theta},\bigtriangledown\cdot\tilde{\bigtriangledown}E^{A,n,\theta}_{c}+\textbf{u}\cdot\bigtriangledown
E^{A,n,\theta}_{c}-\alpha
E^{A,n,\theta}_{c})_{k})+\sum_{k=1}^{n_{el}}\tau_{3}^{\prime}(d_{4},\bigtriangledown\cdot\tilde{\bigtriangledown}E^{A,n,\theta}_{c}\\\
+\textbf{u}\cdot\bigtriangledown E^{A,n,\theta}_{c}-\alpha
E^{A,n,\theta}_{c})_{k}\hskip 199.16928pt\\\
\leq(\mid\tau_{3}^{\prime}\mid+\frac{\tau_{3}^{\prime
2}}{dt-\mid\tau_{3}^{\prime}\mid})(\sum_{k=1}^{n_{el}}(D_{1m}B_{7k}+D_{2m}B_{7k}^{\prime}+D_{u1}B_{8k}+D_{u2}B_{8k}^{\prime}+\mid\alpha\mid
B_{6k}))\|R_{4}^{h,n,\theta}\|\hskip 59.75078pt\\\ \leq\frac{\mid\tau_{3}\mid
T_{0}}{T_{0}-C_{\tau_{3}}}(\sum_{k=1}^{n_{el}}\bar{B}_{4k})\|R_{4}^{h,n,\theta}\|\hskip
5.69054pt(say)\hskip 184.9429pt$ (99)

where $\bar{B}_{4k}$=
$(D_{1m}B_{7k}+D_{2m}B_{7k}^{\prime}+D_{u1}B_{8k}+D_{u2}B_{8k}^{\prime}+\mid\alpha\mid
B_{6k})$.
Finally

$\begin{split}\sum_{k=1}^{n_{el}}\tau_{2}^{\prime}(R_{3}^{h,n,\theta},\bigtriangledown\cdot
E^{A,n,\theta}_{\textbf{u}})_{k}&=\sum_{k=1}^{n_{el}}\tau_{2}^{\prime}(\bigtriangledown\cdot
e^{n,\theta}_{\textbf{u}},\bigtriangledown\cdot
E^{A,n,\theta}_{\textbf{u}})_{k}-\sum_{k=1}^{n_{el}}\tau_{2}^{\prime}(\bigtriangledown\cdot
E^{I,n,\theta}_{\textbf{u}},\bigtriangledown\cdot
E^{A,n,\theta}_{\textbf{u}})_{k}\\\ &\leq
C_{\tau_{2}}(1+\epsilon_{2}+\epsilon_{3})\|e^{n,\theta}_{u1}\|_{1}^{2}+C_{\tau_{2}}(1+\frac{1}{\epsilon_{2}}+\frac{1}{\epsilon_{3}})\|e^{n,\theta}_{u2}\|_{1}^{2}\\\
&\quad
C_{\tau_{2}}h^{2}\sum_{i=1}^{2}(\frac{1+\theta}{2}\|u_{i}^{n+1}\|_{2}+\frac{1-\theta}{2}\|u_{i}^{n}\|_{2})\\\
&\leq
C_{\tau_{2}}\\{(1+\epsilon_{2}+\epsilon_{3})\|e^{n,\theta}_{u1}\|_{1}^{2}+(1+\frac{1}{\epsilon_{2}}+\frac{1}{\epsilon_{3}})\|e^{n,\theta}_{u2}\|_{1}^{2}+C_{2}h^{2}\\}\end{split}$
(100)

Applying assumption (iv) on $\|u_{i}^{n}\|_{2}$ for $n=0,1,2,...,N$ we have
the constant $C_{2}$ in the last line of (100). Now this completes finding
bounds for each term in the $RHS$ of (89). Therefore our next work is to
combine all the results into equation (89). 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

$\|e_{c}^{N}\|^{2}+\sum_{n=0}^{N-1}\int_{t^{n}}^{t^{n+1}}\\{(2\mu_{l}-\epsilon_{1}-2C_{\tau_{2}}(1+\epsilon_{2}+\epsilon_{3}))(\|e_{u1}^{n,\theta}\|^{2}+\|\frac{\partial
e_{u1}^{n,\theta}}{\partial x}\|^{2}+\|\frac{\partial
e_{u1}^{n,\theta}}{\partial y}\|^{2})\\\
+(2\mu_{l}-\epsilon_{1}-2C_{\tau_{2}}(1+\frac{1}{\epsilon_{2}}+\frac{1}{\epsilon_{3}}))(\|e_{u2}^{n,\theta}\|^{2}+\|\frac{\partial
e_{u2}^{n,\theta}}{\partial x}\|^{2}+\|\frac{\partial
e_{u2}^{n,\theta}}{\partial y}\|^{2})+\\\
(2\sigma-\epsilon_{1})\|e_{p}^{n,\theta}\|^{2}+(2D_{\alpha}-\epsilon_{1})(\|e_{c}^{n,\theta}\|^{2}+\|\frac{\partial
e_{c}^{n,\theta}}{\partial x}\|^{2}+\|\frac{\partial
e_{c}^{n,\theta}}{\partial y}\|^{2})\\}dt\\\
\leq\frac{C_{2}h^{2}}{\epsilon_{1}}\sum_{n=0}^{N-1}\int_{t^{n}}^{t^{n+1}}(\|R_{1}^{h,n,\theta}\|^{2}+\|R_{2}^{h,n,\theta}\|^{2}+\|R_{3}^{h,n,\theta}\|^{2}+\|R_{4}^{h,n,\theta}\|^{2})dt+\hskip
48.36967pt\\\
+2\mid\tau_{1}\mid\\{(\sum_{k=1}^{n_{el}}\bar{B}_{1k})\sum_{n=0}^{N-1}\int_{t^{n}}^{t^{n+1}}\|R_{1}^{h,n,\theta}\|+(\sum_{k=1}^{n_{el}}\bar{B}_{2k})\sum_{n=0}^{N-1}\int_{t^{n}}^{t^{n+1}}\|R_{2}^{h,n,\theta}\|\\}dt\\\
+\frac{4\mid\tau_{3}\mid}{T_{0}-C_{\tau_{3}}}(\sum_{k=1}^{n_{el}}(B_{6k}+\bar{B}_{4k})\sum_{n=0}^{N-1}\int_{t^{n}}^{t^{n+1}}\|R_{4}^{h,n,\theta}\|dt+\sum_{n=0}^{N-1}\int_{t^{n}}^{t^{n+1}}C_{\tau_{2}}C_{2}h^{2}\hskip
14.22636pt\\\
+(\sum_{k=1}^{n_{el}}B_{6k})\sum_{n=0}^{N-1}\int_{t^{n}}^{t^{n+1}}\|TE^{n,\theta}\|dt\hskip
170.71652pt\\\ \leq\bar{C}(\textbf{R}^{h})(h^{2}+dt^{2})\hskip 213.39566pt$
(101)

Choose the arbitrary parameters 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 backward Euler
time discretisation scheme and its associated property (18) and the fact that
$\tau_{1},\tau_{3}$ are of order $h^{2}$, we have arrived at the above
$aposteriori$ estimate (101), which does not depend upon exact solution. It
shows that the method is second order accurate in space.

##  4 Coupling of the Stokes-Brinkman/Transport equations through interfaces

This section presents a brief study on implementing stabilized $ASGS$ method
on coupled Stokes-Brinkman/Transport model with interface conditions. Here we
have considered the domain $\Omega$ be partitioned into two sub-domains viz.
$\Omega^{S}$ and $\Omega^{B}$ where the fluid flow in $\Omega^{S}$ is governed
by the Stokes equation and in $\Omega^{B}$ the porous media flow obeys the
Brinkman model. Let $\Gamma$ denote the interface and $\partial\Omega^{S}$ and
$\partial\Omega^{B}$ be the boundaries of $\Omega^{S}$ and $\Omega^{B}$
respectively. Now $\Gamma^{l}=\partial\Omega^{l}\setminus\Gamma$ ($l=S,B$).
Let us first mention here the system of fluid flow and mass transport
equations in $\Omega^{S}$: Find $\textbf{u}^{S}$: $\Omega$ $\times$ (0,T)
$\rightarrow R^{2}$ , $p^{S}$: $\Omega\times$ (0,T) $\rightarrow R$ and
$c^{S}$: $\Omega\times$ (0,T) $\rightarrow R$ such that,

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

and the same set of equations in $\Omega^{B}$ is: Find $\textbf{u}^{B}$:
$\Omega$ $\times$ (0,T) $\rightarrow R^{2}$ , $p^{B}$: $\Omega\times$ (0,T)
$\rightarrow R$ and $c^{B}$: $\Omega\times$ (0,T) $\rightarrow R$ such that,

$\begin{split}-\mu^{B}(c)\Delta\textbf{u}^{B}+\bigtriangledown
p^{B}&=\textbf{f}_{1}^{B}\hskip 5.69054ptin\hskip
5.69054pt\Omega^{B}\times(0,T)\\\
\bigtriangledown\cdot\textbf{u}^{B}&=f_{2}^{B}\hskip 5.69054ptin\hskip
5.69054pt\Omega^{B}\times(0,T)\\\ \phi\frac{\partial c^{B}}{\partial
t}-\bigtriangledown\cdot\tilde{\bigtriangledown}c^{B}+\textbf{u}^{B}\cdot\bigtriangledown
c^{B}+\alpha c^{B}&=g^{B}\hskip 5.69054pton\hskip
5.69054pt\Omega^{S}\times(0,T)\\\ \textbf{u}^{B}&=0\hskip 5.69054pton\hskip
5.69054pt\partial\Omega^{B}\times(0,T)\\\
\textbf{u}^{B}&=\textbf{u}_{0}^{B}\hskip 5.69054ptat\hskip 5.69054ptt=0\\\
\tilde{\bigtriangledown}c^{B}\cdot\textbf{n}&=0\hskip 5.69054ptin\hskip
5.69054pt\partial\Omega^{B}\times(0,T)\\\ c^{B}&=c_{0}^{B}\hskip
5.69054ptat\hskip 5.69054ptt=0\\\ \end{split}$ (103)

where $(\textbf{u}^{S},p^{S})$ and $(\textbf{u}^{B},p^{B})$ are the pairs of
Stokes and Brinkman velocities and pressure respectively. As well $c^{S}$ and
$c^{B}$ denote concentration of the solute in $\Omega^{S}$ and $\Omega^{B}$
respectively. $\mu^{S}(c)$ and $\mu^{B}(c)$ are the Stokes dynamic viscosity
and Brinkman effective viscosity respectively Now the interface conditions on
$\Gamma$ are as follows:

$\textbf{u}^{S}\cdot\textbf{n}^{S}+\textbf{u}^{B}\cdot\textbf{n}^{B}=0$ (104)
$-\mu^{S}(c)\partial_{n}\textbf{u}^{S}\cdot\textbf{n}^{S}+p^{S}=-\mu^{B}(c)\partial_{n}\textbf{u}^{B}\cdot\textbf{n}^{B}+p^{B}$
(105)
$\mu^{S}(c)\partial_{n}\textbf{u}^{S}\cdot\textbf{t}+\frac{\alpha}{\sqrt{\sigma}}\textbf{u}^{S}\cdot\textbf{t}=0$
(106) $c^{S}=c^{B}$ (107)
$\tilde{\bigtriangledown}c^{S}\cdot\textbf{n}^{S}+\tilde{\bigtriangledown}c^{B}\cdot\textbf{n}^{B}=0$
(108)

where $\textbf{n}^{S}$ and $\textbf{n}^{B}$ are the outward normals to
$\Omega^{S}$ and $\Omega^{B}$ respectively. It is quite obvious to observe
that $\textbf{n}^{S}=-\textbf{n}^{B}$. Here $(104),(107),(108)$ represent
continuity conditions of normal velocities and concentration. Where as $(105)$
enforces continuity of normal stresses and $(106)$ is the Beavers-Joseph-
Saffman condition [7].
Let us mention the corresponding spaces for both the sub-problems here: Let
$V^{l}=H_{0}^{1}(\Omega^{l})$ and $Q^{l}=L^{2}(\Omega^{l})$ (for $l=S,B$).
$V_{h}^{l}$ and $Q_{h}^{l}$ be the corresponding finite dimensional subspaces
of $V_{h}^{l}$ and $Q_{h}^{l}$ respectively for $l=S,B$. Following the
previous steps the stabilized formulation for Stokes sub-problem is to find
$\textbf{U}_{h}^{S}=(\textbf{u}_{h}^{S},p_{h}^{S},c_{h}^{S})$: J $\rightarrow$
$V_{h}^{S}\times V_{h}^{S}\times Q_{h}^{S}\times V_{h}^{S}$ such that
$\forall\textbf{V}_{h}^{S}=(\textbf{v}_{h}^{S},q_{h}^{S},d_{h}^{S})$ $\in$
$V_{h}^{S}\times V_{h}^{S}\times Q_{h}^{S}\times V_{h}^{S}$

$(M\partial_{t}\textbf{U}^{S}_{h},\textbf{V}^{S}_{h})+B^{S}_{ASGS}(\textbf{U}^{S}_{h},\textbf{V}^{S}_{h})=L^{S}_{ASGS}(\textbf{V}^{S}_{h})+\int_{\Gamma}(\mu^{S}(c)\partial_{n}\textbf{u}^{S}-p^{S}\textbf{n}^{S})\cdot\textbf{v}_{h}^{S}d\Gamma$
(109)

and for Brinkman sub-problem the stabilized formulation is to find
$\textbf{U}_{h}^{B}=(\textbf{u}_{h}^{B},p_{h}^{B},c_{h}^{B})$: J $\rightarrow$
$V_{h}^{B}\times V_{h}^{B}\times Q_{h}^{B}\times V_{h}^{B}$ such that
$\forall\textbf{V}_{h}^{B}=(\textbf{v}_{h}^{B},q_{h}^{B},d_{h}^{B})$ $\in$
$V_{h}^{B}\times V_{h}^{B}\times Q_{h}^{B}\times V_{h}^{B}$

$(M\partial_{t}\textbf{U}^{B}_{h},\textbf{V}^{B}_{h})+B^{B}_{ASGS}(\textbf{U}^{B}_{h},\textbf{V}^{B}_{h})=L^{B}_{ASGS}(\textbf{V}^{B}_{h})+\int_{\Gamma}(\mu^{B}(c)\partial_{n}\textbf{u}^{B}-p^{B}\textbf{n}^{B})\cdot\textbf{v}_{h}^{B}d\Gamma$
(110)

Now applying the interface conditions it is easy to conclude that the
stabilized $ASGS$ formulation for coupled Stokes-Brinkman/Transport model is
to find $\textbf{U}_{h}^{S}=(\textbf{u}_{h}^{S},p_{h}^{S},c_{h}^{S})$: J
$\rightarrow$ $V_{h}^{S}\times V_{h}^{S}\times Q_{h}^{S}\times V_{h}^{S}$ and
$\textbf{U}_{h}^{B}=(\textbf{u}_{h}^{B},p_{h}^{B},c_{h}^{B})$: J $\rightarrow$
$V_{h}^{B}\times V_{h}^{B}\times Q_{h}^{B}\times V_{h}^{B}$ such that
$\forall\textbf{V}_{h}^{S}=(\textbf{v}_{h}^{S},q_{h}^{S},d_{h}^{S})$ $\in$
$V_{h}^{S}\times V_{h}^{S}\times Q_{h}^{S}\times V_{h}^{S}$ and
$\forall\textbf{V}_{h}^{B}=(\textbf{v}_{h}^{B},q_{h}^{B},d_{h}^{B})$ $\in$
$V_{h}^{B}\times V_{h}^{B}\times Q_{h}^{B}\times V_{h}^{B}$

$(M\partial_{t}\textbf{U}^{S}_{h},\textbf{V}^{S}_{h})+(M\partial_{t}\textbf{U}^{B}_{h},\textbf{V}^{B}_{h})+B^{S}_{ASGS}(\textbf{U}^{S}_{h},\textbf{V}^{S}_{h})+B^{B}_{ASGS}(\textbf{U}^{B}_{h},\textbf{V}^{B}_{h})\\\
+\frac{\alpha}{\sqrt{\sigma}}(\textbf{u}_{h}^{S}\cdot\textbf{t},\textbf{v}_{h}^{S}\cdot\textbf{t})=L^{S}_{ASGS}(\textbf{V}^{S}_{h})+L^{B}_{ASGS}(\textbf{V}^{B}_{h})$
(111)

######  Remark 4.

These terms in $(111)$ do not much differ from the general form of stabilized
$ASGS$ formulation in $(8)$. Only one term consisting of the tangential
component of Stokes velocity is extra in the coupled stabilized formulation
$(111)$. Hence $apriori$ and $aposteriori$ error estimate results will be
almost similar; only coefficients $C^{\prime\prime}$ for $apriori$ and
$\bar{C}(\textbf{R}^{h})$ for $aposteriori$ may sightly differ though the
convergence rates for both that error estimations will be the same.

##  5 Numerical Experiment

In this section we present a comparative study between standard Galerkin
method and stabilized algebraic subgrid scale(ASGS) method as well as we have
verified the convergence rate established theoretically under stabilized
method in the previous sections. We have considered three different models to
work with: Coupled Stokes/Transport Model, Coupled Brinkman/Transport Model
and Coupled Stokes-Brinkman with interface/Transport Model.
For simplicity we have considered bounded square domain $\Omega$= (0,1)
$\times$ (0,1). We have taken continuous piecewise linear finite element(P1)
space into account for approximating velocity, pressure and concentration too.
The expression of concentration dependent viscosity is taken from [20], which
establishes that viscosity of a solvent depends upon concentration of the
solute of a electrolyte solution. The proposed expression for viscosity is
$\mu(c)=0.954e^{27.93\times 0.028c}$.
Let us mention here the exact solutions for all three cases as follows:
$\textbf{u}=(tsin^{2}(\pi x)sin(\pi y)cos(\pi y),-tsin(\pi x)cos(\pi
x)sin^{2}(\pi y))$,
$p=tsin(2\pi x)sin(2\pi y)$ and $c=txy(x-1)(y-1)$
Now in the following we mention the general expressions of the coefficients
involved in the equations:
The reaction coefficients $\alpha=0.01$
The diffusion coefficients: $D_{1}=t^{2}(sin(\pi x))^{4}(sin(2\pi y))^{2}$,
$D_{2}$=$t^{2}(sin(2\pi x))^{2}$ $(sin(\pi y))^{4}$
The stabilization parameters: $\tau_{1}=(4\frac{\mu_{l}}{h^{2}}+\sigma)^{-1}$,
$\tau_{2}=4\mu_{l}h$ and
$\tau_{3}=19(\frac{9}{4h^{2}}+\frac{3}{2h}+\alpha)^{-1}$ , where
$\mu_{l}=0.954e^{27.93\times 0.028\times 0.0625}$
1.Coupled Stokes/ Transport Model In this case $\sigma=0$ and porosity
$\phi=1$. Hence the expression for stabilization parameters are changed
accordingly.
Table 1 and table 2 present the error in $V$ norm (which is standard norm on
the space V, introduced in the section 2.1 and is defined in the initial part
of section 3.3) and order of convergence under Galerkin method and ASGS method
respectively for this case. These tables are clearly showing that both the
methods perform equally well.
2.Coupled Brinkman/ Transport Model In this case $\sigma$ is non-zero and in
particular we consider $\sigma=1$. The value of porosity is taken to be 2. In
Brinkman flow problem we will deal with effective viscosity $\mu_{B}$.
According to [21] the effective viscosity $\mu_{B}$ and viscosity $\mu$ is
related through $\sigma^{2}=\frac{\mu_{B}}{\mu}$. Hence both are considered
same since $\sigma=1$. Here the stabilization parameters take the general
form.
Table 3 and table 4 present the error in $V$ norm and order of convergence
under Galerkin method and ASGS method respectively. The tables represent that
both methods perform equally well.
3.Coupled Stokes-Brinkman/ Transport Model Here $(\textbf{u},p,c)$ take values
$(\textbf{u}^{S},p^{S},c^{S})$ in $\Omega^{S}$ and
$(\textbf{u}^{B},p^{B},c^{B})$ in $\Omega^{B}$. As mentioned in the previous
case $\mu^{S}(c)=\mu^{B}(c)$ for $\sigma=1$ on $\Omega^{B}$. The stabilization
parameters on $\Omega^{S}$: $\tau_{1}^{S}=\frac{h^{2}}{4\mu_{l}}$,
$\tau_{2}^{S}=4\mu_{l}h$ and
$\tau_{3}^{S}=19(\frac{9}{4h^{2}}+\frac{3}{2h}+\alpha)^{-1}$ and on
$\Omega^{B}$: $\tau_{1}^{B}=(4\frac{\mu_{l}}{h^{2}}+\sigma)^{-1}$,
$\tau_{2}^{B}=4\mu_{l}h$ and
$\tau_{3}^{B}=19(\frac{9}{4h^{2}}+\frac{3}{2h}+\alpha)^{-1}$ ,for given
$\mu_{l}=0.954e^{27.93\times 0.028\times 0.0625}$. The value of porosity
$\phi$ is 2 on $\Omega^{B}$.
Table 5 and table 6 show the error in $V$ norm and order of convergence under
Galerkin method and $ASGS$ method respectively for coupled Stokes-Brinkman/
Transport model. These tables represent that stabilized $ASGS$ method performs
well, whereas the convergence rate under the Galerkin method oscillates.

Grid | Error in $H^{1}$ norm | Order of convergence
---|---|---
10$\times$ 10 | 0.950341 |
20 $\times$ 20 | 0.27489 | 1.78959
40 $\times$ 40 | 0.0635241 | 2.11348
80 $\times$ 80 | 0.0190773 | 1.73544
160 $\times$ 160 | 0.0053094 | 1.84524
Table 1: Error and Order of convergence obtained in $V$ norm under Galerkin method for Stokes/Transport Model Grid | Error in $H^{1}$ norm | Order of convergence
---|---|---
10 $\times$ 10 | 0.200567 |
20 $\times$ 20 | 0.0661861 | 1.59948
40 $\times$ 40 | 0.0162986 | 2.02178
80 $\times$ 80 | 0.00434506 | 1.9073
160 $\times$ 160 | 0.00113881 | 1.93185
Table 2: Error and Order of convergence obtained in $V$ norm under ASGS method for Stokes/Transport Model Grid | Error in $H^{1}$ norm | Order of convergence
---|---|---
10 $\times$ 10 | 0.953771 |
20 $\times$ 20 | 0.275123 | 1.79357
40 $\times$ 40 | 0.0635331 | 2.1145
80 $\times$ 80 | 0.0190837 | 1.73516
160 $\times$ 160 | 0.00531838 | 1.84328
Table 3: Error and Order of convergence obtained in $V$ norm under Galerkin method for Brinkman/Transport Model Grid | Error in $H^{1}$ norm | Order of convergence
---|---|---
10 $\times$ 10 | 0.201123 |
20 $\times$ 20 | 0.0661698 | 1.60383
40 $\times$ 40 | 0.0162947 | 2.02177
80 $\times$ 80 | 0.00435187 | 1.90469
160 $\times$ 160 | 0.0011484 | 1.92201
Table 4: Error and Order of convergence obtained in $V$ norm under ASGS method for Brinkman/Transport Model  Grid | Error in $H^{1}$ norm | Order of convergence
---|---|---
10 $\times$ 10 | 0.000603299 |
20 $\times$ 20 | 0.00057062 | 0.080342
40 $\times$ 40 | 0.000139967 | 2.02744
80 $\times$ 80 | 3.04023$e^{-5}$ | 2.20284
160 $\times$ 160 | 2.06788$e^{-5}$ | 0.556032
Table 5: Error and Order of convergence obtained in $V$ norm under Galerkin method for coupled Stokes-Brinkman/Transport Model Grid | Error in $H^{1}$ norm | Order of convergence
---|---|---
10 $\times$ 10 | 0.000479849 |
20 $\times$ 20 | 0.00015364 | 1.64303
40 $\times$ 40 | 4.25356$e^{-5}$ | 1.85281
80 $\times$ 80 | 1.1035$e^{-5}$ | 1.94664
160 $\times$ 160 | 2.80343$e^{-6}$ | 1.97685
Table 6: Error and Order of convergence obtained in $V$ norm under ASGS method
for coupled Stokes-Brinkman/Transport Model

Figure 1: Mesh for 40 $\times$ 40 grid points and Horizontal velocity plot
respectively for both Coupled Stokes/Transport Model and Coupled
Brinkman/Transport Model

Figure 2: Concentration plot for both Coupled Stokes/Transport Model and
Coupled Brinkman/Transport Model

Figure 3: Concentration plot for both Coupled Stokes-Brinkman/Transport Model
with interface

Figure 4: Velocity plot and Velocity concentration plot respectively for both
Coupled Stokes/Transport Model and Coupled Brinkman/Transport Model

Figure 5: Velocity plot and Velocity concentration plot respectively for
Coupled Stokes-Brinkman/Transport Model

Figure 6: Comparison of solution plots at an arbitrary point (0.5,0.5)
(upper); Error plot in $H1$-norm with respect to degrees of freedom (middle);
Order of convergence plot (lower) for Coupled Stokes/Transport Model

Figure 7: Comparison of solution plots at an arbitrary point (0.5,0.5)
(upper); Error plot in $H1$-norm with respect to degrees of freedom (middle);
Order of convergence plot (lower) for Coupled Brinkman/Transport Model

Figure 8: Comparison of solution plots at an arbitrary point (0.5,0.5)
(upper); Error plot in $H1$-norm with respect to degrees of freedom (middle);
Order of convergence plot (lower) for Coupled Stokes-Brinkman/Transport Model
with interface conditions

Figure 9: Exact solution for 40 $\times$ 40 grid points

Figure 10: ASGS solution for 40 $\times$ 40 grid points

######  Remark 5.

The tables are showing that error under $ASGS$ method at each mesh size is
turned out to be lesser than that of Galerkin method and for all three cases
the order of convergence under $ASGS$ method is 2, which justifies
theoretically established result.

######  Remark 6.

It is clear from figure 6,7, and 8 that $ASGS$ solution comparatively more
fast converges to exact solution whereas the Galerkin solution converges
slowly in first two cases and diverges in third case representing coupled
Stokes-Brikman/Transport model.

######  Remark 7.

The error plots in $H1$-norm and order of convergence plot under Galerkin and
$ASGS$ methods in figure 6,7,8 establish the more efficiency of the stabilized
method in compared to Galerkin method. It shows that error under $ASGS$ method
is much lesser than that of Galerkin method at the same mesh size and both are
decreasing for finer mesh.

######  Remark 8.

Figure 9 and 10 present three dimensional view of exact solution and
approximated solution derived under stabilized method. It is easily seen that
the approximated solution is very much alike to the exact one.

##  6 Conclusion

The paper presents $ASGS$ stabilized finite element analysis of two different
aspects of Stokes-Brinkman fluid flow model strongly coupled with unsteady
$VADR$ transport equation; one is unified way of considering the model and
another is coupling system of equations through interface conditions. Whereas
this paper in one hand elaborately derive both apriori and aposteriori error
estimates, on other hand it highlights the way to prove existence and
uniqueness of the solution of variational formulation. 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 apriori error estimation turns out
to be $O(h+h^{2}+dt^{2})$ and for aposteriori it is $O(h^{2}+dt^{2})$ for
backward Euler time discretization method. In numerical experiment section
three cases viz. coupled Stokes/Transport model, coupled Brinkman/Transport
model, coupled Stokes-Brinkman/Transport model with interface conditions, have
been considered to cover all the different aspects of the model and in all of
the three cases stabilized $ASGS$ method presents better performance with
respect to standard Galerkin method.
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] D. Vassilev, I. Yotov, Coupling Stokes-Darcy flow with transport, SIAM J. Sci. Comput. 3661-3684(2009).
  * [2] 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).
  * [3] A. Cesmelioglu, P. Chidyagwai, B. Rivi$\grave{e}$re, Continuous and discontinuous finite element methods for coupled surface-subsurface flow and transport problems.
  * [4] G. Houzeaux, J. Principe, A variational subgrid scale model for transient incompressible flows, International Journal of Computational Fluid Dynamics. 135-152(2008).
  * [5] 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).
  * [6] 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).
  * [7] S. Badia, R. Codina, Unified stabilized finite element formulations for the Stokes and the Darcy problems, SIAM J. NUMER. ANAL, Vol. 47, No. 3, pp. 1971–2000(2009).
  * [8] J. A. Evans, T.J.R. Hughes, Isogeometric Divergence-conforming B-splines for the Darcy-Stokes-Brinkman Equations, ICES REPORT 12-03 (2012).
  * [9] A. Masud, A stabilized mixed finite element method for Darcy–Stokes flow, Int. J. Numer. Meth. Fluids ; 54:665–681(2007).
  * [10] A. Hannukainen, M. Juntunen, R. Stenberg, Computations with finite element methods for the Brinkman problem, Comput Geosci 15:155–166(2011).
  * [11] M. Juntunen, R. Stenberg, Analysis of finite element methods for the Brinkman problem, Calcolo 47: 129–147(2010).
  * [12] X. Xie, J. Xu, G. Xue, Uniformly-stable finite element methods for Darcy-Stokes-Brinkman models, Journal of Computational Mathematics, Vol.26, No.3, 437-455(2008).
  * [13] T. Lévy, Loi de Darcy ou loi de Brinkman? C. R. Acad. Sci. Paris Sér. II Méc. Phys. Chim. Sci. Univers Sci. Terre 292(12), 871–874, Erratum (17):1239 (1981)
  * [14] T. Arbogast, H.L. Lehr, Homogenization of a darcy-stokes system modeling vuggy porous media. Comput. Geosci. 10(3), 291–302 (2006)
  * [15] Y. Cao, M. Gunzburger,F. Hua, X. Wang, Coupled Stokes-Darcy Model With Beavers-Joseph Interface Boundary Condition, Commun. Math. Sci, Vol. 8, No. 1, pp. 1–25(2010).
  * [16] M. Mu, J. Xu, A Two-Grid Method of a Mixed Stokes–Darcy Model for Coupling Fluid Flow with porous media flow, SIAM J. NUMER. ANAL, Vol. 45, No. 5, pp. 1801–1813(2007).
  * [17] M. Cai, M. Mu, A multilevel decoupled method for a mixed Stokes/Darcy model, Journal of Computational and Applied Mathematics, 236, 2452-2465(2012).
  * [18] B. Riviere, Analysis of a Discontinuous Finite Element Method for the Coupled Stokes and Darcy Problems, Journal of Scientific Computing, Volumes 22 and 23,( 2005)
  * [19] Y. Amanbek, M.F. Wheeler, A priori error analysis for transient problems using Enhanced Velocity approach in the discrete-time setting, Journal of Computational and Applied Mathematics 361, 459-471(2019).
  * [20] 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).
  * [21] G. Neale, W. Nader, Practical significance of Brinkman’s extension of Darcy’s law: coupled parallel flows within a channel and a bounding porous, Canad. J. Chem. Eng. 52 475–478(1974).
  * [22] I. Harari, T.J.R. Hughes, Stabilized finite element methods for steady advection-diffusion with production, Comput. Methods Appl. Mech. Engrg. 115 165-191(1994).
  * [23] R. Codina, On stabilized finite element methods for linear systems of convection-diffusion-reaction equations, Comput. Methods Appl. Mech. Engrg. 188 61-82(2000).
  * [24] 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).
  * [25] L. P. Franca, G. Hauke, A. Masud, Revisiting stabilized finite element methods for the advective–diffusive equation, Comput. Methods Appl. Mech. Engrg. 195 1560–1572(2006).
  * [26] G. Hauke, A simple subgrid scale stabilized method for the advection–diffusion-reaction equation, Comput. Methods Appl. Mech. Engrg. 191 2925–2947(2002).
  * [27] A. Sendur, A Comparative Study on Stabilized Finite Element Methods for the Convection-Diffusion-Reaction Problems, Hindawi Journal of Applied Mathematics, Article ID 4259634, 16 pages(2018).
  * [28] T.J.R. Hughes,G. N. Wells, Conservation properties for the Galerkin and stabilised forms of the advection–diffusion and incompressible Navier–Stokes equations, Comput. Methods Appl. Mech. Engrg. 194 1141–1159(2005).
  * [29] 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, Comput. Methods Appl. Mech. Engrg., Vol.32 Issues 1-3 (1982).
  * [30] R. Codina, O. Soto, Approximation of the incompressible Navier–Stokes equations using orthogonal subscale stabilization and pressure segregation on anisotropic finite element meshes, Comput. Methods Appl. Mech. Engrg., Vol 193 Issues 15-16 (2004).
  * [31] J. Blasco, R. Codina, Space and time error estimates for a first order pressure stabilized finite element method for the incompressible NavierStokes equations (2001).
  * [32] L.P.Franca, S.L.Frey, Stabilized finite element methods: II. The incompressible Navier-Stokes equations, Comput. Methods Appl. Mech. Engrg. Vol. 99, Issues 2-3(1992).
  * [33] A. Masud, T.J.R. Hughes, A space-time Galerkin/least-squares finite element formulation of the Navier-Stokes equations for moving domain problems , Comput. Methods Appl. Mech. Engrg. 146 , 91-126(1997).
  * [34] T.J.R. Hughes, Multiscale phenomena: Green’s functions, the Dirichlet-toNeumann formulation, subgrid scale models, bubbles and the origins of stabilized methods, Comput. Methods Appl. Mech. Engrg. 127 (1995) 387-401
  * [35] R. Codina, Comparison of some finite element methods for solving the diffusion-convection-reaction equation, Comput. Methods Appl. Mech. Engrg. 156 185-210 (1998).
  * [36] W.R. Hwang, S. G. Advani, Numerical simulations of Stokes–Brinkman equations for permeability prediction of dual scale fibrous porous media, Phys. Fluids 22, 113101 (2010)
  * [37] P. Benner, S. Dolgov, A. Onwunta, M. Stoll, Low-rank solvers for unsteady Stokes–Brinkman optimal control problem with random data, Comput. Methods Appl. Mech. Engrg. 304 26-54 (2016)