A 2D Hydraulic Simulation Model Including Dynamic Piping and Overtopping Dambreach

: Numerical simulation of unsteady free surface ﬂows using depth averaged equations that consider the presence of initial discontinuities has been often reported for situations dealing with dam break ﬂow. The usual approach is to assume a sudden removal of the gate at the dam location. Additionally, in order to prevent any kind of dam risk in earthen dams, it is very interesting to include the possibility of a progressive dam breach leading to dam overtopping or dam piping so that predictive hydraulic models beneﬁt the global analysis of the water ﬂow. On the other hand, when considering a realistic large domain with complex topography, a ﬁne spatial discretization is mandatory. Hence, the number of grid cells is usually very large and, therefore, it is necessary to use parallelization techniques for the calculation, with the use of Graphic Processing Units (GPU) being one of the most efﬁcient, due to the leveraging of thousands of processors within a single device. The aim of the present work is to describe an efﬁcient GPU-based 2D shallow water ﬂow solver (RiverFlow2D-GPU) supplied with the formulation of internal boundary conditions to represent dynamic dam failure processes. The results obtained indicate that it is able to develop a transient ﬂow analysis taking into account several scenarios. The efﬁciency of the model is proven in two complex domains, leading to > 76 × faster simulations compared with the traditional CPU computation.


Introduction
It is broadly accepted that, even though dams are hydraulic structures that can provide benefits, there exists a risk associated with their failure.When dams fail, the consequences can vary with the extent of the inundation area and the size of the population at risk, but they could be devastating.The possibility to simulate scenarios of dam breach events and the resulting floods with the help of computational tools is essential.It can help to anticipate and reduce threats due to potential dam failure and justifies the search for improved breach analysis tools considering different failure mechanisms.
Among the tools available today for analysing dam failures and their resulting outflow hydrograph, some of the best known are the National Weather Service (NWS) Dam-Break Flood Forecasting Model (DAMBRK), the NWS Simplified Dam-Break Flood Forecasting Model SMPDBK [1] and the NWS FLDWAV [2].The focus is the prediction of the reservoir outflow hydrograph and its routing downstream the dam.The hydrograph depends on the breach characteristics and the reservoir shape and size.Concerning the numerical simulation of the downstream wave propagation given the breach hydrograph, models dealing with unsteady free surface flows using depth averaged equations that consider the presence of initial discontinuities are the most recent [3][4][5][6].It is usual to assume a sudden dam removal but, in order to deal with earthen dams, it is very interesting to include the possibility of a progressive dam breach that can lead to dam over-topping or dam piping.In this scenario, the predictive hydraulic model offers a global analysis of the water flow in dynamic conditions including warning time as sum of the breach initiation time, breach formation time and flood wave travel [7].
It is therefore justified to devote time and effort to the development of accurate and complete advanced simulation models able to incorporate the dynamic breach formation in the context of ambitious wave propagation computations.Recent advances in this sense have led to a new generation of tools based on finite volume techniques for the 2D shallow water equations that can be used for real time flow predictions [8].When considering a realistic large domain with complex topography, a fine spatial discretization is mandatory.Hence, the number of grid cells is usually very large and, therefore, it is necessary to use parallelization techniques for the calculation, with the use of Graphic Processing Units (GPU) being one of the most efficient due to the leveraging of thousands of processors within a single device [9].As the computational efficiency of the numerical method used is limited by the stability restriction, a CUDA-based implementation on GPU is used.This is the core tool of the RiverFlow2D-GPU software developed by the authors and distributed by Hydronia LLC [10][11][12][13].It is able to develop a transient flow analysis taking into account several scenarios.The main objective of the present work is to report on the enhancement of the 2D shallow water flow model RiverFlow2D-GPU by means of a dynamic internal boundary conditions that allow the simulation of the evolution of dam breach processes as well as the propagation of the subsequent overland flow waves.
The outline of the present paper is as follows: First, the governing 2D shallow water equations are recalled together with the formulation chosen to model the dam failure processes including erosion, collapse and overflow.Then, the finite volume scheme used to solve the partial difference system on unstructured triangular cells is outlined.The novel formulation of internal boundary conditions representing several dam breach processes and the strategy adopted to incorporate the presence of the dam and its failure through piping is described.Two test cases are presented to show the model performance.In the absence of good experimental data, synthetic test cases are used.First, a case to simulate a collapse that leads to the transition from piping to dambreach.Then, a representative range of data collected by the USDA is also reproduced to show the efficiency of the model complex domains, leading to >76× faster simulations compared with the traditional CPU computation.

