Monitoring Spatio-Temporal Dynamics in the Eastern Plain Lakes of China Using Long-Term MODIS UNWI Index

: Monitoring the spatio-temporal dynamics of the Eastern Plain Lake (EPL) is vital to the local environment and economy. However, due to the limitations and efﬁciency of traditional image formats in storing and processing large amounts of images and optimal threshold adjustments are often necessary for water/non-water separation based on traditional multi-band/spectral water indexes over large areas and in the long-term, previous studies have either been on a short period or mainly focused on water inundation dynamics of several lakes. To address these issues, a multi-dimensional dataset (MDD) storage format was used to efﬁciently organize more than ~7000 time series composite MODIS images. Furthermore, a universal normalized water index (UNWI) was developed based on full-spectrum information to simplify optimal threshold adjustments. Consequently, the present study analyzed the patterns of spatio-temporal water dynamic patterns and potential driving factors of inundation changes at large lakes (>5 km 2 ) in the EPL during 2000–2020 through MDD and UNWI. In terms of annual inundation patterns, the numbers of lakes that experienced signiﬁcant ( p < 0.05) decreases (17 lakes) and increases (43 lakes) were highest for Class IV lakes among six geographical classes. Variation in intra-annual inundation in Classes I and II is correlated with consumption of chemical fertilizers (CCF), while precipitation accounted for the most change in lake area in Class III. This spatio-temporal analysis of lakes provides a necessary foundation for the sustainable development and continuous investigations of the EPL.


