* Leroy et al. (2021) Leroy, A. K., Schinnerer, E., Hughes, A., et al. 2021, The Astrophysical Journal Supplement Series, 257, 43, doi: 10.3847/1538-4365/ac17f3
  * Levy et al. (2021) Levy, R. C., Bolatto, A. D., Leroy, A. K., et al. 2021, The Astrophysical Journal, 912, 4, doi: 10.3847/1538-4357/abec84
  * Levy et al. (2022) —. 2022, The Astrophysical Journal, 935, 19, doi: 10.3847/1538-4357/ac7b7a
  * Linden et al. (2021) Linden, S. T., Evans, A. S., Larson, K., et al. 2021, The Astrophysical Journal, 923, 278, doi: 10.3847/1538-4357/ac2892
  * Linden et al. (2022) Linden, S. T., Perez, G., Calzetti, D., et al. 2022, The Astrophysical Journal, 935, 166, doi: 10.3847/1538-4357/ac7c07
  * Linden et al. (2023) Linden, S. T., Evans, A. S., Armus, L., et al. 2023, The Astrophysical Journal, 944, L55, doi: 10.3847/2041-8213/acb335
  * Maraston (1998) Maraston, C. 1998, Monthly Notices of the Royal Astronomical Society, 300, 872, doi: 10.1046/j.1365-8711.1998.01947.x
  * Mayya et al. (2008) Mayya, Y. D., Romano, R., Rodríguez-Merino, L. H., et al. 2008, ApJ, 679, 404, doi: 10.1086/587541
  * Mayya et al. (2008) Mayya, Y. D., Romano, R., Rodríguez-Merino, L. H., et al. 2008, The Astrophysical Journal, 679, 404, doi: 10.1086/587541
  * Messa et al. (2018) Messa, M., Adamo, A., Östlin, G., et al. 2018, Monthly Notices of the Royal Astronomical Society, 473, 996, doi: 10.1093/mnras/stx2403
  * Moustakas et al. (2010) Moustakas, J., Kennicutt, Robert C., J., Tremonti, C. A., et al. 2010, ApJS, 190, 233, doi: 10.1088/0067-0049/190/2/233
  * Nguyen Luong et al. (2011) Nguyen Luong, Q., Motte, F., Schuller, F., et al. 2011, Astronomy and Astrophysics, 529, A41, doi: 10.1051/0004-6361/201016271
  * Noeske et al. (2007) Noeske, K. G., Weiner, B. J., Faber, S. M., et al. 2007, The Astrophysical Journal, 660, L43, doi: 10.1086/517926
  * Noll et al. (2009) Noll, S., Burgarella, D., Giovannoli, E., et al. 2009, Astronomy and Astrophysics, 507, 1793, doi: 10.1051/0004-6361/200912497
  * Orozco-Duarte et al. (2022) Orozco-Duarte, R., Wofford, A., Vidal-García, A., et al. 2022, MNRAS, 509, 522, doi: 10.1093/mnras/stab2988
  * Peng et al. (2010) Peng, Y.-j., Lilly, S. J., Kovač, K., et al. 2010, The Astrophysical Journal, 721, 193, doi: 10.1088/0004-637X/721/1/193
  * Popescu et al. (2012) Popescu, B., Hanson, M. M., & Elmegreen, B. G. 2012, The Astrophysical Journal, 751, 122, doi: 10.1088/0004-637X/751/2/122
  * Portegies Zwart et al. (2010) Portegies Zwart, S. F., McMillan, S. L. W., & Gieles, M. 2010, Annual Review of Astronomy and Astrophysics, 48, 431, doi: 10.1146/annurev-astro-081309-130834
  * Pérez et al. (2021) Pérez, G., Messa, M., Calzetti, D., et al. 2021, The Astrophysical Journal, 907, 100, doi: 10.3847/1538-4357/abceba
  * Querejeta et al. (2021) Querejeta, M., Schinnerer, E., Meidt, S., et al. 2021, Astronomy & Astrophysics, 656, A133, doi: 10.1051/0004-6361/202140695
  * Renaud (2018) Renaud, F. 2018, New Astronomy Reviews, 81, 1, doi: 10.1016/j.newar.2018.03.001
  * Rico-Villas et al. (2020) Rico-Villas, F., Martín-Pintado, J., González-Alfonso, E., Martín, S., & Rivilla, V. M. 2020, Monthly Notices of the Royal Astronomical Society, 491, 4573, doi: 10.1093/mnras/stz3347
  * Rodríguez et al. (2023) Rodríguez, M. J., Lee, J. C., Whitmore, B. C., et al. 2023, The Astrophysical Journal, 944, L26, doi: 10.3847/2041-8213/aca653
  * Ryon et al. (2017) Ryon, J. E., Gallagher, J. S., Smith, L. J., et al. 2017, The Astrophysical Journal, 841, 92, doi: 10.3847/1538-4357/aa719e
  * Salim et al. (2007) Salim, S., Rich, R. M., Charlot, S., et al. 2007, The Astrophysical Journal Supplement Series, 173, 267, doi: 10.1086/519218
  * Salim et al. (2016) Salim, S., Lee, J. C., Janowiecki, S., et al. 2016, The Astrophysical Journal Supplement Series, 227, 2, doi: 10.3847/0067-0049/227/1/2
  * Scheuermann et al. (2023) Scheuermann, F., Kreckel, K., Barnes, A. T., et al. 2023, Monthly Notices of the Royal Astronomical Society, 522, 2369, doi: 10.1093/mnras/stad878
  * Schinnerer et al. (2019) Schinnerer, E., Leroy, A., Blanc, G., et al. 2019, The Messenger, 177, 36, doi: 10.18727/0722-6691/5151
  * Schinnerer et al. (2023) Schinnerer, E., Emsellem, E., Henshaw, J. D., et al. 2023, The Astrophysical Journal, 944, L15, doi: 10.3847/2041-8213/acac9e
  * Searle et al. (1980) Searle, L., Wilkinson, A., & Bagnuolo, W. G. 1980, The Astrophysical Journal, 239, 803, doi: 10.1086/158165
  * Sellwood & Masters (2022) Sellwood, J. A., & Masters, K. L. 2022, ARA&A, 60, doi: 10.1146/annurev-astro-052920-104505
  * Sellwood & Wilkinson (1993) Sellwood, J. A., & Wilkinson, A. 1993, Reports on Progress in Physics, 56, 173, doi: 10.1088/0034-4885/56/2/001
  * Silva-Villa & Larsen (2011) Silva-Villa, E., & Larsen, S. S. 2011, Astronomy and Astrophysics, 529, A25, doi: 10.1051/0004-6361/201016206
  * Simonyan & Zisserman (2015) Simonyan, K., & Zisserman, A. 2015, Very Deep Convolutional Networks for Large-Scale Image Recognition, arXiv.  http://arxiv.org/abs/1409.1556
  * Skillman et al. (1996) Skillman, E. D., Kennicutt, Robert C., J., Shields, G. A., & Zaritsky, D. 1996, ApJ, 462, 147, doi: 10.1086/177138
  * Smith et al. (2007) Smith, L. J., Bastian, N., Konstantopoulos, I. S., et al. 2007, The Astrophysical Journal, 667, L145, doi: 10.1086/522120
  * Sormani et al. (2020) Sormani, M. C., Tress, R. G., Glover, S. C. O., et al. 2020, Monthly Notices of the Royal Astronomical Society, 497, 5024, doi: 10.1093/mnras/staa1999
  * Sormani et al. (2023) Sormani, M. C., Barnes, A. T., Sun, J., et al. 2023, Monthly Notices of the Royal Astronomical Society, 523, 2918, doi: 10.1093/mnras/stad1554
  * Stetson (1987) Stetson, P. B. 1987, Publications of the Astronomical Society of the Pacific, 99, 191, doi: 10.1086/131977
  * Sun et al. (2024) Sun, J., He, H., Batschkun, K., et al. 2024, Hidden Gems on a Ring: Infant Massive Clusters and Their Formation Timeline Unveiled by ALMA, HST, and JWST in NGC 3351, doi: 10.48550/arXiv.2401.14453
  * Tang et al. (2023) Tang, J., Grasha, K., & Krumholz, M. R. 2023, Cluster Population Demographics in NGC 628 Derived from Stochastic Population Synthesis Models, doi: 10.48550/arXiv.2301.05912
  * Thilker et al. (2024) Thilker, D., Lee, J., Whitmore, B., Maschmann, D., & Deger, S. 2024
  * Thilker et al. (2022) Thilker, D. A., Whitmore, B. C., Lee, J. C., et al. 2022, Monthly Notices of the Royal Astronomical Society, 509, 4094, doi: 10.1093/mnras/stab3183
  * Thilker et al. (2023) Thilker, D. A., Lee, J. C., Deger, S., et al. 2023, The Astrophysical Journal, 944, L13, doi: 10.3847/2041-8213/acaeac
  * Tress et al. (2020) Tress, R. G., Sormani, M. C., Glover, S. C. O., et al. 2020, Monthly Notices of the Royal Astronomical Society, 499, 4455, doi: 10.1093/mnras/staa3120
  * Turner et al. (2021) Turner, J. A., Dale, D. A., Lee, J. C., et al. 2021, Monthly Notices of the Royal Astronomical Society, 502, 1366, doi: 10.1093/mnras/stab055
  * Turner et al. (2022) Turner, J. A., Dale, D. A., Lilly, J., et al. 2022, MNRAS, 516, 4612, doi: 10.1093/mnras/stac2559
  * van den Bergh & Hagen (1968) van den Bergh, S., & Hagen, G. L. 1968, The Astronomical Journal, 73, 569, doi: 10.1086/110663
  * Vollmer et al. (2004) Vollmer, B., Balkowski, C., Cayatte, V., van Driel, W., & Huchtmeier, W. 2004, Astronomy and Astrophysics, 419, 35, doi: 10.1051/0004-6361:20034552
  * Vázquez & Leitherer (2005) Vázquez, G. A., & Leitherer, C. 2005, The Astrophysical Journal, 621, 695, doi: 10.1086/427866
  * Ward & Kruijssen (2018) Ward, J. L., & Kruijssen, J. M. D. 2018, MNRAS, 475, 5659, doi: 10.1093/mnras/sty117
  * Ward et al. (2020) Ward, J. L., Kruijssen, J. M. D., & Rix, H.-W. 2020, MNRAS, 495, 663, doi: 10.1093/mnras/staa1056
  * Watkins et al. (2023a) Watkins, E. J., Kreckel, K., Groves, B., et al. 2023a, Astronomy and Astrophysics, 676, A67, doi: 10.1051/0004-6361/202346075
  * Watkins et al. (2023b) Watkins, E. J., Barnes, A. T., Henny, K., et al. 2023b, The Astrophysical Journal, 944, L24, doi: 10.3847/2041-8213/aca6e4
  * Wei et al. (2020) Wei, W., Huerta, E. A., Whitmore, B. C., et al. 2020, Monthly Notices of the Royal Astronomical Society, 493, 3178, doi: 10.1093/mnras/staa325
  * Whitmore et al. (2010) Whitmore, B. C., Chandar, R., Schweizer, F., et al. 2010, The Astronomical Journal, 140, 75, doi: 10.1088/0004-6256/140/1/75
  * Whitmore et al. (2021) Whitmore, B. C., Lee, J. C., Chandar, R., et al. 2021, Monthly Notices of the Royal Astronomical Society, 506, 5294, doi: 10.1093/mnras/stab2087
  * Whitmore et al. (2023a) Whitmore, B. C., Chandar, R., Lee, J. C., et al. 2023a, Monthly Notices of the Royal Astronomical Society, 520, 63, doi: 10.1093/mnras/stad098
  * Whitmore et al. (2023b) Whitmore, B. C., Chandar, R., Rodríguez, M. J., et al. 2023b, The Astrophysical Journal, 944, L14, doi: 10.3847/2041-8213/acae94
  * Williams et al. (2024) Williams, T. G., Lee, J. C., Larson, K. L., et al. 2024, PHANGS-JWST: Data Processing Pipeline and First Full Public Data Release, doi: 10.48550/arXiv.2401.15142
  * Wofford et al. (2016) Wofford, A., Charlot, S., Bruzual, G., et al. 2016, Monthly Notices of the Royal Astronomical Society, 457, 4296, doi: 10.1093/mnras/stw150
  * Wright (2020) Wright, N. J. 2020, New A Rev., 90, 101549, doi: 10.1016/j.newar.2020.101549
  * Zaritsky et al. (1994) Zaritsky, D., Kennicutt, Robert C., J., & Huchra, J. P. 1994, ApJ, 420, 87, doi: 10.1086/173544

