Numerical Solution of Two-Dimensional Nonlinear Unsteady Advection-Diffusion-Reaction Equations with Variable Coefficients

. Te advection-difusion-reaction (ADR) equation is a fundamental mathematical model used to describe various processes in many diferent areas of science and engineering. Due to wide applicability of the ADR equation, fnding accurate solution is very important to better understand a physical phenomenon represented by the equation. In this study, a numerical scheme for solving two-dimensional unsteady ADR equations with spatially varying velocity and difusion coefcients is presented. Te equations include nonlinear reaction terms. To discretize the ADR equations, the Crank–Nicolson fnite diference method is employed with a uniform grid. Te resulting nonlinear system of equations is solved using Newton’s method. At each iteration of Newton’s method, the Gauss–Seidel iterative method with sparse matrix computation is utilized to solve the block tridiagonal system and obtain the error correction vector. Te consistency and stability of the numerical scheme are investigated. MATLAB codes are developed to implement this combined numerical approach. Te validation of the scheme is verifed by solving a two-dimensional advection-difusion equation without reaction term. Numerical tests are provided to show the good performances of the proposed numerical scheme in simulation of ADR problems. Te numerical scheme gives accurate results. Te obtained numerical solutions are presented graphically. Te result of this study may provide insights to apply numerical methods in solving comprehensive models of physical phenomena that capture the underlying situations.


Introduction
Te ADR equation is a partial diferential equation used to represent various mathematical models in science and engineering.Te equation describes transport processes driven by advection, difusion, and reaction.Advection involves the collective displacement of species propelled by a fow feld.Difusion, characterized by the movement of species driven by concentration gradients, serves to homogenize concentration distributions over time.Reaction process occurring within the system may either generate or deplete the transported species.In practical scenarios, the coupling of advection, difusion, and reaction processes results in intricate interactions, leading to temporal and spatial alterations within the ADR system.Te ADR equation is used in modeling contaminant transport [1,2], fuid fow [3], mass and heat transfer [4], chemical reaction process [5], difusion across a cell membrane [6], tumor cell expansion [7], population dynamics [8], and so on.
Because of its importance in many felds, obtaining accurate solution for ADR equation has been the interest of numerous researchers.Tere have been a large number of solutions for ADR equations presented in literature using analytic [1,9,10] and numerical [11][12][13][14][15] techniques.Due to the intricate nature of many advection-difusion-reaction problems, numerical techniques are commonly employed to obtain solutions for the ADR equation.Various numerical techniques have been employed in the literature to compute numerical solutions for two-dimensional nonlinear ADR equations.For instance, Mesgarani et al. [16] used radial basis functions to solve time-dependent nonlinear ADR equations with variable coefcients.Tey transformed the ADR equation into system ordinary diferential equations and used the fourth-order Runge-Kutta method to compute solution of the system.Similarly, Ali et al. [17] applied a fnite diference method incorporating with Lucas and Fibonacci polynomials to obtain the solution of one-and two-dimensional nonlinear ADR equations.Ngondiep [18] employed a six-level time-split Leap-Frog/Crank-Nicolson approach to approximate solutions of two-dimensional nonlinear time-dependent ADR equations.He deduced that the numerical scheme is fast, spatial fourth-order convergent, and temporal second-order accurate.Zhang et al. [19] used fourth-order fully implicit fnite diference scheme to solve unsteady nonlinear ADR equations.Tambue [20] presented a fnite volume method combined with exponential time diferencing to solve the ADR equation involving the nonlinear reaction term.
Numerous researchers have solved ADR equations that model real-realistic problems.Yu et al. [21] utilized the homotopy analysis method (HAM) to simulate a contaminant transport model governed by ADR equations, presenting contaminant concentration profles at various time points based on simulation outcomes.Para et al. [12] investigated a water pollution scenario represented by ADR equations, employing a fnite volume method.Teir numerical solution was compared with results obtained using a fnite diference method, revealing satisfactory agreement between the solutions.Gurarslan et al. [22] obtained a numerical solution for contaminant transport described by a one-dimensional advection-difusion equation, employing a fnite diference scheme in space and a Runge-Kutta method in time.Te integrated scheme yielded accurate solutions.Putri et al. [23] conducted numerical simulations of advection-difusion equations representing oxygen demand concentration in waste stabilization ponds using a fnite diference method.Pananu et al. [24] utilized a fnite diference method scheme to investigate a water pollutant dispersion problem.
A Crank-Nicolson method is an implicit fnite diference method for solving partial diferential equations involving time and space derivatives [25].In solving partial diferential equations using the Crank-Nicolson method, both the time and space derivatives are approximated by central diferences.Te method is unconditionally stable and secondorder accurate both in space and time variables for solving the difusion equation.However, the method requires more computing time to solve the resulting system of algebraic equations.Specially, when the method is applied to solve nonlinear and multidimensional problems, linearization techniques like Newton's method [26] are required.Solving the system of nonlinear equations introduced through discretization demands much extra time.Any method solving system of linear equations (matrix-inversion method, Gauss-Seidel method, Lu decomposition method, and so on) can be used to solve the block tridiagonal system to obtain the error correction vector required in Newton's method.
Te objective of this paper is to devise a numerical scheme based on fnite diference approximation to solve ADR equations with all of the following properties: (i) unsteady, (ii) two-dimensional, (iii) nonlinear reaction term, and (iv) spatially varying velocity and difusion coefcients.Te intention is to employ the developed numerical scheme to simulate the distribution of concentration of a contaminant within a fowing fuid.

