Implementation and sensitivity analysis of the Dam-Reservoir OPeration model (DROP v1.0) over Spain

The prediction of water resource evolution is considered to be a major challenge for the coming century, particularly in the context of climate change and increasing demographic pressure. Water resources are directly linked to the continental water cycle, and the main processes modulating changes can be represented by global hydrological models. However, anthropogenic impacts on water resources, and in particular the effects of dams-reservoirs on river flows, are still poorly known and generally neglected in coupled land surface–river routing models. This paper presents a parameterized reservoir model, DROP (Dam-Reservoir OPeration), based on Hanasaki's scheme to compute monthly releases given inflows, water demands and the management purpose. With its significantly anthropized river basins, Spain has been chosen as a study case for which simulated outflows and water storage variations are evaluated against in situ observations over the period 1979–2014. Using a default configuration of the reservoir model, results reveal its positive contribution in representing the seasonal cycle of discharge and storage variation, specifically for large-storage capacity irrigation reservoirs. Based on a bounded version of the Nash–Sutcliffe efficiency (NSE) index, called C_2M, the overall outflow representation is improved by 43 % in the median. For irrigation reservoirs, the improvement rate reaches 80 %. A comprehensive sensitivity analysis of DROP model parameters was conducted based on the performance of C_2M on outflows and volumes using the Sobol method. The results show that the most influential parameter is the threshold coefficient describing the demand-controlled release level. The analysis also reveals the parameters that need to be focused on in order to improve river flow or reservoir water storage modeling by highlighting the difference in the individual effects of the parameters and their interactions depending on whether one focuses on outflows or volume mean seasonal patterns. The results of this generic reservoir scheme show promise for modeling present and future reservoir impacts on the continental hydrology within global land surface–river routing models.

1 Introduction

Dams are used to provide essential services to mankind in terms of economic, environmental and social impacts. They provide water supply for domestic, industrial and irrigation needs, enable hydroelectric power generation and river navigation and prevent extreme hydrological events. There are currently more than 58 700 large dams (heights >15 m) worldwide, with an estimated cumulative storage capacity between 7700 and 8300 km³ (Vörösmarty et al., 2003; Downing et al., 2006; Lehner et al., 2011; ICOLD, 2020). When including millions of smaller dams (∼16.7 M larger than 0.01 ha; Lehner et al., 2011), the total global impounded water may even exceed 10 000 km³ (Chao et al., 2008). This means that reservoirs hold more than 4 times the amount of water stored in rivers (the average annual river water storage ranges from approximately 1200 to 2120 km³), and they account for approximately 20 % of average annual river flow to the oceans (40 000–45 500 km³ yr⁻¹) (Baumgartner and Reichel, 1975; Oki and Kanae, 2006; Syed et al., 2010; Lehner et al., 2011). More than 60 % of the world's largest rivers are fragmented by the construction of dams, which account for 90 % of the flow from these river basins (Dynesius and Nilsson, 1994; Revenga et al., 2005; Grill et al., 2015, 2019).

Several studies have demonstrated the significant impact of reservoirs on river flow regimes at not only local scales, but also at larger regional and global scales: reservoirs impact the magnitude of downstream river flows and alter the temporal pattern of river discharge over the continental surface (Haddeland et al., 2006; Hanasaki et al., 2006; Döll et al., 2009; Biemans et al., 2011; Shin et al., 2019; Gutenson et al., 2020). Through surface evaporation and water exchanges with groundwater, lakes and floodplains, reservoirs not only affect the water budget over land, but also throughout the Earth's horological cycle by interacting with the atmosphere and oceans (Pokhrel et al., 2012; Zhao et al., 2012; Wada et al., 2016; Frederikse et al., 2020). There is therefore an increasingly pressing need to integrate reservoir operations in large-scale land surface and global hydrological models (LSMs-GHMs) to overcome the existing biases in continental water cycle and river flow modeling given the number of highly regulated basins.

Models developed to date which represent reservoir releases at a large scale can be categorized as data-driven and process-based approaches. The first category of models is built on the basis of observed release data, water levels and volumes. These methods range from simplified representations of reservoir operation using linear or multilinear regression (e.g., Young, 1967; Raman and Chandramouli, 1996) to very sophisticated models based on machine learning and artificial intelligence techniques, such as neural-network-based methods (e.g., Maier and Dandy, 2000; Razavi and Karamouz, 2007; Ehsani et al., 2016; Coerver et al., 2018). However, this approach requires specific knowledge of the studied reservoirs and requires access to a large amount of observed data. This approach also remains limited in its applicability for future predictions since the generated operational rules are based on historical data and thus do not take into account future potential socioeconomic and predicted climatic changes in the operation of these reservoirs. Process-based approaches, on the other hand, are based on conceptualizing reservoir responses according to its operational purpose by linking release control to physical processes, such as crop growth and the associated water requirements, or to water and energy demands that can be estimated at the global scale. The representation of dam operations is thus achieved without having to explicitly observe the actual release operations performed on each reservoir (Gutenson et al., 2020). The best-known schemes in this category are those developed by Hanasaki et al. (2006) and Haddeland et al. (2006), which are inflow- and demand-based models as presented by Yassin et al. (2019). These two generic models have been implemented in several global hydrological and water management models (e.g., WaterGAP, Döll et al., 2009; VIC, Haddeland et al., 2006; H08, Hanasaki et al., 2008; PCR-GLOBWB, Van Beek et al., 2011).

Of all the studies which have been carried out with these models, very few have been focused on Iberian Peninsula basins, where the prevailing semi-arid climate leads to a greater necessity to store water in large-capacity reservoirs, which leads to larger reservoir effects on rivers (Batalla et al., 2004; López-Moreno et al., 2009; Lorenzo-Lacruz et al., 2010). Spain is, in fact, among the top 10 dam-building countries, with more than 1064 dams (ICOLD, 2020). The study of Grill et al. (2019), which assessed global river connectivity, revealed that the regulation effect of dams is the dominant pressure source in Spain's rivers, where very high degrees of regulation cause an alteration of the natural river flow regime for its five main rivers.

A large number of studies can be found in the literature focusing on the sensitivity analysis of various models, in various science fields, such as machine learning (e.g., Zouhri et al., 2022) or civil engineering (e.g., Zamanian et al., 2021). Pianosi et al. (2016) provided a classification of the sensitivity analysis methods used in environmental sciences and their benefits. A comprehensive sensitivity analysis, as provided by global methods such as Sobol indices, is essential for a precise understanding of a model (Saltelli et al., 2008; Saltelli, 2013). It can help to improve parameter calibration efficiency and avoid overparameterization (e.g., Shin and Jung, 2022; Tang et al., 2007a, b). It is also an efficient tool to better understand the model structure (Saltelli et al., 2008), its uncertainties (e.g., Pheulpin et al., 2022) and the dominant processes under various conditions (e.g., Huang et al., 2021; Zhang et al., 2013). If understanding and quantifying uncertainties is necessary for hydrological modeling, it is also particularly critical to efficiently weight observations and model states in data assimilation techniques (Liu and Gupta, 2007; Abdolghafoorian and Farhadi, 2016). Finally, for a potential extension to the global scale, in which the model parameters cannot be calibrated and validated using observations (that do not exist or are not accessible in most reservoirs of the world), a full understanding of the model sensitivity to the parameters is crucial, especially when models are used as support tools for decision-making (Herrera et al., 2022). A few sensitivity tests on the Hanasaki scheme can be found in the literature (e.g., Hanasaki et al., 2006; Shin et al., 2019), but all of them only focused on one or two parameters supposed to be the most sensitive, which may not be sufficient for a thorough understanding of the impact of parameter uncertainties (Saltelli and Annoni, 2010). To our knowledge, no comprehensive sensitivity analysis has been conducted on this scheme yet.

