Assessing the Hazard of Flooding from Breaching of The Alacranes Dam in Villa Clara, Cuba

1. Introduction

Floods caused by dam failure or dam breaching represent a hydrological phenomenon of great relevance due to their potential to generate flood peaks significantly greater than those caused by natural events such as intense rainfall [1]. Throughout history, dams have failed due to different events or loading conditions [2]. Generally, most documented dam and levee failures are often linked to exceedances of design water levels according to the United States Army Corps of Engineers [3]. In addition, the scientific community has identified other types of failures, such as piping and internal filtration in the dam inner structure, which results in a piping process [4]. Although dam-break flood hydrographs calculated for these failures usually have lower values than those produced by exceeding the water level, they are still significantly higher than the values generated by discharge through the spillway or the emergency flood drains [4].

Over the past 100 years, there have been around 200 dam failures, resulting in more than 30,000 deaths [2,5]. One example was the failure of the Vaiont Dam in Italy, where 2,600 lives were lost [2]. Another recent dam failure event due to exceeding storage capacity occurred in Derna, Libya, where more than 11,000 people lost their lives [5].

Although a time-based analysis shows that the frequency of large dam failures has decreased by a factor of four over the past 40 years, dam failure events continue to occur with significant frequency [6]. Given the high risk that dam failure poses to nearby communities, studies based on mathematical simulation are often required for risk assessment of new construction. These studies include an assessment of how the new construction would affect river flows, its banks, and floodplains [1,7]. According to Bellos, et al. [8] dam breach analyses are divided into two submodels: (a) the dam breach submodel, which is responsible for producing the flood hydrograph, and (b) the hydrodynamic submodel, which uses the flood hydrograph to determine flood peaks and maximum water depths downstream of the dam.

Studies have developed equations to estimate the size of a dam breach or failure (width, slope, eroded volume, etc.), as well as failure time. These equations were derived from data obtained from earth dams, dams with impermeable cores (e.g., clay, concrete, etc.), and rockfill dams, so they are not directly applicable to concrete dams or earth dams with concrete cores [3]. Notable works in this area include Froehlich David [9], y Froehlich David [10], which performed multiple linear regressions to develop an equation that predicts the peak discharge following an earth dam failure. In 2008, the same presented an improved mathematical expressions based on 74 data points of dam failures (both piping and overtopping), characterizing the phenomenon and analyzing the uncertainty in predicting peak flows and water levels by using Monte Carlo simulations [11]. In addition, the author recommended caution with the results, suggesting their application in situations with a low probability of human fatalities and low potential for property damage. Other authors Xu and Zhang [12], have presented empirical equations to predict dam failure parameters based on 182 failure cases. This study showed the importance of dam erodibility in predicting failure, along with the reservoir shape coefficient and failure mode. In addition, they recommended a multiparameter nonlinear regression model where two case studies are included to demonstrate the applicability of the developed models. Other proposed models include MacDonald y Langridge-Monopolis (1984) and Van Thun y Gillete (1990), which can be found in [13]. Peramuna, et al. [14] reviewed the techniques used in traffic modeling of the flood generated by the phenomenon and explored the different one-dimensional hydrodynamic models (1D) and two-dimensional (2D) hydrodynamic models to simulate the propagation of the flood. USACE [7] explored advantages, disadvantages and differences between these simulations and suggests 2D modeling for some cases, since it can produce better results than 1D modeling.

