Dam Break Modeling in a Cascade of Small Earthen Dams: Case Study of the ˇCižina River in the Czech Republic

: Failures of small dams can pose a serious threat to people and property even if the size of the schemes is relatively low. In many cases, small dams are situated in a cascade along streams, meaning that the failure of the uppermost dam may cause the dams downstream to fail. In this paper, a cascade of three small reservoirs, Lichnov II (14.6 m high), Lichnov III (10 m high), and Pocheˇn (8.5 m high), is the subject of the dam break analyses carried out via various methods such as empirical formulae, analogy, and hydraulic modeling. The dam-break ﬂood routing was simulated using a shallow water ﬂow hydraulic model. The simulations conﬁrm that the attenuation e ﬀ ect of the peak discharge is governed by the ﬂood volume, slope, and morphology of the ﬂoodplain and increases with the distance from the breached dam following an approximately exponential trend. When estimating peak discharge, empirical formulae derived for a single dam break should be applied carefully as they may underestimate the peak outﬂow by up to 10% in the case of a dam cascade. The attenuation volume of small reservoirs is small when compared to the ﬂood volume, meaning that the attenuation of the peak discharge usually varies between 5–10%.


Introduction
Every dam represents a certain threat due to the possibility of its failure, which may lead to a dam-break flood in the area downstream. Statistics show that the most frequent causes of earth-fill dam failures are overtopping and internal erosion [1][2][3]. The safety of small dams, ponds, and dry reservoirs is still topical. In the Czech Republic, there are about 20,000 small dams. During extreme floods in 1997, 2002, and 2010, dozens of collapsed small dams were reported [4]. Therefore, dam break studies, including dam-break flood propagation modeling are also carried out for systems of small dams. An example of the application of dam break modeling to a cascade of small earthen dams is presented in this case study of theČižina River in the Czech Republic.
For Emergency Action Plans, it is necessary to specify the characteristics of the dam-break flood and the extent of the endangered area. The knowledge concerning such flood-prone areas aids in the development of warning systems and evacuation plans. Various hydraulic and sediment transport models can be used for this purpose [5][6][7][8].
Numerous authors have dealt with the issue of dam breach modeling. In the 1970s, Cristofano [1] collected data from historical failures of earth-fill dams occurring up to 1965, based on which he performed a simulation of the onset of erosion. In the early eighties, the BRDAM model [9] was used to simulate the erosion of an earth dam in the case of overtopping or internal erosion. The lessons learned during past extreme floods [2,4,34] show that on average, about ten small dams fail in the Czech Republic during extreme flood events and that such dams are more vulnerable to breaching than large dams, which is primarily due to their insufficient spillway capacity. As a rule, such schemes have a flood attenuation volume, which is relatively small when compared to the volume of the arriving extreme flood. Their fixed spillways usually have a limited capacity due to the low importance of the dams.
This study aims to fill a gap in the research concerning the breaching of small embankment dams in a cascade. It deals with a cascade of three relatively new small embankment dams in theČižina River catchment in the Moravian-Silesian Region of the Czech Republic. One of the aims was to demonstrate the effect of the distance between dams on dam-break flood attenuation. In our case study, the two upper reservoirs are quite close to one another, while the lower one is located more remotely downstream along theČižina River. A complete set of technical data about the dams was available. The lowest dam, Pocheň, has already experienced overtopping twice but it resisted in both cases and did not fail.
Based on the results of the analysis, typical dam-break flood routing characteristics are summarized for a cascade of small dams. The method proposed in the study may be used for the development of Emergency Action Plans containing dam break and flood routing analyses for the affected area below a dam. The results of such analyses contain flood arrival times and maps with flood zones where water velocities and depths are included.

Study Area
The case study covers an area of theČižina River catchment in the Moravian-Silesian Region of the Czech Republic ( Figure 1) that extends for about 103 km 2 .
Water 2020, 12, x FOR PEER REVIEW 3 of 21 such schemes have a flood attenuation volume, which is relatively small when compared to the volume of the arriving extreme flood. Their fixed spillways usually have a limited capacity due to the low importance of the dams. This study aims to fill a gap in the research concerning the breaching of small embankment dams in a cascade. It deals with a cascade of three relatively new small embankment dams in the Čižina River catchment in the Moravian-Silesian Region of the Czech Republic. One of the aims was to demonstrate the effect of the distance between dams on dam-break flood attenuation. In our case study, the two upper reservoirs are quite close to one another, while the lower one is located more remotely downstream along the Čižina River. A complete set of technical data about the dams was available. The lowest dam, Pocheň, has already experienced overtopping twice but it resisted in both cases and did not fail.
Based on the results of the analysis, typical dam-break flood routing characteristics are summarized for a cascade of small dams. The method proposed in the study may be used for the development of Emergency Action Plans containing dam break and flood routing analyses for the affected area below a dam. The results of such analyses contain flood arrival times and maps with flood zones where water velocities and depths are included.

Study Area
The case study covers an area of the Čižina River catchment in the Moravian-Silesian Region of the Czech Republic ( Figure 1) that extends for about 103 km 2 . In the catchment area, three small reservoirs named Lichnov II, Lichnov III, and Pocheň are arranged in a cascade ( Figure 2). The cascade was selected because the distance between the dams differs significantly, allowing the demonstration of the impact of distance on dam-break flood attenuation. The two upper reservoirs are quite close to one another, while the lower one (Pocheň) is  In the catchment area, three small reservoirs named Lichnov II, Lichnov III, and Pocheň are arranged in a cascade ( Figure 2). The cascade was selected because the distance between the dams differs significantly, allowing the demonstration of the impact of distance on dam-break flood attenuation. The two upper reservoirs are quite close to one another, while the lower one (Pocheň) is about 15 km downstream from the Lichnov III dam. A complete set of technical data about the dams was available. The lowest dam, Pocheň, has already experienced overtopping twice, but it resisted in both cases and did not fail.

