Advection-diffusion-reaction equations (ADRE) have been widely used as mathematical models in many areas of science and engineering, for example, in areas such as water pollution, air pollution, molecular diffusion, and chemical engineering. There are now many numerical methods that have been developed to solve linear and nonlinear ADREs. Some examples are as follows. In 2012, Diego et al [1] developed a model consisting of a decoupled system of advection-diffusion-reaction equations, along with the Navier-Stokes equations of incompressible flow, and solved the model by using the finite element method. In 2012, Savovic and Djordjevich [2] proposed an explicit finite difference method to solve a one-dimensional advection-diffusion equation with variable coefficients in semi-infinite media for three dispersion problems. In 2013, Appadu and Gidey [3] proposed two time-splitting procedures that they used to solve a 2D advection-diffusion equation with constant coefficients. In 2013, Bause and Schwegler [4] proposed a higher order finite element approximation for systems of coupled convection-dominated transport equations. In 2015, Mojtabi and Deville [5] proposed a finite element method to solve analytically using separation of variables and numerical solution of a time-dependent one-dimensional linear advection-diffusion equation with Dirichlet homogeneous boundary conditions and an initial sine function. In 2016, Gharehbaghi [6] proposed an explicit and an implicit differential quadrature method to compute a numerical solution of a one-dimensional time-dependent advection-diffusion equation with variable coefficients in a semi-infinite domain. In 2017, Gyrya and Lipnikov [7] proposed an arbitrary order mimetic finite difference discretization for the diffusion equation with non-symmetric tensorial diffusion coefficient in a mixed formulation on general polygonal meshes. In 2017, Bahar and Gurarslan [8] studied the effects of Lie-Trotter and Strang splitting methods on the solution of a one-dimensional advection-diffusion equation. In 2017, Oyjindal and Pochai [9] proposed models for numerical simulations of air pollution measurements near an industrial zone. The numerical experiments consisted of different conditions including normal and controlled emissions. The models were then solved using an explicit finite difference technique and the solutions were compared with the measurements for the controlled and uncontrolled point sources at the monitoring points. In 2018, Suebyat and Pochai [10] developed a three-dimensional air pollution measurement model for a heavy traffic area under a Bangkok sky train platform and used finite difference techniques to approximate the solutions. In 2018, Al-Jawary et al. [11] proposed a semi-analytical technique for finding the exact solutions for different types of Burger’s equations and systems of equations in 1D, 2D, and 3D. Also, the method was applied to solve the diffusion and advection-diffusion equations. In 2018, Kusuma et al. [12] proposed a finite difference (FTCS) method to solve a model of pollution distribution in a street tunnel using two dimensional advection and three dimensional diffusion in a rectangular box domain. In 2018, Lou et al. [13] developed reconstructed Discontinuous Galerkin methods for solving linear advection-diffusion equations on hybrid unstructured grids based on a first-order hyperbolic system formulation. In 2018, Bhatt et al, [14] developed a Krylov subspace approximation-based locally extrapolated exponential time differencing method and studied its accuracy and efficiency for solving three-dimensional nonlinear advection-diffusion-reaction systems. In 2020, Pananu et al. [15] analyzed the convergence of the finite difference method with the implicit forward time central space (FTCS) scheme for the two-dimensional advection-diffusion-reaction equation (ADRE) and applied the scheme to a pollutant dispersion with removal mechanism model in a reservoir. In 2020, Heng and Guodong [16] improved the element-free Galerkin method and used it for solving 3D advection-diffusion problems. In 2020, Cruz-Quintero and Jurado [17] proposed a backstepping design for the boundary control of a reaction-advection-diffusion equation with constant coefficients and Neumann boundary conditions. In 2021, Para et al. [18] proposed the characteristic finite volume method for solving a convection-diffusion problem on two-dimensional triangular grids and compared the accuracy of four piecewise linear reconstruction techniques on structured triangular grids, namely, Frink, Holmes-Connell, Green-Gauss, and least squares methods. In 2021, Hong et al. [19] discussed the numerical solution of 3D unsteady advection-diffusion equations using a meshless numerical scheme with a space-time backward substitution method. In 2021, Irfan and Hidayat [20] proposed a meshless finite difference method with B-splines for the numerical solution of coupled advection-diffusion-reaction problems. In 2021, Shahid et al [21] studied an epidemic type model with advection and diffusion terms for the transmission dynamics of a computer virus model with fixed population density. In 2021, García and Jurado [22] designed an adaptive boundary control for a parabolic type reaction-advection-diffusion PDE under the assumption of unknown parameters for both advection and reaction terms and Robin and Neumann boundary conditions. In 2022, Para et al. [23] proposed an explicit characteristic-based finite volume method (FVM) for the numerical solution of advection-diffusion-reaction equations (ADRE) and applied it to solve some 1D and 2D water pollution problems which can be modeled in terms of ADREs. Then, the FVM results were compared with numerical results obtained using a finite difference method (FDM) with an implicit forward time central space (FTCS) scheme.

