Monitoring Reservoir Drought Dynamics with Landsat and Radar / Lidar Altimetry Time Series in Persistently Cloudy Eastern Brazil

Tropical reservoirs are critical infrastructure for managing drinking and irrigation water and generating hydroelectric power. However, long-term spaceborne monitoring of reservoir storage is challenged by data scarcity from near-persistent cloud cover and drought, which may reduce volumes below those in the observational record. In evaluating our ability to accurately monitor long-term reservoir volume dynamics using spaceborne data and overcome such observational challenges, we integrated optical, lidar, and radar time series to estimate reservoir volume dynamics across 13 reservoirs in eastern Brazil over a 12-year (2003–2014) period affected by historic drought. We (i) used 1560 Landsat images to measure reservoir surface area; (ii) built reservoir-specific regression models relating surface area and elevation from ICESat GLAS and Envisat RA-2 data; (iii) modeled volume changes for each reservoir; and (iv) compared modeled and in situ reservoir volume changes. Regression models had high goodness-of-fit (median RMSE = 0.89 m and r = 0.88) across reservoirs. Even though 88% of an average reservoir’s volume time series was based on modeled area–elevation relationships, we found exceptional agreement (RMSE = 0.31 km3 and r = 0.95) with in situ volume time series, and accurately captured seasonal recharge/depletion dynamics and the drought’s prolonged drawdown. Disagreements in volume dynamics were neither driven by wet/dry season conditions nor reservoir capacity, indicating analytical efficacy across a range of monitoring scenarios.


Introduction
Reservoirs are critical global infrastructure that occupy at least 26 × 10 4 km 2 (0.2%) of global land surface area and contribute 6 × 10 3 km 3 (0.0004%) of global freshwater storage [1].The freshwater held by reservoirs is essential for meeting global demand for drinking water, irrigation water for agriculture [2,3], and hydroelectric power generation [4].Recent research has identified the dual roles of reservoirs in atmospheric carbon dynamics being large-scale sites of both carbon sequestration as well as greenhouse gas emission [5][6][7][8][9][10][11]. Tropical reservoirs only make up 15% (1040) of the term monitoring of reservoir storage throughout drought as well as high precipitation periods is essential to understanding the role of tropical reservoirs on society and environment [8,16]).
Reservoir storage is temporally dynamic and satellite remote sensing time series from optical (e.g., Landsat, MODIS, Sentinel-2) and radar or lidar (e.g., Envisat, GLAS/ICESat) systems support systematic monitoring of changes in reservoir volume at low cost per reservoir [19][20][21][22][23][24][25][26][27].Spaceborne observation of reservoir dynamics are especially valuable for geographically remote or highelevation reservoirs where consistent and long-term monitoring are a challenge, or in regions where the cost of constructing or maintaining a hydrologic gauge network is prohibitively expensive [28][29][30].However, optical remote sensing of tropical reservoir dynamics is complicated by persistent cloud cover [31][32][33].Indeed, as shown in Figure 1, the 178 Brazilian reservoirs included in GRanD are on average obscured by clouds 54% of days between 2003 and 2014.[12].Reservoir counts are summarized at 5° resolution), and percent daily cloud cover is based on MODIS MOD09GA [34] (2003-2014, 1 km resolution, Robinson projection).
Using a case study of 13 reservoirs in eastern Brazil, the goal of this study is to accurately model tropical reservoir volume dynamics from 2003 to 2014 despite near-persistent cloud cover and reservoir depletion due to a historic drought in 2014.To achieve our goal, we estimated surface area  [12].Reservoir counts are summarized at 5 • resolution), and percent daily cloud cover is based on MODIS MOD09GA [34] (2003-2014, 1 km resolution, Robinson projection).
In addition to the challenges of persistent cloud cover, monitoring long-term tropical reservoir dynamics is complicated by sensor spatial resolution, satellite revisit period, and, in some cases, a limited observational record.MODIS has regularly been used to monitor reservoir dynamics at sub-monthly time scales (e.g., [23,25,35]).While spectral unmixing of MODIS imagery has been effective at measuring sub-pixel fractional water area [36,37], the relatively coarse 250 m spatial resolution impedes detection of the reservoir edge and quantification of fine-scale surface area changes (e.g., [38]) with heightened consequences for reservoirs with shallow sloped near-surface bathymetry [25].Landsat's 30 m spatial resolution, over 40 years of global coverage makes it well suited for long-term, high spatial resolution monitoring of reservoir surface water dynamics [39][40][41][42][43]. Unfortunately, Landsat's 16-day revisit period limits the opportunities to mitigate cloud cover or haze compared to 8-or 16-day MODIS temporal composites [44], which further restricts opportunities for pairing surface area estimates with concurrent altimetry elevation observations.While spaceborne radar is broadly robust to atmospheric effects, altimeters on TOPEX/Poseidon and Jason-1, 2, and 3 satellites with 3-10 km spatial resolutions are not well suited for monitoring smaller reservoirs [45,46].More recent platforms such as and Cryosat-2 and Sentinel 3 have higher spatial resolution but lack long-term coverage [24,47,48].However, Envisat RA-2 (Radar Altimeter 2) radar (2002-2012) and ICESat GLAS (Geoscience Laser Altimeter System) lidar (2003-2010) altimeters lend themselves to reservoir monitoring applications given their sub-kilometer resolution and many years of coverage [28,35,[49][50][51][52].Using a case study of 13 reservoirs in eastern Brazil, the goal of this study is to accurately model tropical reservoir volume dynamics from 2003 to 2014 despite near-persistent cloud cover and reservoir depletion due to a historic drought in 2014.To achieve our goal, we estimated surface area and elevation across study reservoirs by integrating Landsat, Envisat RA-2, and ICESat GLAS time series, modeling surface area-elevation relationship, calculating volume changes over the study period, and assessing agreement between modeled and in situ volume time series data per date, month, and year with attention to differences in agreement before and during the drought.

Study Area
Brazil is the largest (8.5 million km 2 ) and most populous country (approximately 209 million residents in 2017) in South America.The country is dominated by tropical rain forest in its west and north, has tropical semideciduous forest along its southeastern Atlantic coast with tropical savannas (cerrado) stretching between the Brazilian Highlands to the Mato Grosso Plateau, and is home to major urban population centers including São Paulo, Rio de Janeiro, and the capital, Brasilia, in the eastern half of the country.Brazil has the largest reserve of renewable surface water in the world, approximately 32% of which is used for agricultural production [53].Beginning in 2012, changing atmospheric circulation patterns, declining rainfall, and increased temperatures culminated in the summer of 2014 being the warmest and driest since 1951 [54] and a series of historically intense droughts that parched extensive cropland, dwindled drinking water supply, and shrank Brazil's rivers and reservoirs [55,56].
During the peak of the drought in the summer of 2014, national hydroelectric production declined by approximately 20% compared to average 2000-2010 levels, and the Cantareira reservoir system relied upon by São Paulo saw nearly 11% reduction in its total capacity [53].Over half of Brazil's largest reservoirs (by hydroelectric production) are located in the drought-affected eastern region of the country.Thirteen reservoirs were included in this study that span southeast Brazil (Figure 2) with at least three dates of observation between Envisat and GLAS.Four study reservoirs-Agua Vermelha, Furnas, Marimbondo, and Tres Marias-with a range in nominal storage capacity from 1.3 to 34.1 km 3 and a diversity of inter-annual storage dynamics were selected for focused data presentation below with the remaining reservoirs' data included in the Appendix A. and elevation across study reservoirs by integrating Landsat, Envisat RA-2, and ICESat GLAS time series, modeling surface area-elevation relationship, calculating volume changes over the study period, and assessing agreement between modeled and in situ volume time series data per date, month, and year with attention to differences in agreement before and during the drought.

Study Area
Brazil is the largest (8.5 million km 2 ) and most populous country (approximately 209 million residents in 2017) in South America.The country is dominated by tropical rain forest in its west and north, has tropical semideciduous forest along its southeastern Atlantic coast with tropical savannas (cerrado) stretching between the Brazilian Highlands to the Mato Grosso Plateau, and is home to major urban population centers including São Paulo, Rio de Janeiro, and the capital, Brasilia, in the eastern half of the country.Brazil has the largest reserve of renewable surface water in the world, approximately 32% of which is used for agricultural production [53].Beginning in 2012, changing atmospheric circulation patterns, declining rainfall, and increased temperatures culminated in the summer of 2014 being the warmest and driest since 1951 [54] and a series of historically intense droughts that parched extensive cropland, dwindled drinking water supply, and shrank Brazil's rivers and reservoirs [55,56].
During the peak of the drought in the summer of 2014, national hydroelectric production declined by approximately 20% compared to average 2000-2010 levels, and the Cantareira reservoir system relied upon by São Paulo saw nearly 11% reduction in its total capacity [53].Over half of Brazil's largest reservoirs (by hydroelectric production) are located in the drought-affected eastern region of the country.Thirteen reservoirs were included in this study that span southeast Brazil (Figure 2) with at least three dates of observation between Envisat and GLAS.Four study reservoirs-Agua Vermelha, Furnas, Marimbondo, and Tres Marias-with a range in nominal storage capacity from 1.3 to 34.1 km 3 and a diversity of inter-annual storage dynamics were selected for focused data presentation below with the remaining reservoirs' data included in the Appendix.  ) of study reservoirs in eastern Brazil [12].Four focus reservoirs whose volume dynamics are illustrated in the manuscript text are identified with a black point on the map and an asterisk (*) in the reservoir list; all reservoirs' dynamics are illustrated in Appendix A.

Modeling Framework
The volumetric modeling framework (Figure 3) for a given reservoir comprised five main objectives: (1) generate reservoir surface area time series; (2) generate surface elevation time series; (3) build and apply a linear regression model relating surface area and elevation; (4) estimate  3 ) of study reservoirs in eastern Brazil [12].Four focus reservoirs whose volume dynamics are illustrated in the manuscript text are identified with a black point on the map and an asterisk (*) in the reservoir list; all reservoirs' dynamics are illustrated in Appendix A.

Modeling Framework
The volumetric modeling framework (Figure 3) for a given reservoir comprised five main objectives: (1) generate reservoir surface area time series; (2) generate surface elevation time series; (3) build and apply a linear regression model relating surface area and elevation; (4) estimate volumetric change over the study period; and (5) assess agreement between modeled and in situ volumes.

Surface Area Time Series Generation
The 13 study reservoirs collectively span 15 Landsat WRS-2 path-row footprints, with three reservoirs stretching across 2-4 footprints, respectively.A total of 1460 Landsat 5, 7, and 8 surface reflectance (SR)-corrected scenes with 20% or less cloud cover were collected from 2003 through 2014; this cloud cover threshold reduced the likelihood of missing data due to cloud cover while still supporting an intra-annually dense time series.We measured surface water time series with the commonly used Modified Normalized Difference Water Index (MNDWI), which is effective at discriminating water from non-water features [57].We included all pixels that have a MDWI value larger than 0.1 as 'water' for each image date.Missing pixels associated with clouds, shadows, and Landsat 7's Scan Line Corrector (SLC) Error data gaps were identified and removed using Landsat's CFmask product.For each image's missing pixels, pixels from images collected within a 36-day window before/after (e.g., Figure 4a-c) were infilled into the image, prioritized by the nearest date of observation (e.g., Figure 4d).Infilling using a 36-day temporal window ensured inclusion of images from up to two along-track (every 16 days) and up to four across-track (7 or 9 days) observations and reduced the amount of missing data within an average reservoir from 8.6% to 1.2% (Figure 5).On some dates, an image's missing pixels could not be directly infilled through temporal compositing; in this case, if these missing pixels were in locations classified as 'water' on at least the median frequency for a given reservoir, the missing pixels were automatically included in the surface water extent.A reservoir's water frequency map showing the total number of days when 'water' was detected (Figure 4e) was also used to generate a consistent spatial extent within which surface water area was measured: pixels classified as 'water' on at least two dates made up an 'ever-water' envelope for this purpose (shown in Figure 4).Finally, each reservoir's surface area time series was clipped by an integrated (i.e., upstream to dam wall) reservoir region that was visually interpreted using very high-resolution imagery hosted by Google Earth.

Surface Area Time Series Generation
The 13 study reservoirs collectively span 15 Landsat WRS-2 path-row footprints, with three reservoirs stretching across 2-4 footprints, respectively.A total of 1460 Landsat 5, 7, and 8 surface reflectance (SR)-corrected scenes with 20% or less cloud cover were collected from 2003 through 2014; this cloud cover threshold reduced the likelihood of missing data due to cloud cover while still supporting an intra-annually dense time series.We measured surface water time series with the commonly used Modified Normalized Difference Water Index (MNDWI), which is effective at discriminating water from non-water features [57].We included all pixels that have a MDWI value larger than 0.1 as 'water' for each image date.Missing pixels associated with clouds, shadows, and Landsat 7's Scan Line Corrector (SLC) Error data gaps were identified and removed using Landsat's CFmask product.For each image's missing pixels, pixels from images collected within a 36-day window before/after (e.g., Figure 4a-c) were infilled into the image, prioritized by the nearest date of observation (e.g., Figure 4d).Infilling using a 36-day temporal window ensured inclusion of images from up to two along-track (every 16 days) and up to four across-track (7 or 9 days) observations and reduced the amount of missing data within an average reservoir from 8.6% to 1.2% (Figure 5).On some dates, an image's missing pixels could not be directly infilled through temporal compositing; in this case, if these missing pixels were in locations classified as 'water' on at least the median frequency for a given reservoir, the missing pixels were automatically included in the surface water extent.A reservoir's water frequency map showing the total number of days when 'water' was detected (Figure 4e) was also used to generate a consistent spatial extent within which surface water area was measured: pixels classified as 'water' on at least two dates made up an 'ever-water' envelope for this purpose (shown in Figure 4).Finally, each reservoir's surface area time series was clipped by an integrated (i.e., upstream to dam wall) reservoir region that was visually interpreted using very high-resolution imagery hosted by Google Earth.

Surface Elevation Time Series Generation
Separate surface water elevation time series were built using Envisat RA-2 (2002-2010) and ICESat GLAS (2003)(2004)(2005)(2006)(2007)(2008)(2009) altimetry data, respectively.The long wavelength of Envisat's RA-2 (Kuband) radar altimetry provides resilience to cloud cover effects while the vertical profiling capabilities of ICESat's GLAS support the discrimination of ground-level conditions in all but the densest cloud cover [31,58].Both RA-2 and GLAS datasets are suitable for reservoir monitoring with along-track

Surface Elevation Time Series Generation
Separate surface water elevation time series were built using Envisat RA-2 (2002-2010) and ICESat GLAS (2003)(2004)(2005)(2006)(2007)(2008)(2009) altimetry data, respectively.The long wavelength of Envisat's RA-2 (Kuband) radar altimetry provides resilience to cloud cover effects while the vertical profiling capabilities of ICESat's GLAS support the discrimination of ground-level conditions in all but the densest cloud cover [31,58].Both RA-2 and GLAS datasets are suitable for reservoir monitoring with along-track
The Envisat RA-2 surface elevation time series was generated using an automated return clustering algorithm detailed in [74], which yielded an average RMSE of 48.9 cm compared to in situ gauge observation data across study reservoirs.ICESat GLAS GLA14 product (Level-2 Global Land Surface Altimetry; [75]) data were filtered to single Gaussian peak returns typical of surface water, saturated returns were corrected or removed using the saturation index, and surface elevation was calculated at the centroid of the single Gaussian peak.All RA-2 and GLAS surface returns collected within 15 days of a Landsat surface area measurement were spatially filtered by the surface area perimeter buffered inwards by 100 m to eliminate potential measurement of near-surface topography.The median elevation over a reservoir's extent on a given date was measured and used to generate reservoir elevation and volume time series.Since GLAS' reference ellipsoid (equatorial radius = 6378.1363km; flattening coefficient = 1/298.257)is offset 70 cm from RA-2's WGS 84 ellipsoid, GLAS elevation data were adjusted to fit the WGS 84 vertical datum following [76].

Surface Area-Elevation Model Generation
Surface area and elevation time series were used to estimate changes in reservoir volumes.However, since the vast majority of surface area measurements lacked same-day measurements of surface elevation, regression models were built for each reservoir to relate surface area measurements to altimeter elevations collected within 15 days of the surface area measurement; linear regression was used following successful regional application (e.g., [23,25,49]).Modeled surface elevation values based on measured surface areas were used to build an elevation time series with the same temporal sampling as the surface area time series.The goodness-of-fit of each reservoir's area-elevation model was measured using the Pearson correlation coefficient (r) and the root mean square error (RMSE).

Volume Time Series Generation
As the volume of a reservoir with unknown bathymetry cannot be remotely measured, the reservoir was rather conceived as a circular cone with a series of water layers as detailed in [77].Reservoir volume on a given date, t, was based on changes in surface area and elevation time series data measured relative to the median 2003-2011 volume (before the drought): where, V med-t represents the relative difference in volume, A med and A t and z med and z t represent surface areas and elevations based on median 2003-2011 values and the date of observation, respectively.The median 2003-2011 in situ volume was measured for each reservoir to serve as a pre-drought baseline, and the relative volume difference for each date of in situ observation was measured.

Comparison between In Situ and Modeled Volumetric Time Series
Changes in modeled reservoir volumes were compared to a median 4383 daily in situ reservoir volume measurements collected from 2003-2014 by the Brazilian Electric Sector available at http: //www.ons.org.br/.However, study reservoir boundaries deviated from (unavailable) management boundaries used by the Brazilian Electric Sector to collect in situ volume and modeled and in situ volume time series could not be directly compared.To mitigate this disparity, both time series were standardized to preserve relative variations with the following equation: where S t and V t are standardized and modeled volume changes for a given reservoir on date, t, respectively, and V and σ are the mean and standard deviation of modeled volume changes over the full time series, respectively.Using standardized (unitless) values, the agreement between modeled and in situ volume dynamics time series were assessed using linear regression, r, and RMSE, and annual and monthly agreement between modeled and in situ data were evaluated.

Surface Area and Elevation Time Series
Over the 12-year study period, a median 216 Landsat image dates were used to generate each reservoir surface area time series (Figures 6 and A1).Prior to the drought's intensification in 2014, reservoirs typically reached their lowest extent in December or January before recharging to the maximum surface area between March and June following the end of the rainy season.Surface water expansion and contraction was highly dynamic with a median annual range (i.e., difference between maximum and minimum) of reservoir surface areas representing 31.3% of the maximum surface area in the pre-drought 2003-2011 period.During the drought period (2012-2014), this range was nominally consistent at 30.7% though the minimum and maximum surface areas declined relative to the pre-drought period by median 4.5% and 3.6%, respectively.Seasonal dynamics shifted under the drought, as well, as early year recharge was typically reduced in 2012 and 2013 and absent by 2014.Reservoirs such as Furnas and Tres Marias contracted during the drought so much that low surface areas (below 900 km 2 and 350 km 2 , respectively) were without precedent in the pre-drought time series.where St and Vt are standardized and modeled volume changes for a given reservoir on date, t, respectively, and  and σ are the mean and standard deviation of modeled volume changes over the full time series, respectively.Using standardized (unitless) values, the agreement between modeled and in situ volume dynamics time series were assessed using linear regression, r, and RMSE, and annual and monthly agreement between modeled and in situ data were evaluated.

Surface Area and Elevation Time Series
Over the 12-year study period, a median 216 Landsat image dates were used to generate each reservoir surface area time series (Figures 6 and A1).Prior to the drought's intensification in 2014, reservoirs typically reached their lowest extent in December or January before recharging to the maximum surface area between March and June following the end of the rainy season.Surface water expansion and contraction was highly dynamic with a median annual range (i.e., difference between maximum and minimum) of reservoir surface areas representing 31.3% of the maximum surface area in the pre-drought 2003-2011 period.During the drought period (2012-2014), this range was nominally consistent at 30.7% though the minimum and maximum surface areas declined relative to the pre-drought period by median 4.5% and 3.6%, respectively.Seasonal dynamics shifted under the drought, as well, as early year recharge was typically reduced in 2012 and 2013 and absent by 2014.Reservoirs such as Furnas and Tres Marias contracted during the drought so much that low surface areas (below 900 km 2 and 350 km 2 , respectively) were without precedent in the pre-drought time series.While each reservoir's seasonal dynamics are apparent, there was often a disproportionate exclusion of cloudy wet season imagery (October-March; Figures 7a and A2).Conversely, since cloud-free observations are more likely in eastern Brazil during winter months, there is a relative abundance of imagery available in July and August when reservoirs tend to be near capacity.Fewer images and generally more cloud cover during the wet season meant fewer near-date images for While each reservoir's seasonal dynamics are apparent, there was often a disproportionate exclusion of cloudy wet season imagery (October-March; Figures 7a and A2).Conversely, since cloud-free observations are more likely in eastern Brazil during winter months, there is a relative abundance of imagery available in July and August when reservoirs tend to be near capacity.Fewer images and generally more cloud cover during the wet season meant fewer near-date images for infilling and compositing, which may yield a larger amount of residual missing data and low biased surface area estimates [33], as well as fewer surface area measurements with which to couple elevation data.While the collection of surface elevation values was not affected by atmospheric conditions in the same way as surface area measurements, having a larger surface area (unobstructed by clouds) meant that more altimetric returns could be collected over the reservoir body (Figures 7b and A3).surface area estimates [33], as well as fewer surface area measurements with which to couple elevation data.While the collection of surface elevation values was not affected by atmospheric conditions in the same way as surface area measurements, having a larger surface area (unobstructed by clouds) meant that more altimetric returns could be collected over the reservoir body (Figures 7b  and A3).

Modeled Surface Area-Elevation Relationships
Linear regression models relating Landsat-derived surface area and altimeter-derived surface elevation measurements were built for the 13 study reservoirs (Figures 8 and A4).The median number of area-elevation pairs for "combined" regression models was 21, with 13 and 17 measurements for GLAS-and RA-2-specific models, respectively.Combined linear models showed high goodness-of-fit with a median RMSE of 0.89 m and r of 0.88 across reservoirs.A poor model may reflect a combination of seasonal variation in surface area and elevation, differences in Landsat and altimetry acquisition dates, variation in bathymetric slope along a reservoir's perimeter and across its depth, as well as lingering cloud cover effects.Though GLAS and RA-2 elevations span a median 74.9% of reservoir surface area range (highlighted in grey in Figures 6 and A1), the elevation distribution is often skewed towards including higher elevations.This skew is most pronounced for reservoirs with fewer dates of observation during the wet season as well as those with a smaller surface area, which limits the area within which altimetry data could be collected.The contraction of surface areas during the drought also limited potential for coupling surface area-elevation values as very low surface areas during the drought had no paired elevation values.

Modeled Surface Area-Elevation Relationships
Linear regression models relating Landsat-derived surface area and altimeter-derived surface elevation measurements were built for the 13 study reservoirs (Figures 8 and A4).The median number of area-elevation pairs for "combined" regression models was 21, with 13 and 17 measurements for GLAS-and RA-2-specific models, respectively.Combined linear models showed high goodness-of-fit with a median RMSE of 0.89 m and r of 0.88 across reservoirs.A poor model may reflect a combination of seasonal variation in surface area and elevation, differences in Landsat and altimetry acquisition dates, variation in bathymetric slope along a reservoir's perimeter and across its depth, as well as lingering cloud cover effects.Though GLAS and RA-2 elevations span a median 74.9% of reservoir surface area range (highlighted in grey in Figures 6 and A1), the elevation distribution is often skewed towards including higher elevations.This skew is most pronounced for reservoirs with fewer dates of observation during the wet season as well as those with a smaller surface area, which limits the area within which altimetry data could be collected.The contraction of surface areas during the drought also limited potential for coupling surface area-elevation values as very low surface areas during the drought had no paired elevation values.GLAS-specific models provided a slightly better fit with median RMSE = 0.68 m and r = 0.94 compared to RA-2's RMSE = 0.94 m and r = 0.90 (see Table A1).The better fit provided by GLAS is expected given its higher vertical accuracy and smaller footprint.For the 8 of 11 reservoirs with both GLAS and RA-2 coverage over the study period, GLAS also captured a larger surface area range despite GLAS having fewer samples for these reservoirs (see Figure A4).Further, GLAS alone provided coverage for two reservoirs, Ilha Solteira and Sao Simao.Despite the driving role of GLAS data in shaping surface area-elevation models, there remains considerable value in combining altimeter measurements into a single area-elevation model.Since GLAS and RA-2 captured reservoirs at different stages over the time series, combining altimetric observations in a single model expanded the seasonality of coverage, which helped compensate for irregularity of observations at a given reservoir.Combining altimetric data also increased the average range of observed surface elevations across reservoirs by at least 1.3 m compared to the elevation range observed by either GLAS or RA-2 alone.This expanded range supports characterizing reservoir dynamics from the drawdown through recharge conditions, which is especially valuable for reservoirs with a large annual range, e.g., Agua Vermelha.3.3.Volume Time Series and In Situ Comparison.
Across reservoirs, the median modeled volume time series is composed of 216 time steps (Figures 9a and A5a).Of these, 88% of time series values were based on reservoir-specific surface area-elevation regression models (i.e., Figure 8) while the remaining 12% were based on direct observations.While some reservoirs such as Agua Vermelha have a clear, consistent seasonality with a pronounced annual recharge and depletion, other reservoirs such as Chavantes or Sao Simao (see Figure A5) are much less dynamic.Linear regression models showed exceptionally high agreement between modeled and in situ volumetric changes with median RMSE = 0.31 km 3 and r = 0.95 (Figure 9b).The low median bias of 0.11 km 3 indicates an underestimation of in situ volume, and, indeed, each reservoir's maximum modeled volume was consistently less than the maximum in situ volume.This divergence likely results from the mismatch between reservoir boundaries used to generate modeled values and boundaries considered in in situ data collection, differing 2003-2011 reference volumes, and spuriously low modeled volumes that most often occur during the wet season.GLAS-specific models provided a slightly better fit with median RMSE = 0.68 m and r = 0.94 compared to RA-2's RMSE = 0.94 m and r = 0.90 (see Table A1).The better fit provided by GLAS is expected given its higher vertical accuracy and smaller footprint.For the 8 of 11 reservoirs with both GLAS and RA-2 coverage over the study period, GLAS also captured a larger surface area range despite GLAS having fewer samples for these reservoirs (see Figure A4).Further, GLAS alone provided coverage for two reservoirs, Ilha Solteira and Sao Simao.Despite the driving role of GLAS data in shaping surface area-elevation models, there remains considerable value in combining altimeter measurements into a single area-elevation model.Since GLAS and RA-2 captured reservoirs at different stages over the time series, combining altimetric observations in a single model expanded the seasonality of coverage, which helped compensate for irregularity of observations at a given reservoir.Combining altimetric data also increased the average range of observed surface elevations across reservoirs by at least 1.3 m compared to the elevation range observed by either GLAS or RA-2 alone.This expanded range supports characterizing reservoir dynamics from the drawdown through recharge conditions, which is especially valuable for reservoirs with a large annual range, e.g., Agua Vermelha.3.3.Volume Time Series and In Situ Comparison.
Across reservoirs, the median modeled volume time series is composed of 216 time steps (Figures 9a and A5a).Of these, 88% of time series values were based on reservoir-specific surface area-elevation regression models (i.e., Figure 8) while the remaining 12% were based on direct observations.While some reservoirs such as Agua Vermelha have a clear, consistent seasonality with a pronounced annual recharge and depletion, other reservoirs such as Chavantes or Sao Simao (see Figure A5) are much less dynamic.Linear regression models showed exceptionally high agreement between modeled and in situ volumetric changes with median RMSE = 0.31 km 3 and r = 0.95 (Figure 9b).The low median bias of 0.11 km 3 indicates an underestimation of in situ volume, and, indeed, each reservoir's maximum modeled volume was consistently less than the maximum in situ volume.This divergence likely results from the mismatch between reservoir boundaries used to generate To examine intra-annual variation in agreement, the mean absolute difference between modeled and in situ volumetric change was estimated for each month (Figures 10 and A6).Wet season months (October-March) showed higher disagreement than winter months (July-August) in 5 of 13 study reservoirs.Since reservoirs tend to have their lowest volumes during wet months of November-January, accurately modeling reservoir volumes during these low volume months would improve modeling of anomalously low volumes during drought periods.Considering inter-annual variation, the mean annual flux (i.e., the range of standardized volumes within a calendar year) of modeled and in situ time series only differed between 0.004 and 0.041 in 2003-2010 across reservoirs (Figure 11).With the drought came a larger disagreement of 0.62 in 2012 but was reduced to a very low disagreement of only 0.04 in 2014.The divergence between modeled and in situ flux in 2011-2013 was likely influenced by fewer Landsat images being collected after Landsat 5 operational imaging ended in November 2011, 17 months before Landsat 8's operational imaging began in April 2013.The agreement in 2014 benefitted from Landsat 8's data collection, a low annual flux due to the drought, and good fitting area-elevation regression models that allowed for estimating elevations beyond the range of historical direct observation.To examine intra-annual variation in agreement, the mean absolute difference between modeled and in situ volumetric change was estimated for each month (Figures 10 and A6).Wet season months (October-March) showed higher disagreement than winter months (July-August) in 5 of 13 study reservoirs.Since reservoirs tend to have their lowest volumes during wet months of November-January, accurately modeling reservoir volumes during these low volume months would improve modeling of anomalously low volumes during drought periods.Considering inter-annual variation, the mean annual flux (i.e., the range of standardized volumes within a calendar year) of modeled and in situ time series only differed between 0.004 and 0.041 in 2003-2010 across reservoirs (Figure 11).With the drought came a larger disagreement of 0.

Discussion
Accurate long-term monitoring of reservoir dynamics achieved by this study is of paramount interest for assessing United Nations Development Program Sustainable Development Goal (SDG) 6 (Clean Water and Sanitation) and 7 (Affordable and Clean Energy) [78], managing national-level water resources, examining reservoir dynamics and hydroelectric power generation (e.g., [3]), assessing the impact of changing precipitation and temperature regimes on reservoir storage [54, 56,79,80], and as input to or validation of hydrodynamic models (e.g., [81][82][83].This study's high agreement between remote sensing and in situ volumetric flux (i.e., RMSE = 0.31 km 3 ) illustrates the potential for effective monitoring in tropical regions with persistent cloud cover or even those affected by historic drought.
Volumetric modeling is an essential tool in monitoring climate effects on global surface water dynamics especially in regions with unavailable or scarce in situ data [43].This study builds on research by [19,20,[22][23][24][25]65,84,85] in pursuing remote sensing-based estimation of reservoir volume dynamics.However, this research alongside [49] extends previous efforts by maintaining high accuracy in modeled volume changes outside the historical range of altimetric observations, in this case due to historic drought.Though the drought began in 2012, the reduction of reservoir surface areas, elevations, and volumes was broadly delayed across study reservoirs until 2013 as is evident in modeled and in situ volumetric change time series (e.g., Figure A5).In 2013, surface elevations at eight of the study reservoirs fell below the observed minimum surface elevation from 2000-2009; by

Discussion
Accurate long-term monitoring of reservoir dynamics achieved by this study is of paramount interest for assessing United Nations Development Program Sustainable Development Goal (SDG) 6 (Clean Water and Sanitation) and 7 (Affordable and Clean Energy) [78], managing national-level water resources, examining reservoir dynamics and hydroelectric power generation (e.g., [3]), assessing the impact of changing precipitation and temperature regimes on reservoir storage [54, 56,79,80], and as input to or validation of hydrodynamic models (e.g., [81][82][83].This study's high agreement between remote sensing and in situ volumetric flux (i.e., RMSE = 0.31 km 3 ) illustrates the potential for effective monitoring in tropical regions with persistent cloud cover or even those affected by historic drought.
Volumetric modeling is an essential tool in monitoring climate effects on global surface water dynamics especially in regions with unavailable or scarce in situ data [43].This study builds on research by [19,20,[22][23][24][25]65,84,85] in pursuing remote sensing-based estimation of reservoir volume dynamics.However, this research alongside [49] extends previous efforts by maintaining high accuracy in modeled volume changes outside the historical range of altimetric observations, in this case due to historic drought.Though the drought began in 2012, the reduction of reservoir surface areas, elevations, and volumes was broadly delayed across study reservoirs until 2013 as is evident in modeled and in situ volumetric change time series (e.g., Figure A5).In 2013, surface elevations at eight of the study reservoirs fell below the observed minimum surface elevation from 2000-2009; by

Discussion
Accurate long-term monitoring of reservoir dynamics achieved by this study is of paramount interest for assessing United Nations Development Program Sustainable Development Goal (SDG) 6 (Clean Water and Sanitation) and 7 (Affordable and Clean Energy) [78], managing national-level water resources, examining reservoir dynamics and hydroelectric power generation (e.g., [3]), assessing the impact of changing precipitation and temperature regimes on reservoir storage [54, 56,79,80], and as input to or validation of hydrodynamic models (e.g., [81][82][83].This study's high agreement between remote sensing and in situ volumetric flux (i.e., RMSE = 0.31 km 3 ) illustrates the potential for effective monitoring in tropical regions with persistent cloud cover or even those affected by historic drought.
Volumetric modeling is an essential tool in monitoring climate effects on global surface water dynamics especially in regions with unavailable or scarce in situ data [43].This study builds on research by [19,20,[22][23][24][25]65,84,85] in pursuing remote sensing-based estimation of reservoir volume dynamics.However, this research alongside [49] extends previous efforts by maintaining high accuracy in modeled volume changes outside the historical range of altimetric observations, in this case due to historic drought.Though the drought began in 2012, the reduction of reservoir surface areas, elevations, and volumes was broadly delayed across study reservoirs until 2013 as is evident in modeled and in situ volumetric change time series (e.g., Figure A5).In 2013, surface elevations at eight of the study reservoirs fell below the observed minimum surface elevation from 2000-2009; by 2014, 12 of the 13 study reservoirs had done so.Despite extrapolating surface area-elevation relationships to capture volumetric change in 2013-2014, the mean standardized modeled and in situ volumetric fluxes (e.g., Figure 11) were in agreement, thereby demonstrating the effectiveness of monitoring.Volume estimate accuracies for Tres Marias, Marimbondo, and others are comparable to those achieved by [19] (i.e., r = 0.99) despite this study's challenges of persistent cloud cover and extreme volumetric loss.
Limited observational frequency remained a challenge in this study.The largest difference between mean annual standardized modeled and in situ volumetric fluxes was 0.62 (unitless) and measured in 2012.This did not result from drought-induced depletion but was rather associated with reduced surface area measurements since only Landsat 7 was operational in 2012.Nonetheless, despite intra-annual variability of surface area measurements, limited observational data, and drought-induced depletion of reservoir volume beyond the range of observations, agreement between modeled and in situ volume dynamics was high.This is perhaps best exemplified by Tres Marias, which lacked any altimetry observations for six months the calendar year yet shows better agreement with in situ volume data (i.e., RMSE = 0.14 km 3 ) than any other reservoir.Future volumetric modeling efforts will benefit from the increased observational frequency by including additional surface area time series from Sentinel-1, Sentinel-2, as well as the upcoming Surface Water and Ocean Topography (SWOT) mission [86].Similarly, while no current hydrodynamic models accurately simulates reservoirs and lakes, recent developments on hydrological modeling have taken into account reservoir operation (e.g., [87][88][89]).Results presented in this study can therefore be used for further model evaluation and improvement of model parameters to better model human activities.
Errors inherent to the volumetric modeling approach include Landsat-derived surface area resolution, altimetry-derived vertical resolution, area-elevation linear regression modeling, and unavailable reservoir boundaries used for in situ data collection.Moreover, referencing the median water frequency coverage for infilling remaining pixels that could not be directly infilled with temporal compositing lead to overestimating surface area and volume on low volume days and underestimating surface area and volume on high volume days; this, in turn, degraded surface area-elevation correlations.While a linear area-elevation relationship that assumes a consistent bathymetric (i.e., hypsographic) profile across all depths and along the entire boundary of a given shoreline has been used effectively by other researchers (e.g., [23]), non-linear models have also been used to model sections or the entirety of a reservoir's hypsographic profile (e.g., [36,85]).A more extensive examination of the relationship between reservoir dynamics and bathymetry would have added valuable information.Unfortunately, there were no available reference bathymetry datasets with which to evaluate this relationship, nor were there surface area validation time series data with which to evaluate the accuracy of surface area dynamics for any reservoir in the study.Including information from a high-resolution bathymetric model would have improved volumetric modeling effort [84].
There are several opportunities for improving the water classification through use of a more advanced classification approach sensitive to reservoir-specific hydrologic conditions such as clarity, turbidity, and bathymetry, e.g., [90].The results of this study point towards more effective monitoring using recently launched sensors such as Sentinel 1 and 2 [91][92][93].The use of 20 m spatial resolution Sentinel 2A and 2B time series data with a nominal 5-day temporal resolution rather than 30 m, 16-day repeat Landsat data have the potential to estimate much more detailed volumetric changes, especially for reservoirs with complex shorelines that cannot be readily captured with Landsat, let alone MODIS imagery.Use of imaging radars such Sentinel 1, SWOT, and RADARSAT-2 [94] with <100 m horizontal resolution and <10 cm vertical resolution as well as passive microwave sensors [95] offer data fusion opportunities that make the most of available optical and passive time series data.

Conclusions
This study presented a remote sensing-based approach for monitoring reservoir dynamics over a 12-year period (2003-2014), which included a historic two-year long drought.Using the case study of regularly clouded eastern Brazil, Landsat 5, 7, 8, ICESat GLAS, and Envisat RA-2 time series data were used to measure surface area, elevation, and volume time series for 13 study reservoirs of varying size and intra-and inter-annual dynamics.Volumetric changes have high overall agreement with in situ volumetric time series data as well as intra-and inter-annual agreement throughout the study period including the drought.The high agreement despite regular cloud obfuscation and drought depletion points the way towards routine, very high-resolution monitoring of reservoir dynamics in regions that lack existing or available in situ monitoring.While this study benefits from reservoir management data, anticipated applications extend to regularly clouded regions without existing or otherwise inaccessible data on reservoir dynamics.A2.Stacked monthly frequency of Landsat images for all 13 study reservoirs.Monthly frequency information for four selected reservoirs is in Figure 7.

Figure 2 .
Figure2.Geographic distribution and storage capacities (km 3 ) of study reservoirs in eastern Brazil[12].Four focus reservoirs whose volume dynamics are illustrated in the manuscript text are identified with a black point on the map and an asterisk (*) in the reservoir list; all reservoirs' dynamics are illustrated in Appendix A.

Figure 2 .
Figure2.Geographic distribution and storage capacities (km 3 ) of study reservoirs in eastern Brazil[12].Four focus reservoirs whose volume dynamics are illustrated in the manuscript text are identified with a black point on the map and an asterisk (*) in the reservoir list; all reservoirs' dynamics are illustrated in Appendix A.
Remote Sens. 2019, 11, x FOR PEER REVIEW 4 of 25volumetric change over the study period; and (5) assess agreement between modeled and in situ volumes.

Figure 3 .
Figure 3. Overview of analytical framework to generate surface area, elevation, and volumetric time series for a given reservoir.Colors shown for each dataset are consistent throughout all figures.

Figure 3 .
Figure 3. Overview of analytical framework to generate surface area, elevation, and volumetric time series for a given reservoir.Colors shown for each dataset are consistent throughout all figures.

Figure 4 .
Figure 4. Example generation of temporally composited, infilled surface water extent maps showing southeastern Chavantes reservoir.(a-c) Surface water coverage from three near-date images (23 October, 8 November, and 24 November 2012) with data gaps.(d) Coverage maps are composited to create an infilled surface water map for 8 November 2012.(e) A time series of infilled water maps are used to generate a water frequency map for 2003-2014.

Figure 5 .
Figure 5. Reductions in average percentage of missing data across all image dates and reservoirs after infilling.

Figure 4 .
Figure 4. Example generation of temporally composited, infilled surface water extent maps showing southeastern Chavantes reservoir.(a-c) Surface water coverage from three near-date images (23 October, 8 November, and 24 November 2012) with data gaps.(d) Coverage maps are composited to create an infilled surface water map for 8 November 2012.(e) A time series of infilled water maps are used to generate a water frequency map for 2003-2014.

Figure 4 .
Figure 4. Example generation of temporally composited, infilled surface water extent maps showing southeastern Chavantes reservoir.(a-c) Surface water coverage from three near-date images (23 October, 8 November, and 24 November 2012) with data gaps.(d) Coverage maps are composited to create an infilled surface water map for 8 November 2012.(e) A time series of infilled water maps are used to generate a water frequency map for 2003-2014.

Figure 5 .
Figure 5. Reductions in average percentage of missing data across all image dates and reservoirs after infilling.

Figure 5 .
Figure 5. Reductions in average percentage of missing data across all image dates and reservoirs after infilling.

Figure 6 .
Figure 6.Landsat-derived surface area time series for four reservoirs with dates of ICESat GLAS (green dot) and Envisat RA-2 (red dot) measurement.The gray background indicates the range of surface areas measured on altimetry data collection dates.Time series for all 13 reservoirs are shown in Figure A1.

Figure 6 .
Figure 6.Landsat-derived surface area time series for four reservoirs with dates of ICESat GLAS (green dot) and Envisat RA-2 (red dot) measurement.The gray background indicates the range of surface areas measured on altimetry data collection dates.Time series for all 13 reservoirs are shown in Figure A1.

Figure 7 .
Figure 7. Monthly frequency of (a) Landsat images and (b) GLAS and RA-2 surface elevation observations for four selected reservoirs from 2003-2014.Monthly frequency of Landsat images and surface elevation observations for all 13 reservoirs are shown in Figures A2 and A3, respectively.

Figure 7 .
Figure 7. Monthly frequency of (a) Landsat images and (b) GLAS and RA-2 surface elevation observations for four selected reservoirs from 2003-2014.Monthly frequency of Landsat images and surface elevation observations for all 13 reservoirs are shown in Figures A2 and A3, respectively.

Figure 8 .
Figure 8. Derived surface area-elevation linear regressions for four selected reservoirs based on Landsat surface area and combined GLAS and RA-2 elevation data.n = number of area-elevation pairs, r = Pearson correlation coefficient, RMSE = root mean square error, m = linear slope, and b = linear intercept.Area-elevation models for all 13 reservoirs are in Figure A4.

Figure 8 .
Figure 8. Derived surface area-elevation linear regressions for four selected reservoirs based on Landsat surface area and combined GLAS and RA-2 elevation data.n = number of area-elevation pairs, r = Pearson correlation coefficient, RMSE = root mean square error, m = linear slope, and b = linear intercept.Area-elevation models for all 13 reservoirs are in Figure A4.

25 Figure 9 .
Figure 9. (a) Time series of difference between standardized reservoir volumes and each reservoir's baseline, pre-drought (i.e., median 2003-2011) volume for modeled (purple) and in situ (green) data, respectively; for dates without altimetry data, regression-based elevation values are used in modeling volume.(b) Linear regressions relating standardized in situ and modeled volumes on dates of mutual observation.Light grey lines in (a) and (b) indicate zero values on respective axes.All volume values are in km 3 .Time series and regressions for all 13 reservoirs are in Figure A5.

Figure 9 .
Figure 9. (a) Time series of difference between standardized reservoir volumes and each reservoir's baseline, pre-drought (i.e., median 2003-2011) volume for modeled (purple) and in situ (green) data, respectively; for dates without altimetry data, regression-based elevation values are used in modeling volume.(b) Linear regressions relating standardized in situ and modeled volumes on dates of mutual observation.Light grey lines in (a) and (b) indicate zero values on respective axes.All volume values are in km 3 .Time series and regressions for all 13 reservoirs are in Figure A5.

25 Figure 10 .
Figure 10.Mean monthly absolute difference between modeled and in situ standardized (unitless) volume changes (2003-2014; red line) with standard deviation range (red field).Monthly comparisons for all 13 reservoirs are in Figure A6.

Figure 11 .
Figure 11.Comparison between mean annual flux of standardized (unitless) modeled (purple line) and in situ (green line) volume dynamics across all study reservoirs with ±1 standard deviation ranges (purple and green fields, respectively).

Figure 10 . 25 Figure 10 .
Figure 10.Mean monthly absolute difference between modeled and in situ standardized (unitless) volume changes (2003-2014; red line) with standard deviation range (red field).Monthly comparisons for all 13 reservoirs are in Figure A6.

Figure 11 .
Figure 11.Comparison between mean annual flux of standardized (unitless) modeled (purple line) and in situ (green line) volume dynamics across all study reservoirs with ±1 standard deviation ranges (purple and green fields, respectively).

Figure 11 .
Figure 11.Comparison between mean annual flux of standardized (unitless) modeled (purple line) and in situ (green line) volume dynamics across all study reservoirs with ±1 standard deviation ranges (purple and green fields, respectively).

Figure A1 .
Figure A1.Landsat-derived surface area time series for all 13 study reservoirs with dates of ICESat GLAS (green dot) and Envisat RA-2 (red dot) measurement.The gray background indicates the range of surface areas measured on altimetry data collection dates.Time series for select four reservoirs are shown in Figure 6.

Figure A1 . 25 Figure A2 .
Figure A1.Landsat-derived surface area time series for all 13 study reservoirs with dates of ICESat GLAS (green dot) and Envisat RA-2 (red dot) measurement.The gray background indicates the range of surface areas measured on altimetry data collection dates.Time series for select four reservoirs are shown in Figure 6.Remote Sens. 2019, 11, x FOR PEER REVIEW 15 of 25

Figure A3 .
Figure A3.Stacked monthly frequency of GLAS and RA-2 surface elevation observations for all 13 study reservoirs.Monthly frequency information for four selected reservoirs is in Figure 7.

Figure A3 .
Figure A3.Stacked monthly frequency of GLAS and RA-2 surface elevation observations for all 13 study reservoirs.Monthly frequency information for four selected reservoirs is in Figure 7.

Figure A4 .
Figure A4.Derived surface area-elevation linear regressions for all 13 study reservoirs based on Landsat surface area and combined GLAS and RA-2 elevation data.n = number of area-elevation pairs, r = Pearson correlation coefficient, RMSE = root mean square error, m = linear slope, and b = linear intercept.Area-elevation models for select four reservoirs are in Figure 8.

Figure A4 .
Figure A4.Derived surface area-elevation linear regressions for all 13 study reservoirs based on Landsat surface area and combined GLAS and RA-2 elevation data.n = number of area-elevation pairs, r = Pearson correlation coefficient, RMSE = root mean square error, m = linear slope, and b = linear intercept.Area-elevation models for select four reservoirs are in Figure 8.

Figure A5 .
Figure A5.(a) Time series of difference between standardized reservoir volumes and each reservoir's baseline, pre-drought (i.e., median 2003-2011) volume for modeled (purple) and in situ (green) data, respectively; for dates without altimetry data, regression-based elevation values are used in modeling volume.(b) Linear regressions relating standardized in situ and modeled volumes on dates of mutual Figure A5.(a) Time series of difference between standardized reservoir volumes and each reservoir's baseline, pre-drought (i.e., median 2003-2011) volume for modeled (purple) and in situ (green) data, respectively; for dates without altimetry data, regression-based elevation values are used in modeling volume.(b) Linear regressions relating standardized in situ and modeled volumes on dates of mutual observation.Light grey lines in (a) and (b) indicate zero values on respective axes.All volume values are in km 3 .Time series and regressions for four select reservoirs are in Figure 9.

Figure A6 .
Figure A6.Mean monthly absolute difference between modeled and in situ standardized (unitless) volume changes (2003-2014; red line) with standard deviation range depicted (red field).Monthly comparisons for four select study reservoirs are in Figure 10.

Figure A6 .
Figure A6.Mean monthly absolute difference between modeled and in situ standardized (unitless) volume changes (2003-2014; red line) with standard deviation range depicted (red field).Monthly comparisons for four select study reservoirs are in Figure 10.

Table A1 .
List of reservoir-specific linear regression model parameters and goodness-of-fit values for GLAS-or RA-2-specific models as well as 'combined' models.n = number of altimeter surface elevation measurements, m = linear slope, b = linear intercept, r = Pearson correlation coefficient, and RMSE = root mean square error.