Introduction
As crucial components of the inland water system, lakes hold indispensable ecological value [1][2][3][4]. In addition to providing a majority of the freshwater resources necessary for the survival and development of humankind, lakes also perform a significant role in maintaining the balance of regional ecological environments via the hydrological cycle, flood control, fisheries, and biodiversity conservation [5][6][7][8]. Lakes are a part of the continental surfaces and highly sensitive to climate and human activities [7,[9][10][11][12][13][14].
In recent decades, the intensification of both global climate change and human activities have caused many lakes to undergo tremendous changes. Water resources security has become a serious problem, threatening nearly 80% of the world population [14]. Not only the number of lakes, but also the area of lakes have experienced statistically significant declines and this phenomenon has appeared in Oklahoma, USA [15], central Asia [16], Alaska, USA and Siberia, Russia [9], and so on. Moreover, the same phenomenon has occurred in saline lakes around the world [17]. Without exception, lakes in China have also substantially changed [10,18]. For example, glaciers melting caused by accelerated climate warming have led to the expansion of lakes on the Qinghai and Tibetan plateaus [11,[19][20][21][22]. Even more drastically, the lakes of the Inner-Mongolian Plateau have decreased by 1223.76 km 2 , while 333 lakes of eastern China with area ≥ 1.0 km 2 have completely disappeared over the past few decades [23]. Therefore, understanding the spatio-spatial dynamics and distribution of lakes in large regions is urgently required to better protect water resource ecosystems and environments [8,24].
Remote sensing methods are unparalleled in their ability to provide large-scale and repeatable Earth observations via real-time images that capture surface information in remote or inaccessible areas, which have become an important method of studying surface water information (size or number) [4,15,[24][25][26][27]. The imagery of various satellites, such as Landsat, Moderate Resolution Image Spectroradiometer (MODIS), and Sentinel, has been used to acquire and analyze water surface inundation [3,4,9,23,[28][29][30][31][32]. A variety of band combinations indexes, such as the Normalized Difference Water Index (NDWI) [33], Modified Normalized Difference Water Index (mNDWI) [34], Automated Water Extraction Index (AWEI), Automated Water Extraction Index shadow (AWEIsh) and Automated Water Extraction Index no shadow (AWEInsh) [25], Land Surface Water Index (LSWI) and Modified Land Surface Water Index (MLSWI) [35,36], Sentinel Water Mask (SWM) [37,38], and Normalized Difference Vegetation Index (NDVI), etc., have been widely used to acquire water-land interface data. Those indexes were used to find sandbars and lakes or reservoir. For example, SWM water index was the best indicator for sandbars detecting in the Lower Vistula River [39], and the first complete study of lake changes in China between the 1960s and 1980s and 2005-2006 was accomplished using more than 24,000 satellite images [10]. In addition, more than 3831 Landsat satellite images were used to collect data for lakes ≥ 1 km 2 in China every five years based on NDWI during 1990-2015 [23]; Tao et al. [11] discussed the changes and their driving forces of lakes in China from the mid-1980s to 2015, while [15] discussed the water body area and number in Oklahoma, USA, using~16,000 Landsat images based on mNDWI, NDVI, and EVI during 1984-2015. Moreover, [40] analyzed the intra-annual trend in Patagonian with different water index: NDWI1, NDWI2, and mNDWI. The MODIS images of NDVI were also used to acquire the area variation of Urmia Lake, Iran, from 2000 and 2011 [41]. At the same time, [42] used the support vector machine model to identify water and non-water area as well as estimated water volume of Lake Victoria during 2003-2017, with long-term MODIS and GRACE images. In addition, the remote sensing images are also used to acquire global inundation maps. For example, [43] produced global inundation maps at a scale of 15 arc seconds based on the topographic downscaling of the SRTM-derived Hydro SHEDS database and a GLC2000 global land cover map. At the same time, all 3.4 million Landsat 5, 7, and 8 scenes from 1999 to 2018 were recently used to generate global surface water maps [8]. Obviously, the optimal threshold for water/non-water separation based on traditional water indexes is key to accurately extracting water area data [23]. However, determining the optimal thresholds for long-term water/non-water separation over large areas poses a substantial challenge. Traditional water indexes are based on multi-band or multi-spectral information, and a large amount of other spectral information has not been utilized, which may be one reason why threshold adjustment would be necessary. Therefore, in the present study, we introduce a full-band index called the Universal Normalized Water Index (UNWI), which aims to make full use of all spectral information, thereby simplifying or reducing the need for threshold adjustment.
The middle and lower reaches of the Yangtze River and the Huai River, along with their associated outflows and freshwater lakes, constitute the main lakes in the EPL region of China [44,45]. The core EPL primarily covers Zhejiang, Jiangsu, Shandong, Hebei, Anhui, Jiangxi, Hubei, and Hunan provinces as well as Shanghai. With 643 lakes ≥ 1 km 2 , accounting for 25.9% of the total area of lakes in China [45,46], the EPL is one of the most biologically diverse wetland systems and the largest freshwater group of lakes on the planet, serving as the most important source of freshwater for human and commercial use in eastern China. Characterized by high urbanization and population density [10,45], it is Remote Sens. 2022, 14, 985 3 of 26 reported that the EPL lakes decreased in area by 7% during the 1960s to 2015 [23]. Due to the influence of human activities and accelerated urbanization [4,10,[47][48][49][50], many of the lakes and associated wetlands in this region have undergone dramatic degradation [3,29,51,52], leading to a series of problems such as reduced shortage of drinking water [48,53,54], drought [13,28,50], vegetation losses [47,51,52], etc. The 90 lakes in the downstream Yangtze Basin decreased by 7.4% from 2000 to 2011 [2], and the annual mean lake areas in Poyang also significantly declined from 2000 to 2010 [4]. However, few reports about the EPL as a whole region are available for recent (2000-2020) long-term water dynamics, such as spatial and temporal patterns, seasonal patterns, and inter-annual variability.
The quantification of trends in spatio-temporal water dynamics is critical to water security and will help elucidate and protect water resource-related ecosystem services; however, processing large amounts of long-term data for an adequate number of lakes can prove difficult and the traditional image data storage and processing methods are not ideal for large areas and long-term remote sensing images. The traditional methods are based on traditional image storage format, such as Geographic Tagged Image File Format (Geo-TIFF) [55] and formats for remote sensing business software (e.g., The Environment for Visualizing Images, ENVI). These traditional data storage formats are based on a three-dimensional cube, which stores 1D-spectral dimension and 2D-spatial dimension information, while long-term series remote sensing data contain four-dimensional information: 1D-spectral dimension, 2D-spatial dimensions, and 1D-temporal dimension, also known as temporal-spatial-spectral dimensional data [56]. Traditional data organization formats face several shortcomings in the storage, management, and analysis of long-term remote sensing data. First, in terms of data management, the need for multiple discrete data files to store and organize time series datasets results in a lack of integrity and consistency and a lack of temporal dimension information. Second, in data analysis, many repeated steps are required to process the data, which is cumbersome and makes it difficult to complete the corresponding work through interactive operations. Lastly, it is necessary to read each node one-by-one from multiple files, tree structures, or linked list structures of each file to obtain the data, resulting in low efficiency and inconvenient operations for the extraction of temporal-spectrum information. Therefore, a multi-dimensional dataset (MDD) format framework is urgently needed to store and manage long-term remote sensing images. The design targets of this present study are as follows:

1.
Build a new framework that can process time-series remote sensing data; 2.
Develop a universal normalized water index (UNWI) for large areas and long-term lake inundation observations that uses full spectrum information; 3.
Use the new framework to record patterns of spatio-temporal dynamics of large lakes in the EPL, China, during 2000-2020, and discuss the forces driving the inundation changes at large lakes in the EPL.

Study Area
The EPL is located in eastern China and contains 18% of the nation's lakes [10,23]. The region contains about~600 lakes larger than 1 km 2 [23], which are mainly distributed in the Yangtze Plain, Huai River Plain, and the canal from Beijing to Hangzhou ( Figure 1) [10,18,57]. These freshwater lakes provide not only freshwater resources for human survival but also serve many ecosystem functions, such as regulating river discharge, irrigation development, water for industry, shipping, etc. [25,29,51,[58][59][60]. More importantly, they also play a significant role in mitigating regional climate change, regulating flood and drought, combating soil erosion and water losses, and protecting biodiversity [4,6,14,43,61,62]. Locations of the 184 study lakes in the EPL region and MODIS-derived annual water frequency, (a) the EPL region and the inset shows the number of lakes with decreasing or increasing trends in annual area for the six classes, (b) the middle of the Yangtze Plain. Lakes for which the annual lake area significantly (p < 0.05) increased or decreased are marked with the symbol "↑" or "↓", respectively. Locations of the 184 study lakes in the EPL region and MODIS-derived annual water frequency, (a) the EPL region and the inset shows the number of lakes with decreasing or increasing trends in annual area for the six classes, (b) the middle of the Yangtze Plain. Lakes for which the annual lake area significantly (p < 0.05) increased or decreased are marked with the symbol "↑" or "↓", respectively.
In the present study, the 184 lakes that were >5 km 2 were selected, and information regarding name, longitude, latitude, and basic region of the lakes was obtained from the China Lake Scientific Database (http://www.lakesci.csdb.cn/, accessed on 14 January 2021) and National Tibetan Plateau Data Center (http://data.tpdc.ac.cn/, accessed on 23 February, 2021). The basic region of each lake was used to extract a subset of MODIS images and was greater in size than the maximum water area [23,63]. In the study of the Yangtze [2,5], the lakes were categorized into three classes based on their discriminant relations with the Yangtze River. Based on the same consideration and strategy, the lakes were categorized into six classes according to the connection between lakes and the Beijing-Hangzhou Grand Reach, Yangtze River, and Huai River in order to better understand the spatial and temporal dynamics of these lakes and their protentional driving forces. Lakes that are freely connected to the Yangtze River are categorized as Class I [2], including Poyang Lake, Dongting Lake, and Shijiu Lake. Lakes that are freely connected to the Huai River were categorized as Class II, while those freely connected to the Beijing-Hangzhou Grand Reach were classified as Class III. The remaining semi-natural lakes (mainly via local flood gates) of the Yangtze River Plain and Huai River Plain were categorized as Class IV and Class V, respectively. The remaining eight lakes were categorized as Class VI).

Datasets
The MODIS surface reflectance data (h27v05, h27v06, h28v05, and h28v06) from Terra (MOD09) and Acqua (MYD09) during 2000-2020 were obtained from the National Aeronautics and Space Administration Land Processes Distribution Active Archive Center (https://ladsweb.nascom.nasa.gov/, accessed on 3 March 2021) and were used to study changes in the spatial and temporal dynamics of the study lakes. Data on lake dynamics were obtained from the MODIS eight-day surface reflectance data, which avoids large errors and negative reflectance partly by selecting the best possible high observation coverage, low viewing angle, and aerosol loading during an eight-day period [5,57], with spatial resolutions of 250 m (MOD09Q1 and MYD09Q1) and 500 m (MOD09A1 and MYD09A1). Band 1 (red spectra of 620-670 nm) and band 2 (near-infrared spectra of 841-876 nm), which are stored at a spatial resolution of 250 m, and bands 3-7 (i.e., 459-479 nm, 545-565 nm, 1230-1250 nm, 1628-1652 nm and 2105-2155 nm), stored at a spatial resolution of 500 m, were re-sampled to 250 m using the sharpen method [4,5]. In addition, MODIS daily surface reflectance data (MOD09GQ and MYD09GQ with a spatial resolution of 250 m; MOD09GA and MOD09GA with a spatial resolution of 500 m) were also downloaded for further validation by higher-resolution (30 m) Landsat images. With the MODIS QA data, influencing factors such as clouds, cloud shadows, or snow can be detected and then excluded from water extraction as well as spatial and temporal dynamic analysis.
The higher resolution (30 m) Landsat Thematic Mapper (TM), and Landsat 8 Operational Land Imager (OLI) were obtained and downloaded from the United States Geological Survey (USGS, https://earthexplorer.usgs.gov/, accessed on 12 March 2021) to validate the water extent maps acquired by MODIS. All Landsat data downloaded from USGS have been geometrically corrected (Level 1T, L1T). The data collected by the Landsat TM and OLI included the surface reflectance generated by the USGS with the Landsat Ecosystem Disturbance Adaptive Processing System algorithm [64] and the Landsat Surface Reflectance Code algorithm [65], respectively.
Precipitation and temperature data were used to examine the potential natural factors that could influence the long-term spatial and temporal dynamics of lakes in the EPL. Long-term annual precipitation data during 2000-2020 were obtained from the Global Precipitation Measurement Integrated Multi-satellite Retrievals (GPM IMERG) Level-3 "Final" product data, which were downloaded from the NASA Goddard Earth Sciences (GES) Data and Information Services Center (DISC) (https://gpm1.gesdisc.eosdis.nasa. gov/data/GPM_L3/GPM_3IMERGM.06/, accessed on 24 March 2021). These products are produced with the IMERG algorithm, which were combined various satellite data [66,67]. With a spatial resolution of 0.1 • × 0.1 • (roughly~10 × 10 km), monthly GPM_3IMERGM data have been shown to perform slightly better than the Tropical Rainfall Measuring Mission data (TRMM 3B43) and are consistent with meteorological field data in terms of long-term trends [68][69][70]. The precipitation conditions of a lake were represented by the total precipitation over each lakes' water area, which was derived from the GPM_3IMERGM data. Annual temperature and wind data were obtained from the National Meteorological Information Center, China (https://data.cma.cn, accessed on 24 March 2021). Similarly, data collected from the closest meteorological station to a lake were used to represent the temperature conditions of each lake.
Several statistical variables were used to account for the potential human activities that might influence the long-term variability of the annual dynamics of lakes in the EPL. These variables included local population (LPO), gross domestic production (GDP), consumption of chemical fertilizers (CCF), irrigation area of cultivated land (IAC), and gross output value of farming, forestry, animal husbandry, and fishery (GFF) [28,57,71]. In addition, only the human activity factors specific to the city or prefecture in which a particular lake was located were used in the analysis.

Methodology
Although extracting water data is easy for a small remote sensing image or for a short time period, e.g., using thresholds for NDVI, NDWI, or MNDWI, the task is quite challenging for 184 different water bodies over 21 years. Through traditional 3D data storage formats, more than 7000 composites images over 21 years would need to be processed and analyzed, requiring enormous computer performance and many repeated steps for processing [72]. Therefore, we performed the necessary long-term time series lake mapping using the MDD structure that stores all band data by time sequence and each spectral cube in the band sequence. First, the long-term time series MODIS images were organized in the MDD structure [56,72]. Second, automated water extraction was conducted, and the final water maps were acquired through visual inspection. Lastly, patterns of the long-term dynamics of the lakes were analyzed. This framework is presented in Figure 2, and the details of each step are described below.

The Multi-Dimensional Dataset: TSB
As one of the MDD formats, the temporal sequential in band (TSB) data structure ( Figure 3) was used to store the long-term MODIS images. A multi-dimensional analysis of remote sensing (MARS) [72] was used to build the TSB. The software can be downloaded from: http://www.geodoi.ac.cn/WebCn/doi.aspx?Id=702 (accessed on 12 June 2021).

The Multi-Dimensional Dataset: TSB
As one of the MDD formats, the temporal sequential in band (TSB) data structure ( Figure 3) was used to store the long-term MODIS images. A multi-dimensional analysis of remote sensing (MARS) [72] was used to build the TSB. The software can be downloaded from: http://www.geodoi.ac.cn/WebCn/doi.aspx?Id=702 (accessed on 12 June 2021).

The Multi-Dimensional Dataset: TSB
As one of the MDD formats, the temporal sequential in band (TSB) data structure ( Figure 3) was used to store the long-term MODIS images. A multi-dimensional analysis of remote sensing (MARS) [72] was used to build the TSB. The software can be downloaded from: http://www.geodoi.ac.cn/WebCn/doi.aspx?Id=702 (accessed on 12 June 2021). As shown in Figure 3, the storage principles of the TSB are that the format stores the first band B1 at time T1 in the order of row-column, the pixel data of each band of T1 are stored sequentially, and all temporal cube data are stored through the above rules. This structure is defined as: where, represents the cube dataset for each time, which can be expressed as: As shown in Figure 3, the storage principles of the TSB are that the format stores the first band B 1 at time T 1 in the order of row-column, the pixel data of each band of T 1 are stored sequentially, and all temporal cube data are stored through the above rules. This structure is defined as: where, F TSB (t) represents the cube dataset for each time, which can be expressed as: Each band is stored in line-by-line order, which can be expressed as: where, t, s, c, and r represent time, spectral band, column, and row, respectively, and the ranges are [1, T], [1, S], [1, C], and [1, R], respectively. A(t, s, c, r) represents the spectral cube of A, the s band, c-th row and the r-line pixel of time t. B TSB represent the data of one band, F TSB represents multiple bands, and P TSB represents a temporal-spatial-spectral dataset.

Water Extraction Algorithms
Zhang et al. [73,74] proposed a sensor-independent remote sensing decomposition that utilizes a full spectrum information algorithm-universal pattern decomposition method (UPDM). This method is based on the fact that the spectral reflectance of each pixel can be decomposed into the standard spectral linear combination of water, vegetation, soil, and yellow leaves. The UPDM decomposition is: where C w , C v , C s , and C 4 are the decomposition coefficients for water, vegetation, soil, and yellow leaves, respectively; and P w , P v , P s , and P 4 are the standard spectral models for sensor band i of water, vegetation, soil, and yellow leaves, respectively. Through this projtection, each pixels of the multi-dimensional remote sensing image were converted into four-dimensional UPDM space. That is, four feature vectors (C w , C v , C s and C 4 ) were extracted through this convertion. Zhang et al. [75] optimized the UPDM decomposition calculation using the following formula: where R represents the reflectance value column vector of the pixel, and n is the value of bands number for remote sensing sensor. For example, the value of n is 7 for MODIS Acqua/Terra sensor. M is the value upon which the normalized values of the continuous spectra of the standard water, vegetation, soil, and yellow leaves are reconstructed into the reflectance values of each band based on the spectral response function of the target sensor. The value of M for the MODIS band can be found in Zhang et al. [75]. C represents the UPDM dimensional parameters of water, vegetation, soil, and yellow leaves. The pure pixels of lake, vegetation, dark mountains, and buildings in the EPL were selected, and the UPDM dimensional space distributions in the EPL are shown in Figure 4. In the UPDM dimensional space, the water pixels exhibit several notable characteristics.
(1) First, the values of C w and C 4 were higher for water pixels than for other objects ( Figure 4a,d), while the value of C w was significantly greater than 0. (2) In addition, C v < 0.04 and C s < 0.08 (p < 0.05), and the values of C v and C s were lower. Based on these characteristics of the UPDM dimensional space, the UNWI was defined as Equation (6): where t represents the time during 2000-2020. Ideally, when the value of C w or C 4 is the largest among C v , C s , and C w > 0, the pixels are considered water. Unfortunately, the environment of inland lakes is very complex, and if a lake is shallow or turbid, the influence of sediment or soil will cause the value of C s to be high. If a lake contains floating algae, the value of C v could also be high. To reduce the influence of sediment, soil or floating algae and enhance the separability of water and land pixels, the UNWI is considered as expressed in Equation (6) by stabilizing the threshold to distinguish water from non-water pixels. The equation results in large negative values for pixels covered by vegetation, dark mountains, or buildings, which have large values of C v and C s (see Figure 5), indicating that 0 could be used as a reasonable threshold to distinguish water from non-water under various environmental conditions. In conclusion, the water pixels were identified as C w > 0, C v < 0.04, C s < 0.08, and UNW I > 0. Using the UNWI method, automated lake maps were acquired, and visual inspection was used to avoid the effects of complex landcover types, mixed MODIS pixels, and noise in the EPL. The final water maps were obtained using the combination of UNWI and visual examination. Cloud-free Landsat images from the same time were used for verification. The accuracy of this method is discussed below.
from non-water under various environmental conditions. In conclusion, the water pixels were identified as 0, < 0.04, < 0.08, and 0 . Using the UNWI method, automated lake maps were acquired, and visual inspection was used to avoid the effects of complex landcover types, mixed MODIS pixels, and noise in the EPL. The final water maps were obtained using the combination of UNWI and visual examination. Cloud-free Landsat images from the same time were used for verification. The accuracy of this method is discussed below.  . Each box-plot shows the location of the maximum, 75th, 50th, 25th percentiles, and the minimum values using horizontal lines. (note, C w , C v , C s , and C 4 are the decomposition coefficients for water, vegetation, soil, and yellow leaves, respectively).

Calculation of Spatial and Temporal Statistics
To reduce the influence of clouds, aerosols, and image noise, the seasonal and annual water bodies were calculated using water body frequencies during the period of 2000-2020, which were calculated using Equation (7). Thresholds of 2/3 and 1/3 were used to identify annual and seasonal water bodies [8,15].
Two methods were used to test whether the changes in lakes area significantly increased or decreased in the long-term: (1) the Mann-Kendall (M-K) trend test [76] was used to determine the pattern distribution of annual lake dynamics in the EPL during 2000-2020; positive and negative z-statistics indicate increasing or decreasing trends, respectively. (2) The second method used linear regression to calculate R 2 values. Positive slopes indicate increasing patterns, while negative slopes indicate decreasing trends. Statistical significance of R 2 was set at p < 0.05.
where Fre water is the frequency of water bodies, ∑ N water is the total number of water bodies, and ∑ N Goodobservation is the total good number during the observation time. Remote Sens. 2022, 14, x FOR PEER REVIEW 10 of 28 Figure 5. UNWI distributions of pure pixels of major land cover types.

Calculation of Spatial and Temporal Statistics
To reduce the influence of clouds, aerosols, and image noise, the seasonal and annual water bodies were calculated using water body frequencies during the period of 2000-2020, which were calculated using Equation (7). Thresholds of 2/3 and 1/3 were used to identify annual and seasonal water bodies [8,15].
Two methods were used to test whether the changes in lakes area significantly increased or decreased in the long-term: (1) the Mann-Kendall (M-K) trend test [76] was used to determine the pattern distribution of annual lake dynamics in the EPL during 2000-2020; positive and negative z-statistics indicate increasing or decreasing trends, respectively. (2) The second method used linear regression to calculate R 2 values. Positive slopes indicate increasing patterns, while negative slopes indicate decreasing trends. Statistical significance of R 2 was set at p < 0.05.
where is the frequency of water bodies, ∑ is the total number of water bodies, and ∑ is the total good number during the observation time. Table 1 presents the verification of water extraction maps using the UNWI method for three lakes: Poyang, Chaohu, and Hongze. The true water maps were acquired through visual discrimination. The values of user's accuracy were all above 90%, with the exception of the image acquired on 5 October 2007. The producer's accuracy was relatively high (>90%). Based on the UNWI method, the average values of the user's and producer's accuracy were 94.52% and 95.65%, respectively. Two pairs of cloud-free Landsat images of Poyang, Chaohu and Hongze were collected and compared to  Table 1 presents the verification of water extraction maps using the UNWI method for three lakes: Poyang, Chaohu, and Hongze. The true water maps were acquired through visual discrimination. The values of user's accuracy were all above 90%, with the exception of the image acquired on 5 October 2007. The producer's accuracy was relatively high (>90%). Based on the UNWI method, the average values of the user's and producer's accuracy were 94.52% and 95.65%, respectively. Two pairs of cloud-free Landsat images of Poyang, Chaohu and Hongze were collected and compared to MODIS daily reflectance images on the same days ( Figure 6). The boundaries were derived from Landsat and acquired through visual discrimination. For example, for the more complex boundary of the water-land interface (see Figure 6a,b), a visual inspection of Poyang indicated that the Landsat boundaries were more detailed than those acquired by MODIS, due to the higher resolution (30 m) of the former. In particular, Hongze (Figure 6c,d), the area derived from MODIS was slightly larger than the Landsat images. By contrast, for relatively simple water-land interface of Chaohu, the boundaries derived from MODIS and Landsat were very similar (Figure 6e,f). The average differences between MODIS and Landsat of each independent pair were 6.3%, 2.0%, and 3.3% for Poyang, Hongze, and Chaohu, respectively. Overall, the water-land interface distributions strongly overlapped and were very similar, indicating that using a threshold of 0 from the UNWI images to delineate the boundaries of the water-land interface was reasonable and acceptable through visual examination and comparison. The spatial resolution of MODIS has a certain influence on the accuracy of lake area extraction. Different level areas of lakes were selected to detect these impacts ( Table 2). The Dongping Lake (C001,~100 km 2 ), Shaobo Lake (C008,~50 km 2 ), Douhu Lake (B002,~5 km 2 ), and Dazonghu Lake (F003,~5 km 2 ) images were acquired on the same day both on MODIS data and Landsat data. The average differences between MODIS and Landsat for each independent pair were 5.69%, 8.84%, 9.27%, and 8.37%, respectively, indicating that MODIS image data affects area accuracy in some smaller lake extraction, but this effect is acceptable (<10%) in the long-term.  In addition, the lake inundation areas were also compared to the NDWI method. The range of NDWI optimal threshold ranges were −0.06 and 0.15 during 2015-2017 based on visual inspection. A total of 103 and 91 pairs of results acquired from MODIS Terra and Aqua are shown in Figure 7. In general, compared to the UNWI method, the R 2 values for Poyang Lake (A002), Taihu Lake (D001), and Chaohu Lake (D143) varied  In addition, the lake inundation areas were also compared to the NDWI method. The range of NDWI optimal threshold ranges were −0.06 and 0.15 during 2015-2017 based on visual inspection. A total of 103 and 91 pairs of results acquired from MODIS Terra and Aqua are shown in Figure 7. In general, compared to the UNWI method, the R 2 values for Poyang Lake (A002), Taihu Lake (D001), and Chaohu Lake (D143) varied from 0.956 to 0.997 for MODIS Terra/Aqua. In addition, the R 2 value for all lake regions (sample sizes of 103 and 91 for MODIS Terra and Aqua, respectively) was 0.998 (Figure 7a,b), indicating that the MODIS water maps acquired using the UNWI method were reasonable and accurate.

Spatial and Temporal Patterns of lake Inundation during the Long-Term Observations
Seasonal water frequency between 2000 and 2020 in the EPL were calculated using Equation (7), and frequency maps were used to assess the seasonal characteristics of each lake inundation area (Figure 8). These frequencies were classified into six levels: 0-0.1, 0.1-0.3, 0.3-0.5, 0.5-0.65, 0.65-0.84, and 0.84-1. The water area of each lake was summed for each meteorological season (March-May, June-August, September-November,and December-February) of each year and used to analyze patterns during 2000-2020. Significant seasonal variation could be observed in almost all lakes located in the upper reaches of the Huai River and the middle and lower reaches of the Yangtze River. In autumn (Figure 8c), the number of lakes with significantly decreased lake inundation area (p < 0.05) was highest, accounting for 34.24% (63/184) of lakes. By contrast, in summer (Figure 8b), 10.33% (19/184) of lakes experienced a significantly increasing trend (p < 0.05), accounting for the largest proportion. This phenomenon may relate to the climate properties of China, characterized by a monsoon climate, high temperature, and rain in summer, with cold and less rain in winter. In particular, the EPL zone is situated in a subtropical monsoon climate zone [77,78], and more precipitation in summer could result in a significantly increasing trend of lake areas.

Spatial and Temporal Patterns of lake Inundation during the Long-Term Observations
Seasonal water frequency between 2000 and 2020 in the EPL were calculated using Equation (7), and frequency maps were used to assess the seasonal characteristics of each lake inundation area (Figure 8). These frequencies were classified into six levels: 0-0.1, 0.1-0.3, 0.3-0.5, 0.5-0.65, 0.65-0.84, and 0.84-1. The water area of each lake was summed for each meteorological season (March-May, June-August, September-November, and December-February) of each year and used to analyze patterns during 2000-2020. Significant seasonal variation could be observed in almost all lakes located in the upper reaches of the Huai River and the middle and lower reaches of the Yangtze River. In autumn (Figure 8c), the number of lakes with significantly decreased lake inundation area (p < 0.05) was highest, accounting for 34.24% (63/184) of lakes. By contrast, in summer (Figure 8b), 10.33% (19/184) of lakes experienced a significantly increasing trend (p < 0.05), accounting for the largest proportion. This phenomenon may relate to the climate properties of China, characterized by a monsoon climate, high temperature, and rain in summer, with cold and less rain in winter. In particular, the EPL zone is situated in a subtropical monsoon climate zone [77,78], and more precipitation in summer could result in a significantly increasing trend of lake areas.  The annual water frequency map is shown in Figure 1, and the annual lake inundation area also exhibited significant increases or decreases (p < 0.05). The annual water frequencies were also classified into six levels: 0-0.1, 0.1-0.3, 0.3-0.5, 0.5-0.65, 0.65-0.84, and 0.84-1, accounting for 27.23%, 15.35%, 10.50%, 7.31%, 11.26%, and 28.36%, respectively. The center of each lake exhibited a higher frequency than the edge. In addition, 13.59% (25/184) of the lakes exhibited significantly decreasing trends during the observation period of 21 years, while 25% (46/184) of the lakes experienced a significantly increasing trend. As shown in Tables 3 and S1, the annual area of A002 (Poyang Lake) showed the largest rate of decrease (−29.09 km 2 year -, 1 r = −0.5, p < 0.05) during 2000-2020, among all 184 study lakes, while the lake D30 (Huangdahu) experienced the larg- Figure 8. Seasonal water frequency of the study lakes in the EPL. Seasonal lake area significantly (p < 0.05) increased or decreased, denoted with the symbol "↑" or "↓", respectively. Water frequencies were calculated for each meteorological season: (a) spring, (b) summer, (c) autumn, and (d) winter.
The annual water frequency map is shown in Figure 1, and the annual lake inundation area also exhibited significant increases or decreases (p < 0.05). The annual water frequencies were also classified into six levels: 0-0.1, 0.1-0.3, 0.3-0.5, 0.5-0.65, 0.65-0.84, and 0.84-1, accounting for 27.23%, 15.35%, 10.50%, 7.31%, 11.26%, and 28.36%, respectively. The center of each lake exhibited a higher frequency than the edge. In addition, 13.59% (25/184) of the lakes exhibited significantly decreasing trends during the observation period of 21 years, while 25% (46/184) of the lakes experienced a significantly increasing trend. As shown in Tables 3 and S1, the annual area of A002 (Poyang Lake) showed the largest rate of decrease (−29.09 km 2 year −1 , r = −0.5, p < 0.05) during 2000-2020, among all 184 study lakes, while the lake D30 (Huangdahu) experienced the largest rate of increase (5.94 km 2 year −1 , r = 0.49, p < 0.05). Among the lakes located downstream of Poyang Lake, 24.73% (23/93) and 11.83% (11/93) exhibited significant (p < 0.05) increases and decreases, accounting for 50% (23/46) and 44% (11/25) of all lakes with significantly (p < 0.05) increasing and decreasing trends, respectively. Table 3. Summary of information (number and location) and changes in lake annual area (average area and standard deviation during 2000-2020). Statistically significant correlations and contributions are annotated with * p < 0.05 and ** p < 0.01. The full table is shown in Table S1. The percentages of lakes of the various classes that exhibited significantly (p < 0.05) decreasing or increasing trends in seasonal inundation area are presented in Table 4. For the different geographical classes, almost all experienced a significantly (p < 0.05) decreasing trend in seasonal inundation area, with the exception of Class I, which showed a decreasing trend in autumn only. The maximum percentages of lakes that exhibited significantly (p < 0.05) decreasing trends in seasonal inundation area for spring, summer, autumn, and winter, were 25% (3/12), 33.33% (4/12), 42.86% (3/7), and 50% (6/12) for lakes in Classes V, V, II, and V, respectively. By contrast, the seasonal patterns for increases in autumn (2.78% of lakes) and winter (10.41%) were similar, and only Class V lakes exhibited significant increases in inundation area. Particularly in summer, the maximum proportion of lakes exhibiting increases occurred in Class VI, accounting for 37.5% (3/8). However, the inter-annual trends for lake inundation area significantly differed among the different classes (upper left of Figure 1). Class IV exhibited the highest number of lakes with significant changes in inter-annual inundation area, with 17 lakes exhibiting significant decreases (p < 0.05) and 43 lakes exhibiting significant increases (p < 0.05) during 2000-2020. Furthermore, the number of lakes with significantly decreasing or increasing trends was the same for Classes V and VI, accounting for two lakes and one lake in each class, respectively. However, the inundation area of Poyang Lake (Class I, A002) and Douhu Lake (Class II, B002) significantly decreased (p < 0.05), with rates of change of −30.54 km 2 year −1 and −0.11 km 2 year −1 , respectively (Figure 9a). At the same time, the inter-annual inundation of lakes C005 (Class III) and E007 (Class V) also showed a significant decreased while lakes D004 (Class IV) and F003 (Class VI) showed a significantly increased trend (Figure 9c). By contrast, no lakes from Classes I or II exhibited significant increases (p < 0.05) in the inundation area (upper left of Figure 1). More importantly, the maximum percentages of lakes with significantly decreasing or increasing trends occurred in Classes I/VI and IV, accounting for 25% (1/4) and 29.86% (43/144), respectively (Figure 9b).
showed a significant decreased while lakes D004 (Class IV) and F003 (Class VI) showed a significantly increased trend (Figure 9c). By contrast, no lakes from Classes I or II exhibited significant increases (p < 0.05) in the inundation area (upper left of Figure 1). More importantly, the maximum percentages of lakes with significantly decreasing or increasing trends occurred in Classes I/VI and IV, accounting for 25% (1/4) and 29.86% (43/144), respectively (Figure 9b).

Driving Forces
Changes in climate factors, such as temperature and precipitation, directly affect the lake inundation area [4,57]. Furthermore, human activities can also drastically impact lake ecosystems [29,51]. However, the climate changes and human activities that affected the inundation area in the EPL during 2000-2020 still remain unclear. To determine the potential factors that lead to variability in intra-annual inundation, correlations between annual lake inundation area and corresponding local climate factors (precipitation, maximum temperature, average temperature, minimum temperature, and wind) and human activities (LPO, GDP, CCF, IAC, and GFF) were calculated for each lake.
Precipitation, GDP, and GFF were the main driving factors affecting local lake area ( Figure 10). For the 46 EPL lakes that experienced significant increases in annual area, lake area was positively correlated with precipitation, GDP, and GFF. By contrast, for the 25 lakes that exhibited significant decreases in annual area, lake area was negatively correlated with GDP, CCE, IAC, and GFF. Moreover, both positive and negative correlations were found between annual area and GDP and CCF, primarily in the lakes of Class IV. In total, 30.43% (56/184) and 23.37% (43/184) of lakes exhibited positive correlations with precipitation and GDP (p < 0.05), respectively, while 13.04% (24/184) and 21.19% (39/184) of lakes exhibited negative correlations with GDP and CCF (p < 0.05), respectively, suggesting that human activities strongly affected water area in the EPL. In addition, precipitation was an important factor affecting lake inundation area in regions with relatively high rainfall. For example, 59 lakes in the EPL exhibited significant (p < 0.05) correlations with annual precipitation, and 81.35% (48/59) of these were in Class IV. Although both positive and negative correlations were observed between annual inundation area and yearly maximum temperature, average temperature and minimum temperature, the relationship between minimum temperature and area was weakest, suggesting that yearly minimum temperature has little effect on the annual area of lakes. Remote Sens. 2022, 14, x FOR PEER REVIEW 21 of 28 (a) (b) Figure 10. Correlation coefficients between annual lake area and their potential driving factors, from 2000 to 2020. The factors represented by rows from top to bottom are precipitation, maximum temperature, average temperature, and minimum temperature, wind, LPO, GDP, CCF, IAC, and GFF. (a) Lake A001-D072 and (b) lake D073-F008. Note: when |r| > 0.369, the correlation is marginally significant at p < 0.1; when |r| > 0.433, the correlation is significant at p < 0.05; and when |r| > 0.549, the correlation is significant at p < 0.01. When r > 0, it indicates a positive correlation; and when r < 0, it represents a negative correlation.

Uncertainties
In the present study, MOD09A1/MYD09A1 images were used to obtain water maps of 184 lakes in the EPL. The quality of MODIS surface reflectance products is one reason for the uncertainty of spatio-temporal monitoring of the dynamics of large lakes. On one hand, each pixel of MOD09A1/MYD09A1 surface reflectance data contains the best possible L2G observations during an eight-day period; thus, large errors, stripe noise, clouds, shadows, and negative reflectance can be avoided in part. On the other hand, QA data with these MODIS surface reflectance products were also used to detect good observation pixels and to exclude factors such as clouds, cloud shadows, and snow. At the same time, banding/gap is also an important factor affecting MODIS data, the strip data were not involved in calculation and analysis. Due to the low spatial resolution of MODIS image data, its accuracy is not very high in some small lakes, and the challenge of low spatial resolution of the MODIS sensor, we currently have few better methods to solve it. In the future, if there are long-sequence, high-frequency, high-spatial-resolution image data, higher-precision results will be obtained, and when the sentinel data accumulates a longer time series, the long sequence, high frequency, and high spatial resolution image data will obtain higher precision results.
Despite the fact that TRMM 3B42 data are critical sources of long-term precipitation products for analyzing the relationships among precipitation, water budget, and inundation area [3][4][5][79][80][81], the spatial resolution of TRMM 3B42 is 0.25° × 0.25°, which is not very suitable for inland lakes. Moreover, this system stopped working and receiving data (a) Lake A001-D072 and (b) lake D073-F008. Note: when |r| > 0.369, the correlation is marginally significant at p < 0.1; when |r| > 0.433, the correlation is significant at p < 0.05; and when |r| > 0.549, the correlation is significant at p < 0.01. When r > 0, it indicates a positive correlation; and when r < 0, it represents a negative correlation.
The most important climate and human factors affecting lake dynamics differed across the geographical classes. For example, CCF was the key driving factor for Classes I and II, for which 75% (3/4) and 42.86% (3/7) of lakes exhibited significant (p < 0.05) correlations between inundation area and CCF, respectively. One noticeable phenomenon was that two lakes in both Classes I and II exhibited significant (p < 0.05) negative correlations with CCF, indicating that human agricultural activities were the main factor causing changes in area for lakes that are freely connected to the Yangtze and Huai Rivers. By contrast, precipitation was the key driving factor affecting changes in area for lakes freely connected to the Beijing-Hangzhou Grand Reach; 33.33% (3/9) of these lakes exhibited significantly (p < 0.05) positive correlations with precipitation in Class III. For Class IV lakes, GDP, precipitation, and CCF were the three most dominant driving factors of changes in lake area, with 38.89% (56/144), 33.33% (48/144), and 32.64% (47/144) of lakes exhibiting significant correlations with GDP, precipitation, and CCF, respectively. In addition, CCF and GDP were the key factors affecting changes in lake area for Classes V and VI, with 41.67% (5/12) and 37.5% (3/8) of lakes, respectively, showing significant correlations with CCF and GDP.
In addition, the construction of the Three Gorges Dam (TGD) may be a primary reason why the area of Poyang Lake (Class I) decreased. The TGD was impounded in 2003, leading to sharply decreased water flow and in turn decreased water area (Figure 9a). Previous studies have demonstrated that the construction of the TGD has led to reductions in the area of Poyang Lake since 2003 [2][3][4].

Uncertainties
In the present study, MOD09A1/MYD09A1 images were used to obtain water maps of 184 lakes in the EPL. The quality of MODIS surface reflectance products is one reason for the uncertainty of spatio-temporal monitoring of the dynamics of large lakes. On one hand, each pixel of MOD09A1/MYD09A1 surface reflectance data contains the best possible L2G observations during an eight-day period; thus, large errors, stripe noise, clouds, shadows, and negative reflectance can be avoided in part. On the other hand, QA data with these MODIS surface reflectance products were also used to detect good observation pixels and to exclude factors such as clouds, cloud shadows, and snow. At the same time, banding/gap is also an important factor affecting MODIS data, the strip data were not involved in calculation and analysis. Due to the low spatial resolution of MODIS image data, its accuracy is not very high in some small lakes, and the challenge of low spatial resolution of the MODIS sensor, we currently have few better methods to solve it. In the future, if there are long-sequence, high-frequency, high-spatial-resolution image data, higher-precision results will be obtained, and when the sentinel data accumulates a longer time series, the long sequence, high frequency, and high spatial resolution image data will obtain higher precision results.
Despite the fact that TRMM 3B42 data are critical sources of long-term precipitation products for analyzing the relationships among precipitation, water budget, and inundation area [3][4][5][79][80][81], the spatial resolution of TRMM 3B42 is 0.25 • × 0.25 • , which is not very suitable for inland lakes. Moreover, this system stopped working and receiving data on 16 April 2015. Fortunately, GPM, TRMM's successor, was launched in February 2014 to provide global-scale precipitation products [69,82,83]. Based on the IMERG algorithm, all microwave satellite sensors and microwave-calibrated infrared (IR) satellites are crosscalibrated, merged and interpolated to generate GPM IMERG products with a spatial resolution of 0.1 • × 0.1 • [66,68,84]. The GPM IMERG products perform better in southern and eastern China [69,85,86]. The accuracy of GPM IMERG products in the EPL also affects the spatial distribution of lake waters. In the present study, annual precipitation data were acquired from monthly GPM IMERG products and then compared to annual precipitation from field meteorological stations to evaluate the precipitation precision of GPM. The correlation coefficients for stations in Yanzhou, Wuxi, Yueyang, and Poyang were calculated ( Figure 11) and were all greater than 0.8, with the exception of Poyang station. The correlations for all data pairs were 0.90, indicating good consistency between annual GPM and meteorological station precipitation. However, a correction algorithm of GPM IMERGE products should be developed in conjunction with meteorological station data in China. At the same time, water level and runoff are other important factors affecting water area, we will analyze the relationship between water level, runoff, and area in a future study.

Implications
The UNWI method is a robust water extraction index, due to its use of full spectrum information, and it can also be applied to Landsat images ( Figure 12). The Landsat OLI images were acquired on 15 April 2020 and 1 October 2014 in Poyang ( Figure 12a) and Hongze (Figure 12c), respectively. The red color in the UNWI images produced good fits with water from the Landsat RGB images by visual examination. The water inundation of Poyang Lake during the dry season acquired by Landsat TM and OLI images via the UNWI method were also compared to the results of Wang et al. [29]. As shown in Figure 13a, the inundation area of Poyang Lake (A002) obtained through UNWI agrees well with the results of the RF classification [29]. Specifically, the R 2 value for all data combined was 0.97, and the slope was 1.0096. These results indicate that water area data obtained using the UNWI method are consistent with the results of the RF classification method. Indeed, Remote Sens. 2022, 14, 985 20 of 26 the UNWI method also works well for Sentinel-2 images, as shown in Figure 13b. The R 2 values for data obtained using the UNWI method and the true inundation area of Nanyi (D006), Shijiu (A001), and Gucheng (D004) Lakes acquired by visual discrimination were 0.99, 0.99, and 0.91, respectively. The R 2 value for all lake regions (84 lakes) was 0.999, and the slope was 1.0033 (Figure 13b), suggesting that the UNWI method can be used successfully for extracting water area based on Landsat and Sentinel-2 images. tial resolution of 0.1° × 0.1° [66,68,84]. The GPM IMERG products perform better in southern and eastern China [69,85,86]. The accuracy of GPM IMERG products in the EPL also affects the spatial distribution of lake waters. In the present study, annual precipitation data were acquired from monthly GPM IMERG products and then compared to annual precipitation from field meteorological stations to evaluate the precipitation precision of GPM. The correlation coefficients for stations in Yanzhou, Wuxi, Yueyang, and Poyang were calculated ( Figure 11) and were all greater than 0.8, with the exception of Poyang station. The correlations for all data pairs were 0.90, indicating good consistency between annual GPM and meteorological station precipitation. However, a correction algorithm of GPM IMERGE products should be developed in conjunction with meteorological station data in China. At the same time, water level and runoff are other important factors affecting water area, we will analyze the relationship between water level, runoff, and area in a future study. Figure 11. Comparisons of meteorological annual precipitation and GPM IMERG precipitation (note, the gray line represents the 1:1 line).

Implications
The UNWI method is a robust water extraction index, due to its use of full spectrum information, and it can also be applied to Landsat images ( Figure 12). The Landsat OLI images were acquired on 15 April 2020 and 1 October 2014 in Poyang ( Figure 12a) and Hongze (Figure 12c), respectively. The red color in the UNWI images produced good fits with water from the Landsat RGB images by visual examination. The water inundation of Poyang Lake during the dry season acquired by Landsat TM and OLI images via the UNWI method were also compared to the results of Wang et al. [29]. As shown in Figure 13a, the inundation area of Poyang Lake (A002) obtained through UNWI agrees well with the results of the RF classification [29]. Specifically, the R 2 value for all data combined was 0.97, and the slope was 1.0096. These results indicate that wa- The MDD format, developed by [87], is an open dataset, and includes five data storage structures: temporal sequential in band (TSB), temporal sequential in pixel (TSP), temporal interleaved by band (TIB), temporal interleaved by pixel (TIP), and temporal interleaved by spectrum (TIS). The TSB format was applied to the extraction of water area over the 21-year study period and was used to store and calculate~7000 composite MODIS images in this study. As part of the MDD format, the TSB format can avoid complex problems caused by traditional image storage when dealing with time series images, since this format has the distinct advantage of grouping all bands of the time series data together and integrating the four-dimensional features involved in remote sensing products, such as 1D-spectral dimension, 2D-spatial dimensions, and 1D-temporal dimension, into a 4-D cube. The five data structures of MDD can be converted into one other according to MARS. Due to differences in physical storage, the different structures have different application scenarios. The TSB format is suitable for band operation and spatial filtering, since it stores spectral cube data by time sequence, and each spectral cube is stored in band sequence. The TSP format is appropriate for extracting spectral data of one or more pixels and for spectral derivative and integral operations. The advantages of the TIB, TIP, and TIS formats include extracting time series data from a single band, processing the time-spectrum data of pixels, and extracting time series data of the pixel spectral curve, respectively. classification method. Indeed, the UNWI method also works well for Sentinel-2 images, as shown in Figure 13b. The R 2 values for data obtained using the UNWI method and the true inundation area of Nanyi (D006), Shijiu (A001), and Gucheng (D004) Lakes acquired by visual discrimination were 0.99, 0.99, and 0.91, respectively. The R 2 value for all lake regions (84 lakes) was 0.999, and the slope was 1.0033 (Figure 13b), suggesting that the UNWI method can be used successfully for extracting water area based on Landsat and Sentinel-2 images.  as shown in Figure 13b. The R values for data obtained using the UNWI method and the true inundation area of Nanyi (D006), Shijiu (A001), and Gucheng (D004) Lakes acquired by visual discrimination were 0.99, 0.99, and 0.91, respectively. The R 2 value for all lake regions (84 lakes) was 0.999, and the slope was 1.0033 (Figure 13b), suggesting that the UNWI method can be used successfully for extracting water area based on Landsat and Sentinel-2 images.   The 21-year inundation frequencies and records of 184 lakes in the EPL are highly valuable for the long-term monitoring and protection of water and wetland ecosystems. Characterized by high levels of urbanization and population density, the EPL region also serves as the most important freshwater source for human and commercial use in eastern China. Therefore, the spatio-temporal monitoring of the dynamics of large lakes in the EPL is crucial to the national economy and peoples' livelihoods. Time series data for lake inundation from 2000 to 2020 provide not only a dynamic database for lake area but also a foundation for water quality monitoring. Based on our findings, lake reclamation and urbanization in the EPL can be quantified in future studies to achieve a better understanding of human-environment issues.

Conclusions
The present study provides a new method for long-term remote sensing monitoring of patterns in spatio-temporal dynamics of water; this method avoids the need for many repeated steps during processing, and it develops a full-band water extraction index, which uses full spectrum information to avoid or reduce optimal threshold adjustments. Compared to Landsat images, the producer's accuracy range of water extraction based on UNWI algorithms was 90.45-98.87%, and the average R 2 value of the lake inundation area obtained using the UNWI versus the NDWI method was 0.97. Based on the MDD framework and UNWI, the inundation conditions of 184 lakes in the EPL during 2000-2020 were accurately mapped through more than 7000 composite eight-day MOD09A1/MYD09A1 images.
Based on our results, significant seasonal variation was observed in almost all lakes located in the upper reaches of the Huai River and the middle and lower reaches of the Yangtze River. The number of lakes exhibiting significant (p < 0.05) decreases or increases was highest in autumn/summer. Meanwhile, during the 21-year observation period, 13.58% (25/184) and 25% (46/184) of lakes significantly decreased and increased in annual inundation area, respectively. Simultaneously, the inter-annual inundation area of 17 Class IV lakes exhibited a significantly decreasing trend (p < 0.05) and 43 Class IV lakes showed a significantly increasing trend (p < 0.05) from 2000-2020, which were the highest percentages for all six geographical classes.
The analysis of climate factors and human activities indicated that precipitation, GDP, and GFF were the main driving factors that led to significant increases in annual area for 46 lakes. CCF was the key driving factor for lakes in Classes I and II. By contrast, precipitation was the key driving factor affecting changes of area in Class III lakes. For Class IV lakes, GDP was the dominant driving factor of changes in lake area, with 38.89% (56/144) exhibiting significant changes. GDP and CCF were also considered key factors affecting changes in area for the lakes of Classes V and VI.
Overall, our findings provide accurate estimates of the spatio-temporal dynamics of large lakes in the EPL, China, between 2000-2020. Seasonal and annual patterns of the lake inundation area were quantified, and potential driving factors of changes in lake area were documented. Our results provide the basis and necessary foundation for the development of water quality monitoring and continuous investigations in the EPL.
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/rs14040985/s1, Table S1: Summary of information (number and location), lake area changes (average area and standard deviation during 2000 ± 2020). Statistically significant correlations and contributions are annotated with * and ** represents p < 0.01 and * represents p < 0.05.