Innovative Risk Assessment Framework for Hydraulic Control of Irrigation Reservoirs´ Breaching

This research introduces an innovative framework aimed at developing a risk assessment to analyse the breaching hydraulic control of non-impounding reservoirs for irrigation purposes, called irrigation reservoirs (IRs). This approach comprises an analytical method based on several empirical formulas where the one that best fits the different geometric characteristics of IR water systems is chosen. Furthermore, a stochastic framework allows for the incorporation of the occurrence probability as a tool to characterize the risk analysis of IRs. This occurrence probability has two components: probability based on the bottom elevation of a final breach and probability based on the failure mode (piping in this case). In risk assessment terms, the ultimate product comprises the maximum hazard probability maps that allow a significant improvement in the representation of the artificial flooding effect. This research was successfully applied in two dimensions, synthetically and realistically, in the Las Porteras and Macías Picavea IR water systems (Spain). This approach may improve the management of this type of hydraulic infrastructure and its surrounding area by reducing the risk of experiencing negative consequences derived from uncontrolled hydraulic breaching.


Introduction
In recent years, flood risk management has been the subject of numerous regulatory developments at many scales. These developments have obliged owners to ensure the safety of their hydraulic infrastructure. Within the scope of Spain's legislation, Royal Decree 9/2008 introduced, for the first time, a new section concerning dam, reservoir, and irrigation reservoir (IR) safety. This fact, together with the publication of Royal Decree 903/2010 on the evaluation and management of flood risk, Directive 2007/60/EC defines a new management framework in which the new infrastructure must be readapted considering its potential danger.
IRs need to be differentiated from dams largely for the following reasons. An IR's location is one of the most important differences from dams because it limits the stage of solicitation, the failure mode and the potentiality of the damage caused. This fact offers notable uncertainty because a priori, the failure zone is unknown. On the other hand, the advantage of an off-river IR location means that the failure hydrograph may be enclosed in the first moments of collapse. Moreover, in the case of IRs, the event associated with the predominant failure mode is clearly localized. The consequences caused by the outflow hydrograph are limited to a smaller spatial context than a natural hydrological process.
The number of IRs in Spain is close to 65,000, mainly concentrated in the Mediterranean area, while the number of large dams in Spain is close to 1,400. Owners of IRs, mainly irrigation communities, complain about their security obligations, which are the same as those of a large dam. These obligations rely on certain requirements, which are based mainly on the movable water height (above 5.00 m) and its corresponding volume (above 100,000 m 3 ). Irrigation communities make claims related to the normative demands inherited from the dam world based on the obvious differences between IRs and dams (Ouyang et al. 2018).
Moreover, specific studies on IR risk assessment are scarce in the scientific literature. There is practically no specific bibliography of IRs. In the case of tailing IRs, the safety analysis has improved, but the effects are not comparable to those of water IRs. Fortunately, it is possible to extrapolate some behaviours from embankment dams to IRs. Embankment dams hold an important position in the scientific literature, especially in safety. A considerable amount of literature has been published on the characterization of embankment dam breaches (Froehlich 1987(Froehlich , 2008MacDonald and Langridge-Monopolis 1984;Von Thun and Gillette 1990). Furthermore, the failure mode of earth dams has been widely studied, especially overtopping (Coleman et al. 2002) and internal erosion (Fell et al. 2003). Additionally, in recent years, different organizations have been working on specific analyses of IRs to improve their safety, such as International Commission on Large Dams (2016) and Ministerio de Agricultura y Pesca, Alimentación y Medio Ambiente (2012). Moreover, the incipient popularization of IRs in Spain is causing the necessary development of scientific studies about their peculiarities. These studies cover several fields, such as the calculation schema used (Soler et al. 2012) or questioning the safety system of IRs (Sánchez 2014;Espejo 2017).
IR water system characterization is firmly associated with the geomorphological variability of the environment. The roughness and terrain slope around the IRs are the most important parameters. The estimation of geomorphological parameters is conditioned by the resolution and accuracy of the digital elevation model (DEM), which, in turn, conditions the routing of the outflow hydrograph produced by an IR failure. The terrain conditions all hydrological processes, and laser imaging detection and ranging (LIDAR) or aerial photogrammetry are perfect tools to cover its representation (Molina et al. 2014). Therefore, the IR water system restricts mitigation measures and the behaviour of the failure hydrograph.
Actions to palliate flooding due to IR breaching can be classified into structural and nonstructural (Gendreau and Gilard 1997). Schanze et al. (2008) proposed ordering mitigation measures according to three fields: construction (structural vs. nonstructural), effect (reduction hazard vs. vulnerability reduction) and functionality (instruments vs. measures). The purpose of this study is to act within the second aspect (effect) on the reduction of dangerousness. Identification of the element of the system that acts under the prism of risk reduction is an essential factor in the modification of the hydrological-hydraulic behaviour of the affected basin drainage (Salazar 2013).
Given the serious scarcity of general risk approaches specifically focused on IRs, this paper mainly aims to establish an innovative risk assessment framework for the hydraulic control of breaching IRs. This will be largely useful to predict the hydraulic operation of IR breaching through numerical simulation (Bladé et al. 2014) in terms of the probability of the occurrence of a certain maximum flow.

Proposed Methodology
This methodology comprises 3 consecutive phases: preprocess, process, and postprocess. The preprocess phase has 3 substages: study of the breach geometric parameter estimation, sectorization and stochastic analysis. The process phase contains the hydraulic simulation substage. Finally, the postprocess phase has 2 substages: maximum hazard probability mapping and mitigation measure assessment for the failure hydrograph ( Fig. 1).

Breach Geometric Parameter Estimation
It is necessary to reverse the current situation and analyse the existing records that allow the characterization of geometric IRs. To delimit the areas of greatest concentration in the legislation, a representative sample of 231 IRs was used for this purpose. This analysis considers the movable water height and its corresponding volume. In this way, the geometric characteristics of IRs are established through the study of the hydraulic response caused by the breach.
Currently, there are three methods of IRs breach analysis, which differ in the physical basis on which they are based: Empirical Models (EM) (Froehlich 2008), Parametric and Analytical Models (PAM) (Al-Rbeawi 2018), and Physical Base Models (Espejo 2017).
First, it is unfeasible to conduct the analysis using PBMs due to the economic cost involved, the required facilities and the team of professionals needed to build the model. On the other hand, performing the analysis with PAMs would allow for the reproduction of the behaviour of the IR embankment. Therefore, the conditional probability of breaking could be estimated given a load scenario. With a complete geotechnical analysis, it would be possible to define a model that represents the response of the IR breaking system. From this model, it would be feasible to extract the probability density function (PDF) of IRs breaking through the geotechnical parameter of an IR embankment. This process is not addressed by the owners of IRs due to the high cost and complexity of geotechnical tests. This fact led us to discard the analysis using PAMs.
The EM, based on numerical simulation, was chosen as the best method of analysis to be adopted for this research. This type of analysis provides the best solution in terms of time and cost, as opposed to the other types of analysis. Furthermore, the results obtained may be extrapolated to a multitude of IRs with the same characteristics.
Therefore, based on the generated record and the EM analytical method, selecting the empirical formula that best fits the IR characteristics and behaviour is possible. This selection is based on an adapted methodological proposal on dam breaches by the State of Colorado (2010) and applied to this research IR sample. In accordance with normative requirements, three application fields are distinguished according to the analysed dam magnitude: minor dams, small dams and large dams (State of Colorado 2010). The aim of these application fields is to recommend an empirical formula that adequately responds to the predictable behaviour of this type of dam. In this way, this recommendation serves as a reference for managers and users of this water infrastructure.
The definitive selection of the formula was associated with a characterization factor called storage intensity (SI, in cubic metres per metre, m 3 /m), which relates the reservoir volume stored (Vw, expressed in cubic metres, m 3 ) with the water depth above the breach (Hw, expressed in metres) (SI = Vw/Hw). This factor, introduced by the State of Colorado (2010), delimits the application fields of some empirical formulas. This formula delimitation is based on the sensitivity to SI variation.
For this purpose, the variables evaluated were breach characterization (final breach width and breach formation time) and failure hydrograph. These variables most influence a hypothetical IR's structural failure (Ahmadisharaf et al. 2016) in terms of flooding. The previous inputs, together with the evolution of these flows through the natural drainage network, define the main study items. This research proposes an IR breaching characterization based on the hydraulic response in the form of a hydrograph generated by structural failure. Among the multitude of existing empirical formulas, those that are best adapted to the existing IR rupture records in a territory were chosen for this analysis. The selected empirical formulas comprise MacDonald-Landgridge-Monopolis's expression (MLM84; MacDonald and Langridge-Monopolis 1984), Von Thun-Gillette's expression (VTG; Von Thun and Gillette 1990), Froehlich's expression (F95;Froehlich 1995), the expression of a Spanish technical guide for classification of dams according to risk potential expression (GUIDE; Ministerio de Agricultura y Pesca, Alimentación y Medio Ambiente 1996), which is a variant of the proposed expression by Froehlich (1987) and Froehlich's expression, F08, (Froehlich 2008).
A study of the causes of earth dam failure carried out by several authors (Aboelata and Bowles 2008;Sills et al. 2008) concluded that the main mode of failure is quality problems in 42% of cases, of which 60% are related to siphoning. Therefore, the most widespread failure mode in this type of solution, "piping" or internal erosion, was chosen (Foster et al. 2000a, b). This study considered siphoning as the only mode of failure. Therefore, to characterize the IRs according to their rupture behaviour, different synthetic types of IRs were constructed, with different movable water heights and various dimensions, covering the geometric variability susceptible to normative obligations.
Therefore, it is possible to construct a numerical model to analyse the hydraulic response of IRs in a wide range of cases. The studied empirical formulas facilitate the geometric characteristics of the breaches for each case. Thus, the equation that best fits the IR characteristics, based on the movable water height and its corresponding volume, was selected. Then, it was possible to select the correct empirical formula for each case of study, from which the breach parameters that feed a model of hydraulic simulation of IR failure were obtained. Therefore, it was possible to develop a breach model to make a comparative analysis of the results for each synthetic IR. This premise was achieved by studying the sensitivity of these results to the variation in movable water height, water volume and associated water load. Obtaining and analysing the failure hydrograph plays a key role in this methodological approach, since according to the selected empirical formula for each geometric type of IR, the resulting failure hydrograph will be different.

Sectorization
One of the greatest uncertainties in the application of IR risk analysis is the location of the breach geometric axis in the numerical simulation. There is no preferential section for the location of the breach geometric axis (Chauhan and Bowles 2004).
The main aim of this methodological phase is to perform the sectorization of the area surrounding the IRs. To achieve this objective, it is necessary to analyse the IR water system. Transfer to reality demands the characterization of the geomorphological environment of the affected basin. This approach aims to ensure a complete analysis of the entire closure dam by sectoring it in its entirety and to propose a single representative type of breach for each sector. This avoids the arbitrariness in its location and ensures that the behaviour of the entire closure dam is analysed. This sectoring requires a thorough knowledge of the surface on which the discharge flows must evolve. Among the possibilities for analysis, geomorphometry has been used, particularly focused on the first vulnerable zone. Technically, this is carried out by delineating the drainage basins.
Certainly, geomorphological characterization is considered a valid and objective tool for identifying the areas susceptible to channelling the failure hydrograph. The morphometric analysis of the digital model studies the relief and orography of the surface; on the other hand, it characterizes it through additional parameters, which are described in Olaya and Conrad (2009). This type of study is based on a focal analysis, which studies the cell values and their environment to obtain a definitive value. To characterize the area surrounding the IRs, it is necessary to define the main geomorphological parameters.
The geometric axis breach is based on the Horton-Strahler stream order (Horton 1945;Strahler 1957), mainly due to its ease of implementation. Therefore, the threshold that defines the sector's channel was established according to the number of cells that feed it. A useful indicator for choosing the representative Horton-Strahler order number is the topographic wetness index (TWI). The TWI, defined by Sörensen et al. (2006), indicates the susceptible zones of accumulating wetness, which helps to explain the involved hydrological processes. This fact helps to determine a sufficiently representative stream order number to avoid the excessive fragmentation of the border zone sectors of the IRs, which would impair the feasible implementation of the methodology.
Later, the location of the geometric axis of IR rupture was calculated. This was performed by delineating an axis on the intersection of the outer slope of the IR with natural terrain. Once this axis was traced, it was possible to obtain the minimum absolute elevation for each sector. Thus, the bottom elevation of the final breach was defined for use in stochastic analysis.

Stochastic Analysis
An innovative approach for a better and more accurate cartographic representation of the area most likely to be dangerous is proposed. The definition of current regulation associated with different return periods for flooding studies may be improved considering this approach. These regulations concerned natural floods, and artificial floods caused by water infrastructure failure were not considered.
This stochastic framework requires the implementation of the Monte Carlo method (Díaz-Emparanza 2002), which requires previous knowledge of the PDF of the object variable (Tung 2011). Due to the lack of data, to characterize this experiment, it was necessary to use a triangular PDF proposed in the scientific literature (Johnson 1997). The determination of triangular PDF parameters, including minimum, maximum and mode, were obtained according to the indications by Ahmadisharaf et al. (2016). This employed a triangular PDF and was based on the bottom elevation of the final breach (Ahmadisharaf et al. 2016).
Hydrologic Engineering Center and Hydrologic Modeling System (HEC-HMS) software, Version 4.2.1, was used to apply the Monte Carlo method. The HEC-HMS model was developed by the US Army Corps of Engineers (US Army Corps of Engineers 2000). In this method, 5,000 simulations were carried out (Díaz-Emparanza 2002) per type defined by its height and for each of the 4 empirical formulas mentioned above. On the one hand, this allowed for the validation of the triangular PDF; on the other hand, it proposed another adjustment to the distribution for the generated flows.
Therefore, there were four sequential objectives in this phase: (1) determining the occurrence probability of a breach considering its final bottom elevation through a triangular PDF; (2) using the parameters obtained in the selected empirical formula of the first phase to simulate IR breaches with the Monte Carlo method; (3) applying the results of the Monte Carlo method in each sector to determine its geomorphometric probability; and (4) analysing the results in terms of maximum generated flows and adjusting a PDF to those flows in the synthetic IRs.

Process: Hydraulic Simulation
To develop probability maps and to propose structural measures that laminate the generated hydrograph, it was necessary to develop a 2D hydraulic model associated with each breach and considered scenario. This hydraulic model was substantiated on the resolution of the 2D shallow water equations (also known as the two-dimensional Saint-Venant equations; Tayfur et al. 1993). In this phase, both synthetic IRs and real cases were analysed, which allowed for the validation of the laws of hydraulic behaviour. The implementation of the model was carried out in free software called IBER (Bladé et al. 2014).
The selection of the DTM considers parameters such as resolution, accuracy, and discretization of the ground system in the hydraulic simulation software. In this case, the LIDAR data of the 2010 Spanish National Plan of Aerial Orthophotography (PNOA) flight were used for the elaboration of the DTM. These data had an initial density of 0.5 points/m 2 and an altimetric accuracy < 0.40 m (Equipo Técnico Nacional 2001). Consequently, a DTM obtained from the PNOA LIDAR data source was considered a reliable basis for the scope of this study. The simulation duration was 7,200 s.
The geometry model created in IBER had a higher resolution than similar studies of dam breaks. The triangles of the numerical mesh used for hydraulic calculation had the following characteristics: 1.0 m minimum side, 50.0 m maximum side and a tolerance (maximum vertical distance between the DTM and the generated geometry) of 0.2 m.
IBER allowed exporting the hydraulic results in the Arc Info ASCII grid format. Therefore, for the postprocess stage, it was useful to obtain maximum maps of both hydraulic variables and hazards.

Post Process
The main aim of this stage was to obtain the maximum hazard probability maps as a reliable tool for the establishment of forthcoming hazard attenuation measures. Under this approach, the results offered by the synthetic IRs were validated regarding the results provided by real IRs. There are several methodological approaches available to validate the results, which can verify and calibrate the parameters of a mathematical model based on real data. Other proposals allow for the evaluation of similar behaviour between different models. In this respect, it was appropriate to validate the data by comparing the results of the synthetic IRs with those from real IRs.

Maximum Hazard Probability Mapping
The results obtained from previous stages made it possible to create a new cartography. This product would improve the management of flood risks caused by the existence of IRs.
Regulations regarding flood risk management proposed sectoring of the territory based on natural floods in return periods of 4, 100 and 500 years. This conception was not adapted for flooding caused by the failure of hydraulic infrastructure. However, this was considered by the scientific community (Di Baldassarre et al. 2010). Proposing a new method for generating risk maps due to the failure of hydraulic infrastructure is considered necessary.
Consequently, a probabilistic characterization of the results for each sector allowed us to produce maximum hazard probability maps. Probabilistic characterization responded to the occurrence probability based on the bottom elevation of the final breach and failure mode.
To achieve the proposed objectives, the following workflow was proposed: 1. Characterization of each sector by the probabilities of geomorphological analysis. 2. Characterization of each sector by the probabilities of failure mode (piping), according to the scientific literature (Foster et al. 2000b). 3. Obtaining the overall probability of each sector as a result of the previous partial probabilities. 4. Transformation of this probability for its total computation with respect to 1 for each sector. 5. Reclassification of maximum hazard rasters obtained from hydraulic simulation. 6. Overlap of reclassified rasters of all sectors. 7. Obtaining the sum of all previous rasters by map algebra. 8. Reclassification of the previous results according to the following criteria: low probability (< 33%), medium probability (33% < P < 67%) and high probability (> 67%). 9. Creation of a maximum hazard probability map of the desired threshold.

Synthetic Implementation
Ten different synthetic IRs with 3 different movable water heights (5.00, 7.50 and 10.00 m) were designed. The slopes of the walls inside and outside the IRs were the same, 2.50H/1 V. The smallest IRs had a quadrangular base of 25.0 × 25.0 m, which gradually increased until 300.0 × 300.0 m in the largest IRs. The synthetic IR types covered all safety regulations for the scope and multitude of SIs. Therefore, the sample IR geometry was characterized by a wide range of synthetic IRs.

Real Implementation
The main aim of this analysis was to obtain the cartography of maximum hazard probabilities. For the development of case studies, sectorization, breach parameters drawn from the empirical formula from the first stage and stochastic analysis were implemented.
To implement this methodology, two types of water systems were chosen that differed in geomorphological characteristics. In this sense, the case studies included the Macías Picavea IRs (Valladolid), with a high altimetric gradient, and the Las Porteras IRs (Ávila), with a low altimetric gradient (Fig. 2).
Macías Picavea IRs had an earthfill wall with a height of 6.75 m, reservoir capacity of 201,623.00 m3 (normal pool level) and an IR crest of 5.0 m. The second studied case, Las Porteras IRs, had a wall height of 7.50 m, reservoir capacity of 397,261.00 m 3 (normal pool level) and crest of 4.5 m. Both case studies were located in the Castilla and Leon region, Spain, in the Duero River Basin. In this sense, hypothetical piping failure was used for both case studies.

Estimation of Breach Geometric Parameters
The sample of 231 IRs was considered representative of the average geometric reality of existing IRs. This record revealed that in the study area, there was a greater density of IRs with approximately a 5 m water height and 20,000 m 3 reservoir volume (minor dams) in the classification of the state of Colorado (2010); however, there was also great variability. Given this context, the empirical formula best suited to changes in movable water height and reservoir volume should be chosen.
The breach progression method and breach development estimators were used to validate the obtained results. The erosion rate (ER, in metres per second, m/s) can be used as a check of the set methods and parameters (State of Colorado 2010). A linear progression was considered valid when 1.6 < ER/Hw < 21.0, where ER = Bavg/tf. Whole breach development was considered valid when Bavg/Hb > 0.60, where Hb is the height, in metres, of the breach and Bavg is the average width of the final breach expressed in metres (m).
The analysis of the results revealed that the maximum and minimum flow rates for the same IRs greatly differed depending on the empirical formula used. In medium-sized IRs (250,000 m 3 ), the flow rate doubled. This is undoubtedly an indicator of the level of uncertainty associated with this analysis and the need for its continued revision.
The VTG formula showed logarithmic behaviour for a movable water height below 5 m. Above this height, the growth was parabolic, softening in the central zone. The comparison between empirical formulas proposed by several authors revealed that the expression presented in Ministerio de Agricultura y Pesca, Alimentación y Medio Ambiente (1996) overestimated flow rates. GUIDE overestimated the flow rate by 35-50% with respect to MLM84 and 20% with respect to F08. The results confirmed that MLM84 is more versatile for small IRs (minor dams). F08 showed the best fit and sensitivity to SI variation.
After applying new validation schemes in terms of formation time, Bowles et al. (2014) considered reasonable formation times associated with the expression tf = 2 Vw/

Fig. 2 Location of case studies
Qp, where Qp is the peak breach flow rate in cubic metres per second (m 3 /s) and tf is the breach formation time in seconds (s). This made it possible to certify that the expressions chosen for each application field were appropriate, such as those with the lowest outflow peak, longer formation time and greater lamination of the hydrograph. Figure 3 shows the empirical formula that best meets these conditions for each IR's normative requirements. A second empirical formula (2 nd choice) to corroborate the initial premise was proposed.

Sectorization
Generated watershed algorithms are based on a drainage network definition and on the number of input cells chosen. In this analysis, the following software was used: SAGA GIS v.5.0 for its level of geomorphometric analysis and Global Mapper v.18.2 for its versatility.
The breach geometric axis location was based on Stream Order 4. Therefore, the threshold was determined according to the number of cells that feed the drainage channel. Uncertainty arises when deciding on a sufficiently representative stream order to avoid an excessive number of basins surrounding the IRs. As mentioned above, the TWI was a great help in recognizing the drainage network in the DTM and avoiding the mentioned uncertainty.
For the Las Porteras IR case study, different cell thresholds were studied for the same stream order: 123,456, 12,345 and 1,234. Once the IR sectors were known, it was necessary to identify the minimum height of each sector (on the intersection of the outer IR slope with natural terrain). There were 4 sectors for Las Porteras IR and 9 sectors for Macias Picavea IR. The breaches' geometric axes should be located on the position of minimum elevations for each sector, and afterwards, stochastic analysis should be applied. Fig. 3 Application field representation of a suitable empirical method for each IR normative requirement

Stochastic Analysis
For each IR, each height and empirical formula, the Monte Carlo framework obtained a distribution of the hydrograph attributes. Additionally, with the implementation of the Monte Carlo method, the occurrence probability was achieved for each breach of case studies. These parameters corresponded to the bottom elevation of the final breach predicted by stochastic analysis, following the triangular PDF.
Once the method was applied, the average, maximum and minimum hydrographs together with the standard deviation were obtained for each case. This allowed for a prediction of the range of results in each case. Then, with the hydrographs obtained and to characterize the breach, the breach bottom elevation results were adjusted to new PDFs. The adjustment was performed again to a triangular PDF validated by 3 methods: Kolmogorov-Smirnov (Massey 1951), Anderson-Darling (Anderson and Darling 1954) and Chi-square (Lancaster 1969). This adjustment was found to be correct for the reality of the IRs. The base elevation of the breach showed a failure rate with exponential growth when approaching the IR wall crest, which denoted randomness in the failure process.
To finally analyse the results obtained in terms of peak outflow generated, for each synthetic IR and height, box plots were constructed. Jointly, mean confidence intervals were also determined. Then, selecting from among more than 50 probability distributions, the one that offered the best fit in the Q-Q plot was developed. Therefore, resulting peak outflows of 60% of the analysed IRs could be adjusted to a Wakeby distribution, which was proposed by Houghton (1978). The Wakeby distribution is a suitable probabilistic model to represent maximum daily or instantaneous flood flows in a region. A smaller percentage of obtained flow rates was fitted to other probability distributions, such as Johnson SB (Campos-Aranda 2015) or Beta-Dagum (Domma and Condino 2017).
Three different behaviours were detected in IRs of similar size for the analysis of generated flow rates. The three smallest IRs had a slight asymmetry of simulated peak flows, which was more pronounced than intermediate and large IRs. This asymmetry was due to higher extreme values than in other cases. Data dispersion dissipated as the IR height increased, as shown in Fig. 4. The three largest IRs had the same tendency of data and thus the same frequency distribution. To conclude, as the size of the synthetic IRs increased, the bias in the peak outflow data gradually decreased (Fig. 4).

Maximum Hazard Probability Maps
To produce the final mapping of the maximum hazard probabilities, the probability data from the stochastic analysis phase were necessary. The occurrence probabilities based on the bottom elevation of the final breach (geomorphometric) were multiplied by the probabilities of the failure mode (causal), and the resulting probability was computed with respect to 1. The corresponding probability was assigned to each sector obtained from the sectorization phase. Finally, three hazard groups were distributed to reclassify the probability raster.
Since the characteristics of the IR water systems of the case studies were dissimilar, the maximum hazard behaviour was different in the two case studies, as spatially represented in Fig. 5. In the Las Porteras Irrigations reservoir case, most maximum hazards were in the IRs close to the infrastructure. In contrast, in the Macias Picavea Irrigation reservoir case, the maximum hazard zone was located in the area where the flow had stabilized. This was F. Espejo et al. due to the steep terrain. This fact demonstrates the importance of the geomorphology of the environment in this type of study. Moreover, the flat areas around the Las Porteras Irrigation Reservoir partially retained the flood generated by IR failure. It would therefore be interesting to apply mitigation measures at the Las Porteras Irrigation Reservoir.

Discussion
An innovative risk assessment framework for the hydraulic control of IR breaching was implemented in this study. Four empirical breach prediction methods, including F08, VTG, MLM84 and GUIDE, were compared to obtain the empirical formula that best adapted to each geometry IR in a territory. A stochastic analysis was used to obtain the occurrence probabilities related to the bottom elevation of the final breach. Sectorization and hydraulic simulation are very useful tools to produce maximum hazard probability maps. The study of the breach geometric parameter estimation was effective. This study characterized the hydraulic operation of IRs and defined the empirical formula that best responded to the geometric reality of each case. However, the restrictions of the method used should always be considered: it only evaluated a failure mode (internal erosion), and the embankment cohesion was considered null.
One of the issues in this study was the lack of IR records. This problem was solved with a representative sample of 231 IRs and the development of synthetic IRs covering the geometric characteristics of the IR water system. This approach was validated by verifying that the hydraulic operation of real case studies (Las Porteras and Macias Picavea IRs) closely approximated the studied synthetic IRs.
Moreover, the infeasibility of tackling an exhaustive geotechnical analysis and building a PAM led to the implementation of an EM. This decision was based on the ease of EM implementation, access to required data and viable extrapolation of results. This EM relies exclusively on geometric parameters.
The comparison between empirical formulas proposed by several authors revealed that the formula presented in some relevant reports overestimates the flow rates. The SI parameter and movable water height were proven to be robust parameters for characterizing the potential risk of the hydraulic infrastructure and delimiting the scope of each formula.

Conclusions
The location of the breach geometric axis of the hypothetical failure was objectively justified after applying the proposed methodology. The IR embankment was divided according to exclusively geomorphological factors. Furthermore, although IR failure is an artificial event, the importance of receiving basin characteristics was revealed. Therefore, the receiving basin conditioned the failure hydrograph; therefore, it was necessary to include its study in the risk analysis of IRs.
It is important to note that mitigation measures will be feasible if geomorphological conditions in the IR environment allow it. The proposed solution should not be implemented in IRs with high-gradient bordering areas.
The stochastic analysis, not addressed thus far in IRs, showed its strengths and weaknesses in the analysis. Implementing the Monte Carlo model allowed for the determination of the breach occurrence probability independent of the probability associated with the failure mode. The probability achieved after applying the Monte Carlo method was based exclusively on geometric factors, since it was associated with the bottom elevation of the final breach. The main weakness in applying this framework of analysis was focused on the utilized triangular PDF. Thus, the option of proposing new PDFs supported by a broad range of case studies remains open for upcoming studies. Furthermore, the Wakeby distribution proved to be the best fit for peak outflow data obtained from stochastic analysis.
Hydraulic simulation was the tool used to produce maximum hazard probability maps using parameters obtained from previous phases. One of the main conclusions of this stage was the sensitivity of the hydraulic results when working with different discretizations of the DTM. This is one of the research lines that remains open and worthy of study.
Maximum hazard probability maps are a new contribution to flood risk managers and land-use planners. The causal occurrence probability and probability associated with the bottom elevation of breaches were used. The proposal presented in this paper improves the knowledge of the territory's response before the existence of this infrastructure with failure risk.
The problem due to the lack of a record of IRs and failed IRs was expounded during the study. Therefore, it would be of interest to systematically address IR registration throughout the territory. This knowledge would make it possible to characterize the IRs reliably, validating the results with geometric methods. This effort would result in better risk management. Thus, the development of a centralized and open access database would allow us to advance the knowledge of failure modes, being able to limit the uncertainties associated with failure mode.
It was stated that the analysis mode used in this study was through the empirical model. The convergence of the three existing analysis modes (EM, PAM and PBM) would be a very interesting development. This development would be useful to determine the laws of behaviour in IR failure.
Acting on the nature of the cartographic source is considered necessary. Thus, LIDAR cartography was used in this work. The next step will be to use other cartographic sources, resolutions, and accuracies to reliably establish the minimum cartographic characteristics needed.
Geometric parameters were considered for the study of the potential risk of IR collapse. There are other uncertain parameters that can be included in the analysis. Therefore, the probabilistic model has considerable room for improvement.
The goal of this whole methodology is the development of a decision support system (DSS) that considers hydraulic and geometric parameters and determines the risk infrastructure for possible hazards. Therefore, the next step of this research is automating the proposal methodology through software development considering free distribution tools.
Author Contributions F.E. and J.L.M. conceived, designed and led the research and paper editing; S.Z. supported the paper editing and the research conceptualization; and R.M. was in charge of the analytical development. The Discussion and Conclusions sections were addressed by all authors, and all authors wrote the paper.
Funding Open Access funding provided thanks to the CRUE-CSIC agreement with Springer Nature. The authors did not receive support from any organization for the submitted work.
Availability of data and materials The datasets generated during and/or analysed during the current study are available from the corresponding author on reasonable request.