Application of Large Time Step TVD High Order Scheme to Shallow Water Equations

: The numerical modeling of actual river floods faces three challenges related to computational efficiency, accuracy, and proper balancing of terms in the governing equations, all of which are discussed in this paper. Herein, a large time step (LTS) scheme is used to improve efficiency, a high order scheme is used to enhance accuracy, and specific treatment of the bed slope term achieves a well-balanced form of shallow water equations. The LTS scheme, originally proposed by LeVeque in 1998, has led to the development of highly efficient computational solvers of the shallow water equations (SWEs). This paper examines use of a total variation diminishing (TVD) high order scheme in conjunction with LTS. We first applied the scheme to the solution of the homogeneous 1D SWEs and obtained satisfactory results for three cases, even though small oscillations neverthe-less occur when the CFL number is very large. The additional source term makes the issue more complicated and can introduce a spurious flow when the term is not correctly handled. Many meth-ods have been developed in traditional differential schemes, but not all are fit for the TVD-LTS scheme; for example, the method of decomposing the source term into simple characteristic waves has proved infeasible. In this paper the TVD-LTS scheme was implemented for the first time for well-balanced SWEs containing bed slope source terms. We found that oscillations were not as suppressed as for the homogeneous SWEs when the TVD-LTS scheme was applied to the three cases of step Riemann problems (SRP) tested for CFL numbers 1 to 10. For free surface flow over a bed hump, the TVD-LTS scheme can only reach CFL number 4 before the solution breaks down.


Introduction
Shallow water equations (SWEs) are depth-integrated continuity and momentum equations for river flow, in which it is assumed that surface waves are long, the depth is relatively small in comparison with the horizontal scale of the flow domain, and the pressure is hydrostatic. In the version considered in the present paper, viscous terms are neglected, which restricts the applicability of the equations but is acceptable for a wide range of flows in open channel hydraulics (see, e.g., Abbott and Minns [1] and Toro [2]). The SWEs can be formulated as a nonlinear hyperbolic system that can be discretized using an upwind method such as the Roe [3], Osher [4], and AUSM [5] schemes. Alternative higher-order schemes such as WENO [6,7] also work by reconstruction, using additional cells next to the interface. However, all the foregoing schemes are constrained by a time step stability condition, such as the Courant-Frederichs-Lewy number CFL < 1, where CFL = umaxΔt/Δx, in which umax is the maximum flow speed, Δt is the time step, and Δx is a spatial increment in the computational mesh. In the 1980s, LeVeque [8,9] proposed the large time step (LTS) method, which led to considerable gains in computational efficiency. In practice, LTS does not conflict with CFL theory in that the interface flux is computed across multiple cells in LTS rather than a single cell in traditional schemes, thus enabling the time step to be increased accordingly. In LTS, more cells take part directly in the computation of fluxes, unlike higher-order schemes that use reconstruction schemes to handle cells adjacent to the interface. Murillo et al. [10] and Morales-Hernández [11][12][13] then applied LTS to solve the SWEs for free surface flows, and similarly, Qian [14] used LTS in solving the Euler equations in aerodynamics. All the foregoing LTS schemes are first order. However, as the CFL number becomes large, spurious oscillations appear in the solution, which do not occur in traditional first-order schemes. Harten introduced the total variation diminishing (TVD) concept [15] that could suppress oscillations for second-or higher-order schemes in the traditional form. Harten [16] then proposed a TVD-LTS scheme that was further developed by Qian and Lee [17].
Godunov-type solvers of the SWEs can generate non-physical flow phenomena from discretized bed slope source terms that are out of balance with the flux gradient terms. Well-balanced schemes overcome such difficulties. LeVeque [18] proposed a quasi-steady wave propagation method that incorporated source terms within high-resolution Godunov methods and obtained accurate results for the propagation of small perturbations. Hubbard and García-Navarro [19] utilized a flux difference splitting method to ensure the flux gradient and source terms were correctly balanced. Zhou et al. [20] suggested a surface gradient method to represent the bed slope source term. Rogers et al. [21] presented an algebraic technique for balancing the flux gradient and source terms for single waterbodies. Liang et al. [22] extended the algebraic technique to multiple waterbodies, with stage and discharge selected as dependent variables.
This paper applies TVD-LTS to solving the SWEs without viscous stresses and expands the scheme to non-homogeneous hyperbolic systems for the first time, to the authors' knowledge. The results indicate that TVD-LTS gives a satisfactory representation of benchmark flows while suppressing spurious oscillations when solving homogeneous SWEs but is not as successful when applied to non-homogeneous SWEs. The paper has the following structure. Section 2 outlines the TVD-LTS solver of the shallow water equations. Section 3 presents results for different benchmark tests: homogeneous cases of two rarefactions; two shock waves; a single rarefaction and a shock wave; non-homogeneous cases of two rarefactions over a step; two shockwaves over a step; a single rarefaction and a bore over a step; and a transcritical flow over a bed hump. Section 4 provides the concluding remarks.