Description of the Small Dams
The main characteristics of the studied dams are described in Table 1. The Lichnov II dam is located about 1 km above the village of Lichnov on Tetřeví Stream, which is a tributary of the Čižina River. The dam was built in 2015. The main purpose of the reservoir is to provide flood protection to the villages downstream. The dam is homogenous and is made of sandy to gravelly clay. The maximum height of the dam above the local terrain is 16 m, and the dam crest is 5 m wide. A longitudinal section along the dam crest is shown in Figure 3, while a typical cross-section is displayed in Figure 4. The dam crest elevation at the left abutment is 0.85 m lower than the dam crest due to the presence of an emergency spillway. This place where overtopping may start, and water flows rapidly along the downstream slope was considered to be the potential point of origin of a breach due to overtopping. At the toe of the right abutment, extensive leakage through disintegrated colluvial soil below the base of the dam was identified via technical surveillance during the first trial filling of the reservoir. It is considered that there is a threat of backward erosion piping and a possible breach at this location.

Description of the Small Dams
The main characteristics of the studied dams are described in Table 1. The Lichnov II dam is located about 1 km above the village of Lichnov on Tetřeví Stream, which is a tributary of theČižina River. The dam was built in 2015. The main purpose of the reservoir is to provide flood protection to the villages downstream. The dam is homogenous and is made of sandy to gravelly clay. The maximum height of the dam above the local terrain is 16 m, and the dam crest is 5 m wide. A longitudinal section along the dam crest is shown in Figure 3, while a typical cross-section is displayed in Figure 4. The dam crest elevation at the left abutment is 0.85 m lower than the dam crest due to the presence of an emergency spillway. This place where overtopping may start, and water flows rapidly along the downstream slope was considered to be the potential point of origin of a breach due to overtopping. At the toe of the right abutment, extensive leakage through disintegrated colluvial soil below the base of the dam was identified via technical surveillance during the first trial filling of the reservoir. It is considered that there is a threat of backward erosion piping and a possible breach at this location. than the dam crest due to the presence of an emergency spillway. This place where overtopping may start, and water flows rapidly along the downstream slope was considered to be the potential point of origin of a breach due to overtopping. At the toe of the right abutment, extensive leakage through disintegrated colluvial soil below the base of the dam was identified via technical surveillance during the first trial filling of the reservoir. It is considered that there is a threat of backward erosion piping and a possible breach at this location.   The length of Tetřeví Stream between the Lichnov II and Lichnov III reservoirs is relatively small, about 500 m. The valley is relatively steep and narrow and is covered by forest and meadows.
The small detention dam Lichnov III is located approximately 400 m upstream of the village of Lichnov. It was built in 2018. Its purpose is to protect the village and to transform incoming flood waves from the Lichnov II dam and from the area adjacent to the reservoir. The dam body is homogeneous. It is made of sandy to gravelly clay and is very similar to that of Lichnov II. The dam height is about 12 m; see Table 1 for its main characteristics and parameters. As a dam-break flood arriving from the Lichnov II dam would be sudden, it is not probable that a gradually evolving piping process would develop, and so only overtopping was considered as a potential failure mode. The place of overtopping was similar to that of Lichnov II, i.e., at the auxiliary crest spillway at the left abutment ( Figure 5). About 600 m downstream of the Lichnov III dam, Tetřeví Stream flows into the Čižina River, which passes through the village of Lichnov. Numerous civil structures such as bridges and culverts are located there, and the area along the 3.5 km long river reach is urbanized with numerous buildings. Further, on the downstream of the village, the Čižina River continues into a region of fields and meadows covered with trees. This section of the Čižina River stretches for about 6 km, after which the river enters the Pocheň reservoir.
The small dam named Pocheň is located on the Čižina River about 5 km upstream of its junction with the larger Opava River. The Pocheň dam was built between 1973 and 1975 as a source of water for agricultural activities in nearby areas. During the extreme flood event in the year 1996, the dam was seriously damaged by overtopping but not breached. Due to this damage, reconstruction of the dam body took place, and the spillway capacity was improved over the following years. The dam height is about 9.5 m, and the dam is constructed from local gravelly clay. The geodetic survey indicated the lowest crest elevation is approximately in the middle of the crest length (see Figure 6). Due to the short duration of the extreme load, a progressive piping failure would not develop, and The length of Tetřeví Stream between the Lichnov II and Lichnov III reservoirs is relatively small, about 500 m. The valley is relatively steep and narrow and is covered by forest and meadows.
The small detention dam Lichnov III is located approximately 400 m upstream of the village of Lichnov. It was built in 2018. Its purpose is to protect the village and to transform incoming flood waves from the Lichnov II dam and from the area adjacent to the reservoir. The dam body is homogeneous. It is made of sandy to gravelly clay and is very similar to that of Lichnov II. The dam height is about 12 m; see Table 1 for its main characteristics and parameters. As a dam-break flood arriving from the Lichnov II dam would be sudden, it is not probable that a gradually evolving piping process would develop, and so only overtopping was considered as a potential failure mode. The place of overtopping was similar to that of Lichnov II, i.e., at the auxiliary crest spillway at the left abutment ( Figure 5).  The length of Tetřeví Stream between the Lichnov II and Lichnov III reservoirs is relatively small, about 500 m. The valley is relatively steep and narrow and is covered by forest and meadows.
The small detention dam Lichnov III is located approximately 400 m upstream of the village of Lichnov. It was built in 2018. Its purpose is to protect the village and to transform incoming flood waves from the Lichnov II dam and from the area adjacent to the reservoir. The dam body is homogeneous. It is made of sandy to gravelly clay and is very similar to that of Lichnov II. The dam height is about 12 m; see Table 1 for its main characteristics and parameters. As a dam-break flood arriving from the Lichnov II dam would be sudden, it is not probable that a gradually evolving piping process would develop, and so only overtopping was considered as a potential failure mode. The place of overtopping was similar to that of Lichnov II, i.e., at the auxiliary crest spillway at the left abutment ( Figure 5). About 600 m downstream of the Lichnov III dam, Tetřeví Stream flows into the Čižina River, which passes through the village of Lichnov. Numerous civil structures such as bridges and culverts are located there, and the area along the 3.5 km long river reach is urbanized with numerous buildings. Further, on the downstream of the village, the Čižina River continues into a region of fields and meadows covered with trees. This section of the Čižina River stretches for about 6 km, after which the river enters the Pocheň reservoir.
The small dam named Pocheň is located on the Čižina River about 5 km upstream of its junction with the larger Opava River. The Pocheň dam was built between 1973 and 1975 as a source of water for agricultural activities in nearby areas. During the extreme flood event in the year 1996, the dam was seriously damaged by overtopping but not breached. Due to this damage, reconstruction of the dam body took place, and the spillway capacity was improved over the following years. The dam height is about 9.5 m, and the dam is constructed from local gravelly clay. The geodetic survey indicated the lowest crest elevation is approximately in the middle of the crest length (see Figure 6). Due to the short duration of the extreme load, a progressive piping failure would not develop, and so only the overtopping failure mode was taken into account. About 600 m downstream of the Lichnov III dam, Tetřeví Stream flows into theČižina River, which passes through the village of Lichnov. Numerous civil structures such as bridges and culverts are located there, and the area along the 3.5 km long river reach is urbanized with numerous buildings. Further, on the downstream of the village, theČižina River continues into a region of fields and meadows covered with trees. This section of theČižina River stretches for about 6 km, after which the river enters the Pocheň reservoir.
The small dam named Pocheň is located on theČižina River about 5 km upstream of its junction with the larger Opava River. The Pocheň dam was built between 1973 and 1975 as a source of water for agricultural activities in nearby areas. During the extreme flood event in the year 1996, the dam was seriously damaged by overtopping but not breached. Due to this damage, reconstruction of the dam body took place, and the spillway capacity was improved over the following years. The dam height is about 9.5 m, and the dam is constructed from local gravelly clay. The geodetic survey indicated the lowest crest elevation is approximately in the middle of the crest length (see Figure 6). Due to the short duration of the extreme load, a progressive piping failure would not develop, and so only the overtopping failure mode was taken into account. Downstream of the Pocheň reservoir, the Čižina River passes through an area of forest, flows below several road and railway bridges, and then enters the Opava River close to the village of Skrochovice.

