NASA’s MODIS/VIIRS Global Water Reservoir Product Suite from Moderate Resolution Remote Sensing Data

: Global reservoir information can not only beneﬁt local water management but can also improve our understanding of the hydrological cycle. This information includes water area, elevation, and storage; evaporation rate and volume values; and other characteristics. However, operational wall-to-wall reservoir storage and evaporation monitoring information is lacking on a global scale. Here we introduce NASA’s new MODIS/VIIRS Global Water Reservoir product suite based on moderate resolution remote sensing data—the Moderate Resolution Imaging Spectroradiometer (MODIS), and the Visible Infrared Imaging Radiometer Suite (VIIRS). This product consists of 8-day (MxD28C2 and VNP28C2) and monthly (MxD28C3 and VNP28C3) measurements for 164 large reservoirs (MxD stands for the product from both Terra (MOD) or Aqua (MYD) satellites). The 8-day product provides area, elevation, and storage values, which were generated by ﬁrst extracting water areas from surface reﬂectance data and then applying the area estimations to the pre-established Area–Elevation (A–E) relationships. These values were then further aggregated to monthly, with the evaporation rate and volume information added. The evaporation rate and volume values were calculated after the Lake Temperature and Evaporation Model (LTEM) using MODIS/VIIRS land surface temperature product and meteorological data from the Global Land Data Assimilation System (GLDAS). Validation results show that the 250 m area classiﬁcations from MODIS agree well with the high-resolution classiﬁcations from Landsat (R 2 = 0.99). Validation of elevation and storage products for twelve Indian reservoirs show good agreement in terms of R 2 values (0.71–0.96 for elevation, and 0.79–0.96 for storage) and normalized root-mean-square error (NRMSE) values (5.08–19.34% for elevation, and 6.39–18.77% for storage). The evaporation rate results for two reservoirs (Lake Nasser and Lake Mead) agree well with in situ measurements (R 2 values of 0.61 and 0.66, and NRMSE values of 16.25% and 21.76%). Furthermore, preliminary results from the VIIRS reservoir product have shown good consistency with the MODIS based product, conﬁrming the continuity of this 20-year product suite. This new global water reservoir product suite can provide valuable information with regard to water-sources-related studies, applications, management, and hydrological modeling and change analysis such as drought monitoring.


Brief Description of the Selected Reservoirs
This product suite provides data for 164 reservoirs (locations are shown in Figure 1) with areas larger than 25 km 2 , which were selected from the Global Reservoir Bathymetry Dataset (GRBD) [33]. It includes 151 man-made reservoirs (2672 km 3 ) and 13 regulated natural lakes (23,811 km 3 ). The total storage capacity of the 151 man-made reservoirs represents 45.82% of the global capacity (in its category) according to the Global Reservoir and Dam Database (GRanD) [34]. Table 1 summarizes the attributes of these 151 reservoirs at the global and continental scales. More details about the individual reservoirs (both man-made reservoirs and regulated natural lakes) and their attributes are provided in Table A1. At the continental scale, it has the largest coverage in Africa (76.61%) and the lowest in Europe (20.59%) in terms of storage capacity.
Sens. 2021, 13, x FOR PEER REVIEW 3 of 29 Meanwhile, some new approaches have been recently developed to estimate evaporation rates and losses from space. For instance, Zhang et al. [14] estimated the monthly evaporation volumes based on pan-derived evaporative rates and Landsat surface areas for more than 200 reservoirs in Texas. Zhao and Gao [28] used the Penman Equation (with the heat storage and fetch effects addressed), and generated a first longterm evaporation data record for over 700 reservoirs in the contiguous United States. Zhao et al. [29] further improved the calculation of the heat storage change term by leveraging MODIS water surface temperature data. Many other approaches have been developed and tested at individual locations [30][31][32].
Despite the development of remotely sensed reservoir datasets, consistent, comprehensive, long-term, and operationally monitored reservoir products are still lacking at the global scale. Therefore, we present here a newly developed NASA global water reservoir product suite (MxD28 and VNP28) from moderate resolution remote sensing data. The suite provides 8-day and monthly measurements of area, elevation, and storage, along with monthly evaporation rates and evaporative volumetric losses for 164 global reservoirs. This paper is organized as follows. Section 2 describes the overview of the reservoir product. Section 3 explains the materials and methods that were used to generate the MODIS products. Section 4 shows the MODIS product validation results and evaluates the consistency between the MODIS and VIIRS results. Section 5 discusses the limitations and potential applications of this reservoir product.