##  Appendix A Additional figures

In this section we show additional figures. Figure 16 presents the correlation
plot between host galaxy distances and relative fraction of characteristic
regions in color-color diagrams as discussed in Section 5. Figures 17 to
26present the spatial distribution plots discussed in Section 6.

Figure 16: Number fraction of C1 and C2 clusters of each galaxy associated
with the main characteristic regions in color-color diagrams found in Section
4.4 as a function of galaxy distance. We show the YCL, the MAP and the OGC in
blue, green and red, respectively. In gray, we show clusters outside the main
regions. We distinguish distance measurements which are estimated from stellar
markers such as Tip of the Red giant Branch (TRGB) or from Cepheid variable
stars are marked with full circles whereas other distant measurements which
are less precise are marked by empty circles. A complete discussion on each
individual distance measurement is provided in Anand et al. (2020) and Anand
et al. (2021). For each panel we show the Pearson correlation coefficient in
the top right. Figure 17: Spatial distributions for ML-classified star
clusters of class 1 and 2 as categorized (Sect. 4.4) into three groups: OGC
(red, top panels), MAP (green, upper middle panels) and YCL (blue, lower
middle panels); plus color-composite images created from the HST U-B-V bands
(bottom panels). The cluster spatial distribution maps are produced by binning
the cluster positions onto a pixel grid which is subsequently convolved with a
Gaussian and normalized to unity. In order to highlight the relation between
young clusters and molecular gas, with magenta lines, we overlay the ALMA
CO(2-1) intensity contours of the 95 percentile on the YCL distribution maps.
In this figure and Figures 18–26, we show the spatial cluster distribution for
all PHANGS–HST galaxies sorted by decreasing $\Delta$MS values (see Figure
13). Figure 18: Continuation of Figure 17 Figure 19: Continuation of Figure 17
Figure 20: Continuation of Figure 17 Figure 21: Continuation of Figure 17
Figure 22: Continuation of Figure 17 Figure 23: Continuation of Figure 17. We
note that for IC 5332 there is no ALMA CO(2-1) detection. Figure 24:
Continuation of Figure 17 Figure 25: Continuation of Figure 17 Figure 26:
Continuation of Figure 17


# Stabilized subgrid multiscale finite element formulation for advection-
diffusion-reaction equation with variable coefficients coupled with Stokes-
Darcy equation

Manisha Chowdhury , B.V. Rathish Kumar  Email addresses:
<EMAIL_ADDRESS>(M. Chowdhury) and<EMAIL_ADDRESS>(B.V.R.
Kumar)

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

###### Abstract

In this paper subgrid multiscale stabilized finite element method for
Advection-Diffusion-Reaction (ADR) equation coupled with Stokes-Darcy flow
problem has been studied. Here the advection velocity involved in ADR equation
obeys Stokes-Darcy flow equation. In this study the approach of algebraic
approximation of stabilization parameter has been considered. Further apriori
error estimation has been elaborately carried out.

Keywords Stokes-Darcy-Brinkman equation $\cdot$ Advection-diffusion-reaction
equation $\cdot$ Subgrid scale method $\cdot$ A priori error estimation

##  1 Introduction