A well-known tool for flood analysis is the Hydrologic Engineering Center's River Analysis System (HEC-RAS) model. It has been used to analyze water flow (Socas, et al. [15], and to determine floodplains, and design hydraulic engineering solutions such as dam breach studies [16,17,18,19,20]. In Albu et al. (2020) HEC-RAS was used to model a dam breach in a mountainous area, highlighting the importance of topographic features in flood propagation. Meanwhile, Marangoz & Anilan (2022) focused on simulating partial dam failures, demonstrating how different dam breaching configurations affect flood dynamics. [21] discusses a management system designed for dam-break hazard mapping using MIKE-21 in complex basin environments. It outlines methodologies for data integration and simulation of dam failure scenarios, which are crucial for generating hazard maps. The research emphasizes the importance of establishing a database for dam-break hazard mapping and spatializing hydrological calculations to create effective flood hazard maps. Ongdas et al. (2020) applied HEC-RAS in an urban context, assessing the impact of dam failure in densely populated areas. The results helped identify critical areas for the implementation of mitigation strategies. This provided valuable information for designing containment and protective structures, highlighting the importance of integrating urban planning with flood management. Another discussion highlights the use of 2D hydraulic models for producing flood hazard maps [22]. Pilotti et al. (2020) explored the model's sensitivity to different hydraulic and terrain parameters, providing a framework to improve simulation accuracy and to reduce uncertainty in results. Thus, the use of HEC-RAS allows for the development of multiple failure scenarios (overtopping or piping), by varying parameters such as the size of the dam breach, the initial conditions of the reservoir, and the topographic characteristics of the valley. These capabilities are vital for risk assessment and hazard management under potential dam failure.

Cuba has faced unique challenges in dam management due to its geography and tropical climate. According to the Cuban National Institute of Hydraulic Resources (INRH in Spanish), a total of 242 dams have been counted, 238 of which are earth dams, and more than 200 of which were built between 1960 and 1980 [23]. Due to frequent slope failures, several authors [24,25,26] have studied the influence of unsaturated soil permeability on slopes, and slope stability during rapid drawdown in earth dams with partially saturated soils. In addition, the occurrence of seismic phenomena on the east side of the island has also been observed by Urquiza-López, et al. [27] in the area, even without events that endanger the mechanical or structural integrity of the dam. Despite the documentation on the problems of Cuban´s earth dams, studies on flood failure and catastrophic flow events are not well documented. These is a need to update methodologies and to provide helpful information to local authorities. The only guideline for calculating time and dam breach parameters for Cuba is the Cuban Standard NC 974 2013. However, this standard does not offer information on the outflow hydrograph, which is essential for assessing the magnitude of the impact of the failure over time. The absence of this information hinders the understanding of risk involving these events, as well as the ability to prepare comprehensive strategies aimed at minimizing the impact of disasters. In addition, the impact of the phenomenon on a regional scale has not been quantified either, although Stucchi, et al. [28] and Socas, González, Marín, Castillo-García, Jiménez, da Silva and González-Rodríguez [15] reported flooding events due to rainfall episodes, these conditions are not similar to dam failure events.

The primary goal of this study was to create a comprehensive hazard map for a potential breach of the Alacranes Dam in Villa Clara, Cuba, using the HEC-RAS software and a 2D downstream flow mode. This study provides useful information for urban planners of how to simulate flood scenarios with a high degree of accuracy and detail for preparing a strategy for early warning systems and evacuation routes to prevent human losses.

2. Materials and Methods

2.1. Procedure for Obtaining Maps for Reservoir Failure

Figure 1 presents the adapted methodological proposal. The Alacranes Dam failure hazard maps were generated using the methodology described by Ferrari, et al. [29], which follows four basic steps. In the first step, the study area was selected (Section 2.2), then a description of the terrain from the reservoir to the boundary conditions at the end of the model, and Manning's roughness coefficient values were obtained (Section 2.3). In the second step and Section 2.4, the failure mechanism (overtopping or piping) and the relationship between the initial hydraulic conditions (boundary conditions) were adopted. Thirteen modeling scenarios were developed: 6 for piping and 7 for overtopping failure mechanism. The third step is applied to the failure model and the flood traffic model (Section 2.5). To obtain the gap opening, the regression models of Froehlich David [10] Froehlich David [11], Van Thun y Gillete USACE [3], Xu and Zhang [12] in their different variants were used, in addition to a physical model (DL BREACH) integrated in the software (Wu [30]. The fourth step obtained the hazard map from using a 2D flood traffic model that included the topographic characteristics of the area (see as Section 3). The depth, speed and hazard maps for breaking the Alacranes Dam used the criteria reported by the State of New South Wales (Australia) and the Department of Planning Industry and Environment [31].

