A Novel Multislope MUSCL Scheme for Solving 2D Shallow Water Equations on Unstructured Grids

Abstract

Within the framework of the two-dimensional cell-centered Godunov-type finite volume (CCFV) method, this paper presents a novel multislope scheme on the basis of the monotone upstream scheme for conservation law (MUSCL) for numerically solving nonlinear shallow water equations on two-dimensional triangular grids. The Riemann states of the considered edge are calculated by an edge-based reconstructing procedure, where a limited scalar slope is employed to prevent potential numerical oscillations. The novel aspect of the new scheme is that it takes advantage of the geometrical characteristics of triangular grids in the reconstructing and limiting procedures, which effectively reduces the cost of computation and provides higher resolution and accuracy compared with classical MUSCL schemes. Seven tests are adopted to verify the scheme, and the results indicate that this scheme is efficient, accurate, robust, and high-resolution, and can be an ideal alternative for solving shallow water problems over uneven and frictional topography.

1. Introduction

Two-dimensional shallow water equations (SWEs) have been increasingly employed to simulate predominantly horizontal, free surface flows over natural topography, for instance, dam-break flow, urban flooding, and tsunami inundation. Analytical solutions of SWEs are only available for a few simplified situations with horizontal or regular bottom and infinite water volume [1,2]. Turning to numerical schemes is necessary in order to keep the balance between accuracy and efficiency.

SWEs are a set of nonlinear hyperbolic-type partial differential equations, and several kinds of methods have been developed to solve them, such as the finite difference method (FDM) [3], the finite element method (FEM) [4], the finite volume method (FVM) [5], and the discontinuous Galerkin (DG) method [6]. The primary advantages of the FVM over the other methods stem from its geometrical flexibility for space discretization and its way of defining flow variables and discretizing convective fluxes [7]. The Godunov-type scheme [8] within the framework of FVM, taking into account the propagation of flow information in resolving a local Riemann problem [9], is probably the most applied method of simulating free surface flows, for its inherent conservation property and stability in mixed flow regimes [10]. There are basically two approaches to derive FVM schemes, the cell-centered finite volume (CCFV) approach and the node-centered finite volume (NCFV) approach. A distinct difference between these two approaches is the relative position between the control volume and the mesh element. Besides, CCFV schemes store solutions at cell centroids, while NCFV schemes store solutions at nodes. In addition, the mixed CCFV and NCFV discretization method, which defines different flow variables separately at centroids with cells as control volumes, or at nodes with the control volumes constructed by the median-dual partition [7,11] with grids, has been proposed and successfully implemented to simulate the bank failure phenomenon [12,13]. This method integrates characteristics of the CCFV and NCFV methods, allowing the slope to be univocally defined in each cell and enabling a straightforward inclusion of the geo-failure operator in the morphodynamic model. The Godunov-type CCFV scheme has been a widely used numerical method for solving SWEs in recent decades [7,14,15,16,17,18,19,20], for its shock-capturing nature and desirable property of automatic conservation. Unfortunately, the conventional Godunov-type CCFV scheme provides only first-order accuracy in space that exhibits relatively large numerical diffusion and follows a low-resolution property. Taking a linear monotone upstream scheme for conservation law (MUSCL) [21] as a means of reconstruction, second-order spatial accuracy can be easily achieved. The Godunov-type CCFV scheme calculates the left and right flow states at a considered edge, assuming a piecewise constant distribution of flow variables in each cell, whereas the MUSCL scheme assumes a piecewise linear function of flow variables in each cell.

The classical MUSCL scheme consists of two steps to achieve the reconstructed values of flow variables. First, neighboring values of the considered cell are used to calculate a predicted slope. Then the slope is limited and gives rise to a vectorial slope, which is used to calculate the values of flow variables on the left and right sides of the considered edge [22]. The MUSCL scheme without stability constraint tends to produce unphysical spurious oscillations near discontinuities due to their unbounded reconstruction procedures. A measurement of local monotonicity, also known as the r-ratio [7,23,24], is usually employed as an argument to limit the slopes and provides oscillation-free results. To this end, combining second-order schemes with selected limiters, such as the Van Leer [25] or Van Albada [26] limiter, results can be obtained with high accuracy, high resolution, and good convergence.

In addition to the monoslope MUSCL scheme presented above, there is another kind of MUSCL scheme, the multislope MUSCL scheme [22]. The former calculates only one vectorial slope per cell, respecting a specified stability constraint to prevent numerical oscillations, whereas the latter calculates three scalar slopes equipped with three specified stability constraints for each edge per cell. The multislope concept is adopted in this work, and a novel multislope MUSCL scheme that takes the mesh geometry into account is proposed. To achieve second-order accuracy in both space and time, the classical two-stage Runge–Kutta time integration scheme is used to update solutions at the new time level. Explicit descriptions of this time integration scheme can be found in [10,27].

The primary purpose of this work is to present a new multislope MUSCL scheme on 2D triangular grids, where a well-designed scalar slope is equipped with an edge-based limiter for each edge. This new scheme is incorporated into the Godunov-type CCFV method to solve SWEs. The rest of this paper is organized as follows: Section 2 introduces the governing equations of SWEs. In particular, a brief description of the eigenstructure of SWEs is stated. Section 3 presents the details of the new multislope MUSCL scheme, including a review of the existing MUSCL-type schemes. In Section 4, some numerical results and comparison analyses over horizontal or uneven topography are provided. A brief conclusion is drawn in Section 5.

2. Governing Equations

When the vertical velocity and displacement of water surface is sufficiently small compared with the frequency of water waves (Figure 1), the governing equations of shallow water flows can be obtained by integrating the Navier–Stokes equation over depth by assuming a hydrostatic pressure distribution, as well as a uniform distribution of the horizontal velocity in the vertical direction and a small bottom slope. Neglecting the diffusion term due to viscosity and turbulence, the Coriolis term, and the surface shear stress term, SWEs can be expressed in a vector form:

$$\frac{\partial \mathbf{U}}{\partial t}+\frac{\partial \mathbf{E}}{\partial x}+\frac{\partial \mathbf{G}}{\partial y}=\mathbf{S}$$

(1)

where $t$ denotes time; $x$ and $y$ are Cartesian coordinates in the horizontal plane; $\mathbf{U}$ represents the vector of conserved variables; $\mathbf{E}$ and $\mathbf{G}$ are the vectors of convective fluxes in $x$ and $y$ directions, respectively; $\mathbf{F}={(\mathbf{E},\mathbf{G})}^{\text{T}}$ is the flux tensor; and $S$ denotes the source term, which can be further divided into the bed slope source term $S_b$ and friction source term $S_f$. These vectors are expressed as:

$$\begin{gathered} U=\left[\begin{array}{c}h\\ hu\\ hv\end{array}\right], E=\left[\begin{array}{c}hu\\ h{u}^2+g{h}^2/2\\ huv\end{array}\right], G=\left[\begin{array}{c}hv\\ huv\\ h{v}^2+g{h}^2/2\end{array}\right], \\ S=S_b+S_f=\left[\begin{array}{c}0\\ -gh\partial z_b/\partial x\\ -gh\partial z_b/\partial y\end{array}\right]+\left[\begin{array}{c}0\\ -C_{f}u\sqrt{u^2+v^2}\\ -C_{f}v\sqrt{u^2+v^2}\end{array}\right] \end{gathered}$$

