Towards the use of satellite remote sensing to validate reservoir storage in global hydrological models: methodology and pilot study in the CONUS

Reservoir parameters such as dam name, location (longitude and latitude), storage capacity (S_c ), and maximum surface area (A_c ) are provided in the ISIMIP3a protocol (Frieler et al 2023). These data are primarily obtained from the Global Reservoir and Dam Database (GRanD) v1.3, which was developed by the Global Water System Project (Lehner et al 2011). The data also include a set of dams provided by Dr Jida Wang from Kansas State University. This collaboration has resulted in a comprehensive database of 7330 dams, either constructed or under construction, spanning the years 286 to 2020. The cumulative global storage capacity of the database is approximately 7000 km³.

As shown in the Introduction, monitoring reservoirs from space is a rapidly evolving study domain. Many studies have focused on one or a few reservoirs and applied the latest techniques for accurate, continuous, and frequent monitoring (e.g. Das et al 2022). Because our final goal was to validate GHMs globally at a decadal scale, we acquired data covering the globe for over 10 years. To the best of our knowledge, there are six datasets that meet this requirement, namely, Database for Hydrological Time Series of Inland Waters (DAHITI) (Schwatke et al 2015), Global Reservoir Surface Area Dataset (GRSAD) (Zhao and Gao 2018), Hydroweb (Crétaux et al 2011), G-REALM (Birkett and Beckley 2010), GloLakes (Hou et al 2022), and Biswas and Hossain (2022). We excluded Hydroweb and G-REALM because the number of reservoirs monitored is one or two orders of magnitude smaller than other datasets. GloLakes is promising but not used because it does not include many dams in CONUS (i.e. the name of only one in seven dams mentioned later was found in the database). We selected DAHITI and GRSAD because their data were easily accessible. It is important to note that in this article, DAHITI represents storage estimates based on satellite altimetry water level data of reservoirs/lakes, and GRSAD represents storage estimates based on optical remote sensing of surface water extents. It should be noted that new datasets are being produced by applying the latest satellite sensors. Among the other potential datasets, the Cooley et al (2021) dataset, which utilizes the ICESat-2 altimetry, is very promising, but we excluded it because it covers only 2018 and later.

2.1.2.1. DAHITI (surface water level time series)

2.1.2.2. GRSAD (surface area time series)

2.1.2.3. Global Reservoir Bathymetry Dataset (GRBD) (bathymetry)

For most reservoirs, as depicted in figure 1, there is a linear relationship between the water surface area (A) and the water surface elevation (h), represented as:

$$\begin{align}{\text{ }}h{\text{ = }}a \times A{\text{ + }}b{\text{,}}\end{align} \tag{ 1 }$$

Zoom In Zoom Out Reset image size

where a and b represent the slope and intercept, respectively, obtained from a linear regression (Gao et al 2012, Busker et al 2019). These parameters are supplied by the GRBD dataset in our study (refer to table S1).

The volume change for a particular period is then calculated as the area of a trapezoid, as described by Gao et al (2012):

$$\begin{align}\Delta S = \frac{{\left( {{h_2} - {h_1}} \right)\cdot\left( {{A_1} + {A_2}} \right)}}{2}\end{align} \tag{ 2 }$$

where ΔS represents the volume change, A₁ and A₂ are surface areas at the start and end of the period, and h₁ and h₂ are their respective water surface elevations.

Gao's method (Gao et al 2012) extended the linear relationship to reservoir storage capacity (S_c ) of the reservoir, resulting in the expression of the corresponding maximum surface area (A_c ) and maximum water surface elevation (h_c ) as:

$$\begin{align}{S_i} = {S_c} - {\text{ }}\frac{{\left( {{h_c} - {h_i}} \right)\cdot\left( {{A_c} + {\text{ }}{A_i}} \right)}}{2},\end{align} \tag{ 3 }$$

