Application of Bed Load Formulations for Dam Failure and Overtopping

The Enhanced HLLC scheme as a robust approximate Riemann solver is used for numerical modeling of three different test cases of mobile bed and stepped mobile bed in dam failure and dam overtopping conditions. The current research has been done in the frame of the finite volume method using shallow water equations along with the Exner equation for sediment continuity. The Ribberink, Wong and Parker formulations have been used for the modelling of bed load movement. A convenient approach based on the Boussinesq hypothesis is deployed for considering turbulence effects in the second case. The affections of stepped and slope condition for the flow bed are considered through a corrected version of the HLLC flux components. Finally, the model is applied for modelling overtopping in the third case. The results of the present model are relatively reasonable by comparing with the experimental data.


Introduction
All kinds of man-made structures such as dam might cause lots of risks depending on the dam location and operating condition. Therefore, Dam failure and overtopping risks are crucial challenges because of their huge effects on the downstream domain. There are three main methods for evaluating of dam failure and overtopping risks containing; studying of the last similar experiences, physical modeling and mathematical modeling. One of the most reasonable ways for studying the mentioned risks is the numerical modelling of waves induced by dam failure and also bed load movement or soil erosion. The numerical results include water surface and bed profiles in the various times after dam breaking or overtopping.
In the recent decades, the finite volume method has been known as one of robust discrimination methods for partial differential equations. After introducing of Godunov's scheme, some considerable activities have been done about waves modeling induced by dam failure and solving Riemann problem. One of the mentioned activities has been done by Harten et al. [1]. They proposed a Riemann solver assuming two separate waves known as HLL Riemann solver. Toro et al. presented a modification of HLL scheme named the HLLC Riemann solver assuming three separated waves considering the shear waves affection [2].
In order to extend HLL scheme to movable bed cases, Fraccarollo et al. studied rapid erosion in a rectangular flume using LHLL scheme as a Godunov-type scheme [3]. Simpson and Castelltort applied modified HLL scheme for solving shallow water equations and sediment continuity equation in the coupled approach using two test cases related to a dam break [4]. The stepped bed condition is one of the important subjects which has been investigated by Vásquez and Leal. They worked on the implementing and modeling of a dam break test case on the stepped mobile bed condition [5]. Castro et al. investigated about one and two-layer flows and also studied the flooding caused by River Mero near the Cecebre dam using the finite volume method and streaming SIMD extensions [6]. Goutier et al. introduced new approximations for the eigenvalues of Saint-Venant and Exner equations, they used HLLC solver in a test case including sloped bed flume under mobile bed conditions [7].
The subject of sediment and water interaction and its affection on wave speed estimations in dam-break problems have been studied by Soares-Frazão and Zech; they worked on a new approximation method aimed at calculating of wave's speeds. Also, they considered the properties of bed material in their research [8]. Then Soares-Frazão et al. worked on the partial dam failure using a long flume under mobile bed conditions; they compared the results of the numerical modelling and experimental test. They also compared the numerical results of twelve different work teams for the same test case [9]. In another study, Talukdar et al. compared a number of significant approaches for sediment movement and related bed load formulations. They investigated about the rate of deviancy between numerical results and experimental data resulted from the various bed load formulations [10]. Liu et al. worked on another viewpoint of dam-break modeling. They represented a two dimensional model for reproducing dam-break flow over wet-dry fronts on irregular topography using a central-upwind scheme and unstructured triangular grids [11].
Also, Liu et al. applied Godunov-type central-upwind scheme for solving shallow water equations and Exner equation in order to model flows over the erodible bed. They used grass bed load formulation to calculate the morphological process [12].
As one of the recent works, Caviedes-Voullième et al. represented a numerical strategy for modeling flow over movable bed conditions based on shallow water equations and Exner equation. Their research presented an applicable approach for sewer flushing analysis [13].
In view of investigating dam overtopping, Schomocker and Hager presented the experimental results of two hydraulic models related to overtopping and breaching dikes [18]. Also, Van Emelen et al. studied about using of various formulations impacts in the dam breaching model [14]. Zhignuo et al. worked on breaching model of an earth dam due to overtopping flow [15]. Also, Saberi and Zenz investigated on the one dimensional and two-dimensional modelling for dam failure due to overtopping flow under various bed slope conditions. They used the revised MacCormack scheme solving the Saint-Venant equation and Exner equation [16].
For modeling of the flow and bed changes and also the erosion of a dike due to overtopping, the shallow water equations along with sediment continuity equation (Exner equation) are solved in a coupled approach in the present research. According to the mentioned background, it is obvious that improvement and extending of shallow water models as robust tools using simple and accurate, applicable methods such as HLLC scheme are still going on, especially in the case of two-phase flows such as the flood-induced by dam failure over movable beds and soil erosion due to dam overtopping. Despite the mentioned researches, it is very important to incorporate the advantages of the last studies and subjoin them the simple turbulence method and also a suitable method considering the stepped and slope bed conditions. Considering sediment properties for wave speed estimation and improving the accuracy of numerical results, the novel wave speed estimation introduced by Soares-Frazão and Zech [8], applied in the current research.