2.2. Location of Study Area

Figure 2a shows the location of the island of Cuba, highlighting the Villa Clara region represented in orange, while Figure 2b shows the Villa Clara province divided by cities (the city of Sagua la Grande is highlighted in red). Figure 2c represents the Digital Terrain Model (DTM) around the Alacranes Dam near the city of Sagua la Grande. The Alacranes Dam is located 7 km southwest of the city of Sagua la Grande and 2 km west of the town of Sitiecito at coordinates E: 590 165, N: 324 733 according to the Lambert NAD 27 Cuba Norte conic projection. The dam was completed in 1972 and retains the waters of the Sagua La Grande River basin, the second largest basin in the country. The 15-km-long and 5-km-wide reservoir, through an extensive channel system, supplies the local chemical industry, as well as irrigates 196 km² of sugarcane in the Armonía area, 110 km² of pastures in Macún, and 150 km² in Sagua la Chica.

The Alacranes Reservoir is the third largest in Cuba, with a storage volume of approximately 350 hm³. The reservoir has a composite section, with a central part (or core) of clay and outer rocky shoulders. Its maximum water level is 36 m and its normal water level is 32.32 m. The intake structure consists of an upstream pressure gallery, a control tower with a segment gate, and a downstream gallery at free flow rate. The guaranteed delivery is 325 hm³ per year. The spillway is a trench with a practical profile crest with a maximum capacity of 2,400 m³/s and ends in an outlet channel carved in the rock until it joins the river. The city of Sagua La Grande is located downstream of the dam (Figure 3) and is at risk of severe flooding if the dam fails or if there is an uncontrolled release of water. Additionally, there are small towns between the city and the coast, parallel to the river, which may also be affected by flooding, given the flat topographical behavior of the mouth toward the Atlantic Ocean. Multiple factors, including both natural and human-induced causes, contribute to potential hazards that can result in life losses. These factors can include earthquakes, storms, climate change, and structural integrity issues from aging that reduce the dam's safety. In addition, downstream communities and cities located directly downstream of the dam are located within the inundation zone. Understanding the interplay of these factors is crucial for effective disaster risk management and for developing strategies to mitigate the potential for life losses. Therefore, this location was selected for developing the province's first hazard map associated with possible reservoir failures.

2.3. Topographic Data

In order to describe the terrain elements that may interact with a flood caused by dam failure, the adoption of a high-resolution DTM is required [29,32]. For the Alacranes Dam, the local DTM developed by the Cuban Geographical Studies Company (GEOCUBA) with a 12.5 m spatial resolution was used (See Figure SI1). Figure 4 represents the local DTM with modifications extracted from the GEOCUBA database. The artificial contour lines within the Alacranes reservoir based on height-volume curves are observed, allowing an approximate knowledge of the spatial variation of the bottom. This methodology is recommended by USACE [3], which explains that in case detailed bathymetric data are not available and a complete unsteady flow path is still desired, the cross-section data can be modified to match the published height-volume curve of the reservoir. The local DTM included those man-made structures present in the area. Among these structures were the so-called "Puertas de Sagua," which consist of two dikes perpendicular to the river, forming a control section (see Figure 3) that allows approximately half of the spillway's maximum flow to be diverted to the coastal plain. These dikes, as indicated by a topographic survey, have been assessed as being in a poor condition. In this study, this topographic survey was used to replicate the configuration of the dikes within the DTM, thus respecting their current condition as indicated by recent topographic measurements.