Brief Description of the Selected Reservoirs
This product suite provides data for 164 reservoirs (locations are shown in Figure 1) with areas larger than 25 km 2 , which were selected from the Global Reservoir Bathymetry Dataset (GRBD) [33]. It includes 151 man-made reservoirs (2672 km 3 ) and 13 regulated natural lakes (23,811 km 3 ). The total storage capacity of the 151 man-made reservoirs represents 45.82% of the global capacity (in its category) according to the Global Reservoir and Dam Database (GRanD) [34]. Table 1 summarizes the attributes of these 151 reservoirs at the global and continental scales. More details about the individual reservoirs (both man-made reservoirs and regulated natural lakes) and their attributes are provided in Table A1. At the continental scale, it has the largest coverage in Africa (76.61%) and the lowest in Europe (20.59%) in terms of storage capacity.    1 The area and storage columns correspond to the values at capacity. "Storage pct." represents the percentage of storage included in our product with respect to the total storage capacity (referring to Global Reservoir and Dam Database (GRanD)) of the corresponding geographical domain.

Reservoir Product Components
The reservoir product suite is organized based on different temporal resolutions: 8-day (MxD28C2 and VNP28C2) and monthly (MxD28C3 and VNP28C3) ( Table 2). The 8-day product includes reservoir area, elevation, and storage values. The monthly product adds reservoir evaporation rate and volume values (in addition to the above-mentioned measurements) for each reservoir.

Input Datasets
The input datasets used for the development of the reservoir product are summarized in Table 3. For MxD28C2, the 8-day Terra/Aqua surface reflectance (MxD09Q1) [35] data were collected for water area extraction. For VNP28C2, we used the 8-day VIIRS/NPP surface reflectance (VNP09H1) [36]. In particular, only the near-infrared (NIR) band was used due to its high spatial resolution (i.e., 250 m for MODIS and 500 m for VIIRS). The NIR band has been commonly utilized for the extraction of water bodies because it is strongly absorbed by water but scarcely absorbed by terrestrial dry soil and vegetation [37]. The Area-Elevation (A-E) relationships were adopted from the GRBD [33], which have proven to be of high quality through validation against in situ data. Then, the 8-day water area estimations were applied to the A-E relationships to derive elevation and storage values (see Section 3.2.1 for more details). Moreover, we used the 8-day day/night land surface temperature (LST) products (MxD21A2 [38] and VNP21A2 [39]) and Global Land Data Assimilation System (GLDAS) [40] meteorological forcing data to estimate the evaporation rates and volumes (see Section 3.2.2 for a more detailed approach).

Methodology
In this section, we use MODIS products to demonstrate the algorithms, and the same algorithms were adopted to generate VIIRS-based products. Detailed descriptions of all of the algorithms involved in the product's development (e.g., the enhancement water classification algorithm and evaporation rate estimation algorithm) are provided in the Algorithm Theoretical Basis Documentation (ATBD) [41], which is available through the product website (https://modis-land.gsfc.nasa.gov/modgwr.html, accessed on 30 December 2020). Since the product suite has yet to be formally released by NASA, the MODIS results presented in this paper were generated by executing the algorithms on the Google Earth Engine (GEE) [42], while the VIIRS results were generated on our local computer cluster (due to lack of VIIRS data on GEE). Figure 2 shows the algorithm for generating the MxD28C2 8-day product. First, the 8-day reservoir area values were extracted from the 250-m Near Infrared (NIR) band, along with the quality assurance (QA) band, of MODIS Terra/Aqua surface reflectance (MxD09Q1) data. During the classification process, an area enhancement algorithm (Zhang et al. [26]) was adopted to minimize the effects of various sources of contamination (e.g., cloud and snow/ice). Then, the area values were applied to the A-E relationship for the given reservoir provided by the GRBD [33] to calculate the corresponding elevation values. Lastly, the reservoir storage was estimated using Equation (1) (Gao et al. [25]):

Algorithm for the 8-Day Products
where V c , A c , and h c represent storage, area, and water elevation values at capacity; and V MODIS , A MODIS , and h MODIS are the estimated storage, area, and water elevation values from MODIS, respectively.

Methodology
In this section, we use MODIS products to demonstrate the algorithms, and the same algorithms were adopted to generate VIIRS-based products. Detailed descriptions of all of the algorithms involved in the product's development (e.g., the enhancement water classification algorithm and evaporation rate estimation algorithm) are provided in the Algorithm Theoretical Basis Documentation (ATBD) [41], which is available through the product website (https://modis-land.gsfc.nasa.gov/modgwr.html, accessed on 30 December 2020). Since the product suite has yet to be formally released by NASA, the MODIS results presented in this paper were generated by executing the algorithms on the Google Earth Engine (GEE) [42], while the VIIRS results were generated on our local computer cluster (due to lack of VIIRS data on GEE). Figure 2 shows the algorithm for generating the MxD28C2 8-day product. First, the 8-day reservoir area values were extracted from the 250-m Near Infrared (NIR) band, along with the quality assurance (QA) band, of MODIS Terra/Aqua surface reflectance (MxD09Q1) data. During the classification process, an area enhancement algorithm (Zhang et al. [26]) was adopted to minimize the effects of various sources of contamination (e.g., cloud and snow/ice). Then, the area values were applied to the A-E relationship for the given reservoir provided by the GRBD [33] to calculate the corresponding elevation values. Lastly, the reservoir storage was estimated using Equation (1) (Gao et al. [25]):

Algorithm for the 8-Day Products
where , , and ℎ represent storage, area, and water elevation values at capacity; and , , and ℎ are the estimated storage, area, and water elevation values from MODIS, respectively. The water area extraction values from MODIS serve as a foundation for elevation and storage estimations. Here we use our experience with the Cahora Bassa reservoir (Mozambique, Africa) as an example to explain the classification and enhancement processes. When classifying this reservoir, we initially buffered the shapefile (provided by HydroLAKES database [43]) outward by 1 km and then further used it to extract the NIR band images. The buffered area covers the maximum water extent, and the following classification and enhancement operations were implemented inside of this area to reduce the computational cost. For each 8-day period, the MxD09Q1 NIR images that overlap The water area extraction values from MODIS serve as a foundation for elevation and storage estimations. Here we use our experience with the Cahora Bassa reservoir (Mozambique, Africa) as an example to explain the classification and enhancement processes. When classifying this reservoir, we initially buffered the shapefile (provided by HydroLAKES database [43]) outward by 1 km and then further used it to extract the NIR band images. The buffered area covers the maximum water extent, and the following classification and enhancement operations were implemented inside of this area to reduce the computational cost. For each 8-day period, the MxD09Q1 NIR images that overlap with the reservoir were Remote Sens. 2021, 13, 565 6 of 29 first selected (Figure 3a). The pixels contaminated by clouds, cloud shadow, and snow/ice (identified using the QA band of MxD09Q1) were labeled as 'NoData' (the blank areas in Figure 3b represent the contaminated pixels). Then, the Otsu thresholding method [44] was used to extract the water area from the NIR image ( Figure 3c). It is obvious that the water area was underestimated due to contaminations (as shown in Figure 3c). Lastly, in combination with the water occurrence (percentile) image provided by the GSW (Figure 3d), we applied the enhancement algorithm developed by Zhang et al. [26] (Figure 3e) to correct the underestimations. The enhanced water areas show clear boundaries (Figure 3f). More details about the framework of MODIS water classification and enhancement can be found in Zhang et al. [26]. Note that it may happen that some extreme area values still exist even after the enhancement process. To eliminate the effects of these outliers, we used three standard deviations to remove them and used an interpolation method to fill in the data gap. This issue will be addressed in a future release of the enhancement algorithm [29]. with the reservoir were first selected ( Figure 3a). The pixels contaminated by clouds, cloud shadow, and snow/ice (identified using the QA band of MxD09Q1) were labeled as 'NoData' (the blank areas in Figure 3b represent the contaminated pixels). Then, the Otsu thresholding method [44] was used to extract the water area from the NIR image ( Figure  3c). It is obvious that the water area was underestimated due to contaminations (as shown in Figure 3c). Lastly, in combination with the water occurrence (percentile) image provided by the GSW (Figure 3d), we applied the enhancement algorithm developed by Zhang et al. [26] (Figure 3e) to correct the underestimations. The enhanced water areas show clear boundaries (Figure 3f). More details about the framework of MODIS water classification and enhancement can be found in Zhang et al. [26]. Note that it may happen that some extreme area values still exist even after the enhancement process. To eliminate the effects of these outliers, we used three standard deviations to remove them and used an interpolation method to fill in the data gap. This issue will be addressed in a future release of the enhancement algorithm [29].

Algorithm for the Monthly Products
The MxD28C3 product includes the evaporation rate and volumetric evaporation loss in addition to the area, elevation, and storage results at monthly temporal resolution. Figure 4 shows the algorithm for generating the MxD28C3 monthly product. The monthly area values were first estimated based on the composite of the 8-day area classifications and then converted into monthly elevation and storage results using the A-E relationship ( Figure 4). In addition, monthly evaporation rates were estimated after the Lake Temperature and Evaporation Model (LTEM) [29] using MODIS LST product (MxD21A2) and meteorological data from GLDAS [40]. Lastly, the monthly evaporative volumetric losses were calculated as the product of evaporation rate and reservoir area values.

Algorithm for the Monthly Products
The MxD28C3 product includes the evaporation rate and volumetric evaporation loss in addition to the area, elevation, and storage results at monthly temporal resolution. Figure 4 shows the algorithm for generating the MxD28C3 monthly product. The monthly area values were first estimated based on the composite of the 8-day area classifications and then converted into monthly elevation and storage results using the A-E relationship ( Figure 4). In addition, monthly evaporation rates were estimated after the Lake Temperature and Evaporation Model (LTEM) [29] using MODIS LST product (MxD21A2) and meteorological data from GLDAS [40]. Lastly, the monthly evaporative volumetric losses were calculated as the product of evaporation rate and reservoir area values. with the reservoir were first selected ( Figure 3a). The pixels contaminated by clouds, cloud shadow, and snow/ice (identified using the QA band of MxD09Q1) were labeled as 'NoData' (the blank areas in Figure 3b represent the contaminated pixels). Then, the Otsu thresholding method [44] was used to extract the water area from the NIR image ( Figure  3c). It is obvious that the water area was underestimated due to contaminations (as shown in Figure 3c). Lastly, in combination with the water occurrence (percentile) image provided by the GSW (Figure 3d), we applied the enhancement algorithm developed by Zhang et al. [26] (Figure 3e) to correct the underestimations. The enhanced water areas show clear boundaries ( Figure 3f). More details about the framework of MODIS water classification and enhancement can be found in Zhang et al. [26]. Note that it may happen that some extreme area values still exist even after the enhancement process. To eliminate the effects of these outliers, we used three standard deviations to remove them and used an interpolation method to fill in the data gap. This issue will be addressed in a future release of the enhancement algorithm [29].

Algorithm for the Monthly Products
The MxD28C3 product includes the evaporation rate and volumetric evaporation loss in addition to the area, elevation, and storage results at monthly temporal resolution. Figure 4 shows the algorithm for generating the MxD28C3 monthly product. The monthly area values were first estimated based on the composite of the 8-day area classifications and then converted into monthly elevation and storage results using the A-E relationship ( Figure 4). In addition, monthly evaporation rates were estimated after the Lake Temperature and Evaporation Model (LTEM) [29] using MODIS LST product (MxD21A2) and meteorological data from GLDAS [40]. Lastly, the monthly evaporative volumetric losses were calculated as the product of evaporation rate and reservoir area values.

This paper primarily focuses on validating the MODIS-based products (Sections 4.1-4.3)
and evaluating the consistencies between VIIRS and MODIS products (Section 4.4). With regard to MODIS, when both 8-day MOD (Terra) and MYD (Aqua) area results are available on the same date, we selected the one with less cloud contamination. Similarly, the validations of the MODIS elevation and storage products also adopted those values associated with the less contaminated areas.

Comparing Water Surface Areas with Landsat-Based Results
At the global scale, long-term in situ reservoir area records are still lacking. Therefore, we compared the MODIS area values with Landsat-based results (at a finer spatial resolution of 30 m) for purposes of area validation. The Landsat monthly reservoir area values for all of the 164 reservoirs were collected from GRSAD between 2000 and 2018 [20]. GRSAD corrected the water area underestimation of the GSW dataset caused by both cloud/shadow/ice contamination and the Landsat-7 scan line corrector failure [20]. Note that the Landsat-based area estimation for a given month was based on the one or two images obtained during that month, while the monthly MODIS area value was derived from the composite of the 8-day classification. However, due to the deficiency of the in situ area values, we used this Landsat-based dataset to validate the overall consistency of the MODIS area products.
According to Figure 5, the MODIS-based area values agree well with the Landsatbased results (with an R 2 value over 0.99). Additionally, the data points are mainly centered on the 1:1 line (slope = 0.99). The disagreements that do exist can be attributed to two sources. The first is because Landsat and MODIS collected data at different dates. If a reservoir experienced a large change within a month, it may have caused a large area discrepancy. The second is related to the low spatial resolution of MODIS, which makes it more susceptible to mixed pixels. This can explain the area underestimations for relatively small reservoirs.

Validating the MODIS Elevation and Storage Products against In Situ Observations
For the elevation and storage validations, we collected in situ daily observations for twelve Indian reservoirs (Ukai, Matatila, Ranapratap Sagar, Gandhi Sagar, Ban Sagar,

Validating the MODIS Elevation and Storage Products against In Situ Observations
For the elevation and storage validations, we collected in situ daily observations for twelve Indian reservoirs (Ukai, Matatila, Ranapratap Sagar, Gandhi Sagar, Ban Sagar, Bargi, Hirakud, Jayakwadi, Sriram Sagar, Nagarjuna Sagar, Yeleru, and Tungabhadra) from the Indian Central Water Commission (http://cwc.gov.in/, accessed on 30 December 2020) between 2000 and 2019. We used the Indian reservoirs for validation purposes because they experience large dynamics that can better evaluate the efficiency of our algorithm.
The validation results of the 8-day MODIS elevation and storage products are shown in Figures 6 and 7, respectively. Overall, the elevation estimations from MODIS agree well with the in situ data ( Next, we validated the monthly elevation and storage products over these twelve Indian reservoirs. The daily in situ elevation and storage values were averaged at a monthly step and were then compared to the monthly MODIS products (Figures A1 and  A2). The validation results show similar patterns as those of the 8-day products but higher

Validating the Evaporation Rate Product against In Situ Observations
The evaporation rate results were validated over two locations-Lake Nasser in Africa and Lake Mead in North America-where high-quality in situ observations are available. The eddy covariance (EC) evaporation rate measurements for Lake Mead between 2010 and 2015 were provided by the United States Geological Survey (USGS) [12]. With regard to Lake Nasser, the evaporation rate estimations were obtained using the Bowen ratio energy budget (BREB) method [45]. Although it is not as accurate as EC observations, the BREB method has been widely utilized to estimate evaporation rates due to its operability and reliability [46]. As shown in Figure 8, the MODIS evaporation rate products have good overall agreement with values obtained via observation (with R 2 values of 0.61 and 0.66, and NRMSE values of 16.25% and 21.76%). Compared to the results that were calculated using the regular Penman equation (without heat storage), our MODIS evaporation rate estimates have shown great improvements in terms of both annual peak values and seasonal variations. For Lake Nasser, the R 2 value increased from 0.30 to 0.61, with the NRMSE decreasing from 26.86% to 16.25%. A greater improvement is observed for Lake Mead in terms of both R 2 (from 0.26 to 0.66) and NRMSE (from 41.06% to 21.76%). These results suggest that better evaporation rate estimates are achieved when considering heat storage within the algorithm. It should be noted that this evaporation rate algorithm was validated over more locations in Zhao et al. [29]. Comparisons of Next, we validated the monthly elevation and storage products over these twelve Indian reservoirs. The daily in situ elevation and storage values were averaged at a monthly step and were then compared to the monthly MODIS products (Figures A1 and A2). The validation results show similar patterns as those of the 8-day products but higher accuracies. This is because the monthly reservoir area values were generated from the composited results of three or four 8-day reservoir areas from MxD28C2, and the composition process further reduced the adverse effects of cloud contamination at the 8-day time step. As shown in Figure A1, the MODIS-based elevations show good consistency with the in situ measured data, with an average R 2 value of 0.90, an average RMSE value of 1.99 m, and an average NRMSE value of 11.33%. Regarding the storage validations ( Figure A2), the results are consistent with those of the elevations, with an average R 2 value of 0.91, an average RMSE value of 0.43 km 3 , and an average NRMSE value of 11.91%.

Validating the Evaporation Rate Product against In Situ Observations
The evaporation rate results were validated over two locations-Lake Nasser in Africa and Lake Mead in North America-where high-quality in situ observations are available. The eddy covariance (EC) evaporation rate measurements for Lake Mead between 2010 and 2015 were provided by the United States Geological Survey (USGS) [12]. With regard to Lake Nasser, the evaporation rate estimations were obtained using the Bowen ratio energy budget (BREB) method [45]. Although it is not as accurate as EC observations, the BREB method has been widely utilized to estimate evaporation rates due to its operability and reliability [46]. As shown in Figure 8, the MODIS evaporation rate products have good overall agreement with values obtained via observation (with R 2 values of 0.61 and 0.66, and NRMSE values of 16.25% and 21.76%). Compared to the results that were calculated using the regular Penman equation (without heat storage), our MODIS evaporation rate estimates have shown great improvements in terms of both annual peak values and seasonal variations. For Lake Nasser, the R 2 value increased from 0.30 to 0.61, with the NRMSE decreasing from 26.86% to 16.25%. A greater improvement is observed for Lake Mead in terms of both R 2 (from 0.26 to 0.66) and NRMSE (from 41.06% to 21.76%). These results suggest that better evaporation rate estimates are achieved when considering heat storage within the algorithm. It should be noted that this evaporation rate algorithm was validated over more locations in Zhao et al. [29]. Comparisons of evaporation rates calculated from different datasets (GLDAS + MxD21 and TerraClimate + MxD11) indicate that the performances are similar, even though the LST and meteorological data are slightly different ( Figure A3).

Consistencies between VIIRS and MODIS Products
MODIS sensors have successfully operated for more than 20 years. The MODIS successor, VIIRS, was launched on 28 October 2011, onboard the Suomi National Polarorbiting Partnership (Suomi NPP). Therefore, VIIRS plays an important role in extending the continuity of this product suite. The reliability of the water area estimation controls the quality of the other products (i.e., the elevation, storage, and evaporation results). Because the VIIRS 8-day reflectance data (VNP09) are not available on GEE and the VIIRS input data are very large and computationally intensive, we selected six Indian reservoirs from which to evaluate the consistency of VIIRS and MODIS data by comparing the 8-day area estimations from January 2012 to August 2020.
As seen in Figure 9, the VIIRS-based area values agree well with those of MODIS, with R 2 values ranging from 0.80 to 0.90. The 8-day elevation and storage values also have good consistency, as they are based on the area estimations. In addition, high-quality monthly products can be generated from 8-day composite results (which are not shown here). Note that the VIIRS results are noisier than those from MODIS because of its coarser spatial resolution (500 m vs. 250 m). We also applied the outlier removal and gap-filling procedures to the area estimations from VIIRS, which led to smoother time series values.

Consistencies between VIIRS and MODIS Products
MODIS sensors have successfully operated for more than 20 years. The MODIS successor, VIIRS, was launched on 28 October 2011, onboard the Suomi National Polarorbiting Partnership (Suomi NPP). Therefore, VIIRS plays an important role in extending the continuity of this product suite. The reliability of the water area estimation controls the quality of the other products (i.e., the elevation, storage, and evaporation results). Because the VIIRS 8-day reflectance data (VNP09) are not available on GEE and the VIIRS input data are very large and computationally intensive, we selected six Indian reservoirs from which to evaluate the consistency of VIIRS and MODIS data by comparing the 8-day area estimations from January 2012 to August 2020.
As seen in Figure 9, the VIIRS-based area values agree well with those of MODIS, with R 2 values ranging from 0.80 to 0.90. The 8-day elevation and storage values also have good consistency, as they are based on the area estimations. In addition, high-quality monthly products can be generated from 8-day composite results (which are not shown here). Note that the VIIRS results are noisier than those from MODIS because of its coarser spatial resolution (500 m vs. 250 m). We also applied the outlier removal and gap-filling procedures to the area estimations from VIIRS, which led to smoother time series values. Additionally, we selected Lake Mead to compare the evaporation rates from January 2012 to April 2015 ( Figure 10). Validation results show that the evaporation rates of VNP28C3 have a higher accuracy than those of MYD28C3 in terms of R 2 (0.76 vs. 0.66) and NRMSE (15.92% vs. 21.76%). Meanwhile, the evaporation rate results suggest a good consistency between VNP28C3 and MYD28C3 (R 2 = 0.92). It is noteworthy that the difference between MODIS and VIIRS evaporation rates is mainly attributed to the LST data. It has been reported that MODIS and VIIRS LST data have good agreement (>0.99) [47]. Therefore, the VIIRS-based evaporation rate values for other reservoirs should also be reliable and have good consistencies with those from MODIS.
Remote Sens. 2021, 13, x FOR PEER REVIEW Figure 10. Comparison of monthly evaporation rates for Lake Mead between MODIS and VIIRS.

Benefits of the MODIS/VIIRS Based Water Elevations as Compared to Radar Altim Products
Global reservoir elevation products are mainly based on measurements col satellite radar altimeters. We compared the MODIS 8-day elevation product w radar altimetry products (Hydroweb and G-REALM) in terms of reservoir cove temporal resolution ( Figure 11). Overall, the MODIS product and radar products agree well with each other (Figure 12). Given that the MODIS and VII products are very similar, Figure 11 also represents the relationship between t Figure 10. Comparison of monthly evaporation rates for Lake Mead between MODIS and VIIRS.

Benefits of the MODIS/VIIRS Based Water Elevations as Compared to Radar Altimetry Products
Global reservoir elevation products are mainly based on measurements collected by satellite radar altimeters. We compared the MODIS 8-day elevation product with two radar altimetry products (Hydroweb and G-REALM) in terms of reservoir coverage and temporal resolution ( Figure 11). Overall, the MODIS product and radar altimetry products agree well with each other (Figure 12). Given that the MODIS and VIIRS 8-day products are very similar, Figure 11 also represents the relationship between the VIIRS product and most radar altimetry products (except that the VIIRS start date is in 2012). Despite the similarities, the MODIS/VIIRS product suite can offer new benefits over traditional radar altimetry products in two ways. First, the MODIS/VIIRS products can provide reservoir water elevations at an 8-day temporal resolution. In comparison, Hydroweb data are at the monthly temporal resolution, and G-REALM data range from 10-day, to 27-day, to 35-day (depending on the sensor). Second, this global water reservoir product suite is characterized by continuous and consistent records since 2000 for all of the 164 reservoirs (some examples are shown in Figure 12

Limitations and Sources of Uncertainties
There are still some limitations that need to be addressed. First, only the relatively large reservoirs are included due to the moderate resolutions of MODIS/VIIRS. For smaller reservoirs, it is more appropriate to take advantage of missions with high spatial resolutions (e.g., Landsat and Sentinel-2). Second, some reservoirs are susceptible to mountain shadows (but this has only small effects on large reservoirs). This overestimation can be eliminated using topography information. For example, Li et al. [48] overlapped water classification maps with Shuttle Radar Topography Mission (SRTM) DEM data to exclude the mountain pixels.
The sources of uncertainties with regard to the reservoir surface area are associated with both the raw image classification using the reflectance product, and the classification enhancement algorithm. The accuracy of the Ostu classification of the NIR images is affected by the mixed pixels (i.e., partially covered by water and partially covered by land) at the reservoir boundaries, as well as by ice over the reservoirs. The reliability of the enhancement algorithm depends on the data quality of both the water occurrence image and the raw water classification. In high-latitude regions, the water occurrence image generally shows small surface area dynamics (i.e., the distribution of occurrence values highly skewed to the left). Thus, the pixels with low occurrence values have relatively large uncertainties. Under such conditions, the enhancement algorithm might generate an incorrect threshold value-typically leading to an overestimation of the surface area.   The reservoir elevation and storage estimation uncertainties include reservoir surface area uncertainties (see above), A-E relationship uncertainties, and reservoir configuration uncertainties. According to Equation (1), the estimated storage will be biased if the characteristics at capacity (storage, area, and elevation) are not accurate. Even when these factors have been correctly documented, the storage capacity may have changed due to sedimentation over time. Since the reservoir elevations are inferred only from areas and A-E relationships, they are not affected by reservoir configuration uncertainties. Moreover, while a non-linear relationship may perform better for some reservoirs, linear relationships were found to perform quite well for most of the reservoirs [33]. The linear A-E relationship represents a second-order polynomial reservoir cross-section, which is more realistic than other cross-section shapes [49].
Sources of evaporation rate uncertainty mainly include forcing data uncertainty and model structure/parameter uncertainty. Specifically, the forcing data used in this study (i.e., GLADS-2.1) are a land-based meteorological record. Although the increased humidity on the reservoir surface is represented by the wind function [28,50], differences in the wind speeds between reservoir and land regions are ignored, which might introduce some uncertainties [51]. In addition, the LTEM and its parameters can also produce uncertainties. For example, the formulation of the light attenuation coefficient (λ_PAR) is simplified. However, λ_PAR is affected by suspended solids, phytoplankton concentration level, and spectral distribution of solar radiation, and thus is constantly changing [52,53].

Limitations and Sources of Uncertainties
There are still some limitations that need to be addressed. First, only the relatively large reservoirs are included due to the moderate resolutions of MODIS/VIIRS. For smaller reservoirs, it is more appropriate to take advantage of missions with high spatial The reservoir volumetric evaporation uncertainty can be attributed to evaporation rate and surface area uncertainties, which have been discussed above.

Future Directions and Potential Applications
Several efforts are planned in the future to improve the current MODIS/VIIRS product suite. To reduce area estimation errors, the edge detection algorithm developed in Zhao et al. [29] will be employed to overcome the limitations of the current image enhancement approach. Furthermore, the design of the Terra (MOD) and Aqua (MYD) combined product components (MCD28C2 and MCD28C3) will allow it to automatically select less contaminated images for the users (thus leading to more accurate area estimates). It is also essential to conduct a comprehensive comparison between the MODIS and VIIRS products, especially given that the MODIS-era is close to an end. This analysis will be carried out after these products are made available through NASA. Additionally, a greater number of large reservoirs can be monitored by developing A-E relationships over locations where radar altimetry and ICESat data are unavailable. For instance, both Zhang and Gao [54] and Li et al. [55] have demonstrated the potential of using DEM data for generating A-E relationships for global reservoirs. Moreover, ICESat-2 (launched in September 2018) [48] and the Surface Water and Ocean Topography (SWOT) mission (planned for launch in 2022) [56] will provide more accurate elevation observations (as compared to radar altimetry and DEM). This should lead to the generation of A-E relationships over many more reservoirs.
The reservoir product can have multiple potential applications. First, it provides critical information for water resources management-especially for reservoirs in transboundary river basins where gauge observations are not shared. Moreover, it can help to calibrate and validate the operation rules incorporated into the global hydrological models that have a reservoir module. This can improve our understanding of the reservoirs' role in the hydrological cycle and water budgets [10,57]. For instance, Shah et al. [58] examined the role of reservoirs in water and energy budgets using the Variable Infiltration Capacity (VIC) model in the Indian sub-continent.
Additionally, time series values of reservoir storage variations at the basin scale can be used to evaluate hydrological droughts. Figure 13 shows the total reservoir storage variations in the Colorado and Murray-Darling river basins. The Colorado River basin has experienced a prolonged hydrological drought since the year 2000, even though the underlying meteorological drought only lasted for the first few years. This hydrological drought has been exacerbated by increased water use [59]. Moreover, the Murray-Darling river basin experienced the longest four-year reservoir-based hydrological drought on record from 2006 to 2010. However, the majority of the prior studies of droughts across the globe [60][61][62][63][64][65][66] considered only streamflow variations in their assessments of hydrological drought. This is because continuous long-term reservoir data with a high temporal resolution are lacking at a global scale (either for use in analysis or for the calibration of global hydrological models). Neglecting the reservoir storage component in hydrological drought assessment can lead to biases and uncertainties, as water is increasingly consumed, stored, and diverted through water management activities (i.e., reservoir operations) [66,67]. Therefore, this product can fill in this crucial research gap and can help stakeholders, water managers, and policymakers to improve existing hydrological drought management practices. temporal resolution are lacking at a global scale (either for use in analysis or for the calibration of global hydrological models). Neglecting the reservoir storage component in hydrological drought assessment can lead to biases and uncertainties, as water is increasingly consumed, stored, and diverted through water management activities (i.e., reservoir operations) [66,67]. Therefore, this product can fill in this crucial research gap and can help stakeholders, water managers, and policymakers to improve existing hydrological drought management practices.

Conclusions
This paper presents the development and validation results of NASA's global water reservoir product suite (MxD28 and VNP28), which includes 8-day and monthly water

Conclusions
This paper presents the development and validation results of NASA's global water reservoir product suite (MxD28 and VNP28), which includes 8-day and monthly water area, elevation, and storage values and monthly evaporation rate and volume values. This product consists of 164 reservoirs that represent nearly half of the total global capacity ( Figure 1 and Table 1). The water area values show good consistency with those from high-resolution Landsat measurements ( Figure 5). Validations against in situ observations over twelve reservoirs suggest a high level of accuracy for the elevation and storage products (Figures 6 and 7). Validations over Lake Nasser and Lake Mead indicate that the evaporation rate product also has a high quality ( Figure 8). Moreover, a comparison of 8-day reservoir area estimations between MODIS and VIIRS ( Figure 9) again shows good consistency, which ensures the continuation of this long-term reservoir product.
Compared with other remotely sensed reservoir products, the advantages of our product can be summarized as follows. First, the product provides consistent measurements for area, elevation, and storage and for evaporation rate and volume values. For most of the current products, the measurements provided are generally for a specific aspect. For example, the radar altimetry products tend to combine observations from multiple satellite missions, which may have data gaps and large uncertainties ( Figure 12). Second, the higher temporal resolution (8-day) can provide timely information useful for water management (e.g., water supply and flood mapping). Although Landsat can provide more accurate area estimations due to its higher spatial resolution, its longer revisit time (16 days) makes it not suitable for timely monitoring purpose. Furthermore, we have demonstrated that reservoir storage variations can be a good indicator of hydrological drought. This reservoir product can help to better understand the drought propagation process from the meteorological to hydrological perspectives.
Acknowledgments: This work has benefitted from the usage of the Google Earth Engine platform and the Texas A&M Supercomputing Facility (http://hprc.tamu.edu, accessed on 30 December 2020). We also thank the editors and three anonymous reviewers for their constructive comments.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript; or in the decision to publish the results.