Convergence Improved Lax-Friedrichs Scheme Based Numerical Schemes and Their Applications in Solving the One-Layer and Two-Layer Shallow-Water Equations

The first-order Lax-Friedrichs (LF) scheme is commonly used in conjunction with other schemes to achieve monotone and stable properties with lower numerical diffusion. Nevertheless, the LF scheme and the schemes devised based on it, for example, the firstorder centered (FORCE) and the slope-limited centered (SLIC) schemes, cannot achieve a time-step-independence solution due to the excessive numerical diffusion at a small time step. In this work, two time-step-convergence improved schemes, the C-FORCE and C-SLIC schemes, are proposed to resolve this problem. The performance of the proposed schemes is validated in solving the one-layer and two-layer shallow-water equations, verifying their capabilities in attaining time-step-independence solutions and showing robustness of them in resolving discontinuities with high-resolution.


Introduction
The Lax-Friedrichs (LF) scheme, also called the Lax method [1], is a classical explicit three-point scheme in solving partial differential equations in, for example, aerodynamics, hydrodynamics, and magnetohydrodynamics [2][3][4].This scheme has many advantages, for example, being stable and monotone up to Courant-Friedrichs-Lewy (CFL) number approaching unity in solving the linear convection problem [4] and being computationally efficient and simple to implement (e.g., compared with a Godunov-type scheme).Due to the excellent monotone and stable properties, the LF scheme, which is of the first-order accuracy in space, has been commonly adopted for a combination use with some higher-order schemes to predict monotone and stable solutions with lower numerical diffusion [5].For instance, Toro and Billett [6] derived the first-order centered (FORCE) scheme, by combining the LF scheme with the two-step, second-order Lax-Wendroff scheme [7].When solving the linear convection problem, compared with the LF scheme, the numerical viscosity of the FORCE scheme is reduced by half [4].The FORCE scheme and its second-order (both in time and in space) extension, the slope-limited centered scheme (SLIC) [8], have been used and verified to perform well in solving the one-layer shallow-water equations (1LSWEs), the Euler equations of compressible gas dynamics, and the magnetohydrodynamic equations [8,9].
Nevertheless, the LF scheme and the schemes devised based on it (e.g., the FORCE and SLIC schemes) have a common shortcoming, which has raised less attention so far.We demonstrate this shortcoming by using the LF scheme to solve the following linear convection problem in the timespace domain [0, ] × [0, ]: Initial condition:  (, 0) =  0 () , Boundary conditions:  (0, ) =  1 () ,  (, ) =  2 () .Here,  is the wave propagation speed and, without loss of generality, is assumed to be positive.An explicit finite-volume discretization for (1a) gives where the superscripts "" and " + 1" denote the values at the old and new time steps, respectively; Δ and Δ denote the time step and space interval, respectively; F denotes the numerical flux and, for the LF scheme, it reads [1] A truncation error analysis shows that the LF scheme in solving (1a), (1b), and (1c), with the discretizations given in (2) and (3), leads to a numerical viscosity of [4] where   is the CFL number and, for the linear convection problem, In Figure 1,  LF is plotted against   .For clear view purpose, the relation in the region   ⩽ 0.01 is plotted on a semilog scale (see Figure 1(a)).For convenience, in Figure 1,  = Δ = 1 is assumed so that   = Δ (see (5)).As seen in Figure 1,  LF increases as Δ decreases.When Δ → 0,  LF → ∞, suggesting that the numerical diffusion is absolutely dominant over the physical convection; this is, of course, unrealistic.In numerical practice, just as usually done via implementing a grid-independence study to choose an appropriate spatial resolution (can be done for the LF scheme as when Δ → 0,  LF → 0), a Δ-independence study is required to choose a reasonable time step.However, due to the rapidly increased numerical diffusion when Δ decreases, the solution computed by using the LF scheme does not converge and thus a Δ-independence study cannot be implemented.Failing to achieve a Δ-independence solution introduces uncertainty in a numerical solution, and the time step can only be determined empirically.Note that, in some situations, for instance, when rapid wet-dry transitions exist (e.g., for simulations in the nearshore regions), or when some other constraints on the time step (e.g., those due to surface tension and fluid viscosity effects [10]) dominate the CFL condition, one may have to use a rather small time step compared with that determined by the CFL condition.In these situations, the LF-based schemes are inappropriate due to the excessive numerical diffusion at a small value of Δ.
We observe that the problem, related to the LF scheme being unable to achieve a Δ-independence solution, can be resolved by the local LF (LLF hereafter) scheme [11].The LLF scheme calculates the numerical flux as where λ denotes the local bound of the wave propagation speed and is determined by and, in solving (1a), (1b), and (1c), λ−1/2 ≡ .In the Appendix, we show that the LLF scheme can be viewed as a firstorder accurate characteristics scheme when solving the linear convection problem.Inserting ( 6) and ( 7) into (2), one attains the following when solving the linear convection problem: which introduces a numerical viscosity of From ( 9) (also see Figure 1), it is seen that  LLF increases with decreasing   .However,  LLF has a maximum and limiting value of 0.5Δ when   → 0. Therefore, the uncontrolled excessive numerical viscosity in the LF scheme is resolved and a Δ-independence study can be performed.
It is worth remarking that, besides the advantage in achieving a Δ-independence solution, the LLF scheme also has the merits (e.g., being stable and monotone) as the LF scheme.
To date, the LLF scheme has been widely used to solve the aerodynamic and magnetohydrodynamic problems [2,3,11] but has not raised much attention in computational hydrodynamics (e.g., in solving the shallow-water equations).The motivation of this work is simple.As the LF-based schemes (e.g., the FORCE and SLIC schemes) share the same drawback as the LF scheme to be unable to achieve a Δindependence solution (note that the numerical viscosity of the FORCE scheme is just half of that defined in (4)), while the LLF scheme overcomes this drawback, it is meaningful to improve the LF-based schemes by a substitution of the LF scheme with the LLF scheme.In this work, the FORCE and SLIC schemes are improved following this idea.For convenience, the modified FORCE (SLIC) scheme will be named as C-FORCE (C-SLIC) scheme hereafter, where the prefix "C" denotes convergence ensured.For illustration purpose and as representatives, the performance of the C-FORCE and C-SLIC schemes will be assessed in solving the 1LSWEs and the two-layer shallow-water equations (2LSWEs), in 1D space and without considering the interfacial shear stress (for the 2LSWEs) and bed friction.
The rest of the paper is organized as follows.In Section 2, a comparison of the performances between the C-FORCE (C-SLIC) scheme and the original FORCE (SLIC) scheme is made in solving the 1LSWEs, and a similar comparison in solving the 2LSWEs is presented in Section 3. Finally, conclusions and discussions are given in Section 4.

