Benchmarking reservoir operation schemes for large-scale hydrological models

There are approximately 62 000 large dams worldwide that significantly alter the hydrological regimes of most major rivers. Despite their importance, reservoirs remain poorly represented in Large-Scale Hydrological Models (LSHMs) due to the complexity of human-driven operations and a widespread lack of observational records. Consequently, reservoir routines in LSHMs must balance structural simplicity with limited data requirements. In this study, we utilize the ResOpsUS dataset to benchmark four reservoir routines of increasing complexity: LISFLOOD, CaMa-Flood, mHM, and STARFIT. We evaluate these routines across 164 reservoirs in the United States and test which target variables are most informative for parameter estimation. Our results indicate that the mHM routine consistently achieves the highest performance; however, its dependence on site-specific demand data limits its applicability at the global scale. In contrast, the CaMa-Flood routine provides a robust compromise, significantly outperforming the linear logic of LISFLOOD and matching the performance of the more complex STARFIT routine. This suggests that increased model complexity does not always yield superior results. Crucially, we find that calibrating to reservoir storage is more informative than calibrating to outflow, as it effectively captures the dynamics of both state variables. This finding paves the way for the use of satellite-derived storage products in the calibration of LSHMs. The findings of this study have been implemented in the upcoming versions of the European and Global Flood Awareness Systems (EFAS v6 and GloFAS v5).

1 Introduction

The increasing global population and expanding economic activities have amplified the pressure on freshwater resources, making human interventions – such as the construction and operation of dams and reservoirs – fundamental elements of the terrestrial water cycle (Biemans et al., 2011; Busker et al., 2019; Sadki et al., 2023; Moreno-Rodenas et al., 2025). These structures are pivotal in fulfilling diverse societal needs, including domestic and industrial water supply, irrigation, hydropower production, and flood control (Mekonnen and Hoekstra, 2012; Abeshu et al., 2023; Sadki et al., 2023; Lehner et al., 2024; Shrestha et al., 2024; Moreno-Rodenas et al., 2025). The global proliferation of dams has been unprecedented over the last century (Lehner et al., 2024; Moreno-Rodenas et al., 2025). Globally, there are approximately 62 000 large dams (heights > 15 m), cataloged by the ICOLD (International Commission on Large Dams) (Moreno-Rodenas et al., 2025), with an estimated cumulative storage capacity ranging between 7420 km³ (GDW v1.0) (Lehner et al., 2024) and 8300 km³ (Sadki et al., 2023). This volume represents roughly 20 % of the average annual river discharge to the oceans and is four times the amount of water stored in river channels worldwide (Sadki et al., 2023).

The immense scale of this intervention profoundly impacts natural hydrology by modifying the magnitude, timing, and duration of streamflows (Shin et al., 2019; Salwey et al., 2024), resulting in the fragmentation of over 60 % of the world's largest rivers (Sadki et al., 2023). In the contiguous United States (CONUS), dams regulate all major rivers, providing a total storage capacity equivalent to roughly 75 % of the mean annual CONUS runoff (Steyaert et al., 2022). While beneficial for human use, these operations cause significant environmental trade-offs, including the homogenization of river dynamics (Scherer and Pfister, 2016; Lehner et al., 2024), alteration of sediment and nutrient transport (Shin et al., 2019; Yassin et al., 2019), and increased open-water evaporation, which can exacerbate water scarcity (Scherer and Pfister, 2016; Yassin et al., 2019; Liu et al., 2026).

Despite their importance, the accurate representation of complex reservoir operations remains a challenge in large-scale hydrological models (LSHMs) and land-surface models (Turner et al., 2021; Sadki et al., 2023; Salwey et al., 2024). Reservoir operations are intrinsically non-linear and governed by complex, human-driven decisions related to water demand, flood risk, and environmental constraints (Coerver et al., 2018; Hanazaki et al., 2022; Steyaert et al., 2022; Salwey et al., 2024; Shrestha et al., 2024). Consequently, models lacking adequate reservoir schemes often exhibit degraded predictive performance, especially in highly regulated basins (Turner et al., 2021; Salwey et al., 2024; Steyaert et al., 2025). Historically, global models have relied on generic operational rules utilizing static characteristics (e.g., capacity, primary purpose) from datasets like the Global Reservoir and Dam (GRanD) database (Lehner et al., 2011). However, these generic approaches often fail to capture the local operating behaviors necessary for simulating realistic daily releases and storage dynamics (Turner et al., 2021). Recent research has therefore pivoted toward exploiting newly available observational data and integrating data-driven techniques, such as machine learning or statistically inferred policies, to better reflect historical decision-making (Coerver et al., 2018; Turner et al., 2021; Steyaert et al., 2025).

This research aims to improve the reservoir module within the OS LISFLOOD hydrological model (Roo et al., 2000; van der Knijff et al., 2010; Burek et al., 2013; JRC, 2026c), a core component of the Global and European Flood Awareness Systems of the Copernicus Emergency Management Services (GloFAS and EFAS) (JRC, 2026a, b). The current LISFLOOD reservoir module employs an operation scheme based on three fixed storage limits (conservation, normal, and flood pools) without explicit classification of reservoir purpose (Zajac et al., 2017). Our objective is to identify a reservoir routine suitable for global application – one that balances portability with efficacy in data-scarce environments.

Specifically, this study presents a systematic benchmark comparing the performance of four distinct reservoir schemes within large-scale hydrological models, all of which are calibrated using in-situ records. We selected schemes from the literature that represent an increasing level of complexity and varying modeling philosophies: the current LISFLOOD storage-based parameterization (Zajac et al., 2017), its evolution as implemented in CaMa-Flood (Hanazaki et al., 2022), a demand-driven approach from mHM (Shrestha et al., 2024), and a seasonality-driven statistical approach, STARFIT (Turner et al., 2021). While not exhaustive, this selection encompasses the primary modeling paradigms currently used in the field. Other models, such as DZTR (Yassin et al., 2019) or SBTS (Shen et al., 2025), could be evaluated in future work. By assessing the capability of these schemes to simulate observed flow regulation and storage dynamics, this research provides critical guidance for the development of human -water interaction modules in the next generation of operational flood systems and LSHMs.

A novel aspect of this study is the exploration of a decoupled calibration strategy, where the reservoir routines are evaluated in isolation from the rest of the catchment hydrology. Rather than calibrating reservoir parameters simultaneously with the rest of the hydrological model against downstream river discharge – a process often plagued by equifinality (Beven and Freer, 2001) – we model reservoir management as a standalone system. To achieve this, we selected reservoirs with observed inflow, outflow, and storage records. This approach ensures that reservoir parameters are sensitive to actual operations rather than compensating for upstream hydrological biases or being masked by gauge observations far downstream.

However, the practical implementation of this strategy faces a significant data bottleneck. While in situ records of reservoir outflow provide the most direct target for calibration, such data are frequently inaccessible or classified as proprietary. Consequently, river discharge measured at gauging stations downstream of the reservoir is often used as a proxy for reservoir operations. This approach inherently assumes that total reservoir outflow is equivalent to river discharge – an assumption that fails when significant volumes are diverted for hydropower, irrigation, or municipal supply. To overcome these limitations, satellite-derived products offer a promising alternative (Schwatke et al., 2015; Pekel et al., 2016; Schwatke et al., 2020; Donchyts et al., 2022; Khandelwal et al., 2022; Hao et al., 2024; Hou et al., 2024). These products offer the potential to monitor volume fluctuations in ungauged regions, making them crucial for future global model calibration and data assimilation (Hanazaki et al., 2022; Shen et al., 2025; Steyaert et al., 2025).