Problem Description
In this study, the ADR equation of the following form is considered [27]: with initial condition and boundary condition where

Discretization of ADR Equation
To apply fnite diference discretization of (4), we divide the spatial domain into uniform meshes of step sizes ∆x and ∆y in x-and y-directions, as shown in Figure 1, and the time domain with step size ∆t.Te grid points are defned as ( In equation ( 4), approximating the time derivative by central diference at point (i, j, k + (1/2)), the derivatives of the difusion coefcients by central diference, the derivatives of C by the averages of central diference at k th and (k + 1) th time levels, and reaction term using the average values of R at k th and (k + 1) th time levels, we get Note that in discretization ( 6), the Crank-Nicolson Method is used.Rearranging equation ( 6), we obtain where International Journal of Mathematics and Mathematical Sciences For j � 2 and i � 2, 3, . . ., N x − 2, j � 3 and i � 2, 3, . . ., N x − 2, . . ., j � N y − 2 and i � 2, 3, . . ., N x − 2, equation (7) yields a system of nonlinear equations with

Newton's Method
For p � 2, 3, . . ., (N x − 2) × (N y − 2), defne equations where the functions f p include variable terms from (k + 1) th time level values and reaction terms and constant terms from the boundary values and ) T be a vector function whose components are the functions in (9).
T be the vector of unknowns in (7).Let X 0 be the initial approximation of X.Using Newton's method, an error correction vector h � (h 1 , h 2 , h 3 , . . ., h (N x − 2)×(N x − 2) ) T such that X � X 0 + h can be obtained iteratively by solving the transformed linear system: where A i,j � (zf i (X 0 )/zX j ).Te iterative process for solving h in (10) at each (k + 1) th time level continues until we obtain X at this time level with the required accuracy.Te computation is performed until f(X) < Tol, where Tol is the tolerance of approximation for Newton's method.Te matrix A in ( 10) is a block tridiagonal matrix.For instance, for sample mesh indicated in Figure 1 with 6 × 6 grid points, A has the following form: To avoid computation with zero entries in matrix A, the Gauss-Seidel iterative method with sparse matrix computation is used to solve the system in (10).A termination criterion ‖h iter+1 − h iter ‖ ∞ < tol is used in the Gauss-Seidel iterative method, where tol is a prescribed number to bound the error.Te values of X at k + 2, k + 3, . .., N t time levels are solved similarly.