where S_i represents the volume of water stored in the reservoir, corresponding to water surface area A_i and water surface elevation h_i . However, h_c records are, in most cases, missing in the reservoir specification inventories. Therefore, the storage estimation from GRSAD data is computed using the linear equation described in equation (1), as follows:

$$\begin{align}{ }{h_c}{\text{ = }}a \times {A_c}{\text{ + }}b.\end{align} \tag{ 4 }$$

Busker's Method (Busker et al 2019) extends the linear relationship toward the minimum storage (i.e. zero storage); thus, the corresponding surface area is also zero, and the water surface elevation is the bed elevation of the reservoir,

$$\begin{align}{S_i} = \frac{{\left( {{h_i} - b} \right)\cdot{A_i}}}{2}.\end{align} \tag{ 5 }$$

Busker's method requires fewer parameters; A_c and h_c are unnecessary. Additionally, the storage volume can be computed using only h or A by substituting the linear A–h relationship (equation (1):

$$\begin{align}Si = \frac{{{{\left( {{h_i} - b} \right)}^2}}}{{2a}} = \frac{{a\cdot{{\left( {{A_i}} \right)}^2}}}{2}{ }{\text{}}.\end{align} \tag{ 6 }$$

Equation (6) is applicable for GRSAD and DAHITI datasets, which contain time series of surface area and elevation, respectively, within the context of our study. For values of 'a' and 'b', we utilized the GRBD database.

Li et al (2023) recently created the global reservoir storage (GRS) dataset which provides a monthly time series of water storage data for 7,245 reservoirs from 1999 to 2018. They convert the lake surface area time series of GRSAD into storage volume by the area-volume relationship derived from two bathymetry datasets. One is the GRBD for 347 reservoirs; the other is a new satellite-based estimation for the remaining 6,898. The essential difference from the storage derived from equation (6) and GRSAD is whether the relationship in bathymetry is linear or non-linear. In addition to the storage derived from GRSAD and DAHITI, we also used GRS in this study.

2.1.4.1. Normalization

The monthly storage time series can be normalized using the following equation:

$$\begin{align}{S_{{\text{norm,}}i}} = \frac{{{S_i} - {\text{min}}\left( S \right)}}{{{\text{max}}\left( {\text{S}} \right) - {\text{min}}\left( {\text{S}} \right)}},\end{align} \tag{ 7 }$$

where min(s) and max(s) represent the minimum and maximum values of the available monthly storage time series, respectively. By normalizing the monthly storage time series, information about the absolute value of reservoir storage is omitted; only the rate of change information is retained, enabling qualitative validation instead of quantitative.

2.1.4.2. Decomposition

The monthly storage time series (S_y,m ) can be decomposed (figure 2) into annual average storage (S_y ), mean annual seasonal variability (hereafter referred to as seasonal variability or ${\bar S_m}$), and residuals (e_y,m ), as follows:

$$\begin{align}{S_{y,m}} = {S_y} + {\text{ }}{\bar S_m} + {e_{y,m}},\end{align} \tag{ 8 }$$

Zoom In Zoom Out Reset image size

S_y denotes a reservoir's annual average storage volume from January to December, computed by averaging the 12 monthly storage values. ${\bar S_m}$ is determined by calculating the mean storage value for each month after subtracting the mean annual storage for that specific year; thus, it represents storage fluctuation within a year due to seasonal factors. e_y,m constitutes the residual storage value after removing both the annual average storage and the seasonal variability, thus representing the storage component not attributable to annual or seasonal variations. Storage decomposition enables splitting the performance of GHMs for representation into annual storage and seasonal variability, which is helpful for reservoirs with carry-over storage capacity (i.e. storage capacity exceeds the mean annual inflow).

2.1.4.3. Evaluation metrics

Two metrics, Pearson's correlation coefficient (r) and Nash–Sutcliffe efficiency (NSE) (see supplementary text S1), were used to validate time series data. Any months corresponding to missing values in either observation or simulation were excluded from the validation process (the percentage of missing values will be reported later).

