Modeling of Breaching Due to Overtopping Flow and Waves Based on Coupled Flow and Sediment Transport

Breaching of earthen or sandy dams/dunes by overtopping flow and waves is a complicated process with strong, unsteady flow, high sediment transport, and rapid bed changes in which the interactions between flow and morphology should not be ignored. This study presents a depth-averaged two-dimensional (2D) coupled flow and sediment transport model to investigate the flow and breaching processes with and without waves. Bed change and variable flow density are included in the flow continuity and momentum equations to consider the impacts of sediment transport. The model adopts the non-equilibrium approach for total-load sediment transport and specifies different repose angles to handle non-cohesive embankment slope avalanching. The equations are solved using an explicit finite volume method on a rectangular grid with the improved Godunov-type central upwind scheme and the nonnegative reconstruction of the water depth method to handle mixed-regime flows near the breach. The model has been tested against two sets of experimental data which show that it well simulates the flow characteristics, bed changes, and sediment transport. It is then applied to analyze flow and morphologic changes by overtopping flow with and without waves. The simulated bed change and breach cross-section shape show a significant difference if waves are considered. Erosion by flow without waves mainly occurs at the breach and is dominated by vertical erosion at the initial stage followed by the lateral erosion. With waves, the flow overtops the entire length of the dune to cause faster erosion along the entire length. Erosion mainly takes place at the upper layer at the initial stage and gradually accelerates as the height of the dune reduces and flow discharge increases, which indicates the simulated results with waves shall be further verified by physical experimental evidence.