A calibration based on satellite products would target reservoir level or storage, instead of outflow (or river discharge as a proxy variable). To assess the feasibility of using satellite products, we compare the performance of the model when calibrated against three distinct targets: (i) reservoir release only, (ii) reservoir storage only, and (iii) a combined objective function utilizing both variables. By evaluating these configurations, this study aims to assess whether satellite-inferred storage can serve as a robust alternative to ground-based outflow records, thereby enabling improved reservoir modeling in ungauged or data-scarce regions.

The benchmarking is supported by the ResOpsUS dataset (Steyaert et al., 2022; Steyaert and Condon, 2024), the first comprehensive, national-scale inventory providing historical daily reservoir operations for 679 major dams across the CONUS. ResOpsUS includes records of reservoir inflow, level, storage and outflow (not all variables for all reservoirs), providing the necessary variables for our benchmark. We leverage this data-rich environment to test the models and calibration targets, establishing a proof-of-concept for operational environments where such granular data is unavailable.

In summary, this study benchmarks four reservoir modeling paradigms across the CONUS domain to address three primary objectives: (i) to identify the most robust reservoir routine for continental or global hydrological models; (ii) to determine the most suitable target variable for calibration: storage, outflow, or both; and (iii) to evaluate the feasibility of using remotely-sensed storage data in the calibration of large-scale models.

2 Data

2.1 Reservoir Selection and Filtering

To evaluate reservoir routines in a data-rich context, we utilized the ResOpsUS dataset (Steyaert et al., 2022; Steyaert and Condon, 2024), which provides daily time series of reservoir operations for 679 reservoirs across the CONUS. From this initial pool, we applied a rigorous filtering process to ensure the suitability of the reservoirs for LSHM:

1. Data Availability: We selected only those reservoirs containing concurrent daily records for all three primary variables: inflow (I), storage (V), and outflow (Q).

2. Physical Dimensions: Reservoirs were required to have a minimum catchment area of 50 km² and a minimum storage capacity (V_max) of 10 hm³.

3. Hydrological Impact: Large-scale hydrological models such as GloFAS or EFAS are primarily focused on modelling river streamflow. For the sake of parsimony, we only want to include in those systems “disruptive” reservoirs – those that significantly alter the downstream flow regime. Shrestha et al. (2024) devised two thresholds that identify disruptive reservoirs using the Degree of Regulation (DOR; Eq. 1a) and the Degree of Disruptivity (DOD; Eq. 1b). Following their recommendations, only reservoirs with DOR≥0.08 and DOD≥0.06 m were selected:

   $$\begin{array}{}\text{(1a)}& {\displaystyle }\mathrm{DOR}={\displaystyle \frac{V_{max}}{{\stackrel{\mathrm{‾}}{I}}_{\mathrm{yr}}}}[-]\text{(1b)}& {\displaystyle }\mathrm{DOD}={\displaystyle \frac{V_{max}}{A_{\mathrm{c}}}}\left[\mathrm{m}\right]\end{array}$$

   where V_max is the reservoir storage capacity [m³], ${\stackrel{\mathrm{‾}}{I}}_{\mathrm{yr}}$ is the average annual inflow [m³ yr⁻¹], and A_c is the catchment area [m²].

4. Record Length: To ensure sufficient data for robust model calibration, we removed reservoirs with time series shorter than four years.

5. Water Balance Integrity: Preliminary analysis revealed large biases between average inflow and outflow in several reservoirs. Such discrepancies may stem from poor data quality or significant consumptive use (e.g., water diversions) not explicitly represented in our modeling framework. To ensure a closed water balance for benchmarking, we removed reservoirs where the absolute bias exceeded 30 %.

Applying these criteria resulted in a final selection of 164 reservoirs. Their geographical distribution, categorized by catchment area and storage capacity, is depicted in Fig. 1.

Figure 1 Distribution of the 164 reservoirs from the ResOpsUS dataset used in the benchmarking. Colors represent the main reservoir use and point size the degree of regulation.

2.2 Dataset Integration and Attribute Enrichment

To satisfy the input requirements of the diverse reservoir routines, we enriched the ResOpsUS records with supplementary hydrometeorological time series, and catchment and reservoir attributes.

As detailed in Sect. 3.1, the routines require time series of direct precipitation and evaporation. These were extracted for each reservoir location from the ERA5 reanalysis dataset (Hersbach et al., 2020).

Static reservoir and dam attributes were primarily sourced from the Global Reservoir and Dam dataset (GRanD) (Lehner et al., 2011). While we included all GRanD attributes in our dataset, four were critical for the simulations: storage capacity (V_max), maximum surface area (A_max), dam crest elevation (Z_max), and dam height (H_dam). A comparison of these attributes against observed time series revealed inconsistencies in some attribute values. To rectify this, we cross-referenced GRanD with the Global Dam Watch (GDW) dataset (Lehner et al., 2024). For instances where the attributes remained inconsistent, we manually updated the records using official data from US authorities, including the U.S. Army Corps of Engineers' National Inventory of Dams (NID), the U.S. Bureau of Reclamation (USBR), and the U.S. Geological Survey (USGS) (US Army Corps of Engineers, 2025; US Bureau of Reclamation, 2025; US Geological Survey, 2025). Similar errors were reported by Shrestha et al. (2024).

For the purpose of parameter regionalization or the development of data-driven models, we further included catchment-scale characteristics derived from GloFAS surface fields (Choulga et al., 2024) and climate indices derived from ERA5 (Hersbach et al., 2020).

We named the enriched dataset Reservoir Operations US and CAtchment and Reservoir Static attributes (ResOpsUS+CARS). It is available via https://doi.org/10.5281/zenodo.15978041 (Casado Rodríguez et al., 2025). The dataset structure follows the standards defined by the CARAVAN initiative (Kratzert et al., 2023), ensuring it is formatted for immediate use in deep learning and large-sample hydrological studies.

The specific time series and static attributes from ResOpsUS+CARS utilized in this study are summarized in Table 1. These include model forcing inputs, variables for calibration and evaluation, and geometric parameters used for the estimation of reservoir area or elevation.

Table 1 Data included in the ResOpsUS+CARS dataset used in this study.

Sources: [1] ResOpsUS; [2] ERA5; [3] GRanD.

3 Methods

3.1 Reservoir models

Each of the reservoir routines evaluated in this study is governed by the water mass balance equation (Eq. 2). At each time step Δt, the model calculates the reservoir storage at the end of the interval (V_t+1) based on the preceding storage state (V_t), the water surface area (A_t), and three forcing variables: inflow (I_t), precipitation (P_t), and evaporation (E_t):

$$\begin{array}{}\text{(2)}& V_{t+\mathrm{1}}=V_t+\left(P_t-E_t\right)\cdot A_t+\left(I_t-Q_{\mathrm{t}}\right)\cdot \mathrm{\Delta }t\end{array}$$

where Q_t is the reservoir outflow (release), the fundamental difference between the four reservoir schemes.

Figure 2 Simplified reservoir shape adopted in this study, including all the model input and output variables.

To represent the reservoir's bathymetry, we assume an inverted half-pyramid geometry (Fig. 2), following the approach of Liebe et al. (2005) and subsequent studies (Shin et al., 2019; Sadki et al., 2023; Shrestha et al., 2024).This geometric approximation allows for the continuous estimation of the water surface area (Eq. 3a), which is required for the mass balance, and the water level (Z_t; Eq. 3b), provided that the dam crest elevation and the dam height are known. We evaluate the validity of this geometric assumption by comparing estimated levels against observed records in a subset of reservoirs where such data is available.