Methods
Both empirical formulae and comparison with real incidents were applied for the first preliminary estimates of the dam break peak discharges. The obtained values provide an idea about the range of the peak outflows and provide a general framework for numerical modeling. Due to the lack of relevant incidents and reliable data from real breached dams, these equations suffer from a great deal of inaccuracy and are mainly useful for the preliminary estimation of parameters, which can then be determined more precisely by numerical models.
The dam break simulation was carried out using the numerical model of breaching due to overtopping and internal erosion (the latter phenomenon at the uppermost dam only), which is normally applied to individual dams. The peak discharges were compared with the results of empirical formulae and with real cases involving the failure of small dams.
Two-dimensional shallow-water flow hydraulic models were applied for the simulation of the dam-break flood propagation into the areas between the dams. The flow and erosion parameters were searched via Monte-Carlo random sampling to achieve the highest peak outflow discharge below the dam.

Empirical Formulae
In order to obtain an initial idea about the peak discharges, published empirical formulae were applied. This approach is the simplest but least accurate way of predicting dam breach parameters. In dam-break simulations, direct calibration of models is usually impossible. One of the methods of verifying numerical models is to compare the modeling results with the outputs of retrospective studies of real incidents and failures. Unlike, e.g., hydrological floods, there is a lack of dam-break calibration data worldwide due to a lack of relevant incidents and reliable data from real breached dams. The empirical formulae discussed in this study (Table 2) are based on such real cases. It is true that statistical analysis is difficult due to the heterogeneity, inaccuracy, and scatter of the dam break data when subjected to regression analysis. However, the published empirical formulae and the analogy method are frequently used in dam break analysis for the comparison of real incidents with the results of numerical analysis.
One of the aims of this study was to demonstrate the difference between a single dam break and a cascade when applying empirical formulae. Table 2. Summary of formulae used for estimating the peak breach outflow Qb. Downstream of the Pocheň reservoir, theČižina River passes through an area of forest, flows below several road and railway bridges, and then enters the Opava River close to the village of Skrochovice.

Methods
Both empirical formulae and comparison with real incidents were applied for the first preliminary estimates of the dam break peak discharges. The obtained values provide an idea about the range of the peak outflows and provide a general framework for numerical modeling. Due to the lack of relevant incidents and reliable data from real breached dams, these equations suffer from a great deal of inaccuracy and are mainly useful for the preliminary estimation of parameters, which can then be determined more precisely by numerical models.
The dam break simulation was carried out using the numerical model of breaching due to overtopping and internal erosion (the latter phenomenon at the uppermost dam only), which is normally applied to individual dams. The peak discharges were compared with the results of empirical formulae and with real cases involving the failure of small dams.
Two-dimensional shallow-water flow hydraulic models were applied for the simulation of the dam-break flood propagation into the areas between the dams. The flow and erosion parameters were searched via Monte-Carlo random sampling to achieve the highest peak outflow discharge below the dam.

