Simulation of advection–diffusion–dispersion equations based on a composite time discretization scheme

In this work, we develop a high-order composite time discretization scheme based on classical collocation and integral deferred correction methods in a backward semi-Lagrangian framework (BSL) to simulate nonlinear advection–diffusion–dispersion problems. The third-order backward differentiation formula and fourth-order finite difference schemes are used in temporal and spatial discretizations, respectively. Additionally, to evaluate function values at non-grid points in BSL, the constrained interpolation profile method is used. Several numerical experiments demonstrate the efficiency of the proposed techniques in terms of accuracy and computation costs, compare with existing departure traceback schemes.


Introduction
This study focuses on the advection-diffusion-dispersion equation described as where u := u or [u, v] T denotes the solution of (1); ν is a positive kinematic viscosity; μ is a dispersive coefficient; J := u x or Many physical phenomena can be described by the advection-diffusion or advectiondispersion equations, including various nonlinear time-dependent partial differential equations (PDEs), such as the Burgers-type and Korteweg-de Vries (KdV) type equations. In particular, the Burgers-type problems are simple models to describe various physical problems, such as turbulence flow, approximate theories of shock waves, vorticity transportation, and dispersion in porous media [9,11,13]. Many powerful methods to simulate Burgers-type problems have been proposed [2, 4, 8, 10, 12, 17, 20-23, 32, 35, 36] and have drawn research interest to find methodologies for nonlinear time-dependent PDEs in many fields, among which the characteristic methods are very attractive [14].
As one of the characteristic methods, the semi-Lagrangian method is one of the most popular numerical methods to simulate flow problems in computational fluid dynamics [1,4,6,7,27,31,33]. Among them, the backward semi-Lagrangian method (BSL) describes the movement of particles in the Lagrangian description and resets the positions of the particles to the Eulerian grid point at each time level to reduce the accumulation and collision of particles. Moreover, it is well known that BSL exhibits good stability, allowing the use of larger temporal steps than the spatial grid size by solving equations implicitly along characteristic curves of particles in the reverse direction to time evolution [3,5,30,31]. Therefore, this paper proposes a robust high-order BSL algorithm for solving nonlinear advection-diffusion-dispersion equations. To do this, we first observe the abstract model (1) in the Lagrangian description, in which individual fluid particles are followed through time. Suppose that the position of a particle at time t is described by a function x * (t). In both the Lagrangian and the Eulerian descriptions, the velocity of the particle at time t is satisfied as follows: Then the model equation (1) can be expressed as about the material derivative. By solving (2) at the same time as (3), we obtain the solution u of the model equation.
In the conventional BSL, several approximate methods are required to solve (3) simultaneously with (2): time and spatial discretization methods for (3), a departure traceback method for nonlinear ordinary differential equation (2), and an interpolation scheme to obtain the solution at non-Eulerian grid points. These approximate methods affect the overall time-to-space accuracy of the BSL. Additionally, a high-order stable temporal discretization method is critical to provide an accurate numerical approximation for the model problem, and it potentially allows for a sparser time step size than lower-order methods. Nevertheless, most current methods focus on developing spatial discretization schemes, since the strong nonlinearity of (3) is coupled with the unknown solution to (2) (or the original problem (1)), which is a significant difficulty for BSLs.
The main contributions of this paper are as follows. We propose an adaptable high-order BSL with good stability to solve nonlinear advection-diffusion-dispersion equations. It is possible to design a high-order departure traceback method, which simultaneously finds the departure points for the three previous time steps in the reverse-flow direction. The departure traceback method combines the methods of classical collocation and error correction, which is called the deferred or defect correction in some literature. The proposed departure traceback method is A-stable, verified by plotting, and determining stability function limits. Also, we use the third-order backward differentiation formula (BDF3) and fourth-order finite difference schemes for diffusion-dispersion problem discretization in the temporal and spatial domains, respectively. Moreover, the constrained interpolation profile (CIP) method [26] is used to evaluate function values at non-grid points. Several numerical experiments demonstrate that the proposed scheme provides third-order accuracy in both temporal and spatial for nonlinear advection-diffusion or -dispersion problems. In terms of computational costs, the proposed scheme is more efficient compared with current methods, including third-order AB (AB3) and fourth-order AM (AM4) [15], since it is not an iterative approach. The proposed scheme allows a more massive time step, which has a particularly significant impact on stiffer initial conditions of smaller viscosity.
The remainder of this paper is structured as follows. Section 2 provides temporal and spatial semi-discretization systems to solve a diffusion-dispersion problem (3). Section 3 designs the departure traceback method to find particle traces along characteristic curves in the reverse time direction by solving the nonlinear initial value problem (IVP) and discusses the stability of the proposed method. Section 4 provides several numerical examples to demonstrate the robustness and adaptability of the proposed scheme as well as the physical behaviors of Burgers-type equations. Finally, Sect. 5 summarizes and concludes the paper.