$$\begin{array}{}\text{(3a)}& {\displaystyle }& {\displaystyle }{A}_t=A_{max}({\displaystyle \frac{V_t}{V_{max}}}{)}^{\mathrm{2}/\mathrm{3}}\text{(3b)}& {\displaystyle }& {\displaystyle }{Z}_t=Z_{max}-H_{\mathrm{dam}}\left[\mathrm{1}-({\displaystyle \frac{V_t}{V_{max}}}{)}^{\mathrm{1}/\mathrm{3}}\right]\end{array}$$

3.1.1 LISFLOOD

The current LISFLOOD reservoir routine (Burek et al., 2013; Zajac et al., 2017) operates as a piecewise linear storage-outflow relationship. The release is determined by the storage zone in which the reservoir resides at any given time step: conservative, normal or flood.

$$\begin{array}{}\text{(4)}& \begin{array}{rl}& Q_{\mathrm{t}}=\\ & \left\{\begin{array}{ll}{Q}_c& \mathrm{if}{V}_t<\mathrm{2}{V}_{\mathrm{c}}\\ Q_c+(Q_n-Q_c)\frac{V_t-\mathrm{2}{V}_{\mathrm{c}}}{V_{\mathrm{n}}-\mathrm{2}{V}_{\mathrm{c}}}& \mathrm{if}\mathrm{2}{V}_{\mathrm{c}}\le V_t<V_{\mathrm{n}}\\ Q_n& \mathrm{if}{V}_{\mathrm{n}}\le V_t<V_{\mathrm{n}}^{\prime }\\ Q_n+(Q_{\mathrm{f}}-Q_n)\frac{V_t-V_{\mathrm{n}}^{\prime }}{V_{\mathrm{f}}-V_{\mathrm{n}}^{\prime }}& \mathrm{if}{V}_{\mathrm{n}}^{\prime }\le V<V_{\mathrm{f}}\\ \begin{array}{rl}max(& \frac{V_t-V_{\mathrm{f}}}{\mathrm{\Delta }t},\\ & min(Q_{\mathrm{f}},max(k{I}_t,Q_n\left)\right))\end{array}& \mathrm{if}{V}_t>V_{\mathrm{f}}\end{array}\right.\end{array}\end{array}$$

where V_c, V_n, $V_{\mathrm{n}}^{\prime }$ and V_f are the conservative storage, the lower and upper bounds of the normal storage zone, and the flood storage limit, respectively. Q_c, Q_n and Q_f denote the outflow values associated with these limits. The release coefficient k allows for larger outflow than inflow when the reservoir exceeds the flood limit; while this was fixed at 1.2 in the original implementation, we treated it as a calibration parameter in this study (Table 2).

Table 2 Calibration parameters in the LISFLOOD reservoir module.

^a Parameter identical in LISFLOOD and CaMa-Flood. ^b Parameter adjusted in a standard LISFLOOD calibration.

In the usual LISFLOOD calibration, only two reservoir parameters were tuned: the normal outflow (Q_n) and the upper limit of the normal storage zone ($V_{\mathrm{n}}^{\prime }$). The remaining breakpoints are kept at their default values to limit the dimensionality of the global calibration, which involves 14 parameters per basin.

However, in this benchmarking study, we allow maximum flexibility to the LISFLOOD routine by calibrating the six parameters listed in Table 2. This expanded parameterization serves two purposes. On one hand, it ensures a fair comparison with more flexible, data-driven schemes. On the other hand, because our decoupled calibration focuses exclusively on reservoir parameters rather than the entire hydrological model, we can afford the increased dimensionality to explore the routine's maximum performance potential.

A fundamental characteristic of this scheme is that it creates a unique (one-to-one) relationship between storage and outflow (Fig. 3). The outflow is strictly a function of storage, independent of seasonality, inflow or water demand. The only exception occurs in the flood zone (V_t≥V_f), where the inclusion of the inflow (I_t) in the release logic allows for deviations from this static behavior. This makes the LISFLOOD scheme an ideal baseline for assessing whether more complex, dynamic routines are necessary to capture observed reservoir behavior.

Figure 3 Storage–outflow relationships for the (a) LISFLOOD and (b) CaMa-Flood reservoir routines. The curves illustrate the calibrated operational rules for the Yellowtail dam (GRanD ID 355). Grey points represent observed daily records, highlighting the spread of historical operations.

3.1.2 CaMa-Flood

The reservoir module in the Catchment-based Macro-scale Floodplain model (CaMa-Flood) (Hanazaki et al., 2022) represents an evolution of the LISFLOOD approach by introducing inflow-dependency into the operational logic. Unlike the static LISFLOOD scheme, CaMa-Flood distinguishes between two operational modes based on the current inflow (I_t).

When I_t is below a threshold (Q_f), the outflow is modeled as a quadratic function of storage. This non-linear behavior restricts releases as the reservoir empties, effectively conserving water for future demand. Conversely, when I_t exceeds Q_f, the system transitions toward a linear reservoir behavior to manage high-flow conditions.

$$\begin{array}{}\text{(5)}& \begin{array}{rl}& Q_{\mathrm{t}}=\\ & \left\{\begin{array}{ll}\left\{Q_n\frac{V_t}{V_{\mathrm{f}}},V_t<V_{\mathrm{c}}\right.& \forall I_t\\ \left\{\begin{array}{ll}{Q}_c+(\frac{V_t-V_{\mathrm{c}}}{V_{\mathrm{e}}-V_{\mathrm{c}}}{)}^{\mathrm{2}}\mathrm{\Delta }Q& V_{\mathrm{c}}\le V_t<V_{\mathrm{e}}\\ Q_{\mathrm{f}}& V_t\ge V_{\mathrm{e}}\end{array}\right.& I_t<Q_{\mathrm{f}}\\ \left\{\begin{array}{ll}{Q}_c+\frac{V_t-V_{\mathrm{c}}}{V_{\mathrm{f}}-V_{\mathrm{c}}}\mathrm{\Delta }Q& V_{\mathrm{c}}\le V_t<V_{\mathrm{f}}\\ Q_{\mathrm{f}}+k\frac{V_t-V_{\mathrm{f}}}{V_{\mathrm{e}}-V_{\mathrm{f}}}(I_t-Q_{\mathrm{f}})& V_{\mathrm{f}}\le V_t<V_{\mathrm{e}}\\ I_t& V_t\ge V_{\mathrm{e}}\end{array}\right.& I_t\ge Q_{\mathrm{f}}\end{array}\right.\\ & \mathrm{where}{Q}_c=Q_n{\displaystyle \frac{V_{\mathrm{c}}}{V_{\mathrm{f}}}}\mathrm{and}\mathrm{\Delta }Q=Q_{\mathrm{f}}-Q_c\end{array}\end{array}$$

where V_e is the emergency storage limit – defined as the upper part of the flood zone – and k is the release coefficient that regulates outflows under flood conditions based on available storage capacity. Hanazaki et al. (2022) defined k as a fixed value derived from the total flood volume (Eq. 6a). Under that formulation, k becomes zero for high-regulation reservoirs, resulting in a constant release even as storage levels approach V_e, effectively neglecting the exhaustion of regulation capacity. To address this unrealistic behavior, we introduced a dynamic calculation of k at each time step based on the instantaneous reservoir volume (V_t, Eq. 6b). This modification ensures that k increases as storage approaches V_e, facilitating higher discharge to mitigate the risk of overtopping.

$$\begin{array}{}\text{(6a)}& {\displaystyle }k=max(\mathrm{1}-{\displaystyle \frac{V_{max}-V_{\mathrm{f}}}{A_{\mathrm{c}}\cdot \mathrm{0.2}}},\mathrm{0})\text{(6b)}& {\displaystyle }k=max(\mathrm{1}-{\displaystyle \frac{V_{max}-V_t}{A_{\mathrm{c}}\cdot \mathrm{0.2}}},\mathrm{0})\end{array}$$

While Hanazaki et al. (2022) utilized default parameters, we calibrated five parameters (Table 3) to ensure the routine has comparable degrees of freedom to the other routines. To maintain consistency with the LISFLOOD benchmark, three parameters (α, δ, ϵ) share identical definition across both models.

Table 3 Calibration parameters in the CaMa-Flood reservoir module.

^* Parameter identical in LISFLOOD and CaMa-Flood.

The primary advantage of the CaMa-Flood routine over LISFLOOD is its ability to represent a non-unique relationship between storage and outflow (Fig. 3). By incorporating I_t into the operational rules, the routine can better capture the observed dispersion in the storage-outflow values, providing a more realistic representation of human-regulated flow variability.

3.1.3 mHM

The reservoir module within the Multiscale Hydrological Model (mHM) (Samaniego et al., 2010; Kumar et al., 2013; Thober et al., 2019) represents a fundamental shift in modeling philosophy compared to the LISFLOOD and CaMa-Flood routines. Rather than treating outflow as a direct function of storage, the mHM routine determines releases primarily based on a water demand time series (D_t). While Shrestha et al. (2024) utilized random forests to predict D_t for individual reservoirs, we implement a seasonal approach to maintain model portability.

In our implementation, we derive an annual demand climatology by computing the mean observed outflow for each day of the year across the available records. To ensure this signal reflects demand rather than total outflow, we apply the water stress ($\stackrel{\mathrm{‾}}{D}/\stackrel{\mathrm{‾}}{I}$) as scaling factor, and smooth the signal using a 28 d moving average. This smoothed climatology serves as the seasonal demand signal D_t for each reservoir, introducing the temporal variability missing in static storage–discharge schemes.

• Reservoir release. The total daily release Q_t is partitioned into a component satisfying the hedged demand (${\widehat{D}}_t$) and a component representing the passthrough of inflow (I_t):

  $$\begin{array}{}\text{(7)}& Q_{\mathrm{t}}=\mathit{\rho }\cdot {\mathit{\kappa }}_t\cdot {\widehat{D}}_t+(\mathrm{1}-\mathit{\rho })\cdot I_t\end{array}$$

  The partition coefficient ρ is a function of the reservoir's degree of regulation (DOR). For highly regulated reservoirs (DOR≥α), the release is entirely demand-driven (ρ=1). In all other cases, the influence of demand relative to inflow is governed by the parameters α and β:

  $$\begin{array}{}\text{(8)}& \mathit{\rho }=min\left(\mathrm{1},{\left({\displaystyle \frac{\mathrm{DOR}}{\mathit{\alpha }}}\right)}^{\mathit{\beta }}\right)\end{array}$$

  The demand component is further moderated by a time-varying coefficient κ_t, which hedges the release based on the current storage level V_t relative to a target normal storage V_n:

  $$\begin{array}{}\text{(9)}& {\mathit{\kappa }}_t=({\displaystyle \frac{V_t}{V_{\mathrm{n}}}}{)}^{\mathit{\lambda }}=({\displaystyle \frac{V_t}{\mathit{\gamma }\cdot V_{max}}}{)}^{\mathit{\lambda }}\end{array}$$

  where γ defines the normal storage fraction and λ controls the sensitivity of the hedging to storage deficits.

• Demand hedging. The hedged demand ${\widehat{D}}_t$ is calculated based on the ratio between average annual demand ($\stackrel{\mathrm{‾}}{D}$) and inflow ($\stackrel{\mathrm{‾}}{I}$). The parameter ω represents the expected annual water surplus and is used to categorize the reservoir's water stress. If there is no water stress ($\stackrel{\mathrm{‾}}{D}/\stackrel{\mathrm{‾}}{I}<\mathrm{1}-\mathit{\omega }$), the daily demand (D_t) is supplied in addition to the average water surplus. If there is stress, the hedged demand is a combination of a fraction of the average inflow and a fraction of the daily demand.

  $$\begin{array}{}\text{(10)}& {\widehat{D}}_t=\left\{\begin{array}{ll}\stackrel{\mathrm{‾}}{I}-\stackrel{\mathrm{‾}}{D}+D_t& \mathrm{if}\frac{\stackrel{\mathrm{‾}}{D}}{\stackrel{\mathrm{‾}}{I}}<\mathrm{1}-\mathit{\omega }\\ \mathit{\omega }\cdot \stackrel{\mathrm{‾}}{I}+(\mathrm{1}-\mathit{\omega })\frac{\stackrel{\mathrm{‾}}{I}}{\stackrel{\mathrm{‾}}{D}}\cdot D_t& \mathrm{otherwise}\end{array}\right.\end{array}$$

• Calibration. Following the procedure established by Shrestha et al. (2024), we calibrated five parameters (Table 4). Other authors used default values of some of the parameters (Hanasaki et al., 2006; Shin et al., 2019; Sadki et al., 2023), which we have used for the simulation with default parameters.

Table 4 Calibration parameters in the mHM reservoir module.

Unlike the rule-based schemes, the mHM routine produces a more realistic storage–outflow relationship, as the release is decoupled from the instantaneous storage volume and instead follows the seasonal cycle of demand. While more flexible, the main limitation of the routine is how to estimate the demand signals.

3.1.4 STARFIT

The Storage Targets And Release Function Inference Tool (STARFIT) framework represents a data-driven paradigm that infers operational rules by learning the seasonality of reservoir storage and outflows (Turner et al., 2021). Unlike previous rule-based models, STARFIT employs harmonic functions to represent seasonal cycles. A key advantage of this formulation is its scale-independence; while the model parameters are typically fitted using weekly data to reduce noise, the harmonic frequency (ω) can be adjusted for daily simulations ($\mathit{\omega }=\mathrm{1}/\mathrm{365}$), ensuring consistency with the temporal resolution of the simulation.

All variables are standardized to facilitate regionalization and comparison across basins (Eq. 11). Storage is converted to reservoir filling, while inflow and outflow are expressed as normalized anomalies relative to the mean annual inflow ($\stackrel{\mathrm{‾}}{I}$):

$$\begin{array}{}\text{(11a)}& {\displaystyle }& {\displaystyle }{\widehat{V}}_t={\displaystyle \frac{V_t}{V_{max}}}\text{(11b)}& {\displaystyle }& {\displaystyle }{\widehat{Q}}_t={\displaystyle \frac{Q_{\mathrm{t}}-\stackrel{\mathrm{‾}}{I}}{\stackrel{\mathrm{‾}}{I}}}\text{(11c)}& {\displaystyle }& {\displaystyle }{\widehat{I}}_t={\displaystyle \frac{I_t-\stackrel{\mathrm{‾}}{I}}{\stackrel{\mathrm{‾}}{I}}}\end{array}$$

• Storage Normal Operating Range (NOR). STARFIT defines the operational target of a reservoir through a Normal Operating Range (NOR), bounded by upper (NOR_up) and lower (NOR_low) seasonal harmonics. These bounds represent the typical filling levels for any given time of the year:

  $$\begin{array}{}\text{(12)}& \begin{array}{rl}{\mathrm{NOR}}_{\mathrm{up},t}& =min(max(A+B\sin \mathrm{2}\mathit{\pi }\mathit{\omega }t+C\cos \mathrm{2}\mathit{\pi }\mathit{\omega }t,\\ & {\widehat{V}}_{min}),{\widehat{V}}_{max})\\ {\mathrm{NOR}}_{\mathrm{low},t}& =min(max(a+b\sin \mathrm{2}\mathit{\pi }\mathit{\omega }t+c\cos \mathrm{2}\mathit{\pi }\mathit{\omega }t,\\ & {\widehat{v}}_{min}),{\widehat{v}}_{max})\end{array}\end{array}$$

  where t is the time (day or week) of the year. The NOR has 10 parameters: six defining the harmonic curves and four capping the maximum and minimum filling (Table 5). To fit these curves, the model uses only the three most extreme (highest and lowest) observed storage for each calendar week (Fig. 4), ensuring the NOR captures the envelope of historical operations.

• Release function. The release policy is conditioned on the reservoir's position relative to the NOR. To maintain the reservoir within its seasonal bounds, STARFIT applies a three-zone operational policy:

  $$\begin{array}{}\text{(13)}& \begin{array}{rl}& Q_{\mathrm{t}}=\\ & \left\{\begin{array}{ll}min(Q_{min}+(Q_{\mathrm{NOR},t}-Q_{min})\frac{{\widehat{V}}_t}{{\mathrm{NOR}}_{\mathrm{low},t}},I_t)& \mathrm{if}{\widehat{V}}_t<{\mathrm{NOR}}_{\mathrm{low},t}\\ max(Q_{min},min(Q_{\mathrm{NOR},t},Q_{max}\left)\right)& \mathrm{if}{\widehat{V}}_t\in {\mathrm{NOR}}_t\\ Q_{\mathrm{NOR},t}+(Q_{max}-Q_{\mathrm{NOR},t})\frac{{\widehat{V}}_t-{\mathrm{NOR}}_{\mathrm{up},t}}{\mathrm{1}-{\mathrm{NOR}}_{\mathrm{up},t}}& \mathrm{if}{\widehat{V}}_t>{\mathrm{NOR}}_{\mathrm{up},t}\end{array}\right.\end{array}\end{array}$$

  We modified the definitions of the outflow for both below and above the NOR, as the original definitions produced in some cases noisy releases. In both cases, we have introduced a linear interpolation between a fixed minimum or maximum outflow (Q_min or Q_max) and the outflow associated with the respective NOR boundary ($Q_{{\mathrm{NOR}}_{\mathrm{low}},t}$ or $Q_{{\mathrm{NOR}}_{\mathrm{up}},t}$). Within the NOR, the standardized outflow (${\widehat{Q}}_t$) is modeled as the sum of an harmonic component (${\widehat{Q}}_{\mathrm{harm},t}$) and a linear component (${\widehat{Q}}_{\mathrm{lin},t}$) that accounts for the current storage and inflow:

  $$\begin{array}{}\text{(14)}& \begin{array}{rl}& {\widehat{Q}}_{\mathrm{NOR},t}={\widehat{Q}}_{\mathrm{harm},t}+{\widehat{Q}}_{\mathrm{lin},t}\\ & {\widehat{Q}}_{\mathrm{harm},t}=d\sin \mathrm{2}\mathit{\pi }\mathit{\omega }t+e\cos \mathrm{2}\mathit{\pi }\mathit{\omega }t\\ & +f\sin \mathrm{4}\mathit{\pi }\mathit{\omega }t+g\cos \mathrm{4}\mathit{\pi }\mathit{\omega }t\\ & {\widehat{Q}}_{\mathrm{lin},t}=h+i{\displaystyle \frac{{\widehat{V}}_t-{\mathrm{NOR}}_{\mathrm{low},t}}{{\mathrm{NOR}}_{\mathrm{up},t}-{\mathrm{NOR}}_{\mathrm{low},t}}}+j\cdot {\widehat{I}}_t\end{array}\end{array}$$

  The release harmonic allows for bimodal release patterns (e.g., two release peaks per year). The linear term is discarded if the linearity is weak (R²<0.2), reverting to a purely seasonal policy.

  The complete release model comprises 7 parameters (Table 5), fitted using the full historical record of weekly releases (Fig. 4). Additionally, the minimum and maximum releases (Q_min and Q_max) are estimated as quantiles of the observed daily inflows; we have changed the original values of these quantiles to 1 % and 99 %, respectively.

Figure 4 The STARFIT model for the Yellowtail dam (GRanD ID 355): (a) seasonal harmonic curves defining the storage normal operating range (NOR) for each calendar week; (b) the seasonal harmonic release; (c) the linear model of release residuals based on available storage (A_st) and standardized inflow.

Table 5 Parameters in the STARFIT reservoir module.

^a Upper case refers to the upper NOR, and lower case to the lower NOR. ^b Percentiles of the observed daily flows.

3.2 Experimental design

To address our research questions, we conducted four distinct simulation experiments designed to evaluate the routines under varying levels of data availability and calibration complexity.

The first run evaluates model performance using default parameters sourced from the literature (see Tables 2, 3, and 4). This experiment assesses the “out-of-the-box” reliability of the routines, which is particularly relevant for continental or global applications where historical operational data are unavailable.

The remaining three experiments involve local parameter optimization targeting different target variables:

• OBJ-Q. Univariate calibration of reservoir outflow. This mirrors the standard calibration procedure for operational systems like EFAS/GloFAS, where model parameters are tuned to match observed streamflow.

• OBJ-V. Univariate calibration of reservoir storage. This experiment evaluates the potential for training routines using satellite-derived storage products, which is vital for ungauged basins.

• OBJ-QV. Bivariate calibration of both storage and outflow. This multi-objective approach analyzes the trade-offs in model fidelity and identifies the information loss associated with single-variable calibration.

Parameter optimization was performed using the Shuffled Complex Evolution (SCEUA) algorithm (Duan et al., 1993, 1994), implemented via the Spotpy Python library (Houska et al., 2015). For each reservoir, the algorithm was executed for up to 1000 iterations using four complexes. Following the recommendations of Shen et al. (2022) for large-sample hydrology, we used the full length of the available daily time series for calibration to capture the widest possible range of hydrological conditions.

Performance was quantified using the modified Kling-Gupta Efficiency (KGE^′) (Kling et al., 2012), which decomposes performance into correlation, bias, and variability components. For the bivariate experiment, we integrated the individual KGE^′ scores into a single multi-objective metric:

$$\begin{array}{}\text{(15a)}& {\displaystyle }{\mathrm{KGE}}^{\prime }=\mathrm{1}-\sqrt{(\mathit{\alpha }-\mathrm{1}{)}^{\mathrm{2}}+(\mathit{\beta }-\mathrm{1}{)}^{\mathrm{2}}+(\mathit{\rho }-\mathrm{1}{)}^{\mathrm{2}}}\text{(15b)}& {\displaystyle }{\mathrm{KGE}}_{\mathrm{bivariate}}^{\prime }=\mathrm{1}-\sqrt{(\mathrm{1}-{\mathrm{KGE}}_Q^{\prime }{)}^{\mathrm{2}}+(\mathrm{1}-{\mathrm{KGE}}_V^{\prime }{)}^{\mathrm{2}}}\end{array}$$

where ρ is the Pearson correlation coefficient, α represents variability (ratio of coefficients of variation), and β is the bias (ratio of means), defined as:

$$\begin{array}{}\text{(16)}& {\displaystyle }\mathit{\beta }={\displaystyle \frac{{\mathit{\mu }}_{\mathrm{sim}}}{{\mathit{\mu }}_{\mathrm{obs}}}}\text{(17)}& {\displaystyle }\mathit{\alpha }={\displaystyle \frac{{\mathrm{CV}}_{\mathrm{sim}}}{{\mathrm{CV}}_{\mathrm{obs}}}}\end{array}$$

While the LISFLOOD, CaMa-Flood, and mHM routines are optimized using the iterative SCEUA algorithm, STARFIT employs a direct statistical fitting procedure. This approach does not require a heuristic search. Instead, the storage model parameters are estimated in a two-step process: an initial ordinary least squares (OLS) regression provides a first guess for the three harmonic parameters, followed by a non-linear optimization using the L-BFGS-B algorithm to finalize all five parameters, including the upper and lower capping values (${\widehat{V}}_{max},{\widehat{V}}_{min}$). In contrast, the parameters for the release model are derived exclusively through OLS regression.

Because the original STARFIT formulation inherently requires observations of both storage (to define the NOR) and outflow (to define the release policy), it is only evaluated within the OBJ-QV bivariate experiment. This ensures a consistent comparison between the statistically-derived STARFIT policies and the calibrated rule-based routines.

4 Results

4.1 Which target variable is most relevant for reservoir calibration?

Figure 5 illustrates the performance of the three rule-based reservoir routines (LISFLOOD, CaMa-Flood, and mHM) across the four experimental configurations. The results highlight a critical dependency between the calibration objective and the model's ability to represent the dual states of discharge and storage.

Figure 5 Comparison of the model performance depending on the variable(s) targeted in the calibration. The panels correspond to the three reservoirs routines that were calibrated: (a) LISFLOOD, (b) CaMa-Flood and (c) mHM. Within each panel, three groups of box plots show the performance in terms of outflow, storage, and both variables together. Four parameter sets are compared – default parameters, calibration of outflow (OBJ-Q), storage (OBJ-V) and both variables (OBJ-QV) – indicated by colors. The dots represent individual reservoirs, indicating the distribution of the overall performance.

The default parameterization (gray) provides a reasonably effective baseline for reservoir outflow across all three routines, with median KGE^′ values of 0.44, 0.52, and 0.60 for LISFLOOD, CaMa-Flood, and mHM, respectively. However, these default parameters fail to capture storage dynamics effectively; median KGE^′ scores for storage drop to 0.30, 0.36, and 0.29, respectively. The significantly larger interquartile range (IQR) in these scores indicates the high variability in storage performance across the 164 reservoirs.

As expected, each routine achieves its peak performance for a specific variable when that variable is the sole target of the calibration. However, there is a clear trade-off when we look at how the models perform on variables they weren't trained on. Calibrating outflow (OBJ-Q) severely degrades storage performance; in LISFLOOD and CaMa-Flood, this configuration yields the lowest median KGE^′ for storage among all experiments; in mHM, while the median is slightly higher than the default parameters, the dispersion between reservoirs is remarkably large. Conversely, calibrating to storage does not result in a similar degradation of outflow performance. In all three routines, OBJ-V maintains outflow KGE^′ values that are superior to the default parameterization.

The bivariate calibration (OBJ-QV) consistently achieves the highest combined performance across all routines. Notably, the performance of the OBJ-QV run is closely mirrored by the storage-only calibration (OBJ-V). These findings suggest that reservoir storage is a more informative variable than outflow for identifying model parameters.

While operational hydrology has traditionally relied on streamflow (outflow) for calibration, our results indicate that storage time series provide a more robust constraint on the internal state of the reservoir. This has significant implications for global hydrology, as it supports the feasibility of using satellite-derived storage products to calibrate reservoir routines in ungauged or data-poor regions (Shen et al., 2025).

4.2 Which is the best-performing reservoir routine?

Figure 6 compares the performance of the four routines across the different experimental configurations. Note that the outflow-only calibration (OBJ-Q) is excluded here for simplicity, as the previous section demonstrated that it is a suboptimal strategy for capturing overall reservoir behavior.

Figure 6 Comparison of the model performance depending on the reservoir routine. The panels correspond to three of the experiments: (a) default parameters, (b) calibration of storage (OBJ-V), (c) bivariate calibration of outflow and storage (OBJ-QV). Within each panel, three groups of box plots show the performance in terms of outflow, storage, and both variables together. Four reservoir routines are compared – LISFLOOD, CaMa-Flood, mHM and STARFIT – indicated by colors. The dots represent individual reservoirs, indicating the distribution of overall performance.

Among the rule-based schemes, the mHM routine consistently achieves the highest performance across nearly all scenarios. Its primary strength lies in its ability to represent seasonal operations through demand-driven logic. The only exception is the uncalibrated (default) simulation of storage, where CaMa-Flood shows a slight advantage. Note that the CaMa-Flood routine was specifically designed for global applications without site-specific tuning (Hanazaki et al., 2022).

The results also highlight the benefits of model complexity: the CaMa-Flood routine, which evolved from the original LISFLOOD scheme by adding quadratic storage logic and inflow-dependency, systematically outperforms LISFLOOD in every experiment. This suggests that even small increases in the physical complexity of the storage-discharge relationship can significantly improve model performance.

Finally, the bivariate calibration results provide the first direct comparison with the STARFIT model. Despite its high parameterization (up to 19 parameters) and seasonal flexibility, STARFIT does not outperform mHM. Its combined performance is generally comparable to CaMa-Flood. This suggests that while seasonal harmonics are valuable, the structured rules in mHM and CaMa-Flood may provide a more robust framework for continuous daily simulations in large-scale hydrological models.

For a detailed analysis of model performance across each routine and variable, Appendix B presents the decomposition of the KGE^′.

4.3 How accurate is the reservoir shape approximation?

To represent the elevation–storage–area relationship in our modeling, we adopted a simplified geometric assumption where all reservoirs are treated as an inverted half-pyramid (Liebe et al., 2005). Under this assumption, the reservoir geometry is defined by the maximum surface area (A_max) and the dam height (H_dam) (Eq. 3). This approximation directly influences the calculation of the instantaneous surface area, which is a critical term for estimating direct precipitation and open-water evaporation.

While in situ observations of reservoir area are often unavailable, we validated this geometric assumption by comparing time series of estimated and observed water elevation. Time series of observed storage were converted into water elevations via Eq. (3b). These estimates were then compared against in situ observations for a subset of 138 reservoirs where such data were available. The performance of this elevation estimation is summarized in Fig. 7.

Figure 7 Performance of the reservoir level estimation assuming an inverted half-pyramid reservoir shape. The swarm plot on the left shows the overall distribution, where each point represents a reservoir and the box plot the quartiles. The map on the right shows the geographical distribution of performance; the dot size indicates storage capacity.

The inverted half-pyramid assumption proved to be a robust approximation for the majority of the study sites. The median KGE^′ for reservoir level estimation is 0.72, with 75 % of the reservoirs achieving a KGE^′ at least 0.44. From the three components of KGE^′, the variability is the main cause of poor performance (median value of 1.25), whereas bias and correlation showed an overall good performance (median values of 1.00 in both cases). The accuracy of the level calculation is highly sensitive to the quality of the static attributes, specifically the crest elevation (Z_max) and dam height (H_dam) (Eq. 3b). As discussed in Sect. 2.2, errors in global databases such as GRanD or GDW can propagate into these elevation estimates; therefore, data cleaning and the use of locally validated attributes are essential for accurate water-level modeling.

5 Discussion

5.1 Balancing model complexity and operational feasibility

This study compared reservoir routines of varying complexity, from the simple storage-dependent logic of LISFLOOD to the demand-driven and statistically-inferred policies of mHM and STARFIT. Our results indicate that model performance generally scales with the number of environmental covariates included in the routine (e.g., inflow or demand).

While mHM emerged as the best-performing routine (Sect. 4.2), its dependence on a demand time series poses a significant operational challenge. In data-rich environments, demand can be inferred via machine learning (Shrestha et al., 2024), but in continental or global systems like EFAS or GloFAS, such data are often unavailable.

A similar constraint applies to STARFIT, although its performance was not markedly superior to the more parsimonious CaMa-Flood routine. A distinct advantage of STARFIT for global applications is its flexibility across different temporal resolutions. The harmonic functions can be fitted at coarse or irregular temporal resolutions (such as those in satellite products) and be executed at daily resolution by adjusting the frequency term (ω). For instance, Steyaert et al. (2025) recently demonstrated a global application of the STARFIT routine within the PCR-GLOBWB2 hydrological model by leveraging satellite-derived storage from the GloLAKES dataset (Hou et al., 2024) to define the storage normal operating range (NOR). To explore this potential, Appendix A provides a comparative analysis of the four reservoir routines when calibrated using satellite-derived storage estimates.

The CaMa-Flood routine has been selected for implementation in the forthcoming major evolution of the open-source LISFLOOD model (OS-LISFLOOD). This update will form the core of the upcoming EFAS v6 and GloFAS v5 operational systems (Salamon et al., 2025). These systems undergo regular updates encompassing updated meteorological forcing, revised surface fields, extended river discharge records and model developments. By evolving the reservoir logic to include inflow-dependency and quadratic storage–outflow relationships, we achieve a measurable gain in performance without increasing the data requirements for calibration. Furthermore, the robust default performance of the CaMa-Flood routine supports global applications in regions where historical records are sparse (Hanazaki et al., 2022).

5.2 Current limitations and the forecast challenge

A notable limitation across the benchmarked routines is the absence of explicit reservoir purpose (e.g., irrigation, hydropower, or flood control …) (Shen et al., 2025). Reservoir operation is fundamentally driven by water demand, which varies significantly by use case. For instance, hydropower reservoirs are often maintained at high levels to maximize hydraulic head, whereas flood control reservoirs are drawn down ahead of wet seasons to maximize storage capacity. Conversely, irrigation reservoirs follow a distinct seasonal filling and release cycle to bridge the gap between wet-season supply and dry-season demand.

While simple schemes like LISFLOOD or CaMa-Flood do not explicitly represent these operational rules, several studies have focused on purpose-specific modeling. This includes routines for water supply in the UK (Salwey et al., 2024), hydropower in France (Baratgin et al., 2024), and non-consumptive uses (Shrestha et al., 2024). Other approaches differentiate between irrigation and non-irrigation schemes (Hanasaki et al., 2006; Sadki et al., 2023) or apply distinct objective functions based on the reservoir's primary use (Haddeland et al., 2006). Although the reservoir routines evaluated in this study do not explicitly consider reservoir purpose, Appendix C shows the model performance disaggregated by the main reservoir use reported in GRanD.

In our comparison, the superior performance of mHM and STARFIT likely stems from their ability to capture these dynamics indirectly. mHM is a demand-driven routine that can infer specific usage patterns directly from demand time series, while STARFIT leverages observed seasonal patterns to learn operational behaviors. The fact that these two models outperformed LISFLOOD and CaMa-Flood, even if by a small margin, suggests that accounting for the temporal variability of human-driven demand is critical for improving global reservoir modeling.

Furthermore, these routines lack a forward-looking component, whereas in practice, reservoir operators rely on meteorological and hydrological forecasts across multiple timescales. While replacing current inflow with forecasted values could address this, it introduces a significant bottleneck for computationally demanding systems. The model would need to simulate the entire upstream catchment for the full forecast horizon before the reservoir release can be determined at any given node. While local or regional systems running simpler lumped models could potentially process these computations sequentially at each reservoir location, this remains a substantial implementation challenge for continental or global distributed systems.

5.3 The potential of satellite products

The scarcity of in situ records remains the primary hurdle for large-scale reservoir modeling (Shen et al., 2025). Satellite-derived products offering reservoir area, elevation, and storage anomalies provide an opportunity to bridge this gap (Zhao and Gao, 2018; Schwatke et al., 2020; Donchyts et al., 2022; Khandelwal et al., 2022; Hao et al., 2024; Hou et al., 2024). Two of the outcomes of this study specifically support the transition toward using satellite products. First, we found that reservoir storage is a more informative calibration target than reservoir outflow, as storage-based calibration yields models that perform well across both state variables (Sect. 4.1). This represents a “lucky coincidence” for the community, as many available satellite products provide direct estimates of reservoir storage or storage variations rather than discharge. Second, our validation of the inverted half-pyramid shape (Sect. 4.3) suggests that even simplified geometric assumptions are sufficient to translate satellite-derived levels into volumes with reasonable accuracy. However, a rigorous analysis of the feasibility of remotely sensed data should consider also spatial and temporal resolution, apart from accuracy.

Recent studies are already leveraging these synergies; for instance, Shen et al. (2025) estimated monthly target storage values for 1178 GRanD reservoirs by combining area time series from Global Water Watch (GWW) (Donchyts et al., 2022) with area-storage curves from GRDL (Hao et al., 2024). This monthly target storage was integrated into an evolved reservoir scheme – based on LISFLOOD and CaMa-Flood – to introduce seasonality into otherwise simple storage–outflow formulations.

The flexibility of STARFIT's harmonic formulation is particularly well-suited for satellite data, as emphasized by Steyaert et al. (2025). Because the frequency term in the harmonic curves allows the model to be fitted at coarser temporal resolutions (e.g., monthly satellite passes) and then run at a daily time step, it overcomes the temporal sampling limitations of many satellite products.

In line with Shen et al. (2025), we have conducted a preliminary analysis of the potential use of the (GWW) dataset (Donchyts et al., 2022), as detailed in Appendix A. Our findings indicate that while satellite-derived storage follows observed historical trends, a persistent bias often exists between in situ and remotely-sensed values. The primary bottleneck for the quality of these products is the availability of accurate elevation-storage-area curves. While several efforts have attempted to provide these curves globally (Yigzaw et al., 2018; Busker et al., 2019; Khazaei et al., 2022), they frequently show shortcomings when validated against in situ data.

5.4 Toward data-driven reservoir modeling

The recent revolution in deep learning has demonstrated that data-driven models can outperform traditional physical models in streamflow prediction (Nearing et al., 2024). To support a similar transition in reservoir modeling, we have formatted the ResOpsUS+CARS dataset following the CARAVAN standard (Kratzert et al., 2023). This extension of the original ResOpsUS dataset (Steyaert et al., 2022) integrates the reservoir operations with meteorological time series, the reservoir attributes from GRanD/GDW, catchment characteristics from GloFAS static maps, and climate indices from ERA5. By consolidating these diverse covariates – encompassing reservoir morphology, catchment physiography, seasonality, meteorology and climate – we provide the necessary information for deep learning models to identify the complex drivers of reservoir release.

Reservoir operations are the product of complex human-natural interactions that process-based routines struggle to capture entirely. While we acknowledge that deep learning offers a path toward superior performance, these models are notoriously data-hungry. This highlights the ongoing necessity of curated, multi-national observational datasets – not just for training, but for the rigorous benchmarking of the hybrid physical-statistical models and the development of global remotely-sensed approaches.

6 Conclusions

This study provided a comprehensive benchmarking of four reservoir routines (LISFLOOD, CaMa-Flood, mHM, and STARFIT) evaluated across 164 reservoirs in the United States using the ResOpsUS and GRanD datasets. By comparing different calibration strategies and model structures, we draw the following conclusions.

Reservoir storage is the most informative variable for calibration. Our results demonstrate a significant asymmetry in model identifiability. While calibrating routines to match observed outflow (the standard practice in many hydrological models) leads to a poor representation of internal reservoir states, calibrating to observed storage yields robust performance for both storage and outflow. This finding supports the shift toward using satellite-derived storage products for model parameterization, particularly in ungauged basins.

mHM offers superior performance, but CaMa-Flood provides the best operational balance. While the demand-driven logic of the mHM routine consistently demonstrated superior performance, its reliance on site-specific demand data constrains its applicability in continental and global systems where such information is often scarce. The CaMa-Flood routine, which will be integrated into the OS-LISFLOOD model, represents an optimal compromise: it provides significant improvements over the original linear storage logic without requiring additional input variables or data-intensive calibration.

The inverted half-pyramid geometric assumption is a robust approximation. Validating reservoir level estimations against in situ observations confirmed that a simplified inverted half-pyramid shape is sufficient for large-scale modeling. This validates the use of maximum area and dam height as sufficient proxies to link storage, area, and elevation, which is critical for estimating evaporative losses and integrating altimetry or reservoir area data.

STARFIT's flexibility is a key asset for remote sensing integration. Although STARFIT did not outperform the rule-based mHM or CaMa-Flood routines in daily simulations, its harmonic formulation is uniquely capable of bridging temporal gaps in observations. Its ability to be fitted at monthly scales and run at daily resolutions makes it a powerful tool for leveraging current and future satellite altimetry missions.

Looking forward, the integration of reservoir modeling into global hydrology must move beyond the “data-poor” paradigm. The creation of standardized, multi-national datasets like the ResOpsUS+CARS is a step toward this goal. By providing these data in the CARAVAN format, we aim to facilitate the development of hybrid models that combine the physical consistency of process-based routines with the predictive power of deep learning. Such advancements will be essential for improving the accuracy of operational flood and drought forecasting systems in an increasingly managed global water cycle.

Appendix A: Evaluation of GWW storage estimates

To assess the feasibility of using Global Water Watch (GWW) (Donchyts et al., 2022) for future global-scale reservoir calibration, we evaluated the accuracy of its volume estimates against our observational sample. Although GWW provides reservoir area for nearly all reservoirs in our dataset, storage data were available for only 12 of them. Figure A1 illustrates the performance and geographical distribution of this subset, comparing the observed storage time series against GWW estimates.

Despite the limited sample size, GWW storage estimates demonstrate relatively high skill, with a median KGE^′ of 0.65 and 75 % of the reservoirs exceeding a score of 0.51. In terms of decomposition, GWW tends to overestimate total storage (median bias of 1.12 with an IQR of 0.13) while slightly underestimating variability (median of 0.91, IQR 0.37). Notably, the temporal dynamics are captured accurately, reflected by a median correlation of 0.88 (IQR 0.10).

To evaluate the feasibility of utilizing satellite products in future applications, we conducted a fifth experiment (OBJ-V-GWW) wherein the reservoir schemes were calibrated against the storage estimates from GWW. This setup mirrors the univariate storage calibration experiment (OBJ-V), with the distinct modification that the optimization target is a satellite-derived product rather than in situ observations. Because the STARFIT formulation requires both storage and outflow time series for its fitting procedure, this specific model was parameterized using observed outflow alongside the satellite-estimated storage. Figure A2 illustrates the performance of this experiment validated against in situ observations (i.e., evaluated against the actual reservoir storage records rather than the satellite data used during optimization).

Calibration using the satellite product consistently yields performance improvements relative to the default parameter baselines. The LISFLOOD and CaMa-Flood routines appear to be particularly sensitive to the constraints of the satellite data. Both models experience a substantial surge in storage simulation performance, with median KGE^′ values increasing by 0.57 and 0.48 for LISFLOOD and CaMa-Flood, respectively. Interestingly, this enhancement does not impact outflow performance in LISFLOOD, whereas it noticeably degrades it in CaMa-Flood (a loss of 0.18 in median KGE^′).

Conversely, the satellite-based calibration yields negligible improvements for the mHM routine (gains of only 0.07 and 0.09 in median KGE^′ for storage and outflow, respectively), resulting in poor overall performance. This outcome contrasts sharply with the original OBJ-V experiment, where univariate in situ storage calibration yielded good performance across all routines, including mHM.

While STARFIT emerges as the highest-performing routine in the OBJ-V-GWW experiment, this comparison is inherently uneven. Because the fitting of this routine required a combination of satellite-estimated storage and observed in situ outflow, its setup is actually a bivariate calibration.

These findings must be interpreted with caution. The evaluation is constrained by a small sample size (n=12) that features reservoirs where the default parameters perform poorly.

Figure A1 Performance of the reservoir storage estimates provided by the Global Water Watch. The swarm plot on the left shows the overall distribution, where each point represents a reservoir and the box plot the quartiles. The map on the right shows the geographical distribution of performance; the dot size indicates storage capacity.

Figure A2 Comparison of model performance across different reservoir routines. The panels represent two experiments: (a) default reservoir parameters, and (b) univariate calibration using reservoir storage estimated by the Global Water Watch (OBJ-V-GWW). Within each panel, three groups of box plots show the performance in terms of outflow, storage, and both variables combined. Four reservoir routines are compared – LISFLOOD, CaMa-Flood, mHM and STARFIT – indicated by colors. The dots represent individual reservoirs, indicating the distribution of overall performance.

Appendix B: Analysis of the KGE^′ components

Figure B1 shows the performance of the multiple simulation runs disaggregated into the individual KGE^′ components: variability (α), bias (β), and correlation (ρ).

The univariate calibration of outflow (OBJ-Q) significantly degrades both the variability and bias of simulated storage. Conversely, the univariate calibration of storage (OBJ-V) only marginally affects the correlation of simulated outflow and has a negligible impact on outflow variability and bias.

The superior performance of the mHM routine in the bivariate calibration (OBJ-QV) is driven by its more accurate representation of outflow variability and correlation, as well as storage variability and bias.

Figure B1 Comparison of the model performance in terms of variability (α), bias (β) and correlation (ρ). The columns correspond to the four experiments: default parameters, calibration of outflow (OBJ-Q), calibration of storage (OBJ-V), and bivariate calibration of outflow and storage (OBJ-QV). Within each panel, two groups of box plots show the performance in terms of outflow and storage. Four reservoir routines are compared – LISFLOOD, CaMa-Flood, mHM and STARFIT – indicated by colors. The dots represent individual reservoirs, indicating the distribution of overall performance.

Appendix C: Performance by reservoir main use

The reservoir routines evaluated in this study are applied uniformly, regardless of the primary purpose of the reservoir. Table C1 analyses how the performance of these routines varies depending on the main reservoir use (as classified in the GRanD database). We exclusively present the values for the bivariate calibration, as it has been shown to yield the highest overall performance and provides the only equitable benchmark for the STARFIT model.

The mHM routine exhibits the highest performance across most categories. This superiority is particularly evident in the outflow simulations, where it is only outperformed by LISFLOOD in the Navigation category. Because mHM is a routine fundamentally driven by water demand time series, it follows logically that it captures release dynamics more effectively than routines primarily governed by available storage. Despite this structural orientation, mHM still emerges as the top-performing routine for storage simulations in three out of the six reservoir use categories. Meanwhile, CaMa-Flood reproduces storage variations slightly better in Flood Control and Water Supply reservoirs. Finally, we note that the lower apparent performance in the Navigation and Recreation categories is likely an artifact of their very small sample sizes.

Table C1 Model performance (in terms of median KGE^′) in the bivariate calibration disaggregated by main reservoir use and reservoir routine. The table shows performance for both target variables: storage and outflow. Bold values indicate the best-performing routine for each reservoir use and target variable, and the Count column indicates the number of reservoirs available in the dataset for that category.