This study proposes a global and parameterized reservoir model, DROP (Dam-Reservoir OPeration), built on the basis of the generic scheme by Hanasaki et al. (2006), and the aim of the study is to implement it in Spain as a case study to represent reservoir releases and to provide a comprehensive understanding of how uncertainties in each of the model parameters are affecting the overall accuracy of its predictions. In this regard, the study aims to set up a comprehensive sensitivity analysis based on Sobol indices in order to quantify the influence of each parameter and to reveal the types of influence that each of them holds by dissociating individual effects from possible interactions. Further work will focus on implementing this scheme in global hydrological models in order to provide a physical representation of reservoirs on large scales.

This paper is organized as follows: Sect. 2 provides a description of the DROP model and a theoretical outline of the sensitivity analysis method. The study area in Spain, the available observational data, the model setup and the sensitivity analysis implementation are described in Sect. 3. Sections 4–5 illustrate and discuss the model's overall performance and the parameter sensitivity results. Conclusions and perspectives of the study are presented in the last section.

2 Methodology

2.1 DROP

The parameterized DROP model has been developed based on the Hanasaki et al. (2006) scheme. The model works at the level of each reservoir individually and is based on the mass balance (as shown below in Eq. 1) of each reservoir to calculate the release at its outlet.

$$\begin{array}{}\text{(1)}& {\displaystyle \frac{\mathrm{d}V}{\mathrm{d}t}}=Q_{\mathrm{in}}-Q_{\mathrm{out}}\end{array}$$

Q_in and Q_out stand, respectively, for the net inflow to the reservoir and its outflow at the outlet. In fact, the net inflow, Q_in, combines different physical processes, as shown in Eq. (2): it includes water inputs from precipitation P, direct runoff R_d representing water flows running off the surrounding ground surface that also feeds the reservoir, and tributary streamflows Q_trib flowing into the reservoir. Evaporation E and groundwater exchange Q_gw are also included in Q_in.

$$\begin{array}{}\text{(2)}& Q_{\mathrm{in}}=(P-E)\times A_{\mathrm{reservoir}}+Q_{\mathrm{trib}}+R_{\mathrm{d}}\pm Q_{\mathrm{gw}},\end{array}$$

where A_reservoir is the reservoir surface area.

In order to simulate dam releases, Q_out, the reservoir model categorizes reservoirs as irrigation and non-irrigation reservoirs, computes the mass balance for each reservoir individually, and calculates releases based on inflow and water demands. A schematic of the DROP model is shown in Fig. 1 with its six parameters (in blue). To give an overview of its overall functioning, the model takes as input inflow and water demands and calculates at a monthly time step dam releases according to the reservoir management purpose and its relative capacity compared to inflow, denoted as c. Operating rules are set following an “operational year”. This type of year has been introduced by Hanasaki et al. (2006), and it differs from hydrological and calendar years. It starts the first month of the water release period and is therefore specific to each reservoir. This specific month, representing the first parameter of the model, will be denoted as m_start in the remaining sections of this article. At the beginning of each operational year, defined by m_start, reservoir storage volume S_init is compared to the ideal filling value S_ideal, defined as a ratio, α, of storage capacity C. Reservoir-simulated releases are impacted by this step as they are retrospectively revised upwards or downwards depending on whether the reservoir has more or less water storage than the ideal rate. Dam monthly releases are then computed following two steps: first, a provisional release is calculated based on water demands and the annual mean inflow. In irrigation reservoirs, two parameters are involved: d_max, setting the control area of the reservoir and thus the water needs to be supplied, and M, defining the minimum release to be provided for environmental requirements. The calculation scheme remains simplistic for other management purposes where releases are set to mean annual inflow. The provisional release is then corrected by incorporating a “demand-controlled release” ratio R, controlled by two parameters c_threshold and b, which accounts for inflow pattern influence in reservoirs with low storage capacities compared to inflow. The DROP model therefore has six parameters, as shown in Fig. 1.

Figure 1 Schematic representation of the DROP model showing its six parameters (in blue): m_start defines the start of the reservoir operational year. α sets the ideal filling rate of the reservoir to be reached at each starting month of the year. The remaining four parameters are used in dam release computation: d_max and M, operating only in irrigation reservoirs, set the control area of the reservoir (and therefore irrigation water demands) and the minimum release to be provided, respectively. c_threshold and b, on the other hand, account for inflow pattern influence in reservoirs with low relative storage capacity to inflow (i.e., run-of-river reservoirs).

The operating rules are detailed below. First, at the beginning of each operational year, an annual release coefficient K_y is computed to determine the filling rate of the reservoir. K_y is a ratio between the reservoir storage at the beginning of the year S_first,y and the long-term target storage S_ideal:

$$\begin{array}{}\text{(3)}& K_{\mathrm{y}}=S_{\mathrm{first},\mathrm{y}}/S_{\mathrm{ideal}},\end{array}$$

where $S_{\mathrm{ideal}}=\mathit{\alpha }\times C$, α being a nondimensional constant, and C is the total storage capacity of the reservoir. S_ideal represents the ideal filling level at the beginning of each year; α was set semi-empirically to 0.85 of the maximum capacity for all reservoirs based on sensitivity tests conducted in Hanasaki et al. (2006). This annual release coefficient is the one used to weight the provisional releases calculated afterwards depending on whether the reservoir has more (K_y>1) or less (K_y<1) water storage than the ideal rate. The underlying aim is to reduce interannual variation in streamflow.

The following steps describe how reservoir releases are computed. The provisional monthly releases, $r_{\mathrm{m}}^{\prime }$, are set depending on the reservoir's main purpose. The scheme is simplified for a non-irrigation reservoir since it constantly releases the mean annual inflow i_mean calculated over the whole simulation period (noted “long-term” in the remaining sections of this article):

$$\begin{array}{}\text{(4)}& r_{\mathrm{m}}^{\prime }=i_{\mathrm{mean}}.\end{array}$$

For an irrigation reservoir, monthly releases are proportional to water demands. They are parameterized here following Shin et al. (2019) as