Semi-discretization system of diffusion-dispersion problem
This section provides the implementation of the temporal and spatial discretization method to obtain a semi-discretization system for (3). We find approximate solutions u and u at t = t m+1 on each spatial grid point by assuming that the approximate values u k i and u k i,j of u(t k , x i ) and u(t k , x i,j ) for k ≤ m and particle departure points x * i,j (t m-κ ), κ = 0, 1, 2, along the characteristic curve arriving at the grid point at target time t = t m+1 have already been calculated at all grid points where a finite set of mesh points {x i } or {x i,j }, (i = 1, 2, . . . , J x , j = 1, 2, . . . , J y ) is equally spaced.
For the two-dimensional (2D) case, we first evaluate (3) at (t, x) = (t m+1 , x i,j ) and then apply BDF3 to approximate the material derivative on the left side of (3). Then we obtain the Helmholtz equation , whereν = ν 6h 11 andμ = μ 6h 11 for a uniform temporal step size h. The one-dimensional (1D) Helmholtz equation can be obtained in the same way as the 2D Helmholtz equation.
• 2D problems identity matrices of size J x = J x -1 and J y = J y -1, respectively; andD (k) x andD (k) y are constructed from D (k) x and D (k) y (k = 2, 3) using interior elements only. (See Appendix A for more details.) Approximate values for particle departure points x * (t m-κ ), (κ = 0, 1, 2) along the characteristic curve arriving at the grid point at target time t = t m+1 are required to obtain full discretization systems for (5) and (6), which is discussed in Sect. 3. Moreover, since x * i,j (t m-κ ), (κ = 0, 1, 2) may not coincide with the grid points, an appropriate interpolation scheme is also required.

Proposed departure traceback scheme
This section develops the proposed third-order composite scheme without iteration to solve the strongly nonlinear IVP (2) by combining error correction [28,30] and classical collocation methods. Let x * i,j (t) be the solution for the following IVP: where x i,j is an arbitrary grid point and u is the solution of the model problem (1).
To find the departure points by solving (2), we first consider the perturbed asymptotically first-order ordinary differential system with initial value ψ i,j (t m+1 ) = 0 2,1 , where 0 r,s denotes the zero matrix of size r × s, and y i,j (t) is a qth-order approximating polynomial for x * i,j (t) obtained by solving (7) that satisfies We choose the first-order polygon y i,j defined as which allows the proposed method to attain up to fourth-order accuracy eventually.
Remark 1 To construct higher-order methods, y i,j (t) can be chosen as an explicit secondorder polynomial, which implies that the proposed method could achieve sixth-order accuracy. To maintain fourth-order and above temporal accuracy overall, a stable method for fourth-order and above must be employed rather than BDf3 to discretize the total time derivative.
The proposed method is then completed by developing a method to solve equation (8) with initial condition ψ i,j (t m+1 ) = 0 2,1 on [t m-2 , t m+1 ]. By applying (22) with (23) to (8), ) T ] T ,ĥ := -3h and I 2 is the identity matrix of size 2. However, since F(t m+1 , ψ i,j (t m+1 )) contains the unknown solution u for (1) at time t = t m+1 to be obtained, an approximation tool is required. Therefore, we use extrapolation to find where We proceed by simplifying (11), and expressing in matrix form where I is the 2 × 2 identity matrix. Finally, by combining (9) with (12) and truncating the asymptotic error term O(h 2+p ), the approximation x * ,m-κ i,j of the three departure points can be approximated as follows:

