Parameter model of peak outflow rate.

## Abstract

Simulation of dam breach process has significant influence on the evaluation of consequence of dam breach flood. In this study, research progresses on the numerical modeling of earth-rock dams’ breach process are summarized, especially the latest research results of the author’s research team in recent years. However, there still has a considerable gap in the versatility of computer software and visualization technology of dam breaching process. It is suggested that more efforts should be made in the future to study the detailed physically based numerical model for core dam and concrete face rockfill dam; further, more attention should be paid to the application of visualization technology in dam breach process simulation. Finally, the universal and friendly visualization computer software that can accurately simulate the dam failure process and flood routing for earth-rock dams is sorely needed.

### Keywords

- earth-rock dam
- numerical model
- computer software
- research progress

## 1. Introduction

China has nearly 100,000 reservoir dams, of which earth-rock dams account for more than 95% [1, 2]. Most of these reservoir dams were built in the 1950s and 1970s. Due to economic and technical conditions at that time, the problem of dangerous reservoirs in China was outstanding [3]. According to statistics [4], from 1954 to 2018, 3541 reservoir dams broke in China. The “75·8” flood occurred in Henan in 1975, which led to the collapse of 2 large reservoirs in Banqiao and Shimantan (Figures 1 and 2), 2 medium-sized reservoirs in Tiangang and Zhugou, and 58 small reservoirs, causing heavy casualties and property losses [5]. In the twenty-first century, with the improvement of the dam safety management level and the comprehensive development of the reservoir’s risk elimination and reinforcement, the number of dam breaks has been significantly reduced, but due to the frequent occurrence of extreme weather events, dams’ breaching still occur frequently. On July 19, 2018, the Zenglongchang Reservoir in Inner Mongolia and the Sheyuegou Reservoir in Xinjiang on August 1, 2018, successively dams’ breaching [4] (Figures 3 and 4).

Therefore, it is necessary to establish a mathematical model and numerical calculation method that reasonably simulates the process of overtopping and seepage failure collapse, improves the prediction accuracy of the flood flow process of earth-rock dam collapse, and provides theoretical and technical support for the evaluation of the consequences of dam collapse and the preparation of emergency plans. This article will briefly introduce the research progresses on the mechanisms and numerical models of earth-rock dams’ breaching, especially the latest research results of the author’s research team in recent years, and make suggestions for future research.

## 2. Study on mathematical model of earth-rock dam break

The mathematical model of earth-rock dams’ breaching is generally divided into three categories [6]: The first category is the parameter model. Most of these models are based on statistical analysis of dam-break case data, and empirical formulas are used to calculate and obtain dam-break-related parameters. Although most models cannot consider the erosion characteristics of damming materials, but the parameter model formula is simple and fast to calculate and is also often used for rapid evaluation of the consequences of dams’ breaching. The second category is a simplified mathematical model based on the mechanism of failure. It is generally assumed that the shape of the fractured breach (rectangular, inverted trapezoidal, triangular, etc.) remains unchanged during the dams’ breaching. The method based on the flow shear stress and the critical shear stress of the dam material or the erosion formula of the dam material is used to calculate the breach development process; the weir flow is used (overtopping dam failure) or pore flow (seepage failure dam breaching) formulas are used to calculate the breach flow. The stability analysis of the breach slope mostly uses the limit equilibrium method; generally, the numerical calculation method based on time step iteration is used to simulate the breach development process and the breach flow process. The advantage of this type of model is that it considers the failure mechanism of earth-rock dams, and the calculation speed is relatively fast, which is the most widely used in the numerical simulation of earth-rock dams dam breaching process. The third category is a detailed mathematical model based on the failure mechanism. In recent years, a series of researches on one-dimensional, average two-dimensional, and three-dimensional mathematical models based on the hydrodynamic dam material erosion equation have made significant progress, which can simulate the dams’ breaching process of earth-rock dams in more detail. In order to deal with the diffuse overtopping flow composed of discontinuous mixed flow states, shock wave capturing methods such as approximate Riemann solution method and total variation declining (TVD) method are generally used, and finite volume method, level set method, and smooth particle hydrodynamic method are used to solve the governing equation. This type of model is a fast-developing simulation method in recent years, but it can only be used for the simulation of the overtopping collapse process of homogeneous dams or landslide dams. It has not been used to simulate the process of seepage and failure of earth-rock dams and the simulation of the process of overtopping failure of other types of earth-rock dams [6].