In order to calibrate the Manning roughness coefficient, a preliminary report of the flood footprint observed on the left bank of the river in the city of Sagua La Grande was used. This record, obtained since the reservoir's creation, provides a valuable reference for adjusting the Manning value. The estimate represented in Figure 5a corresponds to a maximum flow of 1,100 m³/s and is useful for initially calibrating flood behavior downstream of the "Puertas de Sagua," structures that provide partial protection for the city. To calibrate the Manning roughness coefficient, three main steps were followed using Arcement and Schneider [33] and Te Chow and Saldarriaga [34] together with the procedures of Socas, González, Marín, Castillo-García, Jiménez, da Silva and González-Rodríguez [15], Kiwanuka, et al. [35] and Mohamed, et al. [36]. First, initial roughness values were obtained, as represented in Figure 5b, which covers an overview of the study area. Then, modeling was performed for a flow rate of 1,100 m³/s, comparing the results with the actual curve (Figure 5a) until reaching a maximum of R²=0.821, visible in Figure 5c, which focuses on the city. Finally, the Manning´s coefficients were adjusted over the entire area as presented in Figure 5d, generalized for the study area. The calibrated values, between the horizontal and vertical positions of the curve in Figure 5a and the curve modeled in Figure 5b, although specific to the channel and floodplain with an event of 1,100 m³/s, can be extrapolated to the entire area due to the similarity of both topographic and natural characteristics. The results after adjustments remained within the range specified by Arcement and Schneider [33].

To estimate the initial peak discharge of 10,000 m³/s as a mean approximation, the envelope equations were used as described in USACE [3]. Although this value will not be used in the final model, it allows a sensitivity analysis of the roughness coefficient of the study area, using the control points (CP) represented by red circles (Figure 5b, d). Six control points were located at: Sagua Highway Bridge (CP 1, E: 593187, N: 325162); El Triunfo Bridge (CP 2, E: 595557, N: 331276); Plain near Isabela de Sagua (CP 3, E: 595327, N: 336136); Coastal Plain (CP 4, E: 596889, N: 339727); Coastal plain (CP 5, E: 597301, N: 328262) and Coastal plain (CP 6, E: 591616, N: 336873).

Table 1 shows that for a range of ±20% of Manning's roughness coefficient, the flow velocity and depth values vary by less than 20% (see Table SI1). This indicates that, although the Manning´s coefficient is critical for flood routing, for large discharges its influence is not significant when the water depth exceeds 6 m.

2.4. Hydraulic Conditions and Failure Mechanisms

In order to simulate the breach scenarios due to overtopping and piping, simplified boundary conditions were defined. In the case of overtopping failure, the reservoir level was raised by simulating a flood that raised the water level to the crest, roughly representing an inflow hydrograph at the reservoir's tail. In the case of piping failure, the normal water level was maintained, depending solely on internal factors within the dam wall. Given the proximity to the mouth, a downstream sea level of 0 m was established (the use of high and low tides or waves was not considered because they do not occur frequently in that area), allowing for an initial simulation of the flood propagation. To explore a wide range of conditions in the event of an infrastructure failure, 13 different scenarios were modeled. Table 2 shows the different types of breach analyzed and their formulation. The characteristics of each scenario, such as the location and dimensions of the breach, were defined according to Froehlich David [9,11], Von Thun y Gillete [3], y Xu and Zhang [12]. These empirical formulations allowed for estimating parameters of the development of a breach, such as its final length, propagation speed, and shape, based on attributes of the levee or dam involved, such as its height, length, and material.

The classification of hazard levels was based on the NSW and DPIE (2019) criteria, considering the depth and velocity of the water at the assessment point. In this way, it was possible to develop hazard maps for one of the studied scenarios, following methodologies from recent dam failure studies, which proposed a comprehensive approach that combines 2D hydraulic modeling with HEC-RAS and multicriteria assessment of failure scenarios [29,37]. Similar to Paşa, Peker, Hacı and Gülbaz [37], this study prioritized the generation of flood and hazard maps through simulations that consider key parameters such as breach formation time, breach geometry, and initial hydraulic conditions. However, unlike the approach applied to buttress and earth dams in Turkey, this study adapted empirical breach formation models to the context of Cuban earth dams, incorporating unique topographic features of the Villa Clara coastal plain and local Manning´s roughness conditions. Additionally, the analysis was extended to 13 scenarios to evaluate both overtopping and piping failures, replicating the methodological approach employed by Paşa, Peker, Hacı and Gülbaz [37] in the analysis of multiple failure criteria. The scenario that presented the most critical results was selected.

