Analysis of the Effects of Reservoir Operating Scenarios on Downstream Flood Damage Risk Using an Integrated Monte Carlo Modelling Approach

Abstract

The aim of this study is to analyse the effects of reservoir operating scenarios, for flood damage evaluation downstream of a dam, using a Monte Carlo bivariate modelling chain. The proposed methodology involves a stochastic procedure to calculate flood hydrographs and the evaluation of the consequent flood inundation area by applying a 2D hydraulic model. These results are used to estimate the inundation risk and, as consequence, the relative damage evaluation under different water level conditions in an upstream reservoir. The modelling chain can be summarized as follows: single synthetic stochastic rainfall event generation by using a Monte Carlo procedure through a bivariate copulas analysis; synthetic bivariate stochastic inflow hydrograph derivation by using a conceptual fully distributed model starting from synthetic hyetographs above the derived; flood hydrographs routing through the reservoir taking in an account of the initial level in the reservoir; flood inundation mapping by applying a 2D hydraulic simulation and damage evaluation through the use of appropriate depth-damage curves. This allowed for the evaluation of the influence of initial water level on flood risk scenarios. The procedure was applied to the case study of the floodplain downstream from the Castello reservoir, within the Magazzolo river catchment, located in the southwestern part of Sicily (Italy).

1. Introduction

Floods are an environmental hazard that can cause heavy economic, environmental, and social losses. Therefore, their control is an important issue for damage mitigation worldwide. Floodwater storage construction facilities, such as reservoir dams, are some of the most common strategies of flood control through structural measures, which can also provide other benefits to people or local economies, including water supply during dry season, irrigation, recreation, and hydroelectric power [1,2,3].

However, large populations, infrastructures, and properties are located downstream of a dam, and in case of its failure or break, this can pose significant risks. When dams fail or malfunction, they can adversely affect people, their livelihoods, jobs, and businesses. Between 2000 and 2009, more than 200 dam breaks or failures occurred in many countries, causing disastrous effects on downstream areas [4].

Due to climate change and increased urban development, consequences of dam reservoir incorrect operational management have become much higher in the last decades [5,6]. Hence, a correct risk and damage assessment evaluation needs to consider all artificial structures holding water, which represent potential sources of flooding during operational flood control and in the case of possible failure of the structure.

Dam risk management and mitigation has become a high priority of organizations concerned with dams and valley safety, as well as with civil protection procedures. Most of the potential damages and losses occur along the downstream valleys. Past events show this evidence [7], and recent dam safety legislation includes some procedures related to the downstream effects of a dam failure. An effective mitigation of possible hazards, due to a dam accident or incident, clearly imposes an integrated risk management, including both the dam risk control and the valley protection. European Flood Directive 2007/60 requires potential damage evaluation to estimate the magnitude of the consequences of a flood.

In Italy, national laws and several technical guidelines have highlighted the need to set up specific operational dam reservoir management rules in case of downstream flooding risk by preparing appropriate reservoir flood control operation plans [8,9]. Currently, the evaluation and technical approval of projects related to dams is carried out by the National Dams Authority. The current Dam Regulation, relevant not only to the design and construction of new dams and to design of rehabilitation works of existing dams but also to the emergency action management in case of extreme events, is divided into two parts. The first one deals with the formal and administrative procedures, as well as general technical aspects to be followed. The second part is, instead, the “Technical Rules” defining all the technical details to be considered in the dam design and functionality. Emergency Action Plans (EAP) must be set up by local Civil Protection Authorities and coordinated by the Prefecture for various types of risks, especially for the downstream hydraulic risks of the dam. The Technical Rules are, in fact, aimed to ensure that, in case of extreme events, dams assure their operative functionality during the emergency phase. The Directive also establishes updated conditions to activate alert phases for dam safety and management of downstream hydraulic risks, defining the actions to be implemented in these phases.

Estimation of the potential damage, associated with the release of water downstream of a dam if an extreme event occurs, is a fundamental exercise to make decisions in ensuring safety of the properties placed downstream that, in case of flooding, must be reimbursed.

For calculating the potential damage of the properties for a forecasted or measured flood event, obtaining information about areas that would be inundated is a must and could be useful to the decision maker, in real time, to offer flexibility for adapting the water level in the dam to changing realities. Scientific literature emphasizes the importance of the reservoir operation optimization, finalized to define correct release decisions that guarantee not only the more common objectives of water management, such as the hydropower production, a reliable water supply, etc. [5,10,11,12], but the mitigation of downstream flood as well [13,14,15,16].