The third-round framework of the ISIMIP3a is focused on the evaluation and enhancement of impact models within the context of climate change (Frieler et al 2023). As of 9 June 2023, nine models have participated in the global water sector, but only two have completed simulations that include reservoir outputs. In this study, we utilized two GHMs, namely H08 and WaterGAP2 (WGP). H08 and WGP are the GHMs first incorporated reservoirs and participated in all phases of ISIMIP to date. It has been reported that H08 and WGP are almost at the center of the variation of various models in global water balance calculations (Haddeland et al 2011). We used three meteorological forcings that are bias-adjusted combinations of two reanalyzes, beginning in 1979 for ERA5 (European Centre for Medium-Range Weather Forecasts Reanalysis version 5) and W5 × 10⁵ (WFDE5 over land merged with ERA5 over the ocean; the WFDE5 dataset was generated using the WATCH Forcing Data methodology applied to the surface meteorological variables from the ERA5), respectively (Lange et al 2022):

• 1.  
  Global Soil Wetness Project Phase 3 (GSWP3) combined with W5E5 (GSWP3+W5E5, hereafter referred to as GW)
• 2.  
  20th Century Reanalysis version 3 (20CRv3) combined with W5E5 (20CRv3+W5E5, hereafter referred to as CW)
• 3.  
  20CRv3 combined with ERA5 (20CRv3+ERA5, hereafter referred to as CE)

These forcing data are globally available at 0.5° × 0.5° spatial resolution at daily intervals from 1901 to 2019. Combining the two models and three forcings yields six model simulations, with model and forcing names combined (e.g. H08 forced by GW results in H08_GW). Additional details regarding the simulation protocol can be found at https://protocol.isimip.org/ and in the work by Frieler et al (2023).

2.2.1.1. H08 model

2.2.1.2. WaterGAP model

2.2.1.3. Reservoir operation in H08 and WGP

Reservoir storage is always precisely monitored by dam operators, but the long-term time series are seldom published openly. This has been the primary obstacle in global reservoir modeling and analysis in the past. ResOpsUS (Steyaert et al 2022) is an exhaustive dataset containing historical information about reservoir inflows, outflows, and storage time series for 679 major reservoirs across the United States. The data, with daily temporal resolution, enable detailed analysis of reservoir dynamics. However, the temporal coverage varies among reservoirs based on factors such as construction date and data availability. The dataset spans 1930–2020, with the most robust data from 1980 to 2020. Notably, reservoirs in the dataset contain more than half of the total storage capacity of large reservoirs in the U.S., with a minimum storage threshold of 0.1 km³.

The process of identifying common reservoirs across H08, WGP, GRSAD, GRBD, and ResOpsUS datasets is streamlined by the shared use of the GRanD ID. This sharing facilitates the integration and comparison of data across the different datasets. However, DAHITI uses a unique identification system, thereby requiring individual examination of each reservoir for data availability. Accordingly, a meticulous selection procedure was conducted. First, common reservoirs among H08, WGP, GRSAD, GRBD, and ResOpsUS were identified, and there were 22 in total. This considerable shrink in the number of reservoirs is shown in the Venn diagram in figure 3. The shrink is primarily attributed to the availability of the bathymetry data of GRBD (i.e. data are available for 347 reservoirs globally) and the ground observation data of ResOpsUS (i.e. data are available for 679 reservoirs in the CONUS only). Then, they were searched on the DAHITI website. After a comprehensive review of data availability, only seven reservoirs listed in table 1 were found in all datasets and were thus selected as the foundation for analysis. The locations of these reservoirs in the H08 and WGP models within the 0.5° × 0.5° grids are indicated in table S2. Table S1 displays the storage capacity (S_c ) utilized in our study for these reservoirs. Identifying common reservoirs for all datasets is a prerequisite for this study, which evaluates the agreement between satellite products and the performance of GHMs. The starting and ending dates of ground observations for seven reservoirs in the ResOpsUs dataset are shown in table S3.

