Malak Sadki
email: malak.sadki@meteo.fr

Simon Munier

Aaron Boone

Sophie Ricci

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

Keywords: reservoir modeling, large-scale hydrological modeling, sensitivity analysis, Sobol method

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.

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 7,700 km 3 and 8,300 km 3 [START_REF] Vörösmarty | Anthropogenic sediment retention: major global impact from registered river impoundments[END_REF][START_REF] Downing | The global abundance and size distribution of lakes, ponds, and impoundments[END_REF]2011)), the total global impounded water may even exceed 10,000 km 3 [START_REF] Chao | Impact of artificial reservoir water impoundment on global sea level[END_REF]. 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 1,200 to 2,120 km 3 ), and they account for approximately 20% of average annual river flow to the oceans (40,000 -45,500 km 3 /year) [START_REF] Baumgartner | The world water balance: mean annual global, continental and maritime precipitation and run-off[END_REF][START_REF] Oki | Global hydrological cycles and world water resources[END_REF][START_REF] Syed | Satellite-based global-ocean mass balance estimates of interannual variability and emerging trends in continental freshwater discharge[END_REF][START_REF] Lehner | Highresolution mapping of the world's reservoirs and dams for sustainable river-flow management[END_REF]. 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 [START_REF] Dynesius | Fragmentation and flow regulation of river systems in the northern third of the world[END_REF][START_REF] Revenga | Prospects for monitoring freshwater ecosystems towards the 2010 targets[END_REF][START_REF] Grill | An index-based framework for assessing patterns and trends in river fragmentation and flow regulation by global dams at multiple scales[END_REF][START_REF] Grill | Mapping the world's free-flowing rivers[END_REF].

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 [START_REF] Haddeland | Anthropogenic impacts on continental surface water fluxes[END_REF][START_REF] Hanasaki | A reservoir operation scheme for global river routing models[END_REF][START_REF] Döll | Global-scale analysis of river flow alterations due to water withdrawals and reservoirs[END_REF][START_REF] Biemans | Impact of reservoirs on river discharge and irrigation water supply during the 20th century[END_REF][START_REF] Shin | High-resolution modeling of reservoir release and storage dynamics at the continental scale[END_REF][START_REF] Gutenson | Comparison of generalized non-data-driven lake and reservoir routing models for global-scale hydrologic forecasting of reservoir outflow at diurnal time steps[END_REF]. 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 [START_REF] Pokhrel | Model estimates of sea-level change due to anthropogenic impacts on terrestrial water storage[END_REF][START_REF] Zhao | Identifying effective forecast horizon for real-time reservoir operation under a limited inflow forecast[END_REF][START_REF] Wada | High-resolution modeling of human and climate impacts on global water resources[END_REF][START_REF] Frederikse | The causes of sea-level rise since 1900[END_REF]. 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 processbased 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., [START_REF] Young | Finding reservoir operating rules[END_REF][START_REF] Raman | Deriving a general operating policy for reservoirs using neural network[END_REF] to very sophisticated models based on machine learning and artificial intelligence techniques, such as neural-network-based methods (e.g., [START_REF] Maier | Neural networks for the prediction and forecasting of water resources variables: a review of modelling issues and applications[END_REF][START_REF] Razavi | Adaptive neural networks for flood routing in river systems[END_REF][START_REF] Ehsani | A neural network based general reservoir operation scheme[END_REF][START_REF] Coerver | Deduction of reservoir operating rules for application in global hydrological models[END_REF]. 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 [START_REF] Gutenson | Comparison of generalized non-data-driven lake and reservoir routing models for global-scale hydrologic forecasting of reservoir outflow at diurnal time steps[END_REF]. The best-known schemes in this category are those developed by [START_REF] Hanasaki | A reservoir operation scheme for global river routing models[END_REF] and [START_REF] Haddeland | Anthropogenic impacts on continental surface water fluxes[END_REF], which are inflow-and demand-based models as presented by [START_REF] Yassin | Representation and improved parameterization of reservoir operation in hydrological and land-surface models[END_REF]. These two generic models have been implemented in several global hydrological and water management models (e.g., WaterGAP [START_REF] Döll | Global-scale analysis of river flow alterations due to water withdrawals and reservoirs[END_REF]), VIC [START_REF] Haddeland | Anthropogenic impacts on continental surface water fluxes[END_REF], H08 [START_REF] Hanasaki | An integrated model for the assessment of global water resources-Part 1: Model description and input meteorological forcing[END_REF], PCR-GLOBWB [START_REF] Van Beek | Global monthly water stress: 1. Water balance and water availability[END_REF]).

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 [START_REF] Batalla | Reservoir-induced hydrological changes in the Ebro River basin (NE Spain)[END_REF][START_REF] López-Moreno | Dam effects on droughts magnitude and duration in a transboundary basin: The Lower River Tagus, Spain and Portugal[END_REF][START_REF] Lorenzo-Lacruz | The impact of droughts and water management on various hydrological systems in the headwaters of the Tagus River (central Spain)[END_REF]. Spain is, in fact, among the top 10 dam-building countries with more than 1,064 dams (ICOLD, 2020). The study of [START_REF] Grill | Mapping the world's free-flowing rivers[END_REF], 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., [START_REF] Zouhri | Handling the impact of feature uncertainties on SVM: a robust approach based on Sobol sensitivity analysis[END_REF] or civil engineering (e.g., [START_REF] Zamanian | Significant variables for leakage and collapse of buried concrete sewer pipes: A global sensitivity analysis via Bayesian additive regression trees and Sobol'indices[END_REF]. [START_REF] Pianosi | Sensitivity analysis of environmental models: A systematic review with practical workflow[END_REF] 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 [START_REF] Saltelli | Global sensitivity analysis: the primer[END_REF][START_REF] Saltelli | The cautious modeller: craftsmanship without wizardry[END_REF]. It can help to improve parameter calibration efficiency and avoid overparameterization (e.g., [START_REF] Shin | Using a global sensitivity analysis to estimate the appropriate length of calibration period in the presence of high hydrological model uncertainty[END_REF]Tang et al., 2007a, b). It is also an efficient tool to better understand the model structure [START_REF] Saltelli | Global sensitivity analysis: the primer[END_REF], its uncertainties (e.g., [START_REF] Pheulpin | Uncertainty quantification and global sensitivity analysis with dependent inputs parameters: Application to a basic 2D-hydraulic model[END_REF] and the dominant processes under various conditions (e.g., [START_REF] Huang | Global sensitivity analysis for a prediction model of soil solute transfer into surface runoff[END_REF][START_REF] Zhang | Sobol ′ 's sensitivity analysis for a distributed hydrological model of Yichun River Basin, China[END_REF]. 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 [START_REF] Liu | Uncertainty in hydrologic modeling: Toward an integrated data assimilation framework[END_REF][START_REF] Abdolghafoorian | Uncertainty quantification in land surface hydrologic modeling: Toward an integrated variational data assimilation framework[END_REF]. 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 [START_REF] Herrera | Parameter estimation and uncertainty analysis in hydrological modeling[END_REF]. A few sensitivity tests on the Hanasaki scheme can be found in the literature (e.g., [START_REF] Hanasaki | A reservoir operation scheme for global river routing models[END_REF][START_REF] Shin | High-resolution modeling of reservoir release and storage dynamics at the continental scale[END_REF], 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 [START_REF] Hanasaki | A reservoir operation scheme for global river routing models[END_REF], 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.

DROP

The parameterized DROP model has been developed based on the [START_REF] Hanasaki | A reservoir operation scheme for global river routing models[END_REF] scheme. The model works at the level of each reservoir individually and is based on the mass balance (as shown below in Equation 1) of each reservoir to calculate the release at its outlet.

dV dt = Q in -Q out (1)
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 Equation 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 .

Q in = (P -E) × A reservoir + Q trib + R d ± Q gw (2)
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 Figure 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 [START_REF] Hanasaki | A reservoir operation scheme for global river routing models[END_REF] 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 Figure 1.

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 f irst,y and the long-term target storage S ideal :

K y = S f irst,y /S ideal (3)
where S ideal = α×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 [START_REF] Hanasaki | A reservoir operation scheme for global river routing models[END_REF]. This annual release coefficient is the one used to weight the provisional releases calculated afterwards depending on whether the reservoir has more (Ky>1) or less (Ky<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 ′ m , 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):

r ′ m = i mean (4)
For an irrigation reservoir, monthly releases are proportional to water demands. They are parameterized here following [START_REF] Shin | High-resolution modeling of reservoir release and storage dynamics at the continental scale[END_REF] as:

r ′ m =      i mean (M + (1 -M ) d m d mean,y if DPI > 1 -M i mean + d m -d mean,y if DPI ⩽ 1 -M (5)
where d m and d mean,y are, respectively, monthly and annual mean water demands. DPI (demand per inflow) is, as introduced by [START_REF] Shin | High-resolution modeling of reservoir release and storage dynamics at the continental scale[END_REF], 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 [START_REF] Hanasaki | A reservoir operation scheme for global river routing models[END_REF] and [START_REF] Döll | Global-scale analysis of river flow alterations due to water withdrawals and reservoirs[END_REF], and to 0.1 by [START_REF] Biemans | Impact of reservoirs on river discharge and irrigation water supply during the 20th century[END_REF] and [START_REF] Shin | High-resolution modeling of reservoir release and storage dynamics at the continental scale[END_REF].

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 [START_REF] Hanasaki | A reservoir operation scheme for global river routing models[END_REF] (where a crop growth model is considered to calculate the irrigation demand), here the distribution of irrigated areas is based on ECOCLIMAP SG [START_REF] Calvet | L'apport de la télédétection spatiale à la modélisation des surfaces continentales[END_REF][START_REF] Druel | Implementation of a new crop phenology and irrigation scheme in the ISBA land surface model using SURFEX_v8.1[END_REF] which is used by the irrigation module in the ISBA LSM [START_REF] Druel | Implementation of a new crop phenology and irrigation scheme in the ISBA land surface model using SURFEX_v8.1[END_REF] to compute irrigation demands for 5×5 km-resolution grid cells (see section 3). The irrigation water demands are aggregated within the command area for each reservoir. Unlike the [START_REF] Hanasaki | A reservoir operation scheme for global river routing models[END_REF] 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 :

r m =    K y × r ′ m if c ⩾ c threshold (1 -R) × i m + R × K y × r ′ m if 0 ⩽ c < c threshold ( 6
)
where c is the relative capacity of a reservoir and is defined as the ratio between storage capacity C and the long-term mean 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 [START_REF] Hanasaki | A reservoir operation scheme for global river routing models[END_REF] and [START_REF] Biemans | Impact of reservoirs on river discharge and irrigation water supply during the 20th century[END_REF]. [START_REF] Shin | High-resolution modeling of reservoir release and storage dynamics at the continental scale[END_REF], 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 1/α and 1. The monthly releases are also weighted by the annual coefficient, K y from Equation 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 [START_REF] Hanasaki | A reservoir operation scheme for global river routing models[END_REF] and the [START_REF] Shin | High-resolution modeling of reservoir release and storage dynamics at the continental scale[END_REF] 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 equation ( 5) are monthly and annual mean irrigation water demands, respectively. Industrial and domestic water demands are not taken into account, unlike in [START_REF] Hanasaki | A reservoir operation scheme for global river routing models[END_REF] -M , set to 0.1 in the [START_REF] Biemans | Impact of reservoirs on river discharge and irrigation water supply during the 20th century[END_REF] and [START_REF] Shin | High-resolution modeling of reservoir release and storage dynamics at the continental scale[END_REF] schemes and 0.5 in [START_REF] Hanasaki | A reservoir operation scheme for global river routing models[END_REF], is considered as a parameter in this version, keeping the same notation used in [START_REF] Shin | High-resolution modeling of reservoir release and storage dynamics at the continental scale[END_REF].

-A more explicit parameterization of the demand-controlled release ratio R is provided here compared to what is proposed in [START_REF] Shin | High-resolution modeling of reservoir release and storage dynamics at the continental scale[END_REF]. The generalized formulation of R in Equation 7 highlights the last two parameters of the present model version. This enables one to distinguish the dual role of c threshold in Equation 6 from that of the coefficient b, which only describes the transition made between a demand-controlled reservoir and an input-controlled reservoir. [START_REF] Hanasaki | A reservoir operation scheme for global river routing models[END_REF] and [START_REF] Shin | High-resolution modeling of reservoir release and storage dynamics at the continental scale[END_REF] have set (c threshold , b) to (0.5,2) and (1/α,1), respectively.

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

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.

Sensitivity analysis : Sobol's method

The Sobol sensitivity analysis method [START_REF] Sobol | Sensitivity analysis for non-linear mathematical models[END_REF] 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 1 , X 2 , ...X p of size p, the total variance of Y , D(Y ), is decomposed as follows [START_REF] Zhang | Sobol ′ 's sensitivity analysis for a distributed hydrological model of Yichun River Basin, China[END_REF]):

D(Y ) = i D i + i<j D ij + i<j<k D ijk + ... + D 12..p (8)
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 [START_REF] Hoeffding | class of statistics with asymptotically normaldistribution[END_REF].

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 parameters, all orders included. The first-, second-and total-order Sobol sensitivity measure formulas are :

S1 i = D i D (9) S2 ij = D ij D (10) ST i =1 - D ∼i D (11)
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 [START_REF] Herman | SALib: An open-source Python library for Sensitivity Analysis[END_REF]. The parameter samples are generated following quasi-random sequences [START_REF] Saltelli | Making best use of model evaluations to compute sensitivity indices[END_REF]

Dataset

The data series required as inputs to the model and the ones used to validate the outputs are taken from the Spain database. -The Global Reservoir and Dam (GRanD) database: from which the general characteristics of dams are taken [START_REF] Lehner | Highresolution mapping of the world's reservoirs and dams for sustainable river-flow management[END_REF]. Version 1.3 published in 2019 includes 7,320 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 3 in 2016, which exceeds the mean annual streamflow of the eight major rivers of the Iberian Peninsula (55,850 hm 3 /year) [START_REF] Lorenzo-Lacruz | Recent trends in Iberian streamflows (1945-2005)[END_REF]. 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 [START_REF] Druel | Implementation of a new crop phenology and irrigation scheme in the ISBA land surface model using SURFEX_v8.1[END_REF]. 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-day 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í and standard deviation of 441 hm 3 . Using the latitude and longitude information, the chosen reservoirs were located on the new 1/12°-resolution river network derived from the CTRIP river routing model [START_REF] Munier | River network and hydro-geomorphology parametrization for global river routing modelling at 1/12°resolution[END_REF]. 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.

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., [START_REF] Vanderkelen | Evaluating a reservoir parametrization in the vector-based global routing model mizuRoute (v2. 0.1) for Earth system model coupling[END_REF]. In this study, the water balance (Equation 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 [START_REF] Decharme | Recent changes in the ISBA-CTRIP land surface system for use in the CNRM-CM6 climate model and in global off-line hydrological applications[END_REF], evaporation is computed by FLake, a module representing energy balance in lakes [START_REF] Le Moigne | Impact of lake surface temperatures simulated by the FLake scheme in the CNRM-CM5 climate model[END_REF], 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 [START_REF] Munier | River network and hydro-geomorphology parametrization for global river routing modelling at 1/12°resolution[END_REF].

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 (Equation 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 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. 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.

Sensitivity analysis implementation

A sensitivity analysis with respect to the six parameters was conducted on the performance of a Nash-Sutcliffe efficiency (Nash and Sutcliffe, 1970, NSE,) bounded version, called C 2M [START_REF] Mathevet | A bounded version of the Nash-Sutcliffe criterion for better model assessment on large sets of basins[END_REF], on outflows using the Sobol method. In fact, the NSE values in some reservoirs were highly negative for some simulations, thus this metric was not suitable for a variancebased 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). The default values for α, M , c threshold and b are those considered in [START_REF] Hanasaki | A reservoir operation scheme for global river routing models[END_REF]. α-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 tests 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 nonirrigation 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 4,096 was used for this analysis in each category, resulting in 4,096×(2×6+2)=57,344 and 4,096×(2×4+2)=40,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.

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.

Reproducing the flow seasonal shift in irrigation reservoirs

Using the default parameterization with the parameters listed in Table 2, 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),(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, panel (b) in Fig. 4), corresponds to 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. 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)andwhicharef ullydemand-controlled.T hef lowimprovementratesf ortherestof theh 0.5), andindependentof themanagementpurpose, themodel ′ scontributionisalmostnullsincethereservoirsareconsidered"runof -river"andtheinf luenceof theinf lowregimeispredominant.

Results of the sensitivity analysis

The distributions of first-order Sobol indices for each parameter calculated at each of the 107 irrigation and the 108 nonirrigation 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.

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 2).

