Development and Research of a Modiﬁed Upwind Leapfrog Scheme for Solving Transport Problems

: Modeling complex hydrodynamic processes in coastal systems is an important problem of mathematical modeling that cannot be solved analytically. The approximation of convective terms is difﬁcult from the point of view of error reduction. This paper proposes a difference scheme based on a linear combination of the Upwind Leapfrog scheme with 2/3 weight coefﬁcient, and the Standard Leapfrog scheme with 1/3 weight coefﬁcient. The weight coefﬁcients are obtained as a result of solving the problem of minimizing the approximation error. Numerical experiments show the advantage of the developed scheme in comparison with other modiﬁcations of the Upwind Leapfrog scheme in the case when the convective transport prevails over the diffusion one. The proposed difference scheme solves transport problems more effectively than classical difference schemes in the case when the Péclet number falls in the range from 2 to 20. It follows that the considered difference scheme allows hydrodynamic problems to be solved in regions of complex shape effectively.


Introduction
Solving problems of hydrodynamics is an important field of science that allows for the modeling and forecasting of adverse and dangerous natural and anthropogenic phenomena in coastal systems. Examples include run-up phenomena under wind influences, the spread of pollutants, including oil and petroleum products, blooming waters, and the mass death of commercial fish species due to a lack or absence of dissolved oxygen. The problem of modeling unsteady transport processes in a three-dimensional domain of complex shape cannot be solved analytically, and requires the use of numerical methods [1][2][3]. At work [4], the finite volume method on structured non-uniform grids for 1D and 2D shallow water equations is considered.The paper [5] presents a solution to a two-dimensional non-stationary convection-diffusion problem using the boundary element method of radial integration. This algorithm requires only boundary discretization and some internal points, instead of internal cells. The proposed method accuracy and efficiency are researched in this work. Adaptive iterative splitting methods to solve convection-diffusion-reaction problems is proposed in paper [6]. Adaptive techniques are embedded into the splitting approach, when developing this method. The paper [7] presents an upwind difference scheme for interface flux approximations for the finite volume method. The upwind approach is proposed for spatial interface approximation, and the second-order formulation is used for the temporal approximation in this work. The finite volume method on unstructured grids is proposed for transport problems in which convection dominates diffusion [8]. The research [9] considers a high-order finite volume complete flux scheme for a one-dimensional advection-diffusion-reaction equation. In this work, the accuracy increases to the fourth order due to applying Gauss-Legendre quadrature rules to the integral representation. The presented scheme works for both diffusion dominated and advection dominated flow.
A number of papers present reviews and developments of various numerical approaches to linear [10][11][12][13][14][15] and nonlinear [11,13,15] transport, and convection-diffusion problems, as applied to finite element methods and finite difference methods. In the numerical solution of the transport problem, a problem related to the insufficient accuracy of the approximation of convective terms arises. The specialized difference schemes for solving convection-diffusion problems were developed and described in [16]. Upwind Leapfrog schemes were used to solve aeroacoustics problems [17,18]. The monotonicity of difference schemes is researched in [19][20][21], including the Upwind Leapfrog schemes for the transport problem. A number of papers were devoted to the consideration of various modifications of the Upwind Leapfrog scheme, including [22], in which it was proposed to use weight multipliers to eliminate oscillations when solving a nonlinear transport problem. To solve this class of problems, a difference scheme was proposed, which is a linear combination of Standard Leapfrog and Upwind Leapfrog schemes [23]. The considered schemes are included in a linear combination with weight coefficients of 1/3 and 2/3, respectively, which are obtained from the solution of minimizing the order of the approximation error problem [24].
A linear combination of two schemes with similar properties often has better properties than the original schemes, since there is a mutual compensation of approximation errors. The purpose of this work is to determine the range of the grid Péclet number values, in which the proposed scheme, which is a linear combination of the Upwind Leapfrog scheme with 2/3 weight coefficient, and Standard Leapfrog scheme with 1/3 weight coefficient, obtained as a result of solving the problem of minimizing the approximation error, has better accuracy compared to classical schemes, including modifications of the Upwind Leapfrog scheme with limiters. Solutions of test problems of transport in areas of complex shape using the developed scheme are considered in the work [25]. These are problems such as the test problem of the fluid flow between two coaxial circular cylinders at different grid Péclet numbers, and the problem of suspended particles transport during soil dumping in areas with complex geometry. An application of this scheme to solve the problem with nonlinear dispersion wave equations is described by the model Korteweg-de Vries equation considered in [26]. The advantage of the presented difference scheme in comparison with the standard ones in problems with the predominance of convective transport over diffusion is shown.
The article is devoted to the construction and study of a difference scheme, which is a linear combination of the Upwind and Standard Leapfrog schemes. Section 2 discusses the classical schemes of Upwind and Standard Leapfrog through the example of solving the transport problem. Section 3 describes the proposed scheme based on a linear combination of Upwind and Standard Leapfrog schemes. In Section 4, the results of solving a test problem based on the considered scheme are compared with other modifications of the Upwind and Standard Leapfrog schemes. In Sections 5 and 6, the accuracy of solving the model diffusion-convection problem based on the proposed scheme is investigated for different values of the Péclet grid number. In Conclusions, an analysis of the results obtained in the work is given.
A uniform computational grid is introduced: where is the number of nodes in space; L is the characteristic dimension of the computational domain, ω τ = t n | t n = nτ, n = 0, N t , N t τ = T , where τ = t n+1 − t n = const is the time step; N t is the number of time layers; T is the upper bound in time.
The following finite-difference schemes can be used for the numerical solution of the problem (1): • The Upwind Leapfrog scheme: where The Standard Leapfrog scheme: Test problem I. Find the solution to the equation: For the numerical solution of the problem, we introduce a uniform computational grid: Parameters of the calculation grid: time step τ = 0.02 s, space step h = 1 m, the length of the time interval T = 100 s, the length of the space interval L = 100 m. Figure 1 shows the exact solution of the test problem I, and numerical solutions based on the difference schemes (3) and (4), the Upwind Leapfrog difference scheme with TVDlimiters, and the Standard Leapfrog difference scheme with TVD-limiters [27].   Figure 1a,b show that the solutions of test problem I based on difference schemes (3) and (4) have dispersion-induced oscillations in front of the leading and trailing fronts of the wave, respectively. We also note that for the solution based on the difference scheme (3), high harmonics have a "phase velocity" that is higher than the real one, and for the solution based on the scheme (4), it is lower. Comment 2. Difference schemes (3) and (4), due to their non-dissipativity and nonmonotonicity, are not recommended for solving transport problems in the case of discontinuous initial data. Their application to such problems is possible by implementing procedures that monotonize the solution; for example, by adding artificial viscosity. For scheme (3), it is recommended to use a non-linear flow correction (TVD-limiters) [27] based on the maximum principle. Figure 1c,d present the solutions of the test problem I based on the Upwind and Standard Leapfrog schemes with TVD-limiters [27]. Note that the difference scheme obtained as a result of a linear combination of the Upwind Leapfrog [27] and Standard Leapfrog schemes is no longer dissipative, but the dispersion properties of the scheme can be improved and the solution obtained based on a combination of difference schemes on discontinuous functions may have smaller oscillations.