Performance Comparisons in Solving the 1LSWEs
We consider solving the 1LSWEs in this section.A sketch of the one-layer shallow-water system is illustrated in Figure 2(a).

Governing Equations and Numerical
Methods for the 1LSWEs.The 1LSWEs, in a conservative vector form, can be written as with vectors defined by where  denotes time;  denotes the streamwise coordinate;  and   denote the water surface and bed elevations above a horizontal reference level   ; ℎ =  −   denotes the water depth; and   = ℎ denotes the discharge per unit width with  being the depth-averaged velocity in the -direction.
An explicit finite-volume discretization for (10) gives Here, the overbar on a term denotes that this term is evaluated by the temporally evolved values (described later in this section); the source term vector is discretized by using a central difference scheme as [12] the flux vector F is evaluated by the FORCE or SLIC scheme, which reads uniformly where Here, the superscripts "LW" and "LF" denote that the variables are related to the two-step Lax-Wendroff and Lax-Friedrichs fluxes, respectively; the superscripts "L" and "R" denote the reconstructed values at the left and right sides of the cell-interface, respectively.For both the FORCE and SLIC schemes,  = 0.5 and λ ≡ Δ/Δ; for both the C-FORCE and C-SLIC schemes,  is set to be 0.3 empirically and λ is set to be the maximum local characteristic speed as For face value reconstructions, first, the monotonic upstream-centered scheme for conservation laws (MUSCL) [13] is applied; for example, where  denotes the ratio of successive gradients (e.g.,   = (  +1 −    )/(   −   −1 )) and Ψ is the flux limiter.For both the FORCE and C-FORCE schemes, Ψ ≡ 0, and, for both the SLIC and C-SLIC schemes, a minmod slope limiter [14] is employed to ensure numerical stability.Second, the hydrostatic reconstruction approach proposed in [12] is employed to achieve the well-balanced property [15] and nonnegative water depth solutions; in this approach, a singlevalued bed elevation at the cell-interface is defined as Third, the reconstructed face values are temporally evolved as For both the FORCE and C-FORCE schemes, ΔU ≡ 0, and, for both the SLIC and C-SLIC schemes, ΔU is calculated according to Here, the source term vector S   is evaluated similarly to (13) but calculated with the values before being temporally evolved (i.e., U L,R −1/2 instead of U L,R −1/2 ).