The inundated area in downstream areas depends on various factors, such as the height of the dam, its initial water level (IWL), the nature of failure, and the downstream vulnerability. The inundation maps show the areas that are likely to be submerged for different flood hydrographs as a function of the above factors.

When flooding occurs in a medium–big catchment, it is possible to carry out a two-dimensional analysis by running a model, in real time, in a reasonable amount of time to make appropriate decisions. However, in case of flash flooding, the time between the forecasting and the consequent flooding is not enough; hence, the a priori knowledge of the damage associated with a forecasted or measured flood event can be a valid approach.

Data availability is the main problem in flood hydrograph estimation. Even if discharge data are available, the length of the available hydrological series can be too short for any hydrological evaluation. In these cases, flood events may be estimated via rainfall-runoff simulation using observed precipitation data or generated data through stochastic simulation methods as the input. In this latter case, Monte Carlo simulation methods can be used to generate long series of rainfall events defined by their duration, volume, and temporal and spatial storm distribution [17,18,19,20,21].

Moreover, since flood peaks and corresponding flood volumes are variables of the same phenomenon, they should be correlated and, consequently, bivariate statistical analyses, and the theory of the Copula in particular, should be applied [19,20,21,22,23]. In literature, many hydrological applications (i.e., floods, storms, droughts), with bivariate flood frequency analysis by using copulas, are already presented and implemented in many case studies [20,21,22,23,24,25,26,27].

Dam hydrological safety assessment can, hence, be improved by bivariate flood frequency analysis [24,25,27,28]. Again, if long series of flood peak and volume data are available for a dam, the copula approach can be directly applied to the available samples; otherwise, an indirect hydrological simulation has to be carried out. Klein et al. [28], for example, coupled a stochastic rainfall generator and a continuous semi-distributed rainfall-runoff model to generate a synthetic long-term daily discharge time series.

Once the flood hydrograph series are available, a flood propagation model can be implemented to derive the correspondent downstream inundated area and the consequent flood risk scenarios. Damage analysis for different flood risk scenarios is useful to reveal how the reservoir rules operations (in terms of initial reservoir conditions) can, directly, reduce the negative effects of flooding in the downstream floodplain.

Several studies have assessed flood damages at different scales: from local to regional or macro area scales. Such assessments, however, are often limited in evaluating the flood impacts due to the absence of a global database of flood damage functions to translate flood water levels into direct economic damage. Usually, direct flood damage is assessed by using depth–damage curves, which represent the relationship between the flood damage and specific water depths for each asset or for each land use class [29].

In this paper, a new procedure is presented. The potential damage downstream of a dam, in case of flooding, is estimated as a function of the water level into the reservoir and of the forecasted discharge e/o volume to help water authorities in taking actions for its management and reduction. The proposed procedure firstly involves the derivation of single synthetic rainfall events—stochastically derived using a Monte Carlo procedure—through a bivariate copulas analysis; after this, synthetic bivariate stochastic inflow hydrographs are derived by using a conceptual fully distributed model starting from synthetic hyetographs. The flood hydrographs are then routed through the reservoir, taking into account the initial water level in the reservoir; flood inundation mapping through 2D hydraulic simulation and damage evaluation through depth–damage curves is finally carried out. This allows for the evaluation of the influence of initial water level on flood risk scenarios downstream.

The potential of this integrated procedure has been tested for the analysis of different flood management scenarios for the Castello reservoir in Sicily (Italy).

2. Materials and Methods

2.1. Case Study and Data

Castello dam is located in the southwestern part of Sicily (Italy), within the Magazzolo river catchment (Figure 1), and it is managed by the Sicily Regional Water Agency. According to the basic design information, the dam is characterized by a height of 50 m, a total volume of the reservoir of about 26 Mm³, and an impounded lake with an extension of 1.8 km².

The dam is an earth-fill type with a clay core, and it is 792 m long at the top; it was built between 1976 and 1985, and it is one of the largest reservoirs in western Sicily. The upstream catchment of the dam, a sub-catchment of the Magazzolo River, is characterized by an area of 81 km². The Digital Terrain Model (DTM) available for this area, extracted from a 20-m resolution DTM covering the whole Sicilian territory derived by an aerial photogrammetric survey, provides for the analysed sub-catchment at a maximum elevation of 1360 m above sea level (a.s.l.) and a minimum elevation of 250 m a.s.l. at Castello dam section. The normal water level of the reservoir is 293.65 m a.s.l., whereas the maximum water level is 296.65 m a.s.l.

