Potential Dam Breach Analysis and Flood Wave Risk Assessment Using HEC-RAS and Remote Sensing Data: A Multicriteria Approach

Abstract

Dam breach has disastrous consequences for the economy and human lives. Floods are one of the most damaging natural phenomena, and some of the most catastrophic flash floods are related to dam collapses. The goal of the present study is to analyse the impact of a possible failure–collapse on a potentially affected area downstream of the existing Bramianos dam on southern Crete Island. HEC-RAS hydraulic analysis software was used to study the dam breach, the flood wave propagation, and estimate the extent of floods. The analysis was performed using two different relief datasets of the same area: a digital elevation model (DEM) taken from very high-resolution orthophoto images (OPH) of the National Cadastre and Mapping Agency SA and a detailed digital surface model (DSM) extracted from aerial images taken by an unmanned aerial vehicle (UAV). Remote sensing data of the Sentinel-2 satellite and OPH were utilised to create the geographic information system (GIS) layers of a thorough land use/cover classification (LULC) for the potentially flooded area, which was used to assess the impact of the flood wave. Different dam breach and flood scenarios, where the water flows over man-made structures, settlements, and olive tree cultivations, were also examined. The study area is dominated mainly by three geological formations with different hydrogeological characteristics that dictated the positioning and structure of the dam and determine the processes that shape the geomorphology and surface roughness of the floodplain, affecting flow conditions. The results show that the impact of a potential dam break at Bramianos dam is serious, and appropriate management measures should be taken to reduce the risk. The water flow downstream of the collapsed dam depends on the water volume stored in the reservoir. Moreover, the comparison of DSM and DEM cases shows that the detailed DSM may indicate more accurately the surface relief and existing natural obstacles such as vegetation, buildings, and greenhouses, enabling more realistic hydraulic simulation results. Dam breach flood simulations and innovative remote sensing data can provide valuable outcomes for engineers and stakeholders for decision-making and planning in order to confront the consequences of similar incidents worldwide.

1. Introduction

Dams have significantly contributed to the development of our civilisation, providing water for irrigation and domestic and industrial purposes [1,2]. However, due to the massive, stored water volumes, a great threat to the downstream area inevitably exists in the case of a dam break event. A dam breach flood is an abrupt short-term flood; a dam breach is a low-probability, high-risk catastrophe event that can create an abrupt short-term and very destructive flood event, with significant economic and social impacts on both neighbouring and downstream areas [1,2,3,4]. Fundamentally, dam breach is prescribed as the opening formed in the dam body that leads to dam failure, making the stored water to flow uncontrollably towards downstream regions [5].

According to Penman [6] and the British Dam Society [7], dams are divided into concrete dams and embankments, depending on their construction material. Until the 18th century, all dams worldwide were exclusively embankment dams. During the 19th-century, concrete technology and new manufacturing methods were developed, creating new types of dams [8]. These two categories are also divided into subcategories, according to how they are manufactured, their shape, and how they handle and transport the water loads they receive. There are many more categories into which dams are classified, as their shape and, in general, their characteristics vary according to the region and purpose for which they were made [9].

Even with improved engineering knowledge and enhanced construction quality, a complete safety guarantee is not possible, and many accidents occur [10]. In the 20th century, there have been about 200 dam failures worldwide, with significant consequences [11]. These failures have caused multiple disasters in the downstream valleys in terms of both loss of life and widespread damage to infrastructure and public and private property [12,13,14]. According to International Commission on Large Dams (ICOLD) data [15] covering the period from 1900 to 1975 on dams higher than 15 m (but not including dam failures that occurred during construction or disasters due to wars), the cause of dam failures depends on natural hazards, human activities, the dam type (whether embankment or concrete), and the age of the dam [10,16,17]. For concrete dams, the main cause is foundation failure. For embankment dams, leaks are the biggest cause of failure and the consequent overflow of water over the crest. Moreover, the age of a dam is directly related to its risk of failure. Reinforced concrete dams have an expected life span of several decades, while for embankment dams, the life span can be well over 100 years. For each dam, the critical time is the first filling of the reservoir; however, the ageing period that follows, several years after the construction, is also considered to be of high risk. Other possible causes of failure are earthquakes, landslides, extreme weather conditions (e.g., storms), overtopping, piping, foundation/structural faults, failures or damage, and sabotage (e.g., terrorist attacks, hostilities) [5,14,18,19].

Floods are considered the most catastrophic natural phenomenon in terms of fatalities and property loss. Some of the most damaging flash flood events worldwide are associated with dam breaches that release large volumes of water that affect inhabited areas and human property [14,20]. Over the last few decades, a lot of attention has been paid to flood protection and risk management plans due to the increasing frequency, intensity and extent of severe flood events [21]. The European Union adopted Directive 2007/60/EC for the evaluation and control of flood risks. The fundamental aim of this Directive is the creation of a general outline for flood hazard and risk mapping and modelling to mitigate the unfavourable consequences of floods on the environment, human prosperity, and economic actions [22]. The present study essentially attempts to examine the first stage of Directive 2007/60/EU in the case of a flood, following a hypothetical dam break. The Directive does not explicitly mention this type of flood, but it points out that the first stage of the preliminary assessment is aimed at identifying those areas where flood risk is considered to be significant, either today or in the near future.

Hydraulic flood modelling for events that take place as the consequence of dam breaches, introduced during the 1940s, has faced great expansion during the last two decades; a great number of studies using innovative types of software and models have been conducted [14,23]. One of the most frequently used programs for generating and analysing flood simulations is the free software of Hydrologic Engineering Center’s River Analysis System (HEC-RAS) due to its greatly-enhanced capability of simulating natural rivers and canals [14,24,25]. Previously, HEC-RAS was exclusively one-dimensional (1D) [5], but the latest version has the capability of two-dimensional (2D) flow simulation as well as sediment transport simulation. 1D flood models analyse flows that are assumed to take place in a longitudinal direction, such as rivers and confined channels, but are unable to simulate flood wave side diffusion [24,26]. 2D models are used to simulate the floodplain area, with the hypothesis that the water movement in a perpendicular way cannot be neglected in order to represent the range of floods that 1D models cannot provide [24,27].

