Compositing the Minimum NDVI for Daily Water Surface Mapping

: Capturing high frequency water surface dynamics via optical remote sensing is important for understanding hydro-ecological processes over seasonally ﬂooded wetlands. However, it is a di ﬃ cult task due to the presence of clouds on satellite images. This study proposed the MODerate-resolution Imaging Spectroradiometer (MODIS) Normalized Di ﬀ erence Vegetation Index (NDVI) Minimum Value Composite (MinVC) algorithm to generate daily water surface data at a 250-m resolution. The algorithm selected pixelwise minimum values from the combined daily Terra and Aqua MODIS NDVI data within a 15-day moving window. Consisting mainly of cloud and water surface information, the MinVC NDVI data were segmented for water surfaces over the Poyang Lake, China (2000–2017) by using an edge detection model. The water surface mapping result was strongly correlated with the Landsat based result ( R 2 = 0.914, root mean square error, RMSE = 223.7 km 2 ), the cloud free MODIS image based result ( R 2 = 0.824, RMSE = 356.7 km 2 ), the recent Landsat-MODIS image fusion based result ( R 2 = 0.765, RMSE = 403 km 2 ), and the hydrodynamic modeling result ( R 2 = 0.799). Compared to the equivalent eight-day MOD13 NDVI based on the Constraint View-Angle Maximum Value Composite (CV-MVC) algorithm, the daily MinVC NDVI highlighted water bodies by generating spatially homogenous water surface information. Consequently, the algorithm provided spatially and temporally continuous data for calculating water submersion times and trends in water surface area, which contribute to a better understanding of hydro-ecological processes over seasonally ﬂooded wetlands. Within the framework of sensor intercalibration, the algorithm can be extended to incorporate multiple sensor data for improved water surface mapping.


Introduction
A water body is identifiable from satellite images due to its unique spectral features in a wide range of electromagnetic spectrums [1][2][3][4][5]. In the visible and near-infrared (VNIR) bands, the spectral reflectance of liquid water decreases with wavelength. Based on this principle, water surface information has been extracted using various spectral indices, including the Normalized Difference Vegetation Index (NDVI) [6], the Normalized Difference Water Index (NDWI) [7], and the Modified NDWI (MNDWI) [8]. More spectral bands are combined for enhanced water surface signals, such as the Automated Water Extraction Index (AWEI) [9] and the Open Water Likelihood (OWL) algorithm [10]. An increasing number of new algorithms are being developed for detailed and accurate water surface mapping [11][12][13][14]. In contrast, Wolski et al. [15] recently stated that the single short-wave infrared (SWIR) band also worked effectively.
The core idea of water surface mapping is to distinguish water from non-water objects. Due to their fine spatial resolutions, the National Aeronautics and Space Administration (NASA) Landsat result, and the hydrodynamic modeling result. The compositing algorithm and the resulting daily water surface area data are useful for hydro-ecological studies over seasonally flooded wetlands.

Study Area
Poyang Lake (28°22′-29°45′N, 115°47′-116°45′E) is the largest freshwater lake in China and is well known for its ecological significance [53]. The lake receives water from five major river systems within the lake basin and discharges into the Yangtze River through its north outlet ( Figure 1). Belonging to a humid subtropical climate zone, the annual precipitation in this lake basin is 1635.9 mm, and the annual mean surface air temperature was 17.5 °C for 1960-2010. Heavy rainfalls produce great surface discharges to the lake during the wet season from April-June. Rainfall declines sharply from July-September and remains relatively low from October-March. Jointly controlled by river discharges and interactions with the Yangtze River, the lake's inundation area experiences significant seasonal variations, with large areas occurring in wet seasons and small areas in dry seasons. Water surface area can vary from less than 1000 km 2 in winter during the dry season to over 3000 km 2 in summer during the wet season. The climate condition benefits vegetation development in this shallow lake [54]. As a result, a unique lake-vegetation system exists (Figure 1), posing a challenge for accurate lake surface mapping.