A color scale was used, ranging from blue (low risk) to red (extreme risk), with four intermediate classifications that allow for an understanding of the magnitude of the hazard on a large scale. The areas identified as "blind spots" in the two-dimensional hydraulic simulation within the flood hazard map were assigned as "low risk". Although the model could not calculate hydraulic parameters in these areas due to their size and dispersion at very shallow depths, the absence of flood risk, even in the form of stagnant pools of water, cannot be ruled out since that can be areas where water would naturally tend to accumulate. Although the limitations of a two-dimensional model do not allow for an accurate representation of these behaviors, it was prudent to consider the degree of hazard under these circumstances. For the preparation of the study, the assessment of different hydraulic variables at the control points was proposed, such as the total flooded area (TFA), the maximum water depth (MWD), the maximum water velocity (MWV), the flood arrival time (FAT), the reservoir evacuation time (RET), the total time of flood (TTF) and the average water velocity (AWV).

2.5. Numerical Model

The Navier-Stokes equations describe fluid motion in three dimensions. However, in the context of channel and flood modeling, further simplifications are required. A simplified set of equations are the shallow water (SW) equations. To derive the SW equations, certain conditions are assumed: the flow is incompressible, with uniform density, and pressure is considered hydrostatic. Another assumption is that the vertical length scale is much smaller than the horizontal length scales. As a consequence of these assumptions, the vertical velocity is small and the pressure is hydrostatic [3]. Furthermore, the original Navier-Stokes equations are averaged using the Reynolds number to approximate turbulent motion by eddy viscosity [3]. To improve simulation time, a coarse underlying mesh was used, within which finer topographic details can be extracted [16]. HEC-RAS software (version 6.3.0), developed by the United States Army Corps of Engineers (USACE) was used. The software was designed to solve both the full 2D Saint-Venant equation and the 2D diffusion wave equations. This analysis represents a shallow water model, with the diffusion wave approximation of the shallow water equations (Eq. 1 and 2) specified below [3]:

$${\displaystyle \frac{\partial \mathit{H}}{\partial \mathit{t}}}-\nabla ·\mathit{\beta }\nabla \mathit{H}+\mathit{q}=0$$

(1)

$$\mathbf{w}\mathbf{h}\mathbf{e}\mathbf{r}\mathbf{e}:\mathit{\beta }={\displaystyle \frac{{\left(\mathit{R}\left(\mathit{H}\right)\right)}^{5/3}}{\mathit{n}{\left|\nabla \mathit{H}\right|}^{1/2}}}$$

(2)

where H is the water surface elevation, t is time, the differential operator $\nabla$ s a vector of two-dimensional partial derivative operators given by $\nabla =(\partial /\partial x,\partial /\partial y)$, $\beta$ is the momentum coefficient, $\nabla H$ represents the gradient of the water surface elevation, q represents a lateral flow term, R is the hydraulic radius, and n is the Manning roughness coefficient.

For the simulation, all data were imported into HEC-RAS and converted to the Hierarchical Data Format (HDF). Due to the large size of the study area, it was modeled using an unstructured 45 m x 45 m grid on the plain and a smaller grid on the river channel, composed of a total of 314,426 cells downstream of the reservoir, while the dam was modeled as a spillway. This captured all the details of the study area while also ensuring the computational efficiency of the model. In order to design the dam within the HEC-RAS environment, it was necessary to enter the different dimensional parameters required in the module. As indicated in [16], the grid size is directly related to the accuracy of the results when simulating flows. With a high-resolution DEM available, preliminary tests were performed to identify the appropriate cell size. Initially, a 75 m x 75 m mesh was tested and then progressively reduced to 25 m x 25 m. It was found that 45 m x 45 m was the inflection point where the error for the overall volume did not decrease significantly and instead, the computational time began to increase significantly (See Figure SI2).

3. Results