Three rain gauge stations (Bivona, Giuliana, and Mezzojuso), managed by the Sicilian Agrometeorological Service (SIAS–Servizio Informativo Agrometeorologico Siciliano, www.sias.regione.sicilia.it) and located within or around the Magazzolo catchment (Figure 1), have been considered. For these stations, continuous rainfall time series, with a temporal resolution of 10 min, are available, respectively, for the period 2003–2020 for Bivona and Mezzoiuso and from 2005 to 2020 for the Giuliana rain gauge station. Moreover, a hydrometric ultrasonic gauge is installed on the dam crest measuring the lake water level and, indirectly, the flow discharges through the dam spillway.

2.2. The methodology

A Monte Carlo modelling chain (Figure 2) that includes hydrological simulations for synthetic flood hydrograph generation, reservoir routing, and hydraulic simulations for flood propagation and damage estimation has been implemented for the purposes of this study.

In particular, the developed modelling chain considers the following different steps: (a) Monte Carlo generation of ensembles of single synthetic rainfall events by using a stochastic bivariate model based on copulas that need sub-hourly rainfall time series as input; (b) generation of ensembles of single flood hydrographs by using a conceptual, fully distributed model fed with the above generated rainfall events; (c) flood hydrographs routing through the reservoir to obtain outflow hydrographs, for different reservoir conditions, in terms of initial water level (IWL); (d) flood propagation of the obtained outflow hydrographs and flood risk mapping through 2D hydraulic simulation; (e) flood damage modelling based on the results of the flood propagation model and on specific damage-depth curves. The modelling chain was entirely implemented in MATLAB environment with specific codes written for the purposes of the study.

2.2.1. Rainfall Generation Model

Single synthetic rainfall events have been generated by using a statistical model based on copulas theory, presented in detail in [21,30]. The model has a two-module structure: (a) volume-duration module for generation of the duration and the total volume of a single rainfall event; (b) storm temporal pattern module for the generation of the rainfall event profile.

For the model application, it is firstly required to identify independent rainfall events from the sub-hourly rainfall time series available as a function of the inter-event time that separates wet and dry periods. Rainfall time series can be represented as sequences of wet and dry periods where, during the wet period, a certain amount of rainfall is observed, while during the dry period, no rain (or less than a threshold value) is observed.

Between all the extracted rainfall events, it is then necessary to specify and select those considered as hydrologically “significant” for the possible formation of a flood, i.e., the maximum events. Kao and Govindaraju [31] stated how the definition of annual maximum events for multivariate problems is somewhat ambiguous. As matter of fact, extreme rainfall events could be defined as storms that have both high volume and peak intensity. Therefore, the definition of extreme rainfall based on events with annual maximum joint cumulative probability has been considered in this study.

For generating the duration and the total volume of a single rainfall event, a bivariate analysis on the empirical data needs to be performed by applying the well-known theory of copulas. For this application, Frank and Gumbel–Hougard copulas have been considered and adapted to the observed pairs of total volume and duration. Their formulas are respectively:

$$C\left(u,v\right)=-\frac{1}{\theta }ln\left[1+\frac{\left(e^{-\theta ·u}-1\right)\left(e^{-\theta ·v}-1\right)}{e^{-\theta }-1}\right],$$

(1)

$$C\left(u,v\right)=exp\left\{-{\left[{\left(-\ln u\right)}^{\theta }+{\left(-\ln v\right)}^{\theta }\right]}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$\theta $}\right.}\right\},$$

(2)

where θ is the parameter of the two copula functions estimated using the inversion of the Kendall’s coefficient method.

In order to select the copula that best represents the dependence structure of the empirical data, two graphical tools are used; primarily, the goodness of fit is tested by means of the K-plot, as defined by Genest and Rivest [32], where the values of the parametric function K(z) for both copulas are calculated and compared with the nonparametric function Kn(z) derived from the empirical data. The second graphical test is, instead, performed by comparing the level curves (isolines) of the theoretical copulas and those of the empirical ones.