### 2.1 Parametric model

In 1977, Kirkpatrick [7] proposed the first empirical formula for predicting peak outflow *Qp *, and then scholars from various countries proposed a series of models. With the continuous enrichment of dam failure case investigation data and the deepening of research, the dam failure parameter model has gradually evolved from the single-parameter model to a multi-parameter model, and the output results have increased from the original peak outflow of the breach to the final average width of the breach and the duration of the dam and can consider the shape of the dam body, reservoir capacity, dam material characteristics, etc. The peak outflow rate of breach is very important for the evaluation of the consequences of dam breaching. Therefore, domestic and foreign scholars have studied more. The commonly used parameter model of peak outflow rate is shown in Table 1.

Model | Case number | Expression |
---|---|---|

Kirkpatrick (1977) [7] | 19 | Qp = 1.268(hw + 0.3)^{2.5} |

Soil Conservation Service (1981) [8] | 13 | Qp = 16.6hw ^{1.85} |

Hagen (1982) [9] | 6 | Qp = 0.54(hdS)^{0.5} |

Singh and Snorrason (1984) [10] | 28 | Qp = 13.4hd ^{1.89} or Qp = 1.776S ^{0.47} |

MacDonald and Langridge-Monopolis (1984) [11] | 23 | Qp = 1.154(Vwhw )^{0.412} |

Costa (1985) [12] | 31 | Qp = 0.981(hdS)^{0.42} |

Evans (1986) [13] | 29 | Qp = 0.72Vw ^{0.53} |

USBR (1988) [14] | 21 | Qp = 19.1hw ^{1.85} |

Froehlich (1995) [15] | 22 | Qp = 0.607Vw ^{0.295} hw ^{1.24} |

Walder and O’Connor (1997) [16] | 18 | Qp = 0.031 g ^{0.5} Vw ^{0.47} hw ^{0.15} hb ^{0.94} |

Xu and Zhang 1 (2009) [17] | 75 | Qp = 0.175 g ^{0.5} Vw ^{5/6}(hd /hr )^{0.199}(Vw ^{1/3}/hw )^{−1.274} eB _{4} |

Pierce et al. (2010) [18] | 87 | Qp = 0.0176(Vh)^{0.606}or Qp = 0.038 V ^{0.475} h ^{1.09} |

Thornton et al. (2011) [19] | 38 | Qp = 0.1202 L ^{1.7856} or Qp = 0.863 V ^{0.335} hd ^{1.833} Wave ^{−0.663or} Qp = 0.012 V ^{0.493} hd ^{1.205} L ^{0.226} |

Lorenzo and Macchione (2014) [20] | 14 | Qp = 0.321 g ^{0.258}(0.07Vw )^{0.485} hb ^{0.802}(overtopping)Qp = 0.347 g ^{0.263}(0.07Vw )^{0.474} hb ^{−2.151} hw ^{2.992}(seepage failure) |

Hooshyaripor et al. (2014) [21] | 93 | Qp = 0.0212 V ^{0.5429} h ^{0.8713} or Qp = 0.0454 V ^{0.448} h ^{1.156} |

Azimi et al. (2015) [22] | 70 | Qp = 0.0166(gV)^{0.5} h |

Froehlich 2 (2016) [23] | 41 | Qp = 0.0175kMkH (gVwhwhb ^{2}/Wave )^{0.5} |

Mei Shiang et al. (2018) [24] | 154 | Qp = Vwg ^{0.5} hw ^{−0.5}(Vw ^{1/3}/hw )^{−1.58}(hw /hb )^{−0.76}(hd /h _{0})^{0.10}e^{−4.55}(homogeneous dam)Qp = Vwg ^{0.5} hw ^{−0.5}(Vw ^{1/3}/hw )^{−1.51}(hw /hb )^{−1.09}(hd /h _{0})^{−0.12}e^{−3.61} (core-wall dam) |

In 1988, the US Bureau of Reclamation (USBR) [14] proposed the first empirical formula for predicting the final average width of the breach *Bave *, and then scholars from various countries put forward a series of models. The commonly used parameter model of the final average width of the breach is shown in Table 2.

Model | Case number | Expression |
---|---|---|