Linear stability of the proposed departure traceback scheme
In this subsection, to examine the stability of the proposed departure traceback scheme, the scheme is applied to the Dahlquist problem dx * (t) dt = λx * (t) with initial x * (t m+1 ) = x * ,m+1 and an eigenvalue λ. For each departure position at the previous time steps t = t m-κ (κ = 0, 1, 2), stability functions of the proposed departure traceback scheme hold in the following theorem.

Theorem 1
The three departure points of the proposed departure traceback scheme (13) are all A-stable schemes.
Moreover, the existing departure traceback schemes that will be compared with our method in the next section have the stability functions [15] • AB3: Remark 2 Note that while all stability functions of the proposed method are A-stable, the stability functions of AB3 and AM4 are neither A-stable nor L-stable at t = t m-κ , κ = 1, 2.
The impoverished stability constrains the choice of time step size h, when the problem has small viscosity, hence losing the advantages of the BSL.

Numerical experiments of advection-diffusion or -dispersion equation
To assess the efficiency of the proposed method, we simulate six examples. To measure the accuracy, we calculate the computational errors using the maximum norm, L 2 norm and relative L 2 norm, respectively, defined by where u k (x i ) is the exact solution at t = t k and each grid point x i and u k i is an approximation. E ∞ (t k ) and E R2 (t k ) are similarly defined for the 2D case. The numerical convergence for the proposed method, when the exact solution of the problem is available, is calculated as where E k represents E ∞ (t k ) or E R2 (t k ) and τ represents h or x, respectively. For convenience, we label the proposed departure traceback schemes and current schemes as follows.
Remark 3 To find departure points x * i,j (t) at three previous time levels t = t m-κ (κ = 0, 1, 2) arriving on the grid points at target t = t m+1 , AB3 and AM4 are implicit schemes in our BSL, since they simultaneously find the departure points by solving the strongly nonlinear IVP (2) with initial x * i,j (t m+1 ) = x i,j in the reverse direction to time. Therefore, an iteration method, such as the Newton-like iteration scheme, is required for the current departure traceback schemes. We employed the fixed iteration with maximum iteration iter = 20 and tolerance tol = 10 -7 , considering computational costs.

Linear advection equation
This subsection provides experimental results to demonstrate convergence for the proposed BSL for a linear advection problem. To do this, we consider the linear advection problem wherev(t, x) := 1 2 x tan(t) is a known velocity, and the analytical solution of (15) is ρ(t, x) = exp(x 2 cos(t)), x 2 = x 2 + y 2 . Using the semi-Lagrangian formulation, we rewrite (15) as ρ T (t, x) = 0, along the characteristic curve x * (t; x, s) satisfying dx * (t;x,s) dt =v(t, x), which is a different type from the equation discussed previously. Specifically, the distribution function ρ is constant along the characteristic curve, i.e., ρ(t m+1 , x i,j ) = ρ(t m , x * i,j (t m )), and hence, we only need the departure point x * i,j (t) at t = t m . To verify the temporal convergence order of each method, we measure E ∞ and E R2 and check the computational cost (cpu) for a fixed spatial resolution J x = J y = 12 × 100 by varying the temporal step size h from 2 -3 to 2 -9 . The results (the errors versus temporal step size h with a logarithmic scale based on 10) are plotted in Fig. 2 and illustrated in Table 1. The figure shows that the EAC3 and AB3 exhibit third-order convergence rates and EAC4 and AM4 have fourth-order convergence rates. However, EAC3 and EAC4 have greater computational costs than AB3 and AM4, since EAC3 and EAC4 simultaneously compute all three departure points with high accuracy by solving the linear system (12), whereas AB3 and AM4 are designed to calculate with good stability only at t = t m . These results seem that AB3 and AM4 are more efficient than the proposed methods for linear hyperbolic-type problems.