Convergence and Accuracy of the Numerical Scheme
Let L and L D be the diferential operator and the fnite diference operator representing equations ( 4) and ( 6 Te truncation error vanishes as (Δt, Δx, Δy) ⟶ (0, 0, 0).Hence, the fnite diference scheme is consistent with the partial diferential equation.
Let us analyse the stability of the fnite diference scheme (6) using the Von Neumann analysis method for some cases [25].Assume velocity and difusion coefcients are constants and the reaction term is zero.Suppose the discretization of C in (6) at grid point (i, j, k) has the form [13] where , A, β 1 , and β 2 are a constants, and ξ is the amplifcation factor.Substituting (12) in (6) and dividing both sides of equations by Aξ k e Iβ 1 iΔx e Iβ 2 jΔx , we get International Journal of Mathematics and Mathematical Sciences Using the facts that e Iθ � cos θ + Isinθ and e − Iθ � cos θ − Isinθ and solving for ξ from (13), we have where For the stability, we must have |ξ| ≤ 1.From ( 14), Hence, the fnite diference scheme ( 6) is unconditionally stable.According to Lax's equivalence theorem, the fnite diference scheme (6) gives a convergent solution for the considered cases.
In the proposed numerical scheme, the linearized system of equations at each iteration is solved using the Gauss-Seidel Iterative method to obtain error correction vector.Tus, if A is diagonally dominant, the Gauss-Seidel iterative method is convergent.A is diagonally dominant if the following inequality is satisfed for each row in the matrix: where X 0 is the initial guess to be taken for (k + 1) th iteration and m indicates the m th component of X 0 .In (16), the following indexes have to be used: 2 ≤ i ≤ N x − 1 for all terms in the inequality; 2 ≤ j ≤ N y − 1 for the term in the left side and for the frst and second terms in the right side; 2 ≤ j ≤ N y − 2 for the third term; and 3 ≤ j ≤ N y − 1 for the fourth term in the right side.Te accuracy of the numerical scheme is calculated by the absolute maximum error formula given by [ where C k i,j and C(x i , y j , t k ) are the approximate and exact solutions, respectively.Te rate of convergence for the scheme is calculated using the formula as follows: Here, E 1 and E 2 are maximum absolute errors with sets of number of grid points S 1 � N x , Ny, N t   and S 2 � 2N x ,  2Ny, 2N t }, respectively, where the set S 2 is obtained using half of the spatial and temporal step sizes used in S 1 .If the exact solution is not known in (17), the absolute maximum error can be estimated using two approximate solutions as and the rate of convergence can be computed using ( 19) and ( 18) accordingly.

Validation
To validate the convergent, accuracy, and computational efciency of the proposed numerical scheme, the advectiondifusion equation in [29] with constant difusion coefcient Dx � Dy � 1 and velocity components u � cos πy and v � − cos πx is taken.Te equation is with initial condition and boundary conditions An approximate solution of ( 20) is obtained by the numerical scheme with spatial step sizes Δx � Δy � 0.02 and time step size Δt � 0.005 and t � 0.5.A tolerance of 0.000001 is employed for both Newton's method and the Gauss-Seidel iterative method.Figure 2 illustrates the comparison between the approximate and exact solutions.In this computation, the absolute maximum error of the results obtained using the numerical scheme for 0 ≤ x, y ≤ 1 and 0 ≤ t ≤ 0.5 is 7.5275 × 10 − 5 .Furthermore, Figure 3 presents a comparison of the exact and approximate solutions for 0 ≤ t ≤ 1 at x � 0.5 and y � 0.5 and for 0 ≤ x ≤ 1 at y � 0.5 and t � 0.5.Notably, there is a strong agreement between the numerical and exact solutions, as evidenced by the fgures.In addition, Table 1 shows the maximum absolute errors, convergence rates, and CPU time.Te validation process demonstrates that the proposed numerical scheme yields accurate results and can be utilized for addressing similar partial diferential equations.Te scheme requires increased computational time when employing very small spatial step sizes.

Numerical Examples and Discussion
In this section, numerical demonstrations are provided to illustrate the practical utilization of the proposed numerical scheme in solving ADR problems.4), Dx � 0.1x 2 , Dy � 0.1y 2 [1] (see graphs of Dx and Dy in Figure 5), and boundary condition for x, y ∈ zΩ and t > 0 C(x, y, t) � 1, for x � 0 and 0.4 ≤ y ≤ 0.6, 0, otherwise.
Figure 6 displays the solution of Example 1 using the numerical scheme with spatial step sizes Δx � Δy � 0.02 and time step size Δt � 0.001 and at t � 1 and t � 2. In    16) for (a) 0 ≤ t ≤ 1 at x � 0.5 and y � 0.5 and (b) 0 ≤ x ≤ 1 at y � 0.5 and t � 0.5.
International Journal of Mathematics and Mathematical Sciences solving the error correction vector h at each iteration, the matrix A for this problem is diagonally dominant and hence the numerical scheme yields convergent solution.Te maximum absolute error obtained in the computation at t � 2 using ( 16) is 0.1459.It is observed that as the maximum absolute error decreases, we take smaller step sizes.Tere is a rapid change of concentration distribution profle.Te numerical solution of Example 2 at t � 1 and t � 2 is displayed in Figure 8. Te same grid size is used as Example 1 for the computation.Tis example targets to see the efect of difusion coefcient on the contaminant concentration distribution.As we can observe in Figure 8, higher values of the difusion coefcient in the contamination transport give wider coverage area of concentration distribution for the same advection velocity which agrees with situation observed in [12].Te maximum value of the concentration at t � 2 is attained near the entrance gate.
Example 3. Consider the ADR problem in (4) with u(x, y) � e − x cosh y and v(x, y) � e − x sinh y (see the velocity streamline in Figure 4) and having the same difusion coefcients, reaction term, and initial and boundary conditions as Example 1.
Te numerical solution of Example 3 at t � 1 and t � 3 is depicted in Figure 9. Te same grid sizes are used as Example 1 to obtain the solution.For this problem, the matrix A that appears at each iteration in solving error correction vector h is also diagonally dominant.Hence, the numerical scheme is convergent.Te maximum absolute error obtained in the computation is 0.2336.Tere is relatively slow change of concentration distribution profle between t � 1 and t � 3 as it can be observed from the graphs.Te contaminant concentration distribution profle follows the advection velocity streamlines (see Figures 4 and 9).
In Example 4, the reaction term is a third degree polynomial.Te numerical solution of this problem is { } for 0 ≤ x, y ≤ 1 and 0 ≤ t ≤ 1 using ( 16) as indicated in Table 1.For all problems, the maximum absolute error decreases as grid sizes decrease which shows the accuracy and the convergent of the proposed numerical scheme in solving the problems.Figure 4 shows the velocity stream lines of the three variable velocities in the examples.Tese streamlines are plotted using the streamline function of MATLAB.From the numerical results of problems in Examples 1, 3, and 4, it is observed that the contaminant concentration distribution profles follow the fuid fow velocity streamlines which agreed with the condition realized in [1] for small difusion coefcient.Tis reveals the applicability of the proposed numerical scheme in simulations of ADR problems.