Finally, the use of copulas requires the determination of marginal distributions based on univariate data. For this application, the fitting of several statistical distributions is performed by applying the maximum likelihood method, and the best fitted distribution is selected using various criteria, i.e., the Akaike information criterion (AIC), the relative root mean square error (RRMSE) and the Anderson–Darling test [33,34].

For the generation of the rainfall event profile, the storm temporal pattern module is implemented. This approach uses the mass curves concept [35] defined as the normalized cumulative rainfall volume versus normalized time since the event start. According to this definition, each event can be represented by pairs of two dimensionless variables, namely d and v, defined as: (a) d = t/D is the event dimensionless duration obtained by dividing the generic time t by the total rainfall duration D; (b) v = h(t)/V is event dimensionless cumulated rainfall volume obtained as the fraction of the cumulated rainfall depth h(t) at generic time t by the total rainfall volume V.

For the generation of the single synthetic rainfall event, dimensionless hyetograph can be chosen with a random picking from the set of the historical shapes.

2.2.2. Flood Hydrographs Generation Model

Generated rainfall events are used as input for the derivation of the correspondent flood hydrographs entering the reservoir by applying the fully distributed conceptual rainfall-runoff model proposed by Candela et al. [30].

The model has a grid-based structure where each small element (cell) has a distinct hydrological response, treated separately, but incorporates interactions with bordering cells. It is based on the linear kinematic mechanism for the transfer of the effective rainfall from different contributing areas to the catchment outlet and on a nonlinear approach for the rainfall excess calculation. In addition, model input represented by single rainfall events is considered in a fully distributed form.

The linear mechanism is represented by the distributed hydrological response function (UH) in a kinematic form, and it is described through a three-dimensional matrix H(m, n, p) (Figure 3), which represents the space-time distribution of contributing areas (isochrones areas).

The subscripts m = 1, 2 …, X and n = 1, 2, …, Y describe the elements of the matrix, while X and Y are the number of cells in which the catchment is discretized in the directions x and y. The subscript p = 1, 2, …, T counts the number of intervals (with T = Ω/∆t) in which the simulation time Ω is discretized. The matrix element H_m,n,p represents the area of generic cell (m, n) characterized by a concentration time ϑ_c(m,n).

For the evaluation of the concentration time at cell scale, the Wooding formula [36] is implemented in the model:

$${\vartheta }_{c, \left(m,n\right)}=\frac{L_{m,n\to out}^{3/5}}{k_{m,n\to out}^{3/5}·s_{m,n\to out}^{3/10}·r_{m,n}^{2/5}},$$

(3)

where L_m,n→,out [m] is the hydraulic path length between the centroid of the (m, n) cell and the outlet section of the catchment, k_m,n→,out [m^1/3/s] is the Strickler roughness for the same path, s_m,n→,out [m/m] is its slope, and r_m,n [m/s] is the average rainfall intensity for the rainfall event over the (m, n) cell.

The SCS-CN method proposed by USDA Soil Conservation Service [37] has been implemented in this model for the rainfall excess calculation (hydrological losses). The grid-based structure of the model allows direct incorporation of CN spatial distribution maps, taking into account the spatial variation of land use, soil type, and antecedent moisture conditions. The SCS-CN equations are also implemented in the model in a dynamic form [38], as the input rainfall can be variable in time:

$$P_e{}_{\left(q,m,n\right)}=\frac{{\left(P_{\left(q,m,n\right)}-c·S_{m,n}\right)}^2}{\left(P_{\left(q,m,n\right)}+\left(1-c\right)·S_{m,n}\right)},$$

(4)

and

$$S_{m,n}=254\left(\frac{100}{C{N}_{m,n}}-1\right),$$

(5)

Using above equations, a three-dimensional matrix Pe (q, m, n) of the same structure of H matrix, which represents the space-time distribution of excess rainfall, can be obtained.

The subscript q = 1, 2 …, N counts the elements of the matrix to N, which represents the number of steps in which the rainfall event of duration D (with N = D/∆t) is divided.

To compute the direct runoff hydrograph at the catchment outlet, the model solves the discrete convolution equation for a linear system as matrix multiplication:

$$Q=sum\left[diag\left(H\left(m,n,p\right)\times P_e\left(q,m,n\right)\right)\right]$$

(6)

where Q is the flood hydrograph vector with dimension N + T − 1.

