Flood Control Versus Water Conservation in Reservoirs: A New Policy to Allocate Available Storage

: The aim of this study is to contribute to solving conﬂicts that arise in the operation of multipurpose reservoirs when determining maximum conservation levels (MCLs). The speciﬁcation of MCLs in reservoirs that are operated for water supply and ﬂood control may imply a reduction in the volume of water supplied with a pre-deﬁned reliability in the system. The procedure presented in this study consists of the joint optimization of the reservoir yield with a speciﬁc reliability subject to constraints imposed by hydrological dam safety and downstream river safety. We analyzed two di ﬀ erent scenarios by considering constant or variable initial reservoir level prior to extreme ﬂood events. In order to achieve the global optimum conﬁguration of MCLs for each season, we propose the joint optimization of three variables: minimize the maximum reservoir level (return period of 1000 years), minimize the maximum released outﬂow (return period of 500 years) and maximize the reservoir yield with 90% reliability. We applied the methodology to Riaño Dam, jointly operated for irrigation and ﬂood control. Improvements in the maximum reservoir yield (with 90% reliability) increased up to 10.1% with respect to the currently supplied annual demand (545 hm 3 ) for the same level of dam and downstream hydrological safety. The improvement could increase up to 26.8% when compared to deterministic procedures. Moreover, dam stakeholders can select from a set of Pareto-optimal conﬁgurations depending on if their main emphasis is to maintain / increase the hydrological safety, or rather to maintain / increase the reservoir yield.


Introduction
Owing to increasingly risk-averse societies, stronger hydrological safety requirements are being imposed on existing dams in order to fulfil new regulations and prevent dam failures [1]. This problem can be addressed with two different kinds of technical solutions: hard solutions, as the alteration of dam spillways or elevation of dam crest and soft solutions, as the allocation of additional flood control volumes in the reservoir. Hard solutions increase the flood control capacity while maintaining the reservoir storage available for water supply. However, their main drawback is the need to allocate resources for infrastructure works. On the other hand, the implementation of soft solutions is easier and quicker, but can reduce the available volume to supply water with a specific reliability in a water system.
Allocation of flood control volume is addressed by defining a maximum conservation level (MCL), also known as flood-limited water level [2] or flood control level. MCL is the maximum operating level

•
Stochastic assessment of hydrological dam and downstream river safety through return periods related to maximum reservoir levels and maximum outflows.

•
Determination of MCLs that increase/maintain the water yield for a specific reliability while maintaining/improving hydrologic dam safety.
Determination of MCLs accounting for the variability of initial reservoir level prior to flood events.The methodology is illustrated through its application to a gated spillway dam located in Spain.

