Using GRanD Database and Surface Water Data to Constrain Area–Storage Curve of Reservoirs

: Basic information on global reservoirs is well documented in databases such as GRanD (Global Reservoir and Dam) and ICOLD (International Commission on Large Dams). However, though playing a critical role in estimating reservoir storage variations from remote sensing or hydrological models, area–storage curves of reservoirs are not conveniently obtained nor publicly shared. In this paper, we combine the GRanD database and Landsat-based global surface water extent (GSW) data to derive area–storage curves of reservoirs. The reported storage capacity in the GRanD database and water surface area from GSW data were used to constrain the area–storage curve. The proposed method has the potential to derive area–storage curves of reservoirs larger than 1 km 2 archived in the GRanD database. The derived curves are validated with in situ reservoir data collected in US and China, and the results show that in situ records are well captured by the derived curves both in large and small reservoirs with various shapes. The derived area–storage curves could be employed to advance global monitoring or modeling of reservoir storage dynamics.


Introduction
As one of the most important human water resource management activities, artificial reservoir operation, greatly affects water and energy balance both at the global and local scale [1][2][3]. So far, macro-scale simulation-based reservoir operation algorithms [2,[4][5][6][7][8][9][10] have been developed to simulate the storage and release of global reservoirs, but there are still huge gaps in representing the realistic operating rules and reproducing historic reservoir storage dynamics. Due to the unavailability of observation data or characteristic data of most reservoirs, these algorithms are basically calibrated in several selected reservoirs, and then the calibrated uniform modeling parameters are used in all regions, ignoring the heterogeneity of the operating rules. Observed reservoir operation data, i.e., water level, area or storage, could considerably improve the efficiency of these algorithms through calibrating the optimal modeling parameters [11]. However, the in situ observation data is not publicly accessible in most countries and regions mainly due to security concerns, especially in less developed countries. Therefore, remote sensing products, which usually have global coverage, have offered an important source to provide reservoirs storage estimates [12]. Section 2.2. The derived curves are validated with in situ storage records and Landsat-based area estimates or in situ elevation records in the selected reservoirs in the US and China.

Input Datasets
Three input datasets, i.e., GRanD, GSW, and ETOPO1, are used in this study to derive A-S curves. The basic description and their roles in this study are introduced as the following.
The Global Reservoir and Dam (GRanD) database is the most comprehensive georeferenced database which describes characteristics and geospatial distribution of global reservoirs. It is freely available to the public and scientific community. The original spatial coverage polygon and storage capacity of 5542 reservoirs larger than 1 km 2 are extracted from GRanD for use with modification.
The Global Surface Water (GSW) dataset, which classifies each pixel as water or non-water pixel at 30-m and monthly resolution, is generated by Pekel et al. 2016 [23] using three million more scenes from Landsat 5, 7, and 8 between March 1984 and October 2015, though suffering from cloud covers in some regions. Original GRanD polygon of each reservoir, which usually underestimates the maximum coverage of a reservoir, is extended by GSW to be the maximum coverage polygon used later in this study. The annual maximal and minimal water areas of the reservoirs during 1984-2015 are extracted from GSW to validate the A-S curves, since there is no in situ observation of reservoir water area.
One arc-minute (about 1.8 km) global relief model ETOPO1 is intrinsically a digital elevation model (DEM) showing complete global land topography, i.e., the elevation of underwater bedrock or ice surface (around Antarctic and Greenland). It is produced by combining topographic DEM products SRTM30 (1 km) and GLOBE (1 km), shoreline, and bathymetry data sets, if they are available, from various agencies. Though higher-resolution DEM products, such as 250-m GMTED2010, 30-m SRTM and ASTER GDEM, are also globally available and more precise at vertical accuracy, most of the pixels over a reservoir or lake share the same elevation value, which is the detected elevation of water surface instead of underwater bedrock elevation and thus cannot reveal the topography of the reservoir. Recently, Yamazaki et al. 2017 [26] removed the height error of the existing spaceborne DEMs and developed a high-accuracy map of terrain elevations called Multi-Error-Removed Improved-Terrain DEM (MERIT DEM), but it is still not capable of revealing underwater topography and thus not chosen in this study. After examining all the above DEMs, ETOPO1 is found to show the best performance on the issue we focus on, i.e., area-storage curves derivation. Therefore, ETOPO1 is selected as the most proper DEM product, as far as we know, though the absolute vertical error is about as high as 10 m. This disadvantage will be largely compensated by GRanD and GSW with the calibration process, which is introduced later in Section 2.2.
In situ reservoir storage and forebay water elevation records in the US are collected from the Bureau of Reclamation [27] and are used along with remote sensed reservoir water area series [23] to validate the derived curves. Observed area-storage and depth-storage curves of several reservoirs in China are also collected [28,29].

Methods to Derive Area-Storage Curves
The core idea of the method in this study is to sequentially: 1) identify reservoir coverage polygons by extending the original GRanD polygons with the GSW dataset, 2) construct the hypsographic curve of a reservoir from ETOPO1, 3) calculate the cumulative storages with the water area, increasing step by step to obtain the original area-storage curve, and 4) calibrate the original curve with reported storage capacity in the GRanD database. The four steps are explained in detail below.
First, the original coverage polygon of each reservoir is obtained from GRanD and extended to where water was once detected by GSW (30 m) during 1984-2015. The locations of the ten selected reservoirs in this study are shown in Figure 1. This procedure is necessary because most polygons from GRanD are delineated from snapshot imagery in the SRTM Water Body Database (SWBD) [30], which could underestimate the extent of the reservoir and the area at capacity by as much as 65% ( Figure 2 and Table 1). The updated polygons instead of the original GRanD polygons are used in the study. The areas at capacity, i.e., the maximal coverage areas of reservoirs, is determined by the area of the updated polygons. Water 2020, 12, x FOR PEER REVIEW 4 of 14 study. The areas at capacity, i.e., the maximal coverage areas of reservoirs, is determined by the area of the updated polygons.   Table 1 Table 1 Table 1).  Second, the hypsographic relationship between reservoir surface water areas A(D) and the corresponding relative water depths D is derived from ETOPO1 for each reservoir. The elevation of each pixel in the reservoir polygon is extracted from ETOPO1, resampled by the nearest neighbor resampling method at 30-m resolution, which is consistent with the resolution of the updated coverage polygons. Figure 3 shows ETOPO1 DEM values with different pixel colors over the ten selected reservoirs. Large reservoirs have denser ETOPO1 values than small reservoirs, which implies the derived curves of large reservoirs are more reliable while those of small reservoirs may suffer from larger uncertainties. The relative water depth D of one selected pixel, which represents the water depth in this pixel when all and only the pixels no higher than the selected pixel are fully filled with water, is the elevation of the pixel subtracted by the minimum elevation of the reservoir, and is calculated as: where D is the relative water depth (m), E is the elevation value (m) from ETOPO1, and E min is the minimum of all the elevation values in one reservoir. The corresponding water area A(D) of the specific water depth is the sum area of all the pixels no higher than the selected depth and is calculated as: where A(D) is the water area corresponding to water depth D, A cap is the water area at capacity, N is the number of all 30-m pixels in the reservoir, and N d≤D is the number of all pixels with water depth no higher than D. Based on Equation (2), the hypsographic curve of a reservoir can be plotted as the black solid line in Figure 4. We assume the largest water depth over the reservoir exists near the dam, which is observed and recorded as forebay water elevation in reservoirs management.
Water 2020, 12, x FOR PEER REVIEW 6 of 14 Third, the cumulative storage is calculated with increasing water area by every one percent of reservoir coverage area and then, the area-storage curve is obtained. The A-S curves derived in this step are referred as uncalibrated curves since they are directly derived from ETOPO1 and not calibrated or constrained by GRanD database. The increased storage ∆ , as the light green = % * is calculated as: This would be the uncalibrated area-storage curves which suffer great errors mainly originated from the limited absolute vertical accuracy of ETOPO1. to +1 and water depth increasing from to +1 . is the cumulative reservoir storage with water area as . and are respectively the area at capacity and the depth at capacity. Third, the cumulative storage is calculated with increasing water area by every one percent of reservoir coverage area and then, the area-storage curve is obtained. The A-S curves derived in this step are referred as uncalibrated curves since they are directly derived from ETOPO1 and not calibrated or constrained by GRanD database. The increased storage ∆S, as the light green component (approximate right trapezoid) shown in Figure 4, with water area increasing from A i = i% * A cap to A i+1 = (i + 1)% * A cap is calculated as: where I = 0, 1, 2, . . . , 99. One hundred intervals are found adequate in this study. The cumulative storage S(A i ), as the cyan component shown in Figure 4, with water area as A i = i% * A cap is calculated as: This would be the uncalibrated area-storage curves which suffer great errors mainly originated from the limited absolute vertical accuracy of ETOPO1.
Finally, the area-storage curve is linearly calibrated by the ratio of reported capacity CAP in GRanD over calculated capacity S A cap , following: Relative water depth D(A) is defined as the largest water depth over the reservoir when water area is A. The reservoir storage could not be directly inferred through D-S curve since no such depth data are conveniently accessible in both remote sensing products and in situ observation records. Therefore, we define a standardized depth-storage curve (SD-S curve) that depicts the relationship between water depth and storage. The storage records can be inferred from SD-S curves with forebay water elevation, which is assumed as the deepest water over the reservoir. The standardized water depth SD is literally calculated as: where D min and D max are the minimum and maximum of relative water depth values, and SD ranges from 0 to 1. Combining Equations (5) and (6), the SD-S curve can be obtained as: The SD-S curves are used to make indirect validation with in situ reservoir storage and forebay elevation records.