buffering effect and increasingly aligns the simulated release curve (in red) with the seasonal trend cycle of inflow (in orange), as shown in Fig. 7 (a). The buffer role of the reservoir is therefore conditioned by the value of the c threshold .

390

The influence of M on the minimum release is shown in Fig. 7(b). The higher the value of M , the greater the part of mean inflow set as minimum release and the lower is 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 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 dead storage zone or overflowing) in order to bring the reservoir water state to the ideal filling curve over the long term (Fig. 7 (c)). 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 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. 2).

The distributions of total sensitivity indices ST (in grey), alongside with S1 (in white), of each parameter are shown in 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 corresponding 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).

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 (0.01⩽c⩽ 1), a significant part of total variance is due to parameter interactions, apart from b, as all Sobol indices increase significantly from first to 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 Los Bermejales reservoir, for example, a multi-purpose reservoir located in the basin of the Guadalquivir river, south of Spain, which is mainly used for water supply but it is also used to meet irrigation water demands(Fig. 13).

In addition, the model computes releases independently on each reservoir. The cascade of reservoir operations, which can be coordinated with each other, are thus not captured. More specific studies on multi-objective reservoirs [START_REF] Wu | An operation-based scheme for a multiyear and multipurpose reservoir to enhance macroscale hydrologic models[END_REF]; [START_REF] Wang | An improved operation-based reservoir scheme integrated with Variable Infiltration Capacity model for multiyear and multipurpose reservoirs[END_REF] and multi-reservoir systems have been conducted [START_REF] Chang | Optimized cascade reservoir operation considering ice flood control and power generation[END_REF]; [START_REF] Tan | Derivation of optimal joint operating rules for multipurpose multi-reservoir water-supply system[END_REF]; [START_REF] Rougé | Coordination and control-limits in standard representations of multi-reservoir operations in hydrological modeling[END_REF], 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) [START_REF] Moeini | Fuzzy rule-based model for hydropower reservoirs operation[END_REF][START_REF] Feng | Deriving adaptive operating rules of hydropower reservoirs using time-varying parameters generated by the E n KF[END_REF]; [START_REF] Chong | Optimization of hydropower reservoir operation based on hedging policy using Jaya algorithm[END_REF]. 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 [START_REF] Druel | Implementation of a new crop phenology and irrigation scheme in the ISBA land surface model using SURFEX_v8.1[END_REF]. 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 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 [START_REF] Druel | Implementation of a new crop phenology and irrigation scheme in the ISBA land surface model using SURFEX_v8.1[END_REF], 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.

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 [START_REF] Shin | High-resolution modeling of reservoir release and storage dynamics at the continental scale[END_REF] 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 threshold = 1/α) and its optimization in improving the representation of release and stabilizing water storage in the simulated reservoirs, especially those with 0.21 ≤ c ≤ 1.18 (the values of c for which 4c 2 = α c and c = 1/α , α 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.

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 .

Conclusions

In this paper, a global parameterized model, DROP, was reconstructed based on the [START_REF] Hanasaki | A reservoir operation scheme for global river routing models[END_REF] This work is in preparation for the upcoming SWOT wide-swath altimetry mission [START_REF] Biancamaria | The SWOT mission and its capabilities for land hydrology[END_REF], 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 [START_REF] Druel | Implementation of a new crop phenology and irrigation scheme in the ISBA land surface model using SURFEX_v8.1[END_REF] and CTRIP already includes a lake model "MLake" [START_REF] Guinaldo | Parametrization of a lake water dynamics model MLake in the ISBA-CTRIP land surface system (SURFEX v8. 1)[END_REF]) 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 reservoirs 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 (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 influence 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; [START_REF] Voldoire | Evaluation of CMIP6 deck experiments with CNRM-CM6-1[END_REF][START_REF] Séférian | Evaluation of CNRM earth system model, CNRM-ESM2-1: Role of earth system processes in present-day and future climate[END_REF], will provide a major advance in understanding past climate reanalysis and will enable a more realistic representation of future scenarios under climate change. 2).