Governing Equations
The surface flow in the present work is described by means of the depth averaged 2D shallow water system of equations.They can be compactly formulated in terms of the conserved variables [14]: ∂U ∂t where h is the water depth and q x = hu and q y = hv the unit discharges in x and y direction, respectively.The fluxes of mass and momentum are: The bed slopes and friction terms are included on the source term of the system: where bed slopes in both plane directions are: and friction is formulated using the semi-empirical Manning law: where n is Manning's roughness coefficient [15].

Dambreach by Piping-Erosion Model
The novelty presented in this work is the formulation of internal boundary conditions representing several dam breach processes in combination with the 2D shallow water model.The 2D surface flow model is combined with a piping/dam breach evolution model presented in Chen et al. [16].The cross-section of the pipe is assumed to be initially composed of a semicircular arch of radius b/2 at the top with a square of side b at the bottom.Orifice flow and open channel flow equations are chosen to model the discharge for pressurized and free surface flows through the pipe, respectively.

Flow Discharge through the Pipe
If the pipe is fully filled with water, the discharge through can be estimated using the following orifice equation: where: 2 is the cross-sectional area of the pipe; b is the width of the base of the pipe; is the elevation of the pipe centre line; z b is the elevation of the pipe bottom; z s = h + z bed is the elevation of the water surface level; L is the pipe length; R = A/P is the pipe hydraulic radius; P = (3 + 0.5π)b is the pipe wetted perimeter; f = 8gn 2 R −1/3 is the Darcy-Weisbach friction factor of the pipe surface; n = 1 A n d 1/6 50 is Manning's roughness coefficient of the pipe surface; A n = 12 is a constant (Wu [17]).
On the other hand, in the case of a partially filled pipe, roof collapse or overtopping, the discharge is computed by means of a free surface flow equation: where: k sm is a dimensionless submergence correction for tailwater effects ( [18]); β is the breach side slope angle with respect to the horizontal; c r = 1.7 √ m/s is a discharge coefficient for the rectangular part of the breach section [19]; c t = 1.2 √ m/s is a discharge coefficient for the triangular part of the breach section [19].

Pipe Erosion
Regarding the erosion process, a shear stress-based formula is adopted to calculate the erosion rate.As erosion takes place, the full pipe cross section is enlarged along its height and width due to the removal of materials until the collapse of the pipe roof occurs (see Figure 1).The vertical erosion of the pipe is computed using an excess detachment rate relation: where: k d is the measured erosion coefficient at the breach; τ c is the critical stress required to initiate detachment for the material; is the bed shear stress in the pipe surface; ρ w is the water density.
The horizontal erosion rate is assumed to be equal to the vertical erosion rate and, hence, the evolution can be expressed as:

Pipe Roof Collapse
The collapse of the pipe roof is estimated by comparing the weight of the overlying soil and the cohesion of the soil on the two sidewalls of the pipe.The failure planes are assumed to be vertical and, for the sake of simplicity, the collapse is assumed to move downstream instantaneously.The arch finally becomes unstable due to the erosion of the pipe, and the collapse of the soil mass above the arch occurs.The failure of the roof occurs if the top of the eroding pipe (z b + 1.5b) reaches the top of the dam (z crest ).On the other hand, the failure also occurs if the driving force F d exceeds the soil resistant force F r (Figure 2).Hence, the model compares these two forces along the vertical direction.Once the driving force (equal to the weight of the failure part) is larger than the resistant force, the roof above the pipe will collapse: where p is the porosity of the soil material, G s = ρ d /ρ w is the specific gravity of the soil, ρ d is the soil material density and C is the cohesion of the dam fillings.The areas A a , A b and A c are computed as follows: where: L 1 is the dam crest width; α uw is the upward dam slope angle; α dw is the downward dam slope angle.The collapsed pipe roof is assumed to move downstream instantaneously.After the collapse, overtopping failure dominates and the breach flow discharge and the vertical erosion can be estimated using Equations ( 8) and (9), respectively.The relationship between horizontal expansion and vertical undercutting is given by the change in the breach top width ∆B: ∆B = 2∆z b sin β (13) and the change in the breach bottom width ∆b is given by:

Numerical Scheme
The system of Equations ( 2) is mandatory to be solved numerically due to the lack of analytical solution.For that purpose, an explicit upwind finite volume scheme, based on the Riemann Solver of Roe [20], is used in this work.This numerical scheme starts from system (2) which, in compact form, can be expressed as: where E = F, G .Integrating Equation ( 15) in a control volume or cell, Ω, and applying the divergence theorem to the second term: where ∂Ω is the contour of the control volume and n is the outgoing unit vector normal to the Ω volume.By discretizing Equation ( 16) in time and space, the basis of the numerical method in finite volumes is given by: where Ω i is the area of cell i, n is the current time level and the number of neighbouring cells is 3, because triangular cells are used in this work (as in the example in Figure 3).In addition, the fluxes E are evaluated at cell boundaries: where E j is the flux value at cell Ω j , which shares a wall k, of length l k , with cell Ω i with flux value E i .Considering the hyperbolic character of System ( 15), the Jacobian matrix normal to the flow direction, E, can be defined as: The local value of the Jacobian matrix (19) can also be defined at the walls k, being Jn k : where Λk is the diagonal matrix whose non-zero elements are the eigenvalues of the system, λm , and where Pm is the matrix containing the eigenvectors of the system, ẽm , providing this matrix with 3 eigenvalues to the 2D model.Therefore, starting from Expression (17), and using the eigenvalues and eigenvectors of the Jacobian matrix, the updating expression at cell i is obtained: The time step ∆t is calculated dynamically throughout the simulation, following: with 0 < CFL ≤ 1, where CFL is the Courant-Friedrichs-Lewy number, which guarantees stability in the numerical scheme, and where: The time step depends, on the one hand, on the dynamics of the problem to be solved through λm k and, on the other hand, on the cell size chosen for the computational grid, given by l k .Therefore, from Equations ( 22) and (23), it becomes clear that increasing the refinement of the computational leads to smaller time steps, and therefore to a higher computational cost.

Dambreach Flow as Internal Boundary Condition
Expression (21) allows the variables updating scheme at each interior cell of the domain using information that comes from the neighbour cells through the cell edges.When using triangular meshes, there are three contributions to update the three conserved variables.At boundary cells, at least one of the edges does not have a neighbour cell and boundary conditions are necessary to complete the information supplied by the numerical scheme.
The new contribution is that the presence of an internal hydraulic structure can be modelled by means of a mathematical condition defined along an interior line in the computational domain [21,22].This is called the internal boundary condition (IBC), i.e., each pair of cells sharing an edge on that internal line are considered internal boundary cells.These cells are updated using both information from the numerical scheme and from the IBC. Figure 4 shows an example of IBC defined along an internal boundary line, where several pairs of cells (filled with blue) on both sides of the line share an edge.An external law is used to define the module of the discharge through them while the water depth is provided by the numerical scheme updating Expression (21).Assuming the pipe roof collapse has not occurred and the flow direction from left to right, the upstream element is cell L and the downstream is cell R. Water surface elevation levels at the left side z n+1 s,L = (h + z) n+1 L , provided by the numerical scheme (21), are used to evaluate the discharge by means of the external discharge expression.First, the change in the pipe bottom elevation and width due to erosion are computed using ( 9) and ( 10) as: where (τ e − τ c ) n b accounts for the erosive shear stress at the breach evaluated at the time level t n .Then, the driving and the soil resistant forces, F d and F r , respectively, are computed using Equation (11), and the pipe roof collapse condition is checked.Therefore, two cases must be taken into consideration: 1.
If the roof collapse does not occur, the pipe is considered filled with water and pressurized so that the enforced cell discharge in pipe during the next time level is computed using (7) as: where the integrated breach features A, f , L, R are evaluated at time t n+1 ; 2.
If the roof collapse condition is satisfied, the dambreach is assumed open and the enforced cell discharge in pipe during the next time level is computed using (8) as: Note that once the roof collapse occurs, the pipe stability condition does not need to be checked any more, and the dam breach is assumed open.Regardless of whether the pipe is maintained or collapsed, the unit discharge at each cell pair composing the dambreach is assumed normal to the direction of the shared edge nb , and its module is updated as: where W b is the total breach width and l is the length of the shared edge for each internal boundary cell pair.