USBR (1988) [14] | 21 | Bave = 3hw |

Von Thun and Gillette 1 (1990) [25] | 57 | Bave = 2.5 hw + Cb |

Froehlich 2 (1995) [26] | 22 | Bave = 0.1803 K _{0}(Vw )^{0.32}(hb )^{0.19} |

Xu and Zhang 3 (2009) [17] | 75 | Bave = 0.787(hb )(hd /hr )^{0.133}(Vw ^{1/3}/hw )^{0.652} eB _{3} |

Froehlich 4 (2016) [23] | 41 | Bave = 0.27kM (Vw )^{1/3} |

Mei Shiang et al. (2018) [24] | 154 | Bave = hb (Vw ^{1/3}/hw )^{0.84}(hw /hb )^{2.30}(hd /h _{0})^{0.06}e^{−0.90}(homogeneous dam)Bave = hb (Vw ^{1/3}/hw )^{0.55}(hw /hb )^{1.97}(hd /h _{0})^{−0.07}e^{−0.09}(core-wall dam) |

In 1984, MacDonald and Langridge-Monopolis [42] proposed the first empirical formula for predicting the duration of dam failure, and then scholars from various countries proposed a series of models. Commonly used dam-break duration parameter model is shown in Table 3.

Model | Case number | Expression |
---|---|---|

MacDonald and Langridge-Monopolis (1984) [11] | 23 | Tf = 0.0179(0.0261(Vwhw )^{0.769})^{0.364} |

USBR (1988) [14] | 21 | Tf = 0.011Bave |

Froehlich (1995) [26] | 22 | Tf = 0.00254(Vw )^{0.53}(hb )^{−0.9} |

Xu and Zhang 1 (2009) [17] | 75 | Tf = 0.304Tr (hd /hr )^{0.707}(Vw ^{1/3}/hw )^{1.228} eB _{5} |

Froehlich (2016) [23] | 41 | Tf = 63.2(Vw /(ghb ^{2}))^{0.5} |

Mei Shiang et al. 2 (2018) [24] | 154 | Tf = T _{0}(Vw ^{1/3}/hw )^{0.56}(hw /hb )^{−0.85}(hd /h _{0})^{−0.32}e^{−0.20}(homogeneous dam)Tf = T _{0}(Vw ^{1/3}/hw )^{1.52}(hw /hb )^{−11.36}(hd /h _{0})^{−0.43}e^{−1.57}(core-wall dam) |

Due to the difficulties in obtaining the dam-break duration, the relatively low accuracy of the data, and the small number of samples, the dam-break duration model has a large deviation in the calculation of individual cases.

In order to fully consider the dam type, dam breach mode, reservoir characteristics, and breach characteristics, the reservoir capacity (*Vw *) is above the bottom of the breach at the dam break, the water depth (*hw *) above the bottom of the dam at the dam break (*hd *), and the final depth of the rupture (*hb *). For other parameters, the method of statistical regression is used to obtain the results of the peak flow of the breach, the final average width of the breach, and the duration of the dam breach. From the above statistics, it can be seen that the parameter model can simulate the dam-break parameters simply and quickly, which is an efficient and rapid evaluation method, but the parameter model cannot provide the dam-break flood flow process line.

### 2.2 Simplified mathematical model based on failure mechanism

In the 1960s, European and American scholars began to study a simplified mathematical model based on the mechanism of collapse based on hydraulics and sediment transport formulas. This model is also the most widely used mathematical model of earth-rock dams’ breaching. In 1965, from the US Bureau of Reclamation, Cristofano [27] established the first mathematical model of homogeneous dam overtopping failure. Afterward, scholars from various countries proposed a series of mathematical models for simulating earth-rock dam collapse [6, 28]. The most widely used is the NWS BREACH model developed by Fread from the National Weather Service [29]. In recent years, the Nanjing Hydraulic Research Institute and China Institute of Water Resources and Hydropower Research have conducted systematic research work on the mathematical model of earth-rock dams’ breaching, establishing NHRI-DB series and DB-IWHR series dam-break mathematical models, respectively. The commonly used simplified mathematical model of earth-rock dam breaching is shown in Table 4.

Model | Shape of breach | The flow of the breach | Erosion formula | Mechanical analysis | Breach mode | Type of dam |
---|---|---|---|---|---|---|