Governing Equations
The shallow water equations for water flow continuity and momentum along with the Exner equation for bed load continuity are used in the vector form as follows [17]:  Where is the time, ℎ the water depth, and the depth-averaged velocity components in and directions, the gravitational acceleration, and the unit discharge components of the flow, the bed elevation, 0 the bed porosity and , and , the components of the sediment discharges per unit width which are calculated based on Ribberink [7] formulation and Wong and Parker [18] formulation and also parker [6] formulation respectively shown in the Equations 2 to 4 as follows: Where is the relative density of sediment, 50 is a representative grain diameter, * denotes the bed shear stress and * the critical bed shear stress ( * = 0.047). The bed shear stress is calculated as follows: and n is the Manning friction coefficient. The source term vector is calculated as follows [8]: Where 0 and 0 are the bed slope in the x and y directions, respectively, and and are the friction slope components.

The Turbulence Stresses
The Boussinesq turbulence stresses , , , xx xy yx yy     are as follows [21]: Where is the eddy viscosity coefficient, which is calculated by the following relation and replacing Chezy coefficient with Manning roughness [5]: Where n and R are the Manning roughness and hydraulic radius, respectively. The eigenvalues which have been validated by Soares-Frazão and Zech [8] are applied for estimating wave speeds in all test cases.

Modified HLLC Scheme
The HLLC scheme was introduced by Toro [22] is used with some modifications accounting for bed load properties and its effects on wave speed estimation and completing the original HLLC flux with the fourth component as sediment discharge which presents sediment transport rate through a cell interface. The modified HLLC intercell flux is as follows [8]: Where * and * are mass and momentum fluxes in the normal direction respectively, and * is momentum flux in the transverse direction [8]. Therefore the flux expressions are evaluated as follows [22]: Where is the water level and + and − are the maximum and minimum wave speeds at the right and left sides of the cell interfaces.
Finally, the sediment flux in the normal direction is: The waves speed + and − are related to the minimum and maximum wave speeds corresponding to the sediment movement at the interface [22].

Consideration for Stepped Bed Conditions
For extending the method that has been applied by Goutière et al.
[7] method a correction has been used in the second component of the flux as follows: Where g is the gravity acceleration and ℎ , ℎ , and are the water depth on the right and left-hand sides and bed levels in the right and left-hand sides of the interfaces of the cells, respectively.
The system of shallow water and Exner equations are solved by the finite volume method on triangular meshes so the discretization of the Equation 1. is calculated as follows [17]: Where * ( ̅ ) is the flux related to the turbulence stresses and * ( ̅ ) is HLLC-WAF flux and is the velocity vector in the j th interface and represents the rotation vector.