Zoom In Zoom Out Reset image size

Initially, satellite data were compared with ground observations to determine compatibility with evaluations of model simulations. Subsequently, simulated reservoir storage from the two ISIMIP3a models, H08 and WaterGAP, was validated against two satellite datasets, GRSAD and DAHITI. Reservoir storage data were examined in raw, normalized, and decomposed. Refer to table 2 for the data utilized.

We analyzed the data at a monthly interval for two reasons. First, we understand that the latest satellite product can estimate the sub-monthly dynamics of reservoir storage. Even though GHM simulations can be run on a daily scale, the ISIMIP3a simulation protocol foresees only monthly resolution of nearly all variables to avoid excessive data storage and data transfer demand. Second, as Masaki et al (2017) clearly show, GHMs have still observed substantial differences in monthly simulations among models. These differences stand out clearly when the catchment area is relatively small or multiple reservoirs are cascading. For a global simulation at a spatial resolution of half a degree, we considered a monthly resolution largely valid.

The monthly time series of storage volume for the seven selected reservoirs (reservoir storage, S) from two remote sensing datasets were compared with ground-based observations (figure 4). The volumes calculated using satellite data (GRSAD and DAHITI) significantly fluctuated depending on the data source and the calculation method (i.e. Busker's method or Gao's method) utilized for raw storage (figures 4(a)–(g)). Intriguingly, even when the same surface area data from GRSAD were used, the storage estimates varied according to the methodologies adopted (i.e. GRSAD_GRBD, GRSAD_ISIMIP, and GRSAD_Busker, represented by gray lines). However, after normalization, the satellite-derived reservoir volumes aligned well with ground observations (figures 4(h)–(n)).

Zoom In Zoom Out Reset image size

Several factors contribute to the variances in satellite-derived S, utilizing different methods and data. For instance, the difference between S derived from GRSAD_ISIMIP and GRSAD_GRBD (see table S1) can be attributed to the varying reservoir storage capacities used. Intriguingly, S calculated using Busker's method (GRSAD_Busker), which does not consider the maximum storage parameters such as A_c, h_c , and S_c , was closest to the observed storage.

The raw storage volume (S) calculated using DAHITI_Busker and GRSAD_Busker displayed considerable agreement for Hoover, Fort Peck, Toledo Bend, and Coolidge (figures 4(a), (c), (d) and (g)). This agreement is promising because these calculations used entirely different satellite products: surface area imagery and water level altimetry. However, discrepancies were evident for Glen Canyon and Structure 193 (figures 4(b) and (e)). For Glen Canyon, temporal storage variability was lost when surface area data from GRSAD were used but not when surface elevation data from DAHITI were used. With the significant differences in surface area parameters between GRSAD and GRBD, the estimation of the linear A–h relationship for Glen Canyon has limitations (Li et al 2021a). The issue stems from the lake polygon included in the GRanD database. Manual correction was applied to GRBD (i.e. slope a and interception b) but not to GRSAD (i.e. monthly lake area A), resulting in inconsistencies.

Structure 193, also known as Lake Okeechobee, had a shallow average depth of 2.7 m but an extensive surface area (1418.77 km²). Therefore, its A–h relationship differed considerably from the schematic shown in figure 1, such that the A–h relationship had a very small value for 'a' and 'b' was negative (equation (1), table S1).

As seen in the panels for the Wesley and Coolidge Dams, the temporal coverage of DAHITI is rather limited (figures 4(f), (g), (m) and (n)). Because the overlap between DAHITI and ground observation is too short, the normalization does not account for the temporal variations over the long ground observation period. Therefore, the NSE score deteriorates on normalization. For such cases, statistics should be viewed with care.