Cristofano (1965) [27] | Trapezoid | Wide crest weir formula | Cristofano formula | Breach without lateral collapse | Overtopping | Homogeneous |

BRDAM (1981) [30] | Parabolic | Wide crest weir formula, vent flow formula | Schoklitsch formula | Breach without lateral collapse (overtopping), top collapse (seepage) | Overtopping or seepage | Homogeneous |

DAMBRK (1984) [31] | Trapezoid or rectangle | Wide crest weir formula | Even flush | Breach without lateral collapse | Overtopping | Homogeneous |

BEED (1985) [32] | Trapezoid | Wide crest weir formula | Einstein and Brown formula, Meyer-Peter-Mueller formula | Collapse laterally | Overtopping | Homogeneous |

NWS BREACH (1988) [29] | Trapezoid or rectangle | Wide crest weir formula, vent flow formula | Correction Meyer-Peter-Mueller formula | Collapse laterally (overtopping), top collapse (seepage) | Overtopping or seepage | Homogeneous, core wall |

HR BREACH (2002, 2009) [33, 34] | Effective stress method | 1D stable non-uniform weir flow formula | Sediment transport formula or erosion rate formula | Single (two) side erosion, collapse laterally, stability analysis of core wall | Overtopping or seepage | Homogeneous earth dam, core wall |

FIREBIRD (2006) [35] | Trapezoid | Unsteady Saint-Venant equation | Sediment transport formula or erosion rate formula | Collapse laterally | Overtopping | Homogeneous |

WinDAM/SIMBA (2005, 2006, 2010) [36, 37, 38] | Rectangle | Wide crest weir formula | Erosion rate formula | Breach without lateral collapse | Overtopping | Homogeneous |

DLBreach (2013) [40] | Trapezoid | Wide crest weir formula, vent flow formula | Sediment transport formula or erosion rate formula | Single (two) side erosion, collapse laterally, stability analysis of core wall, dam foundation erosion | Overtopping or seepage | Homogeneous, core wall |

Hong Kong University of Science and Technology model [39] | Trapezoid | Wide crest weir formula, vent flow formula | Erosion rate formula | Collapse laterally (overtopping), top collapse (seepage) | Overtopping or seepage | Homogeneous, landslide |

DB-IWHR series dam-break mathematical model of China Institute of Water Resources and Hydropower Research [41, 42, 43] | Trapezoid | Wide crest weir formula | Erosion rate formula | Collapse laterally | Overtopping | Homogeneous, core wall, landslide dam |

NHRI-DB series dam-break mathematical model of Nanjing Hydraulic Research Institute [44, 45, 46, 47, 48] | Trapezoid | Wide crest weir formula, vent flow formula | Sediment transport formula or erosion rate formula | Shearing or dumping of the core wall, panel break, collapse laterally (overtopping), top collapse (seepage) | Overtopping or seepage | Homogeneous, core wall, face dam, landslide |

It can be seen from the above analysis that this type of model is mainly aimed at the two failure modes of earth-rock dam overtopping and seepage failure. By assuming the shape of the breach, different flow calculation formulas and erosion formulas are used to simulate the scouring of the dam material, and different simulation methods are used to analyze the vertical undercut and lateral expansion of the breach. Most of the models use iterative numerical calculation methods based on time steps to simulate the process of dam break and can output the parameters of dam break (such as the flow of the breach, the size of the breach, the water level of the reservoir, etc.) at each time step.

For example, based on the overtopping breach mechanism of the clay-core wall dam, a mathematical model to simulate its breach process is proposed. The model is based on the shape of the dam body and the characteristics of the flood flow to determine the initial scoring position of the downstream slope during erosion. The flow formula of the wide crested weir is used to calculate the rupture flow. The mechanical equilibrium method is used to simulate the tipping and shear failure of the core wall; the model can also consider the erosion of the dam body on one side, the erosion on both sides, and the erosion of the dam foundation and the process of water and soil coupling during dam break.

Based on the mechanism revealed by the model test of the overtopping breach of the homogeneous cohesive earth dam, the author has established a mathematical model that can simulate its collapse process (Figure 5). The specific modules of the model are as follows.