Linear Combination of Upwind and Standard Leapfrog Schemes
Consider the Upwind Leapfrog difference scheme (3) (the case u ≥ 0). Let's use the decomposition of the functions q n±1 i in the Taylor series in the node (i, n) neighborhood: Let's use the decomposition of the functions q n i−1 and q n i in the Taylor series in the node (i − 1/2, n) neighborhood: Taking into account the scheme (1), we have: where c = uτ/h is the Courant number. Consider the Standard Leapfrog difference scheme (4). Let's use the decomposition of the functions q n i±1 into a Taylor series in the node (i, n) neighborhood: Substitute (5) and (11) into (4): With (8), the equality (12) will take the form: Let's consider a linear combination of the Upwind Leapfrog (3) and Standard Leapfrog (4) difference schemes with weight coefficients 2/3 and 1/3, respectively: Taking into account (8), (10), and (13), the difference scheme (14) has a local approxi- The approximation error of difference schemes (3) and (4) is O τ 2 + h 2 , and based on the presented estimate, the approximation error of scheme (14) is O τ 2 + ch 2 , where c = uτ/h is the Courant number. Thus, it is preferable to use scheme (14) for solving problems in the case of small values of the Courant number. Based on Fourier analysis, the increase in the accuracy of the solution when using scheme (14) occurs due to the mutual compensation of errors in the "phase velocity" of the harmonics of the solutions obtained on the basis of schemes (3) and (4). The stability and dispersion properties of a linear combination of the Upwind and Standard Leapfrog difference schemes are discussed in detail in [23,24]. Figure 2 shows the exact solution of test problem I, and a numerical solution based on the scheme (14). The calculated time intervals T were 100 s and 500 s.