Figure 1 .

 1 Figure 1. Schematic representation of the DROP model showing its six parameters (in blue): mstart 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: dmax 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).

  annual inflow water volume (c = C/I 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 :

  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 × (2 × d + 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 data3.1 Spain river basins descriptionSpain has an estimated area of 505,983 km 2 (INE, 2022), representing more than 85% of the Iberian Peninsula (estimated area 583,254 km 2 ,[START_REF] Lorenzo-Lacruz | Recent trends in Iberian streamflows (1945-2005)[END_REF]. It has river basins with sizes ranging from a few square kilometers to more than 80,000 km 2 (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 1,200 dams mainly in the second half of the 20 th 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(IEA, , in 2020,),).

  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 Tran-sition, Spain). The national database includes the location and daily time series of discharge for 1,119 gauge stations and outflows from 347 reservoirs over the period 1900-2014.

  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 steps3.3.1 Cross-referencing global and local datasetsOut 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 3,200 hm 3 (the La Serena dam on the Zujar river in the Gardiana basin) with a mean of 236 hm 3

Figure 2 .

 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.

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

 3 Figure 3. Results of the DROP model contribution (default configuration) in river flow modeling at reservoir outlets : (a) C2M performance metrics of non-reservoir simulation (reference simulation where Qout = Qin) and (b) C2M 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.

Figure 4 .

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

Figure 5 .

 5 Figure 5. Distribution of river flow, C2M , 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, y-axis shows C2M improvement rates.

Figure 6 .

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

Figure 7 .

 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 Table2).

Figure 8 .

 8 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) mstart. The remaining parameters were set at default values (see Table2).