In summary, the raw satellite-based storage time series exhibited considerable uncertainty due to reservoir surface area, the parameters a and b (figure 1), h_c and A_c , and temporal coverage. The success in normalization is mainly due to the proportionate contraction and expansion of the water surface area to different elevations if a significant portion of the total water surface area is monitored. Consequently, a simple normalization enables effective qualitative validation, including signs of change and timing of high/low peaks, of the abilities of hydrological models to simulate reservoir operations.

GRS, one of the latest global satellite storage products, shows remarkably less biased storage estimation, particularly for the Hoover, Glen Canyon, Fort Peck, and Toledo Bend dams (figures 4(a)–(d)). However, there are some systematic discrepancies in the Structure 193, Wesley, and Coolidge dams (figures 4(e)–(g)). After normalization, GRS overlaps with GRSAD (figures 4(h)–(n)). This is not surprising because they commonly used the GRSAD lake area database. Because GRSAD covers a longer temporal duration, hereafter, we analyze only GRSAD.

The normalized time series for satellite-derived and ground-based volumes were decomposed into annual mean storage, seasonal variability, and residuals (figure S1). The correlation coefficient and NSE are displayed in table S6. For reference, the decomposed raw storage (S) is depicted in figure S2.

Overall, satellite-derived decomposed storage components (annual storage, seasonal variability, and residual) consistently compared well with ground-based observation storage components (details in supplementary text S2); correlation (>0.7) and NSE (>0.5) values were high (Moriasi et al 2007). In most cases, annual storage performed prominently among the decomposed components, particularly for GRSAD-based S_norm (table S6).

These satellite-derived components of decomposed normalized monthly storage compared well against their ground observation counterparts and are suitable for validation of model simulations. DAHITI is highly reliable when sufficient, continuous data are available (for instance, data for >5 years). When DAHITI data is unavailable or limited, GRSAD remains a viable (although less robust) alternative. Short-term data (<3 years) and highly discontinuous data, such as Wesley for DAHITI and Coolidge for ground observation, should not be used for validation.

The model simulations were reasonably consistent with satellite-based observations for S_norm,y,m (figures 5(a)–(g)). In particular, simulations for Fort Peck had high correlations with GRSAD and DAHITI. On average, the model performed better when compared with GRSAD than with DAHITI (figures 5(h) and (i)). Particularly, S_norm,y,m for Structure 193 performed well against GRSAD and poorly against DAHITI. Exceptionally, for Toledo Bend (figure 5(d)), performance relative to DAHITI surpassed GRSAD (figures 5(h) and (i)).

Zoom In Zoom Out Reset image size

To illustrate how our proposed method will be finally used, we present a preliminary model and forcing data intercomparison. The performance of WGP ($\bar r$ > 0.5 for 8/12) was superior to the performance of H08 ($\bar r$ > 0.5 for 4/12) (figures 5(h) and (i)). Compared with WGP, H08 was generally more sensitive ($\bar \sigma$) to input forcings (figures 5(a)–(g), S3 and table S7). Among the three forcings, GW ($\bar r$ > 0.5 for 7/12) and CE ($\bar r$ > 0.5 for 7/12) performed better than CW ($\bar r$ > 0.5 for 6/12). A direct comparison of r between CE and GW showed that CE had higher values (CE > GW for 7/12). Noteworthy is the decline in simulation performance since 2005 for unclear reasons. Because most of the DAHITI data included this period, performance relative to GRSAD is generally poor. Considering its long-term consistent coverage, GRSAD demonstrates better consistency with ground observations (i.e. figure 5(h) displays better alignment with figures 5(j) than (i)). Thus, in validating S_norm,y,m for ISIMIP3a, GRSAD is a more reliable evaluation data source than DAHITI when considering the seven reservoirs in this study.

The annual storage simulations were consistent with satellite observations in most cases (figure 6). In particular, Fort Peck and Coolidge simulations demonstrated good agreement with both DAHITI and GRSAD (figures 6(c) and (g)). For most reservoirs, the average correlation coefficient was >0.5 for simulations across two models and three forcings compared with GRSAD (figure 6(h)); this finding was consistent with results from ground observation comparisons (figure 6(j)). We only used GRSAD for further comparisons of GHMs and forcings in this subsection because there were apparent discrepancies between DAHITI and ground observations (figures 6(i)–(j)). This is because more than 50% of DAHITI data covers the period after 2005, when the model's performance deteriorated (section 3.2.1). Because of the difference in the validation period, the results of GRSAD and DAHITI were incomparable.