Empirical Formulae
In order to obtain an initial idea about the peak discharges, published empirical formulae were applied. This approach is the simplest but least accurate way of predicting dam breach parameters. In dam-break simulations, direct calibration of models is usually impossible. One of the methods of verifying numerical models is to compare the modeling results with the outputs of retrospective studies of real incidents and failures. Unlike, e.g., hydrological floods, there is a lack of dam-break calibration data worldwide due to a lack of relevant incidents and reliable data from real breached dams. The empirical formulae discussed in this study (Table 2) are based on such real cases. It is true that statistical analysis is difficult due to the heterogeneity, inaccuracy, and scatter of the dam break data when subjected to regression analysis. However, the published empirical formulae and the analogy method are frequently used in dam break analysis for the comparison of real incidents with the results of numerical analysis.

Author, Reference
Formulae for Q b *) (m 3 /s) The following breach dimensions can be used: Webby (1996) [41] Q Vischer, Hager (1998) [43] *) An explanation of the variables can be seen in the section "List of symbols".
One of the aims of this study was to demonstrate the difference between a single dam break and a cascade when applying empirical formulae.

The Analogy Method
A more reliable peak breach outflow may be derived via the comparison of data from real incidents at dams with similar characteristics (dam type, height, reservoir volume) to the studied dam. The data may be obtained from publications and extensive dam failure databases (both national and international) [2,3,34,42,44,45]. In Table 3, dams similar to those studied are mentioned together with basic relevant characteristics (see also Chapter 2).

Embankment Dam Breach Modeling
A simple parametric model based on hydraulic and transport modules was used for the simulation. The tailor-made code was developed in MATLAB software to allow the dam breaching process to be modified easily according to local conditions and the failure mode.
The main output of the dam break analysis is the outflow hydrograph Q(t), which results from the failure. The following steps were performed before the modeling took place.
A careful investigation was performed with regard to the dam type and structure (material, emergency spillway position, dam surface, type of pavement on the dam crest, etc.), as well as the local geological conditions (subbase material, material discontinuities).
The probable place of failure initiation was located; for an overtopping failure, the place with the minimum dam crest or auxiliary spillway height was chosen. For a piping failure, the location with the weakest subbase providing the highest seepage was chosen.

1.
The shape of the dam break opening was rectangular for overtopping and circular for piping.

2.
The flood hydrograph Q in (t) entering the reservoir was provided by Czech Hydrometeorological Institute.
Water 2020, 12, 2309 8 of 21 3. The bathymetry of the reservoirs was taken from a geodetic survey of the reservoir area. 4.
The initial reservoir water level H 0 (t) was taken from the dam operation manual. It corresponds to the live storage water level.
The following relation between the reservoir inflow, outflow, and a change in water level holds: where V is the reservoir water volume (m 3 ), t is time (s), Q in is the reservoir inflow (m 3 /s), Q b is the breach outflow (m 3 /s) (in the following text, Q bo denotes overtopping discharge and Q bp piping discharge), Q f is the outflow through appurtenant works (outlets, spillways, etc.) in (m 3 /s), and H is the elevation of the water level (m a. s. l.). The outflow through the appurtenant works governed by the reservoir water level was calculated using equations taken from standard dam hydraulics [46]. The equations used and related coefficients were determined according to the arrangement of the hydraulic structures.

Overtopping
In the case of dam overtopping, the breach was approximated using a 1D model. The breach started at the lowest point of the dam crest, where the first overtopping occurred. The overtopping width along the dam crest was suggested as an initial breach width, b 0 . The shape of the breach opening was approximated by a rectangle, and the erosion progress was applied in both the downward and the lateral direction. The water moving along the downstream slope was assumed to exhibit quasi-steady uniform flow [3,13,45] with no submergence from the downstream side. The resistance of the dam body against the surface erosion was evaluated with respect to the water flow velocity at the downstream slope, with non-scouring velocity being used for the evaluation. Parallel uniform gradual backward erosion of the downstream slope was assumed [13]. The angle of the downstream slope and the width of the dam crest were set as input parameters, along with hydraulic parameters such as the discharge coefficient m, Manning roughness n, and the non-scouring velocity v ns and erosion parameters of the soil α 1 and α 2 , together with initial conditions (see below).
The breach outflow during overtopping Q bO was calculated using an overflow formula for a rectangular broad crested weir [47]: where m is the weir coefficient, b is the width of the breach opening (m), g is the acceleration due to gravity (m/s 2 ), and Z is the breach bottom elevation (m a. s. l.). Using the overtopping discharge obtained from Equation (2), the water depth h f and velocity v f along the slope was derived from the Chezy equation [47]: where n is the roughness of the downstream slope, and β is the downstream slope angle. At the moment, when the flow velocity along the slope exceeds the non-scouring velocity v ns (the value is an input parameter), erosion occurs and the breach starts to increase in width and depth [34]: (6) where α 1 and α 2 are erosion parameters. Initial conditions specify the reservoir water level H 0 , the lowest dam crest elevation Z 0 , and initial overtopping width b 0 at the beginning of the solution (t = t 0 ). The set of Equations (1) to (6) was solved numerically using the Newton method with a small-time step ∆t = 1 s.

