# Exploring Well-Posedness and Asymptotic Behavior in an Advection-Diffusion-
Reaction (ADR) Model

Mohammed Elghandouri<EMAIL_ADDRESS>Khalil Ezzinbi
<EMAIL_ADDRESS>Lamiae Saidi<EMAIL_ADDRESS>
###### Abstract

In this paper, the existence, uniqueness, and positivity of solutions, as well
as the asymptotic behavior through a finite fractal dimensional global
attractor for a general Advection-Diffusion-Reaction (ADR) equation, are
investigated. Our findings are innovative, as we employ semigroups and global
attractors theories to achieve these results. Also, an analytical solution of
a two-dimensional Advection-Diffusion Equation is presented. And finally, two
Explicit Finite Difference schemes are used to simulate solutions in the two-
and three-dimensional cases. The numerical simulations are conducted with
predefined initial and Dirichlet boundary conditions.

###### keywords:

Advection–Diffusion–Reaction, Partial Differential Equation, Semigroups
Theory, Existence problems, Mild solutions, Global Attractor, Fractal
dimension, Numerical methods.

###### MSC:

35K57, 47D60, 35A01, 34D45, 65C20.

††journal: Computational and Applied Mathematics

[a]organization=Cadi Ayyad University, Faculty of Sciences Semlalia,
addressline=Departement of Mathematics, postcode=B.P: 2390, city=Marrakesh,
country=Morocco [b]organization=Sorbonne University, EDITE (ED130),
addressline=IRD, UMMISCO, postcode=F-93143, city=Bondy, Paris, country=France

##  1 Introduction

The study of Partial Differential Equations (PDEs) holds a paramount position
in the realm of mathematical analysis, finding applications across a diverse
range of scientific disciplines and engineering fields. Among these, the
Advection–Diffusion–Reaction (ADR) partial differential equations, that have
widespread applications in fluid dynamics, heat or mass transfer, chemical
reaction processes, contaminant transport and population dynamics [2, 9, 12,
13, 19, 25, 27]. They model the temporal evolution of biological species in a
flowing medium such as water or air, or contaminant transport with biological,
chemical or radioactive processes, accounting for advection, diffusion and
reaction processes. Therefore, there exists a wide literature about its
resolution, numerically and analytically, surrounding the existence,
uniqueness, asymptotic behavior, controllability, numerical approximation of
its solutions, etc. In [28], a distinctive classical solution for the single-
step irreversible exothermic Arrhenius-type reaction within an incompressible
fluid subject to Neumann boundary conditions is demonstrated. The authors
provide both proof of the solution’s existence and its uniqueness through the
utilization of a semigroup formulation as well as the application of the
maximum principle. In [22], the authors study the existence and uniqueness of
mild solutions for a fractional reaction-diffusion equation, namely fractional
Fokker–Planck equations. They use semigroup theory to establish their results.
In [17], the authors establish the well-posedness of an initial-boundary value
problem encompassing a broad range of linear time-fractional advection-
diffusion-reaction equations. The analysis hinges on innovative energy
techniques, coupled with the application of a fractional Gronwall inequality
and characteristics of fractional integrals. In [7], the authors explore the
precise characteristics of controllability exhibited by an advection–diffusion
equation within a bounded domain. They achieve this investigation by utilizing
time- and space-dependent velocity fields as controlling factors. Other
relevant works are mentioned in [15, 23, 24, 29, 30]. Nevertheless, it’s
important to note that, in terms of mathematical advancements such as
asymptotic behavior, existence, uniqueness, and positiveness of solutions for
general ADR equations, there hasn’t been significant development, particularly
concerning semigroup methods.

Numerically, there are various methods to solve different partial differential
equations (PDEs) considering initial and boundary conditions like Dirichlet,
Neumann, or Robin. A brief list of some of the numerical methods are: explicit
methods (Finite Difference Methods [5, 10, 18] (Explicit Forward Euler method,
Upwind scheme, the Central Finite Difference scheme, Total Variation
Diminishing (TVD) Schemes [26], etc.) and method of lines); Implicit methods
(implicit backward Euler method, Crank Nicholson, Method of Lines); Iterative
methods (Jacobi, Gauss–Seidel, and relaxation methods); and Finite Element
Methods (Galerkin Method) [16]. However, the numerical resolution of a
nonlinear reaction-advection-diffusion equation while maintaining certain
global physical properties such as mass conservation, positivity, and
causality is complex. The goal is for the numerical solutions to satisfy the
same properties as the exact solution such as positivity, boundedness, or
monotonicity. Most numerical schemes aren’t inherently designed to satisfy
these properties, and they can lead to various numerical artifacts such as
artificial oscillations and numerical dispersion. Explicit schemes are
efficient but suffer from stability issues and require small time steps to
meet the Courant-Friedrichs-Lewy (CFL) condition, thus the overall
computational cost is usually high. Implicit schemes are unconditionally
stable and allow larger time steps but involve solving algebraic systems at
each time step. Semi-implicit schemes strike a balance between explicit and
implicit methods. The choice of which scheme to use depends on the specific
problem and requirements. Implicit schemes are often preferred for stiff
initial value problems due to their stability and larger allowable time steps.
Nevertheless, among the various numerical methods for solving ADR equation,
finite difference method (FDM) seems to be more popular for the ease of
implementation and its simplicity.

Our primary contribution in this work is the comprehensive investigation of
the global existence, uniqueness, and positivity of solutions, along with
exploring the asymptotic behavior of an ADR equation through the application
of semigroup theory. As a result, we establish the existence of a global
attractor while also determining its finite fractal dimension. To implement
our theoretical findings numerically, we employ two fully explicit finite
difference numerical schemes for solving the ADR equation: an explicit
centered scheme in 2-D and an explicit upwind scheme in 3-D.

The organization of this paper is as follows: In Section 2, we present the
model formulation, its origins, and the associated problems. In Section 3, we
review useful concepts that will assist in proving our main results. Section 4
employs semigroup theory to demonstrate the main findings related to the
global existence, uniqueness, and positivity of solutions for our model. The
existence of a global attractor with a finite fractal dimension is established
in Section 5. Section 6 is dedicated to the analytical solution of the
advection-diffusion equation in 2-D, while Section 7 focuses on numerical
simulations using finite difference methods.

##  2 Formulation of the model

###  2.1 Origins and Problematic

The ADR model, which describes the time evolution of chemical species in air,
is characterized by partial differential equations (PDEs) derived from mass
balances. Consider a concentration $c(t,x)$ of a certain chemical species,
with space variable $x\in\Omega\subset\mathbb{R}$ and time $t\geq 0$. Let
$h>0$ be a small number, and consider the average concentration $\bar{c}(t,x)$
in a cell $\Omega(x)=\left[x-\frac{1}{2}h,x+\frac{1}{2}h\right]$,

$\bar{c}(t,x)=\frac{1}{h}\int_{x-h/2}^{x+h/2}c\left(t,y\right)dy=c(t,x)+\frac{1}{24}h^{2}c_{xx}(t,x)+\cdots.$

If the species is carried along by a flowing medium with velocity $u(t,x)$
then the mass conservation law implies that the change of $\bar{c}(t,x)$ per
unit of time is the net balance of inflow and outflow over the cell
boundaries,

$\frac{\partial}{\partial
t}\bar{c}(t,x)=\frac{1}{h}\left[u\left(t,x-\frac{1}{2}h\right)c\left(t,x-\frac{1}{2}h\right)-u\left(t,x+\frac{1}{2}h\right)c\left(t,x+\frac{1}{2}h\right)\right],$

where $u\left(t,x\pm\frac{1}{2}h\right)c\left(t,x\pm\frac{1}{2}h\right)$ are
the mass fluxes over the left and right cell boundaries. Now, if we let
$h\rightarrow 0$, it follows that the concentration satisfies

$\frac{\partial}{\partial t}c(t,x)+\frac{\partial}{\partial
x}(u(t,x)c(t,x))=0.$

This is called an advection equation (or convection equation). In a similar
way, we can consider the effect of diffusion. The change in $\bar{c}(t,x)$ is
induced by gradients in the solution, and the fluxes across the cell
boundaries are
$-k\left(t,x\pm\frac{1}{2}h\right)c_{x}\left(t,x\pm\frac{1}{2}h\right)$, where
$k(t,x)$ represents the diffusion coefficient. The corresponding diffusion
equation is

$\frac{\partial}{\partial t}c(t,x)=\frac{\partial}{\partial
x}\left(k(t,x)\frac{\partial}{\partial x}c(t,x)\right).$

There may also be a change in $c(t,x)$ due to sources, sinks and chemical
reactions, leading to

$\frac{\partial}{\partial t}c(t,x)=f(c(t,x),t,x).$

The overall change in concentration is described by combining these three
effects, leading to the advection-diffusion-reaction equation:

$\frac{\partial}{\partial t}c(t,x)+\frac{\partial}{\partial
x}(u(t,x)c(t,x))=\frac{\partial}{\partial
x}\left(k(t,x)\frac{\partial}{\partial x}c(t,x)\right)+f(c(t,x),t,x).$

We shall consider the equation in a spatial domain
$\Omega\subset\mathbb{R}^{m}$ ($m=1,2\text{ or }3$) with time $t\geq 0$. An
initial profile $c(0,x)$ will be given and we also assume that suitable
boundary conditions are provided. More general, let
$c_{1}(t,x),\cdots,c_{s}(t,x)$ be concentrations of $s$ chemical species, with
spatial variable $x\in\Omega$, and time $t\geq 0$. Then the basic mathematical
equations for transport and reaction are given by the following set of (PDEs):

$\begin{gathered}\frac{\partial}{\partial
t}c_{j}(t,x)+\sum_{i=1}^{m}\frac{\partial}{\partial
x_{i}}\left(u_{i}(t,x)c_{j}(t,x)\right)=\\\
\sum_{i=1}^{m}\frac{\partial}{\partial
x_{i}}\left(k_{i}(t,x)\frac{\partial}{\partial
x_{i}}c_{j}(t,x)\right)\\!+\\!f_{j}\left(c_{1}(x,t),\ldots,c_{s}(x,t),t,x\right),\quad
j=1,2,\ldots,s,\end{gathered}$

with suitable initial and boundary conditions. The quantities $u_{i}$ that
represent the velocities of the transport medium, such as water or air, are
either given in a data archive or computed alongside with a meteorological or
hydrodynamical code. (In such codes Navier-Stokes or shallow water equations
are solved, where again advection-diffusion equations are of primary
importance.) The diffusion coefficients $k_{i}$ are constructed by the
modelers and may include also parameterizations of turbulence. The final term
$f_{j}(c,x,t)$, which gives a coupling between the various species, describes
the nonlinear chemistry together with emissions (sources) and depositions
(sinks). In actual models these equations are augmented with other suitable
subgrid parameterizations and coordinate transformations.

###  2.2 Description of the model

In this work, we undertake a comprehensive exploration of the following ADR
equation:

$\dfrac{\partial c_{j}(t,x)}{\partial
t}=\sum\limits_{i=1}^{m}\left(k_{i}\dfrac{\partial^{2}c_{j}(t,x)}{\partial
x_{i}^{2}}-u_{i}\dfrac{\partial c_{j}(t,x)}{\partial
x_{i}}\right)+R_{j}(t,c(t,x)),$ (1)

for $t\geq 0$, $x=(x_{1},\ldots,x_{m})\in\Omega$, where $\Omega$ is a bounded
set in $\mathbb{R}^{m}$, $m\in\\{1,2,3\\}$, and $j=1,\ldots,s$ with
$s\in\mathbb{N}^{*}$. The coefficients $k_{i}$ and $u_{i}$ are positive
constants, and $c(t,x)=(c_{1}(t,x),\ldots,c_{s}(t,x))$. The function
$R_{j}(t,c(t,x))$ represents the $j^{th}$ chemical reaction between species
and is given by

$R_{j}(t,c(t,x))=\sum\limits_{\kappa=1}^{r}(r_{j\kappa}-l_{j\kappa})h_{\kappa}(t)\prod\limits_{\nu=1}^{s}(c_{\nu}(t,x))^{l_{\nu\kappa}},\hskip
2.84544ptj=1,\ldots,s,\hskip 2.84544ptx\in\Omega,\hskip 2.84544ptt\geq 0.$ (2)

Here, $r\in\mathbb{N}^{*}$ represents the number of chemical reactions between
species, $l_{j\kappa}$ and $r_{j\kappa}$ are non-negative integers describing
the loss and gain of the number of molecules in the $j^{th}$ reaction, and
$h_{\kappa}(t)$ is the rate of reaction, dependent on time $t$ due to external
influences such as temperature and sunlight. Importantly, for each
$\kappa=1,\ldots,r$, $h_{\kappa}(\cdot)$ is almost everywhere continuous, and
there exists $d_{\kappa}>0$ such that $0\leq h_{\kappa}(t)\leq d_{\kappa}$ for
$t\geq 0$. On the boundary of $\Omega$, we use Dirichlet boundary conditions,
i.e.,

$c_{j}(t,x)=0,\hskip 2.84544ptj=1,\ldots,s\hskip 2.84544pt\text{ for }\hskip
2.84544ptx\in\partial\Omega\hskip 2.84544pt\text{ and }\hskip 2.84544ptt\geq
0.$ (3)

As for the initial conditions, we denote them by

$c_{j}(0,x)=c_{j,0}(x),\hskip 2.84544ptj=1,\ldots,s\hskip 2.84544pt\text{ for
}\hskip 2.84544ptx\in\Omega,$ (4)

representing the initial concentrations for the species.

Let $\mathcal{A}:\left(H_{0}^{1}(\Omega)\cap
H^{2}(\Omega)\right)^{s}\to\left(\mathbb{L}^{2}(\Omega)\right)^{s}$ be the map
defined by

$\left(\mathcal{A}z\right)(x)=\sum\limits_{i=1}^{m}\left(k_{i}\dfrac{\partial^{2}z(x)}{\partial
x_{i}^{2}}-u_{i}\dfrac{\partial z(x)}{\partial x_{i}}\right)\text{ for
}z\in\left(H_{0}^{1}(\Omega)\cap H^{2}(\Omega)\right)^{s},\hskip
2.84544ptx\in\Omega.$