This model is based on the shape of the dam body and the characteristics of the flow at the top of the crater to determine the formation position of the “dark ridge.” The traceable erosion formula that can consider the physical and mechanical characteristics of the dam material is used to simulate the movement of the “dark ridge.” The collapse of the dam body: choose a reasonable erosion formula of the dam material to simulate the development of the dam crest and the downstream slope failure, and use the limit equilibrium method to simulate the failure of the collapse slope. The model considers incomplete dam failure and erosion of the dam foundation, as well as erosion on one side and both sides of the dam body.

The flow chart of the model calculation process of the collapse process of the homogeneous earth dam is shown in Figure 6.

There are two major highlights of the NHRI-DB concrete-face dam-break mathematical model [47]: the adoption of total-load nonequilibrium transport equation (Eq. (1)) [72] to simulate the erosion process of sand gravels with a wide range of gradation and the establishment of an analogy to simulate the failure process of each concrete-face slab under various loads during the dam breaching process.

where *t* = time; *x* = longitudinal coordinate; A = cross-sectional flow area in the breach channel; *Ct * = actual total-load sediment concentration; *Ct * = sediment concentration at the equilibrium state; and *L* _{s} = adaptation length characterizing the adjustment of sediment from a nonequilibrium state to equilibrium state.

In the NHRI-DB core dam-break mathematical model [45], a hydraulic method was used to predict the initial scour position for high dam. A time averaged erosion equation was adopted to simulate the backward erosion of dam’s shoulder. The broad-crested weir equation (Eqs. (2) and (3)) [73, 74] was adopted to calculate the breach flow discharge. Furthermore, the sliding or overturning failure was adopted as the key mechanism for the core, which was judged based upon numerical analysis. The calculated results show that the proposed model gives reasonable peak outflow, final breach width, and failure time.

where *Bb *is the bottom width of the breach (m), *H* represents the difference “*zs − zb *” (m), in which *zb *is the elevation of the breach bottom (m), *m* is the side slope (horizontal/vertical) of the breach, *c1 *and *c2 *are the discharge coefficients with values of 1.7 m0.5/s and 1.3 m0.5/s [73], and k_{sm} is the submergence correction factor for tailwater effects on weir outflow.

where *zt *is the tailwater level (m).

The advantage of this type of model is that it can consider the failure mechanism of the earth-rock dam and can use a short calculation time to complete the simulation of the dam-break process; however, most models cannot really consider the water-soil coupling effect during the dam-break process.

### 2.3 Detailed mathematical model based on failure mechanism

In order to fully describe the water-soil coupling effect in the process of dams’ breaching, in recent years, with the improvement of computer performance and the development of sediment science and computational fluid dynamics, a series of nonequilibrium dam material transport theory has emerged based on shallow water hypothetical detailed mathematical model for dam failure [49]. The commonly used detailed mathematical model of earth-rock dams’ breaching is shown in Table 5.

Model | Determination method of breach shape | Flow of breach | Dam material erosion | Mechanical analysis | Calculation method |
---|---|---|---|---|---|

Wang and Bowles (2006) [50] | Scour without sediment motion | Shallow water equations | Erosion rate formula | Three-dimensional collapse laterally | Finite different method |

Faeh (2007) [51] | Two-dimensional Exner equations | Shallow water equations | Traction load and suspended load formula | Collapse laterally | Finite volume method |

Wu et al. (2007, 2012) [52, 53] | One- and two-dimensional nonequilibrium total sand transport equations | General shallow water equations | Total sand transport formula | Collapse laterally | Finite volume method |

Swartenbroekx et al. (2010) [54] | Two-dimensional Exner equations | Shallow water equations | Traction load formula | Collapse laterally | Finite volume method |

Li et al. (2011) [55] | Two-dimensional nonequilibrium sediment transport equations (suspended load) | Shallow water equations | Empirical formulas of sediment carrying rate and Sedimentation rate | Breach without lateral collapse | Finite volume method |

Cao et al. (2011) [56] | Two-dimensional nonequilibrium total sediment transport equations | General shallow water equations | Traction load formula | Collapse laterally | Finite volume method |

Rosatti and Begnudelli (2013) [57, 58] | Two-dimensional mass conservation and energy conservation equations (solid phases) | Shallow water equations (liquid) | Floe concentration formula | Breach without lateral collapse | Finite volume method |