Using this model for the generation of the flood hydrographs requires the calibration of the model parameters. Calibration is necessary for only two parameters: the coefficient c in the SCN-CN rainfall excess calculation formula (Equation (4)) and the Strickler roughness coefficient k in the Wooding formula (Equation (3)). For the other model parameters, i.e., path lengths L, average slopes s and the CN values, a calibration is not required, as their values can be directly extracted from the Digital Elevation Model and the CN map of the catchment. Finally, several metrics have been adopted to verify the goodness of the calibration.

2.2.3. Reservoir Routing Model for Discharged Hydrographs Derivation

The inflow hydrographs derived above have to be routed through the reservoir to obtain the outflow hydrographs from the spillway (peak flows, volumes, and total duration) for different initial conditions in terms of initial water level (IWL).

In level pool routing, the upstream discharge may be expressed explicitly in terms of the downstream discharge and of the channel or reservoir characteristics. Level pool routing is based on the continuity equation:

$$\frac{dW\left(t\right)}{dt}=Q_{in}\left(t\right)-Q_{out}\left[W\left(t\right)\right],$$

(7)

where Q_in(t) is the inflow hydrograph, Q_out[W(t)] is the outflow hydrograph, and W(t) is the reservoir storage. The problems for which this equation is applicable are given by Yevjevich [39]. Equation (7) has been here solved numerically using a fourth-order Runge–Kutta method implemented in a Matlab routine [39]. Particularly, the inflow discharge can be simulated using rainfall-runoff models, and the outflow discharge can be computed using spillway rating curve as follows:

$$Q_{out}\left[W\left(t\right)\right]={\mathrm{C}}_d{L}_e{H}^{3/2},$$

(8)

where C_d is the discharge coefficient assumed equal to 0.385, L_e is the effective length of spillway crest, and 𝐻 is the head on the spillway crest. The effective length of the spillway crest can be computed as follows:

$${\mathit{L}}_{\mathit{e}}=L_n-2\left(N_P{K}_P+K_a\right),$$

(9)

where L_n is the net length of the crest, N_p is the number of the piers, K_p is the pier contraction coefficient, and K_a is the abutment contraction coefficient.

2.2.4. Flood Propagation Modelling Downstream Reservoir

Outflow hydrographs for different reservoir conditions are then used for deriving flood hazard and risk maps downstream.

The Multilevel Flood Propagation 2D (MLFP-2D) model [40,41] has been implemented, in this study, to carry out flood inundation scenarios downstream. It is a hyperbolic model, whose details can be found in Candela and Aronica [42], based on the Saint–Venant equations and where the convective inertial terms are neglected. The model is based on two equations, conservative mass and momentum equations, which depend on the hydraulic resistances that can be calculated in function of the Manning–Strickler parameter. Both equations are solved using a finite-element triangular mesh that needs to be used to carry out model simulations.

Manning’s roughness coefficient is the unique calibration parameter involved in the MLFP-2D propagation model, and for its calibration, the procedure described in detail in Candela and Aronica [42] has been followed.

2.2.5. Flood Damage Evaluation

The evaluation of the total damage caused by the flooding in downstream area has been achieved by means of specific damage–depth curves (damage curves) derived considering that this area is mainly for agricultural use, as shown by the land use classes based on the CORINE Land Cover map for the Magazzolo catchment (Figure 4).

In this case, the damage is linked to the spatial and temporal variability of the flood that can cause strong loss of production and yield. For this reason, potential damage functions have been here obtained based on the methodology proposed by the Joint Research Centre, which considers both “new” costs (replacement) and productivity costs [29], starting from the knowledge of the land use classes of the area of interest.

Following this approach, firstly, the depth–damage functions, expressing the damage in terms of Euros in Purchasing Power parities (PPPs) proposed by Rusmini [43], have been considered (Figure 5). These functions relate nine depths (0, 0.5, 1, 1.5, 2, 3, 4, 5, and 6 or more meters) to the corresponding damage rates (from 0 to 1) for the most important land use classes of the CORINE Land Cover map. Consequently, these curves represent the relationship between water depth in the inundated area and the damage, in percent, for m² for the same purchasing power parity (PPP).

The determination of the PPP is then carried out through the consultation of the National Institute of Agricultural Economics (INEA) reports [44] in order to estimate the average values of agricultural land with the same purchasing power for the study area (

Table 1. Average values per mq of agricultural land with the same purchasing power in Sicily.

Main Agricultural Surface	Average Values per M2 (Euros)
Non irrigable arable land	0.9403
Fruit trees	1.9614
Permanently irrigated land	3.0224
Vineyards	1.5923
Annual and permanent crops	1.3474

).