In this article, we consider the application of a finite difference method based on the implicit forward time central space (FTCS) scheme to the numerical solution of an advection-diffusion-reaction equation of the following form: where $\varphi$ is a scalar quantity, $\mathbf{v}=\mathbf{v}\left(\mathbf{x}\right)$ is a given advection velocity vector, $\epsilon \ge 0$ is a diffusion coefficient, $\kappa$ is a reaction coefficient, $q=q(\mathbf{x},t)$ is a prescribed source term, $\mathbf{x}$ is a position vector and time $t\in [0,T]$, where $T>0$. The present paper is organized as follows. The convergence theory is briefly provided in Section 2. The consistency and stability analysis and numerical results of the method are given in Section 3. Using the method, numerical results and their graphs for the air pollution problems written in terms of the ADRE in Eq. (1) are shown in Section 4. The conclusion of our work is presented in Section 5.

We consider the following initial time, space boundary value problem (IBVP) for the advection-diffusion-reaction equation where $(\mathbf{x},t)\in \mathsf{\Omega }\times (0,T]$, and with the initial time condition and either Dirichlet or Neumann boundary space conditions In Eqs. (2)–(4), T is a finite time, $\partial \mathsf{\Omega }$ is the boundary of the spatial domain $\mathsf{\Omega }$, ${\xi }_D$ and ${\xi }_N$ are given functions, $\mathbf{v}\left(\mathbf{x}\right)$ is an advection vector, $\epsilon >0$ is a diffusion coefficient, $\kappa$ is a reaction constant and $q(\mathbf{x},\mathbf{t})$ is a source term. In addition, L is a linear partial differential operator from a space of continuous functions F to a space of continuous functions H, where the dependent variable $\varphi (\mathbf{x},t)\in F$, and where $L\varphi (\mathbf{x},t)\in H$ and $q(\mathbf{x},t)\in H$.

The first step in developing the finite difference method is to discretize $\mathsf{\Omega }$, $\partial \mathsf{\Omega }$ and $[0,T]$. We assume that there are $N_1$ step sizes of length $\Delta x=\frac{L_1}{N_1}$ in the x direction, $N_2$ step sizes of length $\Delta y=\frac{L_2}{N_2}$ in the y direction and $N_3$ step sizes of length $\Delta z=\frac{L_3}{N_3}$ in the z direction. For simplicity, we will assume that $\Delta x=\Delta y=\Delta z=h$ and therefore $L_1=N_{1}h$, $L_2=N_{2}h$ and $L_3=N_{3}h$. Then the discretized version of $\mathsf{\Omega }$ is the spatial grid where $J_1=\{0,1,2,...,N_1\}$, $J_2=\{0,1,2,...,N_2\}$, $J_3=\{0,1,2,...,N_3\}$. The boundary of the spatial grid is

Similarly, the interval $[0,T]$ is discretized with N step sizes of size $\Delta t=\frac{T}{N}$. The set of discretized times is then

Finally, we define $F_h$ and $H_h$ as discretized function spaces with domain ${\mathsf{\Omega }}_h\times T_K$. We also define a projection ${\mathcal{P}}_h\left(\varphi \right)\in F_h$ of the exact solution $\varphi (\mathbf{x},t)$ of (2) onto ${\mathsf{\Omega }}_h\times T_K$ and a projection ${\mathcal{P}}_h\left(q\right)\in H_h$ of the source term $q(\mathbf{x},t)$ onto ${\mathsf{\Omega }}_h\times T_K$.