Fig. 9 .

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

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

  compared to inflow (DPI ≥ 1-M)becauseinthatcasetheamountof irrigationdemandisnotmoresignif icantinoutf lowcomputing, o costadductionnetwork, basedonP ortogheseet al. (2013)andN everreet al. (2016), toconnecteachirrigationgridcelltoauniqueabst

  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[START_REF] Hanasaki | A reservoir operation scheme for global river routing models[END_REF] 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 GRanD reservoir database (which is global).

Figure 10 .

 10 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

Figure 12 .

 12 Figure 12. Shifts in C2M ' probability distribution function (PDF) depending on c threshold (a) and α (b) range values, within irrigation reservoirs. The unconditional PDFs of the DROP model output C2M , 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.

Figure 13 .

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

Figure 14 .

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

Figure 15 .

 15 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) mstart. The remaining parameters were set at default values (see Table2).

Table 1 .

 1 Description of DROP model parameters.

	Parameter	Description
	α	Refers to the ideal filling rate of the reservoir
	dmax (km)	Defines the control area of the reservoir
	mstart	Refers to the first month of the operational year
	M	Defines the minimum release to be provided for environmental requirements
	c threshold	Threshold of the relative capacity ; influence of the inflow pattern on reservoir release
	b	Influence of the inflow pattern on reservoir release