Zoom In Zoom Out Reset image size

From figure 6, the model and forcing data intercomparison goes as follows. The performances of the two GHMs concerning S_norm,y correlations with satellite data were nearly equivalent; WGP ($\bar r$ > 0.5 for 7/7) was slightly superior to H08 ($\bar r{\text{ }}$ > 0.5 for 6/7) (figure 6(h)). Additionally, H08 displayed a slightly larger standard deviation than WGP (figures 6(a)–(g)), indicating that it had substantially more interannual variability with input forcings. Among the input forcings, CE performed best in terms of the correlation coefficient, followed by GW and CW. WGP_CE correlation coefficients were considerably higher than the correlation coefficients of other GHM-Forcing combinations for most satellites (figure S4).

The GHMs readily captured the interannual variation of reservoir storage.

The model simulations adequately captured the seasonal cycle of reservoir storage in many instances (figure 7). The correlations of simulations with both satellite datasets were generally high, such that many values exceeded 0.5 (21/35 for GRSAD and 24/35 for DAHITI), except for the Hoover and Wesley dams. For instance, the simulated peak timing of the Hoover Dam (April) lagged the satellite products (from January to March) by 1–3 months, resulting in weaker correlations.

Zoom In Zoom Out Reset image size

Both H08 ($\bar r$ > 0.5 for 9/14) and WGP ($\bar r$ > 0.5 for 8/14) performed particularly well in terms of simulating monthly variability in most reservoirs (figures 7(h) and (i)); WGP was superior to H08 (for 9/14 cases). H08 demonstrated less robust performance for the Hoover and Wesley Dams, but WGP displayed relatively strong correlations with observations for these reservoirs. Therefore, WGP demonstrated superior overall performance compared with H08. Moreover, H08 exhibited more significant variability according to input forcings than WGP. Among the input forcings (figures 6(h) and (i)), GW simulations ($\bar r$ > 0.5 for 11/14) outperformed CW ($\bar r$ > 0.5 for 9/14) and CE ($\bar r$ > 0.5 for 8/14). Even in terms of r values, GW performed best for 8/14 cases among the three forcings (figures 6(h) and (i)). This result is consistent with the evaluation relative to ground observations, where 5/7 reservoirs had the highest correlation for GW (figure 7(j)).

A comparison of simulations showed that DAHITI and GRSAD aligned with ground observations; DAHITI demonstrated relatively better consistency. We observed two instances of contradictory outcomes. First, WGP validation for the Fort Peck Dam and the Toledo Bend Dam, where a weak correlation with GRSAD differed from a strong correlation with DAHITI. Cases with very short data availability periods, such as Wesley and Coolidge for DAHITI, should be excluded from validation. In these instances, GRSAD should be used because it can appropriately capture the seasonal variation due to its long-term data availability. However, when DAHITI has sufficient temporal coverage, it outperforms GRSAD. Second, the H08-WGP comparison for Structure 193. In this case, the simulated seasonal cycle of WGP for Structure 193 closely correlated with GRSAD, but the amplitude was considerably smaller. These discrepancies lead to questions regarding the reliability of a single evaluation data source. Therefore, in the absence of ground-based observations, multiple satellite data products and metrics should be used to increase confidence in validation and intercomparison results. When a sufficient number of reliable satellite products were available, it would be possible to calculate the mean and ranges of satellite data ensemble.