Stability and Dispersion Properties of the Difference Scheme Obtained as a Linear Combination of the Upwind and Standard Leapfrog Schemes
Let's use the harmonic method to study the stability of the scheme (14) obtained as a linear combination of the Upwind and Standard Leapfrog schemes. Substitute q n i = ϕ n e jki into scheme (14) (case u ≥ 0): The solution of this quadratic equation with respect to ϕ is written in the form: Assume that the calculation takes the value of the root of a complex number with a non-negative real part (the argument of a complex number belongs to a semi-interval [−π/2, π/2)). Let's consider the case c = 0: From here, we obtain: We studied the values of the function |ϕ 1 (c, k)| ∈ [0, 1] for k ∈ [0, π] and c ∈ [0, 1]. Figure 3 shows the values of the module and the argument of the function ϕ 1 (c, k), depending on the values of the parameters k and c.  Let's study the dispersion properties of scheme (14) in the region of small Courant numbers. Consider the difference between the functions ϕ 1 (c, k) and φ(c, k) = e −jck : (17), and taking into account the fact that grid refinement leads to a linear increase in the number of layers with respect to time, it follows that for small Courant numbers, the error in solving the transport equation based on scheme (14) is proportional to the function:

Results of Test Calculations
Let's compare the calculations based on the scheme (14) with the results obtained using other difference schemes.
First, consider a linear combination of the central difference scheme and the Upwind Leapfrog scheme [27]: Second, consider a two-parameter difference scheme [28,29]: Figure 4 shows the exact solution of the test problem I and numerical solution based on schemes (14), (18), and (19), the Upwind Leapfrog scheme with TVD-limiters, and the Standard Leapfrog scheme with TVD-limiters [27]. The time step is τ = 0.02 s (Figure 4a  The solution error is calculated using the formula: where q i is the exact solution value of the problem in the node i ∈ ω h , andq i is the numerical solution that depends on the time step value with a fixed step size in space [30].

Accuracy of Solving the Heat Conduction Problem
Consider the heat equation with constant coefficients: We assume the necessary smoothness of the functions included in (21)- (23), and the consistency of the initial and boundary conditions.
Find an analytical solution to the heat conduction problem (21) under certain assumptions. We will assume that the functions in (21) can be represented as finite sums, which are expansions in a finite trigonometric basis: where ω = π/L, C It should be noted that in the case of the tabular method of setting u 0 , for example, on a spatial grid, the series will be limited to N by the harmonic, and an interpolation trigonometric polynomial is used to restore the continuous function, where N is the number of discrete values of the function.
Next, we consider functions having a derivative of the order α satisfying the inequality f (α) (x) ≤ K, with a period of 2π. There is an estimation of the residual term of the series (24) for any natural α [31,32]: Substitute the functions u and f into the heat Equation (21): Changing the order of operations of differentiation and summation in (25), and calculating the derivative of a spatial variable: Taking into account the linear independence of the functions sin(ωmx) for different m, (26) can be written in the form: The solution of the scheme (27) will take the form: After the transformations, taking into account the given initial and boundary conditions (22) and (23), we obtain the following representation for the solution function [33]: Difference scheme for the heat equation. For the numerical solution of the problem (21), we use the uniform computational grid (2). Then, the discrete analog of (21) can be written in the form: where q n+σ i = σq n+1 i + (1 − σ)q n i , σ ∈ [0, 1]-weight of the scheme [34]. The necessary stability condition (case σ = 0) obtained on the basis of the harmonic method leads to the following inequality [2,3]: Despite the fact that estimate (31) is a strict limitation for explicit difference schemes, in practice, the time step should be taken even less.

Test problem II. Find the solution of the equation:
Parameters of the calculation grid: the time step τ is in the range from 0.001 s to 10 s, the space step h = 1 m, and the length of the time interval T = 60 s. Figure 7 shows the error of solving test problem II based on the scheme (30). The calculation error is calculated using the Formula (20). The value of the time step τ 0 related to the value τ max = h 2 /2µ (obtained from Equation (31)) is postponed along the horizontal axis (τ 0 = τ/τ max ). In order for the percentage error of the explicit scheme (30) to be equal to 0.01 %, it is necessary to take the value τ 0 = 0.0717; in the case of using the proposed scheme with weights, the parameter is τ 0 = 5.1858.
Test problem III. Consider the problem that arises when modeling the transport of suspended matter in shallow reservoirs. Find a solution to the two-dimensional diffusion equation for a region elongated in one direction: with initial conditions q(x, y, t)| t=0 = (θ(1100 − x) − θ(900 − x))(θ(3 − y) − θ(2 − y)), where 0 ≤ x ≤ L x , 0 ≤ y ≤ L y and boundary conditions are in Dirichlet form.
For the numerical solution of the test problem III, we introduce a uniform computational where h x = x i+1 − x i = const and h y = y j+1 − y j = const are the space steps; N x and N y are the number of nodes in space; L x and L y are the characteristic dimension of the computational domain, ω τ = t n | t n = nτ, n = 0, N t , N t τ = T , where τ = t n+1 − t n = const is the time step; N t is the number of time layers; and T is the upper bound in time.
The discrete analog of test problem III can be written in the form: where q n+σ i,j = σq n+1 i,j + (1 − σ)q n i,j , σ ∈ [0, 1]-weight of the scheme [34]. Parameters of test problem III: µ x = 100 m 2 /s, µ y = 0.5 m 2 /s, L x = 2000 m, L y = 5 m. Parameters of the calculation grid: steps in the space are h x = 100 m and h y = 0.5 m; the time step τ is in the range from 0.001 s to 10 s, and the length of the time interval is T = 600 s. Figure 8 shows solutions to test problem III based on scheme (30). From Figures 7 and 8, we see that in order to achieve the same approximation error, the time step limit for an explicit scheme is significantly stricter than for a scheme with weights. In order for the percentage error of the explicit scheme to be equal to 1 %, it is necessary to take the value τ 0 = 0.01376; in the case of using the scheme (30) with weights, the parameter is τ 0 = 0.34844. Comment 5. The solution obtained on the basis of scheme (30) in the case of σ = 0, i.e., based on the explicit scheme, is stable under the constraint τ ≤ O(h 2 ) [35]. Scheme (30) at σ ≥ 0.5 has no time step restrictions. In practice, when solving problems based on the difference scheme (30) where r is the dimension of space and the parameter ∆ ∼ τ/τ max , where τ is the size of the time step that must be taken so that the calculation error is in an acceptable range, and τ max is the size of the time step obtained from stability conditions for the explicit scheme, while the constraint τ 0 ≤ ∆ is satisfied. In the general case, for solving problems of diffusion-convection, the value of the parameter ∆ is determined on the basis of computational expediency: if the step τ is taken too large, then the calculation error will be large, and if the step τ is small, then the computational costs will be high. For scheme (30), at σ = 0, it is recommended to use the value ∆ = 0.01, and at σ = 0.5, it is recommended to use the value ∆ = 0.3.
Comment 6. Figure 5b shows that in the case of values of the Courant number up to 0.1, the proposed scheme (14) for solving test problem I (convection problem) is more accurate than the Upwind Leapfrog difference scheme with TVD-limiters. From the solutions of the test problem II and the test problem III (diffusion problem) based on scheme (30) at σ = 0 and with respect to Comment 5, a constraint on the time step τ ≤ ∆ × τ max was obtained, where ∆ = 0.01, τ max = h 2 /2µ. Based on the obtained constraints on the value of the Courant number c = |u|τ/h ≤ 0.1 with the limitation on the time step, we obtain ∆|u|h/2µ ≤ 0.1 or Pe ≤ (2 × 0.1)/∆ = 20, where Pe = |u|h/µ is the grid Péclet number [36]. Note that we consider the case of the absence of monotonicity of schemes constructed on the basis of central-difference approximations (the case Pe > 2). Therefore, the proposed approximation of the convective transport operator is effective for 2 ≤ Pe ≤ 20.

Accuracy of Modified Difference Scheme
Consider the nonstationary convection-diffusion equation [37]: with initial and boundary conditions: q(0, x) = q 0 (x), q(t, 0) = q(t, L) = 0. For exact solution of Equation (33), write the finite sum of the trigonometric Fourier series for the function q(x, t) in complex form: where j = √ −1, ω = π/L, m is number of modes, and C m (t) = 2 is complex amplitude of the m-th mode. Substitute (34) in (33): The functions exp (jmωx) are linearly independent; therefore, the last equation can be written in the form: The solution of Equation (35) has the form: Based on the initial and boundary conditions for (33), (34), and (36), the exact solution of Equation (33) can be written in the form: Consider the case u ≥ 0. For a numerical solution of the problem (33), we use the uniform grid (2). The approximation of Equation (33) has the form: where Substitute (39) in (38): The functions exp (jmωhi) are linearly independent; therefore, (40) can be written in the form: When τ → 0, (41) takes the form: Lemma 1. When approximating problem (33) using the difference scheme (38) for each mode of solution function q, the convective exchange rate u and the diffusion exchange rate µ are less than the real values, and differ by and

respectively.
Proof. Scheme (42) can be written in the form: Taking into account (36), the last expression, which is a solution based on the difference scheme (38), corresponds to the solution of the equation Analyze the approximation error of the convective term in space. Let's introduce replacement jωmh = s, when Let's use the decomposition of the functions e ±s in the Maclaurin series: From the Expression (45), we can obtain that scheme (34) approximates the convective term of Equation (33) with the third order of accuracy in space.
Analyze the approximation error of the diffusion term in space: From the Expression (46), we can obtain that the scheme (38) approximates the diffusion term of Equation (33) with the first order of accuracy in space.
Note that the variable r = π/ωmh describes the number of nodes per half wave period, and π > ωmh. Therefore, the accuracy of the solution depends on the number of nodes per half wave period. Figures 9 and 10 show functions and α 2 (r) = 1 − 3(− exp (jπ/r) + 2 − exp (−jπ/r)) π 2 /r 2 (2 + exp (−jπ/r)) , that describe the dependence of the approximation error of the convective and diffusion terms, respectively, by the difference scheme (38) on the number of nodes in comparison with the central difference scheme [1].  Based on Figures 9 and 10, we can conclude that difference scheme (38) is applicable for problems in which convective transport prevails over diffusion.

Approximation of the Convection-Diffusion Problem
To approximate the convection operator in problem (33), we will use the difference operator of the scheme obtained as a result of a linear combination of the Upwind and Standard Leapfrog schemes: Test problem IV. Find a solution to the equation: The exact solution of the test problem IV can be represented as [33]: In Figure 11,     The approximation of the problem (33) based on explicit central difference schemes can be written in form: (48) Figure 13 shows the error function Ψ (in the norm L 1 ) of the numerical solution of the test problem IV (47) and the central difference scheme (48), depending on values of the grid Péclet number. Figure 13 shows the error function Ψ (in the norm L 1 ) that occurs when approximating the solution of test problem IV. The difference scheme (47) and the central difference scheme (48) were used in this case. The value of the approximation error is plotted along the Oy axis, and the value of the grid Péclet number is plotted along the Ox axis. The space step h = 1 m, the time step τ = 0.02, the value of the spatial interval L = 200 m, the value of the time interval T is 100 s, and the range for the diffusion coefficient-from 5 × 10 −3 to 5 m 2 /s, for these graphs.  Figure 13, which demonstrates the calculation error of the test problem IV, shows that the central difference scheme (48) has a higher accuracy compared to the difference scheme (47) for 0 ≤ Pe ≤ 2. Thus, the proposed modification of the Upwind Leapfrog schemes (47) in this paper has a smaller error compared to other considered difference schemes for solving the problems of continuous medium transport described by the diffusion-convection equation in the range of grid Péclet numbers from 2 to 20.

Discussion
The proposed difference scheme is a linear combination of the Upwind Leapfrog scheme with a 2/3 weight coefficient, and the Standard Leapfrog scheme with a 1/3 weight coefficient. The weight coefficients were determined as a result of solving the problem of minimizing the approximation error. This scheme practically does not have grid viscosity when solving convection-diffusion problems; therefore, it allows for more accurate solutions to be obtained in the case of a significant predominance of convective processes over diffusion ones. When modeling hydrodynamic processes in real areas with complex geometry (currents in river channels, estuary zones, along areas of land strongly protruding into the sea-sand bars), large values of grid Péclet numbers arise. The application of the developed vortex-resolving scheme to the problems of hydrodynamics allows small-sized vortexes that occur in coastal systems to be more accurately reproduced.
The proposed difference scheme has an advantage over the classical Upwind and Standard Leapfrog schemes in the case of small Courant numbers. In [24], the approximation error of the developed scheme, which is equal to O ch 2 , is researched, while for classical schemes, the approximation errors are equal to O h 2 . However, for cases of interest from the point of view of applications, the constant c is much less than unity, which implies the advantage of the proposed difference scheme.
In the process of the numerical experiments, we compared the approximation error in the norm of the grid space L 1 for the proposed difference scheme, the scheme that is a linear combination of the Upwind Leapfrog scheme, and the scheme with central differences, as well as the two-parameter difference scheme, which has the third order of accuracy for solving the convection-diffusion problem. The constructed graphs showed the dependence of the errors in solving test problems on the values of the Courant numbers. The study showed that the proposed difference scheme and the scheme based on a linear combination of the Upwind Leapfrog scheme and the central difference scheme have a smaller error compared to other considered schemes. In addition, the developed difference scheme, based on a linear combination of the Upwind and Standard Leapfrog schemes, has an advantage in approximation accuracy compared to other considered schemes for the case of Péclet number values from 2 to 20.