Numerical
Test for the 1LSWEs.The benchmark dambreak problem is considered here, for which the analytic solution is available [16].The simulation domain is 10 m long ( ∈ [0, 10 m]).The initial flow is at rest and the water surface elevation is initially prescribed as Open boundary conditions (i.e., zero gradients of variables at boundaries) are imposed at both ends of the simulation domain.The simulation domain is discretized by equally spaced 1600 mesh cells, determined based on a gridindependence study using the FORCE scheme with   = 0.5.
Here,   is defined as where λ−1/2 is calculated by (17); the minimum is defined over the whole simulation domain.Note that, for this particular test,   is commonly set to be a large value (e.g.,   ≈ 0.5) in literatures [8] and thus the numerical diffusion in the predictions is insignificant.We choose this test here as it is good to demonstrate the different performances of the LF-and LLF-based numerical schemes in achieving Δindependence solutions.
In Figures 3(a) and 3(b), the computed water surface profiles at  = 0.5 s by using the FORCE and C-FORCE schemes with different values of   are displayed.Also shown is the analytic solution of Stoker [16].As seen, for the FORCE scheme, as   decreases, the numerical diffusion increases, causing a smearing of the wave profile; when   ⩽ 0.05, the numerical diffusion becomes significantly large and gradually dominates other processes.It is surprising that, compared with the analytic solution, a larger value of   contributes to a more accurate solution when the FORCE scheme is employed for flux evaluations.However, when the C-FORCE scheme is used for flux evaluations, as can be seen in Figure 3(b), a Δ-independence solution is successfully achieved, which matches the analytic solution with reasonable accuracy.
In Figures 4(a) and 4(b), the computed water surface profiles at  = 0.5 s by using the SLIC and C-SLIC schemes with different values of   are illustrated.As seen, the solution does not converge when the SLIC scheme is used for flux evaluations, though the numerical diffusion is significantly reduced compared with that simulated by the FORCE scheme (cf.Figures 3(a) and 4(a)).The C-SLIC scheme successfully predicts a Δ-independence solution which matches the analytic solution almost exactly (see Figure 4(b)).

Performance Comparisons in Solving the 2LSWEs
In this section, the proposed scheme in solving the 2LSWEs is considered.A sketch of the two-layer shallow-water system is illustrated in Figure 2(b).

Governing Equations and Numerical
Methods for the 2LSWEs.The 2LSWEs are generally written as (10) with vectors defined by [17] where ℎ  ( = 1, 2) denote the layer depths, with  = 1 and  = 2, respectively, referring to the lower-and upper-layers (see Figure 2(b));   is the bed elevation above a reference level   ;   ( = 1, 2) denote the layer surface elevations with respect to   and have the relations  2 =  1 + ℎ 2 and  1 =   + ℎ 1 ;   ( = 1, 2) denote the layer discharges per unit width in the streamwise direction;   =   /ℎ  ( = 1, 2) denote the layer-averaged velocities; and   ( = 1, 2) are the layer densities and  =  2 / 1 is the density ratio.−ℎ 1 (ℎ 2 /) in ( 24) is a nonconservative product term, which has no definition around shocks and discontinuities and is tough to process in numerical discretizations.Following Spinewine et al. [18], and Lu et al. [17], the 2LSWEs are reformulated to make the nonconservative product term vanish.This is done by replacing the equations governing the lower-layer flow motions by a combination of the equations for the two-layer flow system as a whole; namely, the 2LSWEs are recast to (10) with vectors defined by where As the 1LSWEs and 2LSWEs are written in the same conservative form (see (10)), the 2LSWEs are solved similarly as that for the 1LSWEs presented in Section 2.1, with a difference in calculating λ in (17).In solving the 2LSWEs, for both the FORCE and SLIC schemes, λ ≡ Δ/Δ, and, for both the C-FORCE and C-SLIC schemes, λ is set to be the approximately maximum local characteristic speed as [17,18] where ] . (28b)