This study inherited some uncertainties from the previous efforts on satellite monitoring of reservoirs. There were four main issues. First, we assumed a linear A–h relationship for reservoir storage. This relationship is not genuinely linear, particularly when the reservoir is near full or empty. Many existing approaches require knowledge of water surface elevation at storage capacity; such information is not currently available in published global reservoir inventories. Monitoring S_c and h_c from space is challenging because the water level of operational reservoirs has very little opportunity to reach the exact total capacity. The bathymetry approach (e.g. Messager et al 2016) also does not provide information about the water level at total capacity unless other methods specify the highest shoreline. Therefore, significant uncertainties may arise when calculating reservoir storage using these parameters (Gao et al 2012). It should also be noted that A–h relationship is not necessarily constant and subject to change due to exceptional events. Furthermore, the records in the inventories are not necessarily error-free or consistent with information in other inventories. This indicates a need for extensive quality checks among global reservoir inventories. Due to these issues, we abandoned using absolute storage estimates and normalized them. Although we showed that the timing and rate of rise and fall can be validated, a significant limitation of our study was the loss of storage magnitude information. The GRS database (Li et al 2023) adopts a more realistic and non-linear A–h relationship in converting lake area data into storage. As seen in figure 4, GRS outperforms the products with a linear assumption. Improvements in bathymetry estimation are the key to addressing this problem. Second, the limitations of area-based satellite products (i.e. GRSAD). Discrepancies in the consideration of water surface area extents (i.e. distinguishing between reservoir and river), as noted in the case of the Glen Canyon Dam (section 3.1.2), lead to concerns about the reliability of water surface area datasets. Third, there are limitations to altitude-based satellite products (i.e. DAHITI). This study found that DAHITI agrees better with the ground observation than GRSAD, but the advantage is largely deteriorated by its high frequency of missing data. Further technical advancement in data processing is expected for more extensive spatio-temporal coverage of this type of data. Finally, since 2005, there has been a clear discrepancy between simulation results and ground observations. This issue should be further examined from various aspects, including validation with other variables and at different locations.

Model intercomparison aims to identify the superior model among members and the reasons for superiority. The implications earned from our pilot study are summarized as follows.

First, in general, it is challenging to identify the superior model. We learned that the model and forcing performance vary according to the reservoirs. Increasing the number of reservoirs is a prerequisite of model intercomparison to avoid drawing erroneous conclusions by chance. The performance also varies by the spatio-temporal scale. Generally, the discrepancies among forcing data are relatively minor at a continental and monthly scale but exhibit significant variations at a catchment and daily scale. A detailed analysis will be performed in the forthcoming papers with more models and greater spatial coverage.

Second, it is difficult to single out the performance of the reservoir operation scheme because the simulation performance is dependent on the overall hydrological simulation (i.e. the daily water balance calculation carried out by the land surface hydrology and the river routing sub-models in GHMs.) For instance, the superior performance of WGP over H08 could be attributed to the calibration of river discharge (see section 2.2.1), which enhances the overall accuracy of reservoir inflow and outflow estimates in WGP. It could also be attributed to the fact that H08 tends to produce more significant temporal variation in storage than WGP (i.e. the sigma in figures 5, 6 and table S7). The performance metric tends to drop largely for models with the former tendency when the phase was incorrectly simulated. Furthermore, GHM performance cannot be determined by storage volume alone because it is affected by inflows and outflows.

Three additional approaches are considered for comprehensive intercomparison and validation of GHMs. The first approach is a systematic offline intercomparison. The performance of individual reservoir operations algorithms can be compared by applying them on a stand-alone basis at specific reservoirs for which reliable inflow, outflow, and storage data are available. The obvious drawback is that this approach is only applicable to data-rich regions. The second approach is a systematic online intercomparison. By implementing multiple reservoir operation algorithms in each GHM, modelers can identify the impact of the reservoir operation algorithms on their simulations. The disadvantage is that this approach requires intensive coding and more simulation runs. The final approach is a cross-validation. By adding more validation variables, e.g. river discharge upstream and downstream of a reservoir, it is possible to infer whether the simulated inflows and releases are overestimated.