1D Burgers' equation
This section investigates the performance of the proposed method combined with BDF3 for the nonlinear advection-diffusion problem and demonstrates its superiority compared with AB3 and AM4. As an illustrative example, we consider a 1D Burgers' equation with the following initial condition u(0, x) = 2νπ sin(π x) σ +cos(π x) , x ∈ [0, 1] and the Dirichlet boundary condition. The analytic solution for this example is given as [32,35] u(t, x) = 2νπ exp(-π 2 νt) sin(πx) σ + exp(-π 2 νt) cos(πx) , where σ > 1 is a parameter determining the shape of the initial function, and σ and stiffness are inversely proportional. We first numerically investigate temporal convergence for the four methods by measuring the errors E ∞ and E R2 for the smooth initial condition with σ = 100 and ν = 0.1. The experiment is conducted by varying the time step h from 2 -3 to 2 -7 with spatial resolution J x = 400, as shown in Table 2. All four methods, EAC3, EAC4, AB3, and AM4, combined with BDF3 achieve third-order time convergence. To check spatial convergence for the proposed method, simulation is performed until final time t f = 1.0 for various spatial resolutions J x from 2 3 to 2 7 with h = 1/4000 as shown in Table 3. The proposed EAC3 and EAC4 achieve at least third-order spatial accuracy. To demonstrate the superiority of EAC3, E R2 and cpu are obtained via the three methods, EAC3, AB3 and AM4, by using the time step size h = 1/100, 1/200, 1/400, and 1/800 and the spatial resolution J x = 1000 until t = 1.0; we display the cpu versus E R2 in Fig. 3. The figure clearly shows that EAC3 significantly outperforms the other methods, particularly as σ reduces (i.e., increasing initial stiffness or flow velocity). The initial condition becomes steeper and the solution speed increases as σ decreases from 100 to 1.2. Additionally, we measure the E ∞ (t) and E 2 (t) by    EAC3 and other existing methods introduced in [29,35]. The results are listed in Table 4 and show that EAC3 is more accurate than the compared methods presented in [29,35].
To demonstrate the robustness and superiority of the proposed method, we numerically estimate the convergence order of EAC3 and cpu. We measure convergence rates for both solutions u and v simultaneously at t = 1.0 for the same temporal and spatial step sizes h = x from 2 -3 to 2 -7 . As shown in Table 5, EAC3 clearly exhibits the expected third-order convergence rate. To demonstrate the accuracy and superiority of EAC3, we implement the experiment of the example and compare the numerical results obtained by EAC3, AB3, and AM4. We compare the errors, E ∞ and E R2 , obtained by the three methods for different time levels t = 0.5, 1.0, 2.0, and 3.0 with h = 2 -7 and J x = 1024. The results are listed in Table 6 and Fig. 4. EAC3 exhibits significantly superior accuracy, showing   excellent agreement with the exact solutions and slightly superior computational costs compared with AB3 and AM4 (see Fig. 4(b)).