TVD-LTS Scheme for Scalar Case
Consider the homogeneous hyperbolic equation: Harten [16] proposed the following second order TVD-LTS scheme: ( ) Now, we consider the non-homogeneous case: The equation is discretized with the TVD-LTS method: According to Harten's TVD-LTS scheme, (3) is a Lipschitz continuous function, and so: Substituting into Equation (3), and we obtain: Equation (17) becomes: The first part of Equation (14) is second order. The remaining part, is also second order in space, and so Equation (14) is entirely second order in space. The source term is solely a function of space and will not change in time, and so it will not change the property of total variation diminishing, as verified by Harten [16]. Then, it is stable.

TVD-LTS Scheme for Shallow Water Equations
We consider shallow free surface flows in the horizontal dimension and assume hydrostatic pressure. The resulting inviscid shallow water equations (see, e.g., [2]) in the general matrix-vector form are expressed as: For the homogeneous condition (with no bed slope), the vectors of dependent variables, horizontal fluxes, and source terms are: For the non-homogeneous condition with a finite bed slope, the vectors become ( [22]): In Equations (21)-(23), U is the vector of stage and discharge, F(U) is the vector of horizontal fluxes, Sf is the vector of source terms, x is distance, t is time, h is the local water depth, u is mean flow velocity, g is the acceleration due to gravity, η is the stage (i.e., water free surface elevation above a fixed horizontal datum), and z is the bed elevation above the same datum.
In Roe's method, the dependent variables at adjacent cells j and j + 1 are averaged as: Equation (21) describes a nonlinear hyperbolic system for which the Jacobi matrix can be written: After diagonalization: where the characteristic eigenvalues are: The vector difference in U between adjacent cells is transferred to the characteristic space, such that: Equation (21) is then differenced using an explicit upwind scheme to give: Following Harten [16], the numerical flux in TVD-LTS is: Equations (7)- (9) ensure that each grid cell is correctly allocated to the right interface waves ( 0 k  ) when they spread to the left ( 0   ) and to the left interface waves ( 0 k  ) when they spread to the right ( 0   ).

Homogeneous SWEs
The numerical model is tested against three benchmark cases involving different combinations of rarefaction and shock waves propagating in a frictionless, straight, 25 m long one-dimensional channel with a fixed horizontal bed, for which exact solutions are provided by Toro [2]. The domain is discretized into 250 cells, such that ∆x = 0.1 m. The adaptive time step is determined from ∆t = K ∆x/λmax, where K is the CFL number, and λmax is the maximum eigenvalue in Equation (27). The initial time step is 0.01 s, after which its value depends on the instantaneous value of λmax. Throughout the computations, the local flow depth and velocity are updated using Equation (30), in which the final term is set to zero for the flat bed.   Figure 1 presents the distributions of the water free surface and discharge per unit width along the channel for different CFL numbers using the present TVD-LTS model. In Figure 1, the upper-left plots display the distributions obtained for a maximum CFL of unity at time t = 0.70158 s. The predicted water free surface profiles provide a satisfactory match to the solution by Toro's method [2] and Roe's scheme [3], with the head of the leftward rarefaction wave occurring at x = 4.0 m and its tail at x = 12 m and the head of the rightward rarefaction wave at x = 16 m and its tail at x = 22.4 m. The computed and exact mass flux distributions are very similar, except for the slight deviation at the interface between the flow plateau and the rarefaction.

Two Shockwaves in a Horizontal Channel
Next, we consider an idealized hydraulic bore caused by a tide interacting with an opposing river flow in the horizontal channel. The initial conditions listed in Table 2 Figure 2 shows the water free surface and mass flux distributions along the channel, predicted by the TVD-LTS model for different CFL numbers. There is generally close agreement with Toro's [2] exact solution. In Figure 2, the profiles correspond to a rightdirected hydraulic bore at x = 17.2 m and a left-directed hydraulic bore at x = 10.6 m. The upper-left plots display the results obtained for a maximum CFL value of unity.
Oscillations in the TVD-LTS solution become increasingly evident as the CFL number becomes larger.