After running the simulation for the 13 potential dam failure scenarios, the most relevant parameters were extracted in order to emphasize the comparison between scenarios. The most relevant parameters involved (a) TFA, (b) MWD at Sagua Road Bridge (CP1), (c) MWV at Sagua Road Bridge (CP1), (d) FAT from the time the breach occurs, until the flood peak reaches El Triunfo Bridge (CP2), and enters the floodplain downstream to Sagua La Grande, (e) RET; (f) the total time of flooding Sagua La Grande city (TTF), and (g) AWV at El Triunfo Bridge (Table 3).

The results of the comparative analysis between the different dam failure scenarios indicated that the total flooded area varies considerably between the scenarios, with the HEC-RAS physical model (Scenario 13) showing the largest affected area with 604.6 km², while Von Thun and Gillette (1990) A* (Scenario 9) from piping presented the smallest flooded area, with 501.7 km²; this range shows the differences in methodologies and their impact on the magnitude of the flooded areas. Likewise, the maximum water depth in CP1 varied notably with the HEC-RAS physical model reaching 12.47 m, in contrast to the Von Thun and Gillette (1990) B* from piping that reports 8.58 m (Scenario 10). These values indicated the potential severity of the flood depending on the scenario considered. The maximum water velocity and flood arrival time also showed a wide variability of FAT between 1.33 h and 3.33 h, this range is crucial for emergency response planning.

In terms of average water velocity, the HEC-RAS model (Scenario 13) predicted values of 3.21 m/s, implying a significant destructive potential compared to the channel erosion velocities specified in Te Chow and Saldarriaga [34]. Overall, the consistency of results across different methods and models emphasized the importance of considering multiple approaches for an accurate and robust assessment of dam failure hazard. Table 4 shows the results related to dam breach formation such as: average breach width, time to failure, total duration of each simulation, error in estimating the total generated volume, and maximum peak flow during the breach. It was observed that the widths of the breach base vary significantly, from 85 m to 350 m. In the case of the overtopping scenarios, it should be noted that the average breach width value calculated in NC 974-2013 for the Alacranes Dam was 250 m, being within the range of values obtained, and the breach development times ranged from 0.40 h to 11.88 h, while in the NC standard the time obtained was 6.2 h. The overall volume accounting error was generally low, with values from 0.004% to 2.35%. Maximum flows (Q_max) also showed wide variability, with values from 6,860 m³/s to 35,726 m³/s. For all scenarios, the Courant coefficient remained constant between 1 and 0.45. These results reinforce the conclusion of including several types of analysis in these studies.

The 13 scenarios analyzed revealed the exposure of buildings to backwater flooding. Figure 6 shows the saturation of the floodplain with floodwater, highlighting Scenario 13 as the most critical predicted by the HEC-RAS DL Breach physical model. In addition, Figure 6 shows the city of Sagua la Grande, a populated town located downstream the dam. Both averaged and maximum depths were significant (≤10 m) due to the greater destructive potential of flood waves in deeper water layers

A relationship was found between the maximum break discharge variables Qmax and the Total flooded area (TFA) with an R² of 0.90. Despite the limitations from establishing this relationship from only 13 scenarios, the equation that relates to both variables was:

$$\mathbf{T}\mathbf{F}\mathbf{A}=173.9{\mathit{Q}}_{\mathit{m}\mathit{a}\mathit{x}}^{0.1228}$$

(3)

Between flow rates of 7,000 m³/s and 35,000 m³/s, by applying this potential equation it was found that the TFA would range between 500 km² and 600 km². The most critical simulated scenario involved the formation of a 350 m gap in just 0.67 h (Scenario 13).

Figure 6. The Alacranes reservoir flood map for the event of a dam failure obtained for Scenario 13. To the right is the city of Sagua la Grande, with CP 2 marked in the center.

Figure 6. The Alacranes reservoir flood map for the event of a dam failure obtained for Scenario 13. To the right is the city of Sagua la Grande, with CP 2 marked in the center.