We denote by
$u(t)(x)=(u_{1}(t)(x),\ldots,u_{s}(t)(x)):=(c_{1}(t,x),\ldots,c_{s}(t,x))=c(t,x)$,
and let $F(t,u(t))(x)=\left(F_{1}(t,u(t))(x),\ldots,F_{s}(t,u(t))(x)\right)$
for  $t\geq 0$, where

$F_{j}(t,u(t))(x)=R_{j}(t,u(t)(x)),\quad j=1,\ldots,s,\quad x\in\Omega,$ (5)

for $u(t)=(u_{1}(t),\ldots,u_{s}(t))\in D(F)$, where

$D(F):=\left\\{z=(z_{1},\ldots,z_{s})\in\left(\mathbb{L}^{2}(\Omega)\right)^{s}\mid
F(t,z)\in\left(\mathbb{L}^{2}(\Omega)\right)^{s}\text{ for all }t\geq
0\right\\}.$ (6)

Under the above notations, equations (1)–(4) can be expressed in the following
abstract form:

$\left\\{\begin{array}[]{l}u^{\prime}(t)=\mathcal{A}u(t)+F(t,u(t)),\quad t\geq
0\\\ u(0)=u_{0}.\end{array}\right.$ (7)

We then delve into the mathematical properties of equation (7), investigating
conditions that ensure the existence, uniqueness and positivity of solutions,
as well as exploring the asymptotic behavior through a finite fractal-
dimensional global attractor. Additionally, we develop numerical strategies
for effective approximation. Numerical simulations are conducted using the
explicit forward Euler method for the three-dimensional advection-diffusion-
reaction equations with Dirichlet boundary conditions set to zero. Our focus
is on an advection-dominant problem related to the transport of air
pollutants, providing valuable insights into the practical implications of our
theoretical framework.

##  3 Basic working tools

In this section, some necessary definitions, theorems, and lemmas to
demonstrate our main results are recalled. Let $(X,\|\cdot\|_{X})$ be a real
Banach lattice space endowed with the ordering $\geqq$. We define $X_{+}$ as
the set of all elements $x$ in $X$ such that $x\geqq 0$; this set is called
the positive cone of $X$. Additionally, we denote $x^{+}=\max(0,x)$ for $x\in
X$.

###### Definition 3.1

[1, Definition 11.9] A linear operator $A$ on $X$ is called dispersive if for
every $f\in D(A)$, and $\lambda>0$

$\|(\lambda f-Af)^{+}\|_{X}\geq\lambda\|f^{+}\|_{X}$.

The following theorem shows the dispersivity of linear operators in Hilbert
spaces.

###### Theorem 3.2

[1] A linear operator $A$ on a real Hilbert lattice space $H$ is dispersive if
and only if $\langle Af,f^{+}\rangle_{H}\leq 0$ for every $f\in D(A)$.

We state the following definition of $A$-boundedness.

###### Definition 3.3

[1] Let $A$ be a closed linear operator on $X$. A linear operator $B$ with
domain $D(B)$ is called $A$-bounded if $D(A)\subseteq D(B)$ and there exist
$a,b\geq 0$ such that for all $f\in D(A)$ the inequality

$\|Bf\|_{X}\leq a\|Af\|_{X}+b\|f\|_{X},$ (8)

holds. The $A$-bound of $B$ is defined by

$a_{0}=\inf\left\\{a\geq 0\mid\text{there exists $b\geq 0$ such that the
inequality \eqref{In 6} holds}\right\\}$.

The following Theorem will be needed to demonstrate that the linear operator
$\mathcal{A}$ generates a positive $C_{0}$-semigroup of contraction on
$\left(\mathbb{L}^{2}(\Omega)\right)^{s}$.

###### Theorem 3.4

[1, Theorem 13.3.] Let $A$ be the generator of a positive $C_{0}$-semigroup of
contraction on $X$ and $B$ a dispersive and $A$-bounded operator with
$A$-bound $a_{0}<1$. Then $A+B$, defined on $D(A)$, generates a positive
$C_{0}$-semigroup of contraction on $X$.

We recall the following definition of mild solutions for equation (7).

###### Definition 3.5

Let $0<a\leq+\infty$. A function
$u(\cdot):[0,a)\to\left(\mathbb{L}^{2}(\Omega)\right)^{s}$ is called a mild
solution of equation (7) on $[0,a)$ if

$u(t)=S(t)u_{0}+\displaystyle\int_{0}^{t}S(t-\tau)F(\tau,u(\tau))d\tau\text{
for }t\in[0,a),$ (9)

where $(S(t))_{t\geq 0}$ is the $C_{0}$-semigroup generated by $\mathcal{A}$
on $\left(\mathbb{L}^{2}(\Omega)\right)^{s}$.

To demonstrate the existence of solutions for equation (7), we require the
following theorem.

###### Theorem 3.6

[21, Theorem 1.4] Let $f:\mathbb{R}^{+}\times X\to X$ be continuous w.r.t the
first argument and locally Lipschitz continuous w.r.t the second argument. If
$B$ is the infinitesimal generator of a $C_{0}$-semigroup $(T(t))_{t\geq 0}$
on $X$, then for every $v_{0}\in X$, there exists a
$t_{\text{max}}\leq+\infty$ such that the initial value problem:

$\left\\{\begin{array}[]{l}v^{\prime}(t)=Bv(t)+f(t,v(t)),\quad t\geq 0\\\
v(0)=v_{0},\end{array}\right.$ (10)

has a unique mild solution $v$ on $[0,t_{\text{max}}[$. Moreover, if
$t_{\text{max}}<+\infty$, then $\limsup\limits_{t\to
t_{\text{max}}}\|v(t)\|_{X}=+\infty$.

###### Corollary 3.7

[1] Assume that conditions in Theorem 3.6 are satisfied. If $f$ is at most
affine w.r.t the second argument, then $t_{\text{max}}=+\infty$.

###### Remark 3.8

Based on the proofs of Theorems 1.4 and 1.2 in [21], it follows that if $f$ is
almost everywhere continuous and bounded with respect to the first argument,
and locally Lipschitz continuous with respect to the second argument, then the
conclusions of Theorem 3.6 still hold.

The following Lemmas will be needed throughout this work.

###### Lemma 3.9

[1, Proposition A.47] Let $I\subset\mathbb{R}$ be an open interval and $f\in
W^{2,p}_{0}(I;\mathbb{R})$, $p\geq 1$. Then, for any $\epsilon>0$ there is a
constant $C(\epsilon)>0$ such that

$\|f^{{}^{\prime}}\|_{\mathbb{L}^{p}}\leq\epsilon\|f^{{}^{\prime\prime}}\|_{\mathbb{L}^{p}}+C(\epsilon)\|f\|_{\mathbb{L}^{p}}$.

###### Lemma 3.10

[3] Let $v$ be a real, continuous and nonnegative function such that

$v(t)\leq c+\displaystyle\int_{t_{0}}^{t}w(t,\tau)v(\tau)d\tau$ for  $t\geq
t_{0}$,

where $c>0$, $w(t,\tau)$ is continuously differentiable in $t$ and continuous
in $\tau$ with $w(t,\tau)\geq 0$ for $t\geq\tau\geq t_{0}$. Then,

$v(t)\leq
c\exp\left(\displaystyle\int_{t_{0}}^{t}\left[w(\tau,\tau)+\displaystyle\int_{t_{0}}^{\tau}\dfrac{\partial
w(\tau,r)}{\partial\tau}dr\right]d\tau\right)$ for  $t\geq t_{0}$.

##  4 Well-posedness

In this section, we study the existence, uniqueness, and positiveness of the
solution for equation (7). The Hilbert lattice space
$\left(\mathbb{L}^{2}(\Omega)\right)^{s}$ is equipped with the norm
$\|\cdot\|$ defined by

$\|z\|=\left(\sum\limits_{j=1}^{s}\|z_{j}\|^{2}_{\mathbb{L}^{2}(\Omega)}\right)^{\frac{1}{2}}$
for  $z\in\left(\mathbb{L}^{2}(\Omega)\right)^{s}$.

We denote by
$\left(\mathbb{L}^{2}(\Omega)\right)^{s}_{+}=\left\\{z\in\left(\mathbb{L}^{2}(\Omega)\right)^{s}\mid
z_{j}\geq 0,\hskip 2.84544ptj=1,\ldots,s\right\\}$ the positive cone of the
space $\left(\mathbb{L}^{2}(\Omega)\right)^{s}$.

###### Proposition 4.1

$\mathcal{A}$ is an infinitesimal generator of a positive $C_{0}$-semigroup
$(S(t))_{t\geq 0}$ of contraction on
$\left(\mathbb{L}^{2}(\Omega)\right)^{s}$.

###### Proof 1

The linear operator $\mathcal{A}$ can be decomposed as follows:

$\mathcal{A}z=\mathcal{A}_{1}z+\mathcal{A}_{2}z$ for
$z\in\left(H_{0}^{1}(\Omega)\cap H^{2}(\Omega)\right)^{s}$,

where

$(\mathcal{A}_{1}z)(x)=\sum\limits_{i=1}^{m}\left(k_{i}\dfrac{\partial^{2}z(x)}{\partial
x_{i}^{2}}\right)$, $z\in\left(H_{0}^{1}(\Omega)\cap
H^{2}(\Omega)\right)^{s}$, $x\in\Omega$,

and

$(\mathcal{A}_{2}z)(x)=\sum\limits_{i=1}^{m}\left(-u_{i}\dfrac{\partial
z(x)}{\partial x_{i}}\right)$, $z\in\left(H_{0}^{1}(\Omega)\cap
H^{2}(\Omega)\right)^{s}$, $x\in\Omega$.

By Lemma 3.9, we show that $\mathcal{A}_{2}$ is $\mathcal{A}_{1}$-bounded
operator with $\mathcal{A}_{1}$-bounded equal to $0$. Let
$z=(z_{1},\ldots,z_{s})\in\left(H_{0}^{1}(\Omega)\cap
H^{2}(\Omega)\right)^{s}$, $z^{+}=(z^{+}_{1},\ldots,z^{+}_{s})$, and
$z^{-}=(z^{-}_{1},\ldots,z^{-}_{s})$ be such that
$z_{j}^{+}(x)=\max(0,z_{j}(x))$ and $z_{j}^{-}(x)=\max(0,-z_{j}(x))$ for
$j=1\ldots,s$, and $x\in\Omega$. Let $z(x)=z^{+}(x)-z^{-}(x)$, then

$\displaystyle\langle\mathcal{A}_{2}z,z^{+}\rangle_{\left(\mathbb{L}^{2}(\Omega)\right)^{s}}$
$\displaystyle=$
$\displaystyle-\sum\limits_{j=1}^{s}\displaystyle\int_{\Omega}\left(\sum\limits_{i=1}^{m}u_{i}\dfrac{\partial
z_{j}(x)}{\partial x_{i}}z_{j}^{+}(x)\right)dx$ $\displaystyle=$
$\displaystyle\dfrac{1}{2}\sum\limits_{j=1}^{s}\sum\limits_{i=1}^{m}u_{i}\displaystyle\int_{\Omega}\dfrac{\partial}{\partial
x_{i}}\left(z_{j}^{+}(x)\right)^{2}dx$ $\displaystyle\leq$
$\displaystyle\dfrac{m}{2}\max(u_{1},\ldots,u_{m})\sum\limits_{j=1}^{s}\left(z_{j}^{+}(x)\right)^{2}_{\mid_{\partial\Omega}}$
$\displaystyle=$ $\displaystyle 0.$

By Theorem 3.2, we obtain that $\mathcal{A}_{2}$ is dispersive. We know that
$\mathcal{A}_{1}$ is an infinitesimal generator of a positive
$C_{0}$-semigroup of contraction on $\left(\mathbb{L}^{2}(\Omega)\right)^{s}$.
Thus, by Theorem 3.4, we obtain that $\mathcal{A}$ is an infinitesimal
generator of a positive $C_{0}$-semigroup $(S(t))_{t\geq 0}$ of contraction on
$\left(\mathbb{L}^{2}(\Omega)\right)^{s}$.

The following theorem demonstrates the local existence and uniqueness of
solutions for equation (7).

###### Theorem 4.2

Equation (7) has a unique mild solution $u(\cdot)$ satisfying equation (9) on
$[0,t_{\text{max}})$ for some $t_{\text{max}}\leq+\infty$.

###### Proof 2

For now, we fix $u_{0}$ in $(C_{0}^{\infty}(\Omega))^{s}$. The function
$F(t,z)$ is not locally Lipschitz continuous in $z$ uniformly with respect to
$t\geq 0$. Therefore, we cannot directly apply Theorem 3.6 to prove our
result. Instead, we employ the truncation function technique and combines it
with the results of Theorem 3.6 and remark 3.8. The truncation procedure is
applied to $F(t,z)$. For the sake of simplicity, consider a function $G(t,y)$
and any positive integer $n$. Define the sets

$\left\\{\begin{array}[]{l}D_{1}^{n}(G)=\\{(t,y)\mid-n\leq G(t,y)\leq
n\\}\\\\[10.0pt] D_{2}^{n}(G)=\\{(t,y)\mid G(t,y)<-n\\}\\\\[10.0pt]
D_{3}^{n}(G)=\\{(t,y)\mid G(t,y)>n\\}.\end{array}\right.$

Further, define the truncation form of the function $G(t,z)$, denoted by
$G^{n}(t,z)$, as

$G^{n}(t,y)=\left\\{\begin{array}[]{l}G(t,y)\quad\text{for }(t,y)\in
D_{1}^{n}(G)\\\\[10.0pt] -n\quad\text{for }(t,y)\in D_{2}^{n}(G)\\\\[10.0pt]
n\quad\text{for }(t,y)\in D_{3}^{n}(G).\end{array}\right.$

Then, the truncation form of $F$ is defined as
$F^{n}(t,z)=\left(F_{1}^{n}(t,z),\ldots,F_{s}^{n}(t,z)\right)$ for
$(t,z)\in\bigcup\limits_{l=1}^{3}D_{l}^{n}(F)$. It is easily seen that
$F^{n}(t,z)$ is Lipschitz continuous in $z$, uniformly with respect to $t\geq
0$. By Theorem 3.6 and Remark 3.8, we obtain that equation:

$\left\\{\begin{array}[]{l}u^{\prime}(t)=\mathcal{A}u(t)+F^{n}(t,u(t)),\quad
t\geq 0\\\\[10.0pt] u(0)=u_{0},\end{array}\right.$ (11)

has a unique mild solution $u^{n}=(u_{1}^{n},\ldots,u_{s}^{n})$ defined on
$[0,\bar{t}_{\text{max}})$ for some $\bar{t}_{\text{max}}\leq+\infty$ by the
following form:

$u^{n}(t)=S(t)u_{0}+\displaystyle\int_{0}^{t}S(t-s)F^{n}(s,u^{n}(s))\,ds$ for
$t\in[0,\bar{t}_{\text{max}})$.

Let

$k_{n}=\max\limits_{1\leq j\leq
s}\left\\{\|F_{j}^{n}(\cdot,u^{n}(\cdot))\|_{\mathbb{L}^{\infty}(\Omega_{m})}\right\\}$
where  $\Omega_{m}=[0,\bar{t}_{\text{max}})\times\Omega$.

Consider the following evolution equation:

$\left\\{\begin{array}[]{l}\bar{u}_{j}^{\prime}(t)=\mathcal{B}\bar{u}_{j}(t)+F_{j}^{n}(t,\bar{u}(t))-k_{n},\quad
t\geq 0\\\\[10.0pt]
\bar{u}_{j}(0)=u_{j,0}-\|u_{j,0}\|_{\mathbb{L}^{\infty}(\Omega)},\quad
u_{j,0}\in\mathcal{C}_{0}^{\infty}(\Omega),\end{array}\right.$ (12)

for $j=1,\ldots,s$, where

$\left(\mathcal{B}z\right)(x)=\sum\limits_{i=1}^{m}\left(k_{i}\dfrac{\partial^{2}z(x)}{\partial
x_{i}^{2}}-u_{i}\dfrac{\partial z(x)}{\partial x_{i}}\right)\text{ for }z\in
H_{0}^{1}(\Omega)\cap H^{2}(\Omega),\hskip 2.84544ptx\in\Omega.$

Since $\left(\mathcal{B},H_{0}^{1}(\Omega)\cap H^{2}(\Omega)\right)$ generates
a positive $C_{0}$-semigroup $(S_{\mathcal{B}}(t))_{t\geq 0}$, it follows that
equation (12) has a unique mild solution given by

$\bar{u}^{n}_{j}(t)=S_{\mathcal{B}}(t)[u_{j,0}-\|u_{j,0}\|_{\mathbb{L}^{\infty}(\Omega)}]+\displaystyle\int_{0}^{t}S_{\mathcal{B}}(t-s)[F_{j}^{n}(t,\bar{u}^{n}(t))-k_{n}]ds$.

Since $u_{j,0}-\|u_{j,0}\|_{\mathbb{L}^{\infty}(\Omega)}\leq 0$ and
$F_{j}^{n}(t,\bar{u}^{n}(t))-k_{n}\leq 0$, it follows that
$\bar{u}^{n}_{j}(t)\leq 0$, means that $\bar{u}_{j}(t,x)\leq 0$ for
$(t,x)\in\Omega_{m}$. Moreover,

$\bar{u}^{n}_{j}(t,x)=u^{n}_{j}(t,x)-k_{n}t-\|u_{j,0}\|_{\mathbb{L}^{\infty}(\Omega)}$
for  $(t,x)\in\Omega_{m}$,

which implies that $u^{n}_{j}(t,x)\leq
k_{n}t+\|u_{j,0}\|_{\mathbb{L}^{\infty}(\Omega)}$ for  $(t,x)\in\Omega_{m}$.
Similarly, we show that equation:

$\left\\{\begin{array}[]{l}\underline{u}_{j}^{\prime}(t)=\mathcal{A}\underline{u}_{j}(t)+F_{j}^{n}(t,\underline{u}(t))+k_{n},\quad
t\geq 0\\\\[10.0pt]
\underline{u}_{j}(0)=u_{j,0}+\|u_{j,0}\|_{\mathbb{L}^{\infty}(\Omega)},\quad
u_{j,0}\in\mathcal{C}_{0}^{\infty}(\Omega),\end{array}\right.$

for $j=1,\ldots,s$, has a unique mild solution satisfying
$\underline{u}_{j}^{n}(t,x)=u_{j}^{n}(t,x)+k_{n}t+\|u_{j,0}\|_{\mathbb{L}^{\infty}(\Omega)}\geq
0$ for  $(t,x)\in\Omega_{m}$, which implies that $u^{n}_{j}(t,x)\geq-
k_{n}t-\|u_{j,0}\|_{\mathbb{L}^{\infty}(\Omega)}$ for $(t,x)\in\Omega_{m}$.
For

$n\geq 2\max\limits_{1\leq j\leq
s}\left\\{\|u_{j,0}\|_{\mathbb{L}^{\infty}(\Omega)}\right\\}$,

fixed, we can find $t_{\text{max}}\in(0,\bar{t}_{\text{max}})$ (it is
sufficient to chose $t_{\text{max}}\in(0,\bar{t}_{\text{max}})$ such that
$t_{\text{max}}k_{n}<\frac{n}{2}$) such that $-n\leq u^{n}_{j}(t,x)\leq n$ for
$(t,x)\in[0,t_{\text{max}})\times\Omega$. Therefore, it follows from the
definition of $F^{n}(t,u^{n}(t))$ that $F^{n}(t,u^{n}(t))=F(t,u(t))$ for
$t\in[0,t_{\text{max}})$. Thus,

$u(t)=S(t)u_{0}+\displaystyle\int_{0}^{t}S(t-s)F(s,u(s))ds$ for
$t\in[0,t_{\text{max}})$,

is a mild solution of equation (7). Now, for
$u_{0}\in(\mathbb{L}^{2}(\Omega))^{s}$, using the fact that
$(C_{0}^{\infty}(\Omega))^{s}$ is dense in $(\mathbb{L}^{2}(\Omega))^{s}$, we
can find a sequence $(u_{0}^{\nu})_{\nu\geq
0}\subset(C_{0}^{\infty}(\Omega))^{s}$ such that $\|u_{0}-u_{0}^{\nu}\|\to 0$
as $\nu\to+\infty$. Following the same approach presented above, we show that
for each $\nu\geq 0$, there exist
$t^{\nu}_{\text{max}}\in(0,\bar{t}^{\nu}_{\text{max}})$ for some
$\bar{t}^{\nu}_{\text{max}}\leq+\infty$ such that
$t^{\nu}_{\text{max}}k_{n(\nu)}<\frac{n(\nu)}{2}$ for

$n(\nu)\geq 2\max\limits_{1\leq j\leq
s}\left\\{\|u^{\nu}_{j,0}\|_{\mathbb{L}^{\infty}(\Omega)}\right\\}$,
$n(\nu)\in\mathbb{N}^{*}$,

and

$k_{n(\nu)}=\max\limits_{1\leq j\leq
s}\left\\{\|F_{j}^{n(\nu)}(\cdot,u^{n(\nu)}(\cdot))\|_{\mathbb{L}^{\infty}(\Omega_{m})}\right\\}$,
$\Omega_{\nu,m}=[0,\bar{t}^{\nu}_{\text{max}})\times\Omega$,

where $u^{n(\nu)}=\left(u_{1}^{n(\nu)},\ldots,u_{s}^{n(\nu)}\right)$ is the
unique mild solution of equation (11) corresponding to the initial state
$u_{0}^{\nu}$. In a similar way, we show that $-n(\nu)\leq
u^{n(\nu)}_{j}(t,x)\leq n(\nu)$, $j=1,\ldots,s$,
$(t,x)\in[0,t^{\nu}_{\text{max}})\times\Omega$. Then,
$F^{n(\nu)}(t,u^{n(\nu)}(t))=F(t,u^{\nu}(t))$, where

$u^{\nu}(t)=S(t)u^{\nu}_{0}+\displaystyle\int_{0}^{t}S(t-s)F(s,u^{\nu}(s))ds$
for  $t\in[0,t_{\text{max}})$,

is the mild solution of equation (7) corresponding to $u_{0}^{\nu}$. Since
$F(t,u^{\nu}(t))=F^{n(\nu)}(t,u^{n(\nu)}(t))$ is Lipschitz continuous in
$u^{n(\nu)}(t)$, uniformly with respect to $t\geq 0$, using growall Lemma, the
fact that $\|u_{0}-u_{0}^{\nu}\|\to 0$ as $\nu\to+\infty$, and $(S(t))_{t\geq
0}$ is a $C_{0}$-semigroups of contraction, we get that $(u^{\nu}(t))_{\nu\geq
0}$ is a Cauchy sequence in $(\mathbb{L}^{2}(\Omega))^{s}$ that converges to
$u(t)$ the unique mild solution of equation (7) corresponding to $u_{0}$
satisfying (9) on $[0,t_{\text{max}})$, while
$t_{\text{max}}=\lim\limits_{\nu\to+\infty}t^{\nu}_{\text{max}}\leq+\infty$.

From now and throughout the rest of this work, we consider the following
assumption. This assumption means that the system involves only monomolecular
reactions. In other words, the following must be verified: for each
$\kappa\in\\{1,\ldots,r\\}$, $0\leq\sum_{j=1}^{s}l_{j\kappa}\leq 1$ (see
[13]). Since $l_{j\kappa}\in\mathbb{N}$, this is equivalent to:

  1. (H)

For each $\kappa\in\\{1,\ldots,r\\}$, $l_{j\kappa}=0$ for all
$j\in\\{1,\ldots,s\\}$, or there exists a unique $j_{0}\in\\{1,\ldots,s\\}$
such that $l_{j_{0}\kappa}=1$ and $l_{j\kappa}=0$ for all
$j\in\\{1,\ldots,s\\}-\\{j_{0}\\}$.

In such cases, the equation (7) becomes linear in $u$; moreover,
$D(F)=(\mathbb{L}^{2}(\Omega))^{s}$. The following theorem shows the global
existence for equation (7).

###### Theorem 4.3

Assume that (H) holds. Then, equation (7) has a unique global solution defined
on $[0,+\infty[$.

The following Lemma will help us to demonstrate the result stated in Theorem
4.3.

###### Lemma 4.4

Assume that (H) holds, then there exists $\bar{d}>0$ and $\beta\in\\{0,1\\}$
such that $\|F(t,z)\|\leq\bar{d}\|z\|^{\beta}$ for each
$z\in(\mathbb{L}^{2}(\Omega))^{s}$.

###### Proof 3

Let $z=(z_{1},\ldots,z_{s})\in\left(\mathbb{L}^{2}(\Omega)\right)^{s}$. Let
$x\in\Omega$ and $t\in[0,t_{\text{max}})$. We discuss two cases:

i) case 1: if for each $\kappa\in\\{1,\ldots,r\\}$, $l_{j\kappa}=0$ for all
$j\in\\{1,\ldots,s\\}$, then for each $j=1,\ldots,s$, we have

$\displaystyle|F_{j}(t,z)(x)|^{2}$ $\displaystyle\leq$
$\displaystyle\left(\sum\limits_{\kappa=1}^{r}r_{j\kappa}h_{\kappa}(t)\right)^{2}\leq
2(r-1)\sum\limits_{\kappa=1}^{r}\left(r_{j\kappa}d_{\kappa}\right)^{2}.$

This implies that

$\displaystyle\|F_{j}(t,z)\|^{2}$ $\displaystyle\leq$ $\displaystyle
2(r-1)\sum\limits_{\kappa=1}^{r}\left(r_{j\kappa}d_{\kappa}\right)^{2}.$

Thus,

$\displaystyle\|F(t,z)\|$ $\displaystyle\leq$
$\displaystyle\left(2(r-1)\sum\limits_{j=1}^{s}\sum\limits_{\kappa=1}^{r}\left(r_{j\kappa}d_{\kappa}\right)^{2}\right)^{\frac{1}{2}}\|z\|^{\beta}\hskip
2.84544pt\text{ for }\hskip 2.84544pt\beta=0.$

ii) case 2: if for each $\kappa\in\\{1,\ldots,r\\}$, there exists a unique
$j_{0}\in\\{1,\ldots,s\\}$ such that $l_{j_{0}\kappa}=1$ and $l_{j\kappa}=0$
for all $j\in\\{1,\ldots,s\\}-\\{j_{0}\\}$, then for each $j=1,\ldots,s$, we
have

$\displaystyle|F_{j}(t,z)(x)|^{2}$ $\displaystyle\leq$
$\displaystyle\left(\sum\limits_{\kappa=1}^{r}(r_{j\kappa}-l_{j\kappa})h_{\kappa}(t)|z_{j_{0}}(x)|\right)^{2}$
$\displaystyle\leq$ $\displaystyle
2(r-1)\sum\limits_{\kappa=1}^{r}\left[(r_{j\kappa}-l_{j\kappa})h_{\kappa}(t)\right]^{2}|z_{j_{0}}(x)|^{2}.$

This implies that

$\displaystyle\|F_{j}(t,z)\|^{2}$ $\displaystyle\leq$ $\displaystyle
2(r-1)\sum\limits_{\kappa=1}^{r}\left[(r_{j\kappa}-l_{j\kappa})h_{\kappa}(t)\right]^{2}\|z\|^{2}$

Thus,

$\displaystyle\|F(t,z)\|\leq\left(2(r-1)\sum\limits_{j=1}^{s}\sum\limits_{\kappa=1}^{r}\left[\left(r_{j\kappa}-l_{j\kappa}\right)d_{\kappa}\right]^{2}\right)^{\frac{1}{2}}\|z\|^{\beta}\hskip
2.84544pt\text{ for }\hskip 2.84544pt\beta=1.$

In both cases, by denoting

$\bar{d}=\sqrt{2(r-1)}\max\left(\left(\sum\limits_{j=1}^{s}\sum\limits_{\kappa=1}^{r}\left(r_{j\kappa}d_{\kappa}\right)^{2}\right)^{\frac{1}{2}},\left(\sum\limits_{j=1}^{s}\sum\limits_{\kappa=1}^{r}\left[(r_{j\kappa}-l_{j\kappa})d_{\kappa}\right]^{2}\right)^{\frac{1}{2}}\right)$,

we get that for each $z\in(\mathbb{L}^{2}(\Omega))^{s}$,

$\displaystyle\|F(t,z)\|\leq\bar{d}\|z\|^{\beta}\hskip 2.84544pt\text{ for
}\hskip 2.84544pt\beta\in\\{0,1\\}.$

######  Proof 4 (Proof of Theorem 4.3)

According to Theorem 4.2, we know that equation (7) has a unique solution
$u(\cdot)$ defined on $[0,t_{\text{max}})$. To demonstrate that
$t_{\text{max}}=+\infty$, we will employ Lemma 4.4 and discuss two cases:

i) case 1: if $\beta=0$, we obtain that

$\displaystyle\|u(t)\|$ $\displaystyle\leq$
$\displaystyle\|S(t)u_{0}\|+\displaystyle\int_{0}^{t}\|S(t-\tau)F(\tau,u(\tau))\|d\tau$
$\displaystyle\leq$
$\displaystyle\|u_{0}\|+\displaystyle\int_{0}^{t}\|F(\tau,u(\tau))\|d\tau$
$\displaystyle\leq$ $\displaystyle\|u_{0}\|+\bar{d}t$

If $t_{\text{max}}<+\infty$, we get that

$\limsup\limits_{t\to
t_{\text{max}}}\|u(t)\|\leq\|u_{0}\|+\bar{d}t_{\text{max}}<+\infty$.

In light of Theorem 3.6, this leads to a contradiction, and therefore, we
conclude that $t_{\text{max}}=+\infty$.

ii) case 2: if $\beta=1$, then $F$ is at most affine w.r.t the second
argument, using Corollary 3.7, we obtain that $t_{\text{max}}=+\infty$.

The following Theorem shows the positivity of solutions for equation (7) on
$\left(\mathbb{L}^{2}(\Omega)\right)^{s}_{+}$.

###### Theorem 4.5

Assume that (H) holds, then the solutions of equation (7), starting from
$u_{0}\in\left(\mathbb{L}^{2}(\Omega)\right)^{s}_{+}$, remain in
$\left(\mathbb{L}^{2}(\Omega)\right)^{s}_{+}$.

###### Proof 5

We discuss two cases:

i) case 1: if for each $\kappa\in\\{1,\ldots,r\\}$, $l_{j\kappa}=0$ for all
$j\in\\{1,\ldots,s\\}$, then $F(t,z)\equiv F(t)$ doesn’t depend in $z$ and
$F(t)\geq 0$. By Theorem 4.3, we get that equation (7) has a unique solution
$u(\cdot)$ defined on $[0,+\infty[$ by

$u(t)=S(t)u_{0}+\displaystyle\int_{0}^{t}S(t-\tau)F(\tau)d\tau$ for  $t\geq
0$.

Since $(S(t))_{t\geq 0}$ is positive, we obtain that
$u(t)\in\left(\mathbb{L}^{2}(\Omega)\right)^{s}_{+}$.

ii) case 2: if for each $\kappa\in\\{1,\ldots,r\\}$, there exists a unique
$j_{0}\in\\{1,\ldots,s\\}$ such that $l_{j_{0}\kappa}=1$ and $l_{j\kappa}=0$
for all $j\in\\{1,\ldots,s\\}-\\{j_{0}\\}$, then $F(t,z)$ is linear in $z$.
Let $\lambda>0$ be sufficiently large such that $F(t,v)+\lambda v\geq 0$ for
each $v\in\left(\mathbb{L}^{2}(\Omega)\right)^{s}_{+}$. Equation (7) can be
written in the following equivalent form:

$\left\\{\begin{array}[]{l}u^{\prime}(t)=(\mathcal{A}-\lambda
I)u(t)+F(t,u(t))+\lambda u(t),\quad t\geq 0\\\ u(0)=u_{0}.\end{array}\right.$
(13)

The equation (13) has a unique solution $u(\cdot)$ defined on $[0,+\infty[$ by

$u(t)=e^{-\lambda
t}S(t)u_{0}+\displaystyle\int_{0}^{t}e^{-\lambda(t-\tau)}S(t-\tau)\left[F(\tau,u(\tau))+\lambda
u(\tau)\right]d\tau$ for  $t\geq 0$.

The result follows the fact that $(S(t))_{t\geq 0}$ is positive and
$F(t,v)+\lambda v\geq 0$ for each
$v\in\left(\mathbb{L}^{2}(\Omega)\right)^{s}_{+}$.

##  5 Global attractor

Understanding the long-term behaviour of dynamical systems is a fundamental
research challenge. A crucial concept in the study of the behaviour of such
systems is that of global attractors [11]. The fundamental results presented
in the previous section form the cornerstone of our study of the long-term
behaviour of the equation (7) using the theory of global attractors. In this
section, we carry out a complete examination of the long-term dynamics of the
equation (7). We thus establish the existence of a global attractor for our
model.

Let us first revisit some key properties that are derived from the theory of
global attractors and will help us in establishing our results in this
section. For a more comprehensive understanding, readers are encouraged to
refer to [11].

###### Definition 5.1

[11] A semiflow $(\mathcal{U}(t))_{t\geq 0}$ on a complete metric space $X$ is
a one parameter family maps $\mathcal{U}(t):X\to X$, parameter
$t\in\mathbf{R}^{+}$, such that $\mathcal{U}(0)=id_{X}$ ($id_{X}$ is the
identity map on $X$), $\mathcal{U}(t+s)=\mathcal{U}(t)\mathcal{U}(s)$ for
$t,s\geq 0$, and $t\to\mathcal{U}(t)x$ is continuous for all $x\in X$.

###### Definition 5.2

[11] Let $(\mathcal{U}(t))_{t\geq 0}$ be a semiflow on a complete metric space
$(X,d)$.

i)

We say that $(\mathcal{U}(t))_{t\geq 0}$ is bounded, if it takes bounded sets
into bounded sets.

ii)

We say that $B\subset X$ is invariant under $(\mathcal{U}(t))_{t\geq 0}$ if
$\mathcal{U}(t)B=B$ for $t\geq 0$.

iii)

We say that $(\mathcal{U}(t))_{t\geq 0}$ is point dissipative if there is a
bounded set $B\subset X$ which attracts each point of $X$ under
$(\mathcal{U}(t))_{t\geq 0}$.

iv)

A subset $B$ of $X$ is said to attract $C\subset X$ if
$d(\mathcal{U}(t)C,B)\to 0$ as $t\to+\infty$, where

$d(\mathcal{U}(t)C,B):=\inf\\{d(\mathcal{U}(t)x,y)\mid x\in B\text{ and }y\in
C\\}$.

###### Theorem 5.3

[11] Let $(\mathcal{U}(t))_{t\geq 0}$ be a semiflow on a complete space $X$
with the metric $\|\cdot\|_{X}$. Assume that there exists $L>0$ such that for
all $x_{0}\in X$, $\limsup\limits_{t\to+\infty}\|\mathcal{U}(t)x_{0}\|_{X}\leq
L$, then $(\mathcal{U}(t))_{t\geq 0}$ is point dissipative.

###### Definition 5.4

[11] Let $X$ be a complete metric space and $(\mathcal{U}(t))_{t\geq 0}$ be a
semiflow on $X$. Let $\mathcal{M}$ be a subset of $X$. Then, $\mathcal{M}$ is
called a global attractor for $(\mathcal{U}(t))_{t\geq 0}$ in $X$ if
$\mathcal{M}$ is a closed and bounded invariant (under
$(\mathcal{U}(t))_{t\geq 0}$) set of $X$ that attracts every bounded set of
$X$.

The following Theorem shows the existence of global attractors for point
dissipative semiflows.

###### Theorem 5.5

[11, Theorem 4.1.2, page 63] Let $(\mathcal{U}(t))_{t\geq 0}$ be a semiflow on
a complete metric space $X$. Suppose that $\mathcal{U}(t)$ is bounded for each
$t\geq 0$. If $(\mathcal{U}(t))_{t\geq 0}$ is point dissipative, then it has a
connected global attractor.

###  5.1 Why do we need to look for an attractor for equation (7)?

The set of (constant and non-constant) equilibrium points of equation (7) is
given by:

$\mathcal{S}=\left\\{u^{*}=(u_{1}^{*},\ldots,u_{s}^{*})\in
D(\mathcal{A})\mid\text{ such that }\mathcal{A}u^{*}+F(t,u^{*})=0\text{ for
}t\geq 0\right\\}$.

As we can see, these equilibrium points are not isolated. This lack of
isolation means that we cannot apply techniques such as linearization or
Lyapunov’s methods, typically used to study the asymptotic behavior of systems
near isolated equilibrium points. Furthermore, we do not anticipate observing
convergence towards these equilibrium points, which complicates the
understanding of the asymptotic behavior of our system as it approaches these
equilibrium points. Because of these challenges, our primary focus is on
establishing the existence of a finite fractal dimensional global attractor
for equation (7).

###  5.2  Existence of a global attractor for equation (7)

Let us define the family $(\mathcal{U}(t))_{t\geq 0}$ on
$(\mathbb{L}^{2}(\Omega))^{s}_{+}$ as follows: for $t\geq 0$ and
$u_{0}\in(\mathbb{L}^{2}(\Omega))^{s}_{+}$, we set
$\mathcal{U}(t)u_{0}:=u(t,u_{0}):=u(t)$, where $u(\cdot,u_{0})$ represents the
unique solution of equation (7) corresponding to
$u_{0}\in(\mathbb{L}^{2}(\Omega))^{s}_{+}$. It is clear that
$\mathcal{U}(0)u_{0}=u_{0}$, and for each
$u_{0}\in(\mathbb{L}^{2}(\Omega))^{s}_{+}$, the map $t\to\mathcal{U}(t)u_{0}$
is continuous. Additionally, due to the Lipschitz property of the function $F$
with respect to the second argument, we can demonstrate that
$\mathcal{U}(t+\tau)u_{0}=\mathcal{U}(t)\mathcal{U}(\tau)u_{0}$ for all
$t,\tau\geq 0$ and $u_{0}\in(\mathbb{L}^{2}(\Omega))^{s}_{+}$. Consequently,
we conclude that the family $(\mathcal{U}(t))_{t\geq 0}$ constitutes a
semiflow on $(\mathbb{L}^{2}(\Omega))^{s}_{+}$. The following Theorem shows
the boundnesse of the semiflow $(\mathcal{U}(t))_{t\geq 0}$ on
$(\mathbb{L}^{2}(\Omega))^{s}_{+}$.

###### Theorem 5.6

Assume that (H) holds, then $\mathcal{U}(t)$ is bounded for each $t\geq 0$.
Moreover,

$\|\mathcal{U}(t)u_{0}\|\leq e^{\bar{d}t}\left[\|u_{0}\|+1\right]$ for  $t\geq
0$,  and  $u_{0}\in(\mathbb{L}^{2}(\Omega))^{s}_{+}$.

###### Proof 6

We employ Lemma 4.4. Under assumption (H), we have two cases: $\beta=0$ or
$\beta=1$. We recall that

$\mathcal{U}(t)u_{0}=S(t)u_{0}+\displaystyle\int_{0}^{t}S(t-\tau)F(\tau,\mathcal{U}(\tau)u_{0})d\tau$
for  $t\geq 0$.

  1. 1)

Case 1: $\beta=0$. In this case, for $t\geq 0$ and
$u_{0}\in(\mathbb{L}^{2}(\Omega))^{s}_{+}$, by Lemma 4.4, we obtain that

$\displaystyle\|\mathcal{U}(t)u_{0}\|$ $\displaystyle\leq$
$\displaystyle\|u_{0}\|+\displaystyle\int_{0}^{t}\|F(\tau,\mathcal{U}(\tau)u_{0})\|d\tau$
$\displaystyle\leq$ $\displaystyle\|u_{0}\|+\bar{d}t.$

  2. 2)

Case 2: $\beta=1$. In this case, for $t\geq 0$ and
$u_{0}\in(\mathbb{L}^{2}(\Omega))^{s}_{+}$, by Lemma 4.4, we obtain that

$\displaystyle\|\mathcal{U}(t)u_{0}\|$ $\displaystyle\leq$
$\displaystyle\|u_{0}\|+\displaystyle\int_{0}^{t}\bar{d}\|\mathcal{U}(\tau)u_{0}\|d\tau.$

By Gronwall’s lemma, we conclude that

$\displaystyle\|\mathcal{U}(t)u_{0}\|$ $\displaystyle\leq$ $\displaystyle
e^{\bar{d}t}\|u_{0}\|.$

In both cases, we have

$\|\mathcal{U}(t)u_{0}\|\leq e^{\bar{d}t}\left[\|u_{0}\|+1\right]$ for  $t\geq
0$,  and  $u_{0}\in(\mathbb{L}^{2}(\Omega))^{s}_{+}$.

The following result shows the point dissipativeness of
$(\mathcal{U}(t))_{t\geq 0}$ on $(\mathbb{L}^{2}(\Omega))^{s}_{+}$.

###### Theorem 5.7

Assume that (H) holds, the semiflow $\mathcal{U}(t)_{t\geq 0}$ is point
dissipative on $(\mathbb{L}^{2}(\Omega))^{s}_{+}$.

###### Proof 7

As in the proof of Theorem 5.6, assumption (H) imply that $\beta=0$ or
$\beta=1$. Let $u_{0}\in(\mathbb{L}^{2}(\Omega))^{s}_{+}$ be fixed,
$\lambda>1$ be sufficiently large, and $t_{0}:=t_{0}(u_{0})>0$ such that

$\dfrac{e^{(\bar{d}+\lambda)t_{0}}\left[\|u_{0}\|+1\right]}{e^{\lambda t}}\leq
1\hskip 2.84544pt\text{ for }\hskip 2.84544ptt>t_{0}.$ (14)

We recall that equation (7) can take the equivalent form (13), and therefore:

$\mathcal{U}(t)u_{0}=e^{-\lambda
t}S(t)u_{0}\\!+\\!\displaystyle\int_{0}^{t}e^{-\lambda(t-\tau)}S(t-\tau)\left[F(\tau,\mathcal{U}(\tau)u_{0})\\!+\\!\lambda\mathcal{U}(\tau)u_{0}\right]d\tau$,
$t\geq 0$.

i) If $\beta=0$, for $t\geq t_{0}$, using Lemma 4.4, we have

$\displaystyle\|\mathcal{U}(t)u_{0}\|$ $\displaystyle\leq$ $\displaystyle
e^{-\lambda
t}\|u_{0}\|+\displaystyle\int_{0}^{t}e^{-\lambda(t-\tau)}\left[\bar{d}+\lambda\|\mathcal{U}(\tau)u_{0}\|\right]d\tau$
$\displaystyle=$ $\displaystyle e^{-\lambda
t}\|u_{0}\|+\bar{d}\displaystyle\int_{0}^{t}e^{-\lambda(t-\tau)}d\tau+\displaystyle\int_{0}^{t}\lambda
e^{-\lambda(t-\tau)}\|\mathcal{U}(\tau)u_{0}\|d\tau$ $\displaystyle=$
$\displaystyle e^{-\lambda
t}\|u_{0}\|+\frac{\bar{d}}{\lambda}\left[1-e^{-\lambda
t}\right]+\displaystyle\int_{0}^{t}\lambda
e^{-\lambda(t-\tau)}\|\mathcal{U}(\tau)u_{0}\|d\tau.$

By Theorem 5.6, we obtain that

$\displaystyle\|\mathcal{U}(t)u_{0}\|$ $\displaystyle\leq$ $\displaystyle
J_{1,\lambda}(t_{0},t)+\displaystyle\int_{t_{0}}^{t}\lambda
e^{-\lambda(t-\tau)}\|\mathcal{U}(\tau)u_{0}\|d\tau,$

where

$\displaystyle J_{1,\lambda}(t_{0},t)$ $\displaystyle=$ $\displaystyle
e^{-\lambda t}\|u_{0}\|+\frac{\bar{d}}{\lambda}\left[1-e^{-\lambda
t}\right]+\lambda e^{-\lambda
t}\displaystyle\int_{0}^{t_{0}}e^{(\bar{d}+\lambda)\tau}\left[\|u_{0}\|+1\right]d\tau$
$\displaystyle\leq$ $\displaystyle e^{-\lambda
t}\|u_{0}\|+\frac{\bar{d}}{\lambda}\left[1-e^{-\lambda
t}\right]+\frac{\lambda}{\bar{d}+\lambda}e^{-\lambda
t}\left[\|u_{0}\|+1\right](e^{(\bar{d}+\lambda)t_{0}}-1)$ $\displaystyle\leq$
$\displaystyle\dfrac{e^{(\bar{d}+\lambda)t_{0}}\left[\|u_{0}\|+1\right]}{e^{\lambda
t}}+\frac{\bar{d}}{\lambda}.$

By (14), we get that

$\displaystyle\|\mathcal{U}(t)u_{0}\|$ $\displaystyle\leq$
$\displaystyle\left(1+\frac{\bar{d}}{\lambda}\right)+\displaystyle\int_{t_{0}}^{t}w(t,\tau)\|\mathcal{U}(\tau)u_{0}\|d\tau,$

where

$w(t,\tau)=\lambda e^{-\lambda(t-\tau)}$ for  $t\geq\tau>t_{0}$.

By Lemma 3.10, we obtain that

$\displaystyle\|\mathcal{U}(t)u_{0}\|$ $\displaystyle\leq$
$\displaystyle\left(1+\frac{\bar{d}}{\lambda}\right)\exp\left(\displaystyle\int_{t_{0}}^{t}\left[\lambda-\lambda\displaystyle\int_{t_{0}}^{\tau}\lambda
e^{-\lambda(\tau-r)}dr\right]d\tau\right)$ $\displaystyle=$
$\displaystyle\left(1+\frac{\bar{d}}{\lambda}\right)\exp\left(\displaystyle\int_{t_{0}}^{t}\left[\lambda-\lambda(1-e^{-\lambda(\tau-
t_{0})})\right]d\tau\right)$ $\displaystyle=$
$\displaystyle\left(1+\frac{\bar{d}}{\lambda}\right)\exp\left(\displaystyle\int_{t_{0}}^{t}\lambda
e^{-\lambda(\tau-t_{0})}d\tau\right)$ $\displaystyle=$
$\displaystyle\left(1+\frac{\bar{d}}{\lambda}\right)\exp\left(1-e^{-\lambda(t-t_{0})}\right).$

ii) If $\beta=1$, for $t\geq t_{0}$, using Lemma 4.4, we obtain that

$\displaystyle\|\mathcal{U}(t)u_{0}\|$ $\displaystyle\leq$ $\displaystyle
e^{-\lambda
t}\|u_{0}\|+\displaystyle\int_{0}^{t}(\bar{d}+\lambda)e^{-\lambda(t-\tau)}\|\mathcal{U}(\tau)u_{0}\|d\tau.$

By Theorem 5.6, we get that

$\displaystyle\|\mathcal{U}(t)u_{0}\|$ $\displaystyle\leq$ $\displaystyle
J_{2,\lambda}(t_{0},t)+\displaystyle\int_{t_{0}}^{t}(\bar{d}+\lambda)e^{-\lambda(t-\tau)}\|\mathcal{U}(\tau)u_{0}\|d\tau,$

where

$\displaystyle J_{2,\lambda}(t_{0},t)$ $\displaystyle=$ $\displaystyle
e^{-\lambda t}\|u_{0}\|+e^{-\lambda
t}\displaystyle\int_{0}^{t_{0}}(\bar{d}+\lambda)e^{(\lambda+\bar{d})\tau}\left[\|u_{0}\|+1\right]d\tau$
$\displaystyle=$ $\displaystyle e^{-\lambda t}\|u_{0}\|+e^{-\lambda
t}\left[\|u_{0}\|+1\right]\left(e^{(\lambda+\bar{d})t_{0}}-1\right)$
$\displaystyle\leq$
$\displaystyle\dfrac{e^{(\bar{d}+\lambda)t_{0}}\left[\|u_{0}\|+1\right]}{e^{\lambda
t}}.$

By (14), we have

$\displaystyle\|\mathcal{U}(t)u_{0}\|$ $\displaystyle\leq$ $\displaystyle
1+\displaystyle\int_{0}^{t}w(t,\tau)\|\mathcal{U}(\tau)u_{0}\|d\tau,$

where

$w(t,\tau)=(\bar{d}+\lambda)e^{-\lambda(t-\tau)}$ for  $t\geq\tau>t_{0}$.

By Lemma 3.10, we obtain that

$\displaystyle\|\mathcal{U}(t)u_{0}\|$ $\displaystyle\leq$
$\displaystyle\exp\left(\displaystyle\int_{t_{0}}^{t}\left[(\bar{d}+\lambda)-(\bar{d}+\lambda)\displaystyle\int_{t_{0}}^{\tau}\lambda
e^{-\lambda(\tau-r)}dr\right]d\tau\right)$ $\displaystyle=$
$\displaystyle\exp\left(\displaystyle\int_{t_{0}}^{t}\left[(\bar{d}+\lambda)-(\bar{d}+\lambda)(1-e^{-\lambda(\tau-
t_{0})})\right]d\tau\right)$ $\displaystyle=$
$\displaystyle\exp\left(\dfrac{(\bar{d}+\lambda)}{\lambda}(1-e^{-\lambda(t-t_{0})})\right).$

In both cases, we obtain that

$\|\mathcal{U}(t)u_{0}\|\leq(1+\bar{d})e^{(1+\bar{d})}$ for  $t>t_{0}$.

As a consequence,

$\limsup\limits_{t\to+\infty}\|\mathcal{U}(t)u_{0}\|\leq(1+\bar{d})e^{(1+\bar{d})}$.

By Theorem 5.3, we conclude that $(\mathcal{U}(t))_{t\geq 0}$ is point
dissipative on $(\mathbb{L}^{2}(\Omega))^{s}_{+}$.

The main result in this section is the following, which directly follows from
Theorem 5.5, Theorem 5.6, and Theorem 5.7.

###### Theorem 5.8

Assume that (H) holds, then the semiflow $\mathcal{U}(t)_{t\geq 0}$ has a
unique connected global attractor $\mathcal{M}$ on
$(\mathbb{L}^{2}(\Omega))^{s}_{+}$.

###  5.3  Finite fractal dimensional

The main focus of the previous subsection was centered on the presence of the
global attractor $\mathcal{M}$ for the semiflow $(\mathcal{U}(t))_{t\geq 0}$
in the space $(\mathbb{L}^{2}(\Omega))^{s}$. In the realm of infinite-
dimensional dynamical systems, determining the fractal dimension of the global
attractor is of great importance. This significance arises from Mañé’s
Theorem, which states that if the global attractor has a finite fractal
dimension, it becomes possible to reduce the infinite-dimensional dynamical
system to a finite-dimensional counterpart existing within the attractor. Our
goal in this section is to establish the finite fractal dimension property of
the global attractor $\mathcal{M}$ of the semiflow $(\mathcal{U}(t))_{t\geq
0}$.

From [8], we recall the following definition.

###### Definition 5.9

The fractal dimension of $\mathcal{M}$ is defined by

$\dim_{f}(\mathcal{M})=\limsup\limits_{\varepsilon\to
0}ln(N_{\varepsilon})/ln(1/\varepsilon)$,

where $N_{\varepsilon}$ is the minor $n\in\mathbb{N}$ such that there exists
open balls $B_{1},\ldots,B_{n}$ with radius $\varepsilon$ such that
$\mathcal{M}\subset\bigcup\limits_{i=1}^{n}B_{i}$.

To show that $\mathcal{M}$ has a finite fractal dimension, we employ the
following theorem.

###### Theorem 5.10

[6, Theorem 2.1] Let $X$, $Y$, and $Z$ be Banach spaces such that $Y$ is
compactly embedded in $Z$ and let $\mathcal{B}$ be a bounded set of $X$,
invariant under a map $\mathcal{R}$ such that
$\mathcal{R}(\mathcal{B})=\mathcal{B}$. Assume that there exist a map
$\mathcal{T}:\mathcal{B}\to Y$, $\delta<\frac{1}{2}$ ($\delta\geq 0$),
$\mu>0$, and $k>0$ such that

$\|\mathcal{R}(u)-\mathcal{R}(v)\|_{X}\leq\delta\|u-v\|_{X}+\mu\|\mathcal{T}(u)-\mathcal{T}(v)\|_{Z},$

and

$\|\mathcal{T}(u)-\mathcal{T}(v)\|_{Y}\leq k\|u-v\|_{X},$

for all $u,v\in\mathcal{B}$. Under the above assumptions, $\mathcal{B}$ is a
precompact set in $X$ and it has a finite fractal dimension, more precisely,
for any $\nu\in(0,\frac{1}{2}-\delta)$,

$\dim_{f}(\mathcal{B})\leq\frac{\ln\left(N_{\frac{\nu}{k\mu}}^{Z}\left(B_{Y}(0,1)\right)\right)}{\ln\left(\frac{1}{2(\delta+\nu)}\right)},$

where $N_{\frac{\nu}{k\mu}}^{Z}\left(B_{Y}(0,1)\right)$ is the minimal number
of $\frac{\nu}{k\mu}$-balls in $Z$ need to cover $B_{Y}(0,1)$ the open unit
ball in $Y$.

The following Theorem is the main result in this subsection. It shows that the
attractor $\mathcal{M}$ has a finite fractal dimension.

###### Theorem 5.11

Assume that (H) holds, then the global attractor $\mathcal{M}$ of the semiflow
$(\mathcal{U}(t))_{t\geq 0}$ has a finite fractal dimension.

###### Proof 8

Let $u_{0},v_{0}\in\mathcal{M}$, and $\lambda:=\lambda(\mathcal{M})>0$ such
that $e^{-\lambda t}\leq\delta$ for $0<\delta<\frac{1}{2\sqrt{e}}$ and $t>0$.
Writing equation (7) in the equivalent form (13), implies that the semiflow
starting from $u_{0}$ and $v_{0}$ satisfies

$\displaystyle\|\mathcal{U}(t)u_{0}-\mathcal{U}(t)v_{0}\|$ $\displaystyle\leq$
$\displaystyle e^{-\lambda
t}\|u_{0}-v_{0}\|+\displaystyle\int_{0}^{t}e^{-\lambda(t-s)}\|F_{\lambda}(s,\mathcal{U}(s)u_{0})-F_{\lambda}(s,\mathcal{U}(s)v_{0})\|\,ds,$

where $F_{\lambda}(t,v)=F(t,v)+\lambda v$ for
$v\in(\mathbb{L}^{2}(\Omega))^{s}_{+}$. Since $u_{0}$ and $v_{0}$ are in
$\mathcal{M}$ which is invariant under $(\mathcal{U}(t))_{t\geq 0}$, it
follows that there exists $C:=C(\mathcal{M})>0$ such that

$\|F_{\lambda}(s,\mathcal{U}(s)u_{0})-F_{\lambda}(s,\mathcal{U}(s)v_{0})\|\leq(\lambda+C)\|\mathcal{U}(s)u_{0}-\mathcal{U}(s)v_{0}\|.$

Therefore,

$\displaystyle\|\mathcal{U}(t)u_{0}-\mathcal{U}(t)v_{0}\|$ $\displaystyle\leq$
$\displaystyle\delta\|u_{0}-v_{0}\|+\displaystyle\int_{0}^{t}w(t,s)\|\mathcal{U}(s)u_{0}-\mathcal{U}(s)v_{0}\|ds,$

where $w(t,s)=e^{-\lambda(t-s)}(\lambda+C)$ for $t\geq s\geq 0$. By Lemma
3.10, we get the following

$\displaystyle\|\mathcal{U}(t)u_{0}-\mathcal{U}(t)v_{0}\|$ $\displaystyle\leq$
$\displaystyle\delta\exp\left(\displaystyle\int_{0}^{t}\left[\lambda+C-(\lambda+C)e^{-\lambda\tau}\displaystyle\int_{0}^{\tau}\lambda
e^{\lambda r}dr\right]d\tau\right)\|u_{0}-v_{0}\|$ $\displaystyle=$
$\displaystyle\delta\exp\left(\displaystyle\int_{0}^{t}\left[\lambda+C-(\lambda+C)e^{-\lambda\tau}\left(e^{\lambda\tau}-1\right)\right]d\tau\right)\|u_{0}-v_{0}\|$
$\displaystyle=$
$\displaystyle\delta\exp\left((\lambda+C)\displaystyle\int_{0}^{t}e^{-\lambda\tau}d\tau\right)\|u_{0}-v_{0}\|$
$\displaystyle=$
$\displaystyle\delta\exp\left(\dfrac{(\lambda+C)}{\lambda}\left[1-e^{-\lambda
t}\right]\right)\|u_{0}-v_{0}\|.$

Let $t^{*}>0$ such that $(1-e^{-\lambda
t^{*}})\leq\frac{\lambda}{2(\lambda+C)}$. Then,

$\|\mathcal{U}(t^{*})u_{0}-\mathcal{U}(t^{*})v_{0}\|\leq\delta\sqrt{e}\|u_{0}-v_{0}\|\hskip
2.84544pt\text{ for all }\hskip 2.84544ptu_{0},v_{0}\in\mathcal{M}.$

Since $0<\delta\sqrt{e}<\frac{1}{2}$, by applying Theorem 5.10 for
$X=Z=(\mathbb{L}^{2}(\Omega))^{s}$, $Y=(H^{1}(\Omega))^{s}$,
$\mathcal{R}=\mathcal{U}(t^{*})$, $\mathcal{T}=0$, $\mathcal{B}=\mathcal{M}$,
and $\mu=k=1$, we obtain that $\mathcal{M}$ has a finite fractal dimension,
moreover, for any $\nu\in(0,\frac{1}{2}-\delta\sqrt{e})$,

$\dim_{f}(\mathcal{M})\leq\frac{\ln\left(N_{\nu}\left(B_{H^{1}}(0,1)\right)\right)}{\ln\left(\frac{1}{2(\delta\sqrt{e}+\nu)}\right)},$

where $N_{\nu}\left(B_{H^{1}}(0,1)\right)$ is the minimal number of
${\nu}$-balls in $(\mathbb{L}^{2}(\Omega))^{s}$ needed to cover
$B_{H^{1}}(0,1)$ the open unit ball in $(H^{1}(\Omega))^{s}$.

##  6 Analytical solution of an Advection-Diffusion Equation in 2-D using
separation of variables

Finding an analytical solution to the Advection-Diffusion-Reaction (ADR)
equation, with constant advection and diffusion coefficients, Dirichlet
boundary conditions, and an initial condition, which considers the reaction in
its general form, is often infeasible. This is due to the non-linearity added
by the reaction term $R(c)$, the sum of reactions that are products of
concentrations. The analytical solution exists for certain forms of $R(c)$
that can be simplified. In the literature, the Laplace Transform method is
used to solve the ADR equation, with constant advection and diffusion
coefficients, Dirichlet boundary conditions, and an initial condition, with a
constant first-order reaction and a source function [14]. Green’s function is
used when partial differential equations (PDEs) are linear, and the integral
involving the reaction term in the Green function approach converges [20]. The
analytical solution of the Advection-Diffusion equation, with constant
advection and diffusion coefficients, Dirichlet boundary conditions, and an
initial condition, using the separation of variables approach while
considering there is no reaction between species, is revisited in this
section. This analytical solution is one of a subset of ADR equations, with
constant advection and diffusion coefficients, Dirichlet boundary conditions,
and an initial condition.

###  6.1 Analytical solution in 2-D

As presented in [4], the solution of the two-dimension advection-diffusion
equation with constant advection and diffusion coefficient terms and the
following initial and boundary conditions is presented:

$\left\\{\begin{array}[]{l}\dfrac{\partial c}{\partial
t}=k\dfrac{\partial^{2}c}{\partial x^{2}}+k\dfrac{\partial^{2}c}{\partial
y^{2}}-u\dfrac{\partial c}{\partial x}-u\dfrac{\partial c}{\partial
y},\quad\Omega=[0,1]\times[0,1],\\\\[10.0pt]
c_{\mid\partial\Omega}=0,\\\\[10.0pt] c(0,x,y,)=f(x,y).\end{array}\right.$
(15)

The solution is of the separable form

$c(t,x,y)=X(x)Y(y)T(t),$ (16)

substituting into equation (15) and dropping $x,y,t$ terms

$XYT^{\prime}=\left(kX^{\prime\prime}-uX^{\prime}\right)YT+X\left(kY^{\prime\prime}-uY\right)T.$

Rearranging leads to

$\frac{T^{\prime}}{T}=\frac{kX^{\prime\prime}-uX^{\prime}}{X}+\frac{kY^{\prime\prime}-uY^{\prime}}{Y}.$

The variables $t$, $x$ and $y$ are supposed independent, then each function
variable in the above equation is a constant. Meaning this equation is only
possible if each fraction is equal to a number $\lambda$, $\alpha$, and
$\beta$ where $\lambda=\alpha+\beta$. From symmetry, the solutions for the
$x$-dependent differential equation and $y$-dependent differential equation
will be the same with different constants (eigenvalues). The eigenproblem gets
reduced to two first order differential equations with

$\displaystyle T^{\prime}(t)=\lambda T(t),\quad T(0)=1,$ (17) $\displaystyle
kX^{\prime\prime}(x)-uX^{\prime}(x)=\alpha X(x).$ (18)

Solving equation (17), we get

$T(t)=e^{\lambda t}.$ (19)

Below the solution of equation (18) is worked out, and the characteristic
equation for equation (18) is presented:

$\displaystyle kr^{2}-ur-\alpha=0;\hskip 2.84544pt\mu=u^{2}+4k\alpha;\hskip
2.84544ptr_{1}=\frac{u+\sqrt{u^{2}+4k\alpha}}{2k};\hskip
2.84544ptr_{2}=\frac{u-\sqrt{u^{2}+4k\alpha}}{2k}.$ (20)

The solution will depend on the discriminant (20) of the characteristic
equation. The three possible solutions are:

$\displaystyle
X(x)=\left\\{\begin{array}[]{l}Ee^{r_{1}x}+Fe^{r_{2}x},\quad\mu>0\\\
Ce^{r_{1}x}+D,\quad\mu=0\\\
e^{\frac{u}{2d}x}\left[A\sin\left(\frac{\sqrt{-u^{2}-4k\alpha}}{2k}x\right)+B\cos\left(\frac{\sqrt{-u^{2}-4k\alpha}}{2k}x\right)\right],\quad\mu<0.\end{array}\right.$

However the only equation, from which we can construct a non-trivial solution,
is when $\mu<0$ and thus $r$ is imaginary:

$X(x)=e^{\frac{u}{2k}x}\left[A\sin\left(\frac{\sqrt{-u^{2}-4k\alpha}}{2k}x\right)+B\cos\left(\frac{\sqrt{-u^{2}-4k\alpha}}{2k}x\right)\right].$

The next steps determine the solution that satisfies the Dirichlet boundary
conditions then the initial condition, we have:

$\displaystyle c(t,0,y)=0\implies X(0)Y(y)T(t)=0,$ $\displaystyle
c(t,1,y)=0\implies X(1)Y(y)T(t)=0.$ Hence, $\displaystyle X(0)=0\text{ and
}X(1)=0,$ $\displaystyle X(0)=0\implies B=0\hskip 2.84544pt\text{ and }\hskip
2.84544ptX(1)=0\implies\sin\left(\frac{\sqrt{-u^{2}-4k\alpha}}{2k}\right)=0.$
Thus, $\displaystyle\frac{\sqrt{-u^{2}-4k\alpha}}{2k}=m\pi\text{ for
}m=1,2,3\ldots.$

The solution is reduced to:

$\displaystyle X_{m}(x)=A_{m}e^{\frac{u}{2k}x}\sin(m\pi x)\text{ for
}m=1,2,3\ldots,\hskip 2.84544pt\text{ with }\hskip
2.84544pt\alpha_{m}=-km^{2}\pi^{2}-\frac{u^{2}}{4k}.$

Based on the superposition principle, any linear combination $X(x)$ is a
solution of a linear homogeneous problem. Also, the eigenfunctions $X_{m}(x)$
and their eigenvalues will have the same form for $Y(y)$. The product of the
solutions for the $\mathrm{x}$-dependent differential equation and
$y$-dependent differential equation will be of the form of a Fourier Sine
series satisfying Dirichlet’s boundary conditions:

$X(x)Y(y)=\sum_{m=1}^{\infty}\sum_{n=1}^{\infty}A_{m,n}e^{\frac{u}{2k}x}e^{\frac{u}{2k}y}\sin(m\pi
x)\sin(n\pi y).$

Now, remains finding the Fourier coefficient $A_{m,n}$ such that this
expression holds true under the initial condition. Using the orthogonality
properties of the eigenfunctions and by combining equations (16), (19) and the
initial condition, we get:

$c(0,x,y)=f(x,y)=X(x)Y(y)T(0)=X(x)Y(y)$ $\implies
f(x,y)=X(x)Y(y)=\sum_{n=1}^{\infty}\sum_{n=1}^{\infty}A_{m,n}e^{\frac{u}{2k}x}e^{\frac{u}{2k}y}\sin(m\pi
x)\sin(n\pi y).$

Based on Fubini’s theorem where the condition of absolute convergence holds,
integrating, we get that

$\displaystyle\int_{0}^{1}\int_{0}^{1}f(x,y)dxdy=\int_{0}^{1}\int_{0}^{1}\sum_{m=1}^{\infty}\sum_{n=1}^{\infty}A_{m,n}e^{\frac{u}{2k}x}e^{\frac{u}{2k}y}\sin(m\pi
x)\sin(n\pi y)dxdy,$

which implies that

$\displaystyle\int_{0}^{1}\int_{0}^{1}f(x,y)e^{-\frac{u}{2k}x}e^{-\frac{u}{2k}y}\sin(l\pi
x)\sin(k\pi y)dxdy=$
$\displaystyle\sum_{m=1}^{\infty}\sum_{n=1}^{\infty}\int_{0}^{1}\int_{0}^{1}A_{m,n}\sin(m\pi
x)\sin(n\pi y)\sin(l\pi x)\sin(k\pi y)dxdy,$ (21)
$\displaystyle\left\\{\begin{array}[]{l}\text{when }k,l\neq
m,n\implies\eqref{eq_multiply}=0\\\ \text{when
}k,l=m,n\implies\eqref{eq_multiply}=\frac{1}{4}A_{k,l}.\end{array}\right.$

The previous implication occurs when $k,l=m,n$ and stems from:

$\displaystyle\sum_{m=1}^{\infty}\sum_{n=1}^{\infty}\int_{0}^{1}\int_{0}^{1}A_{m,n}\sin(m\pi
x)\sin(n\pi y)\sin(l\pi x)\sin(k\pi y)dxdy$
$\displaystyle=\frac{1}{4}A_{k,l}\int_{0}^{1}\int_{0}^{1}2\sin^{2}\left(k\pi
x\right)2\sin^{2}(l\pi y)dxdy$
$\displaystyle=\frac{1}{4}A_{k,l}\int_{0}^{1}(1-\cos(2k\pi
x))dx\int_{0}^{1}(1-\cos(2k\pi y))dy$
$\displaystyle=\frac{1}{4}A_{k,l}\left[\left(x-\left.\frac{\sin(2k\pi
x)}{2k\pi}\right|_{0}^{1}\right)_{0}\left(y-\left.\frac{\sin(2l\pi
y)}{2k\pi}\right|_{0}^{1}\right)\right]$ $\displaystyle=\frac{1}{4}A_{k,l}.$

Therefore (21) reduces to:

$A_{m,n}=4\int_{0}^{1}\int_{0}^{1}f(x,y)e^{-\frac{u}{2k}x}e^{-\frac{u}{2k}y}\sin(m\pi
x)\sin(n\pi y)dxdy$ (22)

for $m,n=1,2,3\ldots$. Now back to the time-dependent solution:

$\displaystyle T(t)=e^{\lambda
t},\quad\lambda_{m,n}=\alpha_{m}+\beta_{n}=-km^{2}\pi^{2}-\frac{u^{2}}{4k}-kn^{2}\pi^{2}-\frac{u^{2}}{4k}.$

Combining the time-dependent solution with the spatial solution we get the
final solution to the advection-diffusion equation:

$c(t,x,y)=\sum_{m=1}^{\infty}\sum_{n=1}^{\infty}A_{m,n}e^{\left(-km^{2}\pi^{2}-kn^{2}\pi^{2}-\frac{u^{2}}{2k}\right)t}e^{\frac{u}{2k}(x+y)}\sin(m\pi
x)\sin(n\pi y),$

where $A_{m,n}$ is calculated by (22).

###  6.2 Analytical Example

We will specify an advection diffusion equation and initial conditions:

$\left\\{\begin{array}[]{l}\dfrac{\partial c}{\partial
t}=.5\dfrac{\partial^{2}c}{\partial x^{2}}+.5\dfrac{\partial^{2}c}{\partial
y^{2}}-5\dfrac{\partial c}{\partial x}-5\dfrac{\partial c}{\partial
y},\\\\[10.0pt] c(0,x,y)=f(x,y)=\sin(\pi x)\sin(\pi y).\end{array}\right.$

Applying these new advection-diffusion coefficients and initial conditions,
the solutions are:

$\displaystyle
c(t,x,y)=\sum_{m=1}^{\infty}\sum_{n=1}^{\infty}A_{m,n}e^{\left(-.5m^{2}\pi^{2}-.5n^{2}\pi^{2}-25\right)t}e^{5(x+y)}\sin(m\pi
x)\sin(n\pi y),$

where

$A_{m,n}=4\displaystyle\int_{0}^{1}\displaystyle\int_{0}^{1}e^{-5(x+y)}\sin(\pi
x)\sin(\pi y)\sin(m\pi x)\sin(n\pi y)dxdy$, $m,n\geq 1$.

Graphed in Figure 1 are the solutions found using initial conditions and
coefficients. The coefficients of the Fourier Sine series can be calculated by
computer, in this example taking the first 40 terms. We graphed these
solutions using Python 3.10.5 at different times stamps. The analytical
solutions diffuse over time and move away from the origin. In the next
section, we will discuss a numerical simulation of this same example, and
compare the results (convergence and graphs) to these analytical solutions.

Figure 1: Analytical solutions at times 0, 0.05, 0.09, 0.12, with $u=5$ and
$k=0.5$.

##  7 Numerical Simulations with Finite Difference Methods

Numerical methods, such as finite difference, can provide accurate and
efficient approximations to the analytic solution, allowing for the study of
complex systems and realistic scenarios that do not admit simple analytical
solutions. In this section, we will solve numerically the example presented
above, and present a new example including the reaction between Ozone,
$\mathrm{NO}$ and $\mathrm{NO}_{2}$ in clean atmosphere.

###  7.1 Numerical simulations in 2-D using explicit centered difference
method

Considering the same advection-diffusion equation stated in the analytic
example above with the same initial conditions, discretization with the
explicit centered difference method is shown. The time derivative is
approximated using a Forward Difference approximation, the advection and the
diffusion terms are approximated using a Central Difference approximation,
where $i$ and $j$ represent the spatial index and $n$ represents the temporal
index. $c_{i,j}^{n}$ is the solution at grid point $\left(x_{i},y_{j}\right)$
at time $n$:

$\displaystyle\frac{c_{i,j}^{n+1}-c_{i,j}^{n}}{\delta t}=$
$\displaystyle-u\frac{c_{i+1,j}^{n}-c_{i-1,j}^{n}}{2\delta
x}-u\frac{c_{i+1,j}^{n}-c_{i,j-1}^{n}}{2\delta y}$
$\displaystyle+k\frac{c_{i+1,j}^{n}-2c_{i,j}^{n}+c_{i-1,j}^{n}}{\delta
x^{2}}+k\frac{c_{i,j+1}^{n}-2c_{i,j}^{n}+c_{i,j-1}^{n}}{\delta y^{2}}.$

####  7.1.1 Stability

Let’s consider a uniform grid (and advection-diffusion coefficients), and
let’s introduce the Diffusion number and the Cell Peclet number:

$\begin{gathered}R_{x}=R_{y}=k\frac{\delta t}{\delta
x^{2}},P_{x}=\frac{u\delta x}{2k},\text{ and }P_{y}=\frac{u\delta
y}{2k}.\end{gathered}$

The finite difference scheme rearranged gives:

$\displaystyle c_{i,j}^{n+1}=c_{i,j}^{n}+$ $\displaystyle
R_{x}\left(c_{i+1,j}^{n}-2c_{i,j}^{n}+c_{i-1,j}^{n}\right)+R_{y}\left(c_{i,j+1}^{n}-2c_{i,j}^{n}+c_{i,j-1}^{n}\right)$
$\displaystyle-
P_{x}R_{x}\left(c_{i+1,j}^{n}-c_{i-1,j}^{n}\right)-P_{y}R_{y}\left(c_{i,j+1}^{n}-c_{i,j-1}^{n}\right).$

Grouping like terms,

$\displaystyle c_{i,j}^{n+1}=$
$\displaystyle\left(1-2R_{x}-2R_{y}\right)c_{i,j}^{n}+\left(R_{x}-P_{x}R_{x}\right)c_{i+1,j}^{n}+\left(R_{x}+P_{x}R_{x}\right)c_{i-1,j}^{n}$
(23)
$\displaystyle+\left(R_{y}-P_{y}R_{y}\right)c_{i,j+1}^{n}+\left(R_{y}+P_{y}R_{y}\right)c_{i,j-1}^{j}.$

The solution will be stable if the coefficients are all positive. The
following constraints on $R_{x}$ and $P_{x}$ rise:

$1-2R_{x}-2R_{y}=1-4R_{x}>0\quad or\quad R_{x}<\frac{1}{4}\quad and\quad
R_{x}-P_{x}R_{x}>0\quad or\quad P_{x}<1.$

####  7.1.2 Convergence

Next, the convergence of the explicit centered method is analyzed using a
Taylor expansion for each term:

$\displaystyle c_{i,j}^{n+1}=c_{i,j}^{n}+\left(\frac{\partial c}{\partial
t}\right)_{i,j}^{n}\delta t+\left(\frac{\partial^{2}c}{\partial
t^{2}}\right)_{i,j}^{n}\frac{\delta\mathrm{t}^{2}}{2!}+\ldots,$ $\displaystyle
c_{i+1,j}^{j}=c_{i,j}^{n}+\left(\frac{\partial c}{\partial
x}\right)_{i,j}^{n}\delta x+\left(\frac{\partial^{2}c}{\partial
x^{2}}\right)_{i,j}^{n}\frac{\delta
x^{2}}{2!}+\left(\frac{\partial^{3}c}{\partial
x^{3}}\right)_{i,j}^{n}\frac{\delta x^{3}}{3!}\ldots,$ $\displaystyle
c_{i-1,j}^{n}=c_{i,j}^{n}-\left(\frac{\partial c}{\partial
x}\right)_{i,j}^{n}\delta x+\left(\frac{\partial^{2}c}{\partial
x^{2}}\right)_{i,j}^{n}\frac{\delta
x^{2}}{2!}-\left(\frac{\partial^{3}c}{\partial
x^{3}}\right)_{i,j}^{n}\frac{\delta x^{3}}{3!}\ldots,$ $\displaystyle
c_{i,j+1}^{n}=c_{i,j}^{n}+\left(\frac{\partial c}{\partial
y}\right)_{i,j}^{n}\delta y+\left(\frac{\partial^{2}c}{\partial
y^{2}}\right)_{i,j}^{n}\frac{\delta
y^{2}}{2!}+\left(\frac{\partial^{3}c}{\partial
y^{3}}\right)_{i,j}^{n}\frac{\delta y^{3}}{3!}\ldots,$ $\displaystyle
c_{i,j-1}^{n}=c_{i,j}^{n}-\left(\frac{\partial c}{\partial
y}\right)_{i,j}^{n}\delta y+\left(\frac{\partial^{2}c}{\partial
y^{2}}\right)_{i,j}^{n}\frac{\delta
y^{2}}{2!}-\left(\frac{\partial^{3}c}{\partial
y^{3}}\right)_{i,j}^{n}\frac{\delta y^{3}}{3!}\ldots.$

Plugging into (23), simplifying and dropping the subscripts

$\begin{gathered}\frac{\partial c}{\partial t}=k\frac{\partial^{2}c}{\partial
x^{2}}+\mathrm{k}\frac{\partial^{2}c}{\partial y^{2}}-u\frac{\partial
c}{\partial t}-u\frac{\partial c}{\partial x}-\frac{\delta
t}{2!}\left(\frac{\partial^{2}c}{\partial t^{2}}\right)-\frac{u\delta
x^{2}}{3!}\frac{\partial^{3}c}{\partial x^{3}}-\frac{u\delta
y^{2}}{3!}\frac{\partial^{3}c}{\partial y^{3}}.\end{gathered}$

The leading error terms are

$E=-\frac{\delta t}{2!}\left(\frac{\partial^{2}c}{\partial
t^{2}}\right)-\frac{u\delta x^{2}}{3!}\frac{\partial^{3}c}{\partial
x^{3}}-\frac{u\delta y^{2}}{3!}\frac{\partial^{3}c}{\partial y^{3}}.$

So the error is second order in space and first order in time Error
$\sim\delta t,\delta x^{2},\delta y^{2}$. In addition, as $\delta t,\delta
x^{2},\delta y^{2}\rightarrow 0$ the numerical solution converges to the
analytical solution showing consistency.

####  7.1.3 Numerical Example

Using the explicit centered method, relatively small step-sizes in the $x,y$,
and $t$-direction, are chosen, with nx and $ny=46$ the Number of grid points
in $x$\- and $y$-direction. The final resolution is $\delta x=\delta
y=\frac{L}{nx-1}=\frac{L}{ny-1}=0.02222,dt=0.0001,$ in order to verify the
stability conditions, $R_{x}=0.10125<\frac{1}{4}$ and $P_{x}=0.1112<1$. In
Figure 2, the numerical solutions diffuse over time and are moving away from
the origin. As shown in the previous section the numerical solution converges
to the analytical solution, as seen from the similar graphs, which also show
properties of diffusion and dispersion.

####  7.1.4 Error Analysis

Figure 3 shows that the explicit centered numerical scheme performs well at
approximating the analytical solution, considering coefficients that verify
stability conditions of the numerical scheme.

###  7.2 Numerical simulations in 3-D using the Explicit Upwind Scheme

Ozone is the product of two reactions between 4 species in a clean atmosphere,
but in reality this is not the case in the event of pollution, where Ozone is
the result of a chemical reaction between $\mathrm{NO}_{2}$ and $\mathrm{VOC}$
at 99%. The example presented below does not aim to model realistically the
reactions, but aims to demonstrate the mathematical and physical principles by
carrying out simulations in a clean atmosphere. Considering the following
three-dimensional Advection-Diffusion-Reaction equation:

$\frac{\partial c}{\partial t}+u_{x}\frac{\partial c}{\partial
x}+u_{y}\frac{\partial c}{\partial y}+u_{z}\frac{\partial c}{\partial
z}=k_{x}\frac{\partial^{2}c}{\partial
x^{2}}+k_{y}\frac{\partial^{2}c}{\partial
y^{2}}+k_{z}\frac{\partial^{2}c}{\partial z^{2}}+R(c).$

Figure 2: Numerical solutions at times 0, 0.05, 0.09, 0.12, with $u=5$ and
$k=0.5$.

Discretization with the explicit First Order Upwind Scheme is used. The time
derivative is approximated using a Forward Difference approximation, the
advection term is approximated using a One-Sided Difference approximation (for
$u>0$) and the diffusion term is approximated using a Central Difference
approximation, where $i$, $j$ and $k$ represent the spatial index and $n$
represents the temporal index. $c_{i,j,k}^{n}$ is the concentration of each
pollutant at grid point $\left(x_{i},y_{j},z_{k}\right)$ at time $n$:

$\displaystyle c_{i,j,k}^{n+1}$ $\displaystyle\\!\\!\\!=\\!\\!\\!$
$\displaystyle c_{i,j,k}^{n}+\delta
t\left[-u_{x}\frac{c_{i,j,k}^{n}-c_{i-1,j,k}^{n}}{\delta
x}-u_{y}\frac{c_{i,j,k}^{n}-c_{i,j-1,k}^{n}}{\delta
y}-u_{x}\frac{c_{i,j,k}^{n}-c_{i,j,k-1}^{n}}{\delta z}\right.$
$\displaystyle\left.+k_{x}\frac{c_{i+1,j,k}^{n}-2c_{i,j,k}^{n}+c_{i-1,j,k}^{n}}{\delta
x^{2}}+k_{y}\frac{c_{i,j+1,k}^{n}-2c_{i,j,k}^{n}+c_{i,j-1,k}^{n}}{\delta
y^{2}}\right.$
$\displaystyle\left.+k_{z}\frac{c_{i,j,k+1}^{n}-2c_{i,j,k}^{n}+c_{i,j,k-1}^{n}}{\delta
z^{2}}+R(c_{i,j,k}^{n})\right].$

Considering 2 reactions between $\mathrm{NO},\mathrm{NO}_{2},\mathrm{O}_{3}$
and $\mathrm{O}_{2}$:

$\displaystyle\mathrm{NO}_{2}+\mathrm{O}_{2}+hv\stackrel{{\scriptstyle\mathrm{k}_{1}}}{{\longrightarrow}}\mathrm{NO}+\mathrm{O}_{3},$
$\displaystyle\mathrm{NO}+\mathrm{O}_{3}\stackrel{{\scriptstyle\mathrm{k}_{2}}}{{\longrightarrow}}\mathrm{O}_{2}+\mathrm{NO}_{2}.$

These reactions are basic to tropospheric air pollution models. The first
reaction is photochemical and states that $\mathrm{NO}$ and $\mathrm{O}_{3}$
are formed from the photo-dissociation $hv$ caused by solar radiation of
$\mathrm{NO}_{2}$ and $\mathrm{O}_{2}$. This depends on the time of the day
and therefore $k_{l}=k_{l}(t)$. Let the concentrations of the pollutants be
$c_{1}=[NO],c_{2}=\left[\mathrm{NO}_{2}\right]$ and
$c_{3}=\left[\mathrm{O}_{3}\right]$. $\left[O_{2}\right]$ is constant, and
consider a constant source term $\sigma_{2}$ at the cell $(1,1,1)$ simulating
emissions of NO. The speed of reactions are:

$\displaystyle R_{1}(t,c_{1}(t))=f\left(t,c_{1}(t)\right)=\frac{\partial
c_{1}(t)}{\partial t}=k_{1}(t)c_{2}(t)-k_{2}c_{1}(t)c_{3}(t)+\sigma_{2},$
$\displaystyle R_{2}(t,c_{2}(t))=f\left(t,c_{2}(t)\right)=\frac{\partial
c_{2}(t)}{\partial t}=k_{2}c_{1}(t)c_{3}(t)-k_{1}(t)c_{2}(t),$ $\displaystyle
R_{3}(t,c_{3}(t))=f\left(t,c_{3}(t)\right)=\frac{\partial c_{3}(t)}{\partial
t}=k_{1}(t)c_{2}(t)-k_{2}c_{1}(t)c_{3}(t).$

The set of ODEs can be written in system form:

$R(t,c(t))=u^{\prime}(t)=Sg(t,u(t)),\quad c(0)\text{ is given },$

where $S=\left(r_{ij}-l_{ij}\right)$ is an $s\times r$ matrix, of the loss and
gain of the number of molecules $U_{i}$ in the $jth$ reaction, where $s$ is
the number of species, $r$ is the number of reactions,
$g(t,u)=\left(g_{j}(t,u)\right)\in\mathbb{R}^{r}$,

$g(t,u)=\left(\begin{array}[]{c}k_{1}(t)c_{2}(t)\\\ k_{2}c_{1}(t)c_{3}(t)\\\
\end{array}\right),\quad S=\left(\begin{array}[]{rrr}1&-1\\\ -1&1\\\
1&-1\end{array}\right).$

Figure 3: Maximum error between Analytical and Numerical solutions at times 0,
0.05, 0.09, 0.12, with $u=5$, $k=0.5$, and $m=n=40$.

Let $c(0,(1,1,1))=\left(1.3\times 10^{8},5.0\times 10^{11},8.0\times
10^{11}\right)^{T},$ be the initial values of the concentrations of the three
respective pollutants in molecules per $cm^{3}$, $\sigma_{2}=10^{6}$ at the
initial cell $(1,1,1)$ at each time-step $\delta t=1s$ and the reaction
coefficients:

$k_{2}=10^{-16},k_{1}(t)=\begin{cases}10^{-5}e^{7\sec(t)}&\text{ in daytime
}\\\ 10^{-40}&\text{ during nighttime, }\end{cases}$

where

$\begin{gathered}\sec(t)=\left(\sin\left(\frac{1}{16}\pi\left(\bar{t}_{h}-4\right)\right)\right)^{0.2},\bar{t}_{h}=t_{h}-24\lfloor
t_{h}/24\rfloor,t_{h}=t/3600,\end{gathered}$

with daytime between $4$ am and $8$ pm and $\lfloor t_{h}/24\rfloor$ the
largest integer $\leqslant t_{h}/24$. $\sec(t)$ is periodic with a period of
$24\mathrm{~{}h}$ but defined only during daytime. The maximum for
$k_{1}\approx 0.01$ occurs at noon. The solution will be stable with the
following constraint:

$R_{x}=\frac{k\delta t}{\delta x^{2}}\text{ and }P_{x}=\frac{u\delta
x}{k},\quad 2R_{x}+2R_{y}+2R_{z}+P_{x}R_{x}+P_{y}R_{y}+P_{z}R_{z}<1.$

The computational parameters for our simulations are as follows:

$\displaystyle\delta t=1s,$ $\displaystyle\delta x=\delta y=\delta z=10m,$
$\displaystyle\text{Initial conditions}=c((0,(1,1,1))=\left(1.3\times
10^{8},5.0\times 10^{11},8.0\times 10^{11}\right)^{T},$
$\displaystyle\text{Boundary conditions: }c=(0,0,0)^{T}\text{at $x$, $y$ or
$z=0m$ and at $x,y$ or $z=1000m$},$ $\displaystyle u_{x}=u_{y}=u_{z}=1m/s,$
$\displaystyle k_{x}=k_{y}=k_{z}=0.00002m^{2}/s.$

For $\delta x=\delta y=\delta z=10m$, each cell’s volume is $10m^{3}$ we
consider a time step of $1s$, due to the First-Order Upwind scheme stability
constraint known as the CFL (Courant-Friedrichs-Lewy) condition:
$\max\left(\frac{u_{x}\delta t}{\delta x},\frac{u_{y}\delta t}{\delta
y},\frac{u_{z}\delta t}{\delta z}\right)\leq\alpha$, a certain range less than
$1$. Dirichlet boundary conditions are considered at $x,y,z=0$ and $1000m$.
Below are Figures 4-6, the horizontal slice plots at $z=1$ of the distribution
of pollutants $\mathrm{NO}_{2},\mathrm{O}_{3},$ and $\mathrm{NO}$
concentrations at times $0s,20s,50s,100s,200s,300s,400s,500s$, and $600s$ with
$u=1m/s$ and $k=0.00002m^{2}/s$. The plots reveal the concentration
distribution of the dispersed gases across the horizontal plane at $z=1$. They
depict areas of high and low concentration, illustrating regions where the gas
is more or less prevalent. The multiple slices from different time steps
illustrate how the gas concentration evolves over time due to advection,
diffusion, and reactions. Considering $u_{x}=u_{y}=u_{z}=1m/s$ air currents
impact the gas transport. Concentration gradients show how the gas is being
transported from regions of higher concentration to lower concentration due to
advection. Diffusion is apparent in the gradual spreading of the gas
concentration away from its sources. The plots show how the gas diffuses and
disperses. In Figure 4, the maximum concentrations of $NO_{2}$ goes from
$5\times 10^{18}$ molecules$/10m^{3}$ at $t=0s$ to $74.385$ $mol./10m^{3}$ at
$200s$ and $0.078$ $mol./10m^{3}$ at $300s$, to $1.57\times 10^{-}19$
$mol./10m^{3}$ at $600s$. In Figure 5, the maximum concentrations of $O_{3}$
goes from $8\times 10^{18}$ $mol./10m^{3}$ at $t=0s$ to $119.017$
$mol./10m^{3}$ at $200s$ and $0.126$ $mol./10m^{3}$ at $300s$, to $2.52\times
10^{-19}$ $mol./10m^{3}$ at $600s$. In Figure 6, the maximum concentrations of
$NO$ is $1.3\times 10^{15}$ $mol./10m^{3}$ at $t=0s$, then $3.3\times 10^{13}$
at $200s$, at $300s$ and at $600s$. In Figure 6, there is a constant source of
pollutant $NO$ at cell $(1,1,1)$, therefore it is normal that the
concentration plot exhibits a relatively uniform concentration, with distinct
variations around the initial source point due to interactions and dispersion
processes. Close to the source point, higher concentrations can be observed
due to the immediate release of $NO$.

Figure 4: Horizontal slice plots at $z=1$ of $\mathrm{NO}_{2}$ concentrations
at times $0s,20s,50s,100s,200s,300s,400s,500s$, and $600s$ with $u=1m/s$ and
$k=0.00002m^{2}/s$. Figure 5: Horizontal slice plots at $z=1$ of
$\mathrm{O}_{3}$ concentrations at times
$0s,20s,50s,100s,200s,300s,400s,500s$, and $600s$ with $u=1m/s$ and
$k=0.00002m^{2}/s$. Figure 6: Horizontal slice plots at $z=1$ of $\mathrm{NO}$
concentrations at times $0s,20s,50s,100s,200s,300s,400s,500s$, and $600s$ with
$u=1m/s$ and $k=0.00002m^{2}/s$.

###  7.3 Attractor visualization

In Figure 7, a compelling observation arises as we discern the convergence of
all concentrations towards a shared region, accompanied by noticeable
oscillations within its confines. This distinct region, identified as our
attractor, stands out prominently against the backdrop of surrounding data
points. Remarkably, it encapsulates the trajectories of the concentrations of
each species, in each grid cell, in function of the other two species through
time, serving as a central hub for their dynamic behavior. This emphasizes the
pivotal role played by the attractor in dictating the enduring patterns and
long-term dynamics of the system, shedding light on its significant influence
on the overall system behavior.

Figure 7: The long-term evolution of the trajectories of the concentration of
each species in each cell of the grid.

##  8 Conclusion

In this work, several results concerning the existence, uniqueness, and
positivity of solutions as well as the asymptotic behavior through a finite
fractal dimensional global attractor for an Advection-Diffusion-Reaction (ADR)
equation, have been established. To achieve this goal, various results from
semigroup and global attractors theories were used. An analytical solution in
2-D was presented and numerical simulations were conducted using two Explicit
Finite Difference schemes: the Explicit Centered scheme to solve the two-
dimensional Advection-Diffusion equation and the First Order Upwind numerical
scheme to solve the three-dimensional Advection-Diffusion-Reaction equation,
which includes reactions involving $\mathrm{NO}$, $\mathrm{NO}_{2}$,
$\mathrm{O}_{3}$, and $\mathrm{O}_{2}$, with specified initial and boundary
conditions. The resulting plots illustrate the solutions of the two-
dimensional equation and later the concentration distribution of the dispersed
gases across the horizontal plane at $z=1$.

## CRediT Authors Contribution Statement

M. ELGHANDOURI:

Mathematical proofs, manuscript writing, review, and editing.

K. EZZINBI:

Mathematical methodology and supervision, manuscript review.

L. SAIDI:

Analytical solution, numerical simulation, and manuscript writing.

M. ELGHANDOURI & L. SAIDI:

Discussion and attractor visualization.

## Financial disclosure.

The authors assert that they do not possess any identifiable conflicting
financial interests or personal relationships that might have seemed to impact
the research presented in this paper.

## Conflict of interest.

This work does not have any conflicts of interest.

## References

  * [1] A. Batkai, M. K. Fijavz, et A. Rhandi. Positive operator semigroups. Operator Theory: advances and applications, 2017, vol. 257.
  * [2] J. Bear. Hydraulics of Groundwater, Dover, Minneola, 2007.
  * [3] D. D. Bainov and P. S. Simeonov. Integral inequalities and applications, vol. 57, Springer Science and Business Media, 2013.
  * [4] N. M. Buckman. The linear convection-diffusion equation in two dimensions.
  * [5] Burden, R., Faires, J., Burden, A.: Numerical analysis, 10th edn. Cengage Learning Inc., Boston (2015).
  * [6] R. Czaja et M. Efendiev. A note on attractors with finite fractal dimension. Bulletin of the London Mathematical Society, 2008, vol. 40, no 4, https://doi.org/10.1112/blms/bdn044.
  * [7] K. Elamvazhuthi, H. Kuiper, M. Kawski, et al. Bilinear controllability of a class of advection–diffusion–reaction systems. IEEE Transactions on Automatic Control, 2018, vol. 64, no 6, p. 2282-2297, https://doi.org/10.1109/TAC.2018.2885231.
  * [8] M. Efendiev, Finite and Infinite Dimensional Attractors for Evolution Equations of Mathematical Physics, Gakuto International Series. Mathematical Sciences and Applications, vol. 33, Gakkotosho Co., Ltd., Tokyo, 2010.
  * [9] M. T. Van Genuchten, W. J. Alves. Analytical Solutions of the One-Dimensional Convective-Dispersive Solute Transport Equation, U.S. Department of Agriculture, Agricultural Research Service Technical Bulletin, vol. 1661: Government Printing Office; 1982.
  * [10] B. Gustafson, H. Kreiss, and J. Oliger. Time Dependent Problems and Difference Methods, John Wiley & Sons, New York, 1995.
  * [11] J. Hale. Asymptotic behavior of dissipative systems (providence, ri: American mathematical society) go to reference in article (1988).
  * [12] D. K. Hetrick. Dynamics of Nuclear Reactors, University of Chicago, Chicago, 1971.
  * [13] W. H. Hundsdorfer, and G. J. Verwer. Numerical solution of time-dependent advection-diffusion-reaction equations. Vol. 33. Berlin: Springer, 2003.
  * [14] Kim, A.S. Complete analytic solutions for convection-diffusion-reaction-source equations without using an inverse Laplace transform. Sci Rep 10, 8040 (2020), https://doi.org/10.1038/s41598-020-63982-w.
  * [15] D. Lanser et G. J. Verwer. Analysis of operator splitting for advection–diffusion–reaction problems from air pollution modelling. Journal of computational and applied mathematics, 1999, vol. 111, no 1-2, p. 201-216, https://doi.org/10.1016/S0377-0427(99)00143-0.
  * [16] P. E. Lewis and J. P. Ward. The Finite Element Method: Principles and Applications (Addison-Wesley, Wokingham, 1991).
  * [17] W. Mclean, K. Mustapha, R. Ali et al. Well-posedness of time-fractional advection-diffusion-reaction equations. Fractional Calculus and Applied Analysis, 2019, vol. 22, no 4, p. 918-944, https://doi.org/10.1515/fca-2019-0050.
  * [18] N. Mphephu. Numerical solution of 1-D convection-diffusion-reaction equation. M.S. thesis, University of Venda, African Institute for Mathematical Sciences, 2013.
  * [19] J. D. Murray. Mathematical Biology I, Springer-Verlag, Berlin, 2002.
  * [20] Parhizi M., Kilaz G., Ostanek J. K., Jain A., Analytical solution of the convection-diffusion-reaction-source (CDRS) equation using Green’s function technique, International Communications in Heat and Mass Transfer, Volume 131, 2022, 105869, ISSN 0735-1933, https://doi.org/10.1016/j.icheatmasstransfer.2021.105869.
  * [21] A. Pazy. Semigroups of linear operators and applications to partial differential equations. Vol. 44. Springer Science & Business Media, 2012.
  * [22] L. Peng, Y. Zhou. The existence of mild and classical solutions for time fractional Fokker–Planck equations. Monatsh Math 199, 377–410 (2022), https://doi.org/10.1007/s00605-022-01710-4.
  * [23] A. D. Polyanin, V. E. Nazaikinskii. Handbook of Linear Partial Differential Equations for Engineers and Scientist. CRC Press Taylor & Francis Group LLC; 2016.
  * [24] H. Schmidt, C. Derognat, R. Vautard et al. A comparison of simulated and observed ozone mixing ratios for the summer of 1998 in Western Europe. Atmospheric Environment, 2001, vol. 35, no 36, p. 6277-6297, https://doi.org/10.1016/S1352-2310(01)00451-4.
  * [25] T. M. Shih. Numerical Heat Transfer, Springer-Verlag, Berlin, 1984.
  * [26] C. W. Shu. TVD time discretizations, SIAM J. Sci. Stat. Comput., 9 (1988), pp. 1073-1084, https://doi.org/10.1137/090907.
  * [27] J. I. Steinfeld. Atmospheric chemistry and physics: from air pollution to climate change. Environment: Science and Policy for Sustainable Development, 1998, vol. 40, no 7, p. 26-26, https://doi.org/10.1080/00139157.1999.10544295.
  * [28] A. U. Zavaleta, A. L. de Bortoli, M. Thompson. Analysis and simulation for a system of chemical reaction equations with a vortex formulation, Applied Numerical Mathematics, Volume 47, Issues 3–4, 2003, Pages 559-573, https://doi.org/10.1016/S0168-9274(03)00083-7.
  * [29] Z. Zhang, M. V. Tretyakov, B. Rozovskii et al. Wiener chaos versus stochastic collocation methods for linear advection-diffusion-reaction equations with multiplicative white noise. SIAM Journal on Numerical Analysis, 2015, vol. 53, no 1, p. 153-183, https://doi.org/10.1137/130932156.
  * [30] M. F. Wheeler, et C. N. Dawson. An operator-splitting method for advection-diffusion-reaction problems. 1987.