Results
To evaluate the performance of the proposed method, the derived area-storage (A-S) curves and standardized depth-storage (SD-S) curves of ten reservoirs (reservoirs a-f in the US and g-j in China) are selected to make the validation, where the in situ observation records are available. See characteristics of these 10 reservoirs in Figure 2 and Table 1. The areas at capacity of these selected reservoirs vary from 6 to 82 km 2 , and the capacity vary from 51 to 2187 million m 3 , which is larger than the capacity of 93% (6373 of 6862) reservoirs in the GRanD database. This ensures that the validation covers both the small and large reservoirs with various shapes. The coefficient of determination R 2 is used to evaluate the performance. Figure 5 shows the validation of A-S curves of the ten selected reservoirs. The red points in Figure 5a-f represent the observation records of annual maximal or minimal in situ water storage and remote sensed water area during 1984-2015 in the US, while in Figure 5g-j, the red points represent the observed curves of reservoirs collected in China. It can be interpreted from Figure 5 that the derived curves (black solid lines) capture the distribution pattern of red observation points in all 10 reservoirs, with R 2 ranging from 0.78 to 0.98. Reference A-S curves from van Bemmelen et al. 2016 [24] (black dashed lines) serve as comparison to show the performance of derived A-S curves in this study. The reference R 2 (numbers in brackets in Figure 5) between reference curves and observation records are smaller than our results in Figure 5a-f. This implicates the proposed method in this study performs better in most cases, especially in relatively large reservoirs Figure 5e,f, while the reference curves fail to represent the realistic area-storage relationships probably due to smaller geomorphological similarity around the large reservoirs. The reference curves are derived based on the geomorphological similarity according to van Bemmelen et al. 2016 [24]. Therefore, if a large similarity exits between the existing reservoir and the nearby virtual reservoirs, the A-S relationships of the existing reservoir can be well represented by the reference curves. For a large reservoir, it may be difficult to find nearby virtual reservoirs with large similarity, which may lead to the failed representations of the reference curves. Moreover, the selection of the virtual reservoirs is subjective to expert knowledge and can induce large uncertainties. These two reasons contribute to the large uncertainties of the reference curves in some reservoirs, e.g., in Figure 5e,f. The performances of the two methods are comparable in Figure 5g-j. The derived A-S curves of small reservoirs in Figure 5a-c also show good performance, not limited by the requirement suggested in Zhang and Gao (2016) [16]. Figure 6 shows the validation of SD-S curves, which serves as the indirect validation of the derived A-S curves. The red points (seem like continuous lines due to large number of points) in the figure represent the daily in situ storage and standardized forebay water elevation records with a few outliers eliminated during 1984-2015 in Figure 6a-f, while in Figure 6g-j the red points represent the observed curves of reservoirs collected in China. We assume as in Section 2 that the forebay region adjacent to the dam has the deepest water over the reservoir. The black SD-S curves capture the distribution pattern of red observation points in nearly all 10 reservoirs, with R 2 ranging from 0.56 to 0.97. The SD-S curve and the observation diverge greatly from each other with standardized depth around 0.5-0.7 and converge with each other otherwise. The SD-S curves perform well when the standardized depth (SD) approaches 0 or 1 because the curves are well constrained by the minimal and maximal storage capacities provided by the GRanD database. However, when SD is in the range 0.5-0.7, the curves are less constrained and so a large discrepancy occurs. The scale of the discrepancy varies across different reservoirs depending on the uncertainties of input data and the different shapes of reservoirs. Larger biases are observed in SD-S curves ( Figure 6) than in A-S curves ( Figure 5). The accuracy of the derived SD-S curves may be improved if the water elevation information over the reservoirs is accessible. c also show good performance, not limited by the requirement suggested in Zhang and Gao (2016) [16].