Figure 7 shows the velocity map obtained from the simulation of Scenario 13 (see Video SI1). Water flow presented erosive velocities between 3 m/s and 5 m/s near structures such as bridges and protective dikes. The analysis at different CP and at the flood protection structures known as Puertas de Sagua allowed to establish that a dam breach greater than 10,000 m³/s during the peak of the flood would overflow these protective structures. If the cause of a dam failure were piping occurs, there is a probability that the overflow threshold of these dikes could not be reached. However, if this threshold were exceeded, the flood would overflow these protective structures. The coastline would be affected for an estimated length of 105 km, and several settlements and communities within the flood plain would be somehow impacted. It was preliminarily estimated that more than 49,000 population could be exposed to this risk, affecting residents mainly in the city of Sagua La Grande (38,773 inhab.), as well as in the towns of Isabela de Sagua (2,963 inhab.), Sitiocito (3,871 inhab.), La Rosita (1,635 inhab.), Nueva Isabela (980 inhab.), Dos Amigos (289 inhab.), Playa Uvero (200 inhab.), Playa Piñón (24 inhab.) Caharatas (633 inhab.), among other villages and isolated individual homes (population data were obtained from the Cuban Collaborative Encyclopedia (EcuRed) [38].

Figure 8 shows the hazard map developed for the study area with six different color levels defining red as extreme risk, yellow as severe risk, light green as significant risk, dark green as moderate risk, light blue as caution, and dark blue as low risk. The simulation revealed that, in the event of dam failure, the areas adjacent to the reservoir will experience significant flooding, affecting both critical infrastructure and local communities. The areas identified as being at extreme risk are concentrated in the lower valley, where the city of Sagua la Grande, and towns such as Sitiecito, Nueva Isabela, Isabela de Sagua, Dos Amigos, Playa Piñon, and Playa Uvero are located. This finding is consistent with previous studies indicating that topography directly influences water propagation following a structural failure of the dam. It is important to emphasize that the city of Sagua la Grande, with a population of 38,733, would be one of the most affected areas, as it is in the severe to extreme risk classification zone. Likewise, the town closest to the dam, Sitiecito, is in an area with this same hazard classification. Identifying these areas allows for prioritizing mitigation actions, such as creating evacuation routes and implementing early warning systems, which could reduce human loss and property damage.

Table 5 shows the percentage of urban areas near the reservoir along with their flood hazard classification. The city of Sagua La Grande has 91% of its area under a flood hazard classification of severe to extreme flooding risk, which was expected given its location downstream of the reservoir. Similarly, the town of Sitiecito, has 63% of its area classified as extreme risk for flooding. In the case of the town of Isabela de Sagua, since it is located on the coast, the hazard classification ranges from moderate to significant risk. Other populated areas classified as caution and moderate risk for flooding include La Rosita, Playa Uvero, Dos Amigos, and Playa Piñón. Even the town Caharatas, located far from the dam breach, could be affected by the rising water level and humidity of the surrounding soil.

4. Discussion

In this study, the methodology applied in Ferrari, Vacondio and Mignosa [29] was adapted to the local conditions and data limitations of this particular study area in Cuba, starting with hydrological and hydraulic analyses until obtaining the dam failure hazard map similar to those developed by Paşa, Peker, Hacı and Gülbaz [37]. Using the hazard model reported by NSW and DPIE [31], the flood hazard in the area was analyzed and mapped for better understanding of the extent of the affected areas. The depth and velocity values produced by the simulations (Figure 6 and Figure 7) can lead to significant riverbed erosion and landslides. Therefore, the river geometry can change during the flood event as reported in [20,29,32,39]. Furthermore, mud and debris carried by the current can affect the dynamic conditions of the flow, which turns the water into a non-Newtonian fluid [3]. However, the results obtained do not take this into consideration; models with greater capacity and improvements in the implementation of roughness coefficients and/or non-Newtonian models are needed for future work.

Simulation results reported that dam failure affects traffic roads and nearby populated areas such as the city of Sagua La Grande and nearby towns and communities such as Carahatas and Uvero and Isabela beaches in Sagua. Based on the use of rupture models, the most catastrophic event was determined to be the sudden failure of 350 m in 0.67 h (Scenario 13). Although the scenarios related to a piping type dam failure were the least dangerous, simulated flows up to 8,000 m³/s could produce an overtopping of the protective dike at the entrance to the city of Sagua La Grande (See Figure SI4). The results are consistent with the events recorded by USACE [3]; Fattorelli and Fernández [1] and Aureli, Maranzoni and Petaccia [6], which suggested that piping events are less dangerous compared to overtopping dam failure. For the study area, flood depth, flow, and flow velocity were determined to be as high as 12.47 m, 35,727 m³/s, and 5.55 m/s, respectively, at the Sagua La Grande Highway Bridge (CP1) cross-section. Furthermore, results showed that flooding could reach a maximum flow of 35,727 m³/s in 1.13 hours at the same cross-section of the El Triunfo Bridge in the city of Sagua La Grande.

Comparing these results in Table 6 with simulations performed on dams with similar characteristics, it was observed that for dams of similar size and volume, the results of scenarios 1-12 are consistent with other studies [1,2,3]. An example is the Oros Dam in Brazil, with a height of 35.4 m and 700 Mm³, which produced a discharge flow close to 10,000 m³/s. However, Scenario 13 exceeds the values recorded in similar episodes, which is significant since other dams smaller than the Alacranes Dam also increase their peak discharge, mainly influenced by the volume and height as it may occur with the Malpasset Dam in France. These differences can be attributed to several factors including the dam height, crest length, the volume of water stored, data quality, and model accuracy. Furthermore, geographic conditions such as slope, land cover, etc. are other factors that may lead to varying flood depth and velocity as the reported values in studies presented in Table 6.

It is important to consider that simulations are based on hydrological models that, although accurate, cannot predict all variables of water behavior in a real-life scenario [40]. Factors such as soil erosion and climate variability, as well as the presence of atmospheric events, can influence the results. Therefore, follow-up studies are recommended to validate and adjust the model to combinations with intense rainfall, drought, or other initial boundary conditions, including breach parameters.

The use of web-based mapping tools such as DWR and SDSOD [45] to map reservoir breach flooding and flooding hazards in different locations across the state of California is an example of how scientific research can be integrated with state agencies responsible for public safety. Furthermore, future research could explore the effectiveness of different mitigation strategies in specific areas, which would contribute to more informed and effective urban planning. The flood velocity and hazard maps developed in this study indicated exposed areas with high population density. By delimiting these areas and the number of inhabitants affected by the potential of the Alacranes Dam failure this would allow for establishing a relationship between the total flooded area and peak break runoff. While this is an empirical equation and is only valid for use in the study area with limitations, it is a first approximation of a regional study in the province of Villa Clara, Cuba. It is believed that the methodology presented in this research can be replicated in other areas, as a useful tool for preliminary risk assessment for flooding events. This could provide insights and useful information for further implementing a warning system for emergency and evacuation plans for communities living in flood-prone areas.

5. Conclusions

The results of hazard maps for the Alacranes Dam failure scenarios provided valuable information for risk assessment and urban planning by identifying areas effected by potential flooding from dam breaching events. This could provide useful information to develop strategies to minimize the risks associated with dam failure. The simulation of multiple failure scenarios using different mathematical models showed significant variations in parameters such as flooded area, maximum depth and maximum flow velocity, which depended of the type of structural failure of the dam. According to the HEC-RAS physical model, the most critical failure scenario was Scenario 13, with the formation of a 350 m breach in just 0.67 h. The depth, velocity, and maximum discharge could be a significant hazard, endangering important infrastructure and several populated areas downstream of the dam.

For the first time, depth, velocity, and flood hazard maps were generated for the affected areas downstream of the Alacranes Dam in Cuba, allowing the identification of populated areas potentially affected by a potential failure event. Further improvements in methodologies are needed, as well as more detailed studies on the impact of potential dam failures, considering factors such as process dynamics, interaction with channel geometry, and real-time data. These types of studies are useful for emergency planning by identifying exposed areas and potential impacts to minimize the risks associated with dam failures.