During flash floods, the water flow movement and the load transported by the flow generally behave like 3D flow, and the mathematical model should simulate the hydrodynamic process in 3D form in order to fully reflect its features. However, due to the complexity of 3D water flow, many research studies have often simplified the flow of shallow water by averaging the moving elements in the direction of water depth, thus transforming the 3D problem into a 2D problem.

Given the assumption that vertical acceleration can be ignored and the water pressure is distributed as hydrostatic pressure along the direction of water depth, an average value was introduced in the direction of water depth in the Cartesian coordinate system, and viscous stress was ignored. Based on NS equation, the 2D shallow water equation can be derived. ∂ U ∂ t + ∇ ⋅ F = S (1)

In the equation, U is the vector of conserved variables; F is the vector flux; and S is the source term vector. U and F (U) are calculated as follows: U = ( ϕ ϕ u ϕ v ) (2) F = ( ϕ V ϕ u V + 1 2 ϕ 2 i ϕ v V + 1 2 ϕ 2 j ) (3)

In the equations, ϕ = g h is the gravitational potential; u is the flow rate in the x direction; v is the flow rate in the y direction; V = u i + v j represents the velocity vector; g = 9.81 m 2 / s is the acceleration of gravity; and h is the water depth.

In many engineering applications, the source term S includes components such as the geostrophic deflecting force, wind force, and bottom friction. For the simulations performed in this research, the source term S is composed of the bottom slope source term S b and the bottom friction source term S f , namely S = S b + S f , where S b and S f can be expressed as follows: S b = ( 0 g ϕ ∂ H ∂ x g ϕ ∂ H ∂ y ) (4) S f = ( 0 − g ρ τ fx − g ρ τ fy ) (5) τ fx = ρg n 2 u u 2 + v 2 1 R 1 3 (6) τ fy = ρg n 2 v u 2 + v 2 1 R 1 3 (7)

In the equations, H is the difference between the fixed horizontal base level and the bottom elevation; ρ is density; τ f x and τ f y are the friction stresses of the river bed in the x and y directions; n is the roughness coefficient; and R is the hydraulic radius, which is used to approximate the water depth, h, in the open channel.

The shallow water equation is a nonlinear hyperbolic partial differential equation and calculating the grid interface flux is an important step in solving the hyperbolic equation. The first-order algorithm has too much dissipation when it uses the unit average to calculate the flux of the interface, smoothing out what should be a steep shock wave (Roe, 1997). Therefore, higher-order formats, such as the second-order TVD method, are favored by researchers around the world. For higher-order Godunov problems, the second-order format is usually obtained using piecewise linear reconstruction of the local variables in the center of the grid.

Due to the existence of non-zero source terms (bottom slope source term and shear stress source term), the non-homogeneous partial differential Eq. 1 can be separated and decomposed into the homogeneous partial differential, Eq. 8, and the ordinary differential Eq. 9, using only the bottom slope source term and the shear stress source term. The equations were solved in two steps: estimation and correction. Namely, the estimation step with a step size of (∆t)/2 was first used to solve the source term U; then, the solved U was substituted into the correction step with step size of ∆t. ∂ U ∂ t + ∇ ⋅ F = 0 (8) dU dt = S b + S f (9)

When the source term of the bottom slope is 0 and the viscosity of the water is ignored; the water depth h can be used to accurately calculate the gravitational potential gradient δ ϕ i . However, the impact from changes in the river bed surface should be considered during an actual flash flood. Hence, the water depth at the unit interface cannot be accurately determined if the gravitational potential is taken as the conserved term. This is because the depth gradient cannot reproduce true water depth variations. So various numerical methods have been developed to improve the treatment of source terms in the shallow water equations (Pu et al., 2012). The surface gradient upwind method (SGUM) proposed by Hui Pu et al. (2013) and surface gradient method (SGM) proposed by Zhou et al. (2001) proved that the hyperbolic equation is conserved when using surface elevation instead of water depth for gradient calculations. Therefore, the algorithm subtracts the river bottom elevation from the surface elevation η i of each point, instead of using the water depth, h i , to calculate the conservation term, ϕ , thereby ensuring the conservation of the bottom slope source term.

In this paper, to solve the Eq. 8, the numerical solution consists of two steps: prediction and correction, and they are shown, respectively, in Eqs 10, 11. ( AU ) i n + 1 2 = ( AU ) i n − Δ t 2 ( ∑ m = 1 M F ( U m ) n ⋅ L m ) (10) ( AU ) i n + 1 = ( AU ) i n − Δ ⋅ t ( ∑ m = 1 M F ( U m L , U m R ) n + 1 2 ⋅ L m ) (11) where, A is the area of the grid; L m is the side vector of the grid, which is defined as the cell edge length multiplied by the outward-pointing unit normal vector; M is the number of cell edges; F ( U m ) n is the unit interface flux of the m-th side of each grid at the n-th time; F ( U m L , U m R ) n + 1 2 is the m-th edge unit grid interface flux of the correction step after the prediction step.

