**1. Introduction**

Partial differential equations with time-fractional derivatives are suitable to describe complex problems in chaotic dynamics, signal processing, wave propagation, classical and continuum mechanics [1–6]. These equations admit the non-local memory effects; therefore, the fractional models can eliminate the disadvantage of classical models that cannot satisfactorily describe some complex processes such as anomalous diffusion.

In recent years, fractional advection–diffusion equations were intensively studied due to their multiple applications in many practical problems such as transport of pollutants in air and water, contaminating streams in groundwater, migration of pollutants in rivers and seawater, and tracer dispersion in porous media.

Bhrawy and Baleanu [7] developed an efficient Legendre–Gauss–Lobatto method to solve the advection–diffusion equation with a space-fractional Caputo derivative and non-homogeneous initial-boundary conditions. The advection–dispersion equation with a fractional order dependent on time and the Coimbra derivative was numerically solved by Abdelkawy et al. [8]. Zhuang et al. [9] numerically solved the advection–diffusion equation with derivatives of variable fractional order and a nonlinear source. They presented the explicit and implicit Euler approximation, fractional method of line, matrix transfer technique and the extrapolation method. Fazio et al. [10] have analyzed the advection–diffusion equation with the time-fractional Caputo derivative by using a finite difference method with non-uniform grids. They obtained an accuracy method compared with the classical discrete fractional formula for the time-fractional Caputo derivative. Time–space fractional advection–diffusion equations with Riesz derivatives have been numerically investigated by Arshad et al. [11]. A fractional advection–diffusion equation with a nonlinear source has been studied by Jannelli et al. [12]. They obtained a numerical solution and an analytical solution based on the Mittag–Leffler functions. Using the finite volume element method, Badr et al. [13] have studied the fractional advection–diffusion equation with the time-fractional Caputo derivative. Stability of the method has been analyzed. Pimenov [14] has numerically investigated the fractional

advection–diffusion equation with heredity. Singh et al. [15] have obtained analytical solutions for the fractional advection–dispersion equation by employing the q-fractional homotopy analysis transform method. Using the integral transform method, Avci and Yetim [16] have found the fundamental solutions of the fractional advection–diffusion with fractional derivatives of the Atangana–Baleanu type. Fundamental solutions for the Dirichlet problem of the two-dimensional fractional advection–diffusion equation with short-tail memory have been studied by Mirza and Vieru [17]. Using the Laplace transform and the sine–cosine Fourier transform, Mirza et al. [18] have determined analytical solutions for the fractional advection–diffusion equation in the one-dimensional case and with boundary conditions of Robin-type. Interesting fractional diffusive processes have been investigated by Hristov in [19–21]. Povstenko and Klekot [22] and Povstenko and Kyrylych [23] have studied the fractional advection diffusion equation with several boundary conditions. It is important to underline that the time–space fractional advection–diffusion equation could be solved with the Monte–Carlo method as a limit of continuous-time random walks [24].

In this work, the analytical and numerical study of two-dimensional advection–diffusion processes with short-tail memory and with the source concentrated in the symmetry center of the geometric domain is carried out.

The mathematical model is described by a fractional differential equation with the time-fractional Caputo–Fabrizio derivative. An analytical solution of the initial-boundary value problem has been determined by employing the Laplace and finite sine-Fourier transforms. A numerical solution of the studied problem has been determined using finite difference approximations. Numerical simulations for both solutions have been carried out using the software Mathcad and a good agreemen<sup>t</sup> is found.

#### **2. The Statement of the Problem**

Let us consider the two-dimensional advection–diffusion process with memory and the time-dependent concentrated-source in a finite domain:

*D*= (*<sup>x</sup>*, *y*), *a* ≤ *x* ≤ *b*, *c* ≤ *y* ≤ *d*, *a* < *b*, *c* < *d* with the border ∂*D*.

The mathematical model is described by the following time-fractional advection–diffusion equation [7,9]:

$$\begin{cases} \,^c \mathbf{F} D\_t^\alpha u(\mathbf{x}, y, t) + k\_1 \frac{\partial u(\mathbf{x}, y, t)}{\partial \mathbf{x}} + k\_2 \frac{\partial u(\mathbf{x}, y, t)}{\partial y} - \varepsilon\_1 \frac{\partial^2 u(\mathbf{x}, y, t)}{\partial \mathbf{x}^2} - \varepsilon\_2 \frac{\partial^2 u(\mathbf{x}, y, t)}{\partial y^2} = f(\mathbf{x}, y, t), \\ 0 < a \le 1, \quad k\_1, \; k\_2 \ge 0, \; \varepsilon\_1, \; \varepsilon\_2 > 0, \; (\mathbf{x}, y) \in D. \end{cases} \tag{1}$$

where *<sup>u</sup>*(*<sup>x</sup>*, *y*,*<sup>t</sup>*) is the field variable representing the solute concentration, →*v* = *<sup>k</sup>*1<sup>→</sup>*<sup>e</sup> x* + *<sup>k</sup>*2<sup>→</sup>*<sup>e</sup> y* is the constant fluid velocity (the drift velocity), ε1, ε2 are the dispersion coefficients in the x-direction and y-direction, respectively and *f*(*<sup>x</sup>*, *y*, *t*) represents the sources/sinks.

The operator *CFD*<sup>α</sup>*t*(·) denotes the time-fractional Caputo–Fabrizio derivative defined as [25]:

$$\begin{cases} \,^G D\_t^\alpha u(x, y, t) = h(\alpha, t) \* \dot{u}(x, y, t), \\ h(\alpha, t) = \begin{cases} \,^M \frac{\dot{M}(\alpha)}{1 - \alpha} \exp\left(\frac{-\alpha t}{1 - \alpha}\right) & , 0 \le \alpha < 1, \\ \delta(t), \alpha = 1, \end{cases} \end{cases} \tag{2}$$
 
$$\dot{u}(x, y, t) = \frac{\partial u(x, y, t)}{\partial t},$$

where "∗" denotes the convolution product and δ(*t*) is Dirac's distribution. In the above definition, the function*M*(α)is so called the normalization function, which should satisfy conditions*M*(0) = *M*(1) = 1. In the present paper, for simplicity, we will consider for the normalization function *<sup>M</sup>*(α) ≡ 1, α ∈ [0, 1].

The following properties of the Caputo–Fabrizio operator immediately result from Equation (2):

$$\mathbf{u}^{\rm CF} D\_t^0 u(\mathbf{x}, y, t) = u(\mathbf{x}, y, t) - u(\mathbf{x}, y, 0), \\ \mathbf{^{CF}} D\_t^1 u(\mathbf{x}, y, t) = \frac{\partial u(\mathbf{x}, y, t)}{\partial t}, \\ \mathbf{^{CF}} D\_t^d \mathbf{C} = \mathbf{0}, \\ \mathbf{C} = \text{const.} \tag{3}$$

The Laplace transform of the Caputo–Fabrizio derivative is given by:

$$\mathbb{L}\left\{{}^{\text{CF}}D\_{t}^{a}u(\mathbf{x},y,t)\right\} = \mathbb{L}\left\{h(a,t)\*\dot{u}(\mathbf{x},y,t)\right\} = \mathbb{L}\{h(a,t)\}\mathbb{L}\{\dot{u}(\mathbf{x},y,t)\} = \frac{s\mathbb{L}\{u(\mathbf{x},t,t)\} - u(\mathbf{x},y,0)}{(1-a)s+a}.\tag{4}$$

The fractional advection–diffusion equation (Equation (1)), is studied along with following initial and boundary conditions:

$$u(\mathbf{x}, y, \mathbf{0}) = \psi\_0(\mathbf{x}, y), \ (\mathbf{x}, y) \in D,\tag{5}$$

$$
\mu(\mathbf{x}, y, t) = \psi\_1(\mathbf{x}, y, t), \ t > 0, \ (\mathbf{x}, y) \in \partial D, \ t > 0,\tag{6}
$$

Making the transformation:

$$u(\mathbf{x}, y, t) = e^{\oint \mathbf{(x, y)} \cdot \mathbf{(x, y, t)} + \psi\_0(\mathbf{x, y}), \ f\_0(\mathbf{x, y}) = \frac{k\_1 x}{2\varepsilon\_1} + \frac{k\_2 y}{2\varepsilon\_2},\tag{7}$$

Equation (1) becomes:

$$\mathbf{F}^{\rm CF} D\_t^{\rm a} \boldsymbol{\varphi}(\mathbf{x}, \mathbf{y}, t) + \frac{1}{2} \Big( \frac{k\_1^2}{\varepsilon\_1} + \frac{k\_2^2}{\varepsilon\_2} \Big) \boldsymbol{\varphi}(\mathbf{x}, \mathbf{y}, t) - \varepsilon\_1 \frac{\partial^2 \boldsymbol{\varphi}(\mathbf{x}, \mathbf{y}, t)}{\partial \mathbf{x}^2} - \varepsilon\_2 \frac{\partial^2 \boldsymbol{\varphi}(\mathbf{x}, \mathbf{y}, t)}{\partial \mathbf{y}^2} = \mathbf{F}(\mathbf{x}, \mathbf{y}, t) + \mathbf{G}(\mathbf{x}, \mathbf{y}), \tag{8}$$

where

$$\begin{aligned} F(\mathbf{x}, y, t) &= f(\mathbf{x}, y, t) \exp(-f\_0(\mathbf{x}, y)) \\ G(\mathbf{x}, y) &= \left[ \varepsilon\_1 \frac{\partial^2 \psi\_0(\mathbf{x}, y)}{\partial \mathbf{x}^2} + \varepsilon\_2 \frac{\partial^2 \psi\_0(\mathbf{x}, y)}{\partial y^2} - k\_1 \frac{\partial \psi\_0(\mathbf{x}, y)}{\partial \mathbf{x}} - k\_2 \frac{\partial \psi\_0(\mathbf{x}, y)}{\partial y} \right] \exp(-f\_0(\mathbf{x}, y)) .\end{aligned} \tag{9}$$

After transformation (7), the initial and boundary conditions become:

$$
\varphi(\mathbf{x}, y, \mathbf{0}) = \mathbf{0}, \ (\mathbf{x}, y) \in D, \tag{10}
$$

$$\varphi(\mathbf{x}, y, t) = F\_1(\mathbf{x}, y, t) + G\_1(\mathbf{x}, y), \ (\mathbf{x}, y) \in \partial D, \ t > 0,\tag{11}$$

where

$$F\_1(\mathbf{x}, y, t) = \psi\_1(\mathbf{x}, y, t) \exp(-f\_0(\mathbf{x}, y)),\\\ G\_1(\mathbf{x}, y) = -\psi\_0(\mathbf{x}, y) \exp(-f\_0(\mathbf{x}, y)).\tag{12}$$

#### **3. Analytical Solution of the Problem**

The solution of the problem described by Equation (8) and conditions (10) and (11) will be determined by using the integral transform method.

Applying the Laplace transform to Equation (8), using Equation (4) and Equation (10), we obtain the transformed equation:

$$\frac{[1+\beta(1-a)]s+a\beta}{(1-a)s+a}\overline{\varphi}(\mathbf{x},y,\mathbf{s})-\varepsilon\_1\frac{\partial^2\overline{\varphi}(\mathbf{x},y,\mathbf{s})}{\partial\mathbf{x}^2}-\varepsilon\_2\frac{\partial^2\overline{\varphi}(\mathbf{x},y,\mathbf{s})}{\partial y^2} = \overline{F}(\mathbf{x},y,\mathbf{s})+\frac{1}{s}G(\mathbf{x},y),\tag{13}$$

where β = 12 *k*21ε1 + *k*22ε2 and ϕ(*<sup>x</sup>*, *y*,*<sup>s</sup>*) = #∞0 ϕ(*<sup>x</sup>*, *y*, *<sup>t</sup>*)e<sup>−</sup>*st*dt, *<sup>F</sup>*(*<sup>x</sup>*, *y*,*<sup>s</sup>*) = #∞0 *<sup>F</sup>*(*<sup>x</sup>*, *y*, *<sup>t</sup>*)e<sup>−</sup>*st*dt are the Laplace transformsoffunctionsϕ(*<sup>x</sup>*,*y*,*t*),*<sup>F</sup>*(*<sup>x</sup>*,*y*,*t*),respectively.

 Equation(13)hastosatisfytheboundarycondition:

$$\begin{cases} \overline{\boldsymbol{\varphi}}(\mathbf{x}, \mathbf{y}, \mathbf{s}) = \overline{\boldsymbol{\varphi}}\_{1}(\mathbf{x}, \mathbf{y}, \mathbf{s}), (\mathbf{x}, \mathbf{y}) \in \partial \mathcal{D},\\ \overline{\boldsymbol{\varphi}}\_{1}(\mathbf{x}, \mathbf{y}, \mathbf{s}) = \overline{\boldsymbol{\Gamma}}\_{1}(\mathbf{x}, \mathbf{y}, \mathbf{s}) + \frac{1}{s} \mathbf{G}\_{1}(\mathbf{x}, \mathbf{y}). \end{cases} \tag{14}$$

*Symmetry* **2019**, *11*, 879

Now, we apply this to Equation (13), the double sine-Fourier transform with respect to variables *x* and *y* defined as [26]:

$$\widetilde{\overline{\varphi}}(n,k,s) = \int\_{c}^{d} \int\_{a}^{b} \sin\left[\frac{n(\mathbf{x}-a)\pi}{b-a}\right] \sin\left[\frac{k(y-c)\pi}{d-c}\right] \overline{\varphi}(\mathbf{x},y,s) dx dy, \ n, \ k = 1,2,\ldots \tag{15}$$

A direct calculation leads to the following relationships:

$$\int\_{c}^{d} \int\_{a}^{b} \frac{\partial^2 \overline{\varphi}(x, y, s)}{\partial x^2} \sin \left[ \frac{n(x - a)\pi}{b - a} \right] \sin \left[ \frac{k(y - c)\pi}{d - c} \right] dx dy = \overline{H}\_1(n, k, s) - \left( \frac{n\pi}{b - a} \right)^2 \overline{\overline{\varphi}}(n, k, s), \tag{16}$$

where

$$\tilde{H}\_1(n,k,s) = \frac{n\pi}{b-a} \int\_c^d \left[ \overline{\varphi}\_1(a,y,s) - (-1)^n \overline{\varphi}\_1(b,y,s) \right] \sin\left[\frac{k(y-c)\pi}{d-c}\right] dy \tag{17}$$

$$\int\_{a}^{b} \int\_{c}^{d} \frac{\partial^2 \overline{\varphi}(x, y, s)}{\partial y^2} \sin \left[ \frac{n(x - a)\pi}{b - a} \right] \sin \left[ \frac{k(y - c)\pi}{d - c} \right] dy d\mathbf{x} = \widetilde{\Pi}\_2(n, k, s) - \left( \frac{k\pi}{d - c} \right)^2 \widetilde{\overline{\varphi}}(n, k, s) \tag{18}$$

where

$$\widetilde{H}\_2(n,k,s) = \frac{k\pi}{d-c} \int\_a^b \left[\overline{\varphi}\_1(\mathbf{x},c,s) - (-1)^k \overline{\varphi}\_1(\mathbf{x},d,s)\right] \sin\left[\frac{n(\mathbf{x}-a)\pi}{b-a}\right] d\mathbf{x}.\tag{19}$$

Applying the double sine-Fourier transform to Equation (13), using the relationships (16)–(19) and notations 

$$\begin{aligned} a\_{nk} &= 1 + (1 - \alpha) \left[ \beta + \varepsilon\_1 \left( \frac{n\pi}{b - a} \right)^2 + \varepsilon\_2 \left( \frac{k\pi}{d - c} \right)^2 \right] \\ b\_{nk} &= \alpha \left[ \beta + \varepsilon\_1 \left( \frac{n\pi}{b - a} \right)^2 + \varepsilon\_2 \left( \frac{k\pi}{d - c} \right)^2 \right], n, k = 1, 2, \dots, \\ \widetilde{\overline{\Psi}}(n, k, s) &= \widetilde{\overline{F}}(n, k, s) + \frac{1}{s} \widetilde{G}(n, k) + \varepsilon\_1 \widetilde{H}\_1(n, k, s) + \varepsilon\_2 \widetilde{H}\_2(n, k, s), \end{aligned} \tag{20}$$

the transform ϕ(*<sup>n</sup>*, *k*,*<sup>s</sup>*) is written in the following form:

9

$$
\widetilde{\overline{\varphi}}(n,k,s) = \frac{(1-a)s + a}{a\_{nk}s + b\_{nk}} \widetilde{\overline{\Psi}}(n,k,s). \tag{21}
$$

Now, we will write Equation (21) in a suitable form which, after inversion, satisfies the initial and boundary conditions (10) and (11). To do this, we consider the following auxiliary functions:

$$\begin{array}{lcl} \text{v}\_{i}(\boldsymbol{y},t) &= [\boldsymbol{g}\_{i}(\boldsymbol{y},t) - \boldsymbol{g}\_{i}(d,t)]H(\boldsymbol{y}-\boldsymbol{c}) + [\boldsymbol{g}\_{i}(\boldsymbol{y},t) - \boldsymbol{g}\_{i}(\boldsymbol{c},t)]H(d-\boldsymbol{y}) - \\ & - \boldsymbol{g}\_{i}(\boldsymbol{y},t)H(\boldsymbol{y}-\boldsymbol{c})H(d-\boldsymbol{y}); \quad i=1,2, \ \boldsymbol{y} \in [\boldsymbol{c},d], \end{array} \tag{22}$$

$$\begin{array}{lcl}\mathbf{w}\_{i}(\mathbf{x},t) = [h\_{i}(\mathbf{x},t) - h\_{i}(b,t)]H(\mathbf{x}-a) + [h\_{i}(\mathbf{x},t) - h\_{i}(a,t)]H(b-\mathbf{x}) - \\ -h\_{i}(\mathbf{x},t)H(\mathbf{x}-a)H(b-\mathbf{x}); \quad i=1,2, \quad \mathbf{x}\in[a,b],\end{array} \tag{23}$$

where

$$\begin{aligned} g\_1(y,t) &= F\_1(a,y,t) + G\_1(a,y), \; g\_2(y,t) = F\_1(b,y,t) + G\_1(b,y), \\ h\_1(\mathbf{x},t) &= F\_1(\mathbf{x},c,t) + G\_1(\mathbf{x},c), \; h\_2(\mathbf{x},t) = F\_1(\mathbf{x},d,t) + G\_1(\mathbf{x},d), \end{aligned} \tag{24}$$

and

$$H(z) = \frac{\text{sign}(z)[1 + \text{sign}(z)]}{2} = \begin{cases} 0, & z \le 0 \\ 1, & z > 0 \end{cases} \tag{25}$$

is the Heaviside step unit function. With the above functions, we define:

$$\theta\_1(\mathbf{x}, y, t) = \frac{b - \mathbf{x}}{b - a} \mathbf{v}\_1(y, t) + \frac{\mathbf{x} - a}{b - a} \mathbf{v}\_2(y, t) = \frac{\mathbf{x}[\mathbf{v}\_2(y, t) - \mathbf{v}\_1(y, t)]}{b - a} + \frac{b \mathbf{v}\_1(y, t) - a \mathbf{v}\_2(y, t)}{b - a},\tag{26}$$

$$\theta\_2(\mathbf{x}, y, t) = \frac{d - y}{d - c} \mathbf{w}\_1(\mathbf{x}, t) + \frac{y - c}{d - c} \mathbf{w}\_2(\mathbf{x}, t) = \frac{y[\mathbf{w}\_2(\mathbf{x}, t) - \mathbf{w}\_1(\mathbf{x}, t)]}{d - c} + \frac{d \mathbf{w}\_1(\mathbf{x}, t) - c \mathbf{w}\_2(\mathbf{x}, t)}{d - c},\tag{27}$$

$$\begin{array}{ll} \theta\_3(\mathbf{x}, y, \mathbf{t}) = & h\_1(a, t)H(b - \mathbf{x})H(d - y) + h\_1(b, t)H(\mathbf{x} - a)H(d - y) + \\ & h\_2(a, t)H(b - \mathbf{x})H(y - \mathbf{c}) + h\_2(b, t)H(\mathbf{x} - a)H(y - \mathbf{c}), \end{array} \tag{28}$$

The function <sup>θ</sup>(*<sup>x</sup>*, *y*, *t*) defined by

$$
\partial(\mathbf{x}, y, t) = \partial\_1(\mathbf{x}, y, t) + \partial\_2(\mathbf{x}, y, t) + \partial\_3(\mathbf{x}, y, t), \tag{29}
$$

satisfies the following properties:

$$\theta(\theta(\mathbf{z}, y, t) = \operatorname{g}\_1(y, t), \ \theta(b, y, t) = \operatorname{g}\_2(y, t), \ \theta(\mathbf{x}, \mathbf{c}, t) = h\_1(\mathbf{x}, t), \ \theta(\mathbf{x}, d, t) = h\_2(\mathbf{x}, t) \tag{30}$$

Applying the Laplace transform and double sine-Fourier transform to functions (26)–(28), we obtain: 9

$$\begin{split} \overline{\partial}\_{1}(n,k,s) &= \frac{a-(-1)^{n}b}{n\pi}P\_{1k}(s) + \frac{1-(-1)^{n}}{n\pi}P\_{2k}(s), \\ P\_{1k}(s) &= \int \left[\mathbf{\dot{v}}\_{2}(y,s) - \mathbf{\dot{v}}\_{1}(y,s)\right] \sin\left[\frac{k(y-c)\pi}{d-c}y\right] dy, \\ P\_{2k} &= \int \left[\mathbf{\dot{v}}\_{1}(y,s) - a\mathbf{v}\_{2}(y,s)\right] \sin\left[\frac{k(y-c)\pi}{d-c}y\right] dy, \\ \overline{\partial}\_{2}(n,k,s) &= \frac{c-(-1)^{k}d}{k\pi}Q\_{1n}(s) + \frac{1-(-1)^{k}}{k\pi}Q\_{2n}(s), \\ Q\_{1n}(s) &= \int \left[\mathbf{w}\_{2}(x,s) - \mathbf{w}\_{1}(x,s)\right] \sin\left[\frac{n(x-a)\pi}{b-a}\right] dx, \\ Q\_{2n}(s) &= \int \left[d\mathbf{w}\_{1}(x,s) - c\mathbf{w}\_{2}(x,s)\right] \sin\left[\frac{n(x-a)\pi}{b-a}\right] dx, \\ &\qquad\qquad \left(1-(-1)^{\frac{n}{k}}\right)\left(1-(-1)^{\frac{k}{k}}\right) \end{split} \tag{32}$$

$$\widetilde{\overline{\Theta}}\_{3}(n,k,s) = \left[\overline{h}\_{1}(a,s) + \overline{h}\_{1}(b,s) + \overline{h}\_{2}(a,s) + \overline{h}\_{2}(b,s)\right] \left(\frac{1-(-1)^{n}}{n\pi}\right) \left(\frac{1-(-1)^{k}}{k\pi}\right) (b-a)(d-c) \tag{33}$$

respectively,

$$\begin{split} \widetilde{\overline{\Theta}}(n,k,s) &= \frac{a-(-1)^{n}b}{n\pi} P\_{1k}(s) + \frac{1-(-1)^{n}}{n\pi} P\_{2k}(s) + \frac{c-(-1)^{k}d}{k\pi} Q\_{1n}(s) + \frac{1-(-1)^{k}}{k\pi} Q\_{2n}(s) + \\ & \left[ \overline{h}\_{1}(a,s) + \overline{h}\_{1}(b,s) + \overline{h}\_{2}(a,s) + \overline{h}\_{2}(b,s) \right] \left( \frac{1-(-1)^{n}}{n\pi} \right) \left( \frac{1-(-1)^{k}}{k\pi} \right) (b-a) (d-c). \end{split} \tag{34}$$

The function ϕ, given by Equation (21), can be written as:

9

$$
\widetilde{\overline{\boldsymbol{\sigma}}}(\boldsymbol{n},k,\boldsymbol{s}) = \widetilde{\overline{\boldsymbol{\Theta}}}(\boldsymbol{n},k,\boldsymbol{s}) + \frac{1}{a\_{\text{nk}}\mathbf{s} + b\_{\text{nk}}} \widetilde{\overline{\mathbf{\overline{\boldsymbol{\Theta}}}}}\_{1}(\boldsymbol{n},k,\mathbf{s}), \tag{35}
$$

where

$$
\widetilde{\overline{\Psi}}\_1(n,k,s) = [(1-a)s+a]\widetilde{\overline{\Psi}}(n,k,s) - (a\_{nk}s+b\_{nk})\widetilde{\overline{\Theta}}(n,k,s).\tag{36}
$$

Applying the inverse Laplace and Fourier transforms, we obtain the solution

$$\begin{split} \, \, \rho(\mathbf{x}, y, t) &= H(t) \theta(\mathbf{x}, y, t) + \\ \, \, \frac{4}{(b-a)(d-c)} \sum\_{n=1}^{\infty} \sum\_{k=1}^{\infty} \sin \left[ \frac{n(\mathbf{x} - \mathbf{a}) \pi}{b - a} \right] \sin \left[ \frac{k(\mathbf{y} - \mathbf{c}) \pi}{d - c} \right] \int\_{0}^{t} \frac{\exp(-(t - \mathbf{r}) b\_{nk}/a\_{nk})}{a\_{nk}} \, \widetilde{\Psi}\_{1}(\mathbf{r}, k, \mathbf{r}) \, \mathrm{d}\tau. \end{split} \tag{37}$$

Using the property (30), it is easy to see that function (37) satisfies boundary condition (11). The solution for the classical advection equation is obtained by making α = 1 into Equations (20) and (35)–(37).

#### **4. Particular Case: A Time-Dependent Concentrated Source in the Center of Domain** *D*

In this section, for the source *f*(*<sup>x</sup>*, *y*, *t*) and for the initial and boundary conditions, we consider: 

$$\begin{array}{l} f(\mathbf{x}, y, t) = \gamma\_0 \delta \Big(\mathbf{x} - \frac{a + b}{2}\Big) \delta \Big(y - \frac{c + d}{2}\Big) \exp(-pt), \, p \ge 0, \\\ u(\mathbf{x}, y, 0) = \psi\_0(\mathbf{x}, y) = \gamma\_0 \exp(f\_0(\mathbf{x}, y)), \, (\mathbf{x}, y) \in D, \, \gamma\_0 = \text{const.} \\\ u(\mathbf{x}, y, t) = \psi\_1(\mathbf{x}, y, t) = 2\gamma\_0 \exp(f\_0(\mathbf{x}, y)), \, (\mathbf{x}, y) \in \partial D, \, t > 0. \end{array} \tag{38}$$

The physical significance of the above source *f*(*<sup>x</sup>*, *y*, *t*) is that of a source of intensity γ0 exp(−*p<sup>t</sup>*) concentrated in the center *a*+*b* 2 , *c*+*d*2 of the rectangular domain *D*.

Using (38), Equations (9), (11), (12) and (14) lead to:

$$\begin{aligned} F(\mathbf{x}, y, t) &= \gamma\_0 \delta \Big( \mathbf{x} - \frac{a + b}{2} \Big) \delta \Big( y - \frac{c + d}{2} \Big) \exp(-pt - f\_0(\mathbf{x}, y)), \ G(\mathbf{x}, y) = -\frac{\beta \gamma\_0}{2}, \\ F\_1(\mathbf{x}, y, t) &= 2\gamma\_0, \ G\_1(\mathbf{x}, y) = -\gamma\_0, \ \eta\_1(\mathbf{x}, y, t) = \gamma\_0. \end{aligned} \tag{39}$$

By using <sup>ϕ</sup>1(*<sup>x</sup>*, *y*,*<sup>s</sup>*) = γ0*s* in Equations (17) and (19), we have:

$$
\begin{array}{l}
\overline{\overline{H}}\_1(\mathfrak{n},k,\mathfrak{s}) = \frac{\mathcal{V}0}{\mathfrak{s}} \frac{d-c}{b-a} \frac{n[1-(-1)^n][1-(-1)^k]}{k}, \\
\overline{\overline{H}}\_2(\mathfrak{n},k,\mathfrak{s}) = \frac{\mathcal{V}0}{\mathfrak{s}} \frac{b-a}{d-c} \frac{k[1-(-1)^n][1-(-1)^k]}{n}.
\end{array}
\tag{40}
$$

Using the "sifting" property of the Dirac distribution [27]

$$\bigcup\_{a}^{b} f(\mathbf{x})\delta(\mathbf{x}-\mathbf{c})d\mathbf{x} = \begin{cases} \frac{1}{2}f(\mathbf{c}^{+}), \mathbf{c} = a, \\\frac{1}{2}[f(\mathbf{c}^{-}) + f(\mathbf{c}^{+})], \mathbf{c} \in (a, b), \\\frac{1}{2}f(\mathbf{c}^{-}), \mathbf{c} = b, \\\ 0, \mathbf{c} \in (a, b). \end{cases} \tag{41}$$

we find that the double sine-Fourier transform of the functions

$$\begin{cases} \overline{F}(\mathbf{x}, y, \mathbf{s}) = \frac{\gamma \upsilon}{s + p} \delta \left( \mathbf{x} - \frac{a + b}{2} \right) \delta \left( y - \frac{c + d}{2} \right) \exp(-f\_0(\mathbf{x}, y)),\\ G(\mathbf{x}, y) = \frac{-\beta \gamma \upsilon}{2}, \end{cases} \tag{42}$$

are given by

$$\begin{cases} \overline{F}(n,k,s) = \frac{\gamma\_0}{s+p} \sin\left(\frac{\eta \pi}{2}\right) \sin\left(\frac{k\pi}{2}\right) \exp\left(-\frac{k\_1(a+b)}{2x\_1^2} - \frac{k\_2(c+d)}{2x\_2^2}\right),\\ \overline{G}(n,k) = \frac{-\beta\gamma\_0}{2} \frac{(b-a)(d-c)}{kn\pi^2} [1 - (-1)^n][1 - (-1)^k]. \end{cases} \tag{43}$$

Since sin *m*π2 = (−<sup>1</sup>)*<sup>m</sup>*−<sup>1</sup> 2 *im*<sup>+</sup>1, *i*2 = −1, *m* ∈ N, the function 9*<sup>F</sup>*(*<sup>n</sup>*, *k*,*<sup>s</sup>*) is written as

$$\widetilde{\overline{F}}(n,k,s) = \frac{-\gamma\_0}{s+p} \exp\left(-\frac{k\_1(a+b)}{2\varepsilon\_1^2} - \frac{k\_2(c+d)}{2\varepsilon\_2^2}\right) \frac{[1-\left(-1\right)^n][1-\left(-1\right)^k]}{4} i^{n+k}.\tag{44}$$

By using Equations (40) and (43) in the third equation (20), we obtain:

$$\begin{array}{l} \widetilde{\overline{\Psi}}(n,k,s) = \frac{[1-(-1)^{n}][1-(-1)^{k}]}{4} \chi\\ \left[\frac{2\gamma\_{0}}{s}\frac{2\varepsilon\_{1}(n\pi)^{2}(d-c)^{2}+2\varepsilon\_{2}(k\pi)^{2}(b-a)^{2}-\beta(b-a)^{2}(d-c)^{2}}{(k\pi)(n\pi)(b-a)(d-c)} - \frac{\gamma\_{0}\gamma\_{1}i^{n+k}}{s+p}\right],\end{array} \tag{45}$$

*Symmetry* **2019**, *11*, 879

where

$$\gamma\_1 = \exp\left(-\frac{k\_1(a+b)}{4\varepsilon\_1} - \frac{k\_2(c+d)}{4\varepsilon\_2}\right)$$

In the considered particular case, Equations (22)–(24) and (31)–(34) become:

$$\begin{cases} \varrho\_1(y,t) = \varrho\_2(y,t) = \gamma\_{0\prime} v\_1(y,t) = v\_2(y,t) = -\gamma\_0 H(y-c)H(d-y), \\ h\_1(\mathbf{x},t) = h\_2(\mathbf{x},t) = \gamma\_{0\prime} w\_1(\mathbf{x},t) = w\_2(\mathbf{x},t) = -\gamma\_0 H(\mathbf{x}-a)H(b-\mathbf{x}), \end{cases} \tag{46}$$

$$\begin{array}{lcl}\theta\_1(\mathbf{x}, y, t) &=& -\gamma\_0 H(y - c) H(d - y), \,\theta\_2(\mathbf{x}, y, t) = -\gamma\_0 H(\mathbf{x} - a) H(b - \mathbf{x}),\\\theta\_3(\mathbf{x}, y, t) &=& \gamma\_0 H(b - \mathbf{x}) H(d - y) + \gamma\_0 H(\mathbf{x} - a) H(d - y) +\\ & & \gamma\_0 H(b - \mathbf{x}) H(y - c) + \gamma\_0 H(\mathbf{x} - a) H(y - c). \end{array} \tag{47}$$

The transformed Laplace–Fourier of the function (47) are given by:

$$\begin{array}{l} \overline{\overline{\Theta}}\_{1}(n,k,s) = \overline{\overline{\Theta}}\_{2}(n,k,s) = -\frac{1}{s}\gamma\_{0}(b-a)(d-c)\frac{1-(-1)^{n}}{n\pi}\frac{1-(-1)^{k}}{k\pi},\\ \overline{\overline{\Theta}}\_{3}(n,k,s) = \frac{4}{s}\gamma\_{0}(b-a)(d-c)\frac{1-(-1)^{n}}{n\pi}\frac{1-(-1)^{k}}{k\pi},\\ \overline{\overline{\Theta}}(n,k,s) = \frac{2}{s}\gamma\_{0}(b-a)(d-c)\frac{1-(-1)^{n}}{n\pi}\frac{1-(-1)^{k}}{k\pi}. \end{array} \tag{48}$$

Function <sup>Ψ</sup>1(*<sup>n</sup>*, *k*,*<sup>s</sup>*) given by Equation (36) can be written as:

$$\widetilde{\overline{\Psi}}\_1(n,k,s) = \mathbb{C}\_{nk}\frac{1}{s+p} + D\_{nk}\frac{1}{s} + E\_{nk\prime} \tag{49}$$

where

$$\begin{aligned} A\_{nk} &= \left[1 - (-1)^{n}\right] \left[1 - (-1)^{k}\right] \\ B\_{nk} &= \frac{1}{(b-a)(d-c)} \left[2\varepsilon\_{1}(n\pi)^{2}(d-c)^{2} + 2\varepsilon\_{2}(k\pi)^{2}(b-a)^{2} - \beta(b-a)^{2}(d-c)^{2}\right], \\ C\_{nk} &= \frac{1}{2}\gamma\_{0}\gamma\_{1}\left[(1-a)p - a\right]t^{n+k}A\_{nk}, \\ D\_{nk} &= \frac{\gamma\_{0}A\_{nk}}{2(n\pi)(k\pi)} \left[aB\_{nk} - 4(b-a)(d-c)b\_{nk}\right], \\ E\_{nk} &= \frac{A\_{nk}}{4(n\pi)(k\pi)} \left[2(1-a)\gamma\_{0}B\_{nk} - 8\gamma\_{00}(b-a)(d-c)a\_{nk} - (1-a)\gamma\_{0}\gamma\_{1}(n\pi)(k\pi)t^{n+k}\right], \end{aligned} \tag{50}$$

and transform (35) becomes:

9

$$\widetilde{\overline{\varphi}}(n,k,s) = \widetilde{\overline{\Theta}}(n,k,s) + \frac{1}{a\_{nk}} \frac{1}{s + b\_{nk}/a\_{nk}} \left[ \frac{\mathbb{C}\_{nk}}{\text{s} + p} + \frac{D\_{nk}}{\text{s}} + E\_{nk} \right] \tag{51}$$

Applying the inverse Laplace transform and the inverse Fourier transform to function (51), we obtain:

$$\begin{split} \rho(\mathbf{x}, y, t) &= H(t)\theta(\mathbf{x}, y, t) + \frac{4}{(b - a)(d - c)} \sum\_{n = 1}^{\infty} \sum\_{k = 1}^{\infty} \frac{1}{a\_{nk}} \Big\{ E\_{nk} \exp\left(\frac{-b\_{nk}t}{a\_{nk}}\right) + \\ &\int\_{0}^{t} \exp\left(\frac{-b\_{nk}(t - \tau)}{a\_{nk}}\right) (D\_{nk} + \mathbb{C}\_{nk} e^{-p\tau}) d\tau \right\} \sin\left(\frac{n(x - a)\pi}{b - a}\right) \sin\left(\frac{k(y - c)\pi}{d - c}\right), \end{split} \tag{52}$$

where the function <sup>θ</sup>(*<sup>x</sup>*, *y*, *t*) is given by

$$\begin{aligned} \theta(\mathbf{x}, y, t) &= \begin{array}{c} \gamma\_0 [H(\mathbf{x} - a) + H(b - \mathbf{x})][H(y - c) + H(d - y)] - \\ \gamma\_0 H(\mathbf{x} - a)H(b - \mathbf{x}) - \gamma\_0 H(y - c)H(d - y) .\end{array} \end{aligned} \tag{53}$$

Using the property

$$\mathcal{S}(t) = H(t - a\_1) + H(a\_2 - t) - H(t - a\_1)H(a\_2 - t) \equiv 1, \ t \in [a\_1, a\_2],\tag{54}$$

it is found that <sup>θ</sup>(*<sup>a</sup>*, *y*,*<sup>t</sup>*) = <sup>θ</sup>(*b*, *y*,*<sup>t</sup>*) = γ0 for *y* ∈ [*c*, *d*], respectively <sup>θ</sup>(*<sup>x</sup>*, *<sup>c</sup>*,*<sup>t</sup>*) = <sup>θ</sup>(*b*, *d*,*<sup>t</sup>*) = γ0 for *x* ∈ [*a*, *b*]; therefore, function (53) satisfied boundary condition (14). Now, to return to the concentration *<sup>u</sup>*(*<sup>x</sup>*, *y*, *t*), the function ϕ(*<sup>x</sup>*, *y*, *t*) given by Equation (52) will be replaced into transformation (7).

Figures 1–4 have been plotted to highlight the variation of the non-dimensional concentration ϕ(*<sup>x</sup>*, *y*, *t*) versus the fractional parameter α.

**Figure 1.** Profiles of function ϕ(*<sup>x</sup>*, 1.2,*t*) for different values of the fractional parameter α, y = 1.2, t = 0.5 and t = 1.0.

**Figure 2.** Profiles of surface ϕ(*<sup>x</sup>*, *y*, 1) for α = 0.3 and their sections with the planes *y* = 1.2 and *x* = 0.74, respectively.

**Figure 3.** Profiles of surface ϕ(*<sup>x</sup>*, *y*, 1) for α = 0.5 and their sections with the planes *y* = 1.2 and *x* = 0.74, respectively.

**Figure 4.** Profile of kernel *h*(<sup>α</sup>, *t*) for two values of the time t.

To draw these graphs we have used the following values of the parameters: *a* = 0.5, *b* = 1.0, *c* = 1.0, *d* = 1.5, *p* = 0.25, γ0 = 4, ε1 = ε2 = 0.2, *k*1 = 0.25, *k*2 = 0.5.

In Figure 1, we show the plotted curves ϕ(*<sup>x</sup>*, 1.2, 0.5) and ϕ(*<sup>x</sup>*, 1.2, 1.0) for different values of the fractional parameter α ∈ {0.3, 0.5, 0.7, 1.0}. It is observed from Figure 1 that the values of concentration ϕ decrease with the fractional parameter and have a peak (maximum/minimum) in the central area of the geometrical domain *D* = [*a*, *b*] × [*c*, *d*].

In Figures 2 and 3, we show the plotted surfaces ϕ(*<sup>x</sup>*, *y*, 1) for the fractional parameters α = 0.3 and α = 0.5, respectively, along with the sections through planes *y* = 1.2, *x* = 0.74.

For small values of the fractional parameter, the concentration ϕ has maximum values in the central area of the rectangular *D* = [*a*, *b*] × [*c*, *d*]. If the fractional parameter increases, the concentration has a minimum in the same area. This behavior is due to the kernel of the time-fractional Caputo–Fabrizio derivative (see Figure 4) which, for *t* = 0.5, increases for α ≤ 0.55 and decreases for α > 0.55. Therefore, at small values of the fractional parameter, the weight function in the fractional derivative has bigger values and the influence on the concentration is stronger.

## **5. Numerical Solution**

In this section, based on the finite difference approximation of the derivatives, we develop a numerical scheme to the problem formulated in Section 2 of the paper.

We consider Equation (8) along with the initial and boundary conditions (10) and (11), with (*x*, *y*, *t*) ∈ [*a*, *b*] × [*c*, *d*] × [0, *<sup>T</sup>*], *T* > 0. The discrete sets of nodes are defined as:

$$\begin{aligned} \mathbf{x}\_{i} &= a + jh\_{\mathbf{x}}, \ i = 0, 1, \ldots, N\_{1}, h\_{\mathbf{x}} = (b - a) / N\_{1}, \\ \mathbf{y}\_{j} &= c + jh\_{\mathbf{y}}, \ j = 0, 1, \ldots, N\_{2}, h\_{\mathbf{y}} = (d - c) / N\_{2}, \\ \mathbf{t}\_{n} &= nh\_{\mathbf{t}}, n = 0, 1, \ldots, N\_{\prime}, h\_{\mathbf{t}} = T / N. \end{aligned} \tag{55}$$

The numerical estimation of the solution ϕ(*<sup>x</sup>*, *y*, *t*) of Equation (8) in the point (*xi*, *yj*, *tn*) is denoted by ϕ*ni*,*j*= ϕ(*xi*, *yj*, *tn*). The initial and boundary conditions (10) and (11) are written as:

$$\begin{aligned} \varphi\_{i,j}^{0} &= \varphi(\mathbf{x}\_{i}, y\_{j}, t\_{0}) = 0, \; i = 0, 1, \ldots, N\_{1}, \; j = 0, 1, \ldots, N\_{2}, \\ \varphi\_{0,j}^{n} &= \varphi(a, y\_{j}, t\_{n}) = \mathcal{F}\_{1}(a, y\_{j}, t\_{n}) + \mathcal{G}\_{1}(a, y\_{j}) = \mathcal{P}\_{0,j'}^{n}, \; j = 0, 1, \ldots, N\_{2}, \; n = 1, 2, \ldots, N\_{\prime} \\ \varphi\_{N\_{1}, j}^{n} &= \varphi(b, y\_{j}, t\_{n}) = \mathcal{F}\_{1}(b, y\_{j}, t\_{n}) + \mathcal{G}\_{1}(b, y\_{j}) = \mathcal{Q}\_{N\_{1}, j'}^{n}, \; j = 0, 1, \ldots, N\_{2}, \; n = 1, 2, \ldots, N\_{\prime} \\ \varphi\_{i,0}^{n} &= \varphi(\mathbf{x}\_{i}, c, t\_{n}) = \mathcal{F}\_{1}(\mathbf{x}\_{i}, c, t\_{n}) + \mathcal{G}\_{1}(\mathbf{x}\_{i}, c) = \mathcal{R}\_{1,0'}^{n}, \; i = 0, 1, \ldots, N\_{1}, \; n = 1, 2, \ldots, N\_{\prime} \\ \varphi\_{i,N\_{2}}^{n} &= \varphi(\mathbf{x}\_{i}, d, t\_{n}) = \mathcal{F}\_{1}(\mathbf{x}\_{i}, d, t\_{n}) + \mathcal{G}\_{1}(\mathbf{x}\_{i}, d) = \mathcal{S}\_{1, N\_{1}}^{n}, \; i = 0, 1, \ldots, N\_{1}, \; n = 1, 2, \ldots, N\_{\prime} \end{aligned} \tag{56}$$

Using the following approximation for the first time-derivative:

$$\left. \frac{\partial \rho(\mathbf{x}\_i, y\_j, t)}{\partial t} \right|\_{t=t\_k} = \frac{1}{h\_l} [\rho(\mathbf{x}\_i, y\_j, t\_{k+1}) - \rho(\mathbf{x}\_i, y\_j, t\_k)] = \frac{1}{h\_l} [\rho\_{i,j}^{k+1} - \rho\_{i,j}^k], k = 0, 1, \dots, N - 1,\tag{57}$$

the time-fractional Caputo–Fabrizio derivative, in the point *t* = *tn*, can be approximated by:

$$\begin{split} \left. \text{CF}D\_{t}^{\alpha} \boldsymbol{q} (\mathbf{x}, \boldsymbol{y}, t) \right|\_{\left(\mathbf{x}\_{i}, \mathbf{y}\_{j}, t\right)} &= \frac{1}{1 - \alpha} \int\_{0}^{t\_{i}} \frac{\partial \boldsymbol{q} (\mathbf{x}\_{i}, \mathbf{y}\_{j}, t)}{\partial t} \bigg|\_{t = \pi} \exp \left( \frac{-a \left( t\_{i} - \tau \right)}{1 - \alpha} \right) d\tau = \\ \sum\_{k = 0}^{n - 1} \frac{1}{1 - \alpha} \int\_{t\_{k}}^{t\_{k}} \frac{\partial \boldsymbol{q} (\mathbf{x}\_{i}, \mathbf{y}\_{j}, t)}{\partial t} \bigg|\_{t = \pi} \exp \left( \frac{-a \left( t\_{n} - \tau \right)}{1 - \alpha} \right) d\tau = \\ \sum\_{k = 0}^{n - 1} \frac{1}{(1 - \alpha) b\_{i}} \Big( \boldsymbol{q}\_{i,j}^{k + 1} - \boldsymbol{q}\_{i,j}^{k} \bigg) \int\_{t\_{k}}^{t\_{k}} \exp \left( \frac{-a \left( t\_{n} - \tau \right)}{1 - \alpha} \right) d\tau = \sum\_{k = 0}^{n - 1} a\_{n,k} \Big( \boldsymbol{q}\_{i,j}^{k + 1} - \boldsymbol{q}\_{i,j}^{k} \bigg), n = 1, 2, \ldots \\ \end{split} \tag{58}$$

where

$$a\_{n,k} = \frac{1}{a\hbar} \left[ \exp\left(\frac{-a(t\_n - t\_{k+1})}{1 - a}\right) - \exp\left(\frac{-a(t\_n - t\_k)}{1 - a}\right) \right], a \in (0, 1), \\ n = 1, 2, \dots, N, \ k = 0, 1, \dots, (n - 1) \tag{59}$$

It is known [28,29] that Equation (58) approximates the time-fractional Caputo–Fabrizio derivative with an error of order *<sup>O</sup>*(*h*2*t*).

The space derivatives of the second order are approximated by:

$$\begin{aligned} \frac{\partial^2 \varphi(\mathbf{x}, y, t)}{\partial x^2} \Big|\_{ (\mathbf{x}\_i \mathbf{y}\_j, t\_n) } &= \frac{1}{h\_x^2} (\varphi\_{i+1, j}^n - 2\varphi\_{i,j}^n + \varphi\_{i-1,j}^n) \newline i = 1, 2, \dots, N\_1 - 1, j = 0, 1, \dots, N\_2, n = 1, 2, \dots, N, \\\frac{\partial^2 \varphi(\mathbf{x}, y, t)}{\partial y^2} \Big|\_{ (\mathbf{x}\_i \mathbf{y}\_j, t\_n) } &= \frac{1}{h\_y^2} (\varphi\_{i, j+1}^n - 2\varphi\_{i,j}^n + \varphi\_{i, j-1}^n) \newline i = 0, 1, \dots, N\_1, j = 1, 2, \dots, N\_2 - 1, n = 1, 2, \dots, N. \end{aligned} \tag{60}$$

Using the approximation formulae (58) and (60), Equation (8) is written as:

$$\begin{aligned} \sum\_{k=0}^{n-1} a\_{n,k} \{\boldsymbol{\varrho}\_{i,j}^{k+1} - \boldsymbol{\varrho}\_{i,j}^{k}\} + \beta \boldsymbol{\varrho}\_{i,j}^{n} &+ \beta\_1 \{\boldsymbol{\varrho}\_{i+1,j}^{n} - 2\boldsymbol{\varrho}\_{i,j}^{n} + \boldsymbol{\varrho}\_{i-1,j}^{n}\} + \\ \beta\_2 \{\boldsymbol{\varrho}\_{i,j+1}^{n} - 2\boldsymbol{\varrho}\_{i,j}^{n} + \boldsymbol{\varrho}\_{i,j-1}^{n}\} &= H\_{i,j'}^{n}, i = 1,2,\dots,N\_1 - 1, \dots \; j = 1,2,\dots,N\_2 - 1, \; n = 1,2,\dots,N\_\prime \end{aligned} \tag{61}$$

where

$$\beta\_1 = \frac{-\varepsilon\_1}{h\_x^2}, \beta\_2 = \frac{-\varepsilon\_2}{h\_y^2}, H\_{i,j}^{\text{tr}} = F(\mathbf{x}\_{i\prime} y\_{j\prime}, t\_{\text{tr}}) + \text{G}(\mathbf{x}\_{i\prime} y\_{j\prime}) \tag{62}$$

For *n* = 1, Equation (61) becomes:

$$d\_{1,0}q^1\_{\dot{i},\dot{j}} + \beta\_1(q^1\_{1+1,\dot{j}} + q^1\_{\dot{i}-1,\dot{j}}) + \beta\_2(q^1\_{\dot{i},\dot{j}+1} + q^1\_{\dot{i},\dot{j}-1}) = H^1\_{\dot{i},\dot{j}'} \tag{63}$$

where *d*1,0 = *<sup>a</sup>*1,0 + β − <sup>2</sup>(β1 + β2).

> Making *j* = 1 and *i* = 1, 2, ... , *N*1 − 1 into Equation (63), we obtain the following algebraic system:

$$\mathbf{M}\mathbf{X}\_1 + \mathbf{M}\_0 \mathbf{X}\_2 = \mathbf{D}\_{1'}^1 \tag{64}$$

where *M* is a square matrix of order (*<sup>N</sup>*1 − 1) whose elements are given by:

$$\begin{aligned} M\_{0,k} &= d\_{1,0}\delta\_{0,k} + \beta\_1 \delta\_{1,k}, k = 0, 1, \dots, N\_1 - 2, \\ M\_{i,k} &= \beta\_1 \delta\_{i,k+1} + d\_{1,0}\delta\_{i+1,k+1} + \beta\_1 \delta\_{i+2,k+1}, i = 1, 2, \dots, N\_1 - 2, k = k = 0, 1, \dots, N\_1 - 2, \end{aligned} \tag{65}$$

so,

$$\mathbf{M} = \begin{pmatrix} d\_{1,0} & \beta\_1 & 0 & 0 & 0 & \dots & 0 & 0 & 0\\ \beta\_1 & d\_{1,0} & \beta\_1 & 0 & 0 & \dots & 0 & 0 & 0\\ 0 & \beta\_1 & d\_{1,0} & \beta\_1 & 0 & \dots & 0 & 0 & 0\\ \dots & \dots & \dots & \dots & \dots & \dots & \dots & \dots & \dots & \dots\\ 0 & 0 & 0 & 0 & 0 & \dots & \beta\_1 & d\_{1,0} & \beta\_1\\ 0 & 0 & 0 & 0 & 0 & 0 & \dots & \beta\_1 & d\_{1,0} \end{pmatrix} \tag{66}$$

the matrix *M*2 is the diagonal square matrix of order (*<sup>N</sup>*1 − 1) with elements *<sup>M</sup>*1*i*,*<sup>j</sup>* = β2δ*i*,*j*, *i*, *j* = 0, 1, ... , *N*1 − 2, <sup>δ</sup>*i*,*<sup>j</sup>* = 0, *i* - *j*, 1, *i* = *j*, is the Kronecker delta and *X*1, *X*2 *D*11 are the column-matrices.

$$\mathbf{X}\_{1} = \begin{pmatrix} \varrho\_{1,1}^{1} \\ \varrho\_{2,1}^{1} \\ \cdots \\ \varrho\_{N\_{1}-1,1}^{1} \end{pmatrix}, \mathbf{X}\_{2} = \begin{pmatrix} \varrho\_{1,2}^{1} \\ \varrho\_{2,2}^{1} \\ \cdots \\ \varrho\_{N\_{1}-1,2}^{1} \end{pmatrix}, \mathbf{D}\_{1}^{1} = \begin{pmatrix} D\_{11}^{1} \\ D\_{12}^{1} \\ \cdots \\ D\_{1N\_{1}-1}^{1} \end{pmatrix}, \quad D\_{11}^{1} = H\_{1,1}^{1} - \beta\_{1} \mathbf{P}\_{0,1}^{1} - \beta\_{2} \mathbf{R}\_{1,0}^{1},$$

$$D\_{11}^{1} = H\_{1,1}^{1} - \beta\_{2} \mathbf{R}\_{1,0}^{1},$$

$$D\_{1N\_{1}-1}^{1} = H\_{1,1}^{1} - \beta\_{1} \mathbf{D}\_{1,0}^{1},$$

The algebraic system (64) has (*<sup>N</sup>*1 − 1) equations with <sup>2</sup>(*<sup>N</sup>*1 − 1) unknown functions.

Making *j* = 2, 3, ... , *N*2 − 2 and *i* = 1, 2, ... , *N*1 − 1 into Equation (63), we obtain the algebraic system:

$$\mathbf{M}\_0 \mathbf{X}\_{j-1} + \mathbf{M} \mathbf{X}\_j + \mathbf{M}\_0 \mathbf{X}\_{j+1} = D^1\_{j},\tag{68}$$

where

$$\begin{array}{l} D^1\_{j1} = H^1\_{1,j} - \beta\_1 P^1\_{0,j'} \\ D^1\_{\vec{u}} = H^1\_{i,j'} \; i = 2, 3, \dots, N\_1 - 2, \\ D^1\_{jN\_1 - 1} = H^1\_{N\_1 - 1, j} - \beta\_1 Q^1\_{N\_1, j} \end{array} \tag{69}$$

Making *j* = *N*2 − 1 and *i* = 1, 2, ... , *N*1 − 1 into Equation (63), we obtain the algebraic system:

$$\mathbf{M}\_0 \mathbf{X}\_{\text{N}\_2-2} + \mathbf{M} \mathbf{X}\_{\text{N}\_2-1} = \mathbf{D}\_{\text{N}\_2-1}^1 \tag{70}$$

where

$$\begin{array}{l}D^1\_{N\_2-11} = H^1\_{1,N\_2-1} - \beta\_1 P^1\_{0,N\_2-1} - \beta\_2 S^1\_{1,N\_2} \\ D^1\_{N\_2-1i} = H^1\_{i,N\_2-1} - \beta\_2 S^1\_{i,N\_2} \; i = 2,3,\ldots,N\_1 - 2, \\ D^1\_{N\_2-1N\_1-1} = H^1\_{N\_1-1,N\_2-1} - \beta\_1 Q^1\_{N\_1,N\_2-1} - \beta\_2 S^1\_{N\_1-1,N\_2} \end{array} \tag{71}$$

The algebraic systems (64), (68) and (70) can be written as a unique algebraic system of (*<sup>N</sup>*1 − <sup>1</sup>)(*<sup>N</sup>*2 −<sup>1</sup>) equations with (*<sup>N</sup>*1 −<sup>1</sup>)(*<sup>N</sup>*<sup>2</sup> −<sup>1</sup>) unknown functionsϕ<sup>1</sup>*i*,*j*, *i* = 1, 2, ... , *N*1 −1, *j* = 1, 2, ... , *N*2 − 1, namely

$$
\begin{pmatrix}
\mathcal{M} & \mathcal{M}\_0 & 0 & 0 & 0 & \dots & 0 & 0 & 0 \\
\mathcal{M}\_0 & \mathcal{M} & \mathcal{M}\_0 & 0 & 0 & \dots & 0 & 0 & 0 \\
0 & \mathcal{M}\_0 & \mathcal{M} & \mathcal{M}\_0 & 0 & \dots & 0 & 0 & 0 \\
\dots & \dots & \dots & \dots & \dots & \dots & \dots & \dots & \dots \\
0 & 0 & 0 & 0 & 0 & \dots & \mathcal{M}\_0 & \mathcal{M} & \mathcal{M}\_0 \\
0 & 0 & 0 & 0 & 0 & 0 & \dots & \mathcal{M}\_0 & \mathcal{M}
\end{pmatrix}
\begin{pmatrix}
X\_1 \\ X\_2 \\ X\_3 \\ \dots \\ \dots \\ X\_{N\_2 - 2} \\ X\_{N\_2 - 1}
\end{pmatrix} = \begin{pmatrix}
D\_1^1 \\ D\_2^1 \\ D\_3^1 \\ \dots \\ \dots \\ D\_{N\_2 - 1}^1
\end{pmatrix} \tag{72}
$$

The solution of system (72) gives the values of function ϕ(*<sup>x</sup>*, *y*,*<sup>t</sup>*) in points (*xi*, *yj*,*t*1) for *i* = 1, 2, ... , *N*1 − 1 , *j* = 1, 2, ... , *N*2 − 1.

The numerical procedure is continued in a similar way for *n* = 2, 3, ... , *N* when Equation (61) is written in the equivalent form:

$$d\_{\mathbf{u},\mathbf{u}-1}q\_{i,j}^{\mathbf{u}} + \beta\mathbf{i}\left(q\_{i+1,j}^{\mathbf{u}} + q\_{i-1,j}^{\mathbf{u}}\right) + \beta\mathbf{2}\left(q\_{i,j+1}^{\mathbf{u}} + q\_{i,j-1}^{\mathbf{u}}\right) = \mathbf{G}\_{i,j}^{\mathbf{u}}, i = 1,2,\dots,N\_1 - 1, \ j = 1,2,\dots,N\_2 - 1,\tag{73}$$

where

$$\begin{aligned} d\_{n,n-1} &= a\_{n,n-1} + \beta - 2(\beta\_1 + \beta\_2), n = 2, 3, \dots, N\_\prime \\ G\_{i,j}^n &= H\_{i,j}^n + a\_{n,n-1} \varphi\_{i,j}^{n-1} - \sum\_{k=0}^{n-2} a\_{n,k} (\varphi\_{i,j}^{k+1} - \varphi\_{i,j}^k), i = 1, 2, \dots, N\_1 - 1, j = 1, 2, \dots, N\_2 - 1. \end{aligned} \tag{74}$$

Finally, the proposed numerical procedure gives the values of the solution ϕ(*<sup>x</sup>*, *y*, *t*) in grid points (*xi*, *yj*, *tn*), *i* = 1, 2, ... , *N*1 − 1, *j* = 1, 2, ... , *N*2 − 1, *n* = 1, 2, ... , *N*.

To compare the numerical values obtained with the analytical solution, respectively, by the numerical method, we have considered the following particular problem:

$$\begin{aligned} F(\mathbf{x}, y, t) &= 5000 \exp\left(-\frac{k\_1 x}{2x\_1} - \frac{k\_2 y}{2x\_2}\right), G(\mathbf{x}, y) = 0, \ (\mathbf{x}, y, t) \in [0, 0.05] \times [0, 0.05] \times [0, 0.1], \\ F\_1(\mathbf{x}, y, t) &= 0, G\_1(\mathbf{x}, y) = 0, \rho(\mathbf{x}, y, 0) = 0, \\ k\_1 = k\_2 = 0.5, \ \varepsilon\_1 = \varepsilon\_2 = 2, \end{aligned} \tag{75}$$

In this case, the double Fourier transform of the function *<sup>F</sup>*(*<sup>x</sup>*, *y*, *t*) is:

$$\begin{split} \widetilde{F}(n,k) &= \underbrace{\frac{n\pi(b-a)\left[\exp\left(\frac{-ak\_1}{2\tau\_1}\right) - (-1)^n \exp\left(\frac{-kb\_1}{2\tau\_1}\right)\right]}\_{\begin{pmatrix} \left(\frac{k\_1}{2\tau\_1}\right)^2 (b-a)^2 + (n\pi)^2\\ \left(\frac{k\_2}{2\tau\_2}\right)^2 (b-a)^2 + (n\pi)^2 \end{pmatrix}} \times \\ &\underbrace{\frac{k\pi(b-a)\left[\exp\left(\frac{-ck\_2}{2\tau\_2}\right) - (-1)^k \exp\left(\frac{-dk\_2}{2\tau\_2}\right)\right]}\_{\begin{pmatrix} \frac{k\_2}{2\tau\_2} \end{pmatrix}^2 (d-c)^2 + (k\pi)^2}} \end{split} \tag{76}$$

and the analytical solution is given by:

$$\begin{split} \varphi(\mathbf{x}, \mathbf{y}, t) &= \frac{4a}{(b-a)(d-c)} \sum\_{n=1}^{\infty} \sum\_{k=1}^{\infty} \frac{\widetilde{F}(n\boldsymbol{k})}{b\_{nk}} \sin\left[\frac{n(\mathbf{x} - \boldsymbol{a})\pi}{b - a}\right] \sin\left[\frac{k(\mathbf{y} - \boldsymbol{c})\pi}{d - c}\right] + \\ &\frac{4}{(b-a)(d-c)} \sum\_{n=1}^{\infty} \sum\_{k=1}^{\infty} \frac{\widetilde{F}(n\boldsymbol{k})[(1-a)b\_{nk} - a\boldsymbol{a}\_{nk}]}{a\_{nk}b\_{nk}} \sin\left[\frac{n(\mathbf{x} - \boldsymbol{a})\pi}{b - a}\right] \sin\left[\frac{k(\mathbf{y} - \boldsymbol{c})\pi}{d - c}\right] \mathbf{e}^{-\frac{k\_{\mathrm{nl}}t}{d}} .\end{split} \tag{77}$$

The values of function ϕ(*<sup>x</sup>*, *y*, *t*) at the instant *t* = 0.01, obtained with the numerical method and with expression (77) given in the Table 1, are in good agreement.


**Table 1.** Values of function ϕ(*<sup>x</sup>*, *y*, 0.01).

#### **6. Discussion and Conclusions**

A two-dimensional advection–diffusion process with short-tail memory and concentrated source situated in the symmetry center of the geometrical domain has been investigated.

The studied mathematical model is described by a fractional differential equation with a time-fractional Caputo–Fabrizio derivative with exponential kernel.

Analytical solutions of the proposed initial-boundary value problem were obtained using the integral transform method (Laplace and finite Fourier transforms are employed).

Numerical solutions of the studied problem have been determined using finite difference approximations.

Numerical simulations for both analytical and numerical solutions have been carried out using the software Mathcad.

The fractional order of the Caputo–Fabrizio derivative has a significant influence on the pollutant concentration. For small values of the fractional parameter, the concentration has maximum values in the central area of the geometrical domain. If the fractional parameter increases, the concentration has a minimum in the same area. This behavior is due to the memory kernel of the time-fractional Caputo–Fabrizio derivative.

**Author Contributions:** The authors have equally contributed to this paper.

**Funding:** This research received no external funding.

**Conflicts of Interest:** On behalf of all authors, the corresponding author states that there is no conflict of interest.