Conclusions and Discussion
In this paper, we demonstrate an approach to derive area-storage curves of reservoirs using readily available datasets, i.e., the GRanD database and Landsat-based global surface water extent (GSW) dataset. It should be noted that this approach relies heavily on reservoir storage data archived in the GRanD database and surface water area from the Landsat-based global surface water extent (GSW) data. As no bathymetry or geometric information has been used, it can not only be used at the reservoirs archived in the GRanD database. Despite its limitations, one advantage is that the storage data has highest accuracy because it is from ground observation rather than remote sensing estimates. The validation results suggest that the derived curves can capture the distribution pattern of in situ observation records both in large and small reservoirs in the US and China. The approach shows the potential to be applicable to 80% of global reservoirs (with area at capacity larger than 1 km 2 ; totally occupying 98% capacity of global reservoirs in GRanD) including small ones. This would advance global-scale monitoring and modeling of reservoir storage and area variations. A previous study [31] developed a framework for monitoring hydrological droughts using a global, long-term, monthly, remotely sensed reservoir surface area dataset. The method proposed in this study can be applied to derive the area-storage curve of reservoirs archived in the GRanD database. In addition to reservoir surface area monitoring, it will facilitate global reservoir water storage monitoring using the area-storage curves and remotely sensed reservoir surface area. Furthermore, a new indicator based on reservoir water storage rather than reservoir surface area can be defined for global hydrological drought and water stress monitoring. For global scale studies, the accuracy of the derived area-storage curves should be acceptable, given that most of those curves are not publicly shared. The curves can also be employed in other related fields which use area-storage relationship of reservoirs as an important input, e.g., estimating evaporation over reservoirs [32].
Uncertainties and bias in A-S curves come from three main sources: 1.
The coverage polygon of each reservoir. The polygon is based on the original GRanD polygon and expanded by the GSW dataset. The uncertainties in both datasets would propagate to affect the real coverage of the reservoir and in result, affect the area at capacity and the topography delineation from ETOPO1.

