Identifying Emerging Reservoirs along Regulated Rivers Using Multi-Source Remote Sensing Observations

The number of reservoirs is rapidly increasing owing to the growth of the world’s economy and related energy and water needs. Yet, for the vast majority of reservoirs around the world, their locations and related information, especially for newly dammed reservoirs, are not readily available due to financial, political, or legal considerations. This study proposes an automated method of identifying newly dammed reservoirs from time series of MODIS-derived NDWI (normalized difference water index) images. Its main idea lies in the detection of abrupt changes in the NDWI time series that are associated with land-to-water conversion due to the reservoir impoundment. The proposed method is tested in the upper reach of the Yellow River that is severely regulated by constructed reservoirs. Our results show that five newly dammed reservoirs were identified in the test area during 2000–2018. Validated against high-resolution Google Earth imagery, our method is effective to determine both locations of the emerging medium-size reservoirs and the timing of their initial water impoundments. Such information then allows for a refined calculation of the reservoir inundation extents and storage capacities through the combination of higher-resolution Landsat imagery and SRTM DEM. The comparison of our estimated reservoir areas and capacities against documented information further indicates that the integration of multi-mission remote sensing data may provide useful information for understanding reservoir operations and impacts on river discharges. Our method also demonstrates a potential for regional or global inventory of emerging reservoirs, which is crucial to assessing human impacts on river systems and the global water cycle.


Introduction
Reservoirs play an important role in controlling and managing water resources and hence facilitating water supply, flood mitigation, agricultural irrigation, hydroelectric power generation, recreation, and other water uses [1][2][3][4].Globally, water from reservoirs supplies an estimated 30-40% of irrigated areas and contributes 20% of global electricity generation in the form of hydropower [5].Up to today, over half of the world's large rivers are regulated by dams [6].Humans have constructed over 1 million dams worldwide [7] and among them 50,000 are large dams (higher than 15 m) [2,8,9].The number of reservoirs is still increasing markedly owing to the growth of the world's economy and related energy needs.Given their broad distribution and high variety, detailed and updated knowledge of the global inventory of reservoirs over time is critical to hydrological studies, such as assessing the impact of reservoirs on river discharge [10][11][12] and improving climate-and water-related risk management [13][14][15].Despite the recognized importance of reservoirs and dams, global datasets describing their characteristics and geographical distribution are incomplete [9,16,17].For the vast majority of reservoirs around the world, their locations and inherent spatial parameters (including extent, storage, etc.) are either not measured or not readily available due to financial, political, or legal considerations [18][19][20][21][22]. Particularly, newly emerged small reservoirs are of increasing concern to researchers because their cumulative storage may be non-negligible, but they remained poorly inventoried or monitored thus far [23].
Remote sensing has become a promising tool to overcome the difficulty in acquiring the reliable information of reservoir construction and operation [18,24], especially in remote and conflict-torn parts of the world where in-situ monitoring networks are poorly available.This study proposes a novel approach combining MODIS, Landsat imagery and digital elevation models (DEMs) to identify newly-dammed reservoirs and to retrieve information on water extent and storage in any inaccessible region.It consists of three key steps: (i) identification of the presence of new reservoirs by using the breaks for additive season and trend (BFAST) algorithm [25,26] to detect the abrupt land-to-water conversion induced by impoundment from high-frequency MODIS observations since 2000, (ii) accurate delineation of identified reservoir extents based on Landsat imagery, (iii) and estimation of reservoir water storages by combining mapped reservoir boundaries and the bathymetry exposed in the SRTM DEM (as acquired in February 2000 before the reservoirs' water impoundment) [27][28][29][30].The approach is applied to identify newly emerged reservoirs in the upper Yellow River basin in Northwest China, where ground monitoring is impeded by economic conditions and information-sharing policy.Results are validated against available high-resolution historical image series in Google Earth.

Study Area
The study area (Figure 1) is located in the upper reach of the Yellow River, in the northeast of the Qinghai-Tibet Plateau known as "the Roof of the World".With considerable altitudinal drop, it is the origin of many rivers, including the Yellow River, which is about 5464 km long and ranks as the sixth longest river system in the world.In the upper reach of the Yellow River, a large number of dams were constructed for hydropower generation, flood control and agricultural irrigation.For example, two giant reservoirs, the Liujiaxia Reservoir and the Longyangxia Reservoir, were built in 1974 and 1992, respectively, with water storage capacities of 5.7 km 3 and 24 km 3 .

Figure 1.
The study area is located in the upper reach of the Yellow River, with a lot of reservoirs built on it.The overview image shows the location of the study area and the terrains in and around the study area.The background was generated by mosaicking two adjacent Landsat-8 OLI RGB composite images.The locations of newly-built reservoirs/dams in the study area derived from Google Earth imagery are marked in this figure.

Datasets and Preprocessing
In this study, we used the MODIS MOD09A1 surface reflectance product [31], with a 500 m spatial resolution, eight-day composite temporal interval, and high spectral and radiometric resolutions.It allows for reconstructing the time-series water index information to indicate the changes of land to open water at the newly-built artificial reservoirs.To inventory a long period of new reservoirs, we collected the MOD09A1 data from February 2000 to April 2018 at the data center ORNL DAAC (https://daac.ornl.gov/).
Landsat-8 OLI satellite imagery with higher spatial resolution were used to extract accurate reservoirs inundation extents [32][33][34].They provide near infrared and green bands to calculate the normalized difference water index (NDWI) values.All available cloud-free Landsat-8 OLI images acquired from 2013 to 2017 were collected in this study.
USGS C-Band SRTM-1 DEM acquired in February 2000 was used as the reservoir bathymetry reference to estimate the water storage of newly impounded reservoirs [35,36].Its accuracy and RMSE are comparable to the SRTM-3 dataset, while its higher spatial resolution (30 m) makes it possible to estimate the reservoir storage capacities of small water bodies.

General Methodology
Remote sensing has long been a useful tool for detecting land cover changes [37].Particularly, the long-term, consecutive satellite observations can provide comprehensive records on various surface features on Earth from different space and time dimensions, including gradual, seasonal, and ephemeral changes.In general, the construction of new dams/reservoirs leads to a long-term alternation from vegetation cover to water bodies, which can be recorded in satellite-detected spectral signals.Thus, we expect to identify the newly dammed reservoirs by examining the abrupt disturbance existing in low frequency signal (gradual changes) of MODIS-based NDWI time series.The BFAST algorithm has been developed to conduct the signal decomposition, which decomposes time series into trend, seasonal, and remainder components while detecting break points in trend and

Datasets and Preprocessing
In this study, we used the MODIS MOD09A1 surface reflectance product [31], with a 500 m spatial resolution, eight-day composite temporal interval, and high spectral and radiometric resolutions.It allows for reconstructing the time-series water index information to indicate the changes of land to open water at the newly-built artificial reservoirs.To inventory a long period of new reservoirs, we collected the MOD09A1 data from February 2000 to April 2018 at the data center ORNL DAAC (https://daac.ornl.gov/).
Landsat-8 OLI satellite imagery with higher spatial resolution were used to extract accurate reservoirs inundation extents [32][33][34].They provide near infrared and green bands to calculate the normalized difference water index (NDWI) values.All available cloud-free Landsat-8 OLI images acquired from 2013 to 2017 were collected in this study.
USGS C-Band SRTM-1 DEM acquired in February 2000 was used as the reservoir bathymetry reference to estimate the water storage of newly impounded reservoirs [35,36].Its accuracy and RMSE are comparable to the SRTM-3 dataset, while its higher spatial resolution (30 m) makes it possible to estimate the reservoir storage capacities of small water bodies.

General Methodology
Remote sensing has long been a useful tool for detecting land cover changes [37].Particularly, the long-term, consecutive satellite observations can provide comprehensive records on various surface features on Earth from different space and time dimensions, including gradual, seasonal, and ephemeral changes.In general, the construction of new dams/reservoirs leads to a long-term alternation from vegetation cover to water bodies, which can be recorded in satellite-detected spectral signals.Thus, we expect to identify the newly dammed reservoirs by examining the abrupt disturbance existing in low frequency signal (gradual changes) of MODIS-based NDWI time series.The BFAST algorithm has been developed to conduct the signal decomposition, which decomposes time series into trend, seasonal, and remainder components while detecting break points in trend and seasonal components [25,26].The break points in the trend component indicate abrupt changes in the NDWI time series.However, some break points are not related to land-to-water conversions caused by damming activities (thereafter referred to as "false break points" and the contrary referred to as "true break points"), so a scheme of combined threshold setting (see Section 3.2.3)was introduced to filter out the false signals (Figure 2).NDWI time series.However, some break points are not related to land-to-water conversions caused by damming activities (thereafter referred to as "false break points" and the contrary referred to as "true break points"), so a scheme of combined threshold setting (see Section 3.2.3)was introduced to filter out the false signals (Figure 2).
To accurately delineate the inundation extents of the identified new reservoirs, we further derived the NDWI imagery from multi-temporal Landsat-8 OLI images to map the water inundation areas of these reservoirs after impoundments.For each reservoir, the water storage was estimated based on the volume difference between the reference bathymetry as depicted in the SRTM DEM and the level of the reservoir boundary mapped from Landsat images.Considering that the relatively steep slopes along the river valley may induce some uncertainties on reconstructing elevations of reservoir waterline, we assume one half Landsat-pixel shift (about ±15 m) as the error (land-water mixed pixel) of reservoir boundary delineation to estimate the uncertainties of reservoir area and volume estimation by referring to the prior commonly used method [38,39].

Identifying Locations and Time of New Reservoir Emergence
The identification method includes determination of the existence of newly dammed reservoirs and the associated impoundment date by analyzing the abrupt changes of NDWI time series of each pixel (as shown in Figure 2).The effectiveness and performance of BFAST in detecting such changes has been widely validated in the field of vegetation monitoring [25,26,40,41].Several studies also applied the tool to detect abnormal changes of water storage and monitor long-term trend changes To accurately delineate the inundation extents of the identified new reservoirs, we further derived the NDWI imagery from multi-temporal Landsat-8 OLI images to map the water inundation areas of these reservoirs after impoundments.For each reservoir, the water storage was estimated based on the volume difference between the reference bathymetry as depicted in the SRTM DEM and the level of the reservoir boundary mapped from Landsat images.Considering that the relatively steep slopes along the river valley may induce some uncertainties on reconstructing elevations of reservoir waterline, we assume one half Landsat-pixel shift (about ±15 m) as the error (land-water mixed pixel) of reservoir boundary delineation to estimate the uncertainties of reservoir area and volume estimation by referring to the prior commonly used method [38,39].

Identifying Locations and Time of New Reservoir Emergence
The identification method includes determination of the existence of newly dammed reservoirs and the associated impoundment date by analyzing the abrupt changes of NDWI time series of each pixel (as shown in Figure 2).The effectiveness and performance of BFAST in detecting such changes has been widely validated in the field of vegetation monitoring [25,26,40,41].Several studies also applied the tool to detect abnormal changes of water storage and monitor long-term trend changes [42][43][44].NDWI is an effective index to differentiate open water features and other land cover types.Here we used NDWI which normalizes the difference between the green (ρ green ) and near infrared (ρ NIR ) spectral bands as proposed by McFeeters [45].A few pixels in the MOD09A1 satellite images may be in poor quality or missing spectral values, resulting in gaps in the NDWI time series.These missing values were filled by averaging the temporally consecutive values in the NDWI time series.
When conducting the BFAST decomposition algorithm, the effectiveness largely depends on parameter settings.The parameter 'h' represents the minimum interval of the break points divided by the length of the time series.If it was set to 0.1, the minimum time interval of NDWI abrupt change between two breaking points is about two years.This setting can effectively exclude signals of abrupt changes that occurred within one-year cycle of the 18-year time series (see Section 5).The season model was set to 'harmonic', which has proven to be more suitable and robust for detecting abrupt and harmonic changes in time series [26].The value of maximum iteration times was set to 3, which ensures a stable and robust decomposition result.The later algorithm tests suggest that the decomposition iteration times are mostly one or two times and only a small proportion need three times.The number of breaks was set to 1 in order to reduce the possibility of detecting false break points because other small disturbances (e.g., interannual variability of land cover) may also result in the detection of breaking points.By setting the number of break points to one we assume that the newly built reservoir is not removed thereafter.The confidence level was set to 0.05 which gave a satisfactory confidence level of 95%.
Damming activities lead to the conversion from land to water in the reservoir area, which means a sudden rise in the NDWI time series and the BFAST-decomposed trend component after the break point.This phenomenon was utilized to distinguish land-to-water conversion from other kinds of conversion such as from vegetation to bare land.To eliminate such false conversions, we applied a threshold on the mean NDWI of the two segments, namely the mean NDWI value in the trend component before the break point (referred to as T f ormer ) and the mean NDWI value after the break point (T latter ).If T f ormer was smaller than the threshold and T latter was larger than it, this indicates abrupt NDWI increase from the former to the latter periods.If so, this break point was classified as 'true'; otherwise, it was classified as 'false'.To achieve a proper threshold value, a two-step scheme was applied to the study area.In the first step, all detected break points were checked by the aforementioned threshold method.The second step was to delineate the unfiltered pixels in a result raster layer and assess the accuracy of filtering by overlapping the raster layer on cloud-free Landsat-8 images acquired in 2017.Then the first and second steps were repeated with different threshold values until the delineating and filtering results were satisfactory.
Finally, the threshold value was set to around −0.2 based on our tests on the threshold sensitivity, which will be discussed in Section 5.After filtering with it, pixels with break points caused by damming activities in MODIS imagery were delineated in a result raster and the date of the breakpoints was assigned to the raster value.By sampling the plots of original data, decomposed components and breaks generated by BFAST algorithm, the pixels in the study area could be generally classified into three major scenarios.The ID of each pixel was set as its index in MODIS imagery in row by column order.

Scenario One: Land-to-Water Areas with Damming Activities
Damming activities lead to land-to-water transformations and raises of NDWI trend component at the break point, and the values in trend component after the break point would become much higher.It had been proved by the plots of BFAST time series decomposition and break points detection (referred to as "BFAST result plots") sampled in areas of newly dammed reservoirs.As shown in Figure 3, the simultaneously abrupt changes in original NDWI time series and trend component can be detected as land-to-water pixels.The IDs of the samples in Figure 3

Scenario Two: Water Areas without Damming Activities
The water presence in river and lake areas stays relatively stable in contrast to land-to-water areas.Break points may be detected in trend or seasonal components in those areas due to interannual fluctuations, but the values in trend components before the break points in those areas are higher than those in land-to-water areas and close to their counterparts.The IDs of the sampled plots in Figure 4

Scenario Two: Water Areas without Damming Activities
The water presence in river and lake areas stays relatively stable in contrast to land-to-water areas.Break points may be detected in trend or seasonal components in those areas due to interannual fluctuations, but the values in trend components before the break points in those areas are higher than those in land-to-water areas and close to their counterparts.The IDs of the sampled plots in

Scenario Two: Water Areas without Damming Activities
The water presence in river and lake areas stays relatively stable in contrast to land-to-water areas.Break points may be detected in trend or seasonal components in those areas due to interannual fluctuations, but the values in trend components before the break points in those areas are higher than those in land-to-water areas and close to their counterparts.The IDs of the sampled plots in Figure 4

Extracting Reservoir Extents and Estimating Impounded Water Volumes
Reservoir water extents usually vary over time due to their filling and discharging operations.In this study, we mapped multi-temporal reservoir extents from time-series satellite images and derived the median-sized water area to indicate the stable stage of reservoirs.All cloud-free Landsat-8 OLI images acquired from 2013 to 2017 were collected to map reservoir inundation areas.Six to ten Landsat 8 observations were used for each of the identified newly-dammed reservoirs.
The variations of spectral features of different water bodies on different images pose challenges for automated retrieval of reservoir extents.To address the challenges, a two-step NDWI threshold segmentation scheme was applied in this study to extract the water inundation areas [34,[46][47][48][49].In the first step, all potential waterbody pixels are filtered by a loose initial threshold segmentation (NDWI > −0.1) based on calculated NDWI maps of selected images; the second step is locally iterative segmentation for each potential lake entity, and the iteration goal is to obtain stable lake boundary within consecutive segmentations with the least shifted water shoreline.The two-step NDWI threshold segmentation method can efficiently improve the accuracy of waterbody boundary delineation, compared to single-threshold method [47,50,51].Eventually, we mapped the inundation areas from multi-dates Landsat images for each reservoir within manually defined reservoir extents (the evidently widened river segments compared to the upper and downstream reaches).The median-sized extent was chosen to calculate the reservoir area and water storage.
The estimation of impounded water volume is based on the median-sized reservoir extents.By calculating the filling mass above the DEM-based bathymetry reference and under the average altitude of reservoir shorelines, the water storage volume of each identified newly dammed reservoir can be estimated.The surface volume tool in ArcToolbox provided by ESRI was used conduct the estimation and the input DEM data was from the USGS SRTM-1 DEM.

Abrupt Changes in NDWI Time Series Due to Damming Activities
According to the methodology in Section 3.2.2, a total of 14,300 pixels were detected with break points in the study area but only 291 pixels were identified as the potential target reservoir areas based on the procedure of threshold-based exclusion of false break points.By filtering scattered pixels (clustered pixel number < 3), the remaining pixel clusters are identified as five newly dammed reservoirs that are also confirmed in historical high-resolution image series in Google Earth, as shown in Figure 1.By manual inspection in Google Earth, there is a small new reservoir, Dahejia Reservoir, omitted in our automated method.The timings of most of identified pixels, which indicate the timings of land-to-water transformations that happened in these pixels, agree well with factual damming activities indicated by the Google Earth images.

Extracting Reservoir Extents and Estimating Impounded Water Volumes
Reservoir water extents usually vary over time due to their filling and discharging operations.In this study, we mapped multi-temporal reservoir extents from time-series satellite images and derived the median-sized water area to indicate the stable stage of reservoirs.All cloud-free Landsat-8 OLI images acquired from 2013 to 2017 were collected to map reservoir inundation areas.Six to ten Landsat 8 observations were used for each of the identified newly-dammed reservoirs.
The variations of spectral features of different water bodies on different images pose challenges for automated retrieval of reservoir extents.To address the challenges, a two-step NDWI threshold segmentation scheme was applied in this study to extract the water inundation areas [34,[46][47][48][49].In the first step, all potential waterbody pixels are filtered by a loose initial threshold segmentation (NDWI > −0.1) based on calculated NDWI maps of selected images; the second step is locally iterative segmentation for each potential lake entity, and the iteration goal is to obtain stable lake boundary within consecutive segmentations with the least shifted water shoreline.The two-step NDWI threshold segmentation method can efficiently improve the accuracy of waterbody boundary delineation, compared to single-threshold method [47,50,51].Eventually, we mapped the inundation areas from multi-dates Landsat images for each reservoir within manually defined reservoir extents (the evidently widened river segments compared to the upper and downstream reaches).The median-sized extent was chosen to calculate the reservoir area and water storage.
The estimation of impounded water volume is based on the median-sized reservoir extents.By calculating the filling mass above the DEM-based bathymetry reference and under the average altitude of reservoir shorelines, the water storage volume of each identified newly dammed reservoir can be estimated.The surface volume tool in ArcToolbox provided by ESRI was used conduct the estimation and the input DEM data was from the USGS SRTM-1 DEM.

Abrupt Changes in NDWI Time Series Due to Damming Activities
According to the methodology in Section 3.2.2, a total of 14,300 pixels were detected with break points in the study area but only 291 pixels were identified as the potential target reservoir areas based on the procedure of threshold-based exclusion of false break points.By filtering scattered pixels (clustered pixel number < 3), the remaining pixel clusters are identified as five newly dammed reservoirs that are also confirmed in historical high-resolution image series in Google Earth, as shown in Figure 1.By manual inspection in Google Earth, there is a small new reservoir, Dahejia Reservoir, omitted in our automated method.The timings of most of identified pixels, which indicate the timings of land-to-water transformations that happened in these pixels, agree well with factual damming activities indicated by the Google Earth images.

Reservoir Inundation Areas and Water Storages
Based on multi-temporal Landsat-8 images, we extracted median-sized water extent for each newly dammed reservoir, the inundation areas of the five identified reservoirs were superimposed on Landsat background images as shown in Figure 7.The areas of the five newly-dammed reservoirs vary largely.The largest one is Gongboxia Reservoir located in the upper reach, with an area of 21.7 ± 2.6 km 2 .The other four reservoirs range from 2.27 to 6.65 km 2 .
Remote Sens. 2018, 10, x FOR PEER REVIEW 10 of 21 Reservoir (e).They showed the high consistency of locations of pixels detected as true break points and the timing of true break points which correspond to damming dates indicated by Google Earth imagery and documentary materials.The lines in the figures point at the location of dams.

Reservoir Inundation Areas and Water Storages
Based on multi-temporal Landsat-8 images, we extracted median-sized water extent for each newly dammed reservoir, the inundation areas of the five identified reservoirs were superimposed on Landsat background images as shown in Figure 7.The areas of the five newly-dammed reservoirs vary largely.The largest one is Gongboxia Reservoir located in the upper reach, with an area of 21.7 ± 2.6 km 2 .The other four reservoirs range from 2.27 to 6.65 km 2 .Using the median-sized reservoir extents over the SRTM-DEM based reference bathymetry, the estimated water volumes reflect the multi-year average water storages of these newly dammed reservoirs.Documented water storage volumes at the end of 2015 and total reservoir capacities are introduced to compare with remote-sensing-based water impoundment estimation (refer to Table  Using the median-sized reservoir extents over the SRTM-DEM based reference bathymetry, the estimated water volumes reflect the multi-year average water storages of these newly dammed reservoirs.Documented water storage volumes at the end of 2015 and total reservoir capacities are introduced to compare with remote-sensing-based water impoundment estimation (refer to Table 1 note).As shown in Table 1, the estimated reservoir impoundment volumes based on median-sized inundation areas are lower than both the reported reservoir capacity and the documented water storage in the end of 2015.Uncertainties of estimated reservoir volume were reported in Table 1 as described in Methodology.In comparison, the remote sensing-based volume estimations for the three reservoirs Gongboxia Reservoir ((3.99 ± 1.09) × 10 8 m 3 ), Suzhi Reservoir ((3.85 ± 1.00) × 10 7 m 3 ), and Sigouxia Reservoir ((3.33 ± 0.62) × 10 7 m 3 ), are relatively close to the documented reservoir capacities, with a bias of 15-35% approximately.For the other two, the estimated water storages only account for about 20% of reservoir capacity.Several factors might jointly attribute to their differences.First, the actual reservoir water storages with regular release and recharge operations could hardly reach the engineering designed reservoir capacities which indicate the possible maximum water storage volumes.Secondly, water release is usually restricted at the end of each year in case of coming droughts in the spring, while our estimates based on median-sized inundation extent only represent the average status of reservoir water storages.Thirdly, the two reservoirs-Huangfeng Reservoir and Jishixia Reservoir-were dammed in recent years and thus impounded water in small proportions of the reservoir capacity in most cases.Additionally, the remote sensing estimations were based on 30-m resolution DEM and thus cannot avoid the uncertainties for qualifying the water volume for some reservoirs with steep river reaches.

Discussions
The presented methodology, based on BFAST algorithm, provided a cost-effective way to detect newly dammed reservoirs from time series satellite images and their impoundment time.It utilizes massive amounts of free-access remote sensing data with nearly two-decade temporal coverage to detect abrupt trend changes of NDWI time series that are associated with water impoundment of dammed reservoirs during the early 21st century.We also compared the timing of damming activities indicated by Google Earth imageries and true break points indicated that the BFAST-based identification method.The results suggest that the proposed method is able to effectively detect abrupt changes in response to water impoundment.Most false break points can be eliminated with few of them remaining unfiltered due to the inherent drawbacks of threshold-based filtering method, but this does not necessarily affect the judgment of new reservoirs as shown in Figure 8. From Figure 8a-c, it is easy to find that the false break point pixels are mostly far away from the river channel and dispersedly distributed.It could be eliminated by applying the method to a buffer of river channel and isolating scattered pixels by image smoothing.
MODIS quality control flags were not used in this study.In the MOD09A1 product each pixel contains the best possible L2G observation during an 8-day period as selected on the basis of high observation coverage, low view angle, absence of clouds or cloud shadow, and aerosol loading [31], which greatly reduced the number of low quality pixels.From the BFAST results, it can be seen that the remaining proportion of extremely high and low values were decomposed into remainder components and thus does not largely affect the detection of abrupt changes in trend components.
The spatial resolution of MODIS imagery is one primary restriction factor obstructing the accurate identification of newly dammed reservoirs.Due to the relatively low spatial resolution, the effectiveness of BFAST-based time series decomposition and break points detection may be influenced due to the mixed pixels between the land/water interfaces.For example, the pixels in Dahejia Reservoir are all mixed pixels of land and water features (Figure 9), which obstructed the successful identification from time series decomposition and break points detection, rendering no true break points detected among them.
Remote Sens. 2018, 10, x FOR PEER REVIEW 12 of 21 but this does not necessarily affect the judgment of new reservoirs as shown in Figure 8. From Figure 8a-c, it is easy to find that the false break point pixels are mostly far away from the river channel and dispersedly distributed.It could be eliminated by applying the method to a buffer of river channel and isolating scattered pixels by image smoothing.MODIS quality control flags were not used in this study.In the MOD09A1 product each pixel contains the best possible L2G observation during an 8-day period as selected on the basis of high observation coverage, low view angle, absence of clouds or cloud shadow, and aerosol loading [31], which greatly reduced the number of low quality pixels.From the BFAST results, it can be seen that the remaining proportion of extremely high and low values were decomposed into remainder components and thus does not largely affect the detection of abrupt changes in trend components.The spatial resolution of MODIS imagery is one primary restriction factor obstructing the accurate identification of newly dammed reservoirs.Due to the relatively low spatial resolution, the effectiveness of BFAST-based time series decomposition and break points detection may be influenced due to the mixed pixels between the land/water interfaces.For example, the pixels in Dahejia Reservoir are all mixed pixels of land and water features (Figure 9), which obstructed the successful identification from time series decomposition and break points detection, rendering no true break points detected among them.The spatial resolution of MODIS imagery is one primary restriction factor obstructing the accurate identification of newly dammed reservoirs.Due to the relatively low spatial resolution, the effectiveness of BFAST-based time series decomposition and break points detection may be influenced due to the mixed pixels between the land/water interfaces.For example, the pixels in Dahejia Reservoir are all mixed pixels of land and water features (Figure 9), which obstructed the successful identification from time series decomposition and break points detection, rendering no true break points detected among them.
(c) When conducting the BFAST decomposition algorithm, the h parameter that indicates the sensitivity of the period of detected abrupt changes to the whole study period will influence the number and timing of break points (Figure 10a).We randomly sampled five pixels for each identified new reservoir and assessed the performance of BFAST algorithm affected by different settings of h parameter from 1/3 to 1/17.Using h = 1/5, where the minimum length of each trend segment is roughly 3.6 years, abrupt changes were not detected for Gongboxia Reservoir at the time of the recorded damming activities (between 2003 and 2005).Similarly, using h = 1/4, where the minimum length of each trend segment is approximately 4.5 years, abrupt changes cannot be detected effectively in Huangfeng Reservoir at the time of the recorded damming activities (around the year 2014).It suggests that the minimum intervals of breaks were too wide to accurately reflect damming activities over the period of satellite image time series.However, as the value of h decreases, the accuracy of breaks began to increase and became acceptable when it was set to 1/7, and then it would be generally stable till 1/17.The h parameter also has a great impact on the averages of trend components (refer to as "ATC" in abbreviation) before and after break points, which may greatly influence the effectiveness of the afore-mentioned threshold-based pixels with false break points elimination.We sampled 10 pixels in water area and 10 pixels in land area, and ATCs of the sampled 25 reservoir area pixels were also used.Note that most of the 20 pixels were sampled near the river, considering buffer would be introduced to eliminate unfiltered pixels with false break points.Among the land-to-water pixels, when h was set to 1/3, their ATCs before and after break points may be far away from stable values, but when it was set to 1/7, ATCs would be generally stable and when it was set to 1/10 or smaller, every single pixel's ATCs would equal each other, which was shown in Figure 10b,c.ATCs of water area pixels after break points and ATCs of land area pixels before and after break points are generally stable, but ATCs of water area pixels before break points would suddenly drop when h was changed from h = 1/10 to h = 1/11.From the ATC plots of all kinds of pixels in Figures 10b and 11, all of the sampled pixels' ATCs derived from different settings of h parameter varying from 1/3 to 1/17, except three pixels' ATCs when h parameter was set to 1/3, perfectly suit the −0.2 threshold and would be classified into classes they belong to.This indicates that the effectiveness of the threshold which equals to −0.2 is robust enough to be less sensitive to different settings of h parameter.So a likely effective approach to increase the accuracy of breaks detected in newly dammed reservoirs whose timing of damming activities are close to the end or beginning of time series period is to lower the value of h parameter such as one or two-year period (h ~0.10-0.05),and it would not greatly influence the accuracy of newly-dammed reservoir identification.When conducting the BFAST decomposition algorithm, the h parameter that indicates the sensitivity of the period of detected abrupt changes to the whole study period will influence the number and timing of break points (Figure 10a).We randomly sampled five pixels for each identified new reservoir and assessed the performance of BFAST algorithm affected by different settings of h parameter from 1/3 to 1/17.Using h = 1/5, where the minimum length of each trend segment is roughly 3.6 years, abrupt changes were not detected for Gongboxia Reservoir at the time of the recorded damming activities (between 2003 and 2005).Similarly, using h = 1/4, where the minimum length of each trend segment is approximately 4.5 years, abrupt changes cannot be detected effectively in Huangfeng Reservoir at the time of the recorded damming activities (around the year 2014).It suggests that the minimum intervals of breaks were too wide to accurately reflect damming activities over the period of satellite image time series.However, as the value of h decreases, the accuracy of breaks began to increase and became acceptable when it was set to 1/7, and then it would be generally stable till 1/17.The h parameter also has a great impact on the averages of trend components (refer to as "ATC" in abbreviation) before and after break points, which may greatly influence the effectiveness of the afore-mentioned threshold-based pixels with false break points elimination.We sampled 10 pixels in water area and 10 pixels in land area, and ATCs of the sampled 25 reservoir area pixels were also used.Note that most of the 20 pixels were sampled near the river, considering buffer would be introduced to eliminate unfiltered pixels with false break points.Among the land-to-water pixels, when h was set to 1/3, their ATCs before and after break points may be far away from stable values, but when it was set to 1/7, ATCs would be generally stable and when it was set to 1/10 or smaller, every single pixel's ATCs would equal each other, which was shown in Figure 10b and 10c.ATCs of water area pixels after break points and ATCs of land area pixels before and after break points are generally stable, but ATCs of water area pixels before break points would suddenly drop when h was changed from h = 1/10 to h = 1/11.From the ATC plots of all kinds of pixels in Figures 10b and 11, all of the sampled pixels' ATCs derived from different settings of h parameter varying from 1/3 to 1/17, except three pixels' ATCs when h parameter was set to 1/3, perfectly suit the −0.2 threshold and would be classified into classes they belong to.This indicates that the effectiveness of the threshold which equals to −0.2 is robust enough to be less sensitive to different settings of h parameter.So a likely effective approach to increase the accuracy of breaks detected in newly dammed reservoirs whose timing of damming activities are close to the end or beginning of time series period is to lower the value of h parameter such as one or two-year period (h ~ 0.10-0.05),and it would not greatly influence the accuracy of newly-dammed reservoir identification.
(a)  Break points were detected in 14,300 pixels in the study area by BFAST algorithm, but the threshold-based filtering scheme successfully filtered most of the false break points.The threshold values also influence the effectiveness of false break point eliminating.We have tested many different threshold values varying from 0 to −0.3, and results illustrated in Figure 12, with values of −0.10, −0.15, −0.20, −0.25, and −0.30, show the overall trend of influences caused by different values of threshold setting.Higher values of the threshold generated rigorous filtering of false break points, but more pixels of true break points might be mistakenly filtered.By lowering the values of threshold, pixels of true breaks would be safely kept, but more pixels with unfiltered false break points would also remain.When the value of threshold was set to values between −0.10 and −0.15, experiments showed that it is too high to reserve true break points as shown in Figure 12a,b.On the contrary, when the threshold was set between −0.25 and −0.30, terrestrial and vegetation features would not be effectively eliminated as shown in Figure 12d,e.However, threshold around −0.2 effectively reserve most true break points, and the remaining false points could be eliminated by applying stream buffer zones and image smoothing.The sensitivity analysis of the parameter h has also proved that −0.2 was less sensitive to different settings of h parameter in BFAST algorithm, and that threshold value was generated by applying the aforementioned iteration-based threshold generation scheme to the Gongboxia and Suzhi reservoirs, then by applying the generated value to the full study area, its performance was proven to be effective.Break points were detected in 14,300 pixels in the study area by BFAST algorithm, but the threshold-based filtering scheme successfully filtered most of the false break points.The threshold values also influence the effectiveness of false break point eliminating.We have tested many different threshold values varying from 0 to −0.3, and results illustrated in Figure 12, with values of −0.10, −0.15, −0.20, −0.25, and −0.30, show the overall trend of influences caused by different values of threshold setting.Higher values of the threshold generated rigorous filtering of false break points, but more pixels of true break points might be mistakenly filtered.By lowering the values of threshold, pixels of true breaks would be safely kept, but more pixels with unfiltered false break points would also remain.When the value of threshold was set to values between −0.10 and −0.15, experiments showed that it is too high to reserve true break points as shown in Figure 12a,b.On the contrary, when the threshold was set between −0.25 and −0.30, terrestrial and vegetation features would not be effectively eliminated as shown in Figure 12d,e.However, threshold around −0.2 effectively reserve most true break points, and the remaining false points could be eliminated by applying stream buffer zones and image smoothing.The sensitivity analysis of the parameter h has also proved that −0.2 was less sensitive to different settings of h parameter in BFAST algorithm, and that threshold value was generated by applying the aforementioned iteration-based threshold generation scheme to the Gongboxia and Suzhi reservoirs, then by applying the generated value to the full study area, its performance was proven to be effective.

Summary
Knowledge of reservoir inventory over space and time is crucial to assessing the impact of reservoirs, while their updated information and spatial parameters (including location, inundation extent, and water storage) are usually not measured or not readily available due to financial, political, or legal considerations.This study proposes an automated method of BFAST-based identification of newly dammed reservoirs from time series of MODIS-derived NDWI images.The method is tested in the upper reach of the Yellow River which is severely regulated by artificial reservoirs.The results show that five reservoirs newly dammed during 2000-2018 were identified successfully in the study

Summary
Knowledge of reservoir inventory over space and time is crucial to assessing the impact of reservoirs, while their updated information and spatial parameters (including location, inundation extent, and water storage) are usually not measured or not readily available due to financial, political, or legal considerations.This study proposes an automated method of BFAST-based identification of newly dammed reservoirs from time series of MODIS-derived NDWI images.The method is tested in the upper reach of the Yellow River which is severely regulated by artificial reservoirs.The results show that five reservoirs newly dammed during 2000-2018 were identified successfully in the study area.The proposed method is able to effectively derive the reservoir locations and water impoundment time by validation against high-resolution imagery of Google Earth.Only a very small new reservoir is omitted in the detection due to the narrow reservoir reach less than one MODIS pixel.The combination of Landsat imagery and SRTM DEM was used to further estimate the median-sized reservoir inundation areas and water storages in recent years.The comparison of remote sensing-based reservoir area, impoundment time and volume estimations and documented information indicates that remote sensing may provide more practical knowledge of reservoir operations and their impacts on river discharges.Furthermore, the key points including parameters of the BFAST algorithm and threshold setting for excluding false targets are discussed.Another advantage of the BFAST-based damming identification developed in this study is the potential capability of detecting both gradual changes from trend component and phenological changes due to seasonal variation in reservoir water storage.These detected changes provide useful reservoir operation information such as water impoundment and release or adjustment of reservoir regulation policy, and thus assist the assessment of catchment-level hydrological impacts from the construction and operation of artificial reservoirs.

Figure 1 .
Figure 1.The study area is located in the upper reach of the Yellow River, with a lot of reservoirs built on it.The overview image shows the location of the study area and the terrains in and around the study area.The background was generated by mosaicking two adjacent Landsat-8 OLI RGB composite images.The locations of newly-built reservoirs/dams in the study area derived from Google Earth imagery are marked in this figure.

Figure 2 .
Figure 2. Schematic flowchart of the methodology, including identification of newly dammed reservoirs, extraction of high spatial-resolution reservoirs and estimation of reservoir water volume.

Figure 2 .
Figure 2. Schematic flowchart of the methodology, including identification of newly dammed reservoirs, extraction of high spatial-resolution reservoirs and estimation of reservoir water volume.

21 Figure 3 .
Figure 3. Plots of original NDWI time series and BFAST-decomposed seasonal, trend and remainder components and BFAST-detected break points of samples in land-to-water areas.In each plot, the four panels from the top to the bottom are original NDWI time series ("Yt"), seasonal component ("St"), trend component ("Tt") and reminder component ("et"), respectively.The x-axis shows time.The gray vertical dashed lines in the rows of trend components indicate timing of break points within the confidence interval of 95%.

Figure 4 .
Figure 4. Plots of NDWI time series and BFAST-decomposed seasonal, trend, and remainder components and BFAST-detected break points of samples in water areas without damming activities in the period of NDWI time series.These sampled plots showed that the values in trend component before the break points (if detected) are higher than those in land-to-water pixels.

Figure 3 .
Figure 3. Plots of original NDWI time series and BFAST-decomposed seasonal, trend and remainder components and BFAST-detected break points of samples in land-to-water areas.In each plot, the four panels from the top to the bottom are original NDWI time series ("Yt"), seasonal component ("St"), trend component ("Tt") and reminder component ("et"), respectively.The x-axis shows time.The gray vertical dashed lines in the rows of trend components indicate timing of break points within the confidence interval of 95%.

21 Figure 3 .
Figure 3. Plots of original NDWI time series and BFAST-decomposed seasonal, trend and remainder components and BFAST-detected break points of samples in land-to-water areas.In each plot, the four panels from the top to the bottom are original NDWI time series ("Yt"), seasonal component ("St"), trend component ("Tt") and reminder component ("et"), respectively.The x-axis shows time.The gray vertical dashed lines in the rows of trend components indicate timing of break points within the confidence interval of 95%.

Figure 4 .
Figure 4. Plots of NDWI time series and BFAST-decomposed seasonal, trend, and remainder components and BFAST-detected break points of samples in water areas without damming activities in the period of NDWI time series.These sampled plots showed that the values in trend component before the break points (if detected) are higher than those in land-to-water pixels.

Figure 4 .
Figure 4. Plots of NDWI time series and BFAST-decomposed seasonal, trend, and remainder components and BFAST-detected break points of samples in water areas without damming activities in the period of NDWI time series.These sampled plots showed that the values in trend component before the break points (if detected) are higher than those in land-to-water pixels.

Figure 5 .
Figure 5. Plots of NDWI time series and BFAST-decomposed seasonal, trend, and remainder components and BFAST-detected break points of samples in land areas.These sampled plots showed that the values in trend component after the break points (if detected) are lower than those in landto-water pixels.

Figure 5 .
Figure 5. Plots of NDWI time series and BFAST-decomposed seasonal, trend, and remainder components and BFAST-detected break points of samples in land areas.These sampled plots showed that the values in trend component after the break points (if detected) are lower than those in land-to-water pixels.

As shown in Figure 6 ,
high resolution imagery in Google Earth was used to qualitatively evaluate the performances of our extracted locations, extents, and timings of each newly dammed reservoir.Two adjacent scenes of Landsat-8 OLI RGB composite images acquired on 31 July 2017 and 7 August 2017, respectively, were used in the background.Other details about the identified reservoirs were derived from the Yellow River Conservancy Commission of the Ministry of Water Resources official website and government documentations (http://www.qhepb.gov.cn/hjgl/qyzc/,http://www.yellowriver.gov.cn/).Along the flow direction of the Yellow River, there are six reservoirs dammed after 2000 in the study area, which are Gongboxia Reservoir (impounded in 2004), Suzhi Reservoir (2006), Huangfeng Reservoir (2014), Jishixia Reservoir (2011), Dahejia Reservoir (2015, omitted in our detection), and Sigouxia Reservoir (2009).The Liujiaxia Reservoir is also located in the study area, but it was constructed in 1974.Remote Sens. 2018, 10, x FOR PEER REVIEW 8 of 21 As shown in Figure 6, high resolution imagery in Google Earth was used to qualitatively evaluate the performances of our extracted locations, extents, and timings of each newly dammed reservoir.Two adjacent scenes of Landsat-8 OLI RGB composite images acquired on 31 July 2017 and 7 August 2017, respectively, were used in the background.Other details about the identified reservoirs were derived from the Yellow River Conservancy Commission of the Ministry of Water Resources official website and government documentations (http://www.qhepb.gov.cn/hjgl/qyzc/,http://www.yellowriver.gov.cn/).Along the flow direction of the Yellow River, there are six reservoirs dammed after 2000 in the study area, which are Gongboxia Reservoir (impounded in 2004), Suzhi Reservoir (2006), Huangfeng Reservoir (2014), Jishixia Reservoir (2011), Dahejia Reservoir (2015, omitted in our detection), and Sigouxia Reservoir (2009).The Liujiaxia Reservoir is also located in the study area, but it was constructed in 1974.

Figure 8 .
Figure 8. Unfiltered false break points detected in terrestrial areas away from the Yellow River (a), near Liujiaxia reservoir (b), and in the Yellow River (c).From the plots of trend component, what could be discovered is that the values in trend components jumped to a little higher value than the threshold we set.With a little higher value of threshold, these false break points would be effectively eliminated, which is the inherent disadvantage of threshold-based elimination.

Figure 8 .
Figure 8. Unfiltered false break points detected in terrestrial areas away from the Yellow River (a), near Liujiaxia reservoir (b), and in the Yellow River (c).From the plots of trend component, what could be discovered is that the values in trend components jumped to a little higher value than the threshold we set.With a little higher value of threshold, these false break points would be effectively eliminated, which is the inherent disadvantage of threshold-based elimination.

Figure 8 .
Figure 8. Unfiltered false break points detected in terrestrial areas away from the Yellow River (a), near Liujiaxia reservoir (b), and in the Yellow River (c).From the plots of trend component, what could be discovered is that the values in trend components jumped to a little higher value than the threshold we set.With a little higher value of threshold, these false break points would be effectively eliminated, which is the inherent disadvantage of threshold-based elimination.

Figure 9 .
Figure 9. Landsat-8 OLI RGB synthesis image of Dahejia Reservoir with the MOD09A1 pixels in its range and the BFAST decomposition result of the pixel at the dam of it.The spatial extents of each MOD09A1 pixel was transparently shown as a chess board to illustrate that all of the pixels in Dahejia Reservoir are combined of different features, which obstructed the effectiveness of BFAST algorithm.

Figure 10 .Figure 10 .Figure 11 .
Figure 10.Timing of break points detected in the time series (a), averages of trend components before break points (b) and averages of trend components after break points (c), influenced by different settings of h parameter sampled in reservoir area.Values in the x-axis indicate different settings of h parameter, values in the y-axis (a) indicate dates of break points detected by BFAST algorithm while values in y-axis (b and c) indicate averages of trend components decomposed by BFAST algorithm.Each line drawn in each diagram indicates their corresponding values of y-axis influenced by different settings of h parameter.They were drawn in different colors to distinguish different sample pixels.The numbers behind each line in legend indicate the IDs of sample pixels.All of the three diagrams share the same legend to diagram (a).

21 Figure 11 .
Figure 11.Averages of trend components of water area pixels before (left) and after (right) break their point, influenced by different settings of h parameter sampled in water (a) and land (b) area.Values in x-axis indicate different settings of h parameter, while values in y-axis indicate averages of trend components before or after break points.Each line drawn in the diagram indicates average of trend components in each sample pixel influenced by different settings of h parameter.They were drawn in different colors to distinguish different sample pixels.The numbers behind each line in the legend indicate the IDs of sample pixels.The no-data points in the diagrams indicate no break points were detected.

Table 1 .
Total reservoir capacity, water storage volume (V) at the end of 2015, estimated inundation areas, estimated water storage V, and ratios of estimated water V to the total reservoir capacity and the water V in 2015.
Landsat-8 OLI RGB synthesis image of Dahejia Reservoir with the MOD09A1 pixels in its range and the BFAST decomposition result of the pixel at the dam of it.The spatial extents of each MOD09A1 pixel was transparently shown as a chess board to illustrate that all of the pixels in Dahejia Reservoir are combined of different features, which obstructed the effectiveness of BFAST algorithm.