One Rarefaction and One Shockwave in a Horizontal Channel
The third case involves a dam break over a wet bed that is initially quiescent. The initial conditions are summarized in Table 3. Table 3. Initial conditions for a homogeneous case of one rarefaction and one hydraulic bore in a frictionless, horizontal channel.  Figure 3 presents the computed water free surface and mass flux distributions along the channel for different CFL numbers. Here, the exact solution obtained by Toro [2] consists of a right-directed hydraulic bore at x = 20.2 m and a left-directed rarefaction wave, whose head is at x = 2.4 m and tail is at x = 9.9 m. Again, we can see that the predicted distributions provide a reasonable match to the exact solution when CFL is unity. As the CFL number increases, slight oscillations become noticeable at the location between the flat free surface and the shock wave.

Non-Homogeneous SWEs
The next three cases concern the propagation of rarefactions and shocks over an abrupt step in the bed (i.e., a shelf), where the SWEs become non-homogeneous owing to the presence of bed slope terms. In the Godunov-type FVM scheme, each interface between cells is considered as a discontinuity that can only be solved by a Riemann solver, and so these cases (such as the bed step) are benchmark foundations from which the method can be extended to more practical problems. For example, local discontinuities are likely to occur when a continuous river bed is discretized. Exact solutions for the three benchmark cases are given by Rosatti and Begnudelli [23]. The final case considers a transcritical flow at a bed hump, again where the SWEs are non-homogeneous, for which an exact solution is given by Goutal and Maurel [24]. In all cases, the channel has an overall length of 25 m, a constant spatial increment of ∆x = 0.1 m, and an initial time step of 0.01 s (which is altered at every time step thereafter, according to ∆t = K ∆x/λmax). The flow depth and velocity are updated using Equation (30), but with the final term no longer zero at the step.

Two Rarefaction Waves in a Channel Containing a Bed Step
The first test case considers flow in a channel containing a 1 m high abrupt step, with the bed elevation defined as z = 0 for x < 12.5 m and z = 1 m for x ≥ 12.5 m. For diverging flow, the initial depth either side of the shelf is 8.0 and 5.0 m, with corresponding opposing velocities of −2 and 7.1704 m/s (Table 4).  Figure 4 shows the computed and exact solutions for the water free surface and discharge per unit width distributions, obtained using different CFL numbers. In all cases, the free surface profiles consistently contain a left-directed rarefaction wave whose head is at x = 3.1 m and tail is at x = 10 m. A stationary shock occurs at the step. Additionally, a right-directed rarefaction wave is evident, with its head at x = 19 m and its tail at x = 22.4 m. The computed and exact mass flux distributions are in satisfactory agreement, except for a slight oscillation at the interface of the step (x = 12.5 m) in the mass flux from CFL number 1 to CFL number 10.

Two Shockwaves in a Channel Containing a Bed Step
The second test case considers an idealized upstream-propagating tidal bore interacting with a downstream-propagating river flow in the same channel as in the previous case. The upstream depth to the left of the bore is 4.0 m, and the velocity is 4.75 m/s. The downstream depth to the right of the bore is 1.0838 m on top of a 1 m elevated riverbed so that the stage is 2.0838 m and the flow velocity is −2.1854 m/s. Here, the opposing bore and river flows approach the middle interface in opposite directions (Table 5).   Figure 5 shows the computed and analytical solutions for the water free surface and mass flux distributions along the channel, obtained for different CFL numbers. The solution again includes a stationary shock at the front edge of the shelf but is different from the preceding case in that the waves emanating from either side are no longer rarefactions but, instead, are two bores at x = 8 m and x = 21 m, traveling in opposite directions. In the TVD-LTS predictions, tiny oscillations in the mass flux can be discerned at the step (x = 12.5 m) even when the CFL number is only 1. At higher CFL numbers, the oscillations extend progressively to the right of the step (x > 12.5 m), with both the free surface and mass flux exhibiting increasingly intense fluctuations. Oscillations of relatively smaller amplitude also occur to the left of the step (x < 12.5 m).

One Rarefaction and One Shockwave in a Channel Containing a Bed Step
Unlike the previous two cases, this case involves zero flow velocity on either side of the step interface, as indicated in Table 6. The initial depths are 5.0 and 1.0 m to the left and right of the step. Table 6. Initial conditions for a rarefaction and a hydraulic bore propagating over a step.