Recently innovative tools, useful for geospatial analysis, are the geographical information system (GIS) [28,29] and remote sensing (RS) techniques [20,30], which provide valuable information and data manipulation capabilities of land use/cover classification (LULC), streamflow and topography [31]. In the last decade, the use of advanced unmanned aerial systems (UAS), which now have improved spatial and spectral resolution imaging capability, has facilitated the study of various surface phenomena based on the sophisticated products they provide, such as high precision digital surface models (DSMs) [32,33,34,35]. The UAS-based DSM has been demonstrated for applications in different areas such as hydrology [20,36], engineering [37], archaeology [38], or agriculture (to monitor crops) [35,39,40]. Topography is one of the major factors of the overall accuracy of flood inundation mapping, and digital surface models are unique in creating a hydraulic or hydrodynamic model and mapping the flood extent [41]. A high-resolution digital elevation model (DEM) or digital surface model (DSM) is suitable for extracting the details of topography, acquiring water surface elevations and simulating flood inundation and depth [31,33].

Furthermore, the new innovative Copernicus program of the European Space Agency (ESA) offers precise, timely and easily accessible data for the improvement of environmental management. Each mission of Sentinel satellites emphasises a different aspect of Earth observation (EO) [42,43]. Sentinel 2 (S2) optical satellite data has made accessible a new land monitoring and classification system, with 13 spectral bands (of 10 m spatial resolution), that can provide a quantified analysis of LULC type [44,45,46]. LULC features affect the surface characteristics that affect the nature, extent, and speed of surface runoff [46,47,48]. Several studies have been conducted worldwide using EO information to map LULC through diverse methods and datasets [20,49]. The detailed delineation of LULC is important in order to assess the risks related to dam breaks and mitigate the effects of this hazard [25]. The anthropogenic interference on LULC has greatly changed the river hydrology that determines flood hazards and has increased the factors of vulnerability and exposure to disasters and flood hazards [50].

Various studies have been conducted to explore the efficiency of the synergistic use of hydraulic modelling with UAS data and GIS manipulation to assess flood wave behaviour and possible hazards of a dam breach [5,14,31,33,51,52].

The study area includes the Bramianos river’s watershed area and the homonymous dam situated at the southeast part of Crete island, about 5.5 km northwest of the city of Ierapetra. This watershed was selected mainly due to its characteristics, such as intensive agricultural activity and the short distance between the dam and a residential area in the downstream part of the dam. The agricultural land is mainly comprised of many greenhouses covering a large part of the area; these are of high economic value. Furthermore, the coastal village of Gra-Lygia is located at a distance of 2 km; in the case of flash-flooding, this village will be at immediate risk.

The purpose of this study is to investigate the impact of a possible failure of the Bramianos dam, comprising (a) dam-break simulations in both overtopping and piping states; (b) the behaviour of the flood wave that will result from a possible collapse of the dam, applying HEC-RAS 2D on two different elevation profiles—a DEM and a UAS-derived DSM; (c) the calculation of an accurate prediction of the inundation levels (depth and velocity) and the time of flood wave arrival at a given location, utilising HEC-RAS and GIS analysis; (d) a risk assessment according to the downstream LULC affected by the flood wave, using remote sensing data (UAS and EO) and American Society of Civil Engineers (ASCE) methodology.

2. Study Area

2.1. Watershed Characteristics

Crete is the biggest and southernmost island of the Aegean Sea, with an E–W direction. It is a mountainous island of an elongated shape, with a mountain range containing Lefka Ori, Idi and Talea Ori, attaining a maximum height of 2456 m above sea level. Its length is about 260 km along the W–E axis, while its width varies between 60 km in the middle and 15 km at the Ierapetra graben in the east. It covers an area of 8261 km² and has a coastline extending over 1046 km in length. Crete is considered a semiarid region; the average annual precipitation is estimated at approximately 900 mm [53,54,55]. The study area of Bramianos torrent watershed and the constructed dam is located between 35°0′00″–35°6′00″ North and 25°38′30″–25°45′00″ East, in the southeastern part of Crete, in Lassithi prefecture (Figure 1a–c).

According to the Assessment and Management of Flood Risks (EU Directive 2007/60), it belongs to East Crete River Basin “GR41” and the 13th Water District of Greece.

The general orientation of the stream is from north to south; it flows west of the city of Ierapetra, next to the village Gra-Lygia, a seaside village located 5 km west of Ierapetra and 39 km south of Agios Nikolaos. Gra-Lygia is situated at the exit of the basin, with a population of 1528 inhabitants [56].

The basin of Bramianos covers an area of approximately 25.8 km², with minimum and maximum elevations of 30 m and 940 m, respectively, mean annual precipitation of about 650 mm and a mean annual runoff depth of 62 mm (Figure 1b). According to the archive files of the Greek Ministry of Environment and Energy (

Table 1. Flood events that have occurred in the study’s broader area.

Location Name	X-	Y	Date	Flood Characteristics	Type of Damage
Lassithi Prefecture, Sitia	691,561	3,897,932	20 November 1932	Flash Flood	Economic
Heraklion Prefecture	595,200	3,911,100	18 October 1937	Flash Flood	Economic: Property and Rural Land Use; Human loss
Heraklion Prefecture, Yiofiros	600,589	3,909,814	16 January 1974	Flash Flood	Economic
Heraklion Prefecture	595,200	3,911,100	16 January 1994	Flash Flood	Economic: Property and Rural Land Use
Lassithi Prefecture, Plateau	666,507	3,880,867	04 May 2000	No data	Economic: Rural Land Use
Lassithi Prefecture, Metochion	630,218	3,894,324	01–03 December 2001	No data	Economic: Property
Lassithi Prefecture, Ziros	695,070	3,883,325	18 August 2002	No data	Economic: Rural Land Use
Lassithi Prefecture, Kato Choriou	663,175	3,879,616	26 May 2003	No data	Economic: Property
Heraklion Pref., Agias Pelasgias	592,188	3,918,422	06 November 2004	No data	Economic: Property
Heraklion Prefecture, Viannos	629,050	3,879,217	13 November 2010	Flash Flood	Economic
Heraklion Prefecture	595,200	3,911,100	19 November 2020	Flash Flood	Economic: Property and Rural Land Use

), several significant flood events (not associated with a dam break) have occurred in the broader area, mainly due to extreme rainfall episodes [57].