Numerical
Test for the 2LSWEs.The internal dam-break problem is considered here, which is widely adopted to test  the robustness and accuracy of numerical schemes in solving the 2LSWEs [19][20][21].The configuration of this test is the same as that in the one-layer dam-break test presented in Section 2.2, but with an upper-layer superimposed on the lower-layer.The density ratio between the layers is set to  = 0.7.The initial condition for the lower-layer flow is the same as that defined in (22), and the initial surface level of the upper-layer is  2 (,  = 0) = 2 m.The spatial resolution is set to be the same as that in the one-layer dam-break problem test.
As the analytic solution for this test is unavailable, the numerical result computed by using the C-SLIC scheme with   = 0.5 and   = 6400 is used as a reference solution.Here,   is defined by (23) with λ−1/2 calculated by (27);   denotes the number of the mesh cells.In Figure 5, the interface profile of the reference solution at  = 1 s is shown, along with the numerical results computed by Bouchut and de Luna [19] using a time-splitting method to decouple the 2LSWEs into two 1LSWEs with an entropy inequality for controlling the numerical instability (the relevant scheme is abbreviated as BD08 scheme hereafter), by Bouchut and Zeitlin [21] using a source-centered hydrostatic reconstruction scheme (BZ10 scheme), and by Dudzinski and Lukáčová-Medvid' ová [20] using a bicharacteristics theory based finite-volume evolution Galerkin scheme (DL13 scheme).A zoom view in the region  ∈ [4.9 m, 5.7 m] is given on the bottom-right of Figure 5.As seen, spurious oscillations can be easily predicted in simulations for this test.The solution computed by the C-SLIC scheme matches the predictions computed by the BD08 and DL13 schemes well, but without oscillations (see the zoom view in Figure 5).
In Figures 6(a) and 6(b), the computed layer surface profiles at  = 1.0 s by using the FORCE and C-FORCE schemes with different values of   are shown.The relevant results for the SLIC and C-SLIC schemes are presented in Figure 7.Note that, in Figures 6 and 7, the reference solution is obtained by using the C-SLIC scheme with   = 0.5 and   = 6400.As we compare Figure 3 with Figure 6 and Figure 4 with Figure 7, it is found that the effects of   on the results of the two-layer internal dam-break problem and the one-layer dam-break problem are quite similar: the numerical predictions do not converge when the FORCE or SLIC scheme is used, while a Δ-independence solution is achieved by using the C-FORCE or C-SLIC scheme.

Conclusions and Discussions
In this work, we have analyzed the problem of the LF-based schemes (e.g., the FORCE and SLIC schemes) regarding achieving time-step-independence solutions; it is found that a time-step-independence study cannot be performed for these schemes as the numerical solutions do not converge with successive finer time steps, owing to the excessive numerical diffusion at a small time step.
Two convergence improved schemes, the C-FORCE and C-SLIC schemes, have been proposed to overcome the relevant problem related to the LF-based schemes.As representatives, the C-FORCE and C-SLIC schemes have been applied to solve the one-layer and two-layer shallow-water equations.Results show that the proposed schemes successfully achieve Δ-independence solutions.Also, the proposed schemes are verified to be capable of resolving sharp fronts well without oscillations.
The proposed schemes are presented in 1D space in this work.The extension of the proposed schemes to 2D space is straightforward.Also, the use of the C-FORCE and C-SLIC   schemes in solving other hyperbolic equation systems is of great interest and is left for future studies.

Figure 1 :
Figure 1: Numerical viscosity versus   for the LF and LLF schemes in solving the linear convection problem: (a)   ⩽ 0.01 and (b) 0.01 ⩽   ⩽ 1.Note that the result for the LLF scheme is plotted on the right -axis.(a) is plotted on a semilog scale.

Figure 2 :
Figure 2: Sketch of the shallow-water systems: (a) one-layer and (b) two-layer.Not to scale.

Figure 3 :
Figure 3: Computed water surface profiles by using the FORCE scheme (a) and the C-FORCE scheme (b) at  = 0.5 s, for the one-layer dam-break problem test.The analytic solution of Stoker [16] is also shown.

Figure 4 :
Figure 4: Same as Figure 3, but computed by using the SLIC scheme (a) and C-SLIC scheme (b), respectively.

Figure 5 :
Figure 5: Computed interface profile by using the C-SLIC scheme with   = 0.5 and   = 6400 at  = 1 s, along with those predicted by the BC08, BZ10, and DL13 schemes, for the two-layer internal dam-break problem test.A zoom view is shown on the bottom-right.

Figure 6 :
Figure 6: Computed layer surface profiles by using the FORCE scheme (a) and the C-FORCE scheme (b) at  = 1 s, for the two-layer internal dam-break problem test.The reference solution is obtained by using the C-SLIC scheme with   = 0.5 and   = 6400.

Figure 7 :
Figure 7: Same as Figure 6, but computed by using the SLIC scheme (a) and C-SLIC scheme (b), respectively.

Figure 8 :
Figure 8: LLF scheme discretizations on a time-space domain.