Sixth-Order Combined Compact Finite Difference Scheme for the Numerical Solution of One-Dimensional Advection-Diffusion Equation with Variable Parameters

: A high-accuracy numerical method based on a sixth-order combined compact difference scheme and the method of lines approach is proposed for the advection–diffusion transport equation with variable parameters. In this approach, the partial differential equation representing the advection-diffusion equation is converted into many ordinary differential equations. These time-dependent ordinary differential equations are then solved using an explicit fourth order Runge–Kutta method. Three test problems are studied to demonstrate the accuracy of the present methods. Numerical solutions obtained by the proposed method are compared with the analytical solutions and the available numerical solutions given in the literature. In addition to requiring less CPU time, the proposed method produces more accurate and more stable results than the numerical methods given in the literature.


Introduction
The transport of solutes by water takes place in a large variety of environmental, agricultural, and industrial conditions. An accurate prediction of the pollutant transport is crucial to the effective management of these systems [1]. In the deterministic approach with constant transport parameters with respect to time and position the advection-diffusion equation (ADE) is linear and explicit closed-form solutions can generally be derived. Because of the large variability of flow and transport properties in the field, the often transient nature of the flow regime, and the non-ideal nature of applicable initial and boundary conditions, the usefulness of analytical solutions is often limited and numerical methods may be needed [2].
Several numerical solutions of the ADE with variable parameters have been proposed in the literature. A numerical solution for one-dimensional transport of conservative and non-conservative pollutants in an open channel with steady unpolluted lateral inflow uniformly distributed over the whole length of the channel was presented by Ahmad [36]. In Ahmad [36], the flow velocity was considered to be proportional to the distance and the diffusion parameter proportional to the square of velocity to account for the lateral inflow into the channel. The governing equation was converted into an equation with constant parameters using a transformation approach. The transformed equation was then solved using a cubic spline interpolation and the Crank-Nicolson finite difference scheme

Advection-Diffusion Equation
Solute transport through a medium is described using a PDE of the parabolic/hyperbolic type. It is derived on the principle of conservation of mass and Fick's laws of diffusion. This equation is usually known as the ADE. The one-dimensional ADE may be written in subscript notation as [45]: where C is the solute concentration, D is the diffusion parameter and U is the flow velocity at a position x along the longitudinal direction at time t. In Equation (1), D and U can be rewritten as: In the above equation, D 0 and U 0 may be referred to as the initial diffusion parameter and the uniform velocity, respectively. D 0 and U 0 are constants whose dimensions depend Mathematics 2021, 9, 1027 3 of 14 upon the expressions g 1 (x, t) and g 2 (x, t). The initially solute-free state of the semi-infinite medium implies the following initial condition: Because a continuous input concentration is introduced at the origin, whereas the concentration gradient at infinity is assumed to be zero, the following boundary conditions are obtained:

Numerical Method
The numerical solution technique that we consider in this study is based on the sixthorder combined compact finite difference and the MOL approach. Detailed information and discussions about these methods are given in the following sub-sections.

Combined Compact Finite Difference Scheme
Combined compact difference (CCD) schemes are high accuracy methods, where first and second derivatives are evaluated simultaneously utilizing the Hermitian polynomial technique, as discussed in [46][47][48][49][50]. Since these methods provide first and second derivatives simultaneously, one expects an increased accuracy and economy of operations from these methods.
The analysis of the CCD schemes by depicting scale resolution [51,52] has been augmented by properties related to phase and dispersion errors represented by numerical group velocity and phase speed, and a new CCD scheme has been proposed in Sengupta et al. [48], which has been termed as the NCCD scheme.
Consider a domain with N equidistant points with a spacing h, on which a general function f (x, t) is defined. The NCCD scheme is used to simultaneously evaluate the first and second spatial derivatives, indicated by primes as f j , f j , which are evaluated at the spatial location x j , from the following implicit equations for interior nodes, in terms of the function f j values [46,48,49]: For j = 1 and j = 2: : For j = (N − 1) to N: with γ 1 = −0.025 and γ 2 = 0.09, are proposed by Sengupta [53] for better global numerical properties in the domain. Multiplicative constants in Equations (6) and (7), are fixed by matching the Taylor series expansion coefficients up to the sixth-order. Thus, we have a complete linear algebraic system for evaluating the first and second derivatives. Equation (5) are used for j = 1 and Equation (8) are used for j = N. A full-domain spectral analysis [54] for individual derivatives is performed for the CCD schemes to solve non-periodic problems. For the purpose of this analysis, we can write Equations (5)-(8) as: By decoupling these two simultaneous linear algebraic equations we get [48,53]: where More details of this global spectral analysis for CCD scheme is available in Sengupta and Bhaumik [53].