Juez et al. (2013, 2014) [59, 60] | One- and two-dimensional Exner equations | Saint-Venant equations and shallow water equations | 10 different erosion formulas | Breach without lateral collapse | Finite volume method |

Swartenbroekx et al. (2013) [61] | Two-dimensional mass conservation and energy conservation equations (traction load) | Shallow water equations (clean water) | Erosion rate formula | Breach without lateral collapse | Finite volume method |

Guan et al. (2014) [62] | Two-dimensional nonequilibrium sediment transport equations (traction load) | Shallow water equations (pure water) | Traction load formula | Collapse laterally | Finite volume method |

Kesserwani et al. (2014) [63] | Two-dimensional nonequilibrium sediment transport equations (suspended load) | Shallow water equations | Empirical formulas of sediment carrying rate and sedimentation rate | Breach without lateral collapse | Intermittent Galerkin method |

Razavitoosi et al. (2014) [64] | N-S equations (solid phases, non-Newtonian fluid) | N-S equations (liquid, non-Newtonian fluid) | / | Breach without lateral collapse | Smoothed particle hydrodynamics method |

Marsooli and Wu (2015) [65] | Three-dimensional nonequilibrium sediment transport equations | N-S equations | Traction load and suspended load formula | Breach without lateral collapse | Finite volume method and volume of fluid |

Abderrezza et al. (2016) [66] | Two-dimensional Exner equations | Shallow water equations | Traction load formula | Collapse laterally | Finite volume method |

Cantero-Chinchilla et al. (2016) [67] | One-dimensional nonequilibrium sediment transport equations | Saint-Venant equations, vertical momentum equation | Traction load and suspended load formula | Breach without lateral collapse | Finite volume method |

Cristo et al. (2016, 2018) [68, 69] | Two-dimensional mass conservation and energy conservation equations (solid phase) | Shallow water equations (liquid) | Traction load formula | Bed collapse algorithm | Finite volume method |

YAN Zhikun et al. (2019) [70] | Two-dimensional nonequilibrium sediment transport equations | General shallow water equations | Total sand transport formula considering bed slope | Bed collapse algorithm | Finite volume method |

It can be seen from the above statistics that this type of model is mainly based on the continuity equations of water flow (Eq. (4)), momentum equations (Eq. (5)), and energy equations (Eq. (6)), coupled with the sediment movement equation, and the finite volume method and other numerical simulation methods are used to discretely solve the governing equations

In Yan Zhikun’s model [70], based on the continuity equations of water flow, momentum equations, and nonequilibrium sediment transport equations, a planar two-dimensional mathematical model of dam rupture along the depth average is proposed. The sand capacity and the collapse mechanism of the two-dimensional slope during the dam-break process. The fully coupled method is used to convert the hydrodynamic equation and the nonequilibrium sediment transport equation into a shallow water equation with source terms and is based on the finite volume method under a rectangular grid. Discrete processing, using conservative, non-negative water depth numerical reconstruction format to make the model have second-order accuracy in the space–time direction, using HLLC [71] approximate Riemann solver to calculate grid boundary flux, SGM (Surface Gradient Method) format to calculate water surface gradient source terms, semi-implicit format. For the bottom bed friction term, the explicit gradient calculation of the source term of the concentration gradient is used to numerically solve the control equation.

Such models can achieve detailed simulation of the dam-break process, but the calculation speed is slow, and it can only be used for the numerical simulation of the overtopping dams’ breaching. However, this method can fully consider the coupling effect of water-soil coupling in the process of dam failure and can simulate complex boundary conditions, which is the development direction of numerical simulation of earth-rock dams’ breaching.

## 3. Conclusion and suggestions

Earth-rock dams’ breaching mechanism and dam-break process simulation are the foundation of dam-break disaster assessment and emergency response. They involve fluid mechanics, sediment kinematics, soil mechanics, and other disciplines. They are complex water-soil coupling problems. After decades of research and exploration, various mathematical models of dams’ breaching have been developed and made a series of innovative achievements, which provide theoretical support for improving the accuracy of flood disaster prediction of earth-rock dams. It is suggested that in the future, research efforts should be intensified on the mathematical model of the detailed simulation of the earth-rock dam breaching process, focusing on the application of visualization technology in the simulation of the dam-break process and accelerating the development of a universal and friendly simulation of the earth-rock dams’ breaching and the visual calculation of the disaster-causing process software.