We can then write a finite difference scheme (FDS) corresponding to the continuous partial differential equation (2) in the form where $L_h:F_h\to H_h$ is a finite difference operator, corresponding to the partial differential operator L, which acts from the discrete function space $F_h$ to the discrete function space $H_h$. The initial condition in (3) can then be rewritten as the discretized initial condition and, for the Dirichlet case, the boundary conditions in (4) can be rewritten as the discretized boundary conditions

Consider the following initial-boundary value problem in which the above partial differential equation is the component form of (1). Here $\varphi =\varphi (x,y,z,t)$ is the dependent variable and $\mathbf{v}=(u,v,w)$ is the advection velocity where $u,v$ and w are the velocities in the x-, y- and z-directions, respectively.

In this section, we will analyze the convergence of the finite difference method with the implicit FTCS scheme applied to (20). We discretize the domain of the problem: $(x_i,y_j,z_k,t^n)\in {\mathsf{\Omega }}_h\times T_K$, where the spatial grid ${\mathsf{\Omega }}_h$ is defined in (5) and the time grid $T_K=[0,T]$ has points $t^n=n\Delta t$ with $n\in K=\{0,1,2,...,N\}.$

Following the method in Section 2, we use Lax’s equivalence theorem to prove the convergence of the FTCS by proving its consistency and stability. We first show the consistency of the FTCS. The implicit FTCS scheme of Eq. (20) is as follows. where ${\varphi }_{i,j,k}^n={\varphi }^h(x_i,y_j,z_k,t^n)$ and $q_{i,j,k}^n=q(x_i,y_j,z_k,t^n)$.

The truncation error for (22) obtained using the method is Therefore, from (24), we have Taking $(\Delta t,\Delta x,\Delta y,\Delta z)\to 0$, then we have

We will now prove the stability of the finite difference method with the implicit FTCS scheme for numerical solution of (20). In order to illustrate stability of the method, we initially show that condition (19) holds for the homogeneous equation of (22) and then that condition (17) is satisfied for the nonhomogeneous equation (22). For convenience, we set $\tau =\Delta t,h=\Delta x=\Delta y=\Delta z$.

Setting $q_{i,j,k}^n=0$ in equation (22), we obtain the discretized homogeneous version of (22). Assuming a discretized solution of the resulting equation as ${\varphi }_{i,j,k}^n={\lambda }^n{e}^{I(\alpha i+\beta j+\eta k)}$, $I=\sqrt{-1}$[24], we obtain the following equation for the eigenvalue $\lambda$ Then, we can solve (25) for the eigenvalue $\lambda$ as follows. Letting $\gamma =\frac{4\epsilon \tau }{h^2}$ and $\mu =\frac{\tau }{h}$, we then obtain, Consequently, the magnitude of $\lambda$ is Therefore, condition (19) for Von Neumann stability, which is independent of $\tau$ and h, has been proved.

We now prove that condition (17) is satisfied for the nonhomogeneous equation (22) with $q_{i,j,k}^n\ne 0$ by finding a solution ${\varphi }_{i,j,k}^{n+1}$ for the linear system of difference equations (22) and (23) if ${\varphi }_{i,j,k}^n$ is assumed known. Substituting $\Delta x=\Delta y=\Delta z=h$ and $\Delta t=\tau$ in (22) and (23), we obtain After rearranging, we can rewrite (26) in the form and where the step sizes $\tau$ and h can be chosen so that the coefficients $a,b,c,d,e,f,g$ satisfy the conditions

Let $0<p<N_1,0<q<N_2,0<r<N_3$. Then

This completes the proof.

Now we will prove the stability of the implicit FTCS scheme by showing that inequality (17) is satisfied.

It is obvious that if $(i^{*},j^{*},k^{*})\in \overline{J},$ then we have

If $(i^{*},j^{*},k^{*})\in \mathcal{J}=J-\overline{J},$ then replacing $(i,j,k)$ in (26) with $(i^{*},j^{*},k^{*})$ we obtain

Next we make the following hypothesis:

Without loss of generality we assume that ${\varphi }_{i^{*},j^{*},k^{*}}^{n+1}>0$ and further assume that the hypothesis (H1) holds. We can estimate the right side of (34) as