**System Dynamics Modeling for Supporting Drought-Oriented Management of the Jucar River System, Spain**

**Adria Rubio-Martin \* , Manuel Pulido-Velazquez \* , Hector Macian-Sorribes and Alberto Garcia-Prats**

Research Institute of Water and Environmental Engineering (IIAMA), Universitat Politècnica de València, 46022 Valencia, Spain; hecmasor@upv.es (H.M.-S.); agprats@upv.es (A.G.-P.)

**\*** Correspondence: adrumar@upv.es (A.R.-M.); mapuve@hma.upv.es (M.P.-V.)

Received: 25 March 2020; Accepted: 10 May 2020; Published: 15 May 2020

**Abstract:** The management of water in systems where the balance between resources and demands is already precarious can pose a challenge and it can be easily disrupted by drought episodes. Anticipated drought management has proved to be one of the main strategies to reduce their impact. Drought economic, environmental, and social impacts affect different sectors that are often interconnected. There is a need for water management models able to acknowledge the complex interactions between multiple sectors, activities, and variables to study the response of water resource systems to drought management strategies. System dynamics (SD) is a modeling methodology that facilitates the analysis of interactions and feedbacks within and between sectors. Although SD has been applied for water resource management, there is a lack of SD models able to regulate complex water resource systems on a monthly time scale and considering multiple reservoir operating rules, demands, and policies. In this paper, we present an SD model for the strategic planning of drought management in the Jucar River system, incorporating dynamic reservoir operating rules, policies, and drought management strategies triggered by a system state index. The DSS combines features from early warning and information systems, allowing for the simulation of drought strategies, evaluating their economic impact, and exploring new management options in the same environment. The results for the historical period show that drought early management can be beneficial for the performance of the system, monitoring the current state of the system, and activating drought management measures results in a substantial reduction of the economic impact of droughts.

**Keywords:** water management; resources; system dynamics; drought management; drought impacts

#### **1. Introduction**

Drought is a natural hazard and, as such, has to be understood as a natural feature of climate. Whether or not a drought becomes a disaster depends on its social, economic, and environmental impacts [1]. Therefore, the key to understanding drought is to acknowledge its different dimensions. Drought affects both surface and groundwater resources and can lead to reduced water supply for in-home consumption and agricultural and industrial activities. Furthermore, it can deteriorate water quality by rising nitrate, ammonium, and phosphate concentrations, and disturb riparian habitats [2,3]. Agriculture is the most affected sector by droughts, but many other sectors may suffer relevant losses, including energy production, tourism and recreation, transportation, urban water supply, and the environment. Sustained drought can cause social, economic, and energy crises, even leading to migration from affected zones (often rural and agricultural-focused) to other regions or nearby countries [4]. Drought is not the only issue that water resource systems have to face regarding water

availability. Water scarcity refers to continued unsustainable use of water resources and it can be influenced by water management [5]. Increasing water demand due to population growth and the development of the agricultural, energy, and industrial sectors has increased the frequency of water scarcity events that occur when there is a lack of freshwater to meet the demands [6]. Climate change is expected to further aggravate water scarcity because of the increase in drought frequency, severity, and duration [7,8].

There is an increasing concern worldwide about the ineffectiveness of most common drought management practices, largely based on crisis management and on treating symptoms (impacts) rather than the underlying causes associated with them [7]. The European Union has promoted the move from crisis management to drought risk management since 2007 [9]. However, there are gaps in the current water scarcity and droughts policy of the EU, including [10]: conceptual gaps on the understanding of causal relationships between drivers, pressures, status, and impacts; limited data on current and future water demand and availability; policy, governance and implementation gaps regarding measures to increase water supply and to target pressures and impacts caused by droughts.

Drought management plans are tools that aim to reduce the impact of droughts in water resource systems providing a framework for proactive, risk-based management [9]. A coordinated drought plan includes monitoring, early warning and information systems, impact assessment procedures, risk management measures, preparedness plans, and emergency response programs. Without these plans, nations will continue responding drought in a reactive, crisis management mode [7]. A key feature of drought management plans is the use of indices to establish a link between the state of the river basin and the measures to be taken [11]. Drought indices have been developed for assessing drought parameters including intensity, duration, severity, and spatial extent, and are effective tools in the monitoring and management of droughts [2,12]. However, traditional drought indexes often fail at detecting critical events in highly regulated systems, where natural water availability is conditioned by the operation of water infrastructures such as dams, diversions, and pumping wells. Here, ad hoc index formulations are usually adopted based on empirical combinations of several significant hydro-meteorological variables through customized formulations [13]. A system of drought indicators based on levels or thresholds depending upon the degree of water scarcity, and several management actions aiming to mitigate critical situations have been developed in the Jucar River system [11,14]. The creation and institutionalization of multi-sector partnerships have reinforced the development of efficient drought management [15]. To support drought management, scientific approaches including drought characterization, development of risk indicators, and the analysis of economic instruments for risk mitigation are involved in conjunction with the identification, selection, and prioritization of measures to lessen the effects of drought [16]. Decision support systems (DSS) have been developed to study effective drought management strategies, as they are considered one of the most effective tools for integrated water resource management [6]. The use of DSS tools for drought risk management has been increasing during the last decades [17–20]. Studying resource allocation requires the development of DSS able to apply drought management strategies and to dynamically evaluate the status of water resource systems [12]. Multi-criteria decision analysis tools (MCDA) are also oriented to assist the decision-making process in the operation of water resource systems. Nevertheless, a major problem in developing MCDA processes is to understand the risk associated with persistent drought conditions, as risk management involves subjective considerations [6]. The water sector's importance for other sectors requires policies and management strategies that are aware of the potential widespread impacts [21]. Very often, undesired effects can be derived from the execution of drought management strategies. For example, increased groundwater extraction to compensate for the reduction of surface water availability can lower base flows of rivers and streams, and reduce the piezometric level of aquifers [22]. These unexpected consequences can affect river biota, agriculture income, and urban supply in ways that are more damaging or long-lasting in time than the aforementioned drought. Consequently, there is a need for management models able to simulate the complex interactions between different sectors and activities to study the response of water resource systems to drought management strategies.

System dynamics is a theory of system structure and a set of tools for representing complex systems and analyzing their dynamic behavior [23]. This methodology is particularly useful for studying complex water resource systems with interacting elements and policies, whose behavior cannot be easily predicted [24]. The development of system dynamics models to analyze and improve water resource management has a tradition that dates back to the late 1960s. Since then, and thanks to the development of computer technology and user-friendly system dynamics software, all types of qualitative models have been developed for improving system understanding in water resource systems. However, system dynamics have not been yet applied to highly regulated and complex water resource systems for testing drought management strategies with a quantitative approach and integrating a drought early warning system.

The objective of this paper is to develop a decision support system (DSS) based on system dynamics for the efficient drought management of the Jucar River system. The DSS simulates the management of the Jucar multi-reservoir system integrating monthly-defined reservoir operating rules, stream-aquifer interaction and conjunctive use of surface and groundwater, drought management measures (linked to a system state index), and all this taking into account current water demands and allocation criteria. The tool allows studying the effect of policy and management measures in the system, and it serves as a steppingstone towards the understanding of water resource systems as a holistic system. The DSS provides quantitative results comparable to the historical records for the calibration and validation period. The calibrated model facilitates the design, testing, and selection of new drought management strategies. Section 2 introduces the system dynamics modeling method, details some applications of the methodology for the management of water resource systems and describes the Jucar River system case study. Section 2 also introduces and describes the main features of the system dynamics model developed for the case of study. Section 3 shows and discusses the results, first validating the behavior of the model and later discussing the hydrological and economic results for the simulated scenarios. Finally, Section 4 exposes the conclusions.

#### **2. Materials and Methods**

#### *2.1. System Dynamics for Water Resource Systems Modeling*

System dynamics modeling is a methodology of model development that facilities a holistic understanding of water resource systems, as it allows analyzing how different elements of a system relate to one another and permits studying the changing relations within the system when different decisions are included [25,26]. The usual purpose of the analysis of system dynamics is to understand how and why the dynamics of concern are generated and to look for managerial policies that can improve the system performance [27]. In system dynamics, the system structure is determined by the positive and negative relationships between variables, feedback loops, system archetypes, and delays [28,29]. The totality of the relationships between the system components constitutes the system structure, and the system's structure defines its behavior [30]. This methodology focuses on understanding how the physical processes, information flows, and managerial policies interact to create the dynamics of the different variables of interest [31]. To achieve this knowledge, qualitative/conceptual and qualitative/numerical modeling methods are applied.