Materials and Methods
Maximum conservation levels represent the linking variable between flood control operation and water conservation operation of the reservoir (Figure 1). Allocation of flood control volume is addressed by defining a maximum conservation level (MCL), also known as flood-limited water level [2] or flood control level. MCL is the maximum operating level that the reservoir is allowed to reach under regular operation conditions. This reservoir level is below or equal to the maximum normal operating level (MNL) and can vary along the seasons of the year.
MCL is the most significant parameter in the trade-off established between flood control and water supply when increasing flood control volumes by soft solutions [3]. Traditionally, practitioners defining MCL only focused on hydrological dam and downstream safety. They frequently neglected other purposes of the reservoir, such as water supply and the economic consequences derived from loss of water yield reliability [4]. Some authors [3,5,6] have focused on accounting simultaneously for both regular (associated to water supply purposes of the reservoir) and flood control dam operations when defining MCLs. These studies analyzed hydrological dam safety by applying deterministic procedures, in which the return period associated to dam and downstream safety is assumed to be equal to the one associated to the flood event. Several authors [7,8] pointed out that hydrological dam safety and downstream safety should be assessed by analyzing the return periods of the maximum reservoir water levels and maximum outflows respectively.
Another relevant factor is that practitioners usually define MCLs assuming that the reservoir is at the maximum level under normal operating conditions prior to flood arrival [9][10][11]. This hypothesis results in conservative hydrological safety assessments. However, it can reduce the volume of water available to satisfy the demands of the water supply system because it influences the definition of MCLs. Accounting for the variability of the initial reservoir level in a hydrological dam and its downstream safety leads to more realistic results [8,[12][13][14], which consequently can improve the definition of MCLs. Within this study, we propose a stochastic methodology to determine seasonal MCLs. The methodology combines three main innovative aspects: • Stochastic assessment of hydrological dam and downstream river safety through return periods related to maximum reservoir levels and maximum outflows.
• Determination of MCLs that increase/maintain the water yield for a specific reliability while maintaining/improving hydrologic dam safety. Determination of MCLs accounting for the variability of initial reservoir level prior to flood events.The methodology is illustrated through its application to a gated spillway dam located in Spain.

Materials and Methods
Maximum conservation levels represent the linking variable between flood control operation and water conservation operation of the reservoir ( Figure 1).  COD represents the level of the crest of dam, DFL the design flood level, MNL the maximum normal operating level, and Zo the reservoir level prior to the flood event. Volume above MCL corresponds to the flood control operation volume and below MCL corresponds to the water conservation operation volume. We propose a stochastic methodology to obtain the optimal set of seasonal MCLs accounting for both hydrological safety (dam and downstream) and water supply with a specific reliability. Figure 2 shows a scheme relating the main elements and procedures developed and applied in the methodology.
Water 2020, 12, x FOR PEER REVIEW 3 of 13 to the flood control operation volume and below MCL corresponds to the water conservation operation volume.
We propose a stochastic methodology to obtain the optimal set of seasonal MCLs accounting for both hydrological safety (dam and downstream) and water supply with a specific reliability. Figure  2 shows a scheme relating the main elements and procedures developed and applied in the methodology.

Study Set of Seasonal Maximum Conservation Levels
We defined a study set of seasonal MCLs representative of all possible configurations. The number of possible configurations of MCLs varies according to the number of seasons identified (n) and the number of possible maximum conservation levels selected (k) for a proper discretization. As the season of the MCLs matters (it is not the same to have the same MCL in one season or another), the number of possible configurations is k n . The number of seasons is defined as described in Section 2.2.1.

Reservoir Operation Simulations
For each configuration of seasonal MCLs of the study set, we carried out two different reservoir operation simulations: simulation of flood control operation and simulation of water conservation operation.

Simulation of Flood Control Operation
Gabriel-Martin et al. [13] presented a stochastic methodology that enabled us to obtain stochastic inflow hydrographs representative of the observed daily annual floods ( Figure 3). The main steps used to generate the inflow hydrographs were as follows: • Generation of 100,000 pairs of flood duration (D) with their associated maximum annual flood volume (V). Pairs of 100,000 flood event durations were generated following the empirical

Study Set of Seasonal Maximum Conservation Levels
We defined a study set of seasonal MCLs representative of all possible configurations. The number of possible configurations of MCLs varies according to the number of seasons identified (n) and the number of possible maximum conservation levels selected (k) for a proper discretization. As the season of the MCLs matters (it is not the same to have the same MCL in one season or another), the number of possible configurations is k n . The number of seasons is defined as described in Section 2.2.1.

Reservoir Operation Simulations
For each configuration of seasonal MCLs of the study set, we carried out two different reservoir operation simulations: simulation of flood control operation and simulation of water conservation operation.

Simulation of Flood Control Operation
Gabriel-Martin et al. [13] presented a stochastic methodology that enabled us to obtain stochastic inflow hydrographs representative of the observed daily annual floods ( Figure 3). The main steps used to generate the inflow hydrographs were as follows: • Generation of 100,000 pairs of flood duration (D) with their associated maximum annual flood volume (V). Pairs of 100,000 flood event durations were generated following the empirical probability distribution of historical floods. For each element of the 100,000 generated durations, the corresponding hydrograph volume was obtained following the probability distribution of the associated duration within a Monte Carlo framework. basin. By applying the curve number method inversely [15] to the cumulated net precipitation, the value of the cumulated precipitation depth was obtained ( Figure 3a) • Temporal distribution of the 100,000 cumulated precipitation depth values. Each cumulated rainfall depth was distributed temporally by applying an autoregressive moving average (ARMA) (2,2) model [10]. Thus, 100,000 hourly hyetographs were obtained (Figure 3b).

•
Generation of 100,000 hourly-distributed hydrographs. By applying the curve number method [15] and the soil conservation service dimensionless unit hydrograph procedure [15], 100,000 hydrographs were generated, which followed the empirical probability distributions of volume and duration (Figure 3b).
• Generation of 100,000 hourly-distributed hydrographs. By applying the curve number method [15] and the soil conservation service dimensionless unit hydrograph procedure [15], 100,000 hydrographs were generated, which followed the empirical probability distributions of volume and duration (Figure 3b). A detailed description of the method to generate representative hydrographs can be found in Gabriel-Martin et al. [13]. Each hydrograph was associated to one of the seasons defined ( Figure 3c). We identified the distinct seasons by applying a graphical test to the observed daily inflows proposed by Ouarda [16] and Ouarda et al. [17]. This test was based on a peak-over-threshold (POT) analysis. In order to identify the seasons, Ouarda et al. [17] tested different threshold values while assuring independence between the selected POT. We applied the criteria recommended by Lang et al. [18] to define the range of threshold values to be considered and the criteria proposed by the Water Resources Council [19] to assure the independence between two consecutive floods. We plotted the cumulative empirical probability of POT during the year against the time of the year for each threshold value tested. The slope changes within the plot indicated the significant seasons. A detailed description of the method to identify the seasons can be found in Gabriel-Martin et al. [20] and is summarized in Figure  3c. For each configuration of MCLs, we obtained the maximum water level in the reservoir corresponding to a return period of 1000 years (MWRLTR=1000y) and the maximum outflow . The procedure to generate the ensemble of inflow hydrographs (a,b) and its seasonal characterization (c). D represents the flood durations, whereas V represents the maximum annual flood. S1, S2, and S3 represent the different seasons.
A detailed description of the method to generate representative hydrographs can be found in Gabriel-Martin et al. [13]. Each hydrograph was associated to one of the seasons defined ( Figure 3c). We identified the distinct seasons by applying a graphical test to the observed daily inflows proposed by Ouarda [16] and Ouarda et al. [17]. This test was based on a peak-over-threshold (POT) analysis. In order to identify the seasons, Ouarda et al. [17] tested different threshold values while assuring independence between the selected POT. We applied the criteria recommended by Lang et al. [18] to define the range of threshold values to be considered and the criteria proposed by the Water Resources Council [19] to assure the independence between two consecutive floods. We plotted the cumulative empirical probability of POT during the year against the time of the year for each threshold value tested. The slope changes within the plot indicated the significant seasons. A detailed description of the method to identify the seasons can be found in Gabriel-Martin et al. [20] and is summarized in Figure 3c.
For each configuration of MCLs, we obtained the maximum water level in the reservoir corresponding to a return period of 1000 years (MWRL TR=1000y ) and the maximum outflow corresponding to a return period of 500 years (MO TR=500y ). These values correspond to the set of 100,000 seasonal maximum annual inflow hydrographs generated by Gabriel-Martin et al. [13,20].
The analysis was carried out in two different scenarios of the initial reservoir level prior to the flood event (Zo) (Figure 2): • Scenario 1 (Sc.1): Zo is constant and corresponds to the seasonal MCL. The reservoir is assumed to be at its maximum operating level when the maximum annual flood occurs. • Scenario 2 (Sc.2): Zo is variable, and the reservoir can be at any level when the maximum annual flood occurs. Zo is randomly sampled from the cumulative probability distribution of Zo associated to the season of occurrence of the maximum annual flood event. This distribution is obtained from the simulation of the water conservation operation of the reservoir (Section 2.2.2).
The operation of the dam gates was simulated by applying the volumetric evaluation method (VEM), fully described in Giron [21] and Sordo-Ward et al. [11,22].

Simulation of Water Conservation Operation
We considered the yield reliability (Y R ) as the ratio of total volume supplied and the total volume demanded [23][24][25][26] during the period analyzed. For each configuration of MCLs (k), we obtained the reservoir yield with a reliability of 90% (Y R = 90% ) and the cumulative probability distribution of Zo, by developing a monthly water balance model. Reservoir storage is large compared to monthly inflows, and thus the monthly time scale is appropriate. The model manages the dam as an isolated element and applies rules of operation validated in Gabriel-Martin et al. [13]. The water balance considers time series of monthly inflows, environmental flow restrictions, evaporation rates, monthly demand distribution, storage-area-height reservoir curves, and dead storage volume (data extracted from "Duero National Water Master Plan" [27]). The main purpose of the reservoir is irrigation and therefore we adopted a required yield reliability of 90% (Y R = 90% ,) which is adequate for irrigation demands in the region [28].
To identify the maximum amount of water that can be supplied to satisfy a regular demand with a specified reliability, a bipartition method was applied. Excessive values of demands were set (for example, similar to mean monthly runoff) and the simulation was carried out. The deficits were obtained and specified yield reliability requirements were checked. If the specified reliability requirements were not fulfilled, the demand was reduced by half and simulated again. If the specified reliability requirements were satisfied, half of the difference was added and simulated again and so on, until the deficit (or gain) was smaller than a pre-set tolerance (e.g., 0.1 hm 3 /year). In addition, we simulated the operation of the reservoir with the associated mean annual current demand. We repeated the procedure for both Sc.1 and Sc.2 scenarios.

Results Analysis and Solutions Proposal
In order to propose the optimal configurations of MCLs within the case study, as exposed, we selected three main decision variables: MWRL TR = 1000y , MO TR = 500y , and Y R = 90% . We assumed that the MCLs configuration would not fulfil the standards if MWRL TR=1000y was above the design flood level (DFL) and/or MO TR = 500y was greater than the emergency flow (O EMER ). We compared MWRL TR=1000y , MO TR = 500y , and Y R = 90% for all the configurations in the study set of MCLs by identifying configurations that are non-inferior solutions (Pareto framework) in terms of maximum volume of water supplied (with a specific yield reliability of 90%) and hydrological safety.

Determination of Possible MCLs by Applying a Pareto Analysis
The selected variables for conducting the Pareto analysis were MO TR = 500y and Y R = 90% (Figure 4). In this study, we assumed that higher levels in the reservoir (and above MNL) imply greater outflows. This assumption is always fulfilled either for dams with fixed crest spillways or if the VEM is applied in gated spillways (as in this study). It should be noted that this assumption may not hold for all possible specific characteristics of dams, all rules of operation adopted, or all specific values adopted for the hydrological dam and downstream safety. Moreover, in this study, the outflow corresponding to the DFL condition is higher than that corresponding to O EMER , that is, the MO TR = 500y is the most restrictive variable. Therefore, for each scenario (Sc.1 or Sc.2), we had a set of k n MCLs configurations with k n pairs of values MO TR = 500y and Y R = 90% . The purpose was to obtain the configurations of MCLs that minimize the value of MO TR = 500y while maximizing Y R = 90%, which implies a two-objective minimization problem (Equation (1)):  (1) consisted of a set of non-dominated solutions. Therefore, following the procedure proposed in Chong and Zak [29] we conducted the mentioned analysis.
Once the non-dominated solutions were identified for Sc.1 and Sc.2, we eliminated from both scenarios those that did not fulfil the hydrological safety regulation standards (both for MWRL TR = 1000y and/or MO TR = 500y ) and those providing a Y R = 90% lower than the annual demand that is currently satisfied. Afterwards, we compared the proposed solutions, providing dam stakeholders with a set of possible configurations depending on whether their main objective was to increase hydrological dam and downstream river safety or increase the water supply with a specific reliability within the system. corresponding to the DFL condition is higher than that corresponding to OEMER, that is, the MOTR = 500y is the most restrictive variable. Therefore, for each scenario (Sc.1 or Sc.2), we had a set of k n MCLs configurations with k n pairs of values MOTR = 500y and YR = 90%. The purpose was to obtain the configurations of MCLs that minimize the value of MOTR = 500y while maximizing YR = 90%, which implies a two-objective minimization problem (Equation (1)):  (1) consisted of a set of non-dominated solutions. Therefore, following the procedure proposed in Chong and Zak [29] we conducted the mentioned analysis. Once the nondominated solutions were identified for Sc.1 and Sc.2, we eliminated from both scenarios those that did not fulfil the hydrological safety regulation standards (both for MWRLTR = 1000y and/or MOTR = 500y) and those providing a YR = 90% lower than the annual demand that is currently satisfied. Afterwards, we compared the proposed solutions, providing dam stakeholders with a set of possible configurations depending on whether their main objective was to increase hydrological dam and downstream river safety or increase the water supply with a specific reliability within the system.

Case Study
We applied the methodology to Riaño Dam (Table 1). It belongs to the Esla basin water resources system, which is managed by the Duero River Basin Authority (west region of mainland Spain). The main purpose of the reservoir is irrigation. The mean monthly water demands for irrigation/urban water supply are as follows (in hm 3

Case Study
We applied the methodology to Riaño Dam (Table 1). It belongs to the Esla basin water resources system, which is managed by the Duero River Basin Authority (west region of mainland Spain). The main purpose of the reservoir is irrigation. The mean monthly water demands for irrigation/urban water supply are as follows (in hm 3 ): April 10.9/0.08, May 54.5/0.08, June 81.7/0.08, July 158/0.16, August 147.1/0.16, and September 92.6/0.16. From October to March, the demands were 0/0.08 hm 3 . The capacity of the reservoir is 651 hm 3 at MNL (the maximum reservoir level that water might reach under normal operating conditions) [30], with a bottom dead storage of 78 hm 3 . The main characteristics of the Riaño Reservoir and its basin are summarized in Table 1. Riaño Dam has two spillways. The main spillway is controlled by two tainter gates, each eight meters wide and seven meters high. There is a second spillway for emergency purposes. It is a fixed-crest spillway with the crest located at the MNL. Riaño Dam also has two bottom and two intermediate outlets, which we assumed were closed during the floods. Flood damage analyses summarized in the Dam Master Plan concluded that discharges above 700 m 3 /s (O EMER ) could produce damage over urban settlements with more than five inhabitants and infrastructures in the downstream reach. The following data were used to perform the study:

•
Simulation of flood control operation. Besides the flood control structures and dam configuration shown, we used 30 years of unaltered daily flow series data from a gauge located right downstream the Riaño reservoir (from the years 1954 to 1984 and prior to the existence of the dam). With this time series, the 100,000 seasonal synthetic flood hydrographs were generated.

•
Simulation of water conservation operation. We used a monthly time series of naturalized inflows from 1940 to 2013, environmental flow restrictions, evaporation rates, monthly demand distribution, storage-area-height reservoir curves and dead storage volume (all data obtained from the Duero River Basin Management Plan [27]) and the reservoir characteristics previously stated.

Limitations of the Methodology
We applied this methodology to one basin and dam configuration. This might limit the generalization of the results obtained. Furthermore, the water resources management model focused on the regular operation of the dam as an isolated element. This methodology could be extended to take into consideration the interaction with other infrastructures within the system, using suitable water resources management models (e.g., AQUATOOL [31] and WEAP [32]).

Determination of the Study Set of MCLs to be Analyzed
We studied a set of MCLs that ranged from 651 hm 3 (volume at MNL) to 400 hm 3 . For the sake of simplicity, we estimated the flood hydrograph volume of Tr = 5000 years (247 hm 3 ) and defined a maximum flood control volume of 251 hm 3 . We discretized the ranges of reservoir volumes in intervals of 5 hm 3 . Thus, we defined a set of k = 51 possible MCLs per season (associated to a volume in the reservoir of 400, 405, . . . , 645, 651 hm 3 ).
According to Gabriel-Martin et al. [20], three characteristic seasons (regarding maximum annual floods) were identified for the location of Riaño: season 1 (S1) from the beginning of November to the end of January; season 2 (S2) from the beginning February to the end of April; and season 3 (S3) from the beginning of May to the end of October. Therefore, as n = 3 seasons, we had a set of 51 3 = 132,651 configurations of MCLs per scenario analyzed (Sc.1 and Sc.2.)

Simulation of the Water Conservation and Flood Operation of the Dam
Once the configurations of MCLs were defined, for each configuration we obtained the values of MWRL TR = 1000y and MO TR = 500y by simulating the flood operation of the reservoir for Sc.1 and Sc.2. By simulating the water conservation operation of the reservoir, we obtained Y R = 90% . Figure 5 shows the values of MWRL TR = 1000y , MO TR = 500y , and Y R = 90% with respect to the MCLs of each season.
Water 2020, 12, x FOR PEER REVIEW 8 of 13 By simulating the water conservation operation of the reservoir, we obtained YR = 90%. Figure 5 shows the values of MWRLTR = 1000y, MOTR = 500y, and YR = 90% with respect to the MCLs of each season.  Figure 5a shows that variation of MCLs in S1 did not affect YR=90%. This is because the water demands associated with the reservoir in S1 were less than 1% of the annual demand (Da), and, as the irrigation season in the Riaño system extends from May to September (98% of Da), the reservoir was able to recover the reduced volume in S1 within the months previous to irrigation (February, March, and April). Figure 5b,d shows that, for the case of Sc.1, the effects of MCLs were similar in the three seasons in terms of hydrological dam safety (Figure 5b) and downstream river safety ( Figure  5d). However, in Sc. 2 (Figure 5c,e), as regular operation was linked to the flood control operation by the initial reservoir level, MCLs in S1 and S3 did not affect either the hydrological dam safety or  Figure 5a shows that variation of MCLs in S1 did not affect Y R=90% . This is because the water demands associated with the reservoir in S1 were less than 1% of the annual demand (Da), and, as the irrigation season in the Riaño system extends from May to September (98% of Da), the reservoir was able to recover the reduced volume in S1 within the months previous to irrigation (February, March, and April). Figure 5b,d shows that, for the case of Sc.1, the effects of MCLs were similar in the three seasons in terms of hydrological dam safety (Figure 5b) and downstream river safety (Figure 5d). However, in Sc.2 (Figure 5c,e), as regular operation was linked to the flood control operation by the initial reservoir level, MCLs in S1 and S3 did not affect either the hydrological dam safety or downstream river safety. Variations of MCL S2 are the main affection to values MWRL TR = 1000y and MO TR = 500y in Figure 5c,e, respectively.

Solutions Proposal
First, we identified the configurations which did not fulfil hydrological dam safety (MWRL TR = 1000y > DFL) and/or downstream safety (MO TR = 500y > O EMER. ) in both scenarios. A total of 826 configurations (0.6% of 132,651 configurations) had a value of MWRL TR = 1000y higher than DFL in Sc.1 (red dots in Figure 6a), while none of the configurations had MWRL TR = 1000y values higher than the DFL in Sc.2 (Figure 6c). It should be noted that the 826 configurations that did not fulfil hydrological dam safety in Sc.1 also did not fulfil the downstream safety condition. On the other hand, 59,013 configurations (44.5% of the total number of configurations) had a value of MO TR = 500y higher than O EMER. in Sc.1, whereas 4944 (3.7%) in Sc.2 (Figure 6a,c) shows the pair of values Y R = 90% and MO TR = 500y (grey points) for each analyzed configuration in Sc.1 and Sc.2, respectively.
Water 2020, 12, x FOR PEER REVIEW 9 of 13 downstream river safety. Variations of MCLS2 are the main affection to values MWRLTR = 1000y and MOTR = 500y in Figure 5c,e, respectively.

Solutions Proposal
First, we identified the configurations which did not fulfil hydrological dam safety (MWRLTR = 1000y > DFL) and/or downstream safety (MOTR = 500y > OEMER.) in both scenarios. A total of 826 configurations (0.6% of 132,651 configurations) had a value of MWRLTR = 1000y higher than DFL in Sc.1 (red dots in Figure 6a), while none of the configurations had MWRLTR = 1000y values higher than the DFL in Sc.2 (Figure 6c). It should be noted that the 826 configurations that did not fulfil hydrological dam safety in Sc.1 also did not fulfil the downstream safety condition. On the other hand, 59,013 configurations (44.5% of the total number of configurations) had a value of MOTR = 500y higher than OEMER. in Sc.1, whereas 4944 (3.7%) in Sc.2 (Figure 6a,c) shows the pair of values YR = 90% and MOTR = 500y (grey points) for each analyzed configuration in Sc.1 and Sc.2, respectively. In Sc.1, the Pareto-solutions consisted of 277 configurations (Figure 6a, magenta points). One those, 98 configurations satisfied the following: YR = 90% ≥ Da, MOTR = 500y ≤ OEMER, and MWRLTR = 1000y ≤ DFL (proposed solutions, cyan in Figure 6a). For the same analysis in Sc.2, we identified 247 configurations (magenta points in Figure 6c) and 135 (cyan points in Figure 6c), respectively. For both scenarios, the extreme proposed solutions were identified. On one hand, in the case of MOTR = 500y = OEMER = 700 m 3 /s (Figure 6a,c, black point), YR = 90% = 587 hm 3 (for Sc.1) and 600 hm 3 (for Sc.2) representing an improvement (compared to Da) of 7.7% and 10.1%, respectively. On the other hand, in the case of YR = 90% = Da = 545 hm 3 (Figure 6a,c, dark blue point), MOTR = 500y = 452 m 3 /s (for Sc.1) and  Figure 6a,c, respectively. In both scenarios, the highest MCLs were associated with S3, while the lowest was associated with S1.

Comparison between the Proposed Configurations in the Two Scenarios
We compared the limit proposed solutions in both scenarios (Sc.1 and Sc.2; Figure 7). In the case of MO TR = 500y = O EMER = 700 m 3 /s, MWRL TR = 1000y was 1100.6 m.a.s.l. The hydrological dam and downstream river safety were invariant for Sc.1 and Sc.2. However, higher MCLs were obtained for Sc.2, which increased the Y R = 90% of the system. In the case of Y R = 90% = Da = 545 hm 3 , the differences between Sc.2 and Sc.1 for MWRL TR = 1000y and MO TR = 500y were 0.2 m.a.s.l and 86 m 3 /s, respectively. Moreover, the MCL of season two were lower when the variable initial reservoir level was considered. This is because of the increase of MCLs in the other seasons. Despite this, the demand supplied with a reliability of 90% was the same, accounting for lower values of MWRL TR=1000y and MO TR = 500y if variable initial level was considered.
Water 2020, 12, x FOR PEER REVIEW 10 of 13 366 m 3 /s (for Sc.2) representing a 64.2% and 52.3% of OEMER, respectively. It is important to point out that, in the case of using a deterministic procedure focused on hydrological dam and downstream safety (conventional procedure), any studied configuration of MCLs could be a potential solution (grey dots within Figure 6a,c). Thus, if MOTR = 500y = OEMER, the proposed stochastic procedure presented improvements of up to 146 hm 3 (26.8% of Da) compared to the worst regular operation configuration (YR = 90% = 454 hm 3 in Sc.1, Figure 6a). The corresponding configurations of MCLs for the extreme proposed solutions are represented in Figure 6b (Sc.1) and Figure 6d (Sc.2) with the same color scheme as in Figure 6a and 6c, respectively. In both scenarios, the highest MCLs were associated with S3, while the lowest was associated with S1.

Comparison between the Proposed Configurations in the Two Scenarios
We compared the limit proposed solutions in both scenarios (Sc.1 and Sc.2; Figure 7). In the case of MOTR = 500y = OEMER = 700 m 3 /s, MWRLTR = 1000y was 1100.6 m.a.s.l. The hydrological dam and downstream river safety were invariant for Sc.1 and Sc.2. However, higher MCLs were obtained for Sc.2, which increased the YR = 90% of the system. In the case of YR = 90% = Da = 545 hm 3 , the differences between Sc.2 and Sc.1 for MWRLTR = 1000y and MOTR = 500y were 0.2 m.a.s.l and 86 m 3 /s, respectively. Moreover, the MCL of season two were lower when the variable initial reservoir level was considered. This is because of the increase of MCLs in the other seasons. Despite this, the demand supplied with a reliability of 90% was the same, accounting for lower values of MWRLTR=1000y and MOTR = 500y if variable initial level was considered. Within the framework of the present study, accounting for the variability of initial reservoir level implied a reduction of flood control volumes (increase of MCLs) maintaining the risk of overtopping and downstream river safety. Consequently, the volume for satisfying the demands increased, providing a more reliable system in terms of regular operation.

Conclusions
The main conclusions extracted from this study are as follow: Within the framework of the present study, accounting for the variability of initial reservoir level implied a reduction of flood control volumes (increase of MCLs) maintaining the risk of overtopping and downstream river safety. Consequently, the volume for satisfying the demands increased, providing a more reliable system in terms of regular operation.

Conclusions
The main conclusions extracted from this study are as follow: • The use of a stochastic methodology allowed us to assess hydrological dam safety and downstream safety by obtaining the frequency curves of outflow and maximum reservoir water levels, while accounting for the variability in hydrological loads with respect to deterministic procedures. As a drawback, it implied a more complex procedure and computational effort.

•
We proposed a set of 98 non-inferior solutions while considering the initial reservoir level equal to the MCL for each season and 135 possible configurations while considering variable initial reservoir level. From the proposed configurations, dam stakeholders are able to decide which configuration to use depending on whether their preference is to increase dam and downstream hydrological safety or to increase water supply (with a specific reliability) in the water resources system. In the Riaño case study, the presented procedure showed improvements in the regular operation that satisfied an increase of up to 10.1% of the current annual demand of 545 hm 3 (with a reliability of 90%) while maintaining the same level of hydrological dam safety.

•
Accounting for initial reservoir variability resulted in the possibility of supplying an extra demand of 13 hm 3 (2.4% of the current annual demand) compared to the optimal solution without accounting for initial reservoir level variability.

•
The proposed stochastic procedure can improve the results obtained by deterministic procedures, increasing supply up to 26.8% of the current annual demand, from the worst regular operation configuration (not accounting for initial reservoir level variability) to the optimal configuration (accounting for initial reservoir level variability) of MCLs.