Piping
For the piping model, a circular opening with a constant leakage pipe diameter and constant shear stress and hydraulic gradient along the pipe were assumed [19][20][21][22]28,48]. The input parameters were the length of the pipe L, the initial reservoir water level H 0 , the elevation of the pipe outflow H outflow , the critical shear stress along the walls of the pipe τ c , the soil erodibility C e , the initial pipe diameter D 0 , and absolute roughness ∆.
The outflow through the pipe Q bP and average water velocity v pip in the pipe were determined using pipe hydraulics [47]: where α is the kinetic energy (Coriolis) coefficient, L is the length of the pipe (m), D is the pipe diameter (m), λ is the friction loss coefficient, S is the pipe cross-sectional area (m 2 ), and H outflow is the pipe outflow elevation (m a. s. l.).
Using the water velocity in the pipe, the average hydraulic gradient along the pipe I E and shear stress τ acting on the pipe walls were obtained: where ρ is the water density. When the critical shear stress τ c on the pipe wall was exceeded, erosion and an increase in the pipe diameter initiated. The rate of erosion . ε and the change in pipe radius dr were determined according to equations suggested by Wan and Fell [48]: where .
ε is the rate of erosion (kg/s/m 2 ), C e is the coefficient of soil erosion (s/m), ρ b is the bulk density of the eroded soil (kg/m 3 ), and r = D/2 is the pipe radius.
The initial reservoir water level H 0 and initial pipe diameter D 0 are the initial conditions. As in the case of overtopping, the solution of Equations (1) and (7) to (12) was carried out using the Newton method with a very small-time step of ∆t = 1 s.
During the solution, the relation between the pipe diameter and the thickness of the roof of the conduit (overburden) was checked. After the pipe diameter exceeded the overburden thickness, it was assumed that the soil layer above the opening collapsed and the process continued similarly as in the case of overtopping. The width of the breach opening corresponded to the pipe diameter at the instant of the overburden collapse, while the breach bottom was assumed to be at the top of the seepage pipe at the instant of the collapse.

Optimization Procedure
Most of the input parameters in both scenarios are rather uncertain, and realistic intervals of their values were estimated from literature sources and real cases [3,41]. The combination of parameter values from their intervals may result in significant differences in outflow hydrographs and peak outflow discharge. In order to obtain the least favorable result with the highest peak outflow, the "optimization" of the parameters mentioned in Table 4 was performed via the Monte-Carlo sampling procedure [49]. The most influential parameters for the optimization were determined based on the equations governing the breach outflow discharge. The shape of the dam break hydrograph and namely the peak outflow discharge is generally governed dominantly by the size of the breach opening and the water level in the reservoir (Equations (2) and (7)). Crucial variables influencing the shape and size of the breach opening are the erodibility parameters (Equations (4)-(6) for overtopping and Equations (8) and (12) for piping). The reservoir water level is determined via level pool routing through the reservoir, Equation (1), and depends on the reservoir bathymetry V(H), the inflow flood hydrograph Q in , the capacity of the appurtenant works Q f , and the breach outflow Q b . The importance of the mentioned variables was proved by sensitivity and statistical analyses carried out in a previous study [50]. In this study, the members in Equation (1) were not subject to optimization as they were "prescribed" by the Czech Hydrometeorological Institute (inflow hydrograph), by the reservoir bathymetry, and by the capacity of the spillway.
For each dam, about 10 thousand simulations were carried out using parameter values obtained from random sampling. A uniform probability distribution was used for each parameter. As a result, a "sub-optimal" scenario was found after a reasonable amount of computational time with sufficient accuracy when compared to other simplifications adopted (the mathematical model, for example). The scenario with the highest peak outflow was used for further processing.

Dam-Break Flood Routing
The dam-break flood routing between reservoirs was simulated using a two-dimensional (2D) model based on shallow-water flow equations [8]. The solution was performed using the commercial code HEC-RAS v. 5.0.7 [5]. As the model used is standard and the procedures are routine, the model is not described here in detail, and the authors refer those interested to the available literature [5,46]).
The inputs for the calculation were the topology and geometry of the area between the dams taken from a digital model of the terrain (DTM). Civil structures were assumed to be fixed, with no damage suffered due to the flow of water.
The upstream boundary condition was represented by the dam breach flood hydrograph and hydrological flood hydrographs at contributing local streams. The initial condition was set as a stream discharge and water level corresponding to a 20-year return period, which corresponds approximately to the bank-full capacity of the streams.