Time Integration
The method of lines (MOL) is a procedure whereby one or more partial differential equations (PDEs) are reduced to a system of ordinary differential equations (ODEs) in time by approximating the spatial derivatives using standard approaches and then integrating in time using an ODE code [55]. The MOL consists of two steps. In the first step, partial differential equations are converted into a system of time-dependent ordinary differential equations by replacing spatial derivatives with any numerical schemes. In the second step, the resulting system of ODEs is numerically integrated in time using any numerical integration rule [56]. When the MOL approach is applied to Equation (1), the time-dependent initial value problem is established as follows: The first-order system of ODEs given by Equation (22) requires N initial conditions, which are given by Equation (3). The expression θ t, C j is called the right-hand side function at the location x j and θ(t, C) can be written in the vector form as follows: As seen above, we need the boundary conditions in the calculation of derivative terms within θ(t, C) and these are given by Equation (4). The solution of the ODE system given in Equation (22) can be easily obtained by any time-integration scheme. In the MOL approach, the differential matrices (D 1 and D 2 ) and the initial and boundary conditions are first passed to the function which is needed to calculate the numerical values of θ(t, C). Note that, the numerical derivatives in θ(t, C) are calculated with high accuracy using the NCCD scheme. The resulting ODE system is then solved using the explicit fourth order Runge-Kutta (ERK4) scheme. The ERK4 scheme is one of the most popular generalpurpose integrators as it is simple to implement and has good stability characteristics [57]. The formulation of the ERK4 scheme from time step m to m + 1 can be written as:

Numerical Applications
In , and the input parameters (C 0 , D 0 , U 0 ) are the same. The parameters m and ∝, which control only temporal and spatial change, and the total simulation times (t) are different. The main goal of this approach is to observe how the behavior of the problem changes when the velocity (U) and diffusion coefficients (D) are selected as variables. If the velocity and diffusion coefficients are fixed, the second and the third problem is reduced to the first problem.