Test Case 1: Synthetic Dambreach through Piping Process
The purpose of this first test case is to validate the performance of the model at the time of the collapse that leads to the transition from piping to dambreach.The case consists of a reservoir of 7 million cubic metres of water contained by a 25.1 m high dam.At a given instant (t = 3 s), the dam will start to crack due to the piping mechanism.Figure 5 Left shows the topography of both the bottom of the reservoir and the domain through which the water will flow.Figure 5 Right shows the 4000 cells unstructured computational mesh used to discretize the domain.A more refined mesh has not been considered due to the purpose of the case does not require it.Figure 6 shows the distributed numerical results for the water depth at the initial state, t = 0.32 h, t = 0.36 h, t = 0.60 h, t = 1 h and t = 2 h, respectively.The results show a correct and stable temporal update of the water depth both inside and outside the reservoir through the flow through the breach estimated by means of the Equations ( 7) and (8).
On the other hand, Figures 7-9 show the temporal evolution of the breach water discharge Q, top (B) and bottom (b) pipe/breach width and top (z top ) and bottom level (z b ), respectively.In this particular case, the collapse occurs moments after the beginning of the piping process, at t = 0.306 h approximately, due to the upper limit of the breach has reached the crest of the dam (Figure 9). Figure 7 (lower) shows a detail of the abrupt change in breach flow around the collapse time.At t = 1 h approximately, the breach expansion stops due to the limits imposed by the case topography (Figure 8).After that, the breach discharge present minor fluctuations due to the 2D flow present within the reservoir, which causes variations in the upstream water level and, thence, in the breach discharge.As the level of the reservoir decreases, the water flow through the breach also does leading to an hydrograph consistent with those reported in the numerical tests and in the experiments of Chen et al. [16].

Test Case 2: Synthetic Set of Dams
The following test case consists of a collection of synthetic dams with different characteristics and was originally presented in [23].The main goal of this set of cases was to simulate a representative range of United States Department of Agriculture (USDA) small watershed structures by varying the dimensions (height and volume) of the reservoir and dam material composition.In [23], the simulations were performed using the software WinDAM C, developed by USDA Natural Resources Conservation Service (USDA-NRCS) in cooperation with USDA Agricultural Research Service (USDA-ARS) and Kansas State University [24].
The shape of the reservoirs can be described by a hypsometric function [23] as follows: where h represents the reservoir depth, relative to the base of the dam, h d is the maximum depth, V(h) is the volume of water within the reservoir at depth h, and V max the maximum volume.This reservoir shape can be visualized as a cone with vertex at the bottom of the reservoir and the base defined by the water surface (see Figure 10).Regarding the dam geometry, the set of cases includes six different heights and widths for the dam crest.The dimensions are included in Table 1 both using the unit system found in [23] and using the IS of units as used in our calculations.As stated in [23], a single length value of 274 m was set for all the cases as the average of the median lengths of the real world cases.Finally, a ratio of 3:1 was selected for upstream and downstream slopes, according to the US Bureau [25].Three different characterizations of the dam material have been considered (see Table 2), with the erodibility k d and the critical shear stress τ c being the most notable parameters for WinDAM C software [23].Equation ( 7) also requires the median particle size D 50 in mm for predicting the breach water discharge.Regarding the initial conditions, the internal erosion was assumed to be initiated with a void square in cross section and having dimensions of 6.1 cm high and 6.1 cm wide, located at the centreline of the dam and at elevation coinciding with 0.25h d above of dam base.
Figure 11 shows the details of the computational mesh and dam breach line for this test case.The mesh consists of 23,355 triangular cells, which ensures a good representation of the geometry of the dam. Figure 12 shows the numerical results in logarithmic scale obtained for the breach peak discharge in the 18 simulated cases (six geometries × three dam materials).They are labelled as High, Medium and Low, corresponding to the previously defined material characteristics.The results provided by WinDAM C are also shown, extracted from [23].In general, the same linear trends are observed for both models although the two-dimensional model predicts slightly higher values for the breach discharge than WinDAM C, especially for intermediate and high dam height values.The maximum discrepancy between both models appears in the 8 ft tall dam with the medium parameter set, with a relative difference of 75.9%, and for the minimum, also in the 8 ft tall dam, but choosing the high parameter set with a relative difference of 18.2%.Regarding the computational efficiency, the computational costs of the implementation of the model in GPU (two different GPU cards were tested) have been compared with the CPU model for this particular case, using the Low parameter set.The reason for choosing this set of parameters to perform this comparison is the slower evolution of the breach and, therefore, the longer simulation time required to reach the peak breach discharge.Table 3 shows the calculation times for both models, together with the speed-up factor, computed as the gain ratio with respect to the one-core CPU computational time.