h (m)
q (m 2 /s) z (m) Left side 5.0 0 0 Right side 1.0 0 1 Figure 6 presents water free surface and mass flux distributions along the channel obtained for different CFL numbers. In this case, the solution consists of a hydraulic drop at the front face of the step, a right-directed hydraulic bore at x = 20.3 m, and a left-directed rarefaction wave with its head at x = 2.8 m and its tail at x = 8.3 m. The predicted results to the left of the step face (i.e., x < 12.5 m) are better simulated than to the right of the step face, especially for larger values of the CFL number, with small amplitude oscillations occurring at the interface between the rarefaction and the flow plateau. The predicted mass flux at the step face (x = 12.5 m) is lower than the analytical result, even when the CFL number is only 1. Decaying fluctuations in the predicted results occur to the right of the front face of the step (x > 12.5 m), with the amplitude increasing with the CFL number. At the step face itself, the deviation between the predicted and analytical mass flux is at a maximum, independent of the CFL number considered (in the range from 1 to 10).

Transcritical Flow at a Fixed Bed Hump
Goutal and Maurel [24] derived the analytical solution for shallow flow over a hump in the riverbed where the bed level satisfies the following equation: In this case, the channel is 25 m long and carries an upstream inlet discharge of 0.18 m 3 /s. The downstream outlet water level is prescribed to be 0.33 m. Initially, the water level is 0.33 m (see Figure 7). Figure 8 shows the analytical solution at a steady state. In the numerical model, the channel is divided into a regular grid of 250 cells, each with a spatial increment of x = 0.1 m. The initial time step is prescribed a value of 0.01 s, which alters each time step thereafter, according to ∆t = K ∆x/λmax. Figure 9 compares the numerical predictions of the water level and mass flux profiles obtained using Equation (30) against the analytical solution. Reasonable agreement is obtained between the predicted and analytical profiles when the CFL number is 1, 2, 3, and 4, but a converged solution cannot be obtained when the CFL number is larger than 4. The over-prediction of mass flux arises for CFL numbers 1 and 2. Over-and under-predictions are both observed for CFL numbers 3 and 4.    To investigate the efficiency of the TVD-LTS algorithm, the computational time overhead was determined. All computations were undertaken on the same computer, which was configured as follows: Intel Core i9-10900K @ 3.70 GHz with an overclock CPU; Kingston DDR4 128 GB 3200 MHz with XMP opened memory; and Intel 750 series PCIe interface 1.2 TB SSD (also used by the authors in their previous work [25]). Table 7 summarizes the computation effort in terms of the number of computational steps, CPU time, and average run time per computational step (determined as the ratio of CPU time to the number of steps) required for convergence to a steady state (difference of U between two neighboring time steps less than 10 −5 ) for transcritical flow at a hump according to CFL number. It should be noted that the use of the same algorithm on the above computer could lead to slightly different CPU time evaluations because of the auto-arrangement of the CPU load and the use of random memory blocks. In all the tests considered, the discrepancy is in the range of ±0.2 s. The average run time per computational step is 5.73 × 10 −3 s for a CFL number equal to 1. The run time increases to 31.64 × 10 3 s when the CFL number is raised to 4. Unlike the first-order LTS scheme, the TVD-LTS accounts for information from all relevant cells even when their contribution is very small, thus complicating the coefficient (Equation (15)) and leading to an abrupt increase in calculation effort.

Conclusions
Total variation diminishing (TVD) high order large time step schemes have been described for solving homogeneous and non-homogeneous versions of shallow water equations (SWEs), the latter containing bed slope source terms. The resulting model was tested extensively against standard benchmark tests for which standard solutions are available: three tests concerning discontinuous flow in a horizontal channel, three tests of discontinuous flow in a channel with a bed shelf, and one test of transcritical flow at a bed hump. In the first six cases, spurious oscillations develop as the CFL number increases and are most pronounced immediately downstream of the step when it is present. Rarefactions are properly simulated in all cases except for the under-prediction of mass flux at the step interface (x = 12.5 m). In general, the results for the homogeneous SWEs are better predicted than for the non-homogenous SWEs. For transcritical flow at a hump, satisfactory predictions were achieved for CFL numbers up to 4, beyond which the mass flux became overwhelmed by oscillations. The CPU time per simulation time step was found to be much larger than for the conventional first-order Roe scheme, and it increased progressively with increasing CFL number because information from more grid cells had to be accommodated at each time step. In practice, we found that the TVD-LTS scheme required fewer steps but at the cost of greater total computational time, indicating that high order LTS is not economical from a computational perspective.
Author Contributions: Conceptualization, R.X.; methodology, R.X.; writing-original draft preparation, R.X.; writing-review and editing, A.G.L.B. and B.X. All authors have read and agreed to the published version of the manuscript.
Funding: This research was funded by the National Natural Science Foundation of China (51509216, 52079120).

Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Not applicable.