KdV-Burgers' equation
This section investigates the adaptability of the proposed method for the problem with dispersive flow. Consider the following KdV-Burgers equation: with the exact solitary wave solution as given by ( [25]) where a 0 = ν 2 25εη , a 1 = ν 10η , a 2 = 6ν 2 25η ; ε, and ν and η are arbitrary constants. The KdV-Burgers' equation applies to studying weak effects of dispersion, dissipation, and nonlinearity in waves propagating in a liquid-filled elastic tube [34].
To investigate the accuracy and efficiency of the proposed method, we solve the problem with fixed sparse grid sizes h = 0.1 and J x = 10 by varying ε, ν, and η, as shown in Table 7. We measure errors E ∞ and E R2 at t = 0.3, 0.6, and 0.9, and the numerical solutions show very good agreement with the exact solutions. We observe the graphical representation of the time-progressive evolution for the traveling solution with small μ = 0.01, 0.001. Hence, we employ h = 0.02, J x = 800 and h = 0.01, J x = 1000 as shown in Fig. 5. The figure shows  Table 8 Convergence rates of EAC3 for Example 4.5 at t = 0.5 and t = 1.0 with ν = 0.1

2D Burgers' equation
In this subsection, to observe numerical accuracy in 2D nonlinear advection-diffusion equations, we perform the test problem using the 2D unsteady Burgers' equation, with the analytic solution from [37] u(t, Initial and boundary conditions are taken directly from the analytic solution.
We first examine the accuracy of EAC3 for ν = 0.1 by varying the same time step and space grid sizes, h = x = y from 1/50 to 1/400 for two different time levels. As shown in Table 8, EAC3 provides third-order convergence.
We also observe the behaviors of the numerical solutions obtained by EAC3 for different viscosities ν = 0.1, 0.01, and 0.001 at t = 1.0, 2.0, and 3.0 using the time step, and spatial grid sizes h = x = y = 0.01, 0.002, and 0.0005 are used. Figure 6 shows that EAC3 has good performance and works well for small viscosity coefficients. Finally, by comparing the errors and computational cost cpu with the results of AB3, we demonstrate the efficiency of the proposed method at two different time levels t = 0.05 and 1.0. The results displayed in Table 9 show that the proposed method EAC3 provides more accurate solutions with lower computational cost than AB3, clearly indicating that EAC3 is more efficient than AB3.
To observe the behaviors of the numerical solutions, the graphical results of the time evolution for the numerical solutions obtained by EAC3 are observed for small viscosity ν = 0.0001 at three different time levels t = 1.0, 2.0, and 3.0 with h = 0.001 and x = y = 1200. The results depicted in Fig. 7 show good performance with small viscosity ν = 0.0001.
To investigate the accuracy and efficiency of EAC3, we measure the errors of u, v, and cpu for different viscosities ν = 1.0, 0.1, 0.01, and 0.005 for various h and J x = J y at t = 0.5 and 2.0. Table 10 shows that the numerical solutions of EAC3 exhibit very good agreement with the exact solutions for various viscosity coefficients.

Conclusions
In this paper, an adaptable and stable BSL was developed to solve advection-diffusiondispersion problems induced from fluid dynamics by introducing a high-order depar-  ture scheme, which is a A-stable with up to fourth-order accuracy. In particular, the scheme is designed to solve the strongly nonlinear IVP and to have high-order accuracy while having good stability for all approximation of the departure points at three previous times. It is experimentally shown to be more suitable for solving advectiondiffusion-dispersion type problems than other conventional high-order departure traceback methods. Moreover, it is more efficient, compared with other methods, in that the proposed method does not need an iterative technique that is required to solve the strong nonlinear IVP. Through several numerical experiments, it was shown that the proposed scheme is superior to the conventional high-order departure traceback schemes and the latest solvers of Burgers' equation, in terms of accuracy and computational costs.

Appendix A: Finite difference method for space discretization
To implement the proposed BSL, we require fourth-order finite difference formulas of order k (k = 1, 2, 3) to approximate a function's partial derivatives. Therefore, this appendix introduces fourth-order finite difference matrices of size (J x + 1) × (J x + 1) from the discrete partial derivative operators for the Dirichlet boundary condition where the matrix elements were obtained from MATHEMATICA. Using these finite difference matrices [16], partial derivatives ∂ k ∂x k f i or ∂ k ∂x k f i,j (k = 1, 2, 3) of a given function f can be approximated as follows.