Simulations and Results
First, the peak outflow from the dams was estimated using empirical formulae (Table 2) based on the dam and reservoir characteristics (Table 3). A comparison with real dam failures was also performed.
The numerical simulation started from the uppermost dam, Lichnov II. From an analysis of hydrological and empirical data, a check flood inflow hydrograph from the area upstream of Lichnov II was derived by the Czech Hydrometeorological Institute (Figure 7), providing a peak discharge of 50 m 3 /s. This approximately corresponds to the outflow of the largest experienced rainfalls at small river basins in the Czech Republic (more than 200 mm/day) during the floods in the years 1996, 1997, and 2010 at relatively small catchments with areas of 5 to 10 km 2 .  The runoff from the larger catchment area between the reservoirs at the time of the dam break events corresponds to a flood with a 20-year return period. This assumption is also based on past experience from real floods at catchments with areas of 50 to 60 km 2 . This also corresponds to the capacity of related streams in smaller urban areas.
The dam break study proceeded from the upstream reservoir in the downstream direction in five steps containing level pool routing: through the Lichnov II reservoir (simulation of the dam breaching due to both overtopping and piping), -between Lichnov II and Lichnov III, -through the Lichnov III reservoir (simulation of the dam breaching due to overtopping), -between Lichnov III and Pocheň, -through the Pocheň reservoir (simulation of the dam breaching due to overtopping).
In order to achieve the worst-case scenario of dam breaching, in which the maximum peak outflow occurs from the collapsed dams, parameter optimization was carried out. The optimization was conducted using the parameters shown in Table 4, with their realistic intervals adopted from the available literature [2,3,37,40,45,50].
A summary of the final input values, including the final parameters obtained from the optimization procedure, is listed in Table 5. The input values mentioned in Table 5 result in the maximum outflow discharge from the reservoirs during the breach event.  The runoff from the larger catchment area between the reservoirs at the time of the dam break events corresponds to a flood with a 20-year return period. This assumption is also based on past experience from real floods at catchments with areas of 50 to 60 km 2 . This also corresponds to the capacity of related streams in smaller urban areas.
The dam break study proceeded from the upstream reservoir in the downstream direction in five steps containing level pool routing: through the Lichnov II reservoir (simulation of the dam breaching due to both overtopping and piping), -between Lichnov II and Lichnov III, -through the Lichnov III reservoir (simulation of the dam breaching due to overtopping), -between Lichnov III and Pocheň, -through the Pocheň reservoir (simulation of the dam breaching due to overtopping).
In order to achieve the worst-case scenario of dam breaching, in which the maximum peak outflow occurs from the collapsed dams, parameter optimization was carried out. The optimization was conducted using the parameters shown in Table 4, with their realistic intervals adopted from the available literature [2,3,37,40,45,50].
A summary of the final input values, including the final parameters obtained from the optimization procedure, is listed in Table 5. The input values mentioned in Table 5 result in the maximum outflow discharge from the reservoirs during the breach event.  At the Lichnov II dam, both overtopping and piping failure were modeled. Figure 8 shows that the piping initiated approximately at the time of the peak inflow while overtopping occurred about 25 min after the inflow peak. The initialization phase at overtopping was followed by the gradual enlargement of the breach opening up to a final width of 68 m, which was reached about 125 min after the first overtopping as a result of intensive erosion only at the left side towards the central axis of the dam. A comparison of the modeled breach outflow hydrographs (Figure 8) shows that the peak discharge caused by the overtopping failure (802 m 3 /s) was more than 2.5 times higher than that caused by piping (316 m 3 /s). Therefore, the hydrograph from the overtopping failure was adopted for further analysis (Figure 8) of the cascade. This hydrograph provides a quite realistic duration for the breach event (approimately 30 min), which corresponded to that experienced during the breaching of small dams in the Czech Republic [34]. At the Lichnov II dam, both overtopping and piping failure were modeled. Figure 8 shows that the piping initiated approximately at the time of the peak inflow while overtopping occurred about 25 min after the inflow peak. The initialization phase at overtopping was followed by the gradual enlargement of the breach opening up to a final width of 68 m, which was reached about 125 min after the first overtopping as a result of intensive erosion only at the left side towards the central axis of the dam. A comparison of the modeled breach outflow hydrographs (Figure 8) shows that the peak discharge caused by the overtopping failure (802 m 3 /s) was more than 2.5 times higher than that caused by piping (316 m 3 /s). Therefore, the hydrograph from the overtopping failure was adopted for further analysis (Figure 8) of the cascade. This hydrograph provides a quite realistic duration for the breach event (approimately 30 min), which corresponded to that experienced during the breaching of small dams in the Czech Republic [34]. Most of the studied real dam failures exhibited peak discharges of 550 to 680 m 3 /s, except for in the case of the Buffalo Creek dam, which had 1416 m 3 /s ( Table 6). Empirical formulae provide values, which are closer to the higher modeled discharge (485 to 1763 m 3 /s). Most of the studied real dam failures exhibited peak discharges of 550 to 680 m 3 /s, except for in the case of the Buffalo Creek dam, which had 1416 m 3 /s (Table 6). Empirical formulae provide values, which are closer to the higher modeled discharge (485 to 1763 m 3 /s). The flood routing between the Lichnov II and Lichnov III reservoirs along the 500 m long valley reach provided relatively small attenuation of the peak discharge from 802 m 3 /s to approximately 756 m 3 /s with a shift in the peak of about 3 min (Figure 9). The flood routing between the Lichnov II and Lichnov III reservoirs along the 500 m long valley reach provided relatively small attenuation of the peak discharge from 802 m 3 /s to approximately 756 m 3 /s with a shift in the peak of about 3 min (Figure 9). As mentioned above, the Lichnov III dam would be breached primarily by overtopping. The flood arrival was quite rapid (10 min), so the breaching started shortly after the rise in the reservoir water level. About 1 h was needed for the initiation of piping (Figure 8), and, therefore, the failure mode was primarily overtopping.
The modeled peak outflow of 841 m 3 /s from the dam occurred a few minutes after the arrival of the flood from the upper Lichnov II reservoir. The duration of the breach event was about 25 min, and the time at which peak outflow occurred is slightly after the time that peak inflow to the reservoir took place (Figure 10). The maximum breach width was 25 m. From Table 6 it can be seen that the range of values for the peak outflow derived from individual formulae varies from 231 to 1225 m 3 /s; the modeled peak outflow exceeds practically all calculated values (except for those obtained from USACE methodology). This could be attributed to the domino cascade effect that occurs when two dams are in very close proximity. Comparisons with empirical relationships and real dam failures may not be entirely relevant, as the difference between single and cascade dam breach peak discharge As mentioned above, the Lichnov III dam would be breached primarily by overtopping. The flood arrival was quite rapid (10 min), so the breaching started shortly after the rise in the reservoir water level. About 1 h was needed for the initiation of piping (Figure 8), and, therefore, the failure mode was primarily overtopping.
The modeled peak outflow of 841 m 3 /s from the dam occurred a few minutes after the arrival of the flood from the upper Lichnov II reservoir. The duration of the breach event was about 25 min, and the time at which peak outflow occurred is slightly after the time that peak inflow to the reservoir took place (Figure 10). The maximum breach width was 25 m. From Table 6 it can be seen that the range of values for the peak outflow derived from individual formulae varies from 231 to 1225 m 3 /s; the modeled peak outflow exceeds practically all calculated values (except for those obtained from USACE methodology). This could be attributed to the domino cascade effect that occurs when two dams are in very close proximity. Comparisons with empirical relationships and real dam failures may not be entirely relevant, as the difference between single and cascade dam breach peak discharge is about 50-55%. It should be noted that empirical equations correspond to single dam failure and do not take into account the domino effect. Flood routing along the remaining reach of Tetřeví Stream, and further on along theČižina River resulted in a transformed hydrograph entering the Pocheň reservoir. From Figure 10, it could be seen that the peak discharge dropped from 841 m 3 /s to about 330 m 3 /s at the entrance to the Pocheň reservoir due to the attenuation along the 10 km long river valley between the Lichnov III and Pocheň dams ( Figure 2). The time shift in the flood arrival time between both reservoirs is about 1 h.
The input values of the erosion resistance characteristics for the simulation of the Pocheň dam breaching were verified using data from real overtopping events that took place during floods in 1977 and 1996 ( Figure 11). They correspond to a non-scouring velocity of v n ≈ 3.1 m/s, which was set in the calculation (it was not subject to optimization).
The resulting reservoir inflow and outflow hydrographs obtained by the simulation are shown in Figure 12. It can be seen that the peak outflow was delayed for approximately 1 h and 20 min after the culmination at the inflow. This time corresponds to the gradual filling of the flood attenuation reservoir storage, the overtopping of the relatively horizontal crest of the dam, and the slow development of the breach opening in the relatively resistant dam material. It can be seen that the peak inflow of 330 m 3 /s into the reservoir was increased by the breaching to a peak outflow of 557 m 3 /s, which is slightly higher than the average peak discharge obtained via comparative analysis and empirical formulae ( Table 6). The minor difference between the simulated and empirical values was due to the characteristics of the cascade and other influences, including significant flood attenuation between the dams, the improved spillway capacity, and the rather higher dam body resistance against erosion. The resulting reservoir inflow and outflow hydrographs obtained by the simulation are shown in Figure 12. It can be seen that the peak outflow was delayed for approximately 1 h and 20 min after the culmination at the inflow. This time corresponds to the gradual filling of the flood attenuation reservoir storage, the overtopping of the relatively horizontal crest of the dam, and the slow development of the breach opening in the relatively resistant dam material. It can be seen that the peak inflow of 330 m 3 /s into the reservoir was increased by the breaching to a peak outflow of 557 m 3 /s, which is slightly higher than the average peak discharge obtained via comparative analysis and empirical formulae ( Table 6). The minor difference between the simulated and empirical values was due to the characteristics of the cascade and other influences, including significant flood attenuation between the dams, the improved spillway capacity, and the rather higher dam body resistance against erosion. The flood routing simulation continued down to the Opava River, which is located in a wide floodplain with road and railway embankments that create large "reservoirs" when flooded by the dam-break flood. This resulted in a relatively significant reduction of the 557 m 3 /s peak outflow from the Pocheň dam to approximately 310 m 3 /s at the junction with the Opava River ( Figure 12). This peak discharge is lower than the 100-year return period discharge Q100 = 352 m 3 /s of the hydrological flood in the Opava River downstream of the junction with the Čižina. As recommended in Guideline 14 [51], the dam-break flood simulation was terminated at the profile where the peak discharge induced by the dam break dropped below the 100-year hydrological flood discharge.
A summary of the simulations is shown in Figure 13, where the flood hydrographs are plotted starting with the initial hydrological flood that entered the Lichnov II reservoir (on the left) and induced the dam-break flood progression through the system of three reservoirs along the Tetřeví Stream and the Čižina River. In Table 6, the dam breach peak discharges obtained from simulations The flood routing simulation continued down to the Opava River, which is located in a wide floodplain with road and railway embankments that create large "reservoirs" when flooded by the dam-break flood. This resulted in a relatively significant reduction of the 557 m 3 /s peak outflow from the Pocheň dam to approximately 310 m 3 /s at the junction with the Opava River ( Figure 12). This peak discharge is lower than the 100-year return period discharge Q 100 = 352 m 3 /s of the hydrological flood in the Opava River downstream of the junction with theČižina. As recommended in Guideline 14 [51], the dam-break flood simulation was terminated at the profile where the peak discharge induced by the dam break dropped below the 100-year hydrological flood discharge.
A summary of the simulations is shown in Figure 13, where the flood hydrographs are plotted starting with the initial hydrological flood that entered the Lichnov II reservoir (on the left) and induced the dam-break flood progression through the system of three reservoirs along the Tetřeví Stream and theČižina River. In Table 6, the dam breach peak discharges obtained from simulations are compared with the results of the comparisons with real incidents and with data from empirical equations. The comparison indicates that the use of comparisons and empirical equations for a dam break in a cascade may be confusing, especially if the breached dams are very close to one another. These methods may be used for the comparison of the breach peak discharge failures at single dams, but with some caution due to the significant uncertainty involved.