ADE with Constant Parameters
The advection-diffusion equation is solved for a transport domain in which the velocity and the diffusion coefficient are constant, i.e., g 1 (x, t) = 1 and g 2 (x, t) = 1. These conditions describe the propagation of a steep front, which is simultaneously subjected to the diffusion. The analytical solution of the ADE with constant parameters is given by Szymkiewicz [3]: In this example, a numerical solution of an ADE with constant parameters and sharp behavior was performed. The analytical solution of this example has been calculated incorrectly by various researchers [3,20,22]. Therefore, the results of these studies were not used in the comparison. Recently, Irk et al. [58] was conducted a study, by approximating the spatial derivatives with cubic B-spline collocation scheme and extended cubic B-spline collocation. They were adopted the Crank-Nicolson scheme for the time integration. These two schemes will be referred to as CN-CBSC and CN-ECBSC.
L ∞ norm errors of the ERK4-NCCD, CN-CBSC, and CN-ECBSC schemes for different time steps at t = 3000 s are compared in Table 1. In this example, the Peclet number (Pe) was selected as 5. The Courant numbers (Cr) were calculated as 0.01, 0.05, 0.1, 0.2, 0.3, 0.6, and 1, respectively, when the time steps used were ordered from the smallest to the largest. As can be seen from Table 1, the ERK4-NCCD scheme is more accurate and more stable than the CN-CBSC and the CN-ECBSC schemes.  Table 2 shows the numerical results obtained with ERK4-NCCD at t = 3000 s using different time steps. It can be observed from Table 2, the ERK4-NCCD scheme it produces the more accurate solutions for ∆t = 1 s, 5 s, and 10 s and gives acceptable solutions for larger time steps. Table 3 shows the numerical solutions of ERK4-NCCD at t = 2000 s, 4000 s, and 6000 s using the time step ∆t = 1 s. As seen from the Table 3, there is almost no difference between the analytical and numerical solutions. Numerical and analytical solutions of the ADE with constant parameters for ∆t = 1 s is presented in Figure 1. Figure 1 shows that the numerical results obtained for ∆t = 1 s coincide very well with the analytical solution. Table 2 shows the numerical results obtained with ERK4-NCCD at t = 3000 s using different time steps. It can be observed from Table 2, the ERK4-NCCD scheme it produces the more accurate solutions for ∆t = 1 s, 5 s, and 10 s and gives acceptable solutions for larger time steps.
Mathematics 2021, 9, x FOR PEER REVIEW 6 of 15 the largest. As can be seen from Table 1, the ERK4-NCCD scheme is more accurate and more stable than the CN-CBSC and the CN-ECBSC schemes.  Table 2 shows the numerical results obtained with ERK4-NCCD at t = 3000 s using different time steps. It can be observed from Table 2, the ERK4-NCCD scheme it produces the more accurate solutions for ∆ = 1 s, 5 s, and 10 s and gives acceptable solutions for larger time steps. Table 3 shows the numerical solutions of ERK4-NCCD at = 2000 s, 4000 s, and 6000 s using the time step ∆ = 1 s. As seen from the Table 3, there is almost no difference between the analytical and numerical solutions.
Numerical and analytical solutions of the ADE with constant parameters for ∆ = 1 s is presented in Figure 1. Figure 1 shows that the numerical results obtained for ∆ = 1 s coincide very well with the analytical solution. Table 2 shows the numerical results obtained with ERK4-NCCD at t = 3000 s using different time steps. It can be observed from Table 2, the ERK4-NCCD scheme it produces the more accurate solutions for ∆ = 1 s, 5 s, and 10 s and gives acceptable solutions for larger time steps.

ADE with Spatially Variable Parameters
In this example, we assumed that the diffusion coefficient is proportional to the square of the velocity, i.e., g 1 (x, t) = (1+ ∝ x) 2 and g 2 (x, t) = 1+ ∝ x. Under this assumptions, the analytical solution of the ADE with spatially variable parameters is given by Kumar et al. [59] as follows: where β = (U 0 + ∝ D 0 )/ 2 √ D 0 and δ = U 0 /(∝ D 0 ). In this example, the velocity and the diffusion coefficient increase depending on the value of ∝. The numerical results at t = 3000 s for different time steps and ∝= 0.005 m −1 are presented in Table 4. Table 4 shows that the ERK4-NCCD scheme was able to produce acceptable solutions for the maximum time step ∆t = 75 s. When the simulation parameters were increased, there was a slight decrease in the stability of the scheme according to the previous example. Certainly, this inference is also valid for other numerical methods. Figure 2 depicts the analytical and numerical solution of the ADE with spatially dependent parameters at t = 3000 s for ∆t = 1 s and ∝= 0.02 m −1 . These results can be compared with the problem presented above for ∝= 0 and ∆t = 1 s in Figure 1. When we compare Figures 1 and 2, we observe that a change in ∝ leads to a significant change in the solution of the problem. As can be observed in Table 5, this change becomes more evident when we significantly increase ∝. The L ∞ norm errors of the ERK4-NCCD scheme for solving ADE with different ∝ values at t = 3000 s are provided in Table 5.