. Conclusions
In this study, a numerical scheme is employed to obtain the solution of two-dimensional unsteady nonlinear advectiondifusion-reaction equations with velocity and difusion coefcients that vary spatially.By developing and implementing the comprehensive numerical scheme, which integrates the Crank-Nicolson Method, Newton's method, and the Gauss-Seidel iterative method, nonlinear ADR problems are successfully solved.Te absolute maximum errors and convergence rate of the numerical scheme are estimated to investigate the accuracy and efciency of the numerical scheme.Te fndings demonstrate the scheme's capability to accurately approximate solutions for ADR  International Journal of Mathematics and Mathematical Sciences problems and its suitability for simulating related twodimensional nonlinear partial diferential equations.Te scheme requires substantial CPU time during computation when a large number of spatial grid points are utilized.

Example 2 .
Let us take the ADR problem in (4) with diffusion coefcient Dx � 0.1(x + y − x 2 − y 2 ) and Dy � 0.1(1 − xy) (see Graphs of Dx and Dy in Figure 7) having the same velocity components, reaction term, and initial and boundary conditions as Example 1.
Te computed results are illustrated in graphs, and the estimated maximum absolute errors are presented in tabular form.For all examples, Example 1.Consider the ADR problem in (4) with u(x, y) � 0.3 + x, v(x, y) � 0.4 − y (see velocity streamlines in Figure