Transport phenomena, mathematically expressed in terms of Advection-Diffusion-
Reaction(ADR) equation has always been an active area of research in the
fields of Biomedical Engineering, Environmental Sciences, Chemical Engineering
etc. Our previous works [4],[6] have focused on transport equation with
spatially variable diffusion and advection coefficients. In this paper we have
considered only diffusion coefficients as variable and the advection velocity
comes from Stokes-Darcy flow problem. There are few studies [1], [2] available
in literature on coupled Stokes-Darcy transport equation, but none of them has
studied stabilized subgrid multiscale finite element method for both the
equations. Being most general stabilization method now-a-days subgrid method
has been considered more suitable finite element method to be dealt with. It
involves split of the unknown true solution is chosen to be the standard
Galerkin finite element solution and obtaining the unresolvable solution in
terms of the known resolvable one, we will finally arrive at the subgrid
formulation. [4] has derived an expression of the stabilization parameter for
subgrid formulation of ADR equation with variable coefficients and another
study [5] has found out the same for Stokes-Darcy equation. Both these studies
have followed the approach of algebraic approximation of the parameter.
Combining those results we will have the required stabilization parameter for
this coupled problem. Here further we have derived apriori error estimation
the coupled problem.
The paper is organised as: section 2 presents model problem along with weak
formulation and subgrid formulation for the weakly coupled problem. In the
next section we have carried out apriori error estimation for this stabilized
method with the help of interpolation estimates for introducing projection
operators.

##  2 Statement of the problem

Let $\Omega$ be an open bounded domain in $R^{d}$, d=2,3. Here we have worked
with two dimensional model for the sake of simplicity in further derivation,
but it can be extended to three dimensional model straightforward.
For an incompressible solvent fluid the Stokes-Darcy (or Brinkman) equation
representing its flow is to find $u:\Omega\rightarrow R^{2}$ and
$p:\Omega\rightarrow R$ such that

$\begin{split}-\mu\Delta\textbf{u}+\sigma\textbf{u}+\bigtriangledown
p&=\textbf{f}\hskip 5.69054ptin\hskip 5.69054pt\Omega\\\
\bigtriangledown\cdot\textbf{u}&=0\hskip 5.69054ptin\hskip 5.69054pt\Omega\\\
\textbf{u}&=\textbf{0}\hskip 5.69054pton\hskip 5.69054pt\partial\Omega\\\
\end{split}$ (1)

where u= ($u_{1},u_{2}$),p and $\mu$ are velocity, pressure and viscosity of
the fluid, $\sigma$ is the inverse of permeability and f is the body force.
And the ADR equation with spatially variable coefficients representing the
concentration c of transporting solute in the same domain $\Omega$ along with
homogeneous Dirichlet boundary condition is to find c: $\Omega\rightarrow R$
such that,

$\begin{split}-\bigtriangledown\cdot\tilde{\bigtriangledown}c+\textbf{u}\cdot\bigtriangledown
c+\alpha c&=g\hskip 5.69054ptin\hskip 5.69054pt\Omega\\\ c&=0\hskip
5.69054pton\hskip 5.69054pt\Omega\\\ \end{split}$ (2)

where the notation, $\tilde{\bigtriangledown}:=(D_{1}\frac{\partial}{\partial
x},D_{2}\frac{\partial}{\partial y})$
$D_{1},D_{2}$ are spatially variable diffusion coefficients along x-axis and
y-axis respectively, $\alpha$ is the reaction coefficient and g denotes the
source of solute mass. Here we make some assumptions on the coefficients as
follows:
(i) $\mu$, $\sigma$ and $\alpha$ are positive constants.
(ii) The diffusion coefficients $D_{1},D_{2}$ are continuous functions on
bounded domain $\Omega$and hence bounded. Let $D_{1l},D_{2l}$ be upper bounds
of them respectively.
(iii) The body force $\textbf{f}\in(L^{2}(\Omega))^{2}$ and source $g\in
L^{2}(\Omega)$

###  2.1 Weak formulation

Let us consider the standard space $V=H^{1}_{0}(\Omega)$ as the admissible
space for both velocity fields and concentration and the space
$Q=L^{2}(\Omega)$ for pressure.
The weak formulation is to find $(\textbf{u},p)\in V\times V\times Q$ and
$c\in V$ such that

$\begin{split}a_{S}(\textbf{u},\textbf{v})-b(\textbf{v},p)+b(\textbf{u},q)&=l_{S}(\textbf{v})\hskip
5.69054pt\forall\textbf{v}\in V\times V,q\in Q\\\ a_{T}(c,d)&=l_{T}(d)\hskip
5.69054pt\forall d\in V\end{split}$ (3)

where
$a_{S}(\textbf{u},\textbf{v})=\int_{\Omega}\mu\bigtriangledown\textbf{u}:\bigtriangledown\textbf{v}+\sigma\int_{\Omega}\textbf{u}\cdot\textbf{v}$
and $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}(\textbf{v})=\int_{\Omega}\textbf{f}\cdot\textbf{v}$ and
$l_{T}(d)=\int_{\Omega}gd$
The stability of the continuous problem (3) has been shown in [9] and [6].

###  2.2 Subgrid formulation

Let $\Omega$ be discretized into total $n_{el}$ number of sub-domains
$\Omega_{k}$ for k=1,2,…,$n_{el}$ and $h_{k}$ be the diameter of each sub-
domain respectively. Let h= $\underset{k=1,2,...,n_{el}}{max}h_{k}$ and
$\tilde{\Omega}=\bigcup_{k=1}^{n_{e}l}\Omega_{k}$ be the union of interior
elements. Let $V_{h}$ and $Q_{h}$ be the suitable finite dimensional subspaces
of $V$ and $Q$ respectively where
$V_{h}=\\{v\in V:v(\Omega_{k})=\mathcal{P}^{2}(\Omega_{k})\\}$ and
$Q_{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}$.
Now the standard Galerkin finite element formulation of (3) is to find
$(\textbf{u}_{h},p_{h})\in V_{h}\times V_{h}\times Q_{h}$ and $c\in V_{h}$
such that

$\begin{split}a_{S}(\textbf{u}_{h},\textbf{v}_{h})-b(\textbf{v}_{h},p_{h})+b(\textbf{u}_{h},q_{h})&=l_{S}(\textbf{v}_{h})\hskip
5.69054pt\forall\textbf{v}_{h}\in V_{h}\times V_{h},q_{h}\in Q_{h}\\\
a_{T}(c_{h},d_{h})&=l_{T}(d_{h})\hskip 5.69054pt\forall d_{h}\in
V_{h}\end{split}$ (4)

As we earlier mention that for subgrid formulation the finite element solution
will be taken as resolvable scale and we will obtain the subgrid formulation
by expressing the unresolvable scale in terms of known solution. Hence the
subgrid formulation for (3) is to find $(\textbf{u}_{h},p_{h})\in V_{h}\times
V_{h}\times Q_{h}$ and $c\in V_{h}$ such that

$\begin{split}a_{S}(\textbf{u}_{h},\textbf{v}_{h})-b(\textbf{v}_{h},p_{h})+b(\textbf{u}_{h},q_{h})+\int_{\Omega^{\prime}}(-\mathcal{L}_{1}^{*}\textbf{V}_{h})[\tau]\mathcal{L}_{1}\textbf{U}_{h}&=l_{S}(\textbf{v}_{h})+\int_{\Omega^{\prime}}(-\mathcal{L}_{1}^{*}\textbf{V}_{h})[\tau]\textbf{F}\\\
\forall\hskip 2.84526pt\textbf{V}_{h}&=(\textbf{v}_{h},q_{h})\in V_{h}\times
V_{h}\times Q_{h}\\\
a_{T}(c_{h},d_{h})+\int_{\Omega^{\prime}}(-\mathcal{L}_{2}^{*}d_{h})\tau_{3}\mathcal{L}_{2}c_{h}&=l_{T}(d_{h})+\int_{\Omega^{\prime}}(-\mathcal{L}_{2}^{*}d_{h})\tau_{3}g\\\
&\forall\hskip 2.84526ptd_{h}\in V_{h}\end{split}$ (5)

where

$\mathcal{L}_{1}\textbf{U}=\begin{bmatrix}-\mu(c)\Delta\textbf{u}+\sigma\textbf{u}+\bigtriangledown
p\\\ \bigtriangledown\cdot\textbf{u}\\\ \end{bmatrix},its\hskip
2.84526ptadjoint\mathcal{L}_{1}^{*}\textbf{U}=\begin{bmatrix}-\mu(c)\Delta\textbf{u}+\sigma\textbf{u}-\bigtriangledown
p\\\ -\bigtriangledown\cdot\textbf{u}\\\
\end{bmatrix},\textbf{F}=\begin{bmatrix}\textbf{f}\\\ 0\end{bmatrix}$