$$\begin{array}{}\text{(5)}& r_{\mathrm{m}}^{\prime }=\left\{\begin{array}{ll}{i}_{\mathrm{mean}}\left[(M+(\mathrm{1}-M)\frac{d_{\mathrm{m}}}{d_{\mathrm{mean},\mathrm{y}}}\right]& \text{ if }\text{DPI}>\mathrm{1}-M,\\ i_{\mathrm{mean}}+d_{\mathrm{m}}-d_{\mathrm{mean},\mathrm{y}}& \text{ if }\text{DPI}\le \mathrm{1}-M,\end{array}\right.\end{array}$$

where d_m and d_mean,y are, respectively, monthly and annual mean water demands. DPI (demand per inflow) is, as introduced by Shin et al. (2019), the ratio between d_mean,y and i_mean, and M represents the minimum monthly release as a percentage of i_mean. It is set to 0.5 by Hanasaki et al. (2006) and Döll et al. (2009) and to 0.1 by Biemans et al. (2011) and Shin et al. (2019).

When the DPI is above the set threshold (1−M), water needs are considered to be very high and can only be partially fulfilled by the water stored throughout the coming year. The priority is to first ensure a minimum release, M×(i_mean), in order to meet the environmental flow requirements. The remaining part of the annual inflow is released throughout the year on a monthly basis following the sub-annual water demand fluctuation curve. Otherwise, when the DPI is below (1−M), the reservoir releases all the monthly water demand that is needed.

In this scheme, only irrigation demand is considered. Since dams provide water for the downstream demand within a certain distance, a maximum distance, d_max, is a parameter to define, for each reservoir, the irrigated grid cells within the river basin to be supplied and thus delimits a “command area” for each reservoir. In contrast to Hanasaki et al. (2006) (where a crop growth model is considered to calculate the irrigation demand), here the distribution of irrigated areas is based on ECOCLIMAP SG (Calvet and Champeaux, 2020; Druel et al., 2022), which is used by the irrigation module in the ISBA LSM (Druel et al., 2022) to compute irrigation demands for 5×5 km-resolution grid cells (see Sect. 3). The irrigation water demands are aggregated within the command area for each reservoir. Unlike the Hanasaki et al. (2006) scheme, the DROP model entirely separates irrigation from non-irrigation reservoirs and only accounts for irrigation water demands for this category of reservoirs. Moreover, industrial and domestic demands are highly uncertain and hardly accessible. In fact, these consumptions are strongly linked to each country's specific political and economical policies, and extensive records are generally not made public. Some global databases, such as the Food and Agriculture Organization of the United Nations (FAO) global information system AQUASTAT (AQUASTAT, 1994), provide estimates of average annual water withdrawals by activity sector and by water resource publicly accessed at the country level. However, such estimates are still limited in terms of temporal and spatial resolution and would consequently add a non-negligible source of error to the model. They are therefore neglected in this study.

In all possible cases, regardless of the reservoir purpose, the water released over the operational year is equivalent to the long-term mean annual inflow.

The release computed so far is provisional. The real monthly release is calculated as follows:

$$\begin{array}{}\text{(6)}& r_{\mathrm{m}}=\left\{\begin{array}{ll}{K}_{\mathrm{y}}\times r_{\mathrm{m}}^{\prime }& \text{ if }c\ge c_{\mathrm{threshold}},\\ (\mathrm{1}-R)\times i_{\mathrm{m}}+R\times K_{\mathrm{y}}\times r_{\mathrm{m}}^{\prime }& \text{ if }\mathrm{0}\le c<c_{\mathrm{threshold}},\end{array}\right.\end{array}$$

where c is the relative capacity of a reservoir and is defined as the ratio between storage capacity C and the long-term mean annual inflow water volume ($c=C/I_{\mathrm{mean}}$). The parameter R, as introduced by Shin et al. (2019) (and also parameterized as β by Horan et al., 2021), is a so-called “demand-controlled release ratio”. It is parameterized in the current study as

$$\begin{array}{}\text{(7)}& R={\left({\displaystyle \frac{c}{c_{\mathrm{threshold}}}}\right)}^b,\end{array}$$

where R describes the influence of the inflow regime on release for small storage capacity reservoirs. It varies between 0 and 1 and includes two parameters: c_threshold and the b coefficient. In fact, the smaller the reservoir capacity is compared to inflow, the closer it gets to run-of-river dams where release is close to the natural river flow and thus the influence of the inflow annual pattern. Otherwise, when c is above the threshold (large-capacity dams), then R=1 and the release is fully controlled by water demand. c_threshold and b are set, respectively, as 0.5 and 2 in both Hanasaki et al. (2006) and Biemans et al. (2011). Shin et al. (2019), on the other hand, proposed an analytical formula to compute R, which is reproduced in this study, and set c_threshold and b, respectively, to $\mathrm{1}/\mathit{\alpha }$ and 1. The monthly releases are also weighted by the annual coefficient, K_y from Eq. (3), which is calculated at the beginning of the operational year and describes how full the reservoir is relative to the ideal rate. Reservoir releases can therefore be revised upwards or downwards depending on whether the reservoir had more or less water storage than the ideal rate.

The modifications brought to the model from the previous version of Hanasaki et al. (2006) and the Shin et al. (2019) parameterization can be summarized in the following list.

• The starting month of the operational year, which was calculated in previous versions at the level of each reservoir based on observed inflows and dam releases, is considered in this model version as a parameter, denoted as m_start. This overcomes the challenge faced by the old versions when applying the model in ungauged basins. Within our version, m_start is to be set separately for irrigation and non-irrigation reservoirs.

• The same applies to d_max, initially set to a constant value depending on the spatial resolution of the river routing model in which the scheme is being implemented in each study. Here, d_max is defined as a parameter to be set for reservoir command area delimitation.

• d_m and d_mean,y in Eq. (5) are monthly and annual mean irrigation water demands, respectively. Industrial and domestic water demands are not taken into account, unlike in Hanasaki et al. (2006).

• M, set to 0.1 in the Biemans et al. (2011) and Shin et al. (2019) schemes and 0.5 in Hanasaki et al. (2006), is considered as a parameter in this version, keeping the same notation used in Shin et al. (2019).

• A more explicit parameterization of the demand-controlled release ratio R is provided here compared to what is proposed in Shin et al. (2019). The generalized formulation of R in Eq. (7) highlights the last two parameters of the present model version. This enables one to distinguish the dual role of c_threshold in Eq. (6) from that of the coefficient b, which only describes the transition made between a demand-controlled reservoir and an input-controlled reservoir. Hanasaki et al. (2006) and Shin et al. (2019) have set (c_threshold, b) to (0.5,2) and ($\mathrm{1}/\mathit{\alpha }$,1), respectively.

A description of the reservoir model parameters is given in Table 1.

Table 1 Description of the DROP model parameters.

The reservoir volume is derived at each time step from the water balance. Boundary conditions are defined considering two possible scenarios: (i) if the reservoir is full, the excess water is spilled. (ii) When reservoir storage falls below 10 % of the capacity, the reservoir reaches the dead storage zone and water release is prevented.

In order to run the model, reservoir characteristics, such as the storage capacity and main purpose, are needed. The model also requires continuous time series of inflow and water demands to compute releases. In the current study, all of the modeled reservoirs are located in Spain, where the physical characteristics, in situ observations of natural and anthropized flows, and storage volumes are publicly available. Section 3 describes the global and local datasets used herein.

2.2 Sensitivity analysis: Sobol's method

The Sobol sensitivity analysis method (Sobol, 1993) has been widely used in hydrological models in order to identify the parameters that contribute the most to the model output uncertainty. It is a global variance-based approach where, for a chosen variable of interest Y, the total variance is decomposed into fractions attributed to each individual input as well as the interactions between them. If Y=f(X) is a goodness-of-fit metric of the model with X representing the set of parameters $X_{\mathrm{1}},X_{\mathrm{2}},\mathrm{\dots }{X}_p$ of size p, the total variance of Y, D(Y), is decomposed as follows (Zhang et al., 2013):

$$\begin{array}{}\text{(8)}& D\left(Y\right)=\sum _i{D}_i+\sum _{i<j}{D}_{ij}+\sum _{i<j<k}{D}_{ijk}+\mathrm{\dots }+D_{\mathrm{12}..p},\end{array}$$

where D_i is the amount of variance due to a parameter X_i alone, D_ij is the amount of variance arising from the interaction between the parameters X_i and X_j, and so on. Equation (8) is known as Hoeffding decomposition (Hoeffding, 1992).

The sensitivity indices, called Sobol indices, are computed as ratios between the component variances and total variance in order to measure the contribution of each single parameter and each parameter interaction. The first-order Sobol index, S1, captures the sensitivity of Y to each input parameter X_i taken alone. The second-, S2, and higher-order indices describe the contribution of the multiple interactions between parameters. A model with p parameters requires 2^p−1 indices to be evaluated, which rapidly becomes computationally challenging for high values of p. The total-order index, ST, measures the full influence of a parameter by including all of the variance caused by its interactions with the rest of the parameters, all orders included. The first-, second- and total-order Sobol sensitivity measure formulas are

$$\begin{array}{}\text{(9)}& {\displaystyle }S{\mathrm{1}}_i=& {\displaystyle \frac{D_i}{D}},\text{(10)}& {\displaystyle }S{\mathrm{2}}_{ij}=& {\displaystyle \frac{D_{ij}}{D}},\text{(11)}& {\displaystyle }{\text{ST}}_i=& {\displaystyle }\mathrm{1}-{\displaystyle \frac{D_{\sim i}}{D}},\end{array}$$

where D_∼i is the amount of variance due to all of the parameters except for X_i.

Model parameter sampling and Sobol index estimation are performed here using the open-source Python library SALib (Herman, 2017). The parameter samples are generated following quasi-random sequences (Saltelli, 2002) in order to scatter the sample points as uniformly as possible over the parameter space. Following the theorem, the different orders of Sobol indices are estimated from a total number of model runs of $N\times (\mathrm{2}\times d+\mathrm{2})$, where N and d are, respectively, the sample size and the number of parameters. The package also provides confidence intervals of the index estimates at the 95 % confidence level.

3 Study area and data

3.1 Spain river basins description

Spain has an estimated area of 505 983 km² (INE, 2022), representing more than 85 % of the Iberian Peninsula (estimated area 583 254 km², Lorenzo-Lacruz et al., 2012). It has river basins with sizes ranging from a few square kilometers to more than 80 000 km² (i.e., the Ebro basin). The five largest river basins are the Ebro, flowing into the Mediterranean Sea, and the Duero, Tajo, Guardiana and Guadalquivir, flowing into the Atlantic Ocean. They also represent the river basins where the natural river flow pattern is the most significantly altered in the country, with high degrees of regulation resulting from extensive reservoir construction. These different rivers have a Mediterranean hydrological pattern, which is characterized by high flows during the wet season (i.e., fall and winter) and low flows during its dry season (i.e., late spring, summer). This seasonality explains the strong anthropization of the Spanish hydrographic basins and, in particular, the construction of more than 1200 dams mainly in the second half of the 20th century. They are essential for retaining enough water to meet irrigation and domestic water demands as well as for hydropower production, which represents ∼13 % of Spanish electricity generation (IEA, 2022, in 2020,).

3.2 Dataset

The data series required as inputs to the model and the ones used to validate the outputs are taken from the Spain database. Reservoir characteristics are taken from the global GRanD database. Irrigation demands are simulated by the ISBA irrigation module. The needed input data are detailed below.

• Local (Spain) database. In situ observations of natural and anthropized flow and volume data are made publicly available by the Center for Hydrographic Studies of Spain (CEDEX, Ministry of Public Works and Ministry for Ecological Transition, Spain). The national database includes the location and daily time series of discharge for 1119 gauge stations and outflows from 347 reservoirs over the period 1900–2014.

• The Global Reservoir and Dam (GRanD) database, from which the general characteristics of dams are taken (Lehner et al., 2011). Version 1.3 published in 2019 includes 7320 reservoirs and provides the geographical locations of dams as well as attribute information such as the construction year, maximum storage capacity, surface area and set of purposes of reservoirs; 263 of the reservoirs listed in GRanD are located in Spain. Most of them were built from 1955 to 2000, reaching a total storage capacity of 56 480 hm³ in 2016, which exceeds the mean annual streamflow of the eight major rivers of the Iberian Peninsula (55 850 hm³ yr⁻¹) (Lorenzo-Lacruz et al., 2012). Irrigation and hydroelectricity are the most identified purposes of more than half the reservoirs, and water supply comes third. In this paper, all management purposes different from irrigation are grouped in the “other purposes” category.

• Simulated irrigation demands. The irrigation water demands are simulated by the new irrigation scheme implemented in the ISBA LSM (Druel et al., 2022). It uses the ECOCLIMAP-SG land cover classification to identify the areas within the grid cell that can be irrigated. The three main parameters are set in the model to control, respectively, the irrigation triggering (when the plant reaches the wilting point), the period of crop growth where irrigation is possible (between emergence and harvest) and the amount of water used for irrigation. In the default configuration, a predefined amount of water of 30 mm is set for each irrigation and a 7 d minimum return period is fixed between two irrigation operations. These parameters can be user-defined for each vegetation type and in each grid cell. To generate irrigation demands over Spain, we used the 5 km resolution SAFRAN-based meteorological datasets for that country (Quintana-Seguí et al., 2017) that were available over the period 1979–2014 to force the ISBA land surface model. Set to the default configuration, the irrigation module within the LSM computed daily irrigation demands for each grid cell over that time period. At the input of the reservoir model, d_max is set beforehand and delimits for each reservoir a command area made up of a selection of grid cells within the same basin at a lower altitude. The equivalent amount of irrigation water requested from a reservoir in a given month corresponds to the aggregated daily irrigation demands of all crop types within the retained grid cells.

3.3 Model setup: pre-processing steps

3.3.1 Cross-referencing global and local datasets

Out of the 263 reservoirs listed in GRanD, only 216 were kept after cross-referencing the two databases and for which both the characteristics and time series of release and volume could be identified. In fact, two reservoirs were doubly identified in GRanD v1.3 (the IDs were 2882 and 2844) because they were rebuilt and/or renamed; their most recent characteristics are those retained. The 45 remaining reservoirs, located mainly in the northwest and south of Spain, were not identified in the Spanish database since they were built after 2014.

The maximum storage capacity of the chosen reservoirs goes from 9.5 (the San Lorenzo Mongay dam, located on the Segre River in the Ebro basin) to 3200 hm³ (the La Serena dam on the Zujar River in the Gardiana basin), with a mean of 236 hm³ and a standard deviation of 441 hm³. Using the latitude and longitude information, the chosen reservoirs were located on the new $\mathrm{1}/\mathrm{12}$^∘-resolution river network derived from the CTRIP river routing model (Munier and Decharme, 2022). Dam locations were adjusted so as to have comparable drainage surfaces between those given by the GRanD database and those estimated by CTRIP.

An initial analysis conducted on observed river flows upstream and downstream of the reservoirs has identified a common seasonal behavior among those with irrigation, which is that the peak dam release is shifted in time from the natural inflow. This is due to the typical operating mode of these reservoirs, which are designed to retain water arriving upstream during winter (wet season) and release it during summer (dry season) to meet irrigation needs.

3.3.2 Reconstructing inflow

The reservoir scheme requires reservoir net inflows and water demands as inputs to estimate volume variations and outflows. In existing studies, water demands are estimated and net inflows are either modeled by a river routing scheme or estimated from gauge observations in rivers and tributaries upstream the reservoir. Reservoir abstraction is also sometimes accounted for in the reservoir water balance. However, some processes are usually neglected, such as precipitation interception, direct runoff, evaporation or groundwater exchanges. This introduces a bias into the water budget and consequently increases the model uncertainties, especially when inflows are derived from land surface and river routing models (see, e.g., Vanderkelen et al., 2022). In this study, the water balance (Eq. 2) is used in a first step to derive the result of all these components (net inflow) from observations of reservoir volume and outflow, as described below. The advantage compared to previous studies is that it removes the uncertainties related to each of these components, enabling the analysis of the model uncertainties themselves and the capacity of the model to reproduce the reservoir behavior alone, without additional uncertainties or potential compensations between the components of the water budget and the model parameters.

Note that, in a future work, the DROP model aims to be coupled to a series of models that can represent these different processes. In the ISBA-CTRIP land surface–river routing model, for example (Decharme et al., 2019), evaporation is computed by FLake, a module representing energy balance in lakes (Le Moigne et al., 2016), direct runoff is derived from the ISBA land surface model, and inflows from tributaries as well as groundwater exchanges are computed by the CTRIP river routing model (Munier and Decharme, 2022).

For each of the selected reservoirs, the longest continuous common periods of daily observed outflows and volumes were first determined. At this stage, only reservoirs with more than 3-year time series were retained, which leaves 215 reservoirs to be simulated. The net inflows were then derived from outflows and volume variations at the daily scale (Eq. 1). The computed net inflows were then corrected by removing outliers in two steps: first, all peak flows were selected when their maximum value exceeded 5 times the long-term mean. Among the peak values, outliers are then identified when the relative difference of slopes before and after peak flow is less than 10 %. Both thresholds were set empirically. This step differentiates the outliers from the hydraulic behavior of a river in flood recession periods. The outlier is replaced by a linearly interpolated value. The length of the corrected time series goes from 3.5 to 34 years with a median of 23 years. The main purposes, simulation period lengths and relative capacities of the 215 reservoirs simulated are shown in Fig. 2.

Figure 2 Main characteristics of the chosen reservoirs. (a) Classification of their main purpose, (b) distribution of simulation period lengths, and (c) histogram showing decimal log values of their relative capacity.

We note a good distribution of management purposes and relative inflow capacities in the final selected reservoirs. Overall, half the reservoirs are primarily used for irrigation, which is mainly due to the semi-arid climate of the Iberian Peninsula and the high needs of irrigation in the country. The rest of the reservoirs are allocated to hydropower generation, water supply and different other purposes with respective percentages of 29 %, 16 % and 5 % and are grouped in the non-irrigation reservoir category.

3.4 Sensitivity analysis implementation

A sensitivity analysis with respect to the six parameters was conducted on the performance of a Nash–Sutcliffe efficiency (NSE, Nash and Sutcliffe, 1970) bounded version, called C_2M (Mathevet et al., 2006), on outflows using the Sobol method. In fact, the NSE values in some reservoirs were highly negative for some simulations, and thus this metric was not suitable for a variance-based sensitivity analysis method like Sobol's. C_2M is used instead as it is a normalized version of NSE that varies between −1 and 1 and where all negative values are bounded between 0 and −1. Parameter default values, bounds and distributions are listed in Table 2. The parameter distributions were all considered uniform except for relative capacity, for which the distribution is logarithmic to align with the observed pattern on the modeled reservoirs (Fig. 2).

Table 2 Summary of the DROP model parameter default values and feasible ranges for the sensitivity analysis.

The default values for α, M, c_threshold and b are those considered in Hanasaki et al. (2006). α-selected bounds cover a realistic range of ideal reservoir filling rates. The distribution and bound values of the c_threshold parameter are drawn from the relative capacity distribution of the 215 modeled reservoirs (Fig. 2). The b lower limit is 0.5, below which the shift to the demand-controlled state of a reservoir becomes too abrupt once c exceeds c_threshold. The upper limit considered for b is 5, beyond which the transition curve between the behavior of a low relative capacity reservoir and a high relative capacity one becomes unchanged, following the sensitivity test run conducted on this parameter. For d_max, both the default value and the lower and upper limits were set to be consistent with the sizes of Spanish river basins. The operational year starting month, m_start, is set to April for irrigation reservoirs in order to match the beginning of the crop irrigation season considered in the irrigation model. For other reservoirs, m_start is chosen empirically to be May as the default value based on the observed filling curves of the reservoirs, which tend to be at the maximum filling level near May.

The sensitivity analysis was performed on each of the 215 reservoirs separately, distinguishing between irrigation and non-irrigation reservoirs since the number of parameters involved depends on the main purpose of the reservoir (six and four, respectively, as d_max and M are only considered in irrigation reservoirs).

Using Saltelli's quasi-random sampling method (Saltelli et al., 2010), a sample size of 4096 was used for this analysis in each category, resulting in $\mathrm{4096}\times (\mathrm{2}\times \mathrm{6}+\mathrm{2})=\mathrm{57}\mathrm{344}$ and $\mathrm{4096}\times (\mathrm{2}\times \mathrm{4}+\mathrm{2})=\mathrm{40}\mathrm{960}$ model runs for each of the 107 irrigation and 108 non-irrigation reservoirs, respectively. By comparing with smaller sample sizes, the confidence intervals estimated by SALib show that the Sobol indices converged and that the chosen sample size is sufficient to reliably represent the results.

4 Results

This section presents the main results of this study. First, simulation results of the reservoir releases using the default configuration of the model are displayed. Then, a sensitivity analysis of the model parameters is presented.

4.1 Reproducing the flow seasonal shift in irrigation reservoirs

Using the default parameterization with the parameters listed in Table 2, an operating rule is determined for each of the 215 reservoirs, and both outflows and water storage variations are simulated through complete operational years within the overall period 1979–2014. Figure 3 shows the C_2M values at reservoir outlets by evaluating the simulated monthly outflow time series against in situ observations. Panel a shows the results obtained using a reference simulation where rivers are considered in their natural state (Q_out=Q_in). Panel b shows the C_2M improvement rate when considering the DROP model.

Figure 3 Results of the DROP model contribution (default configuration) in river flow modeling at reservoir outlets: (a) C_2M performance metrics of non-reservoir simulation (reference simulation where Q_out=Q_in) and (b) C_2M improvement rates by integrating the reservoir model. Symbols represent the different reservoir management purposes. They come in two different sizes: larger if c≥0.5 and smaller if c<0.5.

Overall, with DROP, the river flow representation is clearly improved at nearly all the reservoirs' locations. The median C_2M index for flows is 0.52, which corresponds to a 43 % improved flow representation when compared to natural river representation. For storage volume representation, the mean correlation is 0.53 with a standard deviation of 0.3. Half of the reservoirs have a correlation greater than 0.63 between observed and simulated storage volumes.

The results reveal the model's positive contribution in representing the seasonal cycle of river flow, specifically for irrigation large-storage capacity reservoirs, as the model reproduces the seasonal shift between inflows and outflows caused by irrigation management rules with reasonable accuracy. For these reservoirs, the correlation between simulated and observed discharge increases from 0.49 (reference simulation) to 0.75 in the median. Regarding storage volumes, correlation reaches 0.74 in the median. As an example, the Gonzalez Lacasa reservoir located in the Ebro basin shows typical results in Fig. 4. In fact, this irrigation reservoir stores incoming water during winter/early spring, which explains a lower outflow than inflow (in blue and orange, respectively, for observed discharges, panels (a) and (b) in Fig. 4) between October and March and releases the water in summer when there is insufficient water to supply all crops' irrigation needs. The period of release, from April to September in this case (shown in the annual cycle, Fig. 4c), corresponds to the crop-growing period and therefore to high irrigation water needs. As a result, the discharge seasonal curve is shifted and the maximum monthly discharge, in this example, is reached in July instead of April. This management scheme is well reproduced by the DROP model, as the simulated outflows (shown in red) align well with the observations (in blue), and consequently so do the volume variations.

Figure 4 Simulation results for the Gonzalez Lacasa irrigation reservoir within the Ebro River basin over the period [1979–1994]. The monthly time series of dam releases and storage volumes are shown in panels (a) and (c), respectively (observations are indicated in blue and simulations in red). Their corresponding mean annual cycles are shown in panels (b) and (d), respectively.

The improvement of the C_2M distribution for each category of dams, in terms of main purpose and relative capacity, is shown in Fig. 5. Note that, for irrigation reservoirs, the improvement rate reaches 80 % in terms of the median and 123 % for those considered of high relative capacity (here c≥0.5) and which are fully demand-controlled. The flow improvement rates for the rest of the high relative capacity reservoirs are dispersed but remain positive at the median despite the simplistic approach of DROP. For low relative capacity reservoirs (here c<0.5), and independent of the management purpose, the model's contribution is almost null since the reservoirs are considered “run-of-river” and the influence of the inflow regime is predominant.

Figure 5 Distribution of river flow, C_2M, improvement rates with the DROP model based on their main purpose (columns) and relative capacity (rows). The x-axis labels show the number of reservoirs in each category, and the y axis shows C_2M improvement rates.

4.2 Results of the sensitivity analysis

The distributions of first-order Sobol indices for each parameter calculated at each of the 107 irrigation and 108 non-irrigation reservoirs are shown in Fig. 6 with a box plot. Overall, it emerges from Fig. 6 that the most influential parameter is c_threshold. In fact, based on the definition of S1, ∼48 % of the total variance in C_2M is attributed to c_threshold alone within irrigation reservoirs. In non-irrigation reservoirs, this parameter accounts for ∼74 % of C_2M variance. The M parameter is ranked second in irrigation reservoirs and accounts for 15 % alone in the median of all of the variance, followed by d_max with an S1 index of 0.03. The parameter controlling the month for which the operational year starts, m_start, alone has very little influence on the overall outflow C_2M variance, although the effect is slightly more noticeable in the non-irrigation reservoirs. α and b are considered the least influential parameters for most reservoirs.

Figure 6 Distribution of first (S1)-order Sobol indices in the modeled reservoirs for each parameter: (a) in irrigation and (b) in non-irrigation reservoirs.

To better illustrate how each parameter individually affects the model outputs, the Gonzalez Lacasa reservoir (GRanD ID 2699; c=0.65) is set as an example. Given the sensitivity analysis results, a screening step is added on c_threshold, M and α separately and the rest of the parameters are set to their default values. The means of monthly outflows over the simulation period (1979–1994) are shown in Fig. 7.

Figure 7 Example of seasonal pattern sensitivity of Gonzalez Lacasa reservoir outflow to three of the DROP model parameters: (a) c_threshold, (b) M and (c) α. The remaining parameters were at their default values (see Table 2).

Regarding the first parameter and for values of c_threshold below the relative capacity of the reservoir (lower than 0.65 for the Gonzalez Lacasa case), the reservoir is considered with high storage capacity, and so monthly releases are completely driven by demand, resulting in a seasonal shift between inflow and outflow and a peak of discharge in July in order to meet the irrigation needs. Conversely, when the c_threshold is higher, the dam is considered to have a low storage capacity, which reduces the buffering effect and increasingly aligns the simulated release curve (in red) with the seasonal trend cycle of inflow (in orange), as shown in Fig. 7a. The buffer role of the reservoir is therefore conditioned by the value of the c_threshold.

The influence of M on the minimum release is shown in Fig. 7b. The higher the value of M, the greater the part of mean inflow set as minimum release and the lower the remaining part of water dedicated to meeting irrigation needs during peak demand. Since the Gonzalez Lacasa reservoir is considered to have a relatively large storage capacity (with c_threshold=0.5 as the default value) among the considered reservoirs, the outflow will still follow the seasonal curve of water demand, but M will influence the release peak by extending or flattening the outflow curve to maintain a minimum release level required through the year. For run-of-river dams, M is irrelevant since the release follows the monthly inflow. d_max also controls the variation of outflows over the year in the same way, but only when the DPI ratio is low. Its impact remains limited beyond the (1−M) threshold, as the release curve is fixed to maintain the required minimum outflow.

The parameter α, on the other hand, does not affect the long-term mean pattern of release. In fact, α only operates on outflows on an inter-annual basis, through K_y, to offset the excess or shortage of stored water from one year to the other, especially when the critical filling zones are reached (reservoir in a dead storage zone or overflowing), in order to bring the reservoir water state to the ideal filling curve over the long term (Fig. 7c). The same behavior is noted for m_start, which is not shown in Fig. 7.

The c_threshold, α and m_start parameters have the same effect on non-irrigation reservoirs. The influence of each of the three parameters on [1979–2013] monthly mean outflows from the “Alcantara II” reservoir (GRanD ID 2800) in the Tagus River basin, which has a relative capacity of 0.6 and is mainly used for electricity generation, is shown in Fig. 8. For this specific reservoir, the simulated period includes several wet years, and there were periods when water flowed over the spillways during winter, which explains the alignment of the outflow with the inflow during this season.

Figure 8 An example of the seasonal pattern sensitivity of the Alcantara II hydroelectricity reservoir outflow to three of the DROP model parameters: (a) c_threshold, (b) α and (c) m_start. The remaining parameters were set at default values (see Table 2).

Figure 9 Distribution of first (S1)- and total (ST)-order Sobol indices in the modeled reservoirs for each parameter: (a) in irrigation and (b) in non-irrigation reservoirs.

The distributions of total sensitivity indices ST (in grey), alongside S1 (in white), of each parameter are shown in Fig. 9. A significant gap is observed between the first- and total-order index distributions. This confirms the non-negligible effect of the parameters' interactions on the output variance involving mainly c_threshold. In the median, based on the ST definition, ∼62 % of C_2M variance in irrigation reservoirs and ∼87 % in non-irrigation reservoirs are attributed exclusively to c_threshold and its interactions with other parameters: it is indeed the most important parameter of the DROP model. α and b, on the other hand, were in the median low values, which makes them the least important parameters.

Figure 10 Distribution of second-order Sobol indices (S2) in the modeled reservoirs for each pair of parameters: (a) in irrigation and (b) in non-irrigation reservoirs.

The distributions of second-order Sobol indices for each parameter couple for the irrigation and non-irrigation reservoirs considered separately are shown in Fig. 10. The results reveal that the parameter with the most interactions overall with other parameters is c_threshold. This is mainly due to its position at the end of the model chain, where provisional releases are corrected. For the rest of the parameters, the interaction between m_start and α is more marked in non-irrigation reservoirs since the scheme is simplistic and the outflow depends mainly on the K_y annual coefficient that is driven by both m_start and α. In irrigation reservoirs, the coefficients d_max and M interact because they control irrigation water demand and the reservoir water storage allocated to it. As they are both involved in the reservoir water balance, they jointly control the outflow. The b coefficient meanwhile interacts only with c_threshold when computing the demand-controlled release ratio R calculation. The total influence of b, mostly resulting from the interaction with c_threshold, is only seen in reservoirs with low relative capacities.

Since c is a defining characteristic of each reservoir and is involved in the outflow computation, the first- and total-order Sobol indices were rearranged according to their reservoir-correspondent c values, here grouped into six ranges to simplify the presentation of the results to better identify the impact on parameter ordering and interactions (Fig. 11).

Figure 11 Distribution of first- (“S1”) and total-order (“ST”) Sobol indices of parameters according to relative capacity: (a, c) in irrigation and (b, d) in non-irrigation reservoirs.

For low values of c (less than 0.01), c_threshold is the only relevant parameter. The total influence of b, on the other hand, is mostly related to its interaction with c_threshold. The remaining parameters are negligible since their first-order Sobol indices are almost zero and the total order is very low. For reservoirs with a medium storage capacity ($\mathrm{0.01}\le c\le \mathrm{1}$), a significant part of total variance is due to parameter interactions, apart from b, as all Sobol indices increase significantly from first to total order. This is even more noticeable when $\mathrm{0.1}\le c\le \mathrm{0.3}$. Above a relative capacity of 0.3, the M parameter in irrigation reservoirs gains in importance and increases interactions with the c_threshold as c is bigger. For non-irrigation reservoirs of the same category, c_threshold takes on more importance as the storage capacity increases (less spreading and very high S1). The remaining parameters lose relative importance. For the very high relative capacity values, α and m_start are almost irrelevant, regardless of the operating objective. Finally, for nearly all reservoirs combined (except those with very small values of c), the b coefficient has almost no influence and can therefore be set to a nominal value.

When using Sobol indices, the representation of model output uncertainty is limited to the variance only, which is not fully representative of all the statistical characteristics (or moments) of the C_2M distribution. Using the distribution function instead provides a complete description of uncertainty in the model output. Here, we evaluate the sensitivity of parameters by assessing their influence on the entire C_2M distribution without reference to a specific moment of the output (Chun et al., 2000; Borgonovo, 2007, moment-independent methods). This method is used as a validation step for the overall conclusions found with the variance-based Sobol method. The deviation of the C_2M distribution function caused by two different parameters c_threshold and α over the irrigation reservoirs is shown in Fig. 12. Here, the same samples generated for Sobol's index calculation were used. First, the unconditional probability distribution function (PDF) of the model output C_2M is obtained when all input parameters are randomly sampled from their distributions. Then, for each input parameter, conditional PDFs of C_2M are computed following different value ranges from the total parameter distribution: here, five ranges are considered. The importance of each parameter is proportional to the magnitude of the conditional PDF deviation from the unconditional one.

Figure 12 Shifts in C_2M' probability distribution function (PDF) depending on c_threshold (a) and α (b) range values, within irrigation reservoirs. The unconditional PDFs of the DROP model output C_2M, obtained when all input parameters are randomly sampled, are displayed in black. The conditional PDFs are shown in color depending on the parameter value range.

The large dispersion of the conditional PDFs (colored) shows the strong impact of c_threshold in the model outlet in Fig. 12a. All PDFs are nearly aligned for α, and thus the value of this parameter has very little significance for C_2M distribution. Moment-independent sensitivity indicators like “Kolmogorov”, “Kuiper” and “Delta” (Borgonovo and Plischke, 2016) were used to measure the deviation from the unconditional PDF, and the results were consistent with the conclusions made with the Sobol index. It is also noted that DROP tends to have very poor performance scores for high values of c_threshold (above 2.76) where all the reservoirs are considered run-of-river and the release would be close to the lines alignments with inflow as shown in Fig. 12. This figure indirectly demonstrates the reservoir model contribution in terms of improving river flow representation in anthropized basins where the magnitude and seasonal flow dynamics are significantly altered by large storage capacity reservoirs, especially those for irrigation purposes.

5 Discussion

5.1 Limits of the DROP model scheme

The DROP model outputs are affected by several uncertainties linked to the model inputs and the algorithm. The parameterized model is based on a generic scheme of reservoir operation, which inevitably implies a simplification in terms of water release. Concerning the representation of operational purposes, the algorithm fails to differentiate between other purposes than irrigation and considers that a constant release is expected by the rest of the reservoirs. In addition, only the main objective of the reservoir is represented, and the releases of the multi-purpose reservoirs are not entirely represented since their management rules are more complex to describe. This explains the poor performance of the model at the level of reservoirs that are used for irrigation, though it is not their main objective: DROP is very simplistic for the rest of its possible purposes, so that the seasonal shift in water discharge is not always well reproduced. This is the case of the reservoir “Los Bermejales”, for example, a multi-purpose reservoir located in the basin of the Guadalquivir River, southern Spain, which is mainly used for water supply, but it is also used to meet irrigation water demands (Fig. 13).

Figure 13 Example of time series (a) and monthly means (b) of simulated releases from Los Bermejales reservoir, a multi-purpose reservoir primarily used for water supply and which also operates to meet irrigation water demands. The reservoir model failed to reproduce observed irrigation releases as it does not consider secondary objectives.

In addition, the model computes releases independently on each reservoir. The cascades of reservoir operations, which can be coordinated with each other, are thus not captured. More specific studies on multi-objective reservoirs (Wu and Chen, 2012; Wang et al., 2019) and multi-reservoir systems have been conducted (Chang et al., 2014; Tan et al., 2017; Rougé et al., 2021), but they all remain complex and are reservoir-specific. These methods were evaluated only at the local scale and are very difficult to extend to a global scale because they need a significant amount of observed input data and require detailed operating rule knowledge.

Regarding the representation of releases from non-irrigation reservoirs more generally, the scheme remains very simplistic since release policy is not driven by physical processes. In fact, operation rules of these types of reservoirs involve complex socioeconomic and political factors that are different in each country. Simulating other management purposes is mainly based on optimization algorithms, as is the case for hydro-power dam releases, for instance, where the objective functions are economically oriented (i.e., to maximize energy production; Moeini et al., 2011; Feng et al., 2017; Chong et al., 2021). These methods remain very reservoir-specific and are currently deemed to be too complex to be applicable at a large scale.

The model provides a relatively good performance in representing irrigation reservoir operations because of its physical approach that links water releases to crop water demands. However, some simplifications are to be noted which could be improved in the future. The irrigation water demand estimation is based on the irrigation scheme in ISBA LSM, which has its own limitations (Druel et al., 2022). Also, in this version of the reservoir scheme, water demands for each irrigation reservoir are reduced to considering pixels that are downstream of the reservoir at a given maximum distance d_max. This creates inconsistencies because water demand is not linked to the reservoir water storage capacity, which leads for some reservoirs to much higher water demands compared to water supply. Moreover, here there is no proportionality rule set between reservoirs with common irrigation grid cells. They are recorded on multiple reservoirs when their command areas are overlapping because the model runs on each reservoir independently. We end up repeatedly counting shared pixels, and this leads to overestimation of water demands at the level of each reservoir. According to results from the sensitivity analysis, parameter d_max alone (ranging from 1 to 250 km) does not have much influence on outflow variance, and this is even more true when irrigation demands are considered excessive compared to inflow (DPI $\ge \mathrm{1}-M$), because in that case the amount of irrigation demand is not more significant in outflow computing; only the seasonal variation defines the reservoir release curve. Zhou et al. (2021) suggested an efficient way to overcome this issue by defining a least-cost adduction network, based on Portoghese et al. (2013) and Neverre et al. (2016), to connect each irrigation grid cell to a unique abstraction point, either a river or a reservoir, by using topographic information, distance and upstream areas of the river abstraction points. This approach has the advantage of considering not only the downstream grid cells, but all the surrounding ones. Most importantly, it ensures that reservoir command areas do not overlap anymore and that each pixel irrigation demand is only counted once. The method is implemented in the routing model of ORCHIDEE (ORganizing Carbon and Hydrology In Dynamic EcosystEms; Nguyen-Quang et al., 2018) and can be easily implemented in other routing models. It is also interesting at this stage to account for all possible sources of water withdrawals, including underground water, canals, but also the abstractions made directly from the reservoirs' storage, in order to have a more realistic representation of reservoir releases.

Another aspect which is not explicitly simulated in the model is water abstraction. In this study, abstractions are taken into account indirectly since the inflow is reconstructed from observations at the inlet of each reservoir. However, once the DROP model is implemented in a hydrological model, the tributary inflows will correspond to river flows simulated by the routing model, and therefore there should be a deterioration of the performance index on discharge with an error spreading along the anthropized rivers. However, by coupling the above with a model that takes irrigation into account, such as the new version of the ISBA land surface model, for example (Druel et al., 2022), the water releases from the reservoir model can be linked to crop irrigation needs, and thus river water withdrawals can be represented as well as those taken directly from the reservoirs.

5.2 Contributions of the sensitivity analysis to a clearer understanding of the DROP model

The sensitivity analysis has revealed the most influential parameters and those that can be set using predefined values without impacting the model output uncertainty distribution. It emerged that c_threshold is the most influential parameter in representing reservoir releases, and this result is consistent with the analysis of Shin et al. (2019) on the role of the release ratio R where c_threshold is involved. Indeed, the aforementioned study showed the positive contribution of this formulation to the R ratio (R=α c, i.e., $c_{\mathrm{threshold}}=\mathrm{1}/\mathit{\alpha }$) and its optimization in improving the representation of release and stabilizing water storage in the simulated reservoirs, especially those with $\mathrm{0.21}\le c\le \mathrm{1.18}$ (the values of c for which 4c²=α c and $c=\mathrm{1}/\mathit{\alpha }$, α being set at 0.85).

Hanasaki et al. (2006) conducted a sensitivity test on the α parameter for a case study of a large relative capacity hydroelectric reservoir (c=2.28). The model was tested with four different values of α (0.65, 0.75, 0.85, and 0.95) and found that this parameter had a low impact on the simulated release but a high impact on the simulated storage, except when the reservoir is full. α showed high sensitivity to releases when events with water passing over spillways were more frequent. Their conclusions concerning the sensitivity to outflows are also in line with the results found here.

Actually, the sensitivity analysis undertaken in this study focused on the model sensitivity in the average representation of releases rather than on the filling levels of reservoirs since the focus of the study is on flow representation and the effect of anthropogenic factors in altering the flow dynamics along the rivers. If we were to focus on water resource availability and water management issues, the sensitivity analysis should also focus on the uncertainty in representing water storage levels in the reservoirs, considering the C_2M over volumes as the variable of interest. In this configuration, as shown in Fig. 14, the parameters α and m_start, alongside c_threshold, will emerge as the most influential parameters with significantly high sensitivity indices of both first and total orders.

Figure 14 Distribution of first- (“S1”) and total-order (“ST”) Sobol indices of parameters based on C_2M variance over volumes: (a, c) in irrigation and (b, d) in non-irrigation reservoirs.

Both α and m_start directly affect the bias in the filling curve representation. In the ideal case where the reservoir state is not on the boundary conditions, the filling curve is vertically shifted following the value of α and the chosen month in order to bring the reservoir volume at the beginning of each operating year in line with the ideal filling rate α C, as shown in Fig. 15, without changing the seasonal release pattern. The larger the seasonality of reservoir water levels, the greater the effect of both parameters, which is very noticeable in large-storage capacity reservoirs used for water supply and hence explains the two parameters' interaction with c_threshold.

Figure 15 Example of seasonal pattern sensitivity of Gonzalez Lacasa (a, b) and Alcantara II (c, d) reservoir volumes to two of the DROP model parameters: (a, c) α and (b, d) m_start. The remaining parameters were set at default values (see Table 2).

6 Conclusions

In this paper, a global parameterized model, DROP, was reconstructed based on the Hanasaki et al. (2006) generic scheme to represent reservoir releases in Spain. Results reveal the positive contribution of the model in representing the seasonal cycle of discharge and storage variation, specifically for irrigation large-storage capacity reservoirs, as the model succeeds in reproducing the seasonal shift between inflows and outflows, improving river flow representation (C_2M improvement rate) by 123 % in the median. The results also provide a further validation of the Hanasaki et al. (2006) formulation on 215 reservoirs in Spain where reliable observed data are readily accessible. While Spain represented an idealized case study in terms of data availability, different remote sensing data will allow the extension of the model to any other river basin over the globe, more specifically to those which are ungauged, since the reservoir model relies only on the GRanD reservoir database (which is global). This work is in preparation for the upcoming SWOT wide-swath altimetry mission (Biancamaria et al., 2016), which will provide the data necessary to make improved global-scale river and reservoir storage and flow estimates. Furthermore, the results highlight the importance of incorporating reservoir operations into large-scale hydrological models for a more realistic representation of river flows and thus the water mass exchange with oceans and the atmosphere. The physical approach of DROP is consistent with that of the ISBA-CTRIP LSM RRM (land surface–river routing model): irrigation demands used as input to the reservoir model can be simulated by the irrigation module recently integrated into ISBA (Druel et al., 2022), and CTRIP already includes a lake model, “MLake” (Guinaldo et al., 2021), that a priori models inland water bodies at a global scale, calculates mass balance and lake outflow at the global scale, and provides the foundation for integrating human reservoir operations. The next step is to implement the DROP model in MLake and create a link between the two anthropization models by coupling the new versions of ISBA (irrigation) and CTRIP (reservoirs).

The sensitivity analysis, based on Sobol's method, was conducted on the C_2M representation of the mean seasonal outflow patterns. The results show that the most important parameter overall is c_threshold. M and d_max are ranked second, in the median, in terms of irrigation reservoir release representation, and their importance is linked to the reservoir relative capacity. α and m_start have less influence on both types of reservoir outflow seasonal dynamics but are important if the focus is on reservoir water storage values. It has also been proven that the significance of the model's parameters varies according to the range of reservoir relative capacities being studied. The results represent an essential step to further improve either river flow modeling or reservoir water storage, through calibration schemes and assimilation of new remote sensing products, by targeting the most influential reservoir model parameters. A supplementary application of this work is provided as a calibration study of the model over Spain (see Supplement). This work provides future studies with a fully generalized parameterization of the reservoir scheme along with a deeper insight into the way each of the parameters influences the model outputs and how their importance changes depending on the reservoir characteristics and the output variable of interest.

Overall, integrating this reservoir model into LSM RRMs, which are in turn coupled to climate and Earth system models (such as the CNRM-CM and CNRM-ESM; Voldoire et al., 2019; Séférian et al., 2019), will provide a major advance in understanding past climate reanalysis and will enable a more realistic representation of future scenarios under climate change.