As the damage is expressed in relative units (per square meter), the calculation can be carried out by averaging the water depths at finite element scale. Specifically, a unique value for each element can be calculated by averaging the three nodal values, and given the element area and the crop category associated with the element, the total damage is obtained.

3. Results

Single synthetic rainfall events have been generated, following the procedure described in Section 2.2.1, on the basis of a sample of 52 historical maximum annual rainfall events extracted as follows.

An inter-event time equal to 5 h was initially adopted to separate wet and dry periods in each sub-hourly rainfall series analysed; then, the events with annual maximum joint cumulative probability were selected for each analysed year for Bivona, Giuliana, and Mezzoiuso rain gauge stations. Finally, as all these rain gauges are in the same hydrologically homogeneous area [45,46], subsequent statistical analyses have been performed by aggregating all the selected events in a unique sample of 52 rainfall events whose main characteristics are reported in

Table 2. Main characteristics of the selected maximum rainfall events.

	Duration (min)	Volume (mm)	Iavg (mm/h)	Imax,30′ (mm/h)
Length of record (years)	17 (2003–2020)
Number of events	52
Max	3870	163.8	26.51	84.80
Min	120	19.2	0.64	9.20
Mean	955.58	60.84	6.12	37.76
Standard deviation	743.96	34.35	5.38	20.01

.

Gumbel–Hougard and Frank copula families were adapted to the observed pairs of total volume and duration. The θ parameter of these two copula functions was estimated using the inversion of Kendall’s coefficient method, and the results so obtained are reported in

Table 3. θ parameter for the analysed copulas.

Copula	θ
Gumbel–Hougard	1.2427
Frank	2.1407

.

To select the copula that best represents the dependence structure of the empirical data, as specified in paragraph 2.2.1, two graphical tools have been used. Firstly, the goodness of fit was tested by means of the K-plot, as defined by Genest and Rivest [32]. In this plot (Figure 6), the values of the parametric function K(z) for both copulas were calculated and compared with the non-parametric function K_n(z) derived from the empirical data.

A second graphical test was performed by comparing the level curves (isolines) of the theoretical copulas and those of the empirical copula (Figure 7).

Both tests confirmed how the Frank copula is well suited to describing the dependence structure between the empirical variables considered, and hence, it has been used for the rainfall generation model.

Moreover, the use of copulas requires the determination of marginal distributions based on univariate data. Therefore, fitting of several statistical distributions (i.e., Exponential, Gamma, Lognormal (LNII), Weibull, and General Extreme Values (GEV)) was considered by applying the maximum likelihood method, and the best fitted distribution was selected using various criteria. Simulations returned LNII distribution as the best marginal distribution for the rainfall volumes (Figure 8a) and Gamma distribution as the best marginal distribution for the rainfall duration (Figure 8b).

Temporal patterns of rainfall (event profile), for each event, were characterized using the mass curves concept, and the dimensionless hyetographs so derived for the selected annual maximum rainfall events (Figure 9a) have been used for the derivation of synthetic rainfall.

With these assumptions, 1500 synthetic rainfall events have been generated. This number was chosen as a trade-off between the statistical significance of the generated variables and the burden of the computational time requested for the simulations. About this, the comparison with the historical events is reported in Figure 9b, showing a good reproducibility of the rainfall main characteristics.

The events so generated have been used as input for the derivation of the correspondent flood hydrographs entering the reservoir by applying the flood hydrographs generation model illustrated in Section 2.2.2. The model was calibrated using the software PEST [47] on the basis of the outflow flood hydrograph through the spillway recorded by the water level gauge installed on the crest of the dam and of the rainfall recorded by the three rain gauges of the measurement network for the event of 25–26 February 2015 (Figure 10a). The parameter estimation software PEST, which is a combination of gradient descent and Newton’s method, implements the Gauss–Levenberg–Marquardt method [48,49] for parameter estimation and uncertainty analysis.

As mentioned before, the coefficient c in the SCN-CN rainfall excess calculation formula (Equation (4)) and Strickler roughness coefficient k in the Wooding formula (Equation (3)) are the only parameters that required calibration. The other model parameters needed to run the model have, instead, been extracted from the Digital Elevation Model (20 m resolution) and the CN map (100 m resolution available for the CN_II value) available for the sub-catchment (Figure 11).