Discussion
The simulation of a dam break in a cascade consists of two basic procedures, namely single dambreak analysis (Sections 3.1 to 3.3), and the analysis of the dam-break-induced flood between the dams (Section 3.4).
Basically, two failure modes-overtopping and piping-should be simulated at the uppermost

Discussion
The simulation of a dam break in a cascade consists of two basic procedures, namely single dambreak analysis (Sections 3.1 to 3.3), and the analysis of the dam-break-induced flood between the dams (Section 3.4).
Basically, two failure modes-overtopping and piping-should be simulated at the uppermost

Discussion
The simulation of a dam break in a cascade consists of two basic procedures, namely single dam-break analysis (Sections 3.1-3.3), and the analysis of the dam-break-induced flood between the dams (Section 3.4).
Basically, two failure modes-overtopping and piping-should be simulated at the uppermost dam. These modes are completely independent; they occur for different reasons at different places in the dam body, and at different times (Figures 3 and 8). The initiation of piping is governed by the threshold hydraulic gradient when critical shear stress is exceeded. Therefore, piping usually starts at the rising limb of the flood hydrograph and progresses during the arrival of the flood peak at the maximum water level in the reservoir. This may happen before dam crest overtopping occurs. At the moment when the "roof" of the pipe becomes thin (here this is assumed to be when it is equal to the "pipe's" diameter), the overburden collapses, and the simulation switches to the overtopping mode at the location of the original piping. This may result in a temporary drop in outflow during piping. The duration of the piping failure until the maximum breach opening developed was more than 1 h ( Figure 8). The worst-case scenario of an overtopping failure occurred at the receding limb of the inflow hydrograph, at the moment when the reservoir was completely full, and critical velocity along the downstream slope was reached.
The optimization procedure showed that the outflow from a reservoir is predominantly governed by the size of the breach opening, and the reservoir water level in relation to the base of the breach opening. Thus, the peak outflow is limited by the location of the first overtopping, the morphology of the dam profile, and the erodibility of the dam material.
With large dams, the impact of the breaching of an upper dam on the failure risk of a downstream one is usually governed by the ratio between the released volume from the upper reservoir and the flood attenuation volume of the lower one. In the case of small dams, the flood attenuation volume is relatively small when compared to the volume of the arriving extreme flood, with the ratio varying between 1:10 and 1:25. The attenuation of peak discharge usually varies between 5-10%. Therefore, the attenuation effect of small reservoirs is not so important for dam break analysis in the case of a cascade of small dams. The capacity of the bottom outlets at small dams usually has a negligible effect on the flood passing through the reservoir and may be neglected.
At lower-lying small dams in a cascade where the reservoir only provides a limited attenuation effect, the rising limb is steep (Figure 9), and the increase in the reservoir water level is quite rapid (5-10 min). In the case of the scheme investigated here, the overtopping and breaching started shortly after the rapid increase in the reservoir water level. About 1 h was needed for the initiation of the piping, and, therefore, the expected failure mode was primarily overtopping. Therefore, the piping failure mode was not taken into account for the downstream dams.
The results of the presented study prove that flood attenuation between reservoirs located in close proximity to one another is practically negligible and that simplifications can be applied (Figure 9) [52]. Flood attenuation increases with the distance from the breached dam and follows an approximately exponential trend. Our results ( Figure 13) correspond well with the results of parametric simulations published in [53], where the longitudinal slope, the dam-break flood volume, the width, and the "roughness" of the channel between the reservoirs are taken into account.
The empirical formulae tested in this study were derived using real collapses of single dams. When used for peak discharge estimation, empirical formulae derived for a single dam break should be applied carefully, as they may underestimate the peak outflow by up to 10% in the case of a dam cascade.