Introduction
The failure of dams/levees often lead to catastrophic floods that cause loss of life and damage to property and agriculture. Overtopping flow has been one of the most common reasons for dam/levee failures, especially when heavy rainfall inland or storm surges along the coastline occur [1]. During such extreme events, wind-generated waves may run up the face of the dam/levee and enhance overtopping, which makes the erosion and breaching process more complicated [2,3]. Consequently, understanding the characteristics of flow, sediment transport, and the dam/levee breach process by overtopping flow and waves, and predicting discharge, breach width, depth, and sediment deposition have great significance for the decision-making in flood prevention and risk analysis.
In recent decades, the dam/levee breach problem has been extensively studied in mechanism, experiment, field, and modeling. A good review and summary can be found in the recent publication of the American Society of Civil Engineers/Environmental & Water Resources Institute (ASCE/EWRI) Task Committee on Dam/Levee Breaching [4]. However, the breaching process and its rapid morphologic changes due to overtopping flow and waves are still little understood due to observing and measuring difficulties in either lab experiments or field investigations. Even fewer results are available for the breach characteristics, such as breach width, depth, section shape, and bed changes, downstream at different time stages to examine the influence of waves during overtopping flows.
Two-dimensional (2D) numerical models solving the nonlinear shallow water equations (SWEs), the sediment transport, and bed change equations have been widely used to mathematically describe free surface flows and their associated morphologic changes over complex topography. Those models are often discretized by finite volume methods or finite difference methods and solved by Roe, Harten, Lax, and van Leer (HLL)/HLLC(C stands for Contact) [5], or TVD (total variation diminishing) schemes [6] in order to simulate the flood shock wave. However, simulations of earthen or sandy breaching processes due to overtopping flow and waves are much more difficult than those over fixed beds because highly unsteady flows, high concentrations of sediment, and rapid changes of bed can occur.
In recent years, a considerable number of studies proposed a coupled flow and sediment modeling concept that takes into account the impacts of bed change and variation of sediment concentration on the flow when simulating flows over movable beds [7][8][9][10][11]. These studies suggested that the flow continuity and momentum equations should include these impacts and the flow, sediment transport, and bed change should be calculated in a coupled approach. Meanwhile, non-equilibrium sediment transport approaches have been recommended to solve the cases of transient flows in which strong erosion and deposition occur [7,[10][11][12], since the traditional assumption of local equilibrium capacity of sediment transport could be problematic. It might be inappropriate to use the sediment entrainment or transport capacity functions developed based on uniform flow conditions for embankment breach flow situations [7,10,13]. Furthermore, the bank erosion or collapse calculation should also be included in the modeling frame to deal with the instability of slopes [13].
To more accurately simulate the breaching and sediment transport processes and analyze their features due to overtopping flow and wave overwash, an improved cell-centered finite volume model based on the Godunov-type central upwind scheme [14] is developed to simulate shallow water flow and sediment transport over arbitrary topography. The model uses a non-equilibrium sediment transport approach and takes into account the effects of sediment concentration and bed change in the flow governing equations. The model solves the flow, sediment transport, and bed change equations in a coupled way based on the structured mesh.
The governing equations and the numerical methods, including the central upwind scheme and nonnegative water depth reconstruction, are presented in Sections 2 and 3, respectively. The model is then tested in Section 4 against measurements and afterwards applied to investigate the flow and erosion processes in the cases of overtopping flow with and without waves in Section 5.

Governing Equations of Flow and Sediment Transport in Breaching Process
During dam/levee breach, the sediment concentrations are high and beds change rapidly. Thus, the interactions between flow, sediment transport, and channel bed cannot be ignored. We take into account the bed changes and density variability of flow and sediment mixture in mass and momentum conservation equations based on the non-equilibrium concept [7,15], which assumes that sediment cannot reach new equilibrium states instantaneously due to the temporal and spatial lags between flow and sediment transport. This assumption is more realistic, especially in cases of strong erosion and deposition. Therefore, considering that the flow density is variable in the vertical direction, the depth-averaged 2D equations governing flow and total load sediment transport processes due to overtopping flow can be expressed as [7] ( ) ( ) ( ) where t is the time; x and y are the longitudinal and lateral coordinates; h is the flow depth; u and v are the flow velocities in x and y directions, 2 2 U u v = + ; b z t ∂ ∂ denotes the rate of change of b z (bed surface elevation above a reference datum); s z is the water level, s z = b z +h; n is the Manning roughness coefficient; g is the gravitational acceleration; t C is the actual volumetric total-load sediment concentration; ρ is the density of the water and sediment mixture in the water column determined by ( ) m z x z y = + ∂ ∂ + ∂ ∂ (considers the lateral erosion along the side slope); t q * is the total-load sediment transport capacity, also called the equilibrium transport rate per unit channel width. Lt is the non-equilibrium adaptation length or saturation length of sediment transport [15][16][17]. The total-load sediment transport capacity and the adaptation length are determined by the same method proposed by Wu [15]. Using the relations ( ) z +h, and Equation (4), the flow density on the left-hand sides of Equations (1)-(3) can be separated from the variables. Thus, Equations (1)-(4) can be rewritten in the following compact form as [15] ( ) where Φ , ( ) F Φ , and ( ) G Φ represent the vectors of unknown variables, fluxes, and source terms, respectively: The bed deformation is determined with Moreover, the sediment transport capacity is determined using the formula of Wu et al. [18] and the effect of gravity on sediment transport over steep slopes is also considered [19]. The settling velocity s ω used in Wu et al.'s [18] formula is determined by the Richardson [20]

Model Discretization
The mathematical model is discretized using a finite volume method based on a rectangular mesh, as shown in Figure 1, in which a non-staggered (collocated) grid system is used. The primary variables h, u, v, and t C , and the bed elevation are defined at cell centers, representing the average values over each cell. The fluxes are calculated at the cell faces. Integrating Equation (5) over the (i, j) control volume and applying the Euler scheme, Equation (5) can be discretized using Green's theorem as follows: where Δt is the time step length, are the cell lengths in x and y directions, 1/2, is the , and , 1/2 n i j+ G is the flux at face ( , 1/ 2) i j + . Si, j represents the source terms evaluated at the cell center. The determination of the fluxes at the cell faces is one of most important procedures to solve Equation (9). Thus, this paper develops an improved mass-balanced 2D shallow water and sediment transport model based on the Godunov-type central upwind scheme [14] and the nonnegative reconstruction of the water depth method [21].

Nonnegative Reconstruction of Riemann State
A correct reconstruction of the Riemann states [6] at the interfaces is the key to determine the interface fluxes accurately. Here the nonnegative reconstruction of the water depth method [21] is adopted to calculate the Riemann states, which have been proved to be accurate to deal with wetting and drying. The Riemann states, state variables, and fluxes defined at the interface can be obtained from the available flow information at the cell center by using a slope-limited linear reconstruction, e.g., left-hand side values at the interface (i + 1/2, j) are computed by and ((hv)x)i, j represent the numerical derivatives of the flow variable at cell (i, j). In order to preserve the stability of the model, a one-parameter family of generalized min mod limiters is used to calculate the numerical derivatives. For example, where θ is the parameter to control numerical viscosity. In general, larger values of θ make the numerical result less dissipative and more oscillatory [6]. θ = 1.3, suggested by Kurganov and Petrova [22], is used in this study. The min mod function is defined as The value of (hx)i, j, ((hu)x)i, j, and ((hv)x)i, j can be obtained in a similar way. Similarly, the right-hand side values at the interface (i + 1/2, j) can be calculated by To obtain the final Riemann states, a single value of the bed elevation suggested by Audusse et al. [23] should be defined as Then the face values of water depth are redefined as Obviously, this method ensures the reconstructed water depth to be nonnegative. However, the reconstructed water depth may be very small or zero and thus lead to large values of the flow velocity, which may cause the model to be unstable. In order to prevent this problem, the regularization technique suggested by Kurganov and Petrova [22] is used to calculate the corresponding velocity components at the interface: where ε is an empirically pre-defined positive number (ε = 10 −6 m in all our simulations). Hence, the Riemann states of the other flow variables can be recalculated by The procedure above obviously does not affect the well-mass-balanced property of the numerical scheme for wet-bed applications. However, the local bed modifications are needed to get the well-balanced solution for a dry-bed application. As shown in Figure 2, the wet cell (i, j) shares the common edge (i + 1/2, j) with the dry cell (i + 1, j), and the bed elevation of the dry cell is higher than the water surface elevation at cell (i, j). The problem of stationary water at the steady state is considered through u ≡ 0, v ≡ 0, and η ≡ const at the wet area. After the aforementioned reconstruction at the cell interface (i + 1/2, j), the components of the Riemann states are given by Therefore, the fluxes at cell interface (i + 1/2, j) are calculated using the bed elevation (zb)i+1/2, j rather than the actual water surface elevation η in cell (i, j). However, fluxes at the cell interface (i − 1/2, j) are estimated using the actual water surface elevation η in cell (i, j), which produces a net flux in cell (i, j) and violates the well-balanced property of the scheme. The local bed modification is thus implemented to overcome this difficulty following Liang [21]. As show in Figure 2, the difference between the actual and fake water at interface (i + 1/2, j) is calculated by The bed elevation and water surface are locally modified by subtracting Δz from the original reconstructed value: Obviously, = η = η = η after the local bed modification and no net fluxes are produced in cell (i, j), which ensures the scheme to be well-balanced for the dry-bed application.

Central Upwind Scheme
Once the left and right variables are estimated using the aforementioned reconstruction method, the interface flux can be calculated using the Godunov-type central upwind scheme as follows [14]: where 1/2, are the reconstructed Riemann states at the left-and right-hand side of interface

Treatments of the Source Terms
The source terms, as shown in Equation (7), can be split into the bed slope terms and friction terms. It is important to discretize the bed slope terms appropriately to ensure the scheme to be well-balanced. Hence, the bed slope terms are approximated as follows [7]: Herein the values of bed elevations are obtained after considering the local bed modification in Equation (19).
In general, a very small water depth may cause numerical instability when a simple, explicit discretization of the friction terms is used. To overcome this problem, implicit treatments [24] and semi-implicit treatments [25] are widely used. In this study, the friction terms are discretized using the semi-implicit treatment [7] given by: where the superscript k represents the kth time level. Substituting Equations (29) and (30) into the time marching Equation (7), the values of the flow variables at the end of the time step can be obtained.

Calculation of Sediment Transport
The total-load sediment transport Equation (4) is solved using Equation (9), while the sediment flux at the cell faces is determined using the upwind scheme (e.g., flux at interface (i + 1/2, j)): The bed deformation in Equation (8) is discretized as In order to avoid instability of sediment transport, flow velocities at the cell centers are evaluated from the intercell fluxes as presented by Wu et al. [7] to determine the bed shear stress and sediment transport capacity t q * .
The non-cohesive slope avalanching algorithm presented by Wu [15] is modified by specifying different angles of repose for the submerged materials and the dry materials above the water surface. The submerged repose angle is about 33 degrees as normally used, while the dry repose angle is about 79 degrees, as suggested by Wu [15], for the test cases presented in Section 4.

Stability Criterion and Boundary Conditions
The current numerical scheme is explicit, and its stability is governed by the Courant-Friedrichs-Lewy (CFL) criterion and additional conditions for sediment transport and bed change computations, as described in the following. In order to preserve the positivity of water depth at each time step, the Courant number has to be kept less than 0.25 [26], reading where NCFL represents the Courant number, and a and b are given by: Moreover, the bed change at each time step should be less than about 10% of the local flow depth [7]. Thus, the CFL restriction for the current model is given as  [27,28]. For wall boundary conditions, the velocity normal to the boundary and the water surface gradient are both set to zero at the boundary. It was shown that this procedure is second-order accurate [29].

Case 1: Tsunami Run-Up onto a Complex Three-Dimensional (3D) Beach
This laboratory experiment was suggested as a benchmark test for numerical models in The Third International Workshop on Long Wave Run-up Models in 2004. It was performed in a 1:400 scale wave tank, which approximated the coastline topography near Monai in Japan. This area was inundated by the 1993 Okushiri tsunami and a maximum runup of 31.7 m was observed during this event. As shown in Figure 3, the tank was 5.5 m long and 3.4 m wide. The incoming wave was generated by the movements of wave paddles and three gauges, Ch5, Ch7, and Ch9 (see Figure 3), were set up to record the time history of water elevation.
For this case, the computational domain is approximated made of 392 × 243 grids with a uniform grid size of 0.014 m. At the beginning, the tank is filled with still water of 0.135 m depth. The left boundary is the wave inflow boundary, in which the time history of the measured surface elevation was specified. The other three side boundaries are set to be closed. The total simulation time is 22.5 s and the Manning coefficient is n = 0.0025, which was suggested by [30]. The computed surface elevation at different times is plotted in Figure 4, while Figure 5 shows the computed water elevation at three gauges compared with the observation data. In general, the movement of the wave is well predicted by the present model. Some discrepancies are found between numerical and experimental results due to the fact that the three-dimensionality of the flow cannot be exactly captured by a depth-averaged 2D model. Nevertheless, the lead wave height (with relative errors less than 2.3%) and arrival times (with relative errors less than 0.5%) are predicted with a good accuracy at these gauges.  (e) (f)

Case 2: One-Dimensional Dam Break Flow over Moveable Bed
The experiment of a one-dimensional dam break flow over a movable bed conducted by Capart and Young [12] was used to test the capability of simulating the sediment erosion and bed changes using the present model. As shown in Figure 6, the experimental flume was 1.3 m in length, 0.2 m in width, and 0.7 m in depth and covered by sediment particles with a depth of 0.05-0.06 m. In this experiment, the sediment particles were light artificial spherical pearls covered with a shiny white coating. The uniform diameter of the particles was 6.1 mm with a density of 1.048 kg/m 3 and a porosity of 0.28. The initial water depth in the upstream was 0.1 m and the downstream was dry. At the beginning of the experiment, the gate was instantly removed to release the water, which generated flood propagation and caused erosion downstream. The present model was used to simulate this process. In the simulation, the computational grid was 0.005 m and the Manning's coefficient was 0.025 sm −1/3 . The adaptation length Lb and adaptation coefficient α were 0.2 m and 2.0, respectively. The settling velocity of sand used in the model was 7.6 cm/s based on the experimental data. The time step was adapted by Equation (37). Figure 7 compares the numerical and experimental water surface and bed elevation at different time steps. As can be seen, the maximum erosion occurred downstream close to the gate, and the eroded area expanded downstream through flood propagation. A hydraulic jump was also observed near the gate of the dam. The erosion magnitudes and wavefront locations are well predicted by the present model. However, the location of the hydraulic jump is less accurately predicted. The main reason can be attributed to the present depth-averaged 2D model, which cannot precisely describe the three-dimensionality of the dam break flow over the movable bed. Overall, the simulated water depth and bed elevation agree well with the measured values. The relative errors of the bed elevation are less than 5.6%. Figure 8 shows the simulated sediment transport rate at different times. The maximum sediment concentration appeared at the front of the flood, where the high velocity can entrain more sediment to move with the water.

Numerical Investigation of 2D Dam Breaching Processes with Overtopping Flow and Waves
The present model has been applied to simulate dam breaching experiments with and without waves in order to analyze the different features of sediment transport and erosion during different overtopping processes.

Breaching Processes Due to Overtopping Flow
The experiment considering flow overtopping without waves was carried out by HR Wallingford in the UK [31]. This non-wave case represents the overflow and erosion due to water level rises. The experiment was conducted in a flume with a dimension of 50 m in length and 10 m in width. The sandy dam, located upstream about 36 m towards the flume entrance, had a height of 0.5 m and an initial crest of 0.2 m in width and 0.02 m in depth at the top. The dam was made of non-cohesive sand with a median diameter of 0.25 mm. Both upstream and downstream slopes were 1:1.7. The initial water level was 0.5 m during the experiment, and the flow then overtopped the dam through the initial breach.
A rectangular grid of 1000 × 200 nodes was used in the simulation. The porosity of the dam was set to 0.36 based on the measurements, and the Manning's coefficient was 0.018. The initial water surface used in the simulation was 0.5 m. In the simulation, the adaptation length was 0.05 m, which was close to the grid space to avoid anti-dune formation in the numerical results [15]. The time step was determined by Equation (37). The inflow discharge of 0.07 m 3 /s was specified at the flume entrance, while a free overall flow condition was applied at the end of the flume. Figure 9 compares the measured and simulated breach flow discharges varying with time. The simulation results generally agreed well with the measured data, with a relative error of peak discharge of less than 8.2%, which shows that the model can provide a good prediction for such a complicated breaching process. Figure 10 shows the bed changes due to flow overtopping at the cross-section of the breach. The erosion at the initial stage was dominated by the vertical erosion and then changed to the lateral erosion, which means the breach became deeper at the beginning, and then became wider.

Breaching Processes Due to Wave Overtopping
To investigate the different behaviors of sediment transport and erosion by wave overtopping, a hypothetical numerical analysis of wave overtopping is carried out in this section for comparison. Unlike the water level steadily rising from the still water surface in the no-wave case above, this case assumes a continuous wave reaches the dam from the upstream. Other than this, the numerical analysis in this section keeps the same conditions as the previous experimental case. In the simulation, a regular wave with a wave height of 0.02 m and period of 1.5 s was applied at the upstream, while the initial water surface was 0.49 m. Therefore, the wave runup can overtop the dam. The computational mesh and other parameters used in the simulation were the same as the case of overtopping flow above, while the wave-adjusted boundary condition was used [32,33]. Figure 12 shows the simulated bed elevations at different time steps during the wave overtopping process. As can be seen, the flow overtopped the entire length of the sandy dune due to the wave attack and, correspondingly, erosion took place along the entire sandy dune in comparison to the breach erosion in Figure 11. Meanwhile, the erosion process was much faster than the case without waves and the initial breach at the top of the dam had little influence on the flow and erosion processes.  Figure 12. Bed changes and water depths with a wave at different times. The simulated water depths are plotted in the colored contours in the 3D topography view: (a) bed elevation and water depth at 22.8s, the sandy dam was eroded along the entire length and sediment was deposited in the downstream; (b) bed elevation and water depth at 28.2 s, the height of the dam was reduced quickly and erosion in the lower part of the dam downslope was increased.

Discussion of Results
To further see the different morphological changes, Figure 13 compares the bed changes at the middle cross-section along the flow direction during both flow and wave overtopping processes. As can be seen in Figure 13a, the downstream slope was eroded first and then the erosion expanded to the upper part during the flow overtopping. However, during the wave overtopping, as shown in Figure  13b, the erosion mainly occurred at the upper part at the initial stage because the wave ran up and overtopped the dune to cause a high flow velocity near the top of the dam. As the height of the dune reduced and the overtopping flow discharge increased, the entire dune was eroded fast as well. This process was significantly different in comparison to the one without the wave since the flow can only pass through the breach as shown in Figure 11.

Conclusions
A depth-averaged 2D coupled flow and sediment model is developed to investigate the earthen dam/levee breach processes due to overtopping flows with and without waves. The effects of sediment concentration and bed change on the flow are considered in the continuity and momentum equations, while the non-equilibrium transport of the total load is used for the sediment transport and bed changes. A specially designed avalanching algorithm is adopted in this model to consider different repose angles for the submerged materials and the dry materials above the water surface. The flow, sediment transport, and bed change equations are solved by the explicit algorithm based on the finite volume method. The model improves the Godunov-type central upwind scheme with nonnegative water depth reconstruction and a slope-limited linear reconstruction to efficiently handle the mixed-regime flows near the breach and the wetting-drying problem. To enhance the performance of the model, a varying time step is used that satisfies the CFL condition and ensures that the bed change at each time step is less than 10% of the local flow depth. The developed model has been tested using two sets of experimental data. The model predicts the water levels, bed changes, breach width, flow patterns near the breach, and the breach flow hydrograph generally well.
The numerical investigation of the breaching process due to overtopping flow with and without waves shows that waves have significant impacts on flow field, breaching width, depth, and the cross-section shape. In the case without waves, the flow through the breach eroded the downstream slope first. The erosion at the breach was dominated by the vertical erosion at the initial stage and then changed to the lateral erosion. In the wave overtopping case, the water overtopped the entire length of the dam due to the wave runup and caused the erosion along the entire length. Meanwhile, the erosion mainly occurred at the upper part at the initial stage. As the height of the dune reduced and flow discharge increased, erosion along the entire dune was faster than that in the no-waves case. However, it should be noted that the numerical results with waves are hypothetical and should be further verified by experiments or filed observations.