2.
The coarse resolution and vertical accuracy of elevation data. The coarse resolution and limited vertical accuracy of the ETOPO1 could cause great uncertainties even errors and restrict its applicability. A finer-resolution (10-60 m) and higher-accuracy (10-25 cm) DEM with water depth information may be derived after SWOT (Surface Water and Ocean Topography mission) [33] to be launched in the year 2021. SWOT would be able to observe lakes with repeated high-resolution elevation measurements, and when combined with high-resolution surface water data, would enable area-storage curves of global lake/reservoirs with better accuracy in the future. For now, however, the proposed approach in this study, which utilizes the readily available ETOPO1, can serve as an alternative approach with satisfied accuracy to provide area-storage curves for the hydrology community. Furthermore, the disadvantages of ETOPO1 are largely compensated in the calibration process by the GRanD and GSW dataset. Although other high-resolution DEMs may have the potential to derive the area-storage curves with some new modified approaches, it is beyond the scope of this study and may be considered and examined in the future. 3.
The reported capacity CAP in the GRanD database and the sedimentation effect. The calibration in this study heavily depends on CAP, which is collected from multiple unguaranteed sources with various qualities. Moreover, in this study (as well as in previously mentioned studies), the storage capacity of the reservoir is assumed not to change through time and thus, a stationary capacity is adopted, which means the sedimentation effect [34] is not considered and may introduce bias in reservoirs with severe sediment deposition. Basson (2009) [35] reported that the annual reservoir sedimentation rate is about 0.5%-1% of global reservoir storage volume. It means the adding loss of storage is approximately 50 km 3 per year worldwide [36]. The highest average reservoir sedimentation rate occurs in arid regions such as the Middle East, Australia, Africa, and western United States [37,38]. Wang and Hu (2009) [39] reported that sedimentation had reduced the reservoir capacity in China by 66%. The previous study also showed that 80% of the useful storage capacity for hydropower production will be lost by 2035, and whereas 70% of the storage volumes used for irrigation will be lost due to sedimentation by 2025 in Asia [40].
As the sediment filling effect is increasingly widespread, many measures have been carried out to preserve reservoir capacity [41,42]. The variation of reservoir storage capacity is difficult to monitor and predict. Wisser et al. (2013) [43] estimated the loss in storage capacity for global large reservoirs from 1901 to 2010, and found that the net reservoir capacity is declining as a result of sedimentation (∼5% compared to the installed capacity). It should be noted that using the fixed reservoir storage capacity in the GRanD database would introduce uncertainty to the area-storage curve estimates. There is an urgent need to update the storage capacity in the database. We also note that the annual change of reservoir capacity is usually small compared to the install capacity. If the storage capacity had been recently updated, the fixed storage capacity may be used to derive the area-storage curve for reservoir water storage monitoring. Local information, e.g., more accurate dynamic capacity of reservoirs, could be used to calibrate the A-S curves on a local scale if higher accuracy is desired. However, for global-scale studies, the proposed approach should work well with acceptable bias.