ADE with Spatially and Temporally Variable Parameters
In this example, we assumed that the diffusion is proportional to square of the velocity parameter as considered second example. In addition, the diffusion and velocity are supposed to vary with temporally in the same proportion as considering g 1 (x, t) = g(m, t) (1+ ∝ x) 2 and g 2 (x, t) = g(m, t)(1+ ∝ x). Under this assumptions, the analytical solution of the ADE with spatially and temporally variable parameters is given by Kumar et al. [59] as following: where , g(m, t) = 1 − sin(mt) and ϑ is a dummy variable. In this example, the velocity and diffusion parameters are slightly increased depending on the m and ∝.
L ∞ norm errors of the ERK4-NCCD scheme for m = 0.1 s −1 and different ∝ values at t = 1500 s can be found in Table 6. Increasing the ∝ value leads to increasing Courant number, Cr = (U∆t)/h. Therefore, we can face with stability problems for explicit time integration methods like ERK4. Referring to Table 6, it is seen that the stability range decreased for ∝= 0.01 m −1 and ∝= 0.02 m −1 . The numerical results of the ERK4-NCCD scheme with ∝= 0.005 m −1 and ∆t = 1 s at t = 1500 s given in Table 7. As can be seen from the table, quite less L 2 and L ∞ values were obtained for all m values. Analytical and numerical solution of the ADE with spatially and temporally variable parameters at t = 2000 s for ∆t = 1 s, and ∝= 0.05 m −1 are presented in Figure 3.
is a dummy variable. In this example, the velocity and diffusion p ters are slightly increased depending on the m and ∝. ∞ norm errors of the ERK4-NCCD scheme for = 0.1 s −1 and different ∝ at = 1500 s can be found in Table 6. Increasing the ∝ value leads to increasing C number, = ( ∆ ) ℎ ⁄ . Therefore, we can face with stability problems for expli integration methods like ERK4. Referring to Table 6, it is seen that the stability ra creased for ∝= 0.01 m −1 and ∝= 0.02 m −1 . The numerical results of the ERK4 scheme with ∝= 0.005 m −1 and ∆ = 1 s at = 1500 s given in Table 7. As can  from the table, quite

Conclusions
In this study, a combined compact finite difference scheme based on the method of lines is proposed for the numerical solution of the solute transport equation with variable parameters. In this approach, the partial differential equation representing the ADE is converted into many ordinary differential equations. These time-dependent ordinary differential equations are then solved using an explicit fourth order Runge-Kutta method. Numerical examples with sharp concentration fronts are presented to demonstrate the effectiveness of the method. In the first example, we compared the numerical results of the ERK4-NCCD scheme with CN-CBSC and CN-ECBSC. It is found that the ERK4-NCCD scheme is more accurate and more stable than the CN-CBSC and the CN-ECBSC schemes. It was observed that the numerical results obtained using small time steps and the analytical results coincided with each other. Moreover, acceptable results were obtained when large time steps were used as well. The ERK4-NCCD scheme gives excellent results for long time-integration, i.e., t = 2000 s, 4000 s, and 6000 s. In second example, the velocity and diffusion coefficients were slightly increased depending on the value of ∝. When the parameters increased, there was a slight decrease in the stability of the scheme compared to the previous example. We observe that a small change in ∝ leads to a significant change in the solution of the problem. This change becomes more evident when we significantly increase the ∝. In the third example, we observed that increasing ∝ value leads to increasing Courant number. Therefore, we can face with stability problems for explicit time integration methods like ERK4. As a result, the proposed method has produced more accurate and more stable results than CN-CBSC and the CN-ECBSC schemes in solving ADE with constant and variable parameters. The proposed scheme produces fairly accurate results for the Cr ≤ 1 and the Pe ≤ 5. It is recommended to use compact upwind schemes for higher Peclet numbers. When the method given in Gurarslan (2014) is followed, the proposed scheme can be easily extended to the two-and three-dimensional ADE.