The computational workload required to solve the equations is relatively large; hence, approximate Riemann solutions in formats such as Roe, Osher, and HLL are often used to reduce the processing time. In this paper, the WAF method and HLL format were used to calculate the interface flux. The WAF interface flux (Cao et al., 2004) can be written as follows: F i + 1 2 = 1 2 ( F i + F i + 1 − 1 2 ∑ k = 1 N c k Δ F i + 1 2 ( k ) ) (12) where, c k is the Courant number corresponding to the velocity S k of the k-th train wave. c k = Δ t S k Δ x (13)

The HLL scheme (Hubbard and Dodd, 2002) is a two-wave train model. The Riemann problem is simplified into three conserved regions separated by two waves, S_L and S_R. The interfaces of the three regions are as shown in Figure 2.

In Figure 2, U L and U R represent the conservation vectors on the left and right sides of the interface; U * represents the intra-interface vector; while S L and S R represent the velocities of waves on the left and right sides of the interface. The HLL scheme can be calculated using the following formula: F i + 1 2 hll = { F L 0 ≤ S L S R F L − S L F R + S L S R ( U R − U L ) S R − S L S L ≤ 0 ≤ S R F R 0 ≥ S R (14) where, S L and S R (Toro, 1992) can be estimated as: S L = min   ( u i − ϕ i , u s − ϕ s ) , S R = max ( u i + 1 + ϕ i + 1 , u s + ϕ s ) (15) u s = 1 2 ( u i + 1 + u i ) + ϕ i − ϕ i + 1 (16) ϕ s = ϕ i + ϕ i + 1 2 + u i − u i + 1 4 (17)

The wave velocity algorithm described above is applicable for wet boundaries. For dry and wet mesh boundaries, Da Silva (2017) proposed a corrected model. When the left river bed is dry, the following calculation can be used: S L = u i − ϕ i ,           S R = u i + 2 ϕ i (18)

When the right river bed is dry: S L = u i + 1 − 2 ϕ i + 1 ,         S R = u i + 1 + ϕ i + 1 (19)

After the homogeneous equations were corrected, the total differential Eq. 9 was decomposed, and the estimation step discrete equation of S b is determined as follows: S b = g ϕ ∂ H ∂ x = ∂ ∂ x ( 1 2 ϕ 2 ) (20) ( AU ) i n + 1 2 = ( AU ) i n + Δ t 2 ∑ m = 1 M ( 1 2 ϕ 2 ) n ⋅ L m (21)

For S f , the algorithm adopted a forward difference scheme to solve the final flux ( A U ) i n + 1 2 of the estimation step; then the final flux was substituted into Eq. 11 before iteration of the correction step. Finally, Δ t was used as the step size to calculate ( A U ) i n + 1 in the total differential equation.

In order to ensure the stability of the numerical scheme, the display algorithm needs to adjust the time step ∆t at each calculation step to reduce the deviation in the propagation of the source term along the characteristic line, thereby ensuring the stability of the algorithm and the accuracy of the numerical solution. The Courant number CFL should satisfy the following condition (LeVeque, 2002): 0 ≤ CFL ≤ 1 (22)

In this paper, the value of CFL was assumed to be 0.9 to ensure the stability of the algorithm. Step ∆t was calculated as per the following formula: Δ t = CFL Δ x max 1 ≤ i ≤ M ( | u i − ϕ i | , | u i + 1 − ϕ i + 1 | ) (23)

The key to calculating the suspended load is to determine the suspended matter content, S, and sediment-carrying capacity, S*, of the river. S was calculated using the finite volume method with a first-order upwind finite volume method. The regions upstream and downstream of the adjacent grids were judged according to the flux in the normal direction of the common grid interface, and the sediment flux at each interface of the grid was calculated. Furthermore, the sediment flux at each interface was integrated at the specified time step to obtain the suspended load content of the grid. This algorithm can unconditionally guarantee the stability of the suspended load content calculations.

There has not been much research on the sediment-carrying capacity of 2D water flow. Therefore, the sediment-carrying capacity of most 2D mathematic models is usually calculated directly using a 1D sediment-carrying capacity formula or a modified formula. Since the suspended load particles are small, the sediment particles were divided into finer groups in the algorithm. The median size of each group of particles is directly used to calculate the sediment-carrying capacity for the suspended load group. The most commonly used formula is the one developed by Zhang Ruijin, which is given below (Tan et al., 2018): S ∗ = k ( u 2 + v 2 ghω ) m (28)

In the above equation, k and m are the empirical constants of the model’s sediment-carrying capacity, which were derived based on the experimental data collected from the Yangtze River, the Yellow River, several reservoirs, and indoor water channels. Hence, the empirical constants have a high degree of universality. Theoretically, Zhang Ruijin’s indexless formula can be coordinated with Bagnold’s energy formula; and by improving Zhang Ruijin’s formula, the comprehensive coefficient K 1 can be obtained as shown below: Φ = K 1 γ m γ s − γ m u 2 + v 2 gh ω (29) K 1 = 0.102 ( u 2 + v 2 gh ω ) 0.6346 − 0.21611 lg ( u 2 + v 2 gh ω ) (30)

In the above formula, γ m and γ s are the bulk densities of clear water and sediment, respectively; and Φ is the sediment-carrying capacity, which is equivalent to the S* term in Eq. 28.

In flash floods, the finest and coarsest sediments can differ by more than 100 times, resulting in differences of thousands of times for the settling velocity. Therefore, the mathematical model should calculate the sediment-carrying capacity per group. The sediment-carrying capacity needs to be considered in the transport of sediments suspended in the water and the erosion of the river bed. When the sediment-carrying capacity is greater than the sediment content in the water, the bed surface is eroded. The sediment movement should be considered when calculating the erosion sediment-carrying capacity. The sediment can only be scoured when the flood flow rate exceeds the sediment starting velocity. Hence, Eq. 29 can be improved as follows: Φ ′ = K 1 γ m γ s − γ m ( u 2 + v 2 − U c ) ( u 2 + v 2 ) gh ω (31)

In the above formula, U c is the initial particle size for the group. The formula from the School of Water Resources and Hydropower Engineering of Wuhan University was applied as follows (Tan et al., 2018): U c = ( D d i ) 0.14 ( 17.6 γ s − γ m γ m d i + 0.000000605 10 + D d i 0.72 ) 1 2 (32)

Therefore, the total carrying rate ϕ i of the grouped sediment group is as follows (taking the i-th sediment group as an example): ϕ i = ϕ i s + ϕ i e = σ P i Φ ( ω s ) + ( 1 − σ ) p bi ′ ϕ ′ ( ω i ) (33) σ = min [ S / Φ ( ω s ) ] (34) ω s = P i ω i (35)

In the above formula, ϕ i s is the sediment-carrying capacity of the water flow used for transporting sediment in the water; ϕ i e is the sediment-carrying capacity of the water flow used for activating the sediment bed; P i is the mass ratio of the suspended sediment with this particle size versus the total suspended sediment content; σ is the weight coefficient; and p b i ' is the effective gradation of the sediment bed.

The key to solving Eq. 33 is to calculate the effective gradation of the sediment bed, which is affected by many factors. First, the scouring energy remaining in the water is fully distributed along the water depth so the height of suspended sediments that can be scoured is limited. When the sediments cannot be suspended, the scouring energy that exists in the local water body cannot be used for scouring. Hence, the small sediment transport invariant used for scouring should be related to the suspended height of the scoured sediments. Moreover, the effective gradation of the sediment bed is also greatly affected by how the coarse bed sediments conceal the finer sediments, the effective action of the current drag along the sediment bed surface, and the action of the shear stress on the bed surface (Yang et al., 2020).

Based on the above factors, the effective gradation of the sediment bed can be expressed as follows: p bi e = η i δ i P i d i ∑ η j δ j P j d j (36)

In the above formula, δ i is the suspended sediment height; and η i is the concealment coefficient. These two values were calculated using the formulas proposed by Lane and Kalinske (1941) and Wu et al. (2000).

The model utilizes the Mayer-Peter formula to calculate the bed load for particle sizes greater than 0.2 mm and the Sharmov’s formula is used to calculate the bed load when particle sizes are less than 0.2 mm.

The Mayer-Peter formula (Chien and Wan, 1999): q bx = 8 ( G − 1 ) G d i 3 U U 2 + V 2 max ( τ ∗ − τ ∗ , c , 0 ) 1.5 (37) q by = 8 ( G − 1 ) G d i 3 V U 2 + V 2 max ( τ ∗ − τ ∗ , c , 0 ) 1.5 (38) τ ∗ = n 2 ( U 2 + V 2 ) 1.5 ( G − 1 ) d i D 1 3 (39)

In the formula, G is the relative weight of sediment to water; τ ∗ is the drag force; τ ∗ , c = 0.047 is the initial drag force; and n is the roughness coefficient.

According to Sharmov’s formula, U O H is critical maximum velocity of the sediment. q bx = 0.95 d i ( U U OH ) 3 ( U − U OH ) ( d i D ) 0.25 (40) q by = 0.95 d i ( V U OH ) 3 ( U − U OH ) ( d i D ) 0.25 (41) U OH = 3.83 d i 1 3 D 1 6 (42)

The river bed deformation calculation includes the elevation changes caused by scouring and deposition of bed and suspended matter. Eq. 26 adopts the time forward difference and integrates the spaces in the grid. Next, the Gaussian deformation was calculated as per the following equation: ∫ Δ S ( 1 − λ ) z b new − z b old Δt ds + ∫ Δ l g → b d l → = ∫ Δ S αω ρ s ( S − S ∗ ) (42)

Namely, z b new = z b old − Δ t Δ S 1 1 − λ ∑ 1 m g ni d l i + Δ t 1 − λ α ω ρ s ( S − S ∗ ) (43)

In the equation, z b old and z b new represent the thickness of the movable layer before and after this step, respectively; Δt is the time step; m is the number of grid interface edges; g n i is the single-width sediment transport rate perpendicular to the element interface, which is positive when pointing outwards, i is the interface serial number, 1∼m; ΔS is the grid; Δt is the time step; d l i is the unit interface length, where i is the unit interface serial number, 1∼m.