More in detail, path lengths L (Figure 11a) and, consequently, average slopes s, have been extracted from the DTM map (Figure 1), whereas the CN values have been derived starting from the CN map available for the moderately wet soil moisture condition (CN_II) (Figure 11b).

Considering that most of the catchments in Sicily are small, with flashy hydrological response and a proneness to flash floods formation, especially when the soil is totally wet (CN_III condition), CN values considered to run the model are those relative to totally wet soil condition (CN_III values), derived from the CN_II values as follows [50]:

$$C{N}_{III}=\frac{23·C{N}_{II}}{10+0.13·C{N}_{II}}$$

(10)

Despite the distributed nature of the Strickler roughness coefficient k, for the sake of simplicity, a “lumped” value was calibrated by considering a spatial average over the whole catchment.

Final results of the model calibration are reported in Figure 10b, where the observed and modelled flood hydrographs are plotted for the optimal values of the two parameters. Their final values are c = 0.155 and k = 19.8 m^1/3/s, respectively.

The calibration efficiency was measured through several metrics; in particular, we selected the Mean Absolute Error (MAE), the Root mean Standard deviation Ratio (RSR), and the Nash-Sutcliffe (NSE) indexes to compare errors between observed and estimated discharge values. As a result, we obtained a MAE equal to 3.647 m³/s, a RSR index equal to 0.283, and a Nash-Sutcliffe index equal to 0.921, with an error in-peak discharge of 0.94% and in-flood volume of −4.96%.

The flood hydrograph generation module was, hence, used to generate 1500 flood hydrographs (Figure 12) entering the reservoir correspondent to the 1500 synthetic rainfall events above generated.

The inflow hydrographs routed through the reservoir, by applying the procedure described in Section 2.2.3, allowed obtaining the outflow hydrographs from the spillway (peak flows, volumes, and total duration) for different initial conditions in terms of Initial Water Level (IWL). For a broader analysis, two specific values were chosen: the highest (IWL1), equal to 293.65 m a.s.l., corresponding to normal water level of the reservoir and the lowest (IWL2), equal to 290.00 m a.s.l., corresponding to the most frequent observed level at the dam.

The outflow hydrographs, for the two different reservoir conditions, were finally used for deriving flood hazard and risk maps downstream from the Castello dam by applying MLFP-2D model and for flood damage evaluation.

The model equations have been solved using a finite-element triangular mesh. The mesh covers a domain area of 1.68 km² downstream, discretized into 25.436 nodes and 49.278 elements. The terrain elevations for the study area were derived starting from a 2-m resolution DTM interpolated from a LIDAR survey available for the floodplain. Manning’s roughness coefficient was the unique calibration parameter involved in the propagation model; particularly, one coefficient for each triangular element can be chosen but, lacking a robust basis for allowing the roughness coefficient to vary, the entire triangular domain was divided into two principal regions—the floodplain area and the river—and for both regions, a calibrated Manning roughness coefficient was considered (0.037 s.m^−1/3 for the river and 0.051 s.m^−1/3 for the floodplain area). As an example, in Figure 13, the domain DTM (Figure 13a) and the flood inundated area for a given reservoir condition (Figure 13b), corresponding to normal water level in the reservoir, are reported.

The evaluation of the total damage caused by the flooding has been carried out by means of the specific damage–depth curves (damage curves) derived, as presented, in Section 2.2.5. Particularly, the damage–depth curves illustrated in Figure 5 have been suitably combined with the average values per m² of agricultural land with the same purchasing power of the area where the study area is located, obtaining the damage curves representing the total damage per m² in function of the water depth for the main agricultural landcover classes of the Magazzolo floodplain (Figure 14).

Finally, for each simulation, the total direct flood damage has been calculated for the two initial reservoir conditions (IWL1 and IWL2) considered.

In Figure 15, the results of the simulations are reported by plotting, with a scale of colour for the total damage, the dots representing the inflow hydrographs (flood peak discharge–flood volume pairs, specifically). The choice of representing these results in relation to the hydrological forcing to the reservoir is due to the technical/practical aspects, which the proposed procedure intends to address.

For these aspects, which will be further discussed, and for a better visualization, the clouds of points have been visualized as smooth surfaces obtained by a 2-D spatial interpolation (Figure 16).