Qualitative modeling (e.g., causal loops diagrams and definition of the positive and negative relationships between variables) improves our conceptual system understanding [29]. This type of modeling is often seen as a propaedeutic step to quantitative modeling, where the behavior of the system and the effects of different intervention policies can be visualized through simulation. Qualitative models can be further developed into quantitative models (Figure 1). This change requires a deep knowledge of the existing physical, analytical, and statistical relationships between the variables of the system. In system dynamics, the relationships between variables can be expressed by linear, non-linear mathematical equations and logical expressions such as IF-THEN statements, to introduce

issue.

management policies and rules. To assess the truthfulness of the quantitative models they are validated by comparing their results to the available historical records. introduce management policies and rules. To assess the truthfulness of the quantitative models they are validated by comparing their results to the available historical records.

*Water* **2020**, *12*, x FOR PEER REVIEW 4 of 20

**Figure 1.** System dynamics modeling framework. **Figure 1.** System dynamics modeling framework.

Traditionally, water resources management models were designed with a one-dimensional optimal engineering approach, performed with little regard for social, environmental, or cultural aspects [32]. However, the increased recognition of complexity and uncertainty has promoted the use of more flexible simulation-based tools such as the ones provided by system dynamics [28]. System dynamics provides tools for the graphical representation of systems, facilitates flexible and transparent modeling, eases the holistic understanding of the problem, captures long-run behavioral patterns and trends, facilitates clear communication of model structure and results, promotes sharing modeling, facilitates sensitivity analysis, and it is suitable for policy assessment and selection [25]. System dynamics modeling environments include Powersim (Powersim Corp., 1993), Simile (Simulistics, 2002), Stella (High Performance Systems, 1992), and Vensim (Ventana Systems, 1996). Nowadays, these environments are able to assist modelers and can handle many variables, delays, and interdependent subsystems, allowing the creation of modular object-oriented models, therefore increasing interchangeability and reusability. Traditionally, water resources management models were designed with a one-dimensional optimal engineering approach, performed with little regard for social, environmental, or cultural aspects [32]. However, the increased recognition of complexity and uncertainty has promoted the use of more flexible simulation-based tools such as the ones provided by system dynamics [28]. System dynamics provides tools for the graphical representation of systems, facilitates flexible and transparent modeling, eases the holistic understanding of the problem, captures long-run behavioral patterns and trends, facilitates clear communication of model structure and results, promotes sharing modeling, facilitates sensitivity analysis, and it is suitable for policy assessment and selection [25]. System dynamics modeling environments include Powersim (Powersim Corp., 1993), Simile (Simulistics, 2002), Stella (High Performance Systems, 1992), and Vensim (Ventana Systems, 1996). Nowadays, these environments are able to assist modelers and can handle many variables, delays, and interdependent subsystems, allowing the creation of modular object-oriented models, therefore increasing interchangeability and reusability.

The application of system dynamics in water resource management has grown since the 90s. Nowadays, we find applications of system dynamics modeling to study a large variety of water resource issues [29]. They range from region-scale models with multiple demands and frequent water scarcity events [33,34], to models coupling surface and groundwater dynamics for a basin [35], flood management or predicting models [36,37], reservoir operation and water supply for multiple water users [38], and the design of water pricing policies [39]. However, system dynamics application to simulate the management of highly regulated water resource systems integrating multiple reservoirs, operating rules, dynamic drought management, groundwater use, and conflicting water demands remains very limited. Yet all these features are required to analyze the issue of drought early warning and management in complex water resource systems. The application of system dynamics in water resource management has grown since the 90s. Nowadays, we find applications of system dynamics modeling to study a large variety of water resource issues [29]. They range from region-scale models with multiple demands and frequent water scarcity events [33,34], to models coupling surface and groundwater dynamics for a basin [35], flood management or predicting models [36,37], reservoir operation and water supply for multiple water users [38], and the design of water pricing policies [39]. However, system dynamics application to simulate the management of highly regulated water resource systems integrating multiple reservoirs, operating rules, dynamic drought management, groundwater use, and conflicting water demands remains very limited. Yet all these features are required to analyze the issue of drought early warning and management in complex water resource systems.

Drought management is a multidimensional concept that includes meteorological, ecological, hydrological, environmental, and socioeconomic perspectives. The development of DSS for improving drought management requires the combination of several models [6]. Coupling and analyzing the interactions between these models is often a difficult issue. System dynamics is a Drought management is a multidimensional concept that includes meteorological, ecological, hydrological, environmental, and socioeconomic perspectives. The development of DSS for improving drought management requires the combination of several models [6]. Coupling and analyzing the interactions between these models is often a difficult issue. System dynamics is a methodology that

methodology that provides a common playground for the interaction of different subsystems and submodels, facilitating the analysis of the existing relationships and providing a holistic view of the provides a common playground for the interaction of different subsystems and submodels, facilitating the analysis of the existing relationships and providing a holistic view of the issue. *Water* **2020**, *12*, x FOR PEER REVIEW 5 of 20

#### *2.2. Case Study: Drought Management in the Jucar River System 2.2. Case Study: Drought Management in the Jucar River System*

The Jucar River system is located in Easter Spain. The system is subjected to a tight equilibrium between total water demand (1505 Mm<sup>3</sup> /year, 2009–2015 period average) and water resource availability (1548 Mm<sup>3</sup> /year) [40]. Agriculture is the largest water use by far (89%), followed by urban (9%) and industrial uses (2%). The Jucar is the main source of urban water supply to the city of Valencia and its metropolitan area (about 1,500,000 inhabitants, third largest municipality in Spain). Water from the Jucar is diverted to the Turia River through a 60 km canal (Canal Jucar-Turia), also used for irrigation of mainly citrus and vegetables (Figure 2). Furthermore, there is an intense water use for irrigation in the lower Jucar, downstream of Tous reservoir, with traditional irrigation districts holding senior water rights dating back to the Middle Ages. Non-consumptive water demands include minimum ecological instream flows and hydropower generation. The Jucar River system is located in Easter Spain. The system is subjected to a tight equilibrium between total water demand (1505 Mm<sup>3</sup> /year, 2009–2015 period average) and water resource availability (1548 Mm³/year) [40]. Agriculture is the largest water use by far (89%), followed by urban (9%) and industrial uses (2%). The Jucar is the main source of urban water supply to the city of Valencia and its metropolitan area (about 1,500,000 inhabitants, third largest municipality in Spain). Water from the Jucar is diverted to the Turia River through a 60 km canal (Canal Jucar-Turia), also used for irrigation of mainly citrus and vegetables (Figure 2). Furthermore, there is an intense water use for irrigation in the lower Jucar, downstream of Tous reservoir, with traditional irrigation districts holding senior water rights dating back to the Middle Ages. Non-consumptive water demands include minimum ecological instream flows and hydropower generation.

**Figure 2.** Main features of the Jucar River system included in the model. **Figure 2.** Main features of the Jucar River system included in the model.

The main surface reservoirs are Alarcon (1112 Mm<sup>3</sup> of capacity), Contreras (463 Mm<sup>3</sup> of useful capacity), and Tous (378 Mm<sup>3</sup> ). The regulation capacity of these reservoirs is mainly multi-annual: Alarcon and Contreras are devoted to consumptive uses, while Tous is mostly used for flood protection. The intense overexploitation of the main groundwater body, the Mancha Oriental aquifer (middle basin, near Albacete), for irrigation since the 1970s has shifted the stream-aquifer interaction between Alarcon and Tous from gaining to losing river, diminishing downstream surface water availability. The sustainable use of this aquifer is one of the challenges in the management of the system [41,42]. During droughts, the Plana de Valencia Sur aquifer, located in the lower basin (downstream of Tous), is used as an alternative water source. The main surface reservoirs are Alarcon (1112 Mm<sup>3</sup> of capacity), Contreras (463 Mm<sup>3</sup> of useful capacity), and Tous (378 Mm<sup>3</sup> ). The regulation capacity of these reservoirs is mainly multi-annual: Alarcon and Contreras are devoted to consumptive uses, while Tous is mostly used for flood protection. The intense overexploitation of the main groundwater body, the Mancha Oriental aquifer (middle basin, near Albacete), for irrigation since the 1970s has shifted the stream-aquifer interaction between Alarcon and Tous from gaining to losing river, diminishing downstream surface water availability. The sustainable use of this aquifer is one of the challenges in the management of the system [41,42]. During droughts, the Plana de Valencia Sur aquifer, located in the lower basin (downstream of Tous), is used as an alternative water source.

Water scarcity, irregular hydrology, and groundwater overdraft result in droughts with significant economic, social, and environmental consequences. This situation is expected to be exacerbated by the impacts of climate and socioeconomic (global) changes and increasing institutional impediments from political disputes among the two main riparian regions, Castilla-La Mancha (upstream; mainly Albacete province) and Valencia (middle and downstream basin). A range of different innovative solutions are considered to face the main water management issues, such as pumping-water right acquisitions during droughts, increasing wastewater reuse, "in lieu" Water scarcity, irregular hydrology, and groundwater overdraft result in droughts with significant economic, social, and environmental consequences. This situation is expected to be exacerbated by the impacts of climate and socioeconomic (global) changes and increasing institutional impediments from political disputes among the two main riparian regions, Castilla-La Mancha (upstream; mainly Albacete province) and Valencia (middle and downstream basin). A range of different innovative solutions are considered to face the main water management issues, such as pumping-water right acquisitions during droughts, increasing wastewater reuse, "in lieu" recharge (providing surplus

recharge (providing surplus surface water to groundwater users, keeping groundwater in storage for

surface water to groundwater users, keeping groundwater in storage for later use), water-saving in agriculture through drip irrigation, new water allocation mechanisms, water banks, water pricing, and irrigated crop drought insurances (among others), which makes this case a real lab for analyzing risk management strategies to cope with drought, extreme events, and climate change [43].

The operation of the system, managed by the Jucar River Basin Authority (Confederacion Hidrografica del Jucar, CHJ), is subject to physical, environmental, and legal constraints. The main physical constraints correspond to the reservoir, river, and canal capacities. The environmental constraints are the minimum flows prescribed in certain river reaches and the inflow requirements of the Albufera wetland. The main legal constraint in the Jucar River system is the Alarcon Agreement, signed between the Spanish Ministry of the Environment and the senior users of the Jucar River—mainly farmers—gathered together in the Unidad Sindical de Usuarios del Jucar (USUJ). The agreement divides Alarcon in two zones by a rule curve. If the water level in Alarcon is above the threshold, water can be freely allocated, but if the storage is below certain value, water in the system is reserved exclusively for the USUJ members. In this case, other water users who want to access water from the Jucar River would have to pay a financial compensation to USUJ. The operators also follow additional criteria to decide the releases during the irrigation period (May–September): not causing undesired spills from Tous (the downstream reservoir), not storing more than 450 Mm<sup>3</sup> in Contreras to avoid stability problems, and not storing more than 72 Mm<sup>3</sup> in Tous at the end of the summer to avoid flood damage during autumn due to intense rainfall events [42].

The Jucar River basin, as most Mediterranean and south-eastern basins of Spain, is very vulnerable to droughts [11]. The recurrence of these events is also an important factor when considering the management of the system, as a high-frequency appearance of droughts do not allow the system to properly recover water storage to face future water-scarcity events. The latest drought periods (1991–1995; 1997–2000 and 2004–2009) were classified as extreme droughts using the SPI index [14]. During the drought period 2005–2008, surface water available for agriculture decreased by up to 40% compared to the average. Because of this, drought emergency wells in the lower basin were activated. Despite of these efforts, the drought caused an important economic impact, especially to agriculture activities. The situation is expected to be exacerbated by the impact of climate change [8,44].

A key feature of drought management plans is the indices that define the different drought stages and trigger mitigation measures. Drought indices should capture the state of the water resource system as a whole, allowing the planner to active measures to reduce its impact. Some of the measures for drought management include conjunctive use of surface and groundwater, awareness campaigns to promote domestic water savings, economic tools, control of the supply to agricultural demands from reservoirs, and water reuse [16]. Traditionally, the management of droughts in the Jucar Basin was regulated as an emergency, and the application of Royal Decrees was necessary to mitigate their impacts [12]. From 2007, drought management in the Jucar River system is regulated by a drought management plan [14,40] that establishes a state index to monitor the system and a set of drought management measures triggered by the different drought stages. This index is calculated using different variables distributed in the area of the river basin, including reservoir storages, groundwater levels, streamflow, precipitation, and reservoir inflows. The state index takes values between 0 and 1, with four system states: normal, pre-alert, alert, and emergency. Then, different drought management measures are applied depending on the system's state index. These measures can be divided into 2 groups, (1) control of water supply for urban and agricultural uses and, (2) increase of water availability by drought emergency wells use and increasing water reuse.

#### *2.3. System Dynamics for the Jucar River System*

The system dynamics model developed for the Jucar River system represents its current management with a monthly time step, including the state index of the system and the management measures linked to this state. The software Vensim Pro [45] has been used for the creation of the model. The Jucar model was divided into 5 subsystems:


based on it.

5. State index: calculates the state index and defines the management measures to take based on it. 5. State index: calculates the state index and defines the management measures to take

**Figure 3.** Subsystem for the general view of the system. **Figure 3.** Subsystem for the general view of the system.

The model incorporates monthly water inflows in 5 sub-basins where data from CEDEX [46], the Spanish institution responsible for collecting and supplying data on civil engineering and water, has been obtained and processed. The main view of the model (Figure 3) captures the water flows through the Jucar system, including water infrastructures and stream-aquifer interaction with the Mancha Oriental aquifer. This structure is based on previous models for the area [42], and provides a general framework to visualize the system's network and to allow the integration of other submodels. The model incorporates a submodel that simulates the current operation of the system (Figure 4), based on historical records and trends of the main variables. The model incorporates monthly water inflows in 5 sub-basins where data from CEDEX [46], the Spanish institution responsible for collecting and supplying data on civil engineering and water, has been obtained and processed. The main view of the model (Figure 3) captures the water flows through the Jucar system, including water infrastructures and stream-aquifer interaction with the Mancha Oriental aquifer. This structure is based on previous models for the area [42], and provides a general framework to visualize the system's network and to allow the integration of other sub-models. The model incorporates a submodel that simulates the current operation of the system (Figure 4), based on historical records and trends of the main variables.

The operating rules of the three reservoirs are defined at the monthly scale, mimicking the operation of the system for the 2003–2013 period, and introducing the constraints that bind the seasonal operation of the Jucar River system. The rules were obtained using fuzzy rule-based systems (FRB), co-developed with the experts from the Operation Office of the Jucar River Basin Authority [42]. A series of workshops and surveys were used to extract the decision-making processes followed in the seasonal operation of the Jucar River system. The implicit operation of the system was encoded into two FRB systems that were validated against historical records on reservoir storages and releases, streamflows, and deliveries to consumptive demands for the 2003–2013 period. The developed FRB were introduced into the SD model through piecewise linear regressions equations (Figure 5). Some flexibility is lost in the process of transforming the FRB rules into linear regressions, as it can be observed in the figure.

reservoirs.

(Figure 4), based on historical records and trends of the main variables.

**Figure 3.** Subsystem for the general view of the system.

The model incorporates monthly water inflows in 5 sub-basins where data from CEDEX [46], the Spanish institution responsible for collecting and supplying data on civil engineering and water, has been obtained and processed. The main view of the model (Figure 3) captures the water flows through the Jucar system, including water infrastructures and stream-aquifer interaction with the Mancha Oriental aquifer. This structure is based on previous models for the area [42], and provides a general framework to visualize the system's network and to allow the integration of other sub-

*Water* **2020**, *12*, x FOR PEER REVIEW 7 of 20

4. Reservoir operation: defines the seasonal operating rules of the system.

management models [40,42].

based on it.

water, and the system deliveries and deficits.

2. Mancha Oriental aquifer: simulates the aquifer using a two-cell embedded multireservoir model, in line with the one used by the CHJ in its water resource

3. Water demands: defines the different monthly water demands, the distribution of

5. State index: calculates the state index and defines the management measures to take

**Figure 4.** Reservoir operating rules subsystem that incorporates the monthly operating rules for each reservoir, variables, and seasonal parameters that determine final releases. (Figure 5). Some flexibility is lost in the process of transforming the FRB rules into linear regressions, as it can be observed in the figure.

**Figure 5.** Graphical representation of the simulated operating rule of Tous reservoir in July. Blue dots represent the values of releases using fuzzy logic. Red crosses show values for the piecewise linear regression introduced into the SD model. **Figure 5.** Graphical representation of the simulated operating rule of Tous reservoir in July. Blue dots represent the values of releases using fuzzy logic. Red crosses show values for the piecewise linear regression introduced into the SD model.

To compensate this loss of flexibility, the obtained rules for Alarcon and Contreras reservoirs were adjusted using seasonal factors (depending on whether it was or not irrigation season) and a scarcity factor different for both winter and summer seasons, to account for differences observed in the management of the system that were not correctly captured by the calculated piece-wise linear equations. Releases from Tous were computed as the minimum value between the downstream demand and the releases calculated by the piecewise linear equations. This implies that the system will not release more water from Tous than needed, minimizing unwanted releases to the sea while still capturing the seasonal behavior provided by the operating rules. The Alarcon Agreement was explicitly introduced into the model's formulation. To compensate this loss of flexibility, the obtained rules for Alarcon and Contreras reservoirs were adjusted using seasonal factors (depending on whether it was or not irrigation season) and a scarcity factor different for both winter and summer seasons, to account for differences observed in the management of the system that were not correctly captured by the calculated piece-wise linear equations. Releases from Tous were computed as the minimum value between the downstream demand and the releases calculated by the piecewise linear equations. This implies that the system will not release more water from Tous than needed, minimizing unwanted releases to the sea while still capturing the seasonal behavior provided by the operating rules. The Alarcon Agreement was explicitly introduced into the model's formulation.

The water demands considered by the model are divided into urban and agricultural demands and were located and compiled from the public information provided by the CHJ [47]. Most of them The water demands considered by the model are divided into urban and agricultural demands and were located and compiled from the public information provided by the CHJ [47]. Most of them

are located downstream Tous, although the model also accounts for the demands located in the

are located downstream Tous, although the model also accounts for the demands located in the middle basin, one of them being a groundwater demand that affects the stream-aquifer interaction. The current operating rules of the system prioritizes water allocation to urban uses. Environmental requirements have been considered as a restriction and are captured by the operating rules of the reservoirs. *Water* **2020**, *12*, x FOR PEER REVIEW 9 of 20

The model simulates stream-aquifer interaction between the Mancha Oriental aquifer and the Jucar River using a two-cell Embedded Multi-reservoir Model (Figure 6) [48]. Its formulation is based on the analytical solution of the stream-aquifer flow equation applied to linear systems, as well as it analogy with the state equation. Groundwater discharge can be expressed as the theoretical sum of an infinite number of linear reservoirs whose discharge is linearly proportional to the stored volume. In normal conditions, a limited number of linear reservoirs is enough to adequately reproduce groundwater discharge. Although the EMM does not calculate spatially-distributed heads and internal groundwater flows, it can provide an accurate representation of stream-aquifer interactions, even in karstic aquifers [49,50] and it is used in some general DSS services for water resource management [42,51]. Groundwater flow is calculated as the integration of the outflow of 2 linear reservoirs in which the discharge is linearly proportional to the volume stored. The EMM built for the Mancha Oriental aquifer represents exclusively the impacts of the anthropic stresses on stream-aquifer interaction, since the natural discharge was already included in the natural inflow time series of the model [42]. The anthropic-induced net recharge corresponds to the agricultural percolation minus groundwater abstractions. As shown by Macian-Sorribes et al., 2017, the calibrated EMM was able to capture well both the over-year trend and the seasonal variation of the historical values. The model simulates stream-aquifer interaction between the Mancha Oriental aquifer and the Jucar River using a two-cell Embedded Multi-reservoir Model (Figure 6) [48]. Its formulation is based on the analytical solution of the stream-aquifer flow equation applied to linear systems, as well as it analogy with the state equation. Groundwater discharge can be expressed as the theoretical sum of an infinite number of linear reservoirs whose discharge is linearly proportional to the stored volume. In normal conditions, a limited number of linear reservoirs is enough to adequately reproduce groundwater discharge. Although the EMM does not calculate spatially-distributed heads and internal groundwater flows, it can provide an accurate representation of stream-aquifer interactions, even in karstic aquifers [49,50] and it is used in some general DSS services for water resource management [42,51]. Groundwater flow is calculated as the integration of the outflow of 2 linear reservoirs in which the discharge is linearly proportional to the volume stored. The EMM built for the Mancha Oriental aquifer represents exclusively the impacts of the anthropic stresses on streamaquifer interaction, since the natural discharge was already included in the natural inflow time series of the model [42]. The anthropic-induced net recharge corresponds to the agricultural percolation minus groundwater abstractions. As shown by Macian-Sorribes et al., 2017, the calibrated EMM was able to capture well both the over-year trend and the seasonal variation of the historical values.

**Figure 6.** Subsystem for the stream aquifer interaction between the Jucar River and the Mancha Oriental aquifer. **Figure 6.** Subsystem for the stream aquifer interaction between the Jucar River and the Mancha Oriental aquifer.

The model also implements a state index subsystem. This subsystem checks the state of the system each time-step during the simulation, as does the state index used by the CHJ on a monthly basis. The equations defining the relationship between past and present system states are taken from the Jucar drought management plan [14,52]. The model also implements a state index subsystem. This subsystem checks the state of the system each time-step during the simulation, as does the state index used by the CHJ on a monthly basis. The equations defining the relationship between past and present system states are taken from the Jucar drought management plan [14,52].

The monthly system state index ( (has the following expression: The monthly system state index (*Si*) has the following expression:

Where

$$\mathbf{S}\_{l} = \ ^{S}\_{\not\equiv} \begin{bmatrix} \mathbf{1} \ \mathbf{1} \\ \mathbf{2} \end{bmatrix} \mathbf{1} \ \begin{bmatrix} 1 + \frac{V\_{l} \ \forall\_{\overline{l}} \ \forall\_{\alpha} \forall\_{\alpha\overline{\nu}}}{V\_{m} \underline{V\_{\alpha\nu\overline{\nu}}} \ \forall\_{\alpha\overline{\nu}} \ \mathbf{0} \end{bmatrix} \mathbf{1} \ \mathbf{V}\_{\not\overline{l}} \ \succeq \mathbf{M}\_{\overline{\text{fl}}\nu}$$

$$\mathbf{S}\_{l} = \begin{array}{c} V\_{l} - \ V\_{min} \\ \mathbf{Z} \{ V\_{av} - V\_{min} \} \end{array} \text{if} \quad V\_{l} < V\_{av}$$

the recorded average, maximum and minimum monthly values of the variable since 1982. In the case of the SD model, the subsystem uses historical data of the average, maximum, and minimum value

is the value of the variable at the beginning of the month *i* and , y are

$$\mathcal{S}\_{\dot{i}} = \frac{V\_{\dot{i}} - V\_{\min}}{2(V\_{av} - V\_{\min})} \text{ if } V\_{\dot{i}} < V\_{av}$$

where *V<sup>i</sup>* is the value of the variable at the beginning of the month *i* and *Vav*, *Vmax* y *Vmin* are the recorded average, maximum and minimum monthly values of the variable since 1982. In the case of the SD model, the subsystem uses historical data of the average, maximum, and minimum value of water storage for each one of the three reservoirs and compares the recorded values to the current state of the system. Although the evaluation of the system state index executed by the water authority for the Jucar River basin takes into account 9 additional variables other than the water storage (including piezometric levels and water inflows), in regulated systems the volume stored in the reservoirs is regarded as a good approximation of the actual status of the whole system [53]. *Water* **2020**, *12*, x FOR PEER REVIEW 10 of 20 of water storage for each one of the three reservoirs and compares the recorded values to the current state of the system. Although the evaluation of the system state index executed by the water authority for the Jucar River basin takes into account 9 additional variables other than the water storage (including piezometric levels and water inflows), in regulated systems the volume stored in the reservoirs is regarded as a good approximation of the actual status of the whole system [53].

The state index subsystem is able to trigger drought management measures depending on the current state of the system (Figure 7). The state index subsystem is able to trigger drought management measures depending on the current state of the system (Figure 7).

**Figure 7.** State index subsystem to calculate the system's state comparing with historical data and **Figure 7.** State index subsystem to calculate the system's state comparing with historical data and incorporating drought management strategies.

incorporating drought management strategies.

system.

**3. Results and Discussion** 

*3.1. Model Evaluation* 

The system state index takes values that range from 0 to 1. Each month, the model transforms the system state index (a floating-point number) to the corresponding integer state (normal, pre-alert, alert, and emergency) applying the thresholds defined by the water authority. Drought management strategies defined in this subsystem are introduced as actions into their respective subsystems using shadow variables. The measures implemented consider both supply and demand side solutions. For instance, when triggered, the variable "Agricultural supply management" is linked to the agricultural supply on the "Water demand, supply and deficit" subsystem applying a restriction of 20% or 40% on the deliveries to the agricultural demands, depending on the state index. "Urban supply management" restricts the water delivered to the urban demand in alert or worse situations by 5%, reproducing the estimated effect of the water-saving awareness campaigns proposed by the water authority [14,52]. "Groundwater extractions" and "Alternative water sources" variables simulate the use of wells and the reuse of wastewater respectively for agricultural supply; the intensity of both actions depends on the monthly state of the system. All the values and management measures The system state index takes values that range from 0 to 1. Each month, the model transforms the system state index (a floating-point number) to the corresponding integer state (normal, pre-alert, alert, and emergency) applying the thresholds defined by the water authority. Drought management strategies defined in this subsystem are introduced as actions into their respective subsystems using shadow variables. The measures implemented consider both supply and demand side solutions. For instance, when triggered, the variable "Agricultural supply management" is linked to the agricultural supply on the "Water demand, supply and deficit" subsystem applying a restriction of 20% or 40% on the deliveries to the agricultural demands, depending on the state index. "Urban supply management" restricts the water delivered to the urban demand in alert or worse situations by 5%, reproducing the estimated effect of the water-saving awareness campaigns proposed by the water authority [14,52]. "Groundwater extractions" and "Alternative water sources" variables simulate the use of wells and the reuse of wastewater respectively for agricultural supply; the intensity of both actions depends on the monthly state of the system. All the values and management measures represented in the state index subsystem are based on the current drought management plan for the system.

The system dynamics model of the Jucar River system was evaluated using the 2003–2013 period. The comparison between the model's results and historical records showed that the model is able to capture the observed operation (Figure 8). Residual plots for the same variables can be found

represented in the state index subsystem are based on the current drought management plan for the

### **3. Results and Discussion**

#### *3.1. Model Evaluation*

scenarios.

The system dynamics model of the Jucar River system was evaluated using the 2003–2013 period. The comparison between the model's results and historical records showed that the model is able to capture the observed operation (Figure 8). Residual plots for the same variables can be found in Appendix A (Figure A1). Total storage was closely reproduced by the model, as can be observed in the plot and in the R-squared index. The Alarcon and Contreras releases were adequately reproduced on a broader view, due to the resemblance of the intra-annual patterns. However, the model results depart from the historical observations in some years. This is due to the fact that the middle basin is modeled in less detail than the lower one. For instance, hydroelectric production has not been included in the middle basin. In any case, storages in Alarcon and Contreras are adequately reproduced (Figure A2) and the overall in-year dynamics of the system was matched, so these deviations do not have a significant impact on the performance of the model. Tous releases results correctly fit the available data. These releases have a major importance for the model since the majority of the surface water demands are located downstream. As for water supply deficits, the simulated values matched the observed data adequately, including the main peaks associated with the 2005–2008 drought, especially during the years when the drought was more severe. Differences between observed and simulation results can also be explained by the fact that the model assumes a constant annual demand for the whole simulation period while, actually, demands changed due to population change, variation in irrigated areas, and shift from gravity to drip irrigation [47]. *Water* **2020**, *12*, x FOR PEER REVIEW 11 of 20 in Appendix A (Figure A1). Total storage was closely reproduced by the model, as can be observed in the plot and in the R-squared index. The Alarcon and Contreras releases were adequately reproduced on a broader view, due to the resemblance of the intra-annual patterns. However, the model results depart from the historical observations in some years. This is due to the fact that the middle basin is modeled in less detail than the lower one. For instance, hydroelectric production has not been included in the middle basin. In any case, storages in Alarcon and Contreras are adequately reproduced (Figure A2) and the overall in-year dynamics of the system was matched, so these deviations do not have a significant impact on the performance of the model. Tous releases results correctly fit the available data. These releases have a major importance for the model since the majority of the surface water demands are located downstream. As for water supply deficits, the simulated values matched the observed data adequately, including the main peaks associated with the 2005–2008 drought, especially during the years when the drought was more severe. Differences between observed and simulation results can also be explained by the fact that the model assumes a constant annual demand for the whole simulation period while, actually, demands changed due to population change, variation in irrigated areas, and shift from gravity to drip irrigation [47].

**Figure 8.** Comparison between observed and simulated values for key variables of the system. **Figure 8.** Comparison between observed and simulated values for key variables of the system.

Releases from Alarcon and Contreras are cumbersome to model because of the uncertainties of the middle basin, changes in downstream demands, and varied criteria of releases and management over the simulated period. Although it would be possible to introduce variable demands into the model, there is no available data to represent the variation of all the demands during the simulation period. Furthermore, although the model incorporates monthly operating rules for the reservoirs based on a fuzzy logic representation of the system operation reported by the managers [42], those rules cannot reproduce discretionary changes in the operation of the system during the simulation period. Releases from Alarcon and Contreras are cumbersome to model because of the uncertainties of the middle basin, changes in downstream demands, and varied criteria of releases and management over the simulated period. Although it would be possible to introduce variable demands into the model, there is no available data to represent the variation of all the demands during the simulation period. Furthermore, although the model incorporates monthly operating rules for the reservoirs based on a fuzzy logic representation of the system operation reported by the managers [42], those rules cannot reproduce discretionary changes in the operation of the system during the simulation period.

Once verified that the developed model matches adequately the historical behavior of the Jucar River system, further simulations were launched to test different management assumptions and

*3.2. System State Index and Drought Management Strategies* 

Once verified that the developed model matches adequately the historical behavior of the Jucar River system, further simulations were launched to test different management assumptions and scenarios.

#### *3.2. System State Index and Drought Management Strategies*

The SD model has been applied to study the interaction between the previously indicated drought management strategies and other variables of the system. A comparison between simulations with and without the drought management strategies introduced into the management in 2007 was performed. Results obtained when applying the drought management measures show improvements for the state index of the system and for the system's total water storage (Figure 9). *Water* **2020**, *12*, x FOR PEER REVIEW 12 of 20 simulations with and without the drought management strategies introduced into the management in 2007 was performed. Results obtained when applying the drought management measures show improvements for the state index of the system and for the system's total water storage (Figure 9).

**Figure 9.** (**a**) State index and (**b**) total storage with and without drought management strategies. **Figure 9.** (**a**) State index and (**b**) total storage with and without drought management strategies.

The system state index benefits from applying the drought management strategies defined in the state index subsystem. Thanks to them, the system state does not drop into an emergency state during the 2005–2008 drought. It also recovers earlier from the alert state during that drought, and it enters the prealert stage months before than the scenario with no drought management measures. The system state index benefits from applying the drought management strategies defined in the state index subsystem. Thanks to them, the system state does not drop into an emergency state during the 2005–2008 drought. It also recovers earlier from the alert state during that drought, and it enters the prealert stage months before than the scenario with no drought management measures.

drought managed model, even when the drought is over (from 2010 onwards). According to the model, water storage in reservoirs is increased significantly when drought management measures are applied. The difference is up to almost 100 Mm³ during October 2008. This is the result of the 94

management strategies taken in anticipation thanks to the state index and the four threshold levels defined. The anticipated management also allows to reduce the system vulnerability by 62% in comparison with the scenario without drought management and considering vulnerability as the ratio between total water supply deficit and the number of failures to meet the demands during the whole period. A reduction in vulnerability means that the average water shortage is lower, although

After the system enters a normal state, it is worth pointing out that the state index is higher for the

After the system enters a normal state, it is worth pointing out that the state index is higher for the drought managed model, even when the drought is over (from 2010 onwards). According to the model, water storage in reservoirs is increased significantly when drought management measures are applied. The difference is up to almost 100 Mm<sup>3</sup> during October 2008. This is the result of the management strategies taken in anticipation thanks to the state index and the four threshold levels defined. The anticipated management also allows to reduce the system vulnerability by 62% in comparison with the scenario without drought management and considering vulnerability as the ratio between total water supply deficit and the number of failures to meet the demands during the whole period. A reduction in vulnerability means that the average water shortage is lower, although the frequency of these shortages may increase. These drought management measures entail the use of drought emergency wells for water abstraction within a maximum of 98 Mm<sup>3</sup> /year (Figure 10) following the plan defined by the water authority [14,52]. *Water* **2020**, *12*, x FOR PEER REVIEW 13 of 20 the frequency of these shortages may increase. These drought management measures entail the use of drought emergency wells for water abstraction within a maximum of 98 Mm<sup>3</sup> /year (Figure 10) following the plan defined by the water authority [14,52]. Water pumping from drought emergency wells located in the lower basin compensates the reduced surface water supply and alleviate the drought impact on agriculture. These groundwater abstractions are activated when the system falls into the alert state, and water abstraction scale up above 8 Mm³/month if the emergency state is reached (Figure 10).

**Figure 10.** Water abstraction from emergency wells during drought compared to the system state index. **Figure 10.** Water abstraction from emergency wells during drought compared to the system state index.

*3.3. Economic Impact of Droughts*  Results show that the total reservoir storage of the basin improves when drought management measures are applied. It is to expect that the gained storage will benefit the early recovery of the system allowing for more regular deliveries to agricultural demands. Indeed, it is possible to calculate the Water pumping from drought emergency wells located in the lower basin compensates the reduced surface water supply and alleviate the drought impact on agriculture. These groundwater abstractions are activated when the system falls into the alert state, and water abstraction scale up above 8 Mm<sup>3</sup> /month if the emergency state is reached (Figure 10).

#### economic losses associated with the mismanagement of droughts for the 2003–2009 period (Figure 11). *3.3. Economic Impact of Droughts*

the 2000's drought.

Results show that the total reservoir storage of the basin improves when drought management measures are applied. It is to expect that the gained storage will benefit the early recovery of the system allowing for more regular deliveries to agricultural demands. Indeed, it is possible to calculate the economic losses associated with the mismanagement of droughts for the 2003–2009 period (Figure 11).

Economic losses were calculated by economically characterizing the monthly demands of the system defined as targets [47] using demand curves or functions obtained by Positive Mathematical Programming (PMP) [54] for the different agricultural demands [55]. Benefits were obtained as the integration of the demand function between zero and the level of supply. It can be observed (Figure 11) that economic losses concentrate on the drought period (2005–2008), particularly when the system state index stays in alert for several months (2006–2008). During the irrigation season in drought periods is when economic losses rise due to water scarcity. The fact that, as defined by the drought management strategies subsystem, in alert and emergency states the water supply for agriculture is reduced by up to 40% its original demand could be thought of as detrimental for agricultural interests. However, according to the simulations, the water saved helps a faster recovery of the system, guarantees urban water supply, and reduces the long-term impact of droughts. In the model, the economic impact of the

**Figure 11.** Estimation of economic losses for agriculture compared to the system state index during

Economic losses were calculated by economically characterizing the monthly demands of the system defined as targets [47] using demand curves or functions obtained by Positive Mathematical Programming (PMP) [54] for the different agricultural demands [55]. Benefits were obtained as the

2005–2008 drought was reduced from 89 M€ to 29 M€ thanks to the drought management strategies implemented. Due to conjunctive use of superficial and groundwater, agricultural activities suffer lower impact even considering the significant restrictions they suffer during the alert and emergency states. When the amount of available water is scarce, using groundwater to supply crops under deficit irrigation guarantees the survival of the plantations and minimizes economic losses. *3.3. Economic Impact of Droughts*  Results show that the total reservoir storage of the basin improves when drought management measures are applied. It is to expect that the gained storage will benefit the early recovery of the system allowing for more regular deliveries to agricultural demands. Indeed, it is possible to calculate the economic losses associated with the mismanagement of droughts for the 2003–2009 period (Figure 11).

*Water* **2020**, *12*, x FOR PEER REVIEW 13 of 20

the frequency of these shortages may increase. These drought management measures entail the use

Water pumping from drought emergency wells located in the lower basin compensates the reduced surface water supply and alleviate the drought impact on agriculture. These groundwater abstractions are activated when the system falls into the alert state, and water abstraction scale up

/year (Figure 10)

of drought emergency wells for water abstraction within a maximum of 98 Mm<sup>3</sup>

following the plan defined by the water authority [14,52].

above 8 Mm³/month if the emergency state is reached (Figure 10).

**Figure 11.** Estimation of economic losses for agriculture compared to the system state index during the 2000's drought. **Figure 11.** Estimation of economic losses for agriculture compared to the system state index during the 2000's drought.

#### Economic losses were calculated by economically characterizing the monthly demands of the **4. Conclusions**

system defined as targets [47] using demand curves or functions obtained by Positive Mathematical Programming (PMP) [54] for the different agricultural demands [55]. Benefits were obtained as the This paper presents a system dynamics DSS for drought management of the Jucar River system, taking into account the combination of a state index and several drought management strategies. The resulting DSS showed the potential of system dynamics for simulating the management of multi-reservoir systems, integrating monthly-defined operating rules for the reservoirs, stream-aquifer interaction, conflicting water demands, and drought management strategies. The model adequately reproduces the operation of the system and is able to produce accurate quantitative results, as shown by the comparison with the historical records.

The DSS takes advantage of the holistic concept that drives the methodology and incorporates components from different disciplines (hydrology, economics, social sciences, laws, etc.) into its modular structure. The state index subsystem is an example of how it is possible to integrate policies and management strategies into a water resource model using a system dynamics approach. Likewise, water policy or legislation has been incorporated into the model—e.g., the Alarcon agreement.

The DSS opens up the possibility of analyzing different drought management strategies and assessing the interactions, feedbacks, and impacts within and between multiple sectors and variables.

Results showed that drought management strategies have a net positive effect in the Jucar River system from both the economic (agriculture) and the water management perspective. The defined measures lowered agricultural losses for the 2005–2008 drought period and increased the amount of stored water during drought allowing the faster recovery of the system. Although the model provides quantitative results similar to the historical data available, the main goal of a system dynamics model is neither to forecast nor to optimize, but studying patterns, trends, and interactions between different variables of the model [24]. Modeling and dynamically simulating the change in water resources over time provides a scientifically defensible basis for proactive management strategies, enhancing our prospects to maximize the adaptive capacity of the system as a whole [29].

Moreover, the same methodology used to study drought management strategies can be applied to study the impact of different realities and inputs into the system. The DSS model developed for the Jucar River system uses a quantitative approach for its simulation. Consequently, it requires numeric data and well-tuned equations to capture the behavior of the system in detail. Qualitative variables and inputs can also be implemented in this kind of model. Qualitative modeling often introduces "soft" variables to study the general patterns of behavior of the model, rather than precise numbers [56]. In this case, qualitative modeling can be restricted to new subsystems for the testing of different non-easily quantifiable hypothesis.

The model herein presented was successfully developed for the Jucar case study and it could be replicated in any basin or system where enough information and data are available. The development of quantitative system dynamics models requires the use of a large volume of data coming from different fields (from hydrological to economic and reservoir data) as well as a deep understanding of the system structure and behavior. Very often, the most complex issue of this type of model is the development of the monthly operating rules for the reservoirs. In this case, the final rules were inferred using fuzzy logic, but additional tests showed that it is possible to simulate the operation of the system using other approaches and calibrating the rules with the historical records for the releases and water storage of the reservoirs. Although the model is able to reproduce the stream aquifer interaction between the Jucar River and the Mancha Oriental aquifer, it simulates neither groundwater heads nor aquifer storage. Groundwater head specifically is a determinant factor for the Mancha Oriental aquifer, as it has suffered continuous drops in groundwater levels due to intense pumping since the early 1970s. To assess the effect of drought policies on groundwater levels, it would be necessary to apply a detailed groundwater model, such as finite-difference model, coupling it with the system dynamics model either through scripting, wrapping, or spreadsheet coupling [57].

The model developed using system dynamics for the Jucar River system has the potential to grow and increase its scope by integrating new dynamics that can modify the behavior of the whole system. Future lines of work include linking the agricultural demand subsystem and a land-use subsystem, which would allow for introducing changes in agricultural land use based on economic benefit from previous years and on changes in land-use policies. System dynamics provides an excellent framework to study trade-offs that land use changes can introduce in specific sectors and communities [58]. Furthermore, it is already possible to activate population growths or losses over time to study how changes in urban demand can affect the system. These functionalities are required to test the effect of different climate change narratives within the next decades, which is also a future line of research to explore.

**Author Contributions:** Conceptualization, A.R.-M. and M.P.-V.; methodology, H.M.-S. and M.P.-V.; software, A.R.-M. and H.M.-S.; validation, A.R.-M. and H.M.-S.; formal analysis, A.R.-M. and H.M.-S.; investigation, A.R.-M., H.M.-S. and M.P.-V.; resources, A.G.-P. and M.P.-V.; data curation, A.R.-M. and H.M.-S.; writing—original draft preparation, A.R.-M.; writing—review and editing, A.G.-P., H.M.-S. and M.P.-V.; visualization, A.G.-P., A.R.-M., H.M.-S. and M.P.-V.; supervision, A.G.-P. and M.P.-V.; project administration, A.G.-P. and M.P.-V.; funding acquisition, A.G.-P. and M.P.-V. All authors have read and agreed to the published version of the manuscript.

**Funding:** The data used in this study was obtained from the references included. We acknowledge the European Research Area for Climate Services consortium (ER4CS) and the Agencia Estatal de Investigación for their financial support to this research under the INNOVA project (Grant Agreement: 690462; PCIN-2017-066). This study has also been partially funded by the ADAPTAMED project (RTI2018-101483-B-I00) from the Ministerio de Ciencia, Innovación y Universidades (MICIU) of Spain.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Appendix A Appendix A**

*Water* **2020**, *12*, x FOR PEER REVIEW 16 of 20

**Figure A1.** Residual plots of the variables presented in Figure 8. **Figure A1.** Residual plots of the variables presented in Figure 8.

The residuals for the variables of Alarcon & Contreras releases, Tous releases, and water supply deficit show a lack of general pattern and are distributed pretty symmetrically around the 0 line. Total storage, however, shows a pattern that was already observed in Figure 8: the model tends to store more water at the beginning of the decade and during the drought period, and it storages less water towards the late period. Several reasons have been given to explain this pattern. As most water resource management models, stationary conditions have been assumed for water demand and reservoir operation during the whole period. However, in reality, water demand and the operation of the reservoirs was changing during the 10-year period. There is not available data to correctly represent the variation of all the water demands, but we know that the demand at the beginning of the decade was greater than during the last years, due to changes in regulation and the improvement of control. We have assumed an average water demand based on the available data. This may explain in part why the model has more water than the observed at the beginning (in reality, the water demand was greater than the introduced) and less at the end (the water demand introduced is greater in the model). The same trend can be observed in Figure A2. Regarding the impact of the operating rules in the results, the rules are based on interviews and analysis performed in collaboration with the decision-makers [42] and are, in some regard, influenced by the knowledge gained during the decade simulated in our model. In reality the logic behind the operation of the reservoirs was evolving and changing during the whole period. The residuals for the variables of Alarcon & Contreras releases, Tous releases, and water supply deficit show a lack of general pattern and are distributed pretty symmetrically around the 0 line. Total storage, however, shows a pattern that was already observed in Figure 8: the model tends to store more water at the beginning of the decade and during the drought period, and it storages less water towards the late period. Several reasons have been given to explain this pattern. As most water resource management models, stationary conditions have been assumed for water demand and reservoir operation during the whole period. However, in reality, water demand and the operation of the reservoirs was changing during the 10-year period. There is not available data to correctly represent the variation of all the water demands, but we know that the demand at the beginning of the decade was greater than during the last years, due to changes in regulation and the improvement of control. We have assumed an average water demand based on the available data. This may explain in part why the model has more water than the observed at the beginning (in reality, the water demand was greater than the introduced) and less at the end (the water demand introduced is greater in the model). The same trend can be observed in Figure A2. Regarding the impact of the operating rules in the results, the rules are based on interviews and analysis performed in collaboration with the decision-makers [42] and are, in some regard, influenced by the knowledge gained during the decade simulated in our model. In reality the logic behind the operation of the reservoirs was evolving and changing during the whole period.

*Water* **2020**, *12*, x FOR PEER REVIEW 17 of 20

**Figure A2.** Observed vs simulated storage in Alarcon and Contreras reservoirs. **Figure A2.** Observed vs simulated storage in Alarcon and Contreras reservoirs.

#### **References References**


European Commission: Brussels, Belgium, 2012.


Iacovides, I.; Ashton, V. *Gap Analysis of the Water Scarcity and Droughts Policy in the EU. SWD 380 Final,*


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

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

### **Ivan Gabriel-Martin , Alvaro Sordo-Ward \* , David Santillán and Luis Garrote**

Civil Engineering Department: Hydraulics, Energy and Environment, Universidad Politécnica de Madrid, ES 28040 Madrid, Spain; i.gabriel@alumnos.upm.es (I.G.-M.); david.santillan@upm.es (D.S.); l.garrote@upm.es (L.G.)

**\*** Correspondence: alvaro.sordo.ward@upm.es

Received: 6 February 2020; Accepted: 26 March 2020; Published: 1 April 2020

**Abstract:** The aim of this study is to contribute to solving conflicts that arise in the operation of multipurpose reservoirs when determining maximum conservation levels (MCLs). The specification of MCLs in reservoirs that are operated for water supply and flood control may imply a reduction in the volume of water supplied with a pre-defined reliability in the system. The procedure presented in this study consists of the joint optimization of the reservoir yield with a specific reliability subject to constraints imposed by hydrological dam safety and downstream river safety. We analyzed two different scenarios by considering constant or variable initial reservoir level prior to extreme flood events. In order to achieve the global optimum configuration 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 outflow (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 flood 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<sup>3</sup> ) 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 configurations depending on if their main emphasis is to maintain/increase the hydrological safety, or rather to maintain/increase the reservoir yield.

**Keywords:** hydrological dam safety; initial reservoir level; maximum conservation level; water conservation volume; flood control volume; yield reliability; regular operation; stochastic methodology

#### **1. 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

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–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–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: •


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.

#### **2. Materials and Methods**

Maximum conservation levels represent the linking variable between flood control operation and water conservation operation of the reservoir (Figure 1).

**Figure 1.** Conceptual scheme representing the role of maximum conservation levels (MCLs). 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 belowMCL 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.

**Figure 2.** Scheme of the methodology proposed.

### *2.1. 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.

### *2.2. 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.

#### 2.2.1. 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:


•

•

•

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)


**Figure 3.** The procedure to generate the ensemble of inflow hydrographs (**a**,**b**) and its seasonal **Figure 3.** 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 (MWRLTR=1000y) and the maximum outflow corresponding to a return period of 500 years (MOTR=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):


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].

#### 2.2.2. Simulation of Water Conservation Operation

We considered the yield reliability (YR) as the ratio of total volume supplied and the total volume demanded [23–26] during the period analyzed. For each configuration of MCLs (k), we obtained the reservoir yield with a reliability of 90% (Y<sup>R</sup> <sup>=</sup> 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<sup>R</sup> <sup>=</sup> 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<sup>3</sup> /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.

#### *2.3. 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: MWRLTR <sup>=</sup> 1000y, MOTR <sup>=</sup> 500y, and Y<sup>R</sup> <sup>=</sup> 90%. We assumed that the MCLs configuration would not fulfil the standards if MWRLTR=1000y was above the design flood level (DFL) and/or MOTR <sup>=</sup> 500y was greater than the emergency flow (OEMER). We compared MWRLTR=1000y, MOTR <sup>=</sup> 500y, and Y<sup>R</sup> <sup>=</sup> 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 MOTR <sup>=</sup> 500y and Y<sup>R</sup> <sup>=</sup> 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 OEMER, that is, the MOTR <sup>=</sup> 500y is the most restrictive variable. Therefore, for each scenario (Sc.1 or Sc.2), we had a set of k <sup>n</sup> MCLs configurations with k <sup>n</sup> pairs of values MOTR <sup>=</sup> 500y and Y<sup>R</sup> <sup>=</sup> 90%. The purpose was to obtain the configurations of MCLs that minimize the value of MOTR <sup>=</sup> 500y while maximizing Y<sup>R</sup> <sup>=</sup> 90%, which implies a two-objective minimization problem (Equation (1)):

$$\text{Min} \langle \mathbf{f}\_1(\mathbf{x\_i}), \mathbf{f}\_2(\mathbf{x\_i}) \rangle,\tag{1}$$

in which x<sup>i</sup> = [MCLS1, MCLS2, . . . , MCLSn], f1(x<sup>i</sup> ) = [MOTR=500yi ] and f2(x<sup>i</sup> ) = [–YR=90%i ]; being i = 1,2 . . . , k n . The solution of Equation (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 MWRLTR <sup>=</sup> 1000y and/or MOTR <sup>=</sup> 500y) and those providing a Y<sup>R</sup> <sup>=</sup> 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.

**Figure 4.** Scheme of the procedure for solutions proposal.

#### *2.4. 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<sup>3</sup> ): 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<sup>3</sup> . The capacity of the reservoir is 651 hm<sup>3</sup> at MNL (the maximum reservoir level that water might reach under normal operating conditions) [30], with a bottom dead storage of 78 hm<sup>3</sup> . The main characteristics of the Riaño Reservoir and its basin are summarized in Table 1.


**Table 1.** Main characteristics of the Riaño basin, dam, and reservoir.

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<sup>3</sup> /s (OEMER) 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:


#### *2.5. 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]).

#### **3. Results and Discussion**

#### *3.1. Determination of the Study Set of MCLs to be Analyzed*

We studied a set of MCLs that ranged from 651 hm<sup>3</sup> (volume at MNL) to 400 hm<sup>3</sup> . For the sake of simplicity, we estimated the flood hydrograph volume of Tr = 5000 years (247 hm<sup>3</sup> ) and defined a maximum flood control volume of 251 hm<sup>3</sup> . We discretized the ranges of reservoir volumes in intervals of 5 hm<sup>3</sup> . 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<sup>3</sup> ).

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<sup>3</sup> = 132,651 configurations of MCLs per scenario analyzed (Sc.1 and Sc.2.)

#### *3.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 MWRLTR <sup>=</sup> 1000y and MOTR <sup>=</sup> 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<sup>R</sup> <sup>=</sup> 90%. Figure 5 shows the values of MWRLTR <sup>=</sup> 1000y, MOTR <sup>=</sup> 500y, and Y<sup>R</sup> <sup>=</sup> 90% with respect to the MCLs of each season.

**Figure 5.** Representations of the 132,651 values of MWRLTR <sup>=</sup> 1000y, MOTR <sup>=</sup> 500y, and Y<sup>R</sup> <sup>=</sup> 90% for each of the seasons (S1, S2, and S3). MCLS1 , MCLS2 , and MCLS3 represent the maximum conservation levels in S1, S2, and S3 respectively. (**a**) Representation of the values YR=90%, with a color bar, which are the same in Sc.1 (initial reservoir level equal to MCL) and Sc.2 (variable initial reservoir level). (**b**,**c**) Representation of the values MWRLTR=1000y with a color bar in Sc.1 (**b**) and Sc.2 (**c**). (**d**,**e**) Representation of the values MOTR <sup>=</sup> 500y with a color bar in Sc.1 (**d**) and Sc.2 (**e**).

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 downstream river safety. Variations of MCLS2 are the main affection to values MWRLTR <sup>=</sup> 1000y and MOTR <sup>=</sup> 500y in Figure 5c,e, respectively.

#### *3.3. Solutions Proposal*

First, we identified the configurations which did not fulfil hydrological dam safety (MWRLTR <sup>=</sup> 1000y > DFL) and/or downstream safety (MOTR <sup>=</sup> 500y > OEMER. ) in both scenarios. A total of 826 configurations (0.6% of 132,651 configurations) had a value of MWRLTR <sup>=</sup> 1000y higher than DFL in Sc.1 (red dots in Figure 6a), while none of the configurations had MWRLTR <sup>=</sup> 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 <sup>=</sup> 500y higher than OEMER. in Sc.1, whereas 4944 (3.7%) in Sc.2 (Figure 6a,c) shows the pair of values Y<sup>R</sup> <sup>=</sup> 90% and MOTR <sup>=</sup> 500y (grey points) for each analyzed configuration in Sc.1 and Sc.2, respectively.

**Figure 6.** (**a**,**c**) Grey dots represent pair of values Y<sup>R</sup> <sup>=</sup> 90% and MOTR <sup>=</sup> 500y for Sc. 1 and Sc.2, respectively. Red dots (**a**) represent the configurations in which MWRLTR <sup>=</sup> 1000y > DFL. Magenta dots indicate the pareto front. Cyan dots represent the proposed solutions. The black dot represents the solution that maximizes the water volume supplied with a reliability of 90% and MOTR <sup>=</sup> 500y = OEMER. The blue dot represents the solution that minimizes the maximum outflow released (being Y<sup>R</sup> <sup>=</sup> 90% = Da). (**b**,**d**) MCLs for each month/season which corresponds to the proposed solutions previously selected (blue and black lines correspond to blue and black dots in (**a**) and (**c**)), for Sc.1 and Sc.2, respectively. The red dashed line represents the maximum normal level (MNL), the red dashed–dotted line represents the design flood level (DFL), and the red continuous line represents the crest of dam (COD).

≥ Da, MO ≤ O ≤ In Sc.1, the Pareto-solutions consisted of 277 configurations (Figure 6a, magenta points). One those, 98 configurations satisfied the following: Y<sup>R</sup> <sup>=</sup> 90% ≥ Da, MOTR <sup>=</sup> 500y ≤ OEMER, and MWRLTR <sup>=</sup> 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 <sup>=</sup> 500y = OEMER = 700 m<sup>3</sup> /s (Figure 6a,c, black point), Y<sup>R</sup> <sup>=</sup> 90% = 587 hm<sup>3</sup> (for Sc.1) and 600 hm<sup>3</sup>

(for Sc.2) representing an improvement (compared to Da) of 7.7% and 10.1%, respectively. On the other hand, in the case of Y<sup>R</sup> <sup>=</sup> 90% = D<sup>a</sup> = 545 hm<sup>3</sup> (Figure 6a,c, dark blue point), MOTR <sup>=</sup> 500y = 452 m<sup>3</sup> /s (for Sc.1) and 366 m<sup>3</sup> /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 <sup>=</sup> 500y = OEMER, the proposed stochastic procedure presented improvements of up to 146 hm<sup>3</sup> (26.8% of Da) compared to the worst regular operation configuration (Y<sup>R</sup> <sup>=</sup> 90% = 454 hm<sup>3</sup> 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,c, respectively. In both scenarios, the highest MCLs were associated with S3, while the lowest was associated with S1.

#### *3.4. 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 <sup>=</sup> 500y = OEMER = 700 m<sup>3</sup> /s, MWRLTR <sup>=</sup> 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<sup>R</sup> <sup>=</sup> 90% of the system. In the case of Y<sup>R</sup> <sup>=</sup> 90% = Da = 545 hm<sup>3</sup> , the differences between Sc.2 and Sc.1 for MWRLTR <sup>=</sup> 1000y and MOTR <sup>=</sup> 500y were 0.2 m.a.s.l and 86 m<sup>3</sup> /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 <sup>=</sup> 500y if variable initial level was considered.

**Figure 7.** Spider plot comparing the value of the maximum conservation levels for different seasons (MCLsx), MWRLTR <sup>=</sup> 1000y, MOTR <sup>=</sup> 500y, and Y<sup>R</sup> <sup>=</sup> 90%. In black, solutions that maximize the volume supplied with MOTR <sup>=</sup> 500y = OEMER. In blue, the solutions that supply the current annual demand minimize the maximum outflow released. The continuous line corresponds to Sc.1 and the dashed line to Sc.2.

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.

### **4. Conclusions**

The main conclusions extracted from this study are as follow:


**Author Contributions:** Conceptualization, I.G.-M., A.S.-W., D.S. and L.G.; Data curation, I.G.-M. and A.S.-W.; Formal analysis, I.G.-M., A.S.-W., D.S. and L.G.; Funding acquisition, A.S.-W.; Investigation, I.G.-M., A.S.-W., D.S. and L.G.; Methodology, I.G.-M., A.S.-W., D.S. and L.G.; Resources, A.S.-W.; Software, I.G.-M.; Visualization, I.G.-M.; Writing – original draft, I.G.-M.; Writing – review & editing, A.S.-W., D.S. and L.G. All authors have read and agreed to the published version of the manuscript.

#### **Funding:** Universidad Politécnica de Madrid: VJIDOCUPM19AFSW.

**Acknowledgments:** The authors acknowledge the computer resources and technical assistance provided by the Centro de Supercomputación y Visualización de Madrid (CeSViMa) and the funds from Universidad Politécnica de Madrid in the framework of their Program "Ayudas para contratos predoctorales para la realización del doctorado en sus escuelas, facultad, centro e institutos de I+D+i" and "ayudas a proyectos de I+D de investigadores posdoctorales".

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