$\mathcal{L}_{2}c=-\bigtriangledown\cdot\widetilde{\bigtriangledown}c+\overline{u}\cdot\bigtriangledown
c+\alpha c$ and its $adjoint$
$\mathcal{L}_{2}^{*}c=-\bigtriangledown\cdot\widetilde{\bigtriangledown}c-\overline{u}\cdot\bigtriangledown
c+\alpha c$
the stabilization parameters[],[] are obtained as

$[\tau]=diag(\tau_{1},\tau_{1},\tau_{2})=\begin{bmatrix}(c_{1}\frac{\mu}{h^{2}}+\sigma)^{-1}&0&0\\\
0&(c_{1}\frac{\mu}{h^{2}}+\sigma)^{-1}&0\\\ 0&0&c_{2}\mu\end{bmatrix}and\hskip
2.84526pt\tau_{3}=(\frac{9D}{4h^{2}}+\frac{3U}{2h}+\alpha)^{-1}$

##  3 Error estimation

In this section we are going to derive apriori estimate in $V$ norm that is in
standard $H^{1}$ norm. Before that we introduce projection operators and carry
out error splitting as follows:
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})$. Let us introduce the projection operator corresponding
to each component as the following,
(i)For any $\textbf{u}\in V\times V$ let there exist an interpolation
$P^{h}_{\textbf{u}}:V\times V\longrightarrow V_{h}\times V_{h}$ satisfying
$b(\textbf{u}-P^{h}_{\textbf{u}}\textbf{u},q_{h})=0$ $\forall q_{h}\in Q_{h}$
and component wise satisfy $L^{2}$-orthogonality condition i.e. for any
$u_{i}\in V$ $(u_{i}-P^{h}_{u_{i}}u_{i},v_{h})=0$ $\forall v_{h}\in V_{s}^{h}$
for i=1,2.
(ii) Let $P^{h}_{p}:Q\longrightarrow Q_{h}$ be the $L^{2}$ orthogonal
projection given by
$\int_{\Omega}(p-P^{h}_{p}p)q_{h}=0$ $\forall p\in Q,\hskip 2.84526pt\forall
q_{h}\in Q_{h}$
(iii) Similarly let $P^{h}_{c}:V\longrightarrow V_{h}$ satisfy
$\int_{\Omega}(c-P^{h}_{c}c)v_{h}=0$ $\forall c\in V,\hskip 2.84526pt\forall
v_{h}\in V_{h}$
Now each components 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}-P^{h}_{u1}u_{1})+(P^{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}$
Interpolation estimates [7, 8] : for any true solution $u$ with regularity
upto (m+1)

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

where l 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. We will use $\|\cdot\|$ instead of $\|\cdot\|_{0}$ to
denote $L^{2}(\Omega)$ norm.
Now we carry out apriori error estimate in two parts. In first part we will
bound auxiliary error and in the later part using the result obtained in first
part we will find final form of apriori estimate.

######  Theorem 1.

Auxiliary error estimate: For velocity $\textbf{u}_{h}=(u_{1h},u_{2h})$,
pressure $p_{h}$ and concentration $c_{h}$ belonging to $V_{h}\times
V_{h},Q_{h},V_{h}$ respectively satisfying (5), assume sufficient regularity
of exact solution in equations (1)-(2). Then there exists constants $C_{1}$
and $C_{2}$, depending upon u,p and c respectively, such that

$\begin{split}\|E^{A}_{u1}\|^{2}_{V}+\|E^{A}_{u2}\|^{2}_{V}+\|E^{A}_{p}\|_{Q}^{2}&\leq
C_{1}(\textbf{u},p)h^{2}\\\ \|E^{A}_{c}\|_{V}^{2}&\leq
C_{2}(c)h^{2}\end{split}$ (7)

###### Proof.

We divide the proof into two parts: first part contains estimation of
auxiliary error for velocity and pressure components separately and in the
later part we will find bound for auxiliary error corresponding to
concentration. Now for true solution (5) becomes

$\begin{split}a_{S}(\textbf{u},\textbf{v}_{h})-b(\textbf{v}_{h},p)+b(\textbf{u},q_{h})&=l_{S}(\textbf{v}_{h})\hskip
5.69054pt\forall\textbf{v}_{h}\in V_{h}\times V_{h},q_{h}\in Q_{h}\\\
a_{T}(c,d_{h})&=l_{T}(d_{h})\hskip 5.69054pt\forall d_{h}\in V_{h}\end{split}$
(8)

First part: Subtracting first equation of (8) from that of (5) we will have,
$\forall\textbf{v}_{h}\in V_{h}\times V_{h},q_{h}\in Q_{h}$

$a_{S}(\textbf{u}-\textbf{u}_{h},\textbf{v}_{h})-b(\textbf{v}_{h},p-p_{h})+b(\textbf{u}-\textbf{u}_{h},q_{h})=\int_{\Omega^{\prime}}(\mathcal{L}_{1}^{*}\textbf{V}_{h})[\tau](\textbf{F}-\mathcal{L}_{1}\textbf{U}_{h})$
(9)

After explicitly writing the bilinear forms and using the error splitting we
have the above equation as follows

$\mu\int_{\Omega}\bigtriangledown
E^{I}_{\textbf{u}}:\bigtriangledown\textbf{v}_{h}+\mu\int_{\Omega}\bigtriangledown
E^{A}_{\textbf{u}}:\bigtriangledown\textbf{v}_{h}+\sigma\int_{\Omega}E^{I}_{\textbf{u}}\cdot\textbf{v}_{h}+\sigma\int_{\Omega}E^{A}_{\textbf{u}}\cdot\textbf{v}_{h}-\int_{\Omega}(\bigtriangledown\cdot\textbf{v}_{h})E^{I}_{p}-\\\
\int_{\Omega}(\bigtriangledown\cdot\textbf{v}_{h})E^{A}_{p}+\int_{\Omega}(\bigtriangledown\cdot
E^{I}_{\textbf{u}})q_{h}+\int_{\Omega}(\bigtriangledown\cdot
E^{A}_{\textbf{u}})q_{h}\\\ =\int_{\Omega^{\prime}}\tau_{1}(\mu\Delta
v_{1h}-\sigma v_{1h}+\frac{\partial q_{h}}{\partial
x})(-\mu\Delta(E^{I}_{u1}+E^{A}_{u1})+\sigma(E^{I}_{u1}+E^{A}_{u1})+\frac{\partial
E^{I}_{p}}{\partial x}+\frac{\partial E^{A}_{p}}{\partial x})+\\\
\quad\int_{\Omega^{\prime}}\tau_{1}(\mu\Delta v_{2h}-\sigma
v_{2h}+\frac{\partial q_{h}}{\partial
y})(-\mu\Delta(E^{I}_{u2}+E^{A}_{u2})+\sigma(E^{I}_{u2}+E^{A}_{u2})+\frac{\partial
E^{I}_{p}}{\partial y}+\frac{\partial E^{A}_{p}}{\partial y})+\\\
\quad\int_{\Omega}\tau_{2}(\bigtriangledown\cdot\textbf{v}_{h})(\bigtriangledown\cdot
E^{I}_{\textbf{u}})+\int_{\Omega}\tau_{2}(\bigtriangledown\cdot\textbf{v}_{h})(\bigtriangledown\cdot
E^{A}_{\textbf{u}})\hskip 113.81102pt$ (10)

Since this holds for all $\textbf{v}_{h}\in V_{h}\times V_{h}$ and $q_{h}\in
Q_{h}$, replacing them by $E^{A}_{\textbf{u}}$ and $E^{A}_{p}$ respectively
and using properties of projection operators, we finally have

$\mu\int_{\Omega}\bigtriangledown E^{A}_{\textbf{u}}:\bigtriangledown
E^{A}_{\textbf{u}}+\sigma\int_{\Omega}E^{A}_{\textbf{u}}\cdot
E^{A}_{\textbf{u}}=\int_{\Omega^{\prime}}\tau_{1}(\mu\Delta E^{A}_{u1}-\sigma
E^{A}_{u1}+\frac{\partial E^{A}_{p}}{\partial x})(-\mu\Delta E^{I}_{u1}+\sigma
E^{I}_{u1}+\\\ \quad\frac{\partial E^{I}_{p}}{\partial
x})+\int_{\Omega^{\prime}}\tau_{1}(\mu\Delta E^{A}_{u1}-\sigma
E^{A}_{u1}+\frac{\partial E^{A}_{p}}{\partial x})(-\mu\Delta E^{A}_{u1}+\sigma
E^{A}_{u1}+\frac{\partial E^{A}_{p}}{\partial
x})+\int_{\Omega^{\prime}}\tau_{1}(\mu\Delta E^{A}_{u2}-\sigma E^{A}_{u2}\\\
\quad+\frac{\partial E^{A}_{p}}{\partial y})(-\mu\Delta E^{I}_{u2}+\sigma
E^{I}_{u2}+\frac{\partial E^{I}_{p}}{\partial
y})+\int_{\Omega^{\prime}}\tau_{1}(\mu\Delta E^{A}_{u2}-\sigma
E^{A}_{u2}+\frac{\partial E^{A}_{p}}{\partial y})(-\mu\Delta E^{A}_{u2}+\sigma
E^{A}_{u2}+\frac{\partial E^{A}_{p}}{\partial y})\\\
\quad+\int_{\Omega^{\prime}}\tau_{2}(\bigtriangledown\cdot
E^{A}_{\textbf{u}})(\bigtriangledown\cdot
E^{I}_{\textbf{u}})+\int_{\Omega^{\prime}}\tau_{2}(\bigtriangledown\cdot
E^{A}_{\textbf{u}})(\bigtriangledown\cdot
E^{A}_{\textbf{u}})-\mu\int_{\Omega}\bigtriangledown
E^{I}_{\textbf{u}}:\bigtriangledown
E^{A}_{\textbf{u}}-\int_{\Omega}(\bigtriangledown\cdot
E^{I}_{\textbf{u}})E^{A}_{p}\\\
=S_{1}+S_{2}+S_{3}+S_{4}+S_{5}+S_{6}+S_{7}+S_{8}\hskip 5.69054pt(say)\hskip
99.58464pt$ (11)

Now we estimate each of the above terms separately. Before that here we like
to make an important observation: by the virtue of chosen finite element
spaces it can be clearly said that every element belonging to $V_{h}$ and
$Q_{h}$ and their derivatives all are bounded functions over each sub-domain
$\Omega_{k}$. Therefore let us consider the positive numbers
$M_{1k},M_{2k},M_{3k},M_{4k},M_{5k},M_{1k}^{\prime},M_{2k}^{\prime},M_{3k}^{\prime},M_{4k}^{\prime},M_{5k}^{\prime}$
as bounds for $E^{A}_{u1},\Delta E^{A}_{u1},\frac{\partial E^{A}_{p}}{\partial
x},\frac{\partial E^{A}_{u1}}{\partial x},\frac{\partial E^{A}_{u1}}{\partial
y},E^{A}_{u2},\Delta E^{A}_{u2},\frac{\partial E^{A}_{p}}{\partial
y},\frac{\partial E^{A}_{u2}}{\partial x},\frac{\partial E^{A}_{u2}}{\partial
y}$ respectively.
To bound $T_{1}$ simply multiplying the terms, then applying Cauchy-Schwarz
inequality over each of them and finally using bounds for $E^{A}s$ we will
have

$\begin{split}S_{1}&\leq\mid\tau_{1}\mid\\{(\sum_{k=1}^{n_{el}}(\mu^{2}M_{2k}+\sigma\mid\mu\mid
M_{1k}+\mid\mu\mid
M_{3k}))\|E^{I}_{u1}\|_{2}+(\sum_{k=1}^{n_{el}}(\sigma\mid\mu\mid M_{2k}+\\\
&\quad\sigma^{2}M_{1k}+\mid\sigma\mid
M_{3k}))\|E^{I}_{u1}\|+(\sum_{k=1}^{n_{el}}(\mid\mu\mid M_{2k}+\sigma
M_{1k}+M_{3k}))\|E^{I}_{p}\|_{1}\\}\\\ &\leq
C\mid\tau_{1}\mid\\{(\sum_{k=1}^{n_{el}}(\mu^{2}M_{2k}+\sigma\mid\mu\mid
M_{1k}+\mid\mu\mid
M_{3k}))\|u_{1}\|_{2}+h^{2}(\sum_{k=1}^{n_{el}}(\sigma\mid\mu\mid M_{2k}\\\
&\quad+\sigma^{2}M_{1k}+\mid\sigma\mid
M_{3k}))\|u_{1}\|_{2}+(\sum_{k=1}^{n_{el}}(\mid\mu\mid M_{2k}+\sigma
M_{1k}+M_{3k}))\|p\|_{1}\\}\end{split}$ (12)

Applying similar arguments we have the following bounds

$\begin{split}S_{2}&\leq\mid\tau_{1}\mid\sum_{k=1}^{n_{el}}(\mu^{2}M_{2k}^{2}+2\sigma\mid\mu\mid
M_{1k}M_{2k}+\sigma^{2}M_{1k}^{2}+M_{3k}^{2})\\\ S_{3}&\leq
C\mid\tau_{1}\mid\\{(\sum_{k=1}^{n_{el}}(\mu^{2}M_{2k}^{\prime}+\sigma\mid\mu\mid
M_{1k}^{\prime}+\mid\mu\mid
M_{3k}^{\prime}))\|u_{2}\|_{2}+h^{2}(\sum_{k=1}^{n_{el}}(\sigma\mid\mu\mid
M_{2k}^{\prime}\\\ &\quad+\sigma^{2}M_{1k}^{\prime}+\mid\sigma\mid
M_{3k}^{\prime}))\|u_{2}\|_{2}+(\sum_{k=1}^{n_{el}}(\mid\mu\mid
M_{2k}^{\prime}+\sigma M_{1k}^{\prime}+M_{3k}^{\prime}))\|p\|_{1}\\}\\\
S_{4}&\leq\mid\tau_{1}\mid\sum_{k=1}^{n_{el}}(\mu^{2}M_{2k}^{\prime
2}+2\sigma\mid\mu\mid M_{1k}^{\prime}M_{2k}^{\prime}+\sigma^{2}M_{1k}^{\prime
2}+M_{3k}^{\prime 2})\\\ S_{5}&\leq\mid\tau_{2}\mid
Ch\sum_{k=1}^{n_{el}}(M_{4k}\|u_{1}\|_{2}+M_{5k}^{\prime}\|u_{2}\|_{2})\\\
S_{6}&\leq\mid\tau_{2}\mid\sum_{k=1}^{n_{el}}(M_{4k}^{2}+M_{5k}^{\prime
2})\end{split}$ (13)

Applying Cauchy-Schwarz and Young’s inequality we get

$\begin{split}S_{7}&=-\mu\int_{\Omega}(\frac{\partial E^{I}_{u1}}{\partial
x}\frac{\partial E^{A}_{u1}}{\partial x}+\frac{\partial E^{I}_{u1}}{\partial
y}\frac{\partial E^{A}_{u1}}{\partial y}+\frac{\partial E^{I}_{u2}}{\partial
x}\frac{\partial E^{A}_{u2}}{\partial x}+\frac{\partial E^{I}_{u2}}{\partial
y}\frac{\partial E^{A}_{u2}}{\partial y})\\\ \end{split}$ (14)
$\begin{split}&\leq\mid\mu\mid\\{\frac{1}{2\epsilon_{1}}(\|\frac{\partial
E^{I}_{u1}}{\partial x}\|^{2}+\|\frac{\partial E^{I}_{u1}}{\partial
y}\|^{2}+\|\frac{\partial E^{I}_{u2}}{\partial x}\|^{2}+\|\frac{\partial
E^{I}_{u2}}{\partial y}\|^{2})+\frac{\epsilon_{1}}{2}(\|\frac{\partial
E^{A}_{u1}}{\partial x}\|^{2}+\\\ &\quad\|\frac{\partial E^{A}_{u1}}{\partial
y}\|^{2}+\|\frac{\partial E^{A}_{u2}}{\partial x}\|^{2}+\|\frac{\partial
E^{A}_{u2}}{\partial y}\|^{2})\\}\\\
&\leq\mid\mu\mid\\{C\frac{h^{2}}{2\epsilon_{1}}(\|u_{1}\|_{2}^{2}+\|u_{2}\|_{2}^{2})+\frac{\epsilon_{1}}{2}(\|\frac{\partial
E^{A}_{u1}}{\partial x}\|^{2}+\|\frac{\partial E^{A}_{u1}}{\partial
y}\|^{2}+\|\frac{\partial E^{A}_{u2}}{\partial x}\|^{2}+\|\frac{\partial
E^{A}_{u2}}{\partial y}\|^{2})\\}\end{split}$ (15)

Applying similar argument

$\begin{split}S_{8}=\int_{\Omega}(\frac{\partial E^{I}_{u1}}{\partial
x}+\frac{\partial E^{I}_{u2}}{\partial y})E^{A}_{p}&\leq(\|\frac{\partial
E^{I}_{u1}}{\partial x}\|+\|\frac{\partial E^{I}_{u2}}{\partial
y}\|)\|E^{A}_{p}\|\\\ &\leq
C\frac{h^{2}}{\epsilon_{2}}(\|u_{1}\|_{2}^{2}+\|u_{2}\|_{2}^{2})+\frac{\epsilon_{2}}{2}\|E^{A}_{p}\|^{2}\end{split}$
(16)

This completes finding bounds for each of the terms in the right hand side of
(11). Now we put all the results obtained above into the equation (11) and
take out those terms, which are common with the terms of $LHS$, in the $LHS$
from the $RHS$ of the equation. Consequently in the $RHS$ we remain with the
terms multiplied by $h^{2}$ as the stabilization parameters are of order
$h^{2}$. Now we consider $\epsilon_{1},\epsilon_{2}$ in such a manner that all
the coefficients in $LHS$ can be made positive. After that taking minimum of
all those coefficients and dividing both the sides by that we finally arrive
at the completion of the proof as follows:

$\|E^{A}_{u1}\|_{V}^{2}+\|E^{A}_{u2}\|_{V}^{2}+\|E^{A}_{p}\|_{Q}^{2}\leq
C(\textbf{u},p)h^{2}$ (17)

Second part: Subtracting second equation of (8) from that of (5) we will have,
$\forall d_{h}\in V_{h}$

$a_{T}(c-c_{h},d_{h})=\int_{\Omega^{\prime}}(\mathcal{L}_{2}^{*}d_{h})\tau_{3}(g-\mathcal{L}_{2}c_{h})=\int_{\Omega^{\prime}}(\mathcal{L}_{2}^{*}d_{h})\tau_{3}\mathcal{L}_{2}(c-c_{h})$
(18)

After explicitly writing the bilinear forms and using the error splitting we
have the above equation as follows

$\int_{\Omega}\tilde{\bigtriangledown}E^{I}_{c}\cdot\bigtriangledown
d_{h}+\int_{\Omega}\tilde{\bigtriangledown}E^{A}_{c}\cdot\bigtriangledown
d_{h}+\int_{\Omega}d_{h}\textbf{u}\cdot\bigtriangledown
E^{I}_{c}+\int_{\Omega}d_{h}\textbf{u}\cdot\bigtriangledown
E^{A}_{c}+\alpha\int_{\Omega}E^{I}_{c}d_{h}+\\\
\quad\alpha\int_{\Omega}E^{A}_{c}d_{h}=\int_{\Omega^{\prime}}(-\bigtriangledown\cdot\widetilde{\bigtriangledown}d_{h}-\overline{u}\cdot\bigtriangledown
d_{h}+\alpha
d_{h})\tau_{3}(-\bigtriangledown\cdot\widetilde{\bigtriangledown}E^{I}_{c}+\overline{u}\cdot\bigtriangledown
E^{I}_{c}+\alpha E^{I}_{c})\\\
\int_{\Omega^{\prime}}(-\bigtriangledown\cdot\widetilde{\bigtriangledown}d_{h}-\overline{u}\cdot\bigtriangledown
d_{h}+\alpha
d_{h})\tau_{3}(-\bigtriangledown\cdot\widetilde{\bigtriangledown}E^{A}_{c}+\overline{u}\cdot\bigtriangledown
E^{A}_{c}+\alpha E^{A}_{c})$ (19)

Since this holds for all $d_{h}\in V_{h}$, replacing it by $E^{A}_{c}$ in the
above equation and using properties of projection operators, we finally have

$\int_{\Omega}\tilde{\bigtriangledown}E^{A}_{c}\cdot\bigtriangledown
E^{A}_{c}+\alpha\int_{\Omega}E^{A}_{c}E^{A}_{c}=\int_{\Omega^{\prime}}(-\bigtriangledown\cdot\widetilde{\bigtriangledown}E^{A}_{c}-\overline{u}\cdot\bigtriangledown
E^{A}_{c}+\alpha
E^{A}_{c})\tau_{3}(-\bigtriangledown\cdot\widetilde{\bigtriangledown}E^{I}_{c}+\\\
\quad\overline{u}\cdot\bigtriangledown E^{I}_{c}+\alpha
E^{I}_{c})+\int_{\Omega^{\prime}}(-\bigtriangledown\cdot\widetilde{\bigtriangledown}E^{A}_{c}-\overline{u}\cdot\bigtriangledown
E^{A}_{c}+\alpha
E^{A}_{c})\tau_{3}(-\bigtriangledown\cdot\widetilde{\bigtriangledown}E^{A}_{c}+\overline{u}\cdot\bigtriangledown
E^{A}_{c}+\alpha E^{A}_{c})\\\
\quad-\int_{\Omega}\tilde{\bigtriangledown}E^{I}_{c}\cdot\bigtriangledown
E^{A}_{c}-\int_{\Omega}E^{A}_{c}\textbf{u}\cdot\bigtriangledown
E^{I}_{c}-\int_{\Omega}E^{A}_{c}\textbf{u}\cdot\bigtriangledown E^{A}_{c}\\\
=T_{1}+T_{2}+T_{3}+T_{4}+T_{5}\hskip 2.84526pt(say)\hskip 170.71652pt$ (20)

We will follow the same procedure as the previous part. Under the same
observation let us consider the positive numbers
$M_{6k},M_{7k},M_{8k},M_{7k}^{\prime},M_{8k}^{\prime}$ as bounds for
$E^{A}_{c},\frac{\partial E^{A}_{c}}{\partial
x},\frac{\partial^{2}E^{A}_{c}}{\partial x^{2}},\frac{\partial
E^{A}_{c}}{\partial x},\frac{\partial E^{A}_{c}}{\partial y}$ respectively
over each sub-domain.
First simply multiplying the terms of $T_{1}$, then applying Cauchy-Schwarz
inequality over each of them and finally using bounds for $E^{A}s$ we can find
estimate of $T_{1}$ as follows:

$\begin{split}T_{1}&\leq\mid\tau_{3}\mid(\sum_{k=1}^{n_{el}}(D_{1l}M_{8k}+D_{2l}M_{8k}^{\prime}+D_{u1}M_{7k}+D_{u2}M_{7k}^{\prime}+\alpha
M_{6k}))\\{(D_{1l}+D_{2l})\\\
&\quad\|E^{I}_{c}\|_{2}+(\bar{D}_{u1}+\bar{D}_{u2})\|E^{I}_{c}\|_{1}+\alpha\|E^{I}_{c}\|\\}\\\
&\leq\mid\tau_{3}\mid
C(\sum_{k=1}^{n_{el}}(D_{1l}M_{8k}+D_{2l}M_{8k}^{\prime}+D_{u1}M_{7k}+D_{u2}M_{7k}^{\prime}+\alpha
M_{6k}))\\{(D_{1l}+D_{2l})\\\
&\quad+h(\bar{D}_{u1}+\bar{D}_{u2})+h^{2}\alpha\\}\|c\|_{2}\end{split}$ (21)

Applying similar arguments we have the following bounds

$\begin{split}T_{2}&\leq\mid\tau_{3}\mid\sum_{k=1}^{n_{el}}\\{(D_{1l}M_{8k}+D_{2l}M_{8k}^{\prime}+D_{u1}M_{7k}+D_{u2}M_{7k}^{\prime}+\alpha
M_{6k})(D_{1l}M_{8k}+\\\ &\quad
D_{2l}M_{8k}^{\prime}+\bar{D}_{u1}M_{7k}+\bar{D}_{u2}M_{7k}^{\prime}+\alpha
M_{6k})\\}\end{split}$ (22)

For estimating next terms we are going to use Cauchy-Schwarz and Young’s
inequalities as follows:

$\begin{split}T_{3}&=-\int_{\Omega}(D_{1}\frac{\partial E^{I}_{c}}{\partial
x}\frac{\partial E^{A}_{c}}{\partial x}+D_{2}\frac{\partial
E^{I}_{c}}{\partial y}\frac{\partial E^{A}_{c}}{\partial y})\\\
&\leq\frac{D_{1l}}{2\epsilon_{3}}\|E^{I}_{c}\|_{1}^{2}+\frac{D_{2l}\epsilon_{3}}{2}(\|\frac{\partial
E^{I}_{c}}{\partial x}\|^{2}+\|\frac{\partial E^{I}_{c}}{\partial y}\|^{2})\\\
&\leq\frac{D_{1l}}{2\epsilon_{3}}h^{2}\|c\|_{2}^{2}+\frac{D_{2l}\epsilon_{3}}{2}(\|\frac{\partial
E^{I}_{c}}{\partial x}\|^{2}+\|\frac{\partial E^{I}_{c}}{\partial y}\|^{2})\\\
\end{split}$ (23) $\begin{split}T_{4}&=-\int_{\Omega}(u_{1}\frac{\partial
E^{I}_{c}}{\partial x}+u_{2}\frac{\partial E^{I}_{c}}{\partial y})E^{A}_{c}\\\
&\leq\frac{(u_{1l}+u_{2l})}{2\epsilon_{4}}\|E^{I}_{c}\|_{1}^{2}+\frac{(u_{1l}+u_{2l})}{2}\epsilon_{4}\|E^{A}_{c}\|^{2}\\\
&\leq\frac{(u_{1l}+u_{2l})}{2\epsilon_{4}}h^{2}\|c\|_{2}^{2}+\frac{(u_{1l}+u_{2l})}{2}\epsilon_{4}\|E^{A}_{c}\|^{2}\end{split}$
(24)

Similarly

$\begin{split}T_{5}&=-\int_{\Omega}(u_{1}\frac{\partial E^{A}_{c}}{\partial
x}+u_{2}\frac{\partial E^{A}_{c}}{\partial y})E^{A}_{c}\\\
&\leq\frac{(u_{1l}+u_{2l})}{2\epsilon_{4}}(\|\frac{\partial
E^{A}_{c}}{\partial x}\|^{2}+\|\frac{\partial E^{A}_{c}}{\partial
y}\|^{2})+\frac{(u_{1l}+u_{2l})}{2}\epsilon_{4}\|E^{A}_{c}\|^{2}\end{split}$
(25)

Now combining all the results obtained above into the equation (20) and
proceeding as explained at the end of the first part we finally get

$\|\frac{\partial E^{A}_{c}}{\partial x}\|^{2}+\|\frac{\partial
E^{A}_{c}}{\partial y}\|^{2}+\|E^{A}_{c}\|^{2}\leq C_{2}(c)h^{2}$ (26)

This completes the proof. ∎

######  Theorem 2.

Apriori error estimate: Assuming the same condition as in the previous theorem

$\begin{split}\|u_{1}-u_{1h}\|^{2}_{V}+\|u_{2}-u_{2h}\|^{2}_{V}+\|p-p_{h}\|_{Q}^{2}+\|c-c_{h}\|_{V}^{2}&\leq
C^{\prime}(\textbf{u},p,c)h^{2}\\\ \end{split}$ (27)

where $C^{\prime}$ is a constant depending upon $\textbf{u},p,c$.

###### Proof.

We prove this result by first applying triangle inequality and then
interpolation estimates and the results obtained in the previous theorem as
the following

$\|u_{1}-u_{1h}\|^{2}_{V}+\|u_{2}-u_{2h}\|^{2}_{V}+\|p-p_{h}\|_{Q}^{2}+\|c-c_{h}\|_{V}^{2}\\\
=\|E^{I}_{u1}+E^{A}_{u1}\|_{V}^{2}+\|E^{I}_{u2}+E^{A}_{u2}\|_{V}^{2}+\|E^{I}_{c}+E^{A}_{c}\|_{V}^{2}+\|E^{I}_{c}+E^{A}_{c}\|_{V}^{2}\hskip
56.9055pt\\\
\leq\bar{C}(\|E^{I}_{u1}\|_{V}^{2}+\|E^{I}_{u2}\|_{V}^{2}+\|E^{I}_{p}\|_{V}^{2}+\|E^{I}_{c}\|_{V}^{2}+\|E^{A}_{u1}\|_{V}^{2}+\|E^{A}_{u2}\|_{V}^{2}+\|E^{A}_{p}\|_{V}^{2}+\|E^{A}_{c}\|_{V}^{2})\\\
\leq C^{\prime}(\textbf{u},p,c)h^{2}\hskip 284.52756pt$ (28)

This completes the proof. ∎

######  Remark 1.

This error estimation is not computable as it depends upon the true solution,
but it implies convergence of the stabilized method and the order of
convergence for this method is 2.

##  4 Numerical Experiment

In this section through a test case based on hydrological importance [10], we
have shown subgrid multiscale stabilization method is performing better than
Galerkin method for small diffusion for this weakly coupled Stokes-Darcy/
transport equation. We have considered a source of pollutant dispersed into an
incompressible fluid in a simple bounded square domain
$\Omega=(0,1)\times(0,1)$. The advection velocity comes from the solution of
Stokes-Darcy flow problem and the diffusion coefficients are taken as
$D_{1}=10^{-7}(1+0.02x)^{2}$, $D_{2}=10^{-8}(1+0.02y)^{2}$ and the reaction
term $\alpha=10$. Table 1 represents the comparison between Galerkin method
and SGS method for small diffusion.
Again for higher diffusion coefficients compared to the reaction term both the
Galerkin method and SGS method perform well. Let us consider the values of
diffusion coefficients as $D_{1}=(1+0.02x)^{2}$, $D_{2}=0.1(1+0.02y)^{2}$ and
reaction term $\alpha=0.001$. Table 2 presents the comparison of both Galerkin
method and SGS method for diffusion dominated case.

| Galerkin Method | Galerkin Method | SGS Method | SGS Method
---|---|---|---|---
Mesh size | Error | Order of convergence | Error | Order of convergence
10 | 0.00049606 |  | 0.000265901 |
20 | 0.000152662 | 1.70017 | 6.36986 $e^{-5}$ | 2.06155
40 | 3.94719 $e^{-5}$ | 1.95144 | 1.68468 $e^{-5}$ | 1.91879
80 | 1.43178 $e^{-5}$ | 1.46302 | 4.73744 $e^{-6}$ | 1.83029
160 | 4.79682 $e^{-6}$ | 1.57766 | 5.4902 $e^{-7}$ | 1.91276
320 | 1.38697 $e^{-6}$ | 1.79015 | 4.45223 $e^{-8}$ | 1.89875
Table 1:  Comparison between Galerkin method and SGS method for small diffusion | Galerkin Method | Galerkin Method | SGS Method | SGS Method
---|---|---|---|---
Mesh size | Error | Order of convergence | Error | Order of convergence
10 | 0.000571444 |  | 0.000574645 |
20 | 0.000139961 | 2.02959 | 0.000141085 | 2.0261
40 | 3.20151 $e^{-5}$ | 2.1282 | 3.22881 $e^{-5}$ | 2.12749
80 | 7.74389 $e^{-6}$ | 2.04762 | 7.82478 $e^{-6}$ | 2.04488
160 | 2.02045 $e^{-6}$ | 1.93838 | 2.03639 $e^{-6}$ | 1.94204
320 | 4.69837 $e^{-7}$ | 2.10444 | 4.75272 $e^{-7}$ | 2.09919
Table 2:  Comparison between Galerkin method and SGS method for diffusion
dominated case

######  Remark 2.

The numerical experiment has proved theoretically established result and shown
that the order of convergence under SGS method for both the cases is 2 whereas
for small diffusion Galerkin method oscillates and for diffusion dominated
flow it behaves well as stabilized method.

##  5 Conclusion

In this paper we have studied subgrid multiscale stabilization formulation of
weakly coupled Stokes-Darcy and transport equation. As well as we have carried
out apriori error estimation, which implies the convergence of the method and
finally the numerical experiment verifies the theoretically established
result.

## References

  * [1] D. Vassilev, I. Yotov, Coupling Stokes-Darcy flow with transport, SIAM J. Sci. Comput.(2009) 3661-3684.
  * [2] A. Cesmelioglu, P. Chidyagwai, B. Rivi$\grave{e}$re, Continuous and discontinuous finite element methods for coupled surface-subsurface flow and transport problems.
  * [3] G. Houzeaux, J. Principe, A variational subgrid scale model for transient incompressible flows, International Journal of Computational Fluid Dynamics. (2008) 135-152.
  * [4] M. Chowdhury, B.V.R. Kumar, On subgrid multiscale stabilized finite element method for advection-diffusion-reaction equation with variable coefficients.(pre-prints), arXiv: 1812.06316 [math.AP] 15 Dec. 2018.
  * [5] R. Codina, A stabilized finite element method for generalized stationary incompressible flows, Comput. Methods Appl. Mech. Engrg. 190 (2001) 2681-2706.
  * [6] B.V.R. Kumar, M. Chowdhury, A priori and a posteriori error estimation for finite element approximation of advection-diffusion-reaction equation with spatially variable coefficients.(pre-prints), arXiv: 1811.05283v1 [math.AP] 13 Nov. 2018.
  * [7] P.G. Ciarlet, The Finite Element Method for Elliptic Problems (North-Holland, amsterdam, 1978).
  * [8] Cl. Strang and G. Fix, An Analysis of the Finite Element Method (Prentice-Hall, Englewood Cliffs, NJ, 1973).
  * [9] S. Badia and R. Codina, Unified stabilized finite element formulations for the Stokes and the Darcy problems, SIAM Journal on Numerical Analysis,Vol. 47 (2009), 1971-2000
  * [10] A.Sanskrityayn, V.P.Singh, V.K.Bharati and N.Kumar, Analytical solution of two-dimensional advection-dispersion equation with spatio-temporal coefficients for point sources in an infinite medium using Green’s function method, Environ. Fluid Mech. 18(2018) 739-757.


* [DLZ12] Ugo Dal Lago and Margherita Zorzi.  Probabilistic operational semantics for the lambda calculus.  RAIRO–Theoretical Informatics and Applications, 46(3):413–450, 2012.  doi:10.1051/ita/2012012.
  * [DP95] Ugo Deliguoro and Adolfo Piperno.  Nondeterministic extensions of untyped $\lambda$-calculus.  Information and Computation, 122(2):149–177, 1995.  doi:10.1006/inco.1995.1145.
  * [DPHW05] Alessandra Di Pierro, Chris Hankin, and Herbert Wiklicky.  Probabilistic lambda-calculus and quantitative program analysis.  Journal of Logic and Computation, 15(2):159–179, 2005.  doi:10.1093/logcom/exi008.
  * [Fil03] Jean-Christoph Filliâtre.  Verification of non-functional programs using interpretations in type theory.  Journal of Functional Programming, 13(4):709–745, 2003.  doi:10.1017/S095679680200446X.
  * [GK05] Philipp Gerhardy and Ulrich Kohlenbach.  Extracting Herbrand disjunctions by functional interpretation.  Archive for Mathematical Logic, 44:633–644, 2005.  doi:10.1007/s00153-005-0275-1.
  * [Hoa69] Charles Antony Richard Hoare.  An axiomatic basis for computer programming.  Communications of the ACM, 12(10):576–580, 1969.  doi:10.1145/363235.363259.
  * [JP89] Claire Jones and Gordon Plotkin.  A probabilistic powerdomain of evaluations.  In Proceedings of Logic in Computer Science (LICS’89), pages 186–195. IEEE Press, 1989.
  * [Kre59] Georg Kreisel.  Interpretation of analysis by means of functionals of finite type.  In A. Heyting, editor, Constructivity in Mathematics, pages 101–128. North-Holland, Amsterdam, 1959.
  * [Kri09] Jean-Louis Krivine.  Realizability in classical logic in interactive models of computation and program behaviour.  Panoramas et synthéses, 27, 2009.
  * [McC10] Guy McCusker.  A graph model for imperative computation.  Logical Methods in Computer Science, 6(1:2), 2010.  doi:doi.org/10.2168/LMCS-6(1:2)2010.
  * [NAMB07] Aleksandar Nanevski, Amal Ahmed, Greg Morrisett, and Lars Birkedal.  Abstract predicates and mutable ADTs in Hoare type theory.  In Proceedings of the European Symposium on Programming (ESOP’07), volume 4421 of LNCS, pages 189–204, 2007.  doi:10.1007/978-3-540-71316-6_14.
  * [OP15a] Paulo Oliva and Thomas Powell.  A constructive interpretation of Ramsey’s theorem via the product of selection functions.  Mathematical Structure in Computer Science, 25(8):1755–1778, 2015\.  doi:10.1017/S0960129513000340.
  * [OP15b] Paulo Oliva and Thomas Powell.  A game-theoretic computational interpretation of proofs in classical analysis.  In Reinhard Kahle and Michael Rathjen, editors, Gentzen’s Centenary: The Quest for Consistency, pages 501–531. Springer, 2015.  doi:10.1007/978-3-319-10103-3_18.
  * [PCW05] Iman Poernomo, John N. Crossley, and Martin Wirsing.  Adapting Proofs-as-Programs.  Monographs in Computer Science. Springer, 2005.  doi:10.1007/0-387-28183-5.
  * [Pow16] Thomas Powell.  Gödel’s functional interpretation and the concept of learning.  In Proceedings of Logic in Computer Science (LICS ’16), pages 136–145. ACM, 2016.  doi:10.1145/2933575.2933605.
  * [Pow18] Thomas Powell.  A functional interpretation with state.  In Proceedings of Logic in Computer Science (LICS ’18), pages 839–848. ACM, 2018.  doi:10.1145/3209108.3209134.
  * [Pow20] Thomas Powell.  Well quasi-orders and the functional interpretation.  In Peter Schuster, Monika Seisenberger, and Andreas Weiermann, editors, Well Quasi-Orders in Computation, Logic, Language and Reasoning, volume 53 of Trends in Logic, pages 221–269. Springer, 2020\.  doi:10.1007/978-3-030-30229-0_9.
  * [PSW22] Thomas Powell, Peter Schuster, and Franziskus Wiesnet.  A universal algorithm for krull’s theorem.  Information and Computation, 287:104761, 2022.  doi:10.1016/j.ic.2021.104761.
  * [Red96] Uday Reddy.  Global state considered unnecessary.  Lisp and Symbolic Computation, 9:7–76, 1996.
  * [Rey02] John C. Reynolds.  Separation logic: a logic for shared mutable data structures.  In Proceedings of Logic in Computer Science (LICS’02), pages 55–74. IEEE, 2002.  doi:10.1109/LICS.2002.1029817.
  * [SD78] Nasser Saheb-Djaromi.  Probabilistic LCF.  In Proceedings of International Symposium on Mathematical Foundations of Computer Science (MFCS’78), volume 64 of LNCS, pages 442–451, 1978.
  * [Spe62] Clifford Spector.  Provably recursive functionals of analysis: a consistency proof of analysis by an extension of principles in current intuitionistic mathematics.  In F. D. E. Dekker, editor, Recursive Function Theory: Proc. Symposia in Pure Mathematics, volume 5, pages 1–27. American Mathematical Society, 1962.
  * [Tro73] Anne S. Troelstra.  Metamathematical Investigation of Intuitionistic Arithmetic and Analysis, volume 344 of Lecture Notes in Mathematics.  Springer-Verlag, 1973.  doi:10.1007/BFb0066739.
  * [VVB04] Joost Vennekens, Sofie Verbaeten, and Maurice Brunooghe.  Logic programs with annotated disjunctions.  In ICLP 2004: Logic Programming, volume 3132 of Lecture Notes in Computer Science, pages 431–445, 2004.  doi:10.1007/978-3-540-27775-0_30.
  * [YHB07] Nobuko Yoshida, Kohei Honda, and Martin Berger.  Logical reasoning for higher-order functions with local state.  In Proceedings of Foundations of Software Science and Computation Structures (FOSSACS’07), volume 44213 of LNCS, pages 361–377, 2007.  doi:10.1007/978-3-540-71389-0_26.