4. Discussion

By looking more in detail into Figure 15, a strong correlation between the total damage and the flood peaks–volume pairs, in the case of IWL1 condition, is revealed, while this correlation is definitively weaker (or likely absent) in the case IWL2 condition.

This behaviour appears to be reasonable in the view of the two specific conditions: for IWL1 condition, the reduction effect due to the reservoir volume is lower than the IWL2 condition and, hence, all the flood hydrographs are routed, discharged downstream from the dam, and the original correlation structure between hydrographs characteristics (peak discharges and volumes) is essentially preserved. In the IWL2 condition, instead, some hydrographs are retained and not routed through the reservoir. This results in less events producing flood inundations and, generally, lower values of total damages.

Interpolated 2-D “damage surfaces” (Figure 16), corresponding to the two IWL conditions, can be used for quantifying the expected damage downstream from the dam for a given input hydrograph by simply entering the plot with a given pair of flood peak-volume.

In both cases, the higher damage level can be observed only when both variables (flood peaks and volume) show the major values. Again, this evidence is clearer for the IWL1 condition than IWL2, due to the fact of the influence of the reservoir (flood control) on the routed hydrographs.

In detail, these surfaces clearly show these specific results, whereas the regions characterized by higher values of total damage are more extended in IWL1 condition than IWL2 and vice versa. Regions characterized by lower values of total damage are more extended in IWL2 condition than IWL1.

To better clarify these statements, consider two distinct hydrological scenarios for the reservoir, i.e., the 50-yrs and 100-yrs return time input hydrographs, with each characterized by a specific pair of peak discharge and volume, obtained by bivariate analysis of the generated sample. The corresponding values are: 653.0 m³/s and 7.92 Mm³ for 50-yrs return time; 824.0 m³/s and 8.99 Mm³ for 100-yrs return time.

Now, by entering the surface plot with these two pair, it is possible to evaluate the damage, which is equal to 811,862 Euro for IWL1 condition and 560,765 Euro for IWL2 condition, both for 100-yrs return time. Similar results can be obtained for 50-yrs return time, as shown in

Table 4. Total damage in Euro for different hydrological input and IWLs.

Return Time	Qmax/Vtot	IWL1 (293.65 m)	IWL2 (290.00 m)
50 years	653.0 m3/s	€ 735,738	€ 354,563
	7.92 Mm3
100 years	824.0 m3/s	€ 811,862	€ 560,765
	8.99 Mm3

.

Hence, lowering the reservoir level produces a difference of about 30.9% in total direct damage downstream from the dam for 100-yrs return time and 51.8% for 50-yrs return time.

As matter of fact, the comparison between the results for the two IWL conditions shows how it is possible to have a direct insight on the role which the reservoir plays in protecting the downstream floodplain and to prove its capability to control inundation and reduce direct flood damage.

Damage surfaces allow us to quantify the impact of the different reservoir conditions on the direct flood damage downstream, and they can help the flood risk managers in setting up specific operation rules for the reservoir in order to mitigate the impact of extreme hydro meteorological events.

5. Conclusions

In this study, an efficient and reliable Monte Carlo modelling chain for producing flood risk scenarios downstream to a dam, as a function of the initial water level in the reservoir, has been presented. The methodology allows us to quantify the total direct flood damage in a floodplain downstream of a dam through the generation of ensembles of single rainfall events by using a copula based model, of flood hydrographs by using a conceptual fully distributed rainfall-runoff model, as well as of discharged hydrographs from a reservoir routing and a two-dimensional hydraulic flood propagation model.

Flood damage scenarios have been produced for the Magazzolo river floodplain, downstream the Castello dam in Sicily, using 1500 stochastic extreme rainfall events. Damage analysis for the different flood damage scenarios revealed how the reservoir rules operation (in terms of initial reservoir conditions or IWL) strongly influences the effects of flooding in the downstream floodplain. In particular, the performed analysis reveals how there is a strong correlation between the total damage and the flood peaks–volume pairs in the case of normal water level conditions, while the correlation is definitively weaker in the case of a more frequent water level.

Finally, the a priori knowledge of the possible damage associated with a forecasted or measured flood event, as a function of the water level of the reservoir, can be very useful in case of flood warning to help the Water Authorities to quantify the potential downstream damage and to make appropriate decisions, i.e., lowering the reservoir level for the mitigation of the damage downstream.