(2)

where $h$ represents water depth; $u$ and $v$ are depth-averaged velocities in $x$ and $y$ directions, respectively; $g$ is gravitational acceleration; $z_b$ is the bottom elevation relative to the datum plane, as demonstrated in Figure 1; and $C_f$ is the bed roughness coefficient evaluated by:

$$C_f=g{n}^2/h^{1/3}$$

(3)

where n denotes the Manning coefficient.

Nontrivial steady solutions to SWEs are physically prevalent in shallow water phenomena as the time partial derivative term in Equation (1) vanishes and the SWE is deduced to $\nabla \cdot \mathbf{F}=\mathbf{S}$. One of the most important and interesting steady solutions is the “lake at rest” one:

$$u=v=0,h+z_b=Const.$$

(4)

If the water surface remains constant over the domain and the flow velocities vanish, the bed source term and the flow depth gradient term in Equation (2) reach balance automatically. However, it is not always the case, because these two terms are discretized separately. Numerical difficulties may arise in situations when they should exactly balance with each other [28].

The flux Jacobian matrix of SWEs (Equations (1) and (2)) is:

$$J_n=\frac{\partial \mathbf{E}}{\partial \mathbf{U}}{n}_x+\frac{\partial \mathbf{G}}{\partial \mathbf{U}}{n}_y=\left[\begin{array}{ccc}0& n_x& n_y\\ gh{n}_x-u(u{n}_x+v{n}_y)& 2u{n}_x+v{n}_y& u{n}_y\\ gh{n}_y-v(u{n}_x+v{n}_y)& v{n}_x& u{n}_x+2v{n}_y\end{array}\right]$$

(5)

where $n_x$ and $n_y$ are, respectively, two components of the unit vector in the $x$ and $y$ directions.

The eigenstructure of the flux Jacobian matrix is:

$$J_n=R_n{\Lambda }_n{L}_n$$

(6)

where $R_n$ and $L_n$ represent the right and left eigenvector matrices, respectively, and ${\Lambda }_n$ is the diagonal matrix of the eigenvalues of $J_n$, given by:

$${\Lambda }_n=\left[\begin{array}{ccc}{\lambda }_1& 0& 0\\ 0& {\lambda }_2& 0\\ 0& 0& {\lambda }_3\end{array}\right]$$

(7)

where ${\lambda }_1$, ${\lambda }_2$, and ${\lambda }_3$ denote three eigenvalues of $J_n$ and ${\lambda }_1=u{n}_x+v{n}_y-\sqrt{gh}$, ${\lambda }_2=u{n}_x+v{n}_y$, and ${\lambda }_3=u{n}_x+v{n}_y+\sqrt{gh}$. Obviously, these three eigenvalues are all real and distinct from each other, provided the water depth is positive, so SWEs confirm strict hyperbolicity.

3. Numerical Method

3.1. FVM on Triangular Grids

Two-dimensional grids used for numerical schemes range from structured to highly unstructured. Unstructured triangular grids are adopted in this work for their extreme flexibility with complex flow fields. Values of flow variables are stored at cell centroids, and a piecewise constant distribution is assumed over the computation domain, which conforms to the CCFV scheme. As illustrated in Figure 2, $T_i$ denotes the centroid of cell i; $N_1$, $N_2$ and $N_3$ are three nodes anticlockwise around $T_i$; j is one of the three neighbors of cell i sharing a common edge e; and $M_{ij}$ is the midpoint of edge e where a local Riemann problem is constructed.

SWEs (Equations (1) and (2)) can be integrated over cell i:

$${\displaystyle \underset{{\Omega }_i}{\int }\frac{\partial U}{\partial t}}d\Omega +{\displaystyle \underset{{\Omega }_i}{\int }\nabla \cdot \left(E,G\right)}d\Omega ={\displaystyle \underset{{\Omega }_i}{\int }S}d\Omega$$

(8)

where ${\Omega }_i$ is the control volume, as shown in Figure 2.

Applying Gauss’s divergence theory to the spatial derivative term in Equation (8) yields:

$${\mathrm{A}}_{\mathrm{i}}\frac{\partial U_i}{\partial t}+{\displaystyle \underset{{\Gamma }_i}{\oint }(E{n}_x+G{n}_y)}d\Gamma ={\displaystyle \underset{{\Omega }_i}{\int }S}d\Omega$$

(9)

where $A_i$ and ${\Gamma }_i$ denote, respectively, the area and boundary of cell i and ${\mathbf{U}}_i$ represents the cell-averaged values of conserved variables at a given time.

Assuming a uniform distribution of fluxes along each edge, which are equal to the values at the midpoint of the edge, the surface integral in Equation (9) is redefined and a semidiscrete scheme is given by:

$$\frac{\partial U_i}{\partial t}+\frac{1}{{\mathrm{A}}_{\mathrm{i}}}{\displaystyle \sum _{k=1}^N{F}_{ik}(U)}{l}_{ik}=\frac{1}{{\mathrm{A}}_{\mathrm{i}}}{\displaystyle \underset{{\Omega }_i}{\int }S}d\Omega$$

(10)

where N is the number of edges for cell i; $l_{ik}$ is the length of edge k; and $F_{ik}$ is the outward normal flux vector across edge k and can be estimated with the help of Riemann solvers. The expression of $F_{ik}$ is:

$$F_{ik}(U)={\mathbf{E}}_{ik}(U)n_{fx}+{\mathbf{G}}_{ik}(U)n_{fy}=\left[\begin{array}{c}hu{n}_{fx}+hv{n}_{fy}\\ (h{u}^2+g{h}^2/2)n_{fx}+huv{n}_{fy}\\ huv{n}_{fx}+(h{v}^2+g{h}^2/2)n_{fy}\end{array}\right]$$

(11)

where $n_{fx}$ and $n_{fy}$ are two components of the unit outward normal vector at the midpoint of edge f.

3.2. Approximated Riemann Solver

The Harten–Lax–van Leer–contact (HLLC) approximated Riemann solver [29] is used in the construction of numerical fluxes. Other solvers can be adopted as well, but are not presented in this work. The HLLC approximated Riemann solver restores the missing rarefaction wave by some estimates, for instance, using the Roe average velocity as the speed of the middle wave, which is very efficient compared with the exact Riemann solver, while behaving robustly. The numerical flux vector normal to the considered edge is:

$$F_{}^{HLLC}({\mathbf{U}}_L,{\mathbf{U}}_R)=\{\begin{array}{l}{F}_{L}if{S}_L\ge 0\\ F_{*L}if{S}_L<0\le S_{*}\\ F_{*R}if{S}_{*}<0<S_R\\ F_{R}if{S}_R\le 0\end{array}$$

(12)

where ${\mathbf{U}}_L$ and ${\mathbf{U}}_R$ are the Riemann states to the left and right of the considered edge; $F_L$ can be calculated by Equation (11) using ${\mathbf{U}}_L$ as the independent variable, and the similar operation for $F_R$; $S_L$ and $S_R$ are the left and right wave speeds; $S_{*}$ denotes the contact wave speed; and $F_{*L}$ and $F_{*R}$ are the numerical flux vectors to the left and right of $S_{*}$, as shown in Figure 3.

The values of the wave speeds $S_L$, $S_R$ and $S_{*}$ can be evaluated by [30]:

$$\begin{gathered} S_L=\{\begin{array}{ll}{u}_R^{\perp }-2\sqrt{g{h}_R}& if{h}_L=0\\ \min (u_L^{\perp }-\sqrt{g{h}_L},u_{*}^{\perp }-\sqrt{g{h}_{*}})& if{h}_L>0\end{array} \\ {S}_R=\{\begin{array}{ll}{u}_L^{\perp }+2\sqrt{g{h}_L}& if{h}_R=0\\ \max (u_R^{\perp }+\sqrt{g{h}_R},u_{*}^{\perp }+\sqrt{g{h}_{*}})& if{h}_R>0\end{array} \\ {S}_{*}=\frac{S_L{h}_R(u_R^{\perp }-S_R)-S_R{h}_L(u_L^{\perp }-S_L)}{h_R(u_R^{\perp }-S_R)-h_L(u_L^{\perp }-S_L)} \end{gathered}$$

(13)

where the superscript $\perp$ denotes the quantities of flow variables normal to the considered edge, for example, $u^{\perp }=u{n}_{fx}+v{n}_{fy}$, and $h_{*}$ and $u_{*}^{\perp }$ represent, respectively, the water depth and flow velocity normal to the considered edge in the star region and can be evaluated based on the two-rarefaction assumption [31]:

$$\begin{gathered} h^{*}={\left[\left(u_L^{\perp }-u_R^{\perp }\right)/4+\left(\sqrt{g{h}_L}+\sqrt{g{h}_R}\right)/2\right]}^2/g \\ {u}_{*}^{\perp }=\left(u_L^{\perp }+u_R^{\perp }\right)/2+\left(\sqrt{g{h}_L}-\sqrt{g{h}_R}\right) \end{gathered}$$

(14)

$F_{*L}$ and $F_{*R}$ can be calculated in a simplified version based on the HLL approximated Riemann solver [32]:

$$F_{*L}=\left[\begin{array}{c}{\mathrm{F}}_{*}^{(1)}\\ {\mathrm{F}}_{*}^{(2)}{n}_{fx}-u_L^{\text{P}}{\mathrm{F}}_{*}^{(1)}{n}_{fy}\\ {\mathrm{F}}_{*}^{(2)}{n}_{fy}+u_L^{\text{P}}{\mathrm{F}}_{*}^{(1)}{n}_{fx}\end{array}\right]F_{*R}=\left[\begin{array}{c}{\mathrm{F}}_{*}^{(1)}\\ {\mathrm{F}}_{*}^{(2)}{n}_{fx}-u_R^{\text{P}}{\mathrm{F}}_{*}^{(1)}{n}_{fy}\\ {\mathrm{F}}_{*}^{(2)}{n}_{fy}+u_R^{\text{P}}{\mathrm{F}}_{*}^{(1)}{n}_{fx}\end{array}\right]$$

(15)

where the superscript $\text{P}$ denotes the quantities of flow variables tangential to the considered edge, for example, $u^{\text{P}}=-u{n}_{fy}+v{n}_{fx}$, and ${\mathrm{F}}_{*}^{(1)}$ and ${\mathrm{F}}_{*}^{(2)}$ are two components of the numerical flux vector calculated by the HLL approximated Riemann solver:

$$F_{*}=\frac{S_R{F}_L^{\perp }-S_L{F}_R^{\perp }+S_R{S}_L\Delta U^{\perp }}{S_R-S_L}$$

(16)

where $F_L$, $F_R$ and $\Delta U$ are all column vectors with two components, $F_L^{\perp }=F({\mathbf{U}}_L^{\perp })$, $F_R^{\perp }=F({\mathbf{U}}_R^{\perp })$, and $\Delta U={\mathbf{U}}_R^{\perp }-{\mathbf{U}}_L^{\perp }$, respectively.

3.3. Treatment of Source Term

Simulating shallow water flows over irregular topography demands appropriate treatment of the bed slope source term to achieve well-balanced property. A bed slope source term discretization method developed by Hou [33] is adopted in this work. A brief introduction is given below.

The integral form of the bed slope source term in cell i can be rewritten as:

$${\displaystyle \underset{{\Omega }_i}{\int }{S}_b}d\Omega ={\displaystyle \underset{{\Gamma }_i}{\oint }{F}^{S_b}}(U)nd\Gamma ={\displaystyle \sum _{k=1}^N{F}_{ik}^{S_b}}(U)n_{fk}{l}_{ik}$$

(17)

Integrating Equation (17) with the first-order time marching scheme of SWEs yields:

$${\mathbf{U}}_i^{n+1}={\mathbf{U}}_i^n-\frac{\Delta t^n}{{\mathrm{A}}_{\mathrm{i}}}{\left({\displaystyle \sum _{k=1}^N{F}_{ik}^{S_b}}(U)n_{fk}{l}_{ik}-{\displaystyle \sum _{k=1}^N{F}_{ik}(U)}{l}_{ik}\right)}_n$$

(18)

where the bed slope source term has been transformed into the sum of fluxes across the edges, expressed as:

$$F_{}^{S_b}{n}_f=\left[\begin{array}{c}0\\ -g{n}_{fx}(h_f^L+h_L)(z_{bf}-z_{bL})/2\\ -g{n}_{fy}(h_f^L+h_L)(z_{bf}-z_{bL})/2\end{array}\right]$$

(19)

where $h_L$ and $z_{bL}$ are, respectively, the water depth and bottom level at the centroid of cell i; and $h_f^L$ and $z_{bf}$ represent the reconstructed values of the water depth and the bottom level obtained by Equations (43) and (41), respectively.

With respect to the friction source term, the splitting implicit method proposed in [20] is adopted in this work to properly handle the friction source term:

$${\mathbf{U}}^{n+1}={\mathbf{U}}^n+\Delta t{\overline{S}}_f^n$$

(20)

where the continuity equation is omitted and ${\overline{S}}_f$ is the implicit friction source term expressed by:

$${\overline{\mathrm{S}}}_{fx}^n={\left(\frac{{\mathrm{S}}_{fx}}{1-\Delta t\left(\partial {\mathrm{S}}_{fx}/\partial (hu)\right)}\right)}^n,{\overline{\mathrm{S}}}_{fy}^n={\left(\frac{{\mathrm{S}}_{fy}}{1-\Delta t\left(\partial {\mathrm{S}}_{fy}/\partial \left(hv\right)\right)}\right)}^n$$

(21)

where ${\overline{\mathrm{S}}}_{fx}$ and ${\overline{\mathrm{S}}}_{fy}$ are limited by a criterion:

$${\overline{\mathrm{S}}}_f\{\begin{array}{l}\ge -U/\Delta tifU\ge 0\\ \le -U/\Delta tifU\le 0\end{array}$$

(22)

3.4. Time Integration Scheme

To achieve second-order accuracy in time, the two-stage Runge–Kutta time integration scheme is adopted in this work to update solutions to the new time step. The detailed description of this time integration scheme is given by:

$${\mathbf{U}}_i^{n+1}=\left[\left(U_i^n+U_i^{n^{*}}\right)+K\left(U_i^{n^{*}}\right)\right]/2$$

(23)

where

$${\mathbf{U}}_i^{n^{*}}={\mathbf{U}}_i^n+K\left(U_i^n\right)$$

(24)

where $K\left(U_i^n\right)$ is evaluated by:

$$K\left(U_i^n\right)=\frac{\Delta t^n}{A_i}\left({\displaystyle \underset{{\Omega }_i}{\int }S}d\Omega -{\displaystyle \sum _{k=1}^3{F}_{ik}(U_i^n)}{l}_{ik}\right)$$

(25)

Stability Criteria

The time step used in an explicit time integration scheme must be limited by some stability conditions, otherwise errors could be produced and convergence difficulties could arise. The Courant–Friedrichs–Lewy (CFL) condition is used in this work to estimate the time step:

$$\Delta t=\mathrm{CFL}\cdot \min {(\frac{R}{\sqrt{u^2+v^2}+\sqrt{gh}})}_{i}i=1,2,\dots ,N_t$$

(26)

where CFL is the Courant number, specified in the range of (0, 1] (CFL = 0.5 is used in this work); R denotes the minimum distance from the centroid to the boundary of cell i; u, v and h are the values of flow variables located at the cell centroid; and $N_t$ is the total number of cells.

3.5. Second-Order Spatial Reconstructing and Limiting Method

To achieve high-order accuracy in space, the MUSCL scheme introduced by Van Leer [21] has been a widely used numerical technique for its inherent simplicity and adaptation capacity [22]. Numerical results show that the MUSCL scheme can effectively reduce diffusion in the original Godunov-type scheme [34], leading to improved resolution of water surface where large gradient or discontinuity occurs. Specifically, there are two types of MUSCL schemes, depending on how the values of flow variables are reconstructed: the monoslope and multislope MUSCL schemes. The primary difference between these two schemes is that the former calculates only one vectorial slope per cell, whereas the latter calculates three specific scalar slopes for each edge per cell. It is notable that the monoslope MUSCL scheme involves an additional step for selecting the best limiter [35]. Correspondingly, for the multislope MUSCL scheme, scalar slopes are calculated in the same loop of evaluating numerical fluxes, which apparently leads to a more straightforward and efficient scheme. The multislope MUSCL scheme appears to be more robust when dealing with complex shallow water problems and more accurate than the monoslope one [34]. Moreover, the edge-based limiters chatter far less than the cell-based limiters, and better converging performance can be achieved. A review of these two schemes is presented below. Furthermore, a novel multislope MUSCL scheme is detailed.

3.5.1. A Review of Monoslope MUSCL Scheme

The basic idea of the monoslope MUSCL scheme is briefly reviewed in this subsection by introducing the limited central difference (LCD) scheme [27], which is probably the most widely used MUSCL-type scheme to date. The LCD scheme consists of two stages: the reconstructing stage and the limiting stage. As illustrated in Figure 2, U denotes one of the components of the flow variable vector U and $T_i{M}_{ij}$ denotes the direction vector from $T_i$ to $M_{ij}$. The value of $U_L^M$ can be extrapolated from the centroid $T_i$ with a limited slope:

$$U_L^M=U_i+{\phi }_i\nabla U_i\cdot T_i{M}_{ij}$$

(27)

where $\nabla U_i$ represents the gradient of U in cell i, which can be evaluated by neighbors’ values with the help of the least squares method [36] or Green–Gauss method [37]; ${\phi }_i$ denotes the limiting function:

$${\phi }_i=\underset{k=1,2,3}{\min }{\alpha }^k$$

(28)

where the expression of ${\alpha }^k(k=1,2,3)$ can be found in [27].

3.5.2. Review of Multislope MUSCL Scheme

A recent representative multislope MUSCL scheme proposed by Touze [34] is introduced in the following, which essentially extends the original multislope MUSCL schemes [17,22], and based on that a review of the multislope MUSCL scheme is presented.

A wide stencil is illustrated in Figure 4, which consists of cell i and its neighbors sharing at least a common node; $T_{i1}^{-}$ and $T_{i2}^{-}$ represent the two most upwind neighboring centroids relative to $T_i$ in the direction of $T_i{M}_{ij}$; $T_{i1}^{+}$ and $T_{i2}^{+}$ represent the two most downwind neighboring centroids relative to $T_i$ in the direction of $T_i{M}_{ij}$; and $C_i^{-}$ and $C_i^{+}$ are the intersections of $T_i{M}_{ij}$ with $T_{i1}^{-}{T}_{i2}^{-}$ and $T_{i1}^{+}{T}_{i2}^{+}$, respectively. These two points are exactly the upwind and downwind points relative to $T_i$ in the direction of $T_i{M}_{ij}$. The highlight as well as the challenge of Touze’s scheme is to construct the upwind and backwind points $C_i^{-}$ and $C_i^{+}$, where the values of flow variables are used to calculate the upwind and downwind slopes, $p_{ij}^{-}$ and $p_{ij}^{+}$. These two slopes are then employed as arguments in the limiting function to calculate the directional scalar slope for the considered edge.

To sum up, the principle of the multislope MUSCL scheme is to construct the upwind and downwind slopes. Despite the merits of better robustness, higher accuracy, and better converging performance compared with the monoslope scheme, the drawback of the multislope MUSCL scheme is also obvious: a wider stencil is necessary for evaluating the slopes, which makes the scheme not straightforward because of the procedure for selecting the most upwind and downwind points, so additional computation cost is inevitable. Hence, a more straightforward and efficient numerical scheme is proposed in this work that can efficiently avoid the complicated procedure, while at the same time, the merits of the multislope MUSCL scheme are perfectly inherited.

3.5.3. A Novel Multislope MUSCL Scheme

Based on the basic idea of the multislope MUSCL scheme proposed by Touze, a novel multislope MUSCL scheme equipped with the edge-based Van Albada slope limiter, which takes into account the geometric characteristics of triangular grids, is presented below.

Considering the CCFV scheme adopted in this work, the first thing to calculate is the values of flow variables located at cell nodes, which can be evaluated from neighboring centroids with the method of inverse distance weighting (IDW) [38]:

$$U_{N_k}={\displaystyle {\sum }_{i=1}^{n_c}\left({\omega }_i{U}_{T_i}\right)}/{\displaystyle {\sum }_{i=1}^{n_c}{\omega }_i}$$

(29)

where $U_{N_k}$ denotes the value of the flow variable U located at $N_k$; $n_c$ is the total number of neighboring cells of $N_k$; $U_{T_i}$ denotes the value of the flow variable U located at $T_i$, which is one of the neighboring centroids of $N_k$; and ${\omega }_i$ denotes the weighting coefficient related to the distance between $N_k$ and $T_i$:

$${\omega }_i=1/\left|T_i{N}_k\right|$$

(30)

Two different cases should be taken into consideration, as shown in Figure 5. The criterion is straightforward and independent with the shapes of triangular grids:

$$\{\begin{array}{l}\mathrm{Case}\left(\mathrm{a}\right)if{M}_{ij}{T}_j\times M_{ij}{T}_i\le 0\\ \mathrm{Case}\left(\mathrm{b}\right)if{M}_{ij}{T}_j\times M_{ij}{T}_i>0\end{array}$$

(31)

where the coordinates of the midpoint $M_{ij}$ can be obtained by those of $N_2$ and $N_3$. In case (a), the intersection $C_i^{+}$ of $M_{ij}{T}_i$ and cell j locates at the edge ${N^{\prime }}_2{N}_3$, whereas in case (b), the intersection $C_i^{+}$ locates at the edge $N_2{N^{\prime }}_2$. For the sake of clarity, the following steps are only for case (a).

Similar to Touze’s scheme, cell j can be seen as the most downwind cell relative to cell i, and the intersection $C_i^{+}$ can be regarded as the downwind point in the direction of $T_i{M}_{ij}$. The value of $U_{ij}^{+}$ located at $C_i^{+}$ is evaluated by:

$$U_{ij}^{+}={\nu }_{i1}{U}_{{N^{\prime }}_2}+{\nu }_{i2}{U}_{N_3}$$

(32)

where ${\nu }_{i1}$ and ${\nu }_{i2}$ denote the weighting coefficients with respect to the neighboring nodes of $C_i^{+}$ , that is, ${N^{\prime }}_2$ and $N_3$ in case (a), which can be evaluated by:

$${\nu }_{i1}=\left|C_i^{+}{N}_3\right|/\left|{N^{\prime }}_2{N}_3\right|, {\nu }_{i2}=\left|C_i^{+}{N^{\prime }}_2\right|/\left|{N^{\prime }}_2{N}_3\right|$$

(33)

As shown in Figure 5, the line connecting $T_i$ and $M_{ij}$ also passes through the node $N_1$, which is selected as the upwind point in the direction of $T_i{M}_{ij}$. Considering that the value located at $N_1$ is known, the upwind and downwind slopes for the considered edge are almost ready and given by:

$$p_{ij}^{-}=\left(U_{T_i}-U_{N_1}\right)/\left|T_i{N}_1\right|, p_{ij}^{+}=\left(U_{ij}^{+}-U_{T_i}\right)/\left|T_i{C}_i^{+}\right|$$

(34)

where $U_{N_1}$ denotes the value of U located at $N_1$.

The modified Van Albada limiting function ${\phi }^{VA}$ (Equation (35)) [17,26,39] is employed here to turn the upwind and downwind slopes into a limited slope, which is then used to calculate the left and right Riemann states:

$${\phi }^{VA}\left(a,b\right)=\{\begin{array}{ll}\zeta a+(1-\zeta )b& ifab>0\\ 0& ifab\le 0\end{array}$$

(35)

where $\zeta$ ranges between 0 and 1 and is defined by Equation (36), and $a$ and $b$ denote, respectively, the upwind and downwind slopes; and

$$\zeta (a,b)=\left(b^2+e\right)/\left(a^2+b^2+2e\right)$$

(36)

where e is a constant used to prevent division by zero and $e={10}^{-6}$ is adopted in this work.

The flow states on the considered edge, for example, the left flow state of the flow variable U located at $M_{ij}$, can be calculated by:

$$U_L^M=U_{T_i}+{\phi }^{VA}(p_{ij}^{-},p_{ij}^{+})\left|T_i{M}_{ij}\right|$$

(37)

Similarly, the flow states on the other edges can be obtained by repeating the above steps.

In order to preserve the C-property, only $\eta$, $h$, $hu$ and $hv$ are reconstructed [14], which gives rise to the left and right Riemann states, i.e., ${\eta }_L$, $h_L$, ${(hu)}_L$, ${(hv)}_L$, ${\eta }_R$, $h_R$, ${(hu)}_R$ and ${(hv)}_R$. The bottom level and velocity components in x and y directions are subsequently obtained by:

$$z_b=\eta -h,u=(hu)/h,v=(hv)/h$$

(38)

A criterion proposed by Hou [33] is adopted here to identify problematic cells next to the wet-dry interface:

$$h_L\le \min (\left|z_{bL}-z_{bi}\right|,0.25h_i)\text{ or }{h}_i\le {\epsilon }_{wd}$$

(39)

where $h_L$ and $z_{bL}$ denote, respectively, the water depth and bottom level to the left of the considered edge; $h_i$ and $z_{bi}$ denote, respectively, the water depth and bottom level at the centroid of cell i; and ${\epsilon }_{wd}$ is the threshold of water depth used to distinguish wet and dry cells, and ${\epsilon }_{wd}={10}^{-6}\mathrm{m}$ [40,41] is adopted in this work.

3.6. Treatment of the Wet-Dry Interface

The second-order MUSCL scheme degenerates to first-order accordingly once the criterion in (Equation (39)) is met, that is, the left and right Riemann states are respectively constant with cell-averaged values. By doing so, abnormal reconstructed values are effectively avoided. Furthermore, numerical schemes integrated with an appropriate Riemann solver can properly handle the transience of flow state from wet to dry zones, as long as the topography is horizontal or slowly varying. However, realistic topography is quite inhomogeneous: adverse or abruptly steep slopes are common in cells next to the wet-dry interface. Water depth gets smaller and smaller as flow advances, leading to extremely high velocity, which is unacceptable for simulation. Hence, particular numerical techniques must be used in the representation of wet-dry interface. The classic hydrostatic reconstruction method [14] and the local bed modification technique [18,42] are combined in this work to deal with the wetting and drying problem, stated as follows.

A single value of the bottom level at the common interface is defined by:

$$z_b=\max (z_{bL},z_{bR})$$

(40)

Considering those sensitive cells over adverse slopes, a bed modification strategy is carried out to hold the stopping flow condition:

$$z_b\leftarrow z_b-\max (0,z_b-{\eta }_L)$$

(41)

The water surface elevation on the right side of the common interface is modified accordingly:

$${\eta }_R\leftarrow {\eta }_R+(z_b-z_{bR})$$

(42)

Then the left- and right-side water depths are reevaluated as:

$$h_L={\eta }_L-z_b,h_R={\eta }_R-z_b$$

(43)

After these steps, ${(\Delta z_b)}_{LR}=-{(\Delta h)}_{LR}=0$ is obtained at the wet-dry interface, which is in agreement with the equilibrium condition proposed by Brufau [42]. The water surface elevation over adverse slopes is equal to η_L, so spurious flow motion is effectively prevented and the well-balanced property is perfectly preserved.

3.7. Boundary Conditions

Boundary conditions provide information about how the water along boundaries responds to external disturbances. In order to achieve an optimal simulation for shallow water flows, an adequate treatment of boundary conditions should be considered. Basically, there are two approaches to implementing boundary conditions: using ghost cells and directly calculating fluxes across boundaries. In this work, the latter is adopted by solving Riemann problems along boundaries. Boundaries need to be divided into two distinct sections: the liquid boundary and the solid boundary. The former can be further subdivided into four different cases.

3.7.1. Liquid Boundary

There are four cases to take into account for the liquid boundary, depending on the flow velocity normal to boundaries $u^{\perp }$ and the local Froude number $Fr$, expressed as $Fr=\left|u^{\perp }\right|/\sqrt{gh}$ [43]:

Supercritical inflow ($Fr>1$):

$$u_R^{\perp }=u_c^{\perp }, u_R^{\text{P}}=0, h_R=h_c$$

(44)

Supercritical outflow ($Fr>1$):

$$u_R^{\perp }=u_L^{\perp }, u_R^{\text{P}}=u_L^{\text{P}}, h_R=h_L$$

(45)

Subcritical inflow ($0<Fr\le 1$):

$$u_R^{\perp }=u_c^{\perp }, u_R^{\text{P}}=0, h_R={\left[\sqrt{h_L}-(u_R^{\perp }-u_L^{\perp })/2\sqrt{g}\right]}^2$$

(46)

Subcritical outflow ($0<Fr\le 1$):

$$u_R^{\perp }=u_L^{\perp }+2\sqrt{g}(\sqrt{h_L}-\sqrt{h_R}), u_R^{\text{P}}=u_L^{\text{P}}, h_R=h_c$$

(47)

In order to calculate the numerical fluxes through boundaries, a reverse transforming operation is applied to velocities as soon as $u_R^{\perp }$ and $u_R^{\text{P}}$ have been obtained:

$$u^R=u_R^{\perp }{n}_{fx}-u_R^{\text{P}}{n}_{fy}, v^R=u_R^{\text{P}}{n}_{fx}+u_R^{\perp }{n}_{fy}$$

(48)

where $u_R^{\text{P}}$ is equal to $u_L^{\text{P}}$ along boundaries.

Once the flow states to the right of boundaries are obtained, a local Riemann problem can be set up and then settled by the HLLC Riemann solver, where the left flow states of the Riemann problem can be evaluated by averaging those at the two adjacent nodes located on the boundary edges, since they have been obtained from neighboring centroids with the IDW method.

3.7.2. Solid Boundary

Normal velocity is ignored in this situation, and the value of water depth is taken from the internal computational domain, that is, $h_R=h_L$. The normal fluxes at the midpoint of the considered edge along solid boundaries can be evaluated by a simplified form of Equation (11):

$$F_{ik}^{solid}=\left[\begin{array}{c}0\\ (g{h}_R^2/2)n_{fx}\\ (g{h}_R^2/2)n_{fy}\end{array}\right]$$

(49)

4. Verification and Application

4.1. Stationary Flow with Wet-Dry Interface

A well-balanced scheme should correctly compute the solutions to a steady state and maintain the stationary state throughout the simulation for a quiescent shallow water problem. A simple test is employed here, as also used in [44,45], to verify the well-balanced property of the new scheme. A hump is defined in an 8000 m long square basin with topography given by:

$$z_b\left(x,y\right)=\max \left(0,2000-0.00032\left({\left(x-4000\right)}^2+{\left(y-4000\right)}^2\right)\right)$$

(50)

where the solid boundary condition is imposed on vertical walls. The initial water surface is given by 1000 m, so the hump is partly submerged. The computational domain is discretized into 15,360 Delaunay triangles, and the simulation is run until t = 1500 s.

A variable $R^n(U)$ is defined to quantify the global relative error of solutions following Zhou [45]:

$$R^n(U)=\sqrt{{\displaystyle \sum _{N_w}\left(A_i\cdot {\left(\frac{U_i^n-U_i^{n-1}}{U_i^n}\right)}^2\right)}/{\displaystyle \sum _{N_w}{A}_i}}$$

(51)

where ${\mathbf{U}}_i^n$ and ${\mathbf{U}}_i^{n-1}$ denote, respectively, the vectors of flow variables at time level n and n − 1 for cell i.

Figure 6 gives the numerical results of water surface and unit width discharges in x and y directions. The stationary state is visible at the end of the simulation, since the numerical water surface elevation almost coincides with the analytical result. With regard to the numerical unit width discharge, a little wriggle can be observed in Figure 6b in terms of $q_x$. As a matter of fact, the global relative error is quite small, as illustrated in Figure 7, where both $R(q_x)$ and $R(q_y)$ hold rather small values (lower than ${10}^{-30}$) throughout the simulation; such a result is acceptable and the well-balanced property is well preserved.

4.2. Potential Flow over Uneven Bottom

As proposed by Ricchiuto [46], the steady flow with divergence-free velocities over uneven bottom is simulated here to verify the capacity of the new scheme to preserve monotonicity and converge to steady-state solutions in the presence of irregular topography. This test case was performed in a smooth square domain, $[-1\mathrm{m},1\mathrm{m}]\times [-1\mathrm{m},1\mathrm{m}]$, where the velocity field is obtained by the harmonic function:

$$\psi \left(x,y\right)=xy$$

(52)

The velocities in $x$ and $y$ directions are expressed as the partial derivatives of $\psi \left(x,y\right)$:

$$u\left(x,y\right)=\frac{\partial \psi }{\partial y}=x,v\left(x,y\right)=-\frac{\partial \psi }{\partial x}=-y$$

(53)

The analytical water depth and bottom level of the test case are calculated by:

$$h\left(x,y\right)=1.5+\psi ,z_b\left(x,y\right)=\left(30-\left(x^2+y^2\right)/2\right)/g-\psi -1.5$$

(54)

Considering the Froude number inside the computational domain is always less than 1, the north and south edges are used as subcritical inlet boundaries and the east and west edges are used as subcritical outlet boundaries during the simulation. The computational domain is discretized into 6341 Delaunay triangles, as shown in Figure 8. The simulation is run with the initial conditions calculated by Equations (53) and (54) until t = 50 s.

Figure 9a,b illustrate, respectively, the numerical results for water depth and flow velocities computed at the end of the simulation by the new scheme proposed in this work. It can be observed that the new scheme is able to converge to the steady-state results in terms of water depth and flow velocities on Delaunay triangles at the end of the simulation. The contours of the water depth are symmetric about two diagonal lines of the computation domain, which agrees well with the analytical results, since the analytical water depth (Equation (54)) matches a hyperbolic-paraboloid distribution. In addition, the velocity field of this test is also well predicted at the end of the simulation, because the distribution of flow directions is very regular, as shown in Figure 9b, and obviously there are no unphysical high velocities inside the computational domain. The numerical results confirm that monotonicity in the limitation procedure is guaranteed, and the new scheme can achieve convergent results in the case of irregular topography and nonzero velocities.

4.3. Water Sloshing in a Parabolic Basin

The test proposed by Thacker [1] on water oscillation in a smooth 2D parabolic basin, as illustrated in Figure 10, has been a classical benchmark for evaluating the accuracy of and capacity to deal with the wetting and drying problem of numerical schemes [17,47,48] and is employed in this subsection. There is a parallel test, developed in [49], that takes the effect of friction source term into consideration and is widely used by many other researchers [33,50]. The bottom level of the parabolic basin is defined by:

$$z_b(x,y)=\frac{h_0}{a^2}(x^2+y^2)\left(-\mathrm{10,000}\mathrm{m}\le x,y\le \mathrm{10,000}\mathrm{m}\right)$$

(55)

where $h_0$ denotes the water depth at the center of the basin and $a$ denotes the distance from the center to the shoreline with zero elevation; the analytical solutions of velocities and water depth are:

$$\begin{gathered} u=-\sigma \omega \sin (\omega t) \\ v=-\sigma \omega \cos (\omega t) \\ h\left(x,y,t\right)=\max \left(0,\frac{\sigma h_0}{a^2}\left(2x\cos \left(\omega t\right)+2y\sin \left(\omega t\right)-\sigma \right)+h_0-z_b\left(x,y\right)\right) \end{gathered}$$

(56)

where σ is a constant denoting the amplitude of the water motion and ω denotes the frequency of the water motion and is calculated by $\sqrt{2g{h}_0}/a$; values of the parameters used in this test are $a=8025.5\mathrm{m}$, $h_0=10\mathrm{m}$, $\sigma =802.55\mathrm{m}$, $\omega =1.7453\times {10}^{-3}{\mathrm{s}}^{-1}$, and $T=3600\mathrm{s}$.

The computational domain is discretized into four scaled meshes with sizes of ∆L = 55.18 m, 108.49 m, 138.58 m and 205.02 m, where ∆L is the averaged grid size, calculated by ${\displaystyle \sum _{i=1}^n\sqrt{A_i}}/n$. For each scaled mesh, the number of Delaunay triangles is 128,418, 33,038, 20,516 and 9304, respectively. The analytical solutions at t = 0 calculated by Equation (56) are used as the initial conditions. In addition, the solid boundary condition is employed in this test. This simulation is carried out for three periods, i.e., t = 10,800 s.

The comparison between the analytical and numerical results of water surface along the x-direction symmetric axis at t = T, T + T/3, T + T/2 and 2T + T/6 for the smallest scaled mesh is depicted in Figure 11. Compared with the analytical solutions, the numerical water surface and shoreline show good agreement with the exact results, meaning the new scheme can accurately capture the moving boundaries and handle the wetting and drying problem well. Moreover, two components of flow velocity at the point (−5000, −100) versus time are also depicted in Figure 12, demonstrating the difference between the numerical and analytical flow velocities, which is almost negligible.

Figure 13 shows the relative error of the total volume of water during the simulation over time. The relative error ranges from $-5.42\times {10}^{-17}$ to $-5.43\times {10}^{-19}$, indicating a rather small deviation from the initial total water volume, and the new scheme proposed in this work ensures mass conservation. The errors in $L_1$ and $L_2$ norms, reflecting the deviation between numerical and analytical solutions, are widely used to investigate the convergence performance of numerical schemes and are also used in this work. Water depth and two components of unit width discharge at t = 3T versus the grid size ∆L in logarithm scale are plotted in Figure 14. It is obvious that, for all flow variables, the convergence rate in errors in both $L_1$ and $L_2$ norms is greater than 1.6. The accuracy is acceptable, considering the second-order reconstruction procedure has been switched to first-order in some problematic cells near the wet-dry interface.

4.4. Steady Flow over Frictional Uneven Bottom

A two-dimensional steady-state flow over uneven bottom considering friction effects, originally proposed by Murillo [51] and also adopted in [7,10,52], is employed here to test the new scheme when perturbations by friction source term are involved. The analytical water depth and bottom level satisfy Equations (57) and (58):

$$h(x,y)=a+q_{x}x+q_{y}y$$

(57)

where a is constant and equal to 0.5, and

$$z_b(x,y)=-\frac{1}{2g}\frac{\left(q_x^2+q_y^2\right)+2g{h}^3}{h^2}+\frac{3}{7}\frac{\left|q_x\right|n^2\sqrt{2}}{h^{7/3}}$$

(58)

The computational domain is a square of 10 m. The unit width discharges in Equation (57) are assumed to be constant over the domain and $q_x=q_y=0.1{\mathrm{m}}^2/\mathrm{s}$. The Manning coefficient n throughout the domain is also a constant value taken as 0.03 $\mathrm{s}/{\mathrm{m}}^{1/3}$. As shown in Figure 15, in order to study the adaptability and robustness of the new scheme on various irregular triangular grids, four types of triangular meshes are adopted in this simulation under the same initial and boundary conditions: the Delaunay mesh, the orthogonal-I mesh, the orthogonal-II mesh, and the distorted mesh. For each type of mesh, the number of discrete triangles is 1020, 800, 800 and 1111, respectively, and the number of nodes is 551, 441, 441 and 602, respectively. Once again, the global relative error $R^n(U)$ defined in Equation (51) is employed in this test to quantify the convergent performance of the new scheme on different types of meshes, and $R(h)<{10}^{-8}$ is used as the criterion of convergent results.

Before computation, the initial water surface is set to zero and the still water condition $u=v=0$ is assumed over the domain. The constant unit width discharges $q_x=q_y=0.1{\mathrm{m}}^2/\mathrm{s}$ are imposed on the west and south boundaries. The east and north boundary conditions are given by the water depth calculated by Equation (57). The simulation runs on four types of meshes for 25 s with a time step of 0.005 s.

Figure 16 illustrates the global relative errors of the new scheme in terms of water depth for different types of meshes. In these four series of runs, the new scheme is able to converge to the steady state on all the mesh types at the end of the simulation ($R(h)<{10}^{-8}$), even in the case of distorted mesh. Figure 17 displays the exact and numerical water surface at the end of the simulation along the diagonal section (Figure 17a) and horizontal-symmetry section (Figure 17b) of the computational domain on the distorted mesh. The comparisons show good agreement between the numerical and exact water surface even on the distorted grid (the difference is rather small, on the order of ${10}^{-5}\mathrm{m}$, for both recorded sections). In summary, this test demonstrates not only the validity of the friction source term treatment introduced in Section 3.3, but also the flexibility of the new multislope MUSCL scheme for irregularly shaped grids even in the presence of obtuse triangles.

4.5. Flow in a Compound Channel

This test considers a commonly seen scenario in natural rivers: an asymmetric compound channel, as sketched in Figure 18. The physical experiment was conducted by Rajaratnam [53]. The channel was 18.3 m long and 1.22 m wide. The most significant characteristic of this channel was the distinctly different velocity distributions in the main channel and the floodplain, since a high step separates the main channel and the floodplain. Hence, the test simulated here is used to verify the capability of the new scheme to give an accurate prediction of different flow regimes caused by rapidly varying topography. For further comparison, three kinds of MUSCL-type schemes are used: the LCD scheme, Touze’s scheme, and the new scheme proposed in this work. Detailed parameters are listed in

Table 1. Data of the physical experiment.

D (cm)	d (cm)	h (cm)	B (cm)	b (cm)	Q (m3/s)	Sw × 103	Fr
11.28	1.52	9.75	71.1	50.8	0.027	0.45	0.37

.

The computational domain is discretized into 17,164 Delaunay triangles. In order to sharpen the resolution of the topography equipped with an abrupt step, the computation domain is discretized using a longitudinal break line that exactly coincides with the line dividing the main channel and the floodplain in the horizontal plane. Triangles are distributed on both sides of the break line, resulting in the break line being discretized with triangle edges. The inlet boundary condition is given by a uniform velocity and the outlet boundary condition is given by a constant water surface elevation. In addition, the channel is supposed to be smooth, so the friction source term is negligible.

Figure 19 plots a comparison of the measured and simulated longitudinal velocities along the cross-section located 9.15 m from the inlet boundary. Three schemes can reproduce the dramatic decrease in streamwise velocity from the main channel to the floodplain. Correspondingly, a large gradient of streamwise velocity occurred at the interface of the main channel and the floodplain (the dotted line in Figure 19) due to an abrupt change of flow depth, consistent with observations in the published literature [54,55]. However, all the schemes exhibit some sort of numerical diffusion near the step, which leads the numerical results to deviate from the measured data. It means the two-dimensional shallow water model has an inherent limitation in dealing with this kind of case with sharp changing topography. It is also obvious that the numerical diffusions drag the flow momentum in the main channel to the floodplain and narrow the velocity gap between these two regions.

It is worth noting that compared with the LCD scheme, both multislope schemes perform better in capturing the large gradient of velocity, which is also stated by Delis [17] and Hubbard [27]. Moreover, slight differences can be found between the numerical results of Touze’s scheme and the new scheme proposed in this work. Similar accuracies are achieved by these two schemes and parallel large-slope capturing capabilities can be seen in Figure 19. However, the new scheme is more efficient than Touze’s scheme and the LCD scheme (4.353% and 4.976% runtime, respectively, saved in this test). This is because the new scheme avoids the addition loop for selecting an appropriate limiter in the LCD scheme and makes full use of the geometric characteristic of triangular grids, so the upwind point used for calculating the upwind slope is already available, which yields a more straightforward and efficient algorithm.

4.6. Flow in a Contracting-Expanding Channel over Uneven Bottom

Parshall flume, a widely used device for field measurement in open channels, as shown in Figure 20, is employed to verify the applicability of the new scheme for varying flow regimes.

The corresponding experiment was initially conducted by Ye [56]. The flow discharge is 0.0145 m³/s at the inlet boundary, and the outlet boundary conditions do not influence the flow regime upstream. The water surface along the path is presented in Figure 21a. Zone 1 denotes the converging section where the cross-section contracts over a flat bottom to create a critical depth; zone 2 denotes the throat section over a sloping bottom where supercritical flow occurs; and zone 3 denotes the diverging section over an adverse slope bottom where the flow is supercritical. The flow regime along the flume varies from subcritical to critical and eventually reaches supercritical conditions downstream of the throat section. Complexly varying flow regimes present a great challenge to the accuracy and stability of the simulation. A robust numerical scheme should provide an accurate prediction for different flow regimes. The parameters of this experiment are listed in

Table 2. Values of parameters of the experiment.

A (m)	B (m)	C (m)	D (m)	E (m)	F (m)	G (m)	M (m)	N (m)	T (m)	W (m)
0.389	0.381	0.305	0.305	0.267	0.152	0.483	2.54	0.086	0.889	0.152

.

For numerical computation, the domain is discretized into 10,567 Delaunay triangles. The inlet boundary condition is given by u = 0.4 m/s and h = 0.12 m. The outlet boundary condition is not imposed since the supercritical flow regime. The initial condition is given by h = 0.05 m. In addition, the channel is supposed to be smooth.

Figure 21a depicts the distribution of water surface along the Parshall flume, including the numerical result computed by the scheme proposed in this work and the one provided in [56]. It can be seen that there is excellent agreement of the water surface profile between the numerical and measured results along the flume. Compared with the numerical results in [56], the numerical result in this work is obviously closer to the experimental data regardless of flow regime. This phenomenon is more pronounced for Ye’s numerical results in the supercritical flow regime downstream, where a large deviation from the experimental data occurs. Figure 21b illustrates the flow field of the numerical x-directional velocity calculated by the new scheme, which is shown to agree closely with the experimental data, although the velocity downstream is slightly greater than that of the experimental data. Since the frictional resistance will more or less affect the measured results, the difference of velocity downstream in the flume between the numerical results and measured data is acceptable. It is concluded that this new scheme can provide satisfactory solutions for varying flow regimes.

4.7. A Symmetric 2D Riemann Problem

The last test was proposed by Brio [57] and was widely adopted in other studies [17,58] to verify the shock-capturing capability of numerical schemes. This test aims to simulate the flow pattern arising from two streams intersecting at a straight angle over a frictionless square domain. As presented in Figure 22, the domain is divided into four subregions of equal size, where the bottom levels are all zero. The initial conditions are discontinuous and listed as follows:

The interactions between different regions result in a complex two-dimensional wave pattern, which provides an excellent basis for verification of numerical schemes in terms of shock-capturing ability. The computational domain is discretized into 15,360 Delaunay triangles, and the new multislope MUSCL scheme is used to predict wave propagation. The initial flow conditions are given by the data in

Table 3. Initial conditions for the 2D Riemann problem.

Region (#)	Coordinates (m)	Water Depth (m)	u (m/s)	v (m/s)
1	x ≤ 100, y ≤ 100	1	10	10
2	x > 100, y ≤ 100	1	0	10
3	x ≤ 100, y > 100	1	10	0
4	x > 100, y > 100	10	0	0

and the boundary conditions are not prescribed, to avoid water reflection.

Figure 23 demonstrates the 3D view of the simulated water surface at t = 4 s in a square subregion $[50\mathrm{m},150\mathrm{m}]\times [50\mathrm{m},150\mathrm{m}]$. It is clear that the complex wave propagation pattern is well reproduced. The interactions between regions 1 and 2 yield two shock waves propagating along the x axis in opposite directions. A symmetric wave propagation pattern occurs as a result of the interactions between regions 1 and 3. The interactions between regions 4 and 2 yield a shock wave propagating along the negative y axis and a refraction wave along the positive y axis. A symmetric wave propagation pattern occurs as a result of the interactions between regions 4 and 3, where both a shock wave and a refraction wave arise. The interactions between regions 4 and 1 together with the interactions of the shock waves as mentioned above result in a complex wave propagation pattern.

Figure 24 illustrates the contours of water surface elevation (equidistance = 0.1 m) and the vector field of flow velocities calculated by the new scheme. The contours of the simulated water surface (Figure 24a) and the vector field of flow velocities (Figure 24b) are symmetric about the bottom-left to top-right diagonal line of the domain. All the shock waves are perfectly captured, caused by the interactions between neighboring regions, indicating the excellent shock-capturing capability of the new scheme.

5. Conclusions

A new multislope MUSCL scheme on triangular grids is proposed in this paper within the framework of the Godunov-type CCFV method for simulating 2D shallow water flows over uneven topography. The values of flow variables on both sides of the considered edge are calculated by an edge-based reconstructing and limiting procedure where the upwind slope is easily obtained, since the upwind point is fixed at the node corresponding to the considered edge. Moreover, the downwind slope is also easy to evaluate by means of a straightforward criterion. The algorithm is simplified in the reconstructing and limiting procedures compared with the conventional LCD scheme [27] and other existing multislope MUSCL schemes [17,34], so computation cost is obviously reduced without sacrificing accuracy and robustness.

Seven tests were carried out to verify this scheme. The results show that lower time cost for reaching convergent results is achieved compared with the classical monoslope and multislope MUSCL schemes. In addition, it is also shown that this new scheme is more flexible on triangular grids and more straightforward to implement for taking account of the geometrical characteristics of triangular grids without sacrificing accuracy, robustness, and resolution. Furthermore, this new multislope reconstruction method is independent of HLLC approximated Riemann solvers and SWEs, so it can be integrated with numerical schemes based on other hyperbolic-type equations and Riemann solvers. To sum up, the new scheme provides an ideal alternative to the existing monoslope and multislope MUSCL schemes for numerically solving SWEs and other hyperbolic-type equations.