Conclusions
The paper describes the analysis of a dam break process in a cascade of small dams, including the simulation of flood routing along the adjacent streams. The aim is to fill a gap in the existing research concerning the breaching of small embankment dams in a cascade. This is a specific type of event which may differ from a dam break in a cascade of large dams, where the reservoir attenuation volume may play a significant role. The case study concerns an area along Tetřeví Stream and theČižina River in the Moravian-Silesian Region of the Czech Republic. The area contains three small dams: Lichnov II, Lichnov III, and Pocheň. A simplified parametric model for piping and overtopping erosion with parameter optimization was used for the dam break simulations, while a standard 2D unsteady shallow-water flow model was applied for the flood in the areas downstream of each reservoir. As no calibration data were available, the comparison of the simulation results was carried out with the use of empirical formulae and by means of comparison with real small dam failures.
It is recommended that the solution for a dam cascade be divided into separate sequences corresponding to the dam breaching process and flood routing along the streams downstream of the dams. The location of the first overtopping and piping event should be identified based on the dam geometry, the arrangement of appurtenant works (namely the auxiliary spillway), observed leaks, etc. It is desirable to carry out at least a simplified optimization procedure, which guarantees that a combination of parameters providing the worst possible realistic scenario for the simulated event is obtained. When applying empirical formulae or comparison with real dam failures to a dam cascade, one should be aware that these methods may underestimate the peak discharge by up to 10%.
All flood hydrographs ( Figure 13) showing arrival times and flood zones with water depth ( Figure 14) provide important information for warning and emergency plans, and for evacuation procedures.
Further research should focus on the comprehensive probabilistic analysis of dam failures and flood routing between reservoirs in order to create a basis for the development of risk maps for endangered areas.

Conflicts of Interest:
The authors declare no conflict of interest.