Conclusions
In this work, the improvement of RiverFlow2D-GPU, a surface flow model based on 2D Shallow Water equations with a dambreach piping utility, has been presented.The overflow dambreach situation is also considered in the model.The dam is defined as an internal boundary condition where the operation conditions take into account all the possible hydraulic scenarios.The model is able to develop a transient water flow all over the computational domain, including the reservoir and the downstream valley at the same time.Every time step, the discharge through the dam is updated according to the flow conditions in order to obtain a perfect tracking of the dambreach evolution.Two applications at different spatial scales have been shown for validation of the model and to test the efficiency of the GPU implementation.The first case was used to validate the performance of the model at the time of the collapse that leads to the transition from piping to dambreach.The results are consistent with those reported in the numerical tests and in the experiments.Additionally, both CPU and GPU models provide exactly the same numerical results.A very good agreement was shown as well as good GPU efficiency.The second case consists of a collection of synthetic dams with different characteristics.The main goal of this set of cases was to simulate a representative range of United States Department of Agriculture (USDA) small watershed structures by varying the dimensions (height and volume) of the reservoir and dam material composition.The results provided by WinDAM follow the same linear trends observed in our model in general.The main conclusion of the present work is the successful enhancement of the 2D shallow water flow model RiverFlow2D-GPU by means of dynamic internal boundary conditions that allow the simulation of the evolution of dam breach processes, as well as the propagation of the subsequent overland flow waves.Additionally, the results obtained in this work support the generalized paradigm shift that is taking place within the scope of hydraulic simulation.The use of GPUs against CPUs, especially in computational meshes with a large number of cells, provides significant gains in simulation efficiency.In particular, to the authors' knowledge, no previous works consider either GPU computing for a 2D model including a dynamical internal boundary condition of the dam breach type.

Figure 1 .
Figure 1.Cross-section of the expansion due to piping process before the dam collapse (left) and trapezoidal breach evolution after the dam collapse (right).

Figure 2 .
Figure 2. Schematic diagram of the piping situation.

Figure 3 .
Figure 3. Diagram of the cells in a two-dimensional case with triangular cells.

Figure 5 .
Figure 5. Bed level (a) and detail of the dam geometry together with the triangular computational mesh (b).

Figure 6 .Figure 7 .
Figure 6.Temporal evolution of the surface water depth h after the dam breach at the initial state, t = 0.32 h, t = 0.36 h, t = 0.60 h, t = 1 h and t = 2 h, respectively.

Figure 8 .
Figure 8. Temporal evolution of the top (B) and bottom (b) width of the breach for a piping process starting at time t initial = 0.3 h.Full simulation time (a) and detail at the collapse time (b).

Figure 9 .
Figure 9. Temporal evolution of the top (z top ) and bottom level (z b ) of the breach for a piping process starting at time time t initial = 0.3 h.

Figure 10 .
Figure 10.Three-dimensional view of the Test Case 2 geometry.

Figure 11 .
Figure 11.Details of the computational mesh and dam breach line for Test Case 2.

Figure 12 .
Figure 12.Peak breach discharge as a function of dam height for synthetic sets (continuous line) and comparison with the results provided by WinDam C (dashed line).

Table 1 .
Set of dam heights and widths for Test Case 2.

Table 2 .
Set of dam material parameters for Test Case 2.