Table 2 .

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

	Parameter	Default value Min value Max value Distribution
	α	0.85	0.6	0.95	Uniform
	dmax (km)	100	1	250	Uniform
	mstart (irrigation;other)	4;5	1	12	Discrete uniform
	M	0.5	0	1	Uniform
	c threshold	0.5	0.001	20	Logarithmic
	b	2	0.5	5	Uniform

Region in France.

Acknowledgements. We would like to thank Dr. Pere Quintana-Segui from the Ebro Observatory for providing the SAFRAN-based meteorological dataset for Spain which was used to force the ISBA LSM and compute irrigation demands. We also thank him for providing pre-processed products of the in situ natural and anthropized flow and volume observations over Spain.

This study is part of Malak SADKI's thesis work, co-funded by the French National Center for Space Studies (CNES) and the Occitania

Code and data availability. The DROP model and the sensitivity analysis codes are available on Zenodo. Post-processing codes are also available. All information can be found in the following repository : (https://doi.org/10.5281/zenodo.6389405; [START_REF] Sadki | Implementation and sensitivity analysis of a Dam-Reservoir OPeration model (DROP v1.0) over Spain -Supplement[END_REF].

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][1980][1981][1982][1983][1984][1985][1986][1987][1988][1989][1990][1991][1992][1993][1994] are shown in Fig. 7.

Regarding the first parameter, and for values of c threshold below the relative capacity of the reservoir (lower than 0.65 for the total order. This is even more noticeable when 0.1 ≤ c ≤ 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 [START_REF] Chun | An uncertainty importance measure using a distance metric for the change in a cumulative distribution function[END_REF]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.

The large dispersion of the conditional PDFs (colored) shows the strong impact of c threshold in the model outlet in Fig. 12 (a). All PDFs are nearly aligned for α, and thus the value of this parameter has very little significance on C 2M distribution.

Moment-independent sensitivity indicators like "Kolmogorov", "Kuiper" and "Delta" [START_REF] Borgonovo | Sensitivity analysis: a review of recent advances[END_REF] 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 up 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.

Discussion

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 Author contributions. MS, SM, AB and SR designed the study and determined the methodology. MS and SM developed the reservoir model and the sensitivity analysis algorithm. MS performed the analysis and wrote the original draft. All authors contributed to the editing and review of the paper.

Competing interests. The authors declare that they have no conflict of interest.