Boundary Conditions
The computational domain is assumed which discretized by N cells, so that cells = 1, … , lie within the computational domain. Also, the MTh cell is assumed to be the boundary cell. For applying boundary conditions one fictitious cell next to each boundary is needed. For the left boundary the fictitious cell is denoted by i = M-1 and for the right boundary, it is denoted by = + 1. Two types of boundary conditions are discussed [22].
Transmissive boundary conditions are given by following relations: Reflective boundary conditions are given by following relations:   (19) Where h is the water depth, u is the velocity component in the or direction considering boundary direction, ℎ is the sediment depth [22] In the current research, the target domain was discretized by triangular meshes as the first step, then for each intersection of a triangular mesh, the shallow water equations conjuncted to the Exner equation have been solved using the HLLC Riemann solver for obtaining Fluxes. Furthermore, extra terms for each intersection were added exerting turbulence effects along with improvement accounting for the stepped bed condition. Finally, all equations have been solved for upper time step considering boundary conditions in the framework of the finite volume method.

Application of the Present Model
For validating the present model, three different test cases including movable, stepped bed and overtopping conditions are chosen as follows:

The First Test Case: Ideal Dam Break Over Movable Bed
This test has been conducted by Spinewine and Zech [23]. A rectangular and horizontal flume was used in the test case. The length of the flume is 2.5 and its width is 0.1 also the height of the sidewall is 0.35 . The test initial conditions, are presented in Figure 1. A fully saturated sediment layer whose thickness is 5 extended all over the flume bed. A sluice gate is installed in the middle of the flume as a simulated dam. The gate opening time is predicted to be less than 0.05 second. The simulated sediment layer is composed of small PVC pellets whose properties are presented in Table 1. Where is relative density, is mean diameter of cylindrical pellets, w is sediment settling velocity.

Figure 1. Initial conditions for the second test case [23]
The shape of the bed material is cylindrical whose diameter and height are 3.2 and 2.8 , respectively, so their diameters have been assumed 3.5 approximately [23]. Also, a weir whose crest level is the top level of the bed layer is installed on the downstream end of the flume. The water height on the upstream side of the gate is 0.1 . The Wong and Parker [19] formulation has been applied for the modelling of bed load movement in the current test case. According to the results presented in Figure 2, in view of bed erosion and deposition in the different times of the test, the simulated bed levels near = 1.25 are in good convergence with experimental data of bed levels but there is a little difference in the simulation of erosion. In view of sediment deposition, the maximum levels of the bed rising in all time of the test approximately are well simulated.

Sensitivity Analysis
Sensitivity analysis has been done in the present test case. For example, the computational domain is divided into 15372, 23045, 30734 triangular cells. For evaluating the numerical results at the = 1 in view of water level simulation, 0.015 , 0.012 and 0.011 are obtained as the correspondent RMSEs after running the model. Therefore, considering the time of running the model and accuracy simultaneously, finally the domain area is divided into 23045 triangular cells in the first test case.

Model Results
In view of water surface simulation, there are a little diversity between the numerical results and experimental data in all four graphs. The total shapes of water surface are in good convergence with the experimental data. The Root-Mean-Square Errors (RMSE) of the numerical results compared to experimental data are represented in Table 2.

The Second Test Case: Ideal Dam Break Over Stepped Bed
The initial conditions of this test represented in Figure 3, includes a horizontal flume with the stepped mobile bed condition. The flume length is 19.2 and its width and height are 0.5 and 0.7 , respectively [24]. A gate as a simulated dam is installed in the middle of the flume. The opening time of the gate is between 0.1 and 0.2 . The initial conditions of the two selected tests which have been conducted in the same mentioned flume named TS.25 and TS.28 are presented in Table 3 [24]. The bed covered with sand whose properties presented in Table 4 as follows:

Sensitivity Analysis
Similar to the test case1, sensitivity analysis is done in the present test case. The computational domain is divided into 17632, 24292 and 34169 triangular cells. At the = 1 in view of water level simulation, 0.041 , 0.036 and 0.035 are obtained as the correspondent RMSEs after running the model. Therefore, considering the time of running the model and accuracy of the results simultaneously, finally the domain is divided into 24292 triangles in the second test case.

Model Results
The Ribberink [18] formulation has been applied for the modelling of bed load movement in the current test case. Comparing the present model results and the experimental data of Leal et al. [24] for Ts.25 is illustrated in Figure 4a and 4b at 1 and 4 seconds. The good agreement between the numerical results and experimental data are obtained in view of the water surface simulation. Although the water surface resulted from the present model is a bit lower than the experimental data of Leal et al. [24]. In view of bed level simulation, there are some differences near bed step and the general form of bed simulation is acceptable. Comparing of the present model results and experimental data of Leal et al. [24] for Ts.28 which are illustrated in Figure 4c and 4d at 1 and 4 seconds after dam break. It can be observed, the model results are in good agreement with the experimental data both in view of water and bed levels. Also, as it is illustrated in Figure 4d, the bed levels difference near = 10 at the 4 th second is regenerated well with a longitudinal lag of 1 . The RMSEs of the present model results compared to the experimental data presented in Table 5.

The Third Test Case: Earth Dam Erosion Due to Overtopping
The selected test includes the dike erosion, which conducted by Schmocker and Hager [25] in the rectangular channel whose length, width and height are 8, 0.4 and 0.7 , respectively. The bottom of the channel whose section represented in Figure 5. is horizontal and uniform inflow is generated at the channel intake. The width of the channel was adjusted to 0.2 . The dike in the selected test was located one meter from the channel intake. The shape of the dike was trapezoidal and it has been built with the non-cohesive material without any core. Seepage through the dike is impossible because of installing bottom drainage. The channel is divided to 1042 triangular mesh represented in Figure 6. The constant parameters of the test are as follows: Dike height: = 0.2 , Dike width: = 0.2 , the length of crest: = 0.1 , Dike slope: 0 = 1: 2, the length of Dike: = 0.9 . Also The dike material properties presented in Table 6 as follows:

Sensitivity Analysis
The computational domain is divided into 870, 1042, 1472 triangles, after running the model, 0.007 , 0.005 and 0.005 are obtained as the correspondent RMSEs for water level simulation at the = 100 . According to the results of sensitivity analysis, 1042 triangular cells are selected in this test.

Model Results
The parker [20] formulation has been used to model soil movement along with Goutiere et al.
[7] correction for HLLC Riemann solver in the stepped bed condition which has been extended to the current test case condition according to dike upstream and downstream slopes.
As it is shown in Figure 7, Considering the first order accuracy of enhanced HLLC Riemann solver which has been employed in the present model along with the simple and convenient Boussinesq turbulence model, the numerical results are acceptable. According to the calculated RMSE shown in Table 7 especially at the times of = 50 and = 100 (Figures 7e  and 7f), the simulation of the dike erosion has been successful.

Conclusion
In the current research three different laboratory tests in the various conditions such as mobile bed, stepped mobile bed and overtopping over a dike, have been simulated by the enhanced HLLC Riemann solver in the frame of the finite volume method. Three different bed load formulations named Wong and Parker [19] formulation and Ribberink [18] formulation and Parker [20] formulation have been used for the first, second and third test cases, respectively. According to the results of the current research, using a convenient turbulence model besides an extended and applicable correction method for the stepped and sloped bed conditions has been successful. Considering the first order accuracy of the original HLLC scheme, the applied enhancements for the mentioned scheme make it more robust and effective at least for three test cases which have been applied in the current research. The convergence between the numerical results and experimental data both in view of the bed level simulation and water surface simulation are relatively reasonable according to the mentioned discussion in the current research and with due attention to the calculated RMSEs.