Long-Term Surface Water Dynamics Analysis Based on Landsat Imagery and the Google Earth Engine Platform: A Case Study in the Middle Yangtze River Basin

: Dynamics of surface water is of great significance to understand the impacts of global changes and human activities on water resources. Remote sensing provides many advantages in monitoring surface water; however, in large scale, the efficiency of traditional remote sensing methods is extremely low because these methods consume a high amount of manpower, storage, and computing resources. In this paper


Introduction
Water resources are critical in promoting sustainable development, as they support human communities, maintain the functions of ecosystems, and ensure economic growth [1].Surface water is a key indicator of water resources.As a land cover type, it plays an important role in climate regulation, biogeochemical cycling, and surface energy balance, among many others [2].In recent decades, many countries, especially developing countries, have experienced rapid urbanization [3].Changes in surface water caused by human activities strongly affect surface temperature, soil moisture, biological diversity, ecosystem functioning, and even human wellbeing [4][5][6].Therefore, monitoring the dynamics of surface water is of great importance for natural environmental health and sustainable economic development [7].
Research in remote sensing of surface water started before 1970 [8,9].Studies then emerged in large numbers [10][11][12][13][14][15].Around 2000, accompanied by the rapid development of remote sensing satellites, several efficient water indexes were proposed to identify surface water coverage, such as the Normalized Difference Water Index (NDWI) [16] and the modified Normalized Difference Water Index (mNDWI) [17].A variety of passive and active remote sensors with visible and microwave bands have been used to estimate inundation area and delineate water boundaries [18], and such sensors include Moderate Resolution Imaging Spectrometer (MODIS) [7,19], the Landsat Thematic Mapper (TM) [2,6,11,[20][21][22][23][24][25], Synthetic Aperture Radar (SAR) [15,26,27], SPOT [28,29], and IKONOS [30].Up to now, remote sensing has helped acquire valuable information by providing huge amounts of images that cover the Earth's surface over a period of 40 years [31].However, traditionally, data acquisition and storage, obscure file formats, and multitudes of geospatial data processing frameworks are significant obstacles to take full advantage of these images, especially in large-scale and long-term applications [32].Recently, a free cloud computing platform called the Google Earth Engine (GEE) [33] has been used to store and process large volumes of remote sensing images.The GEE provides a programming and graphic interface to a copy of remote sensing images, as well as the power of dedicated cloud storage and computational hardware [34].It has been used in large-scale applications, including mapping vegetation cover [35][36][37], settlement and population [34,38], and the detection of boundaries of urban areas [3].The GEE also has advantages in terms of exploring the long-term dynamics of surface water in large areas.Global seasonal and permanent surface water extent was calculated by recording the months and years when water was present between 1984 and 2015, and the results are available in the GEE platform, named as Joint Research Centre (JRC) Yearly Water Classification History (v1.0) [6].Earth's surface water change over the past 30 years was also analyzed through the Deltares Aqua Monitor [39] with assistance from Google Earth Engine [40].However, the automatic method for a global scale lacks suitability for a regional scale.The diversity of regional features will lead to inconsistency in the accuracy of the results.Related research about the dynamics of regional surface water and analyses of their driving forces are rare.
The Yangtze River is the longest river in China and the third longest river in the world, it plays a key role in the social and economic development of China [41].The Middle Yangtze River Basin (MYRB) is located in the middle of China.Appropriate climate and geographical conditions have made the MYRB an important food production area in China.With "the Development Plan for the City Cluster along the Middle Yangtze River Basin" proposed as a national strategy [42], the MYRB will become an important economic growth pole for China.The MYRB contains a complex system of surface water, including marshes, multiple river channels, and thousands of lakes [43].Over the last century, the environment in the MYRB has been changed a lot due to increased human activities [43].The MYRB has a large population and rapid economic development, meanwhile frequent floods and changeable hydrological conditions threaten the economy and ecology.Therefore, accurate and rapid monitoring of the long-term dynamics of the surface water in the MYRB is important to the regional economy and ecosystem [41].
Thus, the objectives of this study are (1) to propose a new method of quickly estimating annual surface water, using the existing mature surface water index, (2) to apply the proposed method to the variation in surface water in the MYRB from 1990 to 2017, based on the GEE platform, and (3) to analyze the causes of surface water variation in the MYRB.The approach will provide a new perspective on surface water research, and the results of this study will assist in the remote-sensing-based evaluation of water resource distribution and dynamics.Furthermore, this study also facilitates the development of appropriate water resource conservation and management strategies.

Study Area
The MYRB is located in the center of the whole Yangtze River Basin (Figure 1a) and covers an area of about 680,000 km 2 [44].It is composed of four sub-basins: the Hanjiang River Section (HjRS), Middle Mainstream Section (MMS), Dongting Lake Section (DtLS), and the Poyang Lake Section (PyLS) [45][46][47] (Figure 1b).The MYRB experiences a humid continental climate, with humid summers and cold winters [48].The altitude ranges from −146 to 3562 m, with an average value of 627 m (Figure 1b) [49].There are a number of lakes and reservoirs, including the two largest freshwater lakes in China: Dongting Lake and Poyang Lake [46].The region is rich in biodiversity with 11,770 km 2 nature reserves [50], over 400 species of fish and over 7000 species of vegetation.The large population and rapid economic development in the MYRB have led to drastic changes in surface water.

Study Area
The MYRB is located in the center of the whole Yangtze River Basin (Figure 1a) and covers an area of about 680,000 km 2 [44].It is composed of four sub-basins: the Hanjiang River Section (HjRS), Middle Mainstream Section (MMS), Dongting Lake Section (DtLS), and the Poyang Lake Section (PyLS) [45][46][47] (Figure 1b).The MYRB experiences a humid continental climate, with humid summers and cold winters [48].The altitude ranges from −146 to 3562 m, with an average value of 627 m (Figure 1b) [49].There are a number of lakes and reservoirs, including the two largest freshwater lakes in China: Dongting Lake and Poyang Lake [46].The region is rich in biodiversity with 11,770 km 2 nature reserves [50], over 400 species of fish and over 7000 species of vegetation.The large population and rapid economic development in the MYRB have led to drastic changes in surface water.

Data Preparation
Surface water changed frequently within every year.To acquire more information about its extent, we used a dense time series of Landsat images of the study area.Landsat images with a spatial resolution of 30 m are suitable for monitoring regional land covers.Complete coverage of the study area is achieved with 55 tiles of Landsat Worldwide Reference System (WRS) paths/rows.In this study, a total of 2343 United States Geological Survey (USGS) Landsat calibrated top-of-atmosphere (TOA) reflectance, orthorectified images were acquired for the years of 1990, 2000, 2010, and 2017, respectively, and were available on the GEE as image collections.Detailed statistical information about the Landsat Thematic Mapper (TM), and Operational Land Imager (OLI) data used in this study is given in Table 1.The spatial distribution of total observation counts in the year of 2010 in the study area is presented in Figure 3A.However, mean annual cloud frequencies of the study area over 2000-2014 were between 56.49% and 89.81%, estimated through the Global 1-km Cloud Cover [51].Such a high probability of cloud coverage has significant impacts on image interpretation.As a result, it is difficult to generate a wall-to-wall map of the study area.To get the utmost out of all the acquired images, the index "cloud score" was calculated for every pixel in every image by the built-in algorithm

Data Preparation
Surface water changed frequently within every year.To acquire more information about its extent, we used a dense time series of Landsat images of the study area.Landsat images with a spatial resolution of 30 m are suitable for monitoring regional land covers.Complete coverage of the study area is achieved with 55 tiles of Landsat Worldwide Reference System (WRS) paths/rows.In this study, a total of 2343 United States Geological Survey (USGS) Landsat calibrated top-of-atmosphere (TOA) reflectance, orthorectified images were acquired for the years of 1990, 2000, 2010, and 2017, respectively, and were available on the GEE as image collections.Detailed statistical information about the Landsat Thematic Mapper (TM), and Operational Land Imager (OLI) data used in this study is given in Table 1.The spatial distribution of total observation counts in the year of 2010 in the study area is presented in Figure 3A.However, mean annual cloud frequencies of the study area over 2000-2014 were between 56.49% and 89.81%, estimated through the Global 1-km Cloud Cover [51].Such a high probability of cloud coverage has significant impacts on image interpretation.As a result, it is difficult to generate a wall-to-wall map of the study area.To get the utmost out of all the acquired images, the index "cloud score" was calculated for every pixel in every image by the built-in algorithm "simpleCloudScore" of the GEE platform [52].The source code of the "simpleCloudScore" can be found at [53].Pixels with a low cloud score index, called "good observation", were then used.The "simpleCloudScore" algorithm determines a cloud-likelihood score in the range [0, 100] by using a combination of the Normalized Difference Snow Index (NDSI) and the brightness and temperature from the Landsat TOA reflectance imagery [54].The simpleCloudScore algorithm cannot be used as a robust cloud detector, and its main purpose is to compare multiple appearances at the same point for relative cloud likelihood [54].For this study, a cloud score threshold of 20 was used based on the visual interpretation of Landsat images [55] to give adequate performance at detecting clouds, and pixels with cloud scores less than 20 were defined as "good observation", and pixels showing a cloud score higher than 20 were masked out.The GEE Javascript codes for processing the Landsat image collections are available as Supplementary Materials.Pixels in good observation images account for 51.10% of all pixels in total observation images at the study area in 2010.The spatial distribution of good observation counts in the year of 2010 in the study area is presented in Figure 3B.A stratified random sampling strategy was utilized to select reference samples.For the periods of 2000, 2010, and 2017, we collected reference data from high-resolution images available on Google Earth.For the period of 1990, we selected samples from historic topographic maps available in the library of Wuhan University.Because the focus of this study is the spatio-temporal dynamics of surface water, the reference samples consisted of two categories (water and non-water) and four sub-categories (water, vegetation, built-up and bare areas, and shadow).Detailed information of reference samples is given in Table 2.Among these points, 70% of them were gathered as training samples during classification.The remaining points were used to verify the accuracy of the classification results.

Methodology
Since there is great variation in surface water within a given year, in order to obtain abundant cover information about surface water, the minimal surface water and maximal surface water in each period were extracted.In this study, the minimal surface water means permanent surface water, and the maximal surface water means all the seasonal and permanent surface water.The yearly minimal and maximal surface water were generated by creating the greenest and wettest pixel composite images respectively for each period, the details are described in the following section.A flowchart of extracting the surface water extent in this study is presented in Figure 4.It mainly consists of two parts: (1) removing cloudy pixels with the index "cloud score", then adding the Normalized Difference Vegetation Index (NDVI), NDWI, and mNDWI into the bands of images, and removing obvious non-water pixels by NDVI masking, and NDWI masking procedures, and (2) further extracting surface water pixels, including image classification and manual noise removal.A detailed description of these procedures will be provided in the following sub-sections.

Removing Obvious Non-Water Pixels by Normalized Difference Indexes
Three widely used vegetation indices-the NDVI [56], the NDWI [16], and mNDWI [17]-were calculated for each imagery and added to the image collections, using Equations ( 1 where ρ N IR , respects the values of near infrared band, ρ red respects the values of red band, ρ green respects the values of green band, and ρ SW IR respects the values of shortwave infrared band in the Landsat Images.The internal function "simpleCloudScore" implemented in the GEE was used for cloud masking, as described in Section 2.2.The qualityMosaic method, embedded in the GEE, sets each pixel in the composite based on which image in the collection has a maximum value for the specified band.In this study, the NDVI band was selected to create a qualityMosaic image that was a composite of the "greenest" pixels and called greenest image [57].The greenest image depicted the greatest extent Remote Sens. 2018, 10, 1635 7 of 18 of vegetation observed in the studied periods with minimal surface water extent [58].The NDWI band was selected to create a qualityMosaic image which was composite of the "wettest" pixels and called wettest image.The wettest image depicting the maximal possible surface water in the studied periods.The greenest and wettest images of 2017 in the DtLS sub-basin were compared and are shown in Figure 5.   Since the water pixels have relatively lower NDVI values and higher NDWI values, we can acquire a preliminary surface water map by masking out all pixels that have higher NDVI values and lower NDWI values.In this study, the threshold value of NDVI greater than or equal to 0.2 and NDWI less than or equal to −0.3 were set to remove most non-water pixels, then improve the efficiency and accuracy of the classification algorithm.The threshold was established through checking the effects of repeated experiments with different values.In the GEE, the procedure was easily accomplished by the code "updateMask(ndvi.lt(0.2) and (ndwi.gte(−0.3)))"("ndvi" and "ndwi" are new bands added to the image by the former step).Through this operation, pixels with obvious vegetation or non-water were removed.
However, it was found that the NDWI cannot efficiently suppress the signal from built-up surfaces and shadows, and using the rule (NDVI >= 0.2 and NDWI <= −0.3) does not accurately enable discriminating parts of built-up surfaces and shadows from water pixels [59,60], so the RF classification model was used to obtain more accurate surface water coverage in the next step.

Further Extracting Surface Water Pixels by Random Forest (RF)
After masking with the NDVI and the NDWI, the rest of the imagery can be divided into three classes: (1) surface water, (2) built-up areas, which also includes bare ground, roads, and building facilities, and (3) shadows of mountains, buildings, and clouds.In order to further extract the surface water, the RF model was chosen to execute classification.The RF model is an ensemble, nonparametric modeling approach [38] that grows a "forest" of individual classification or regression trees and improves upon bagging [61] by using the best of a random selection of predictors at each node in each tree [62,63].
In the random forest model, the number of classification trees and the number of features at each node for spilling are vital to the results [64].RF classifiers with 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 15, 20, 25, and 30 classification trees were applied, and the optimal item was chosen based on the overall accuracy and Kappa accuracy.The overall accuracy implies the overall degree of agreement in the error matrix [65].The Kappa coefficient expresses how much better the classification is compared with a random classification [66].As shown in Figure 6, 20 is the optimal number of classification trees in our experiment, for classification result gains a slight improvement in performance when the number is set above 20.Moreover, it has been suggested in previous works that the number of features of classification trees should be the square root of the total number of features [67,68].In this study, the three calculated normalized difference indices, combined with the original spectral bands, were used to improve confidence [69].Hence, we set the number of features as 4. Since the water pixels have relatively lower NDVI values and higher NDWI values, we can acquire a preliminary surface water map by masking out all pixels that have higher NDVI values and lower NDWI values.In this study, the threshold value of NDVI greater than or equal to 0.2 and NDWI less than or equal to −0.3 were set to remove most non-water pixels, then improve the efficiency and accuracy of the classification algorithm.The threshold was established through checking the effects of repeated experiments with different values.In the GEE, the procedure was easily accomplished by the code "updateMask(ndvi.lt(0.2) and (ndwi.gte(−0.3)))"("ndvi" and "ndwi" are new bands added to the image by the former step).Through this operation, pixels with obvious vegetation or non-water were removed.
However, it was found that the NDWI cannot efficiently suppress the signal from built-up surfaces and shadows, and using the rule (NDVI >= 0.2 and NDWI <= −0.3) does not accurately enable discriminating parts of built-up surfaces and shadows from water pixels [59,60], so the RF classification model was used to obtain more accurate surface water coverage in the next step.

Further Extracting Surface Water Pixels by Random Forest (RF)
After masking with the NDVI and the NDWI, the rest of the imagery can be divided into three classes: (1) surface water, (2) built-up areas, which also includes bare ground, roads, and building facilities, and (3) shadows of mountains, buildings, and clouds.In order to further extract the surface water, the RF model was chosen to execute classification.The RF model is an ensemble, nonparametric modeling approach [38] that grows a "forest" of individual classification or regression trees and improves upon bagging [61] by using the best of a random selection of predictors at each node in each tree [62,63].
In the random forest model, the number of classification trees and the number of features at each node for spilling are vital to the results [64].RF classifiers with 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 15, 20, 25, and 30 classification trees were applied, and the optimal item was chosen based on the overall accuracy and Kappa accuracy.The overall accuracy implies the overall degree of agreement in the error matrix [65].The Kappa coefficient expresses how much better the classification is compared with a random classification [66].As shown in Figure 6, 20 is the optimal number of classification trees in our experiment, for classification result gains a slight improvement in performance when the number is set above 20.Moreover, it has been suggested in previous works that the number of features of classification trees should be the square root of the total number of features [67,68].In this study, the three calculated normalized difference indices, combined with the original spectral bands, were used to improve confidence [69].Hence, we set the number of features as 4. To promote the visual interpretation, a false color composite of the TM Bands 7 (swir2), 5 (swir1), and 3 (green) and a composite of OLI Bands 7 (swir2), 6 (swir1), and 4 (green) were generated.The best band color combination for detecting surface water was dark blue.
In order to reduce noise in the classification result, we applied manual noise removal as implemented in ESRI ArcMap [70] in post-processing.Trivial cells with an area of less than 5400 m 2 were removed manually to improve the homogeneity of the classified raster.
Based on the masked with cloud score, the NDVI and the NDWI, the quality mosaicked with the NDVI or the NDWI, the RF classification, and the manual correction, yearly minimal and maximal surface water coverages were extracted for 1990, 2000, 2010, and 2017.

Accuracy Assessment
Considering that the minimal and maximal surface water was extracted by the same method and original image collections, we assumed that classification accuracies of these two maps are similar.Additionally, the 3424 validation samples described in Table 2 were all selected from permanent water surface.Therefore, in this study, only the accuracies of minimal surface water were validated.As shown in Table 3, error matrixes of the classification results in 1990, 2000, 2010, and 2017 were calculated, respectively.The surface water maps in these periods all had high accuracies with overall accuracies of 86-93%, specifically, 86% in 1990, 88% in 2000, 90% in 2010 and 93% in 2017.The surface water category had producer accuracies over 84% and user accuracies over 86%.Errors were mainly caused by misclassification of water surfaces with shadows.The results of the accuracy assessment show that surface water extraction results are consistent with the ground truth.To promote the visual interpretation, a false color composite of the TM Bands 7 (swir2), 5 (swir1), and 3 (green) and a composite of OLI Bands 7 (swir2), 6 (swir1), and 4 (green) were generated.The best band color combination for detecting surface water was dark blue.
In order to reduce noise in the classification result, we applied manual noise removal as implemented in ESRI ArcMap [70] in post-processing.Trivial cells with an area of less than 5400 m 2 were removed manually to improve the homogeneity of the classified raster.
Based on the masked with cloud score, the NDVI and the NDWI, the quality mosaicked with the NDVI or the NDWI, the RF classification, and the manual correction, yearly minimal and maximal surface water coverages were extracted for 1990, 2000, 2010, and 2017.

Accuracy Assessment
Considering that the minimal and maximal surface water was extracted by the same method and original image collections, we assumed that classification accuracies of these two maps are similar.Additionally, the 3424 validation samples described in Table 2 were all selected from permanent water surface.Therefore, in this study, only the accuracies of minimal surface water were validated.As shown in Table 3, error matrixes of the classification results in 1990, 2000, 2010, and 2017 were calculated, respectively.The surface water maps in these periods all had high accuracies with overall accuracies of 86-93%, specifically, 86% in 1990, 88% in 2000, 90% in 2010 and 93% in 2017.The surface water category had producer accuracies over 84% and user accuracies over 86%.Errors were mainly caused by misclassification of water surfaces with shadows.The results of the accuracy assessment show that surface water extraction results are consistent with the ground truth.

Area and Distrubution of Surface Water in the MYRB from 1990 to 2017
The study region encompasses over 680,000 km 2 divided into two categories: surface water and non-water.The minimal and maximal surface water extents for each of the studied years are given in Figures 7 and 8.The yearly minimal surface water extents in 1990, 2000, 2010, and 2017 are 14,751.23km 2 , 14,403.48 km 2 , 13,601.48 km 2 and 15,697.42km 2 , respectively.Over the same time period, the yearly maximal surface water extents are 18,174.76km 2 , 20,671.83km 2 , 19,097.73 km 2 and 18,235.95km 2 , respectively.The total areas of the small surface water and large surface water, whose areas larger than 50 km 2 in the MYRB, are tabulated in Table 4 for each of the studied periods.There was a surge in yearly maximal surface water extent between 1990 and 2000, which may be due to the catastrophic floods of the Yangtze River in 1998.

Area and Distrubution of Surface Water in the MYRB from 1990 to 2017
The study region encompasses over 680,000 km 2 divided into two categories: surface water and non-water.The minimal and maximal surface water extents for each of the studied years are given in Figures 7 and 8.The yearly minimal surface water extents in 1990, 2000, 2010, and 2017 are 14,751.23km 2 , 14,403.48 km 2 , 13,601.48 km 2 and 15,697.42km 2 , respectively.Over the same time period, the yearly maximal surface water extents are 18,174.76km 2 , 20,671.83km 2 , 19,097.73 km 2 and 18,235.95km 2 , respectively.The total areas of the small surface water and large surface water, whose areas larger than 50 km 2 in the MYRB, are tabulated in Table 4 for each of the studied periods.There was a surge in yearly maximal surface water extent between 1990 and 2000, which may be due to the catastrophic floods of the Yangtze River in 1998.In order to quantify specific areal extents of surface water in the four sub-basins of the MYRB.Areas of surface water in these sub-basins were calculated in the study periods, shown as Table 5 and Figure 9.In order to quantify specific areal extents of surface water in the four sub-basins of the MYRB.Areas of surface water in these sub-basins were calculated in the study periods, shown as Table 5 and Figure 9.

Comparison with JRC Yearly Water Classification History
In this section, our surface water results were compared with the JRC Yearly Water Classification History (v1.0) [6].The JRC Yearly Water Classification History (v1.0) dataset contains maps of the location and temporal distribution of surface water from 1984 to 2015.The dataset was generated using 3,066,102 scenes from Landsat 5, 7, and 8 acquired between 16 March 1984 and 10 October 2015.Each pixel was individually classified into water/non-water using an expert system.By means of classification of the seasonality of water throughout the year, surface water was divided into seasonal surface water and permanent surface water.Our yearly minimal and maximal surface water was compared with the permanent surface water and total surface water (seasonal surface water and permanent surface water) of JRC results, respectively.Two typical regions from two sub-basins (MMS and DtLS) were selected to demonstrate the similarities and differences between our results and JRC Global Surface Water, shown as Figure 10.In these figures, it can be seen that the contours and shapes of the surface water obtained by the two methods are comparable.JRC Global Surface Water defined seasonal and permanent surface water by the appearance of surface water in each month through processing every image, and the results are comprehensive.However, the method proposed in this paper is simpler and faster, and the seasonal surface water can be obtained by subtracting the minimum surface water from the maximum surface water.

Comparison with JRC Yearly Water Classification History
In this section, our surface water results were compared with the JRC Yearly Water Classification History (v1.0) [6].The JRC Yearly Water Classification History (v1.0) dataset contains maps of the location and temporal distribution of surface water from 1984 to 2015.The dataset was generated using 3,066,102 scenes from Landsat 5, 7, and 8 acquired between 16 March 1984 and 10 October 2015.Each pixel was individually classified into water/non-water using an expert system.By means of classification of the seasonality of water throughout the year, surface water was divided into seasonal surface water and permanent surface water.Our yearly minimal and maximal surface water was compared with the permanent surface water and total surface water (seasonal surface water and permanent surface water) of JRC results, respectively.Two typical regions from two sub-basins (MMS and DtLS) were selected to demonstrate the similarities and differences between our results and JRC Global Surface Water, shown as Figure 10.In these figures, it can be seen that the contours and shapes of the surface water obtained by the two methods are comparable.JRC Global Surface Water defined seasonal and permanent surface water by the appearance of surface water in each month through processing every image, and the results are comprehensive.However, the method proposed in this paper is simpler and faster, and the seasonal surface water can be obtained by subtracting the minimum surface water from the maximum surface water.

Surface Water Changes Influenced by Natural and Human Factors
The yearly surface water extent in the MYRB fluctuated over the past 27 years, influenced by both natural and human factors.The maximal annual precipitation occurred in 1998 and 2010 (shown in Figure 2), which is consistent with the maximal surface water in 2000 and 2010.Particularly, the

Surface Water Changes Influenced by Natural and Human Factors
The yearly surface water extent in the MYRB fluctuated over the past 27 years, influenced by both natural and human factors.The maximal annual precipitation occurred in 1998 and 2010 (shown in Figure 2), which is consistent with the maximal surface water in 2000 and 2010.Particularly, the yearly maximal surface water extent of the MYRB had dramatically changed between 1990 and 2000 because of the catastrophic floods in 1998.The dynamics of yearly maximal surface water extent in the four sub-basins were similar to the overall situation of the MYRB.Nevertheless, as the annual permanent water, minimal surface water mainly reflects the surface water conditions of lakes, reservoirs, and rivers, and was less affected by the annual precipitation.On the aspect of human activities, both the population and Gross Domestic Product (GDP) in the MYRB increased rapidly in recent decades.Especially the GDP has grown sharply since the late 1990s [71], accompanying with urbanization and industrial production, resulting in that part of the seasonal surface water has become the built-up areas [72].However, the large water bodies of yearly minimal and maximal surface water extent clearly changed by 38.14% and 42.92%, respectively, during 2010-2017, mainly due to China's increasing emphasis on environmental protection and water conservation, especially the protection of key lakes and reservoirs.Meanwhile, relatively small water bodies disappeared over the course of urbanization, keeping the balance of the overall surface water extent.The yearly minimal surface water extent of HjRS (2460.83km 2   9A), these diversities were partially owing to the activities for water resource conservation directed by the local governments.

Advantages and Limitations of Using GEE
For this work, the GEE provided a robust platform and processing environment, providing free accessibility to a series of Landsat image data and to the rapid assessment of surface water dynamics in the MYRB.The GEE platform is designed to allow petabyte-scale data processing, analysis, and visualization [38].The GEE contains more than 200 common datasets, more than 5 million images, including the Landsat scenes, MODIS collections, Sentinel-2 images, and many other remote-sensing images and vector-based datasets.It also provides an integrated environment including cloud-computing and classical algorithms.Compared with traditional image processing tools, such as ENVI, the GEE can process massive images quickly and in batches.In addition, its processing power is constructed based on Google Cloud, not subject to time and place.Scientists can share data and codes via URLs based on the GEE's convenient mechanism.The platform serves as a uniform access interface to a variety of data types, free from image bands, projection, resolution, or bit-depth, making multi-sensor analysis feasible [38].
However, there are a few limitations to using the GEE to monitor land cover.First, if further analysis is required, it is time-consuming to download large-scale classification results, since the pixel-based map contains a high amount of small patches that increase the data volume.Second, there is a lack of entire control of the GEE platform, and the raw images are on line, not in-house.Programs conducted by government agencies may be apprehensive about the accessibility in further research.In addition, the algorithm provided by GEE is not fully explained in some cases.

Conclusions
In this paper, a method of extracting yearly minimal and maximal surface water based on the GEE was proposed, which utilized the algorithm of estimating cloud cover for pixels, images mosaicing based on feature index and random forest classification.And a comparison with existing datasets was present.Using the GEE and Landsat images, temporal changes in the extent of surface water in the Middle Yangtze River Basin were identified through obtaining the minimal and maximal surface water extent in the year of 1990, 2000, 2010 and 2017.Specifically, the areal extent of the surface water in the MYRB showed that yearly minimal surface water had changed to 14,403.48km 2 in 2000 from 14,751.23 km 2 in 1990, and then changed to 13,601.48km 2 in 2010 and changed to 15,697.42 km 2 in 2017.Yearly maximal surface water areas had changed from 18,174.76 km 2 to 20,671.83km 2 first, then consecutively changed to 19,097.73 km 2 and 18,235.95km 2 .The surface water dynamics in the four sub-basins of the MYRB and natural and human factors resulting in these changes were further discussed.The result of this study enhances an understanding of the implications of surface water variations in the developing countries and offers a critical reference for future monitoring and restoration efforts to utilize water resources in the MYRB.
In the proposed method, the pixels with values exceeding certain thresholds of cloud Score, NDVI and NDWI were removed, eliminating a large number of obvious non-water pixels quickly, thus improved the efficiency of the subsequent classification, but also brought some uncertainty, for example water pixels under the cloud may be directly removed.
Although the GEE has gradually become a platform for large-scale remote sensing research, its potential has not been fully explored.This study presents a quick and feasible method for detecting yearly minimal and maximal surface water extents using the GEE.The methodology, combined with the resources available through the GEE, can be easily extended to other domains with similar problems.

Figure 1 .Figure 1 .
Figure 1.An overview of the Middle Yangtze River Basin (MYRB): (a) the location of the MYRB in China; (b) the topography and division of the study area.The average annual precipitation and average temperature of the study area during 1980 to 2016 are shown in Figure 2. Raw data were downloaded from the National Meteorological Information Center (http://data.cma.cn).

Figure 2 .
Figure 2. Average annual precipitation and average temperature changes in the MYRB.

Figure 2 .
Figure 2. Average annual precipitation and average temperature changes in the MYRB.

19 Figure 3 .
Figure 3. Spatial distribution of total observation counts and good observation counts over the MYRB in study periods: (A) total observation in 2017; (B) good observation in 2017; (C) total observation in 2010; (D) good observation in 2010.

Figure 3 .
Figure 3. Spatial distribution of total observation counts and good observation counts over the MYRB in study periods: (A) total observation in 2017; (B) good observation in 2017; (C) total observation in 2010; (D) good observation in 2010.

Figure 4 .
Figure 4. Flowchart of surface water calculation performed in this study.Figure 4. Flowchart of surface water calculation performed in this study.

Figure 4 . 19 Figure 5 .
Figure 4. Flowchart of surface water calculation performed in this study.Figure 4. Flowchart of surface water calculation performed in this study.

Figure 5 .
Figure 5.The greenest and wettest images in the Dongting Lake Section (DtLS) sub-basin in 2017.

Figure 6 .
Figure 6.Impacts of different numbers of trees in the Random Forest (RF) based on the overall accuracy and Kappa accuracy.

Figure 6 .
Figure 6.Impacts of different numbers of trees in the Random Forest (RF) based on the overall accuracy and Kappa accuracy.

Figure 7 .
Figure 7. Yearly minimal and maximal surface water extent.Figure 7. Yearly minimal and maximal surface water extent.

Figure 7 .
Figure 7. Yearly minimal and maximal surface water extent.Figure 7. Yearly minimal and maximal surface water extent.

Figure 9 .BFigure 9 .
Figure 9. Dynamics of surface water extents in four sub-basins based on minimal and maximal extent.(A) Areal changes in minimal surface water extent in all four sub-basins.(B) Area changes in maximal surface water extent in all four sub-basins.

Figure 10 .
Figure 10.Comparisons between our results and Joint Research Centre (JRC) Yearly Water Classification History (Minimum represents minimal surface water, Maximum represents maximal surface water, JRC_PW represents permanent surface water of the JRC Yearly Water Classification History, and JRC_All represents the sum of the annual seasonal surface water and permanent surface water.).(A) Two results in the typical area of the Middle Mainstream Section (MMS) in 2000.(B) Two results in the typical area of the MMS in 2010.(C) Two results in the typical area of the Dongting Lake Section (DtLS) in 2000.(D) Two results in the typical area of the DtLS in 2010.

Figure 10 .
Figure 10.Comparisons between our results and Joint Research Centre (JRC) Yearly Water Classification History (Minimum represents minimal surface water, Maximum represents maximal surface water, JRC_PW represents permanent surface water of the JRC Yearly Water Classification History, and JRC_All represents the sum of the annual seasonal surface water and permanent surface water.).(A) Two results in the typical area of the Middle Mainstream Section (MMS) in 2000.(B) Two results in the typical area of the MMS in 2010.(C) Two results in the typical area of the Dongting Lake Section (DtLS) in 2000.(D) Two results in the typical area of the DtLS in 2010.

Table 1 .
Properties of image collections selected for this study.

Table 1 .
Properties of image collections selected for this study.

Table 2 .
Detail information of reference data.

Table 2 .
Detail information of reference data.

Table 3 .
Results of error matrixes of land cover in 1990, 2000, 2010, and 2017.The value after the symbol "±" represents the margin of error with confidence level 95%.
Notes: the best classification results were marked with bold font, and the worst classification results were marked with italic font.

Table 3 .
Results of error matrixes of land cover in 1990, 2000, 2010, and 2017.The value after the symbol "±" represents the margin of error with confidence level 95%.
Notes: the best classification results were marked with bold font, and the worst classification results were marked with italic font.

Table 4 .
Areas and change rates of large surface water and others of minimal and maximal surface water extent in the study periods.

Table 4 .
Areas and change rates of large surface water and others of minimal and maximal surface water extent in the study periods.

Table 5 .
Minimal and Maximum Surface water extent in the four sub-basins of the MYRB in the study periods.
in 1990, 2361.32 km 2 in 2000, 2293.17km 2 in 2010 and 2170.22 km 2 in 2017) and PyLS (4653.30km 2 in 1990, 4618.39 km 2 in 2000, 4357.58 km 2 in 2010 and 3906.24km 2 in 2017) showed continuous changes, but the yearly minimal surface water extent of MMS and DtLS changed from 3408.31 km 2 to 4768.98 km 2 and from 3551.97 km 2 to 4204.66 km 2 respectively in 2010-2017 (shown in Table 5 and Figure