2.2. Dam Characteristics

The Bramianos dam was constructed in the early 1980s according to the Final Design of 1974 and started operating in 1986. The dam’s lake covers an area of 105 ha and has a capacity of 15 × 10⁶ m³; it is the second-largest wetland of southern Greece after Potami dam at Amari valley. The reservoir of the dam serves the irrigation needs of the great valley of Ierapetra and is enhanced by the waters of Kalamafkianos and Myrtos torrents and the waters of the Malavra springs. It is an embankment dam; the main body consists of clayey–sandy silt. The upstream rockfill is made of crushed quarry stones. The upstream slope of the dam body is 2.6/1.0 (Horizontal/Vertical), and the downstream slope is 2.4/1.0 [58].

According to the sections of the dam design, the dam consists of upstream transitional zones, namely, coarse-grained and fine-grained filters with a total thickness of 1.20 m between the dam body upstream rockfill and the main body (clayey–sandy silt). The height of the dam (from the natural bed of the torrent) is 43.5 m, the length and width of the crest are 540 and 9 m, respectively, the elevation of the crest and the impermeable zone are 75.5 and 74.5 m, respectively, and the volume of the dam is 1.27 × 10⁶ m³ (Figure 2). There is also a spillway consisting of the inlet works, a 200 m long box cross-section conduit of 3.5 × 3.0 m; the outlet works on the Diavata torrent. The inlet works are in the western boundary of the reservoir, about 400 m north of the dam. The spillway is neglected in the case of a dam breach because the maximum flow that it can deliver is too low compared to the breach discharge (Figure 1c) [58,59].

The land use/cover of the dam’s downstream area is mainly composed of greenhouses, olive groves and artificial infrastructures such as the village of Gra-Lygia and the road network.

2.3. Geological and Hydrolithological Characteristics

The island of Crete is situated along the southernmost bend of the Hellenic Arc, where active extension occurs at the overriding plate of the Africa–Eurasia retreating subduction [53,60,61]. It consists of pre-alpine, alpine and post-alpine formations [62,63,64] and is characterised by active tectonics. Both onshore and offshore faults in Crete display significant morphological features, and their strike is both arc-parallel and arc-normal. The tectonics of Crete are dominated by east–west-trending normal faults inherited from pre-alpine tectonics. The active normal faults have a north–south direction and are characterised by sharp fault scarps [63,64,65]. The Ierapetra fault zone is a normal fault, located about 10 km eastwards of the study area. It traverses the entire width of eastern Crete (>20 km), from the northern to the southern coast, with a NNE–SSW direction [54,66]. Along its southern coast and mainly at the eastern part of the island, several uplifted marine terraces (whose total length is almost 24 km on both sides of the city of Ierapetra) have been described [53].

The geology of Crete can be separated into an alpine basement and sediments of Neogene and Quaternary age. The pre-Neogene rock succession of nappes forms a pile and lies between two metamorphic belts. Neogene and Plio-Pleistocene sediments cover approximately one-third of the island and overlay the alpine formations. They present great dissimilarities in age and lithology and are composed mainly of marls, conglomerates, and marly limestones. Quaternary sediments lie over all the older formations and consist of terrestrial, marine–brackish deposits, sand, pebbles, boulders and clays, loose or slightly consolidated [67,68,69].

The Bramianos basin geologically consists of Makrylia and Kalamavka formations that belong to the Early Tortonian Stage and the Ammoudares formation of the Late Tortonian/Early Messinian Stage, both of Late Miocene of the Neogene Period (11.6–5.3 million years ago). The Makrylia formation is characterised by alterations of fossiliferous blush-grey marls within the middle part and a variable amount of poorly lithified graded turbidity sands. The Ammoudares formation consists of alternating sandy, bioclastic limestones as well as yellow-grey homogeneous and laminated marls, while the Kalamavka formation contains a succession of up to 600 m, with a gravelly lower part and a sandy and marly upper part (Figure 3) [70,71].

As far as the hydrogeological properties of the basin are concerned, the Neogene formation includes three main series of rocks with different characteristics. The first series consists of 5 layers of calcareous sandstones between which silty marls intervene. Its total thickness is about 72 m, and it is highly permeable (K = 10⁻³ m/s). The second is a transitional series, has a variable permeability (K = 10⁻⁵ m/s) and consists of a calcareous sandstone layer. Over that layer, there are two horizons of marls and silty marls, 6 and 10 m thick, respectively, and below it, there is a 9-m thick layer of silty marls. Its total thickness is 35 m. In the third series, sandstones, siltstones and limestones of low thickness intervene, which are a few to tens of meters apart, with a changing structure and thickness of over 200 m. Between the layers, there is no connection and water intakes appear in some places following the inclination of the layers. The layers sink under the sea level; therefore, the effect of calcareous sandstones and layers on water losses is low (K < 10⁻⁷ m/s).

Based on the results extracted from borehole drillings during the study of the dam construction, the cementation of the calcareous sandstones was possible but quite expensive. Therefore, the crest of the dam was designed further north and upstream of the boundary of the Ammoudares formation, within the limits of the Makrylia formation, where the water permeability of marls is considered lower or even impermeable.

3. Dam Breach Mechanisms and Flood Wave Impact Overview

The procedure for estimating a potential dam breach mechanism and behaviour and the flood wave consequences requires the collection of information on dam characteristics, such as the type of construction, reservoir volume, and height and length in order to define and evaluate the impacts of a dam failure. Furthermore, the dam failure scenario, including type of disruption, time of year, time of day and weather conditions, plays a crucial role with regards to flood wave behaviour and effect.

Numerous studies of dam failures have indicated that overtopping and piping are the two major failure modes of earth and rockfill dams [18,33,72,73,74]. The other causes of dam collapse are seismic vibrations, slop sliding, and sabotage incidents, generally leading to the occurrence of the above overtopping and piping scenarios or less severe scenarios [75,76]. Overtopping failure is defined as water flowing over the dam crest, which leads to the erosion of the channel along the downstream face of the dam. Dam crest overtopping can occur in extreme flood events, in combination with mismanagement of the spillway, landslides of unstable slopes of the reservoir and the trapping of trees or other objects in the spillway, resulting in the insufficient operation of the dam. Piping is an internal erosion process that leads to the occurrence of spot leakages on the body of the dam or its foundation and can be classified into backwards erosion, interior erosion, tunnelling, suffusion and heave [18,72,73]. Therefore, these two dam-failure causes are selected and simulated in this study. For each scenario that is examined, the theoretical background of the mathematical models, the description of the data and the results (as exported from the applied models) are developed and presented.

The following step concerns the identification of the flooded area, which depends on shape and size of the basin, landscape, dam height, reservoir volume, LULC and dam failure scenario. The water depth and velocity constitute the most important factors to understand the water flow and flood procedure, which depends on the failure scenario, distance from the dam, and the morphometry of the inundation area. They directly affect the potential for damage to structures, loss of life and impact on the environment. As the depth of water inside structures increases, the damage increases. Nevertheless, even shallow water flowing at high velocity can considerably damage a structure, especially when the water carries debris. Excessive water velocity can lead to increased erosion and the loss of environmental assets [77].

The final stage regards the identification of the flood wave’s direct or indirect impact, determining the level of damage to structures, infrastructure and potential economic losses or loss of life, which depend on the depth and velocity of the water. The level of damage is not certainly an indicator of the potential for loss of life. Aspects that can affect the probability of loss of life include the distance between the dam and the infrastructure–settlements as well as the existence of a warning system to alert the community of an actual dam failure in order to evacuate the area on time.

4. Materials and Methods

4.1. Orthophotos and Sentinel-2 Data—LULC Mapping

For detailed delineation of the basin’s downstream area LULC, remote sensing data were obtained, such as the very high resolution (1 m) orthophoto dataset acquired in 2014 from the National Cadastre and Mapping Agency SA. Moreover, a Sentinel-2 image, cloud-free and atmospherically corrected (Level 2A), was collected free of charge via the ESA portal (https://scihub.copernicus.eu/; acquisition date: 18 May 2018). The UAV images and the orthophotos were utilised to distinguish the LULC types and validate and update the S2-based classified map. S2 data preprocessing involved image resampling (10 m, since S2 images have different pixel sizes), applying a subset (given the size of the study area, there was no need to process the entire image; thus, it was clipped for processing reasons to a more manageable size) and masking (pixels outside the basin limits were excluded). Finally, a stack (collocation) of 10 spectral bands for each date was executed by excluding Bands 1, 9, and 10 [20,62].

Subsequently, the classification of the data was performed. The S2 image was classified utilising the random forest (RF) classifier, a nonparametric regression tree machine learning algorithm [78]. In general, the RF classification method demonstrates good simplification by aggregating the decisions of individual trees. It can also handle large numbers of features (input datasets), given the inherent feature selection in the model generation process. Finally, it is comparatively fast to train and, by extension, set to predict [78,79,80].

The postclassification process included the processes of filtering (majority filter) to remove isolated pixels, smoothing the ragged class boundaries and generalising classified output by removing isolated regions. Moreover, using the high accuracy images of the UAV, cross-checked with the orthophotos, the misclassified spots were minimised [81]. The classification training samples (80 points) for the classification process (70%, 56 samples) and the accuracy assessment (30%, 24 samples) were collected using farmers’ declared parcels (from the year 2018), the orthophoto images and Google Earth images (of the same year). Using the first training sample set, five classes were determined: settlements, artificial areas, bare soils and rocks, olive groves, greenhouses, and forest–shrubs.

4.2. Unmanned Aerial System (UAS)—Digital Surface Model

The unmanned aerial system (UAS) that was used consisted of the SenseFly eBee unmanned aerial fixed-wing vehicle that is able to cover up to 12 km² in a single automated flight and can acquire images with a ground sampling distance (GSD) of as low as 1.5 cm per pixel. Additionally, a Canon S110 RGB camera of 12-megapixel resolution and a ground control base and a conventional portable computer from which the flights could be planned and controlled with E-Motion software were used. For aerial photography, 10 flights of a total time of 4.5 h were conducted, and 2395 images were taken in order to cover a total area of 338 ha [32,35].

The lateral and longitudinal overlap between the images was programmed at 60% and 80%, respectively, and the flight altitude was defined so as to achieve a spatial resolution of between 5 and 8 cm/pixel. The postflight processing of the aerial images was done by Pix4D software using ‘Structure for Motion’ (SfM), which is a photogrammetric range imaging technique for the reconstruction of three-dimensional structures from a series of two-dimensional images coupled with local motion signals. Through SfM, a detailed 3D point cloud and a DSM (spatial resolution of 5 cm/pixel) of the surface of the area of interest were produced after the photogrammetric processing of the aerial images by applying triangulation and stereo-pairs matching [32] (Figure 4a,b,d). The DSM represents the highest elevation of the terrain and, therefore, follows the vegetation canopy, the roofs of buildings, greenhouses, and the ground when there is nothing else on top of it. The DSM was used as input to the HEC-RAS hydraulic simulation model to investigate whether it improves the final simulation result, compared to the digital elevation model (DEM) [82,83].

4.3. Digital Elevation Model

The DEM was derived from the stereo pairs of very high accuracy photogrammetric aerial images of the National Cadastre and Mapping Agency SA. It is a numerical representation of the Earth’s surface that contains actual height points, representing the topography, as well as the method to calculate elevations between the height points. Typically, a DEM describes the surface of the Earth through a regular rectangular grid, which can be visualised as a wire grid framework draped over the surface of the Earth. Each grid intersection is defined by two components, i.e., the coordinate reference grid location (y, x) and the associated height (H). The x and y coordinates are referenced to the datum and coordinate systems, while the height component is referenced to the mean sea level (Figure 4c).

4.4. Simulation Models of Dam Failure Mechanisms

The main values for overtopping and piping modelling failure mechanisms are related to the geometrical characteristics of the breach formed, which are (a) mean width (Bav), (b) bottom width (Wb), (c) height (hb), and (d) the height of the water from the bottom to the maximum free surface level (hw). These parameters were calculated utilising the empirical equations of HEC-RAS and validated with the corresponding developed data from the Special Secretariat for Water (Greek Ministry of Environment and Energy) [84].

HEC-RAS software was used to calculate the development of the breach. There are two breaching approaches to choose from—user-entered data (UED) and simplified physical—with regards to the way that the input data are manipulated. For our purposes, the UED methodology was used for estimating geometrical and time characteristics, which incorporates regression equations based on observations of data from real failures. Most of the regression equations were developed from a small subset of the existing database of dam failures (108 historic dam breaches listed in the US Bureau of Reclamation report). The dams included in the analysis are a mixture of homogenous earthen dams and zoned earthen dams [85]. Thus, the pick of the appropriate regression equation should be made with caution.

According to the analyses of the US Hydrologic Engineering Center (HEC) report [86], there are several available techniques (which were taken under consideration in our study) such as (a) Froehlich and David 1995a, 1995b, 2008, (b) MacDonald and Langridge-Monopolis 1984, (c) Von Thun and Gillette 1990, and (d) Xu and Zhang 2009 [85,87,88,89,90]. These studies provide detailed information regarding the selection of the most suitable model. The use of many different sets of equations to obtain a range of values for each parameter and to evaluate the corresponding output hydrographs produced during the failure has been suggested. The equation sets have been produced from a great variety of constructed projects of different scales.

In the present study, Froehlich’s method [88] was selected as the most recent and representative method of this research after considering the HEC report and other evaluation research [91]. According to this method, Froehlich updated his breach equations based on the addition of new data, having utilised 74 earthen, zoned earthen, earthen with a core wall (i.e., clay), and rockfill datasets to develop a set of equations to predict average breach width, side slopes, and failure time. The height of the dam and the volume of the water at breach time, used for his regression analysis, ranged from 3.05 to 92.96 m (with 93% <30 m and 81% <15 m) and 0.0139 to 660.0 m³ × 10⁶ (with 86% <25 m³ × 10⁶ and 82% <15 m³ × 10⁶), respectively. The regression equations for average breach width and failure time are presented as Equations (1) and (2).

$${\mathrm{B}}_{\mathrm{ave}}=0.27{\mathrm{K}}_{\mathrm{o}}{\mathrm{V}}_{\mathrm{w}}^{0.32}{\mathrm{h}}_{\mathrm{b}}^{0.04}$$

(1)

$${\mathrm{t}}_{\mathrm{f}}=63.2\sqrt{\frac{{\mathrm{V}}_{\mathrm{w}}}{{\mathrm{gh}}_{\mathrm{b}}^2}}$$

(2)

where B_ave is the average breach width (m), K_o is a constant (1.3 for overtopping failures, 1.0 for piping), V_w is the reservoir volume at the time of failure (m³), h_b is the height of the final breach (m), g is the gravitational acceleration (9.807 m/s²), and t_f is the breach formation time (s). According to the results of this analysis, the average side slopes should be 1.0 H/1 V for overtopping failures and 0.7 H/1 V for piping and seepage.

4.5. Hydraulic Model Analysis

The data required to solve the Bramianos dam break scenarios on HEC-RAS software are as follows: (a) geometric features along the valley downstream of the dam, (b) morphology of the area, (c) hydraulic characteristics of the valley (roughness, contractions, expansions), (d) areas of interest for the assessment of the impact of the flood wave, and (e) initial and boundary conditions of the two-dimensional routing area in terms of flow and level hydrographs.

The method used to route the inflowing flood hydrograph through the reservoir was level pool routing because of the wide and short reservoir shape. The connection of the reservoir and the downstream 2D area was made by SA/2D area connection, which represents the dam body.

Manning’s values determined the surface roughness and were assigned according to the different land cover types, as classified by the remote sensing data. The appropriate determination of Manning’s values regarding the existence of buildings, greenhouses, fences, and vegetation is crucial to providing a more realistic simulation of water flow characteristics [92].

In the present study, the simulation was based on a 6-h-rainfall at basin-level with a returning period of T = 10,000 and the assumption that the initial level of the reservoir is equal to the elevation of the spillway (+73 m). The inflow supply to the reservoir has a maximum of 407 m³/s and a total volume of 2.7 hm³. The capacity of the reservoir is quite large (for a level of +73, it is 15 hm³) in comparison to the volume of the incoming flood; thus, the generated size of the flood wave is more affected by the existing volume of water in the reservoir at the start of the break.

For downstream flood wave propagation, the mathematical two-dimensional flow model of the HEC-RAS 2D was used. This model solves flooding using either the complete set of Saint-Venant equations or the two-dimensional diffusion wave (DW) equations. The first option allows for the more detailed implementation in a variety of problems, and it is frequently used for the downstream flood propagation of dam breaks due to the greater accuracy that is achieved in situations characterised by significant changes in flow characteristics in terms of time and space. On the other hand, the second option offers faster solutions, and, at the same time, shows optimal numerical stability results. The software-integrated algorithm uses an irrational finite-volume format that allows for a larger computational step over rational methods. The format of finite volumes is characterised by increased stability and reliability of the solutions and realistically deals with the phenomena of “drying” and “wetting” of the cells of the flow domain. In this study, we chose to use the DW approximation for downstream flood wave propagation due to its stability and the exact representation of the flooded areas. However, a scenario using the complete set of Saint-Venant equations was also implemented and compared with the results of the DW approximation to verify that the accuracy obtained with the selected approach was sufficient for the current application. The comparison was made for the overtopping scenario using the DEM.

More specifically, the 2D DW computational method allows the computations to run faster and with greater stability. Most 2D modelling situations, such as flood modelling, can be accurately modelled using this solver, where inertial forces tend to dominate frictional and other forces. The DW computational method can be used in situations where the flow is mainly driven by gravity and friction, and fluid acceleration is monotonic and smooth (i.e., no waves). The general form of Equation (3) is as follows:

$$\frac{\partial \mathrm{h}}{\partial \mathrm{h}}+\mathrm{C}\left(\mathrm{h}\right)\frac{\partial \mathrm{h}}{\partial \mathrm{x}}+\mathrm{D}\left(\mathrm{h}\right)\frac{{\partial }^2\mathrm{h}}{\partial {\mathrm{x}}^2}=0$$

(3)

where C (h) > 0 is the transport velocity of the wave propagation, while D (h) is the diffusion coefficient and h is the flow depth.

The main limitation of the model concerns the size of the grid in relation to the resolution time step. The time step must be selected in such a way that a maximum limit is preserved by the user, while, at the same time, the CFL criterion is met (u Δt/Δx ≤ 1, where u is the velocity of the wave propagation and Δx, Δt are the spatial and temporal discretisations, respectively).

In the present study, the two scenarios of dam failure (piping and overtopping) were simulated with two different terrain models. A digital elevation model (DEM) taken from the National Cadaster and a digital surface model (DSM) taken with the use of an unmanned aerial system (UAS). The results produced were maps of the flooded area, flow depth, flow velocity and risk. Finally, the results were compared in terms of scenarios and relief. The DEM refers to the geometric information of soil morphology, while the DSM also refers to all the surfaces found on that soil (anthropogenic elements such as buildings, greenhouses, infrastructure, natural vegetation, and crops), which were factored into the calculations by our model (Figure 4c,d).

4.6. Flood Hazard Assessment

The intensity of a flood is considered a widely accepted indicator for flood hazard estimation [93]. Intensity refers to the attributes of a hazard that causes damage and is commonly calculated utilising water depth and flow velocity (Equation (4)) [93,94].

Flood intensity = water depth × flow velocity,

(4)

Several criteria for potential-damage assessments have been developed worldwide, such as those of the American Society of Civil Engineers (ASCE, 1952) [95] and the US Bureau of Reclamation (USBR, 1988) [96], which were also used by the Greek Ministry of Environment and Energy for flood risk management plans [97]. In the present study, the ASCE criterion was examined, which considers the values obtained by the multiplication of water depth and flow velocity, derived from the hydraulic model. The effects of the collapse of the examined dam were classified into two categories: (a) loss of human lives and (b) material damage. Regarding the loss of human lives, we accept the criterion of ASCE, according to which depth × velocity >2.10 m²/s endangers human lives.

Table 2. Flood hazard classification according to the American Society of Civil Engineers (ASCE) criterion.

Flood Hazard Classification	Flood Intensity Values (Water Depth × Flow Velocity, m2/s)
Low–Medium	<2.1
High	2.1–3
Very High	>3

shows the flood hazard classification according to the ASCE criterion, with the corresponding flood intensity.

5. Results

5.1. LULC Results and Validation

The accuracy assessment of the LULC classification map was accomplished using the second validation set. The results were assessed using overall accuracy (OA) [98] and the Kappa coefficient [99]. The OA and the Kappa coefficient were estimated at 95.83% and 0.89, respectively. The results indicate high accuracy, especially in the cultivated plain area.

From the classification process, it was revealed that the LULC of the Bramianos basin downstream area, which covers a region of 6.09 km², contains 1.61 km² olive groves (26.4%), 1.21 km² settlements (Gra-Lygia, 19.87%), 1.13 km² greenhouses (18.56%), 0.88 km² forest-shrubs (14.45%) and 1.26 km² artificial land, bare rocks and soils (20.69%) (Figure 5). This area constitutes the most vulnerable area that will be directly affected by a probable dam breach flood.

5.2. Floodextents, Depths, Velocities and Arrival Time Estimation for the Two Elevation Profiles and the Two Dam Failure Mechanisms

The comparison of the results obtained by the complete set of Saint-Venant equations and the diffusion wave approximation for the overtopping scenario using the DEM revealed that the accuracy achieved with DW approximation is comparable to that, and it can be considered sufficient for the current application. Specifically, as can be seen in Figure 6a, the maximum flood extents obtained using the two methods are almost identical and have a difference of less than 3%. Furthermore, the estimated discharge hydrographs (e.g., Figure 6b) obtained in both cases in several examined cross-sections are also almost identical (R² greater than 0.99, peak discharge difference less than 1%), indicating that the accuracy of the DW approximation, as implemented in HEC-RAS 2D model, is adequate for the propagation of the flood wave downstream of the dam break, also considering the various uncertainties caused by the required assumptions in similar applications in terms of hydrology, geology, and flow domain characteristics. It should be noted that the DW approximation was used only for downstream propagation and not for the dam breach characteristics. Accordingly, the DW approximation was used in the following analysis due to its lower computational requirements.

After implementing the simulation for the two breach scenarios, the most relevant parameters were extracted to highlight their differences. The most relevant parameters involved flood extents, maximum water depths, flow velocities, and the arrival time from the moment the breach occurs until the first flood wave occurs and enters the Bramianos floodplain area.

From the relevant created hydrographs of the dam’s downstream areas (Figure 7), it is apparent that the overtopping scenario shows a much higher flow rate, equal to 15,550 m³/s over 11,540 m³/s, and a longer arrival time than the piping scenario, 05.00 h over 04.19 h, respectively.

The two dam breach situations were also compared with the inlet hydrograph of T = 10,000 in the reservoir, which constitutes a very low value of the peaks and the volume of overtopping and piping cases. During the processing and analysis of the data on HEC-RAS using the DEM and DSM terrain forms, the flood extent, flow depth, flow velocity and time of arrival maps were extracted for the overtopping and piping scenarios and the results were compared, as shown in Figure 8, Figure 9 and Figure 10.

The flood extent layer comprises the inundation boundary containing the areas covered by the flood during the event. According to the overtopping scenario, the inundated areas were estimated to be about 177.8 and 185.8 ha using the DEM and DSM datasets, respectively. For the piping scenario, the affected areas were calculated to be about 169.2 ha and 172.5 ha, respectively. Thus, the overtopping scenario represents 5% and 7% larger extents of the inundated areas than the piping scenario when using DEM and DSM, respectively. The estimated total flood extent increased by 4.5% when the DSM was utilised instead of the DEM in the overtopping scenario.

Mean, standard deviation and maximum depth values were assessed by evaluating the zonal statistics of the inundated pixels. In the case of the overtopping scenario using the DSM, they were 3.81, 2.86 and 13.47 m, respectively, and, for the piping scenario using DSM, they were 3.43, 2.57 and 12.92 m, respectively. For the DEM dataset, all values were slightly lower than those of the DSM for the overtopping and piping scenarios (Figure 8 and Figure 9).

Another important result of the 2D hydraulic modelling process was the flow velocity estimation. In the case of the DSM for the overtopping scenario, the mean, standard deviation and maximum values of flow velocity were 3.26, 3.4 and 31.7 m/s, while for the piping scenario, they were 3.0, 3.03 and 31.83 m/s, respectively. For the DEM dataset, all values were little lower than those of the DSM for the overtopping and piping scenarios (Figure 8 and Figure 9).

Finally, the arrival time of the flood to the farthest point of the area was longer in the DSM than the DEM for both the overtopping (0.82 h versus 0.55 h) and piping scenarios (0.83 h versus 0.72 h). The two different breach scenarios for the case of the DSM display equal arrival times and, as far as it concerns the DEM dataset, the overtopping scenario displays a shorter arrival time (Figure 10).

5.3. Flooded Areas

The flooded areas were assessed through the interrelation of the LULC map derived from the remote sensing classified data and the flood extent maps derived from the overtopping and piping scenarios. The flooded areas, according to the overtopping scenario, were 177.83 and 185.47 ha using DEM and DSM data, respectively. The flooded areas, according to the piping scenario, were, correspondingly, 169.21 and 172.53 ha. Thus, it is clear that the difference of the flooded areas varies from 1.92% in the piping scenario to 4.12% in the overtopping scenario. Moreover, in the case of the more accurate DSM data, there is a difference of 6.98% among the two scenarios, while, with the use of DEM data, a difference of 4.85% was observed. In

Table 3. Flooded areas for the various LULC types in the affected area.

LULC	DEM		DSM
	Overtopping (I) (ha)	Piping (II) (ha)	Overtopping (II) (ha)	Piping (II) (ha)
Settlements	31.57	30.43	28.90	26.73
Artificial, Bare soils and rocks	8.70	8.23	9.35	8.37
Olives	67.03	64.34	73.00	69.11
Greenhouses	40.99	38.32	41.71	38.46
Forest, Shrubs	29.54	27.89	32.51	29.86
Total	177.83	169.21	185.47	172.53

, the flooded areas are represented regarding the two elevation datasets and the two hydraulic scenarios. From this estimation, it can be noted that the percentage of total flooded areas varies approximately from 15.49% (DSM-overtopping) to 17.98% (DEM-piping) regards the Gra-Lygia settlement. The most affected cover concerns the olives groves, varying from 37.69% (DEM-overtopping) to 40.06% (DSM-piping), and greenhouses, varying from 22.29% (DSM-piping) to 23.05% (DEM-overtopping).

5.4. Hazard Relation to Depth and Velocity Results

Exploiting the extracted results, the outcomes of water depth and flow velocity multiplication were calculated and classified according to the ASCE [95] flood intensity hazard criterion, which is a reliable and generally accepted index for flood hazard estimation. Consequently, the corresponding geospatial hazard maps were extracted (Figure 11a–d). Flood intensity values that are higher than 2.10 m²/s and water depths that are greater than 1 m can endanger human lives and cause extended damage.

Table 4. Areas belonging to the different categories of the ASCE hazard criterion for the various LULC types in the affected area.

	DEM			DSM
ASCE	<2.1	2.1~3	>3	<2.1	2.1~3	>3
Overtopping	30.07	6.64	140.09	39.79	8.82	136.86
Settlements	7.59	1.89	22.09	7.08	1.99	19.83
Artificial, Bare soils and rocks	1.3	0.48	6.92	1.62	0.78	6.95
Olives	8.08	1.71	57.24	13.11	2.54	57.35
Greenhouses	9.55	1.68	29.76	11.22	2.07	28.42
Forest-Shrubs	4.18	0.88	24.48	6.76	1.44	24.31
Piping	33.87	5.96	129.38	42.6	8.46	121.47
Settlements	9.02	1.73	19.68	9.42	1.74	15.57
Artificial, Bare soils and rocks	1.68	0.39	6.16	1.39	0.44	6.54
Olives	9.57	1.88	52.89	13.56	2.57	52.98
Greenhouses	9.52	1.12	27.68	11.92	2.56	23.98
Forest-Shrubs	4.08	0.84	22.97	6.31	1.15	22.40

shows the percentages of the affected LULC areas, according to the ASCE criterion, when applying the dam breach scenarios and the elevation models. More than 78% of the affected areas belong to high and very-high hazard classes. From these, approximately 32% are in regard to olive groves, and 11% are in regard to the threatened area of the village of Gra-Lygia.

6. Discussion

From the comprehensive analysis of the results, it becomes noticeable that the overtopping scenario is more adverse than the piping one, both in terms of the peak and the time of breach event completion.

Moreover, the analysis of the affected areas reveals that the predicted flooded area calculations of the DEM and the DSM, in the case of both dam breach scenarios, are profoundly different. The topographic and surface characteristics of an area, particularly for local-scale events, have a significant role in water flow and hydraulic modelling [100]. This is easily understood if we consider that the theoretical routing of water between large volumes such as buildings or greenhouses lead it to longer distances. Since these volumes are the actual current state of the surface downstream of the dam, it is important to consider them as significant contributing data in future studies. Thus, the extraction of very precise surface characteristics of an area using UAS or other sources such as laser scanning (light detection and ranging (LiDAR)) may improve the prediction of flood behaviour considerably [50,101]. However, during a flood of this magnitude, part of the anthropogenic elements, such as greenhouses and the natural vegetation, represented in the DSM may be destroyed and drift away. Accordingly, the optimum representation of topographic characteristics in the flow domain is challenging and requires further study on which features should be retained depending on flow depth and velocity.

The use of the DSM instead of the DEM increases the depth values in both overtopping and piping scenarios, probably due to the more detailed representation of the coarse surface by the very detailed model, with much higher spatial resolution, which is also able to reveal the spaces between the high buildings and the greenhouses and the “urban gorges” from where the water channelling leads to greater water depths. On the other hand, the velocity values for both cases (overtopping and piping) using the DSM increases inside the riverbed, as the detected width (by the more precise DSM) is narrower. However, in this case, the velocity is significantly lower in the plain area because of the improved assessment of the water flow obstacles, such as the existing infrastructure and greenhouses. The water depth and velocity show higher values in the overtopping scenario than the piping scenario when using either the DEM or the DSM. The water arrival time decreases considerably with the use of the DSM because of the model characteristics that provide extremely detailed and high accuracy surface relief, helping in more precise calculations. The arrival time for both scenarios (overtopping, piping) looks equal when the DSM is utilised (Figure 10).

Furthermore, the flooded areas are almost equal in all scenarios, with the overtopping scenario corresponding to a slightly larger area than the piping scenario in the case of the DEM. In the case of the DSM, the flooded areas in both scenarios increase, due almost certainly to the many physical (e.g., vegetation) and anthropogenic (e.g., greenhouses, buildings) barriers involved in the relief of the DSM, which force the water to pass through them, reaching and covering further and larger areas. Along the coastline, the flood reaches a width of about 2 km in all cases, which is probably dictated by the topography of the affected area.

Although the very high spatial resolution of UAS data contributes significantly to the effectiveness of such studies, they create big data that may require very demanding processing, revealing challenges and opportunities for several aspects of state-of-the-art parallel processing systems [102]. However, the very high spatial resolution of less than half a metre offers the opportunity to include in the studies small volumes (rocks, bushes, various piles), which, as a whole, may affect the scenario of the routing significantly.

It must be considered that in order to achieve more accurate results in the future, the model should also include other important DSM data parameters that directly affect flood behaviour, such as the resistance and durability of structures, the existence of vegetation and other obstacles that the water encounters during its passage, water permeability, the roughness of artificial surfaces, the type of vegetation, and objects that drift easily (such as vehicles, small trees or light structures), which may cause mass downhill accumulation.

Remote sensing data are effective tools for detailed LULC classification and hazard-related applications like flooding, which require precise mapping [103,104]. The open-access and continuous data of ESA’s Copernicus program of Sentinel satellites offer the possibility for comprehensive mapping of Earth’s surface [105,106] and the systematic exploitation of more operational frameworks required for the monitoring of floods [107,108].

Finally, the hazard analysis of depth and volume parameters of water waves has revealed that, in our case, a very big percentage of the downstream residential area and the human infrastructure and activities (greenhouses) belong to the very-high risk category of the ASCE hazard criterion, which increases the potential impact on both human losses and severe economic consequences dramatically. Motivated by this, it is obvious that effective and comprehensive decisions should be taken in the future, starting from understanding the extent of potential hazards and risks, to identify and evaluate every similar vulnerable area by analyzing the local conditions of each case that could pose a potential threat to people, properties, critical services and economic activities [109]. Therefore, the usefulness of the implementation of such studies and their derived results is related to the fact that these studies can be used as a guide for other similar areas that could undergo comparable implications from analogous catastrophic phenomena.

7. Conclusions

In the present study, the dam break simulations in both overtopping and piping scenarios, using HEC-RAS software to compare the differences of using elevation (DEM) and surface (DSM) models, were assessed. Subsequently, a hazard analysis was performed to determine potential threats in the vulnerable area downstream of the Bramianos dam in Crete.

Breach dimensions, peak discharges and flood wave characteristics (such as variation of flow velocity and water depth) that were calculated based on the hydraulic model for the examined scenarios were compared. The results indicate the substantial effect of different break types on the hydraulic characteristics of a dam-break flood wave in the downstream areas, with the overtopping scenario being the more adverse event. The outcomes indicate that the extent of water flow downstream of the collapsed dam depends largely on the water volume stored in the reservoir. Moreover, the comparison of DSM and DEM cases shows that the DSM represents the surface relief and the existing natural obstacles more accurately, providing a more precise depiction of the existing barriers, such as vegetation, buildings, and greenhouses, and, hence, more realistic hydraulic simulation results. However, further study on the optimum representation of topographic characteristics of the flow domain, which depend on flow depth and velocity, is required.

As a concluding remark, hazard resistance to natural disasters needs to be further reinforced because the value of human life and assets are of great significance. The evaluation of dam failure occurrence probability and the possible effects, using dam-breach flood simulations and innovative remote sensing data, can deliver better knowledge and lead to adequate preparation to confront such flood incidences. The protection of society from a given type of hazard is of primary importance, and being able to assess the spatial distribution of the threatening phenomena and their intensities is critical to this endeavour.