Using the Minimum NDVI to Highlight Water Surfaces
Water bodies generally have negative NDVI values. Although the negative values vary as a result of differences in water turbidity and other mixtures, NDVI is an effective indicator of water bodies [6,36,44,48,49,[55][56][57]. In contrast, vegetation generally has positive NDVI values. The contrasting NDVI behaviors form the basis for vegetation-water separation in the lake-vegetation system. Clouds are common in humid areas, especially in wet seasons. In the VNIR bands, spectral responses of thick clouds are generally flat, corresponding to NDVI values approaching zero. Attenuated by thin clouds, signals reflected from the underlying surfaces are blurred and the contrast of the NDVI image is reduced. Specifically, the NDVI value decreases over vegetation and increases over water bodies. Given the close-to-zero cloud NDVI values, the sign of NDVI values will be kept for cloud contaminated water and vegetation pixels.
Built-up areas, shadows, and snow may not have a large impact on water surface mapping in the Poyang Lake area. Within the lake area, almost no pixels show spectral features of artificial structures at the 250-m spatial scale due to sporadic villages and dense vegetation cover. The effects of shadows and snow can also be neglected in this flat and low-latitude area. Therefore, MODIS images consist mainly of water bodies, vegetation, and clouds (thick and thin) over the lake area. The exceptions may include some mudflat and sand banks [11,53] during the water recession period when vegetation has not yet regenerated. With an intermediate NDVI value between thick clouds and water bodies, these wet elements might be mistaken as water bodies by using image segmentation methods [58]. Despite this, the NDVI minima more likely correspond to a water pixel because water bodies always have lower NDVI values. Except for consecutive days of thick clouds, cloud free NDVI images can be generated by using the Minimum Value Composite (MinVC) algorithm (Table 1). In wet seasons with more cloudy days, a longer compositing period might increase the possibility of generating totally cloud free NDVI. The effects of cloud shadow are similar to cloud effects and can also be reduced by the MinVC algorithm. MODIS is the key instrument onboard Terra and Aqua satellites that have operated successfully since their launches in December 1999 and May 2002. The complementary a.m. and p.m. satellites provide twice-daily global daytime observations at 250-m, 500-m, and 1-km resolutions. MODIS data in the VNIR bands have been well calibrated, corrected for atmospheric effects, and released as 250-m surface (SUR) reflectance products (Terra: MOD09Q1 and Aqua: MYD09Q1, the latest Collection 061 version used in this study, hereinafter referred to as MOD09). NDVI data are calculated by using the red band (620-670 nm) and the near-infrared band (841-876 nm) reflectance data, and composited into a 16-day product with the CV-MVC algorithm [59]. The product is released as MOD13Q1 (Terra) and MYD13Q1 (Aqua), intended to highlight vegetation information. In each year, MOD13Q1 starts from DOY001-353 (DOY means Day of Year), and MYD13Q1 from DOY008-361. A combination Remote Sens. 2020, 12, 700 5 of 21 of MOD13Q1 and MYD13Q1 produces an equivalent 8-day NDVI dataset (hereinafter referred to as MOD13). In this study, all MOD13 data (the latest Collection 006 version) were used as of 2017. The dataset has masked out some water bodies by assigning a defaulted NDVI value of −0.3.

Generating the MinVC NDVI Data
Daily NDVI data were first calculated by using the 250-m MOD09Q1 (2000-2017) and MYD09Q1 (2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017) surface reflectance data in the red band (620-670 nm) and the near-infrared band (841-876 nm). Then, the pixelwise minimum NDVI values were determined by using the combined Terra and Aqua MODIS NDVI data within a 15-day moving window. The 15-day compositing period inherited all advantages of the 16-day compositing period used to produce the MOD13 NDVI product. The odd number meant the same days (i.e., one week) before and after a target day. More specifically, the ith-day MinVC NDVI image was created by using the daily MODIS NDVI images within the day i-7 and day i+7. Before the Aqua era, note that only Terra data were used. Different from the MODIS NDVI CV-MVC algorithm, this algorithm generated temporally continuous NDVI data. For example, the composite data at DOY1, 2010, also made use of MODIS NDVI data from DOY359-365, 2009.
In addition to the baseline algorithm, MinVC NDVI variants were generated by using a different compositing period or data source. First, a shorter 7-day moving window was used. The 15-day and 7-day periods were compared with respect to their performances on cloud removal. Second, the top-of-atmosphere (TOA), instead of SUR reflectance, was used. In this case, the MOD02QKM (Terra) and MYD02QKM (Aqua) products (the latest Collection 061 version) were used. The data were initially calibrated to TOA reflectance, and then followed the processing steps for SUR reflectance data. The two data sources were compared to investigate atmospheric effects on the MinVC NDVI data. In sum, the MinVC NDVI data had four variants, including the 15-day SUR reflectance based (baseline), the 7-day SUR reflectance based, the 15-day TOA reflectance based, and the 7-day TOA reflectance based. In this study, we focused on the baseline MinVC NDVI data, because the MOD13 NDVI product was generated by using SUR reflectance data with a similar compositing period.

Water Surface Mapping
An active contour model ("snake detector") was used to segment the daily MinVC NDVI images. Active contour models are developed within the framework of energy minimization. Generally, these models include an energy function and a terminal function [60]. The energy function is composed of two parts, an internal one to control the smoothness of the curve (evolving boundary) and an external one to attract the curve towards the desired boundary. The external energy can be provided by but not limited to the image gradient. Chen and Vese [61] propose an active contour model based on the Mumford-Shah segmentation technique and the level-set method. The model has an energy function where C is the evolving boundary, L is the length of C, inside(C) and outside(C) are the areas inside and outside C, with mean image values of c 1 and c 2 , and u 0 is the image value at coordinate (x, y). Minimizing Equation (1) generates the desired boundary C 0 . The model does not necessarily depend on image gradient, works well for noisy images, and is less dependent on the position of initial conditions [61]. These features are positive for water surface mapping over shallow lakes, for example, the Poyang Lake. First, these lakes possess varying degrees of water turbidity, especially near river estuaries. The heterogeneous and noisy lake surfaces can be discerned by using this model. Second, the lake bottomland is topographically complicated, and the transition between water and land phases is frequent as a result of hydrological variability in the Poyang Lake. The model can deal with such topological changes in the lake surface by using identical initial conditions. In this study, the model was applied to the daily MinVC NDVI and the MOD13 Remote Sens. 2020, 12, 700 6 of 21 NDVI. To implement the model, 40×40 uniformly distributed circles at a radius of 1.5 km were used as initial conditions. These processing steps were completed on the MATLAB R2015a (MathWorks Inc.) platform.

Water Surface Data Comparison and Validation
Water surface area data derived from MinVC NDVI were comprehensively evaluated. First, the composting algorithms were evaluated for their effectiveness to highlight water surface information (7-day vs. 15-day; SUR vs. TOA reflectance) (see Section 3.1). Second, the MinVC NDVI based result was validated by reference to the 30-m Landsat based result (NDWI threshold method), which was used as ground truth (see Figure 5 in Section 3.2). Within an image area including a land-water boundary, the NDWI histogram was analyzed to determine a single threshold value for land-water separation. Third, the MinVC NDVI based result was compared to the daily MIKE21 hydrodynamic modeling result [62][63][64][65] (see Figure 6 in Section 3.2). To drive the MIKE21 model, runoff data were obtained from gauged stations, and runoffs from the ungauged area were calculated by using a linear extrapolation of the gauged runoffs [63]. The total runoffs were specified as the upstream boundary conditions at the outlets of the five catchments (defined by the five major rivers in Figure 1). Daily water stages at the Hukou hydrological station (near the lake outlet) were used to define the downstream boundary condition. Values of other parameters followed those used in [64]. Due to the limitation of input data, the model was only implemented from 2000-2012. For hydrodynamic simulation, the whole lake area was divided as 20450 triangular elements, ranging from 0.007-0.605 km 2 . Those elements with water depth >0.4 m were counted as inundated, and the areas of inundated elements were summed as total water surface area [64]. The modeling result has been validated by using gauged water stage and river runoff data [64] and has been recently used for investigating lake surface dynamics [65]. Fourth, the MinVC NDVI based result was compared to other MODIS based results, including the cloud free MODIS NDVI based result (using the NDVI based thresholding method), the MOD13 NDVI based result (using the same active contour model), and the 30-m Landsat-MOD13 NDVI fusion based result [37,66] (see Section 3.3). For all comparisons, the coefficient of determination (R 2 ), root mean square error (RMSE), regression slope, and intercept were used as statistical metrics. To sum up, Table 2 shows the list of all water surface area data.

Calculation of Water Submersion Time and Trend in Water Surface Area
The spatial distribution and temporal evolution of water surfaces were analyzed by using the MinVC NDVI based water surface area data. In each year from 2000-2017, daily water masks were generated by assigning the value 1 for water surfaces and the value 0 for non-water surfaces. Spatially, the annual maps of water submersion time (or inundation frequency) were produced by averaging daily water surfaces within each year [38,42,53]. The values of water submersion time may vary from 0%-100%. The map of annual mean water submersion time (2000-2017) was calculated by averaging the 18-year maps. Temporally, time series water surface area data were presented, and the temporal trends were calculated. The differences among multiple water surface area data and the potential causal factors were investigated. To sum up, Figure 2 shows the flowchart of this study.

MinVC NDVI Performance on Highlighting Water Surfaces
The MinVC NDVI highlighted water surfaces by suppressing cloud effects. Figure 3 shows MOD13 NDVI and four MinVC NDVI composites (7-day and 15-day; TOA and SUR) in DOY001, 2003 (dry season) and DOY177, 2003 (wet season). MOD13 NDVI showed more details on vegetation distribution over land in both seasons. Over water surfaces, some pixels were assigned the defaulted NDVI value (-0.3) and others even showed a typical vegetation spectral feature (Figure 3a,f), which might have an impact on water surface mapping. By contrast, the four NDVI composites had closeto-zero values over land and negative values over water. Both land and water areas had spatially homogenous NDVI values. The 15-day MinVC NDVI had more negative NDVI values and exhibited water surfaces than the 7-day MinVC NDVI. The latter might still be obscured by clouds, especially in wet seasons (Figure 3g,i). By using more than doubled MODIS images, the 15-day MinVC NDVI reduced the risk being affected by clouds. No marked differences were observed between the TOA and SUR reflectance based NDVI composites, except that the latter had overall lower values. As expected, large negative NDVI values were observed occasionally for the latter due to inaccurate atmospheric correction.

MinVC NDVI Performance on Highlighting Water Surfaces
The MinVC NDVI highlighted water surfaces by suppressing cloud effects. Figure 3 shows MOD13 NDVI and four MinVC NDVI composites (7-day and 15-day; TOA and SUR) in DOY001, 2003 (dry season) and DOY177, 2003 (wet season). MOD13 NDVI showed more details on vegetation distribution over land in both seasons. Over water surfaces, some pixels were assigned the defaulted NDVI value (-0.3) and others even showed a typical vegetation spectral feature (Figure 3a,f), which might have an impact on water surface mapping. By contrast, the four NDVI composites had close-to-zero values over land and negative values over water. Both land and water areas had spatially homogenous NDVI values. The 15-day MinVC NDVI had more negative NDVI values and exhibited water surfaces than the 7-day MinVC NDVI. The latter might still be obscured by clouds, especially in wet seasons (Figure 3g,i). By using more than doubled MODIS images, the 15-day MinVC NDVI reduced the risk being affected by clouds. No marked differences were observed between the TOA and SUR reflectance based NDVI composites, except that the latter had overall lower values. As expected, large negative NDVI values were observed occasionally for the latter due to inaccurate atmospheric correction.
The compositing period and data source had an impact on water surface mapping results. Figure 4 shows concurrent MODIS NDVI and four MinVC NDVI composites at DOY175, 2016, (wet season) and DOY351, 2016 (dry season), when cloud free Landsat images are available within the composting period (indicating the availability of cloud free MODIS images). With at least one cloud free MODIS image, the 7-day and 15-day NDVI composites generated similar lake surface area data in the wet season. The conclusion was true for both the TOA reflectance based and the SUR reflectance based NDVI composites. In the dry season, the 15-day MinVC NDVI generated larger water surface area data than the 7-day counterpart. Compared to the Landsat based result, the baseline algorithm (15-day SUR reflectance) might underestimate water surface areas in wet seasons and perform better in dry seasons. In addition, the cloud free MODIS NDVI based result also underestimated lake surface areas relative to the Landsat based result, which might be due to scaling issues. It seemed that TOA NDVI outperformed SUR NDVI for lake surface mapping, especially in wet seasons.
Remote Sens. 2020, 12, 700 8 of 22 The compositing period and data source had an impact on water surface mapping results. Figure 4 shows concurrent MODIS NDVI and four MinVC NDVI composites at DOY175, 2016, (wet season) and DOY351, 2016 (dry season), when cloud free Landsat images are available within the composting period (indicating the availability of cloud free MODIS images). With at least one cloud free MODIS image, the 7-day and 15-day NDVI composites generated similar lake surface area data in the wet season. The conclusion was true for both the TOA reflectance based and the SUR reflectance based NDVI composites. In the dry season, the 15-day MinVC NDVI generated larger water surface area data than the 7-day counterpart. Compared to the Landsat based result, the baseline algorithm (15-day SUR reflectance) might underestimate water surface areas in wet seasons and perform better in dry seasons. In addition, the cloud free MODIS NDVI based result also underestimated lake surface areas relative to the Landsat based result, which might be due to scaling issues. It seemed that TOA NDVI outperformed SUR NDVI for lake surface mapping, especially in wet seasons.

Accuracy of MinVC NDVI Based Water Surface Area Data
The baseline MinVC NDVI based water surface area data were in good agreement with both the 30-m Landsat based result (R 2 = 0.914, Figure 5) and the daily hydrodynamic modeling result (R 2 = 0.799, Figure 6). The former elucidated the absolute accuracy (RMSE = 223.7 km 2 ) and the latter

Accuracy of MinVC NDVI Based Water Surface Area Data
The baseline MinVC NDVI based water surface area data were in good agreement with both the 30-m Landsat based result (R 2 = 0.914, Figure 5) and the daily hydrodynamic modeling result (R 2 = 0.799, Figure 6). The former elucidated the absolute accuracy (RMSE = 223.7 km 2 ) and the latter elucidated the temporal consistency of the derived water surface area data. Our result might underestimate water surface area in wet seasons (corresponding to large water surfaces in Figure 5) and slightly overestimate water surface area in dry seasons (corresponding to small water surfaces in Figure 5). The MIKE21 hydrodynamic model had known artifacts (e.g., coarse grid size; topography data not updated) in modeling water regimes over topographically flat lake bottoms (e.g., small lakes in the western lake areas). These were part of the reasons why large discrepancies were observed between the MinVC NDVI based and the MIKE21 based water surface area data (RMSE = 560.7 km 2 , Figure 6).

Performance of MODIS Based Water Surface Area Data
The baseline MinVC NDVI based water surface area data were comparable to other MODIS based results ( Figure 7). As previously illustrated in Figure 4, the SUR reflectance based MinVC NDVI might generate smaller water surface area data than the TOA reflectance based MinVC NDVI,

Performance of MODIS Based Water Surface Area Data
The baseline MinVC NDVI based water surface area data were comparable to other MODIS based results ( Figure 7). As previously illustrated in Figure 4, the SUR reflectance based MinVC NDVI might generate smaller water surface area data than the TOA reflectance based MinVC NDVI, especially in wet seasons (Figure 7a). Despite this, the two results were in good agreement with each

Performance of MODIS Based Water Surface Area Data
The baseline MinVC NDVI based water surface area data were comparable to other MODIS based results ( Figure 7). As previously illustrated in Figure 4, the SUR reflectance based MinVC NDVI might generate smaller water surface area data than the TOA reflectance based MinVC NDVI, especially in wet seasons (Figure 7a). Despite this, the two results were in good agreement with each other (R 2 = 0.833). Compared to the cloud free MODIS NDVI based result, our result might overestimate water surface areas in dry seasons and underestimate water surface areas in wet seasons (Figure 7b). The two results were also highly consistent with each other (R 2 = 0.824). For NDVI composite product, MOD13 NDVI generated overall larger water surface area data relative to the MinVC NDVI, especially for small water surfaces (Figure 7c). Integrated with information from Landsat data, the MOD13 NDVI was more powerful in water surface delineation. As a result, the Landsat-MOD13 NDVI generated a more consistent result with the MinVC NDVI (Figure 7d). The former had a systematical positive bias value of~211.7 km 2 , and the regression slope value was very close to 1.
Remote Sens. 2020, 12, 700 11 of 22 overestimate water surface areas in dry seasons and underestimate water surface areas in wet seasons (Figure 7b). The two results were also highly consistent with each other (R 2 = 0.824). For NDVI composite product, MOD13 NDVI generated overall larger water surface area data relative to the MinVC NDVI, especially for small water surfaces (Figure 7c). Integrated with information from Landsat data, the MOD13 NDVI was more powerful in water surface delineation. As a result, the Landsat-MOD13 NDVI generated a more consistent result with the MinVC NDVI (Figure 7d). The former had a systematical positive bias value of ~211.7 km 2 , and the regression slope value was very close to 1. The baseline MinVC NDVI based result might overestimate water surface areas in dry seasons, in particular in rapid water surface expansion and recession periods. A comparison of the MinVC NDVI based and the cloud free MODIS NDVI based results (Figure 8) showed that positive differences (MinVC relative to cloud free) were large (up to 400 km 2 ) in March and April (normally corresponding to rapid water surface expansion period) and in October and November (normally corresponding to rapid water surface recession period). Negative differences were observed in June-August, which might be caused by incomplete cloud removal due to permanent cloud coverage in wet seasons.  The baseline MinVC NDVI based result might overestimate water surface areas in dry seasons, in particular in rapid water surface expansion and recession periods. A comparison of the MinVC NDVI based and the cloud free MODIS NDVI based results (Figure 8) showed that positive differences (MinVC relative to cloud free) were large (up to 400 km 2 ) in March and April (normally corresponding to rapid water surface expansion period) and in October and November (normally corresponding to rapid water surface recession period). Negative differences were observed in June-August, which might be caused by incomplete cloud removal due to permanent cloud coverage in wet seasons.

Spatial Patterns of Lake Water Surface
The baseline MinVC NDVI based result provided temporally continuous water surface area data, based on which the evolution of the flood/drought event can be investigated. Figure 9 illustrates monthly water surface distribution for the Poyang Lake in 2010, a known typical flooding year. In April, the lake surface expanded rapidly to a large area and stayed until the end of September. The long-lasting large water surfaces made the year 2010 as the most severe flooding year in the recent decade. It is well known that transitions between flooding and drought events are common in the Poyang Lake. Figure 10 illustrates monthly water surface distribution in the next year, the most severe drought year in the recent decade. The lake surface remained at a small area until the end of May. With a relatively large area in June (still smaller than the lake area in June 2010), the water surface began to shrink until the end of year. During the drought year, the small isolated lakes in the western lake area were almost dried, which might have caused significant ecological consequences. In contrast, these lakes can be observed during the flooding year, especially from October-December, 2010 ( Figure 9). An average of these daily water surface maps corresponded to an overall description of water submersion time in the Poyang Lake ( Figure 11).

Spatial Patterns of Lake Water Surface
The baseline MinVC NDVI based result provided temporally continuous water surface area data, based on which the evolution of the flood/drought event can be investigated. Figure 9 illustrates monthly water surface distribution for the Poyang Lake in 2010, a known typical flooding year. In April, the lake surface expanded rapidly to a large area and stayed until the end of September. The long-lasting large water surfaces made the year 2010 as the most severe flooding year in the recent decade. It is well known that transitions between flooding and drought events are common in the Poyang Lake. Figure 10 illustrates monthly water surface distribution in the next year, the most severe drought year in the recent decade. The lake surface remained at a small area until the end of May. With a relatively large area in June (still smaller than the lake area in June 2010), the water surface began to shrink until the end of year. During the drought year, the small isolated lakes in the western lake area were almost dried, which might have caused significant ecological consequences. In contrast, these lakes can be observed during the flooding year, especially from October-December, 2010 ( Figure 9). An average of these daily water surface maps corresponded to an overall description of water submersion time in the Poyang Lake ( Figure 11).
Dry and wet years were determined from the annual maps of water submersion time ( Figure 12). The most severe drought event in 2011 was characterized by extremely short water submersion time all over the lake area except for some isolated reservoirs (Figure 12l). The known drought events in 2006 and 2013 were also manifested on the water submersion maps (Figure 12g,n), yet with lesser strength. Jointly controlled by meteorological (e.g., precipitation) and hydrological (inflow and outflow) conditions, smaller water surfaces may occur in some but not all months. Therefore, the maps might not show a markedly short water submersion time (e.g., 2006). Nevertheless, the wet years (2010, 2012, and 2016) were generally associated with long submersion time (Figure 12k,m,q). It seemed that the 2010 flooding event was the most severe in the recent decade, followed by the 2016 and the 2012 flooding events. Note that the results in 2000-2002 were less reliable because Aqua MODIS images were not included for NDVI compositing until July 2002. Similar results have been reported in [38,67,68].
May. With a relatively large area in June (still smaller than the lake area in June 2010), the water surface began to shrink until the end of year. During the drought year, the small isolated lakes in the western lake area were almost dried, which might have caused significant ecological consequences. In contrast, these lakes can be observed during the flooding year, especially from October-December, 2010 ( Figure 9). An average of these daily water surface maps corresponded to an overall description of water submersion time in the Poyang Lake ( Figure 11).      Dry and wet years were determined from the annual maps of water submersion time ( Figure  12). The most severe drought event in 2011 was characterized by extremely short water submersion time all over the lake area except for some isolated reservoirs (Figure 12l). The known drought events in 2006 and 2013 were also manifested on the water submersion maps (Figure 12g,n), yet with lesser strength. Jointly controlled by meteorological (e.g., precipitation) and hydrological (inflow and outflow) conditions, smaller water surfaces may occur in some but not all months. Therefore, the maps might not show a markedly short water submersion time (e.g., 2006). Nevertheless, the wet years (2010, 2012, and 2016) were generally associated with long submersion time (Figure 12k

Temporal Variations in Lake Water Surface Area
Temporally continuous NDVI generated opposite trends in water surface area to the unevenly distributed cloud free NDVI. Within a common period in 2003-2016, water surface area increased by 9.3 km 2 ·a -1 as derived using the MODIS MinVC NDVI (daily resolution), and increased by 8.2 km 2 ·a -1 Figure 12. Annual maps of water submersion time for the Poyang Lake in 2000-2017 derived from the baseline MinVC NDVI based water surface area data.

Temporal Variations in Lake Water Surface Area
Temporally continuous NDVI generated opposite trends in water surface area to the unevenly distributed cloud free NDVI. Within a common period in 2003-2016, water surface area increased by 9.3 km 2 ·a −1 as derived using the MODIS MinVC NDVI (daily resolution), and increased by 8.2 km 2 ·a −1 as derived using the Landsat-MOD13 NDVI (8-day resolution) (Figure 13). Large water surfaces in the recent wet years (2010, 2012, and 2016) might have reversed the declining trends in water surface area during the early 2000s. In contrast, the cloud free NDVI still revealed a declining trend at a rate of 8.2 km 2 ·a −1 . The uneven temporal distribution of cloud free MODIS images, more in dry seasons and less in wet seasons (Figure 14a), might account for the opposite trends. Although long-lasting large water surfaces might have occurred in recent wet years, the number of totally cloud free images decreased as a result of continuous cloudy and rainy days. Figure 14b shows that the cloud free NDVI based water surface area data followed a Gumbel distribution (R 2 = 0.919) with a most likely water surface area of 865.5 km 2 .
Remote Sens. 2020, 12, 700 15 of 22 area during the early 2000s. In contrast, the cloud free NDVI still revealed a declining trend at a rate of 8.2 km 2 ·a -1 . The uneven temporal distribution of cloud free MODIS images, more in dry seasons and less in wet seasons (Figure 14a), might account for the opposite trends. Although long-lasting large water surfaces might have occurred in recent wet years, the number of totally cloud free images decreased as a result of continuous cloudy and rainy days. Figure 14b shows that the cloud free NDVI based water surface area data followed a Gumbel distribution (R 2 = 0.919) with a most likely water surface area of 865.5 km 2 .

Added value of the MinVC NDVI for water surface mapping
Cloudy and rainy days lead to less cloud free images in wet seasons when water surfaces are, in the meantime, undergoing rapid expansion. High frequency lake surface data are needed in response to this contradiction. In addition, temporally continuous data are required to generate a reliable trend in water surface area. Fortunately, SAR images (e.g., Sentinel-1 at a 10-m resolution) are now capable of high frequency and high-resolution water surface mapping, although limited data are available for tracking long-term water surface area changes (e.g., dating back to 1980s) compared to optical remote sensing data. The passive microwave radiometers and geostationary sensors offer candidate data for all-sky or high frequency water surface monitoring; however, their spatial resolutions are insufficient for most inland water bodies [3,50]. MODIS sensors onboard Terra and Aqua satellites well balance the requirements of spatial and temporal resolutions, and play an important role in water surface mapping. There is a general tendency to exploit the entire MODIS data archive, from both Terra and Aqua, for added information in change detection studies [49,51]. These data are important for tracking rapidly changing natural processes at daily scales, including water surface changes. In this context, this study proposes a simple yet effective method for daily water surface mapping.
The NDVI MinVC algorithm suppresses cloud effects by using complementary cloud free or thin cloud covered information on bi-monthly Terra and Aqua MODIS images. There are at least 30 Terra and Aqua MODIS images within the 15-day moving window, and the number increases toward highlatitude regions, increasing the probability to generate totally cloud free NDVI images. In addition to cloud removal, the algorithm can deal with other issues that complicate water surface mapping. Two of them are sun glint and water turbidity, both leading to increased NDVI values of water pixels. Generally caused by transient hydrological (e.g., large river runoff as a result of heavy rainfall) and meteorological (e.g., strong wind) events, these effects can be reduced by selecting the NDVI minima. The algorithm even improves the effective spatial resolution of water surface maps. Because MODIS is a whiskbroom scanner, the 250-m resolution in nadir-viewing direction can be reduced to over 1km across track and 500-m along track. Fortunately, our algorithm inclines to select nadir-viewing water pixels instead of land-water mixed pixels that generally have larger NDVI values.

Added Value of the MinVC NDVI for Water Surface Mapping
Cloudy and rainy days lead to less cloud free images in wet seasons when water surfaces are, in the meantime, undergoing rapid expansion. High frequency lake surface data are needed in response to this contradiction. In addition, temporally continuous data are required to generate a reliable trend in water surface area. Fortunately, SAR images (e.g., Sentinel-1 at a 10-m resolution) are now capable of high frequency and high-resolution water surface mapping, although limited data are available for tracking long-term water surface area changes (e.g., dating back to 1980s) compared to optical remote sensing data. The passive microwave radiometers and geostationary sensors offer candidate data for all-sky or high frequency water surface monitoring; however, their spatial resolutions are insufficient for most inland water bodies [3,50]. MODIS sensors onboard Terra and Aqua satellites well balance the requirements of spatial and temporal resolutions, and play an important role in water surface mapping. There is a general tendency to exploit the entire MODIS data archive, from both Terra and Aqua, for added information in change detection studies [49,51]. These data are important for tracking rapidly changing natural processes at daily scales, including water surface changes. In this context, this study proposes a simple yet effective method for daily water surface mapping.
The NDVI MinVC algorithm suppresses cloud effects by using complementary cloud free or thin cloud covered information on bi-monthly Terra and Aqua MODIS images. There are at least 30 Terra and Aqua MODIS images within the 15-day moving window, and the number increases toward high-latitude regions, increasing the probability to generate totally cloud free NDVI images. In addition to cloud removal, the algorithm can deal with other issues that complicate water surface mapping. Two of them are sun glint and water turbidity, both leading to increased NDVI values of water pixels. Generally caused by transient hydrological (e.g., large river runoff as a result of heavy rainfall) and meteorological (e.g., strong wind) events, these effects can be reduced by selecting the NDVI minima. The algorithm even improves the effective spatial resolution of water surface maps. Because MODIS is a whiskbroom scanner, the 250-m resolution in nadir-viewing direction can be reduced to over 1-km across track and 500-m along track. Fortunately, our algorithm inclines to select nadir-viewing water pixels instead of land-water mixed pixels that generally have larger NDVI values.

Causes of Differences among Multiple Lake Surface Data
The 30-m Landsat images can accurately delineate water surfaces at the days of image acquisition. At a 250-m resolution, the cloud free MODIS NDVI generates high frequency (relative to Landsat, 6-day effective temporal interval, unevenly distributed) lake surface information. Concurrent daily Terra and Aqua images further corroborate the derived water surface area data, which serve as a solid basis for algorithm validation. By reference to the cloud free NDVI, the MinVC NDVI overestimates water surface areas for small lake surfaces and underestimates water surface areas for large lake surfaces. By selecting minimum NDVI values, the NDVI MinVC algorithm is inclined to include more wet pixels in dry seasons. Generally, multiple scenes of cloud free or thin cloud covered images are available in dry seasons. The algorithm tends to reserve more pixels in the second half of the compositing period (if any) at the water rising stage, while reserving more pixels in the first half of the composition period (if any) at the water recession stage. This situation is more serious at rapid water rising and recession stages. Although lake surface areas tend to be overestimated in the transitional periods, the temporal pattern might be kept by using the NDVI MinVC algorithm. In this circumstance, the 15-day moving window functions like a filter that selects the maximum value in the 15-day compositing period. Because of the 15-day compositing period, the maximum phase difference between our result and the actual lake surface change can be one week in advance (having cloud free image with the largest lake surface) or one week lagging behind (having cloud free image with the smallest lake surface). The accurate phase difference, however, varies depending on the availability of cloud free MODIS images. Meanwhile, the occurrence of wet sand and mudflat (both showing negative NDVI values) can be confused with water surfaces at the water recession stage when vegetation has not yet regenerated. In wet seasons, however, residual cloud effects might contribute to the underestimation of lake surface areas. Given the good consistency between the MinVC NDVI based and the Landsat-MOD13 NDVI based results (Figures 7d and 13a,b), the issues of overestimation and underestimation are likely shared by image compositing based methods.
Various factors may determine the differences among multiple water surface area data. Intended to highlight vegetation information, the MOD13 NDVI seems to reduce NDVI contrast over the lake area. Therefore, it generates overall larger lake surface areas compared to the MinVC NDVI by using the same active contour model. Nevertheless, the Landsat-MOD13 fusion technique decomposes the original MOD13 NDVI using the Landsat NDVI, which largely improves the ability of MOD13 NDVI for water surface mapping. The MIKE21 model also generates overall larger water surface areas than the MinVC NDVI. Tan et al. [37] showed the poor performance of MIKE21 to model water inundation in topographically flat lake areas. Zhang et al. [65] also found that the modeling results overestimated water surface areas that were determined from cloud free MODIS images, which is in agreement with our conclusion. The differences might also be caused by topographical changes in the Poyang Lake. Topography data collected in 1998 were used to drive the model, which cannot reveal the effects of changing topography on lake surface areas in recent decades [62]. Specifically, sand mining near the north lake outlet has been reported to increase the lake's discharge to the Yangtze River and probably caused the decline in lake surface areas, particularly in dry seasons [69].

Improvements of Current NDVI MinVC Algorithm
Long compositing period (i.e., 15-day) has an impact on water surface mapping. Currently, the baseline algorithm uses no prior knowledge and assumes a relaxed uniform compositing period. The main objective is to derive cloud free NDVI images in wet seasons as much as possible. However, the treatment also exacerbates the situation of lake surface area overestimation in dry seasons. To reduce such effects, MODIS cloud mask data can be used to determine time-varying compositing periods. It is safe to assume that the compositing period is inversely proportional to the overall cloud amount.
More specifically, no cloud coverage (totally cloud free) means a daily compositing period, and larger cloud amount corresponds to a longer compositing period. Because the mean temporal interval of cloud free MODIS images is~six days, such time-varying compositing periods are reasonable and are expected to perform better than the fixed 15-day compositing period. With regard to data source, TOA instead of surface reflectance is recommended for NDVI MinVC composting, although minor differences are observed between results derived from the two composites.
A combination of multiple remote sensing data enhances the possibility of cloud free NDVI data. Within a day, satellite images acquired at a different time from Terra and Aqua MODIS images (10:30 and 13:30 Local Time) can provide complementary cloud free information. Such sensors include, but are not limited to, the ENVISAT MERIS (10:00 Local Time, 2002-2012), the Sentinel-3 OLCI (10:00 Local Time, 2016-present), and the NOAA/MetOp AVHRR (overpassing times varying dependent on orbital drift, 1970s-present). Prior to the MODIS era, the AVHRR data can be used alone to generate MinVC NDVI. Several contemporary a.m. and p.m. orbiting AVHRR sensors may overpass the earth surface at different times, especially when considering the orbital drift issues for earlier sensors. By applying the MinVC algorithm to those dense NDVI images, cloud effects can be further reduced. Due to different spectral response functions (relative spectral responses), however, the inter-sensor NDVI differences should be adjusted in advance before any compositing. In most cases, a linear equation is qualified for a reliable spectral adjustment [70][71][72]. Within the framework of NDVI intercalibration, moderate-resolution and high-resolution images can be combined for water surface mapping. Note that the proposed algorithm can be accommodated to other indices that are sensitive to water surfaces.
High-resolution satellite images are required for improved water surface mapping. The relevant sensors include the NASA Landsat series sensors, the ESA Sentinel-2 Multispectral Imager (MSI), and the Chinese HJ/GF series sensors. These sensors data can be employed in at least two aspects. First, a synergy of high-resolution images (e.g., Landsat-8 OLI and Sentinel-2 MSI as in [73,74]) provide an improved reference dataset for validation purposes. Second, high-resolution satellite images can be integrated with the MinVC NDVI to derive long-term high-resolution water surface area data (e.g., [36,37]). These data are expected to benefit monitoring small lakes with an area even in the order of several MODIS pixels [66].

Conclusions
This study proposed a simple MODIS NDVI based MinVC algorithm to generate daily water surface area data over seasonally flooded wetlands by using the combined daily 250-m Terra and Aqua MODIS NDVI images. The baseline algorithm first created daily Minimum Value Composite NDVI by selecting pixelwise NDVI minima within a 15-day moving window. Consisting mainly of water pixels over water and cloud pixels over land, the daily MinVC NDVI data were then segmented for water surfaces from the cloud background. The derived water surface result compared well with the Landsat based, the cloud free MODIS image based, the hydrodynamic modeling based, and the Landat-MOD13 image fusion based results, and outperformed the MOD13 NDVI based result. Our result provided high frequency and reliable descriptions of water surface dynamics in both space and time. The proposed algorithm has the potential to make full use of fragmented cloud free information from the entire MODIS VNIR data archive for water surface mapping of seasonally flooded wetlands. Within the framework of satellite intercalibration, the algorithm can integrate multiple satellite senor data, both moderate-resolution and high-resolution, for improved water surface mapping.
Author Contributions: Conceptualization, X.F. and Y.L.; methodology, X.F.; validation, G.W. and X.Z.; writing-original draft preparation, X.F.; writing-review and editing, Y.L. All authors have read and agree to the published version of the manuscript.
Funding: This study was funded by the National Natural Science Foundation of China (NSFC), grant number. 41701414.