Dynamic Monitoring of Surface Water Area during 1989–2019 in the Hetao Plain Using Landsat Data in Google Earth Engine

: The spatio-temporal change of the surface water is very important to agricultural, economic, and social development in the Hetao Plain, as well as the structure and function of the ecosystem. To understand the long-term changes of the surface water area in the Hetao Plain, we used all available Landsat images (7534 scenes) and adopted the modiﬁed Normalized Di ﬀ erence Water Index (mNDWI), Enhanced Vegetation Index (EVI), and Normalized Di ﬀ erence Vegetation Index (NDVI) to map the open-surface water from 1989 to 2019 in the Google Earth Engine (GEE) cloud platform. We further analyzed precipitation, temperature, and irrigated area, revealing the impact of climate change and human activities on long-term surface water changes. The results show the following. (1) In the last 31 years, the maximum, seasonal, and annual average water body area values in the Hetao Plain have exhibited a downward trend. Meanwhile, the number of maximum, seasonal, and permanent water bodies displayed a signiﬁcant upward trend. (2) The variation of the surface water area in the Hetao Plain is mainly a ﬀ ected by the maximum water body area, while the variation of the water body number is mainly a ﬀ ected by the number of minimum water bodies. (3) Precipitation has statistically signiﬁcant positive e ﬀ ects on the water body area and water body number, which has statistically signiﬁcant negative e ﬀ ects with temperature and irrigation. The ﬁndings of this study can be used to help the policy-makers and farmers understand changing water resources and its driving mechanism and provide a reference for water resources management, agricultural irrigation, and protection.


Introduction
As an important part of the land-water cycle, surface water resources play an important role in promoting national economic development and maintaining the balance of terrestrial and aquatic ecosystems, agriculture, and the ecological environment [1]. In terms of global climate change, all kinds of water resources show obvious characteristics of temporal and spatial differentiation [2][3][4][5]. It is necessary to strengthen the dynamic monitoring and investigation of water resources, especially for arid and semi-arid areas. The Yellow River Basin is an important ecological barrier and economic At the same time, the ability of cloud computing has been significantly improved in the recent years, and it has shown great application potential in large-scale land cover mapping. For example, the Google Earth Engine (GEE), with high-performance parallel computing capabilities, massive remote sensing, and geospatial data, is free to use [25,49]. It facilitates the use of all images archived on the platform for global and national land cover remote sensing mapping and land cover change monitoring [10,[34][35][36][50][51][52][53][54]. Based on the GEE platform, 7534 remote sensing images of Landsat Thematic Mapper (TM), Landsat Enhanced Thematic Mapper (ETM+), and Landsat Operational Land Imager (OLI) captured from 1989 to 2019 were obtained by JavaScript programming operation. The method of combining the modified Normalized Difference Water Index (mNDWI), Enhanced Vegetation Index (EVI), and Normalized Difference Vegetation Index (NDVI) was used to extract the accuracy of surface water extraction. The temporal and spatial distribution pattern of surface water in the Hetao Plain in the recent 31 years was obtained, and four water types were generated based on the water body frequency: maximum water body, seasonal water body, permanent water body, and annual average water body [23,55,56]. The interannual variation trend of the surface water body area and water body number of four water types was analyzed. Finally, the relationship between the changing trend of surface water and the major driven factors was discussed.

Study Area
The Hetao Plain is located in the middle reaches of the Yellow River Basin (40 • to 41 • 18 N, 106 • 21 to 111 • 50 E), north of the Great Wall and south of the Yinshan Mountains, from Helan Mountain in the west and east to Hohhot in the east (Figure 1). The Hetao Plain is formed by the alluvium of the Yellow River, and some of its tributaries, extending from east to west along the Yellow River, with a length of about 591.6 km, a width of 20-90 km from north to south, and an area of 23,575 km 2 , are fanned out in an arc. The administrative departments belong to Alxa League, Bayannaoer, Erdos, Baotou, and Hohhot in the Inner Mongolia Autonomous Region, and the geomorphology of the Hetao Plain is mainly alluvial Plain, the average slope is 1.8º ( Figure S1). The elevation of the study area is high in the west and low in the east. The Hetao Plain is located in an arid and semi-arid region. It is hot and humid in the summer, rainy and hot in the same period, and cold and dry in the winter, and there is a large temperature difference between the day and night. From 1989 to 2019, the annual average temperature is 8.9 • C, and the temperature change decreases gradually from west to east. The annual cumulative precipitation is 272 mm, and the trend increases from west to east. The lowest annual cumulative precipitation was 173 mm in 1991, and the highest annual cumulative precipitation value of 454 mm occurred in 2012. Due to a large amount of evaporation and lack of precipitation, the agricultural water consumption in the Hetao Plain mainly comes from the Yellow River, which is a typical irrigated agricultural region. The main crops are spring wheat, corn, and sunflower.
Water 2020, 12, x FOR PEER REVIEW 3 of 22 Google Earth Engine (GEE), with high-performance parallel computing capabilities, massive remote sensing, and geospatial data, is free to use [25,49]. It facilitates the use of all images archived on the platform for global and national land cover remote sensing mapping and land cover change monitoring [10,[34][35][36][50][51][52][53][54]. Based on the GEE platform, 7534 remote sensing images of Landsat Thematic Mapper (TM), Landsat Enhanced Thematic Mapper (ETM+), and Landsat Operational Land Imager (OLI) captured from 1989 to 2019 were obtained by JavaScript programming operation. The method of combining the modified Normalized Difference Water Index (mNDWI), Enhanced Vegetation Index (EVI), and Normalized Difference Vegetation Index (NDVI) was used to extract the accuracy of surface water extraction. The temporal and spatial distribution pattern of surface water in the Hetao Plain in the recent 31 years was obtained, and four water types were generated based on the water body frequency: maximum water body, seasonal water body, permanent water body, and annual average water body [23,55,56]. The interannual variation trend of the surface water body area and water body number of four water types was analyzed. Finally, the relationship between the changing trend of surface water and the major driven factors was discussed.

Study Area
The Hetao Plain is located in the middle reaches of the Yellow River Basin (40° to 41°18′ N, 106°21′ to 111°50′ E), north of the Great Wall and south of the Yinshan Mountains, from Helan Mountain in the west and east to Hohhot in the east (Figure 1). The Hetao Plain is formed by the alluvium of the Yellow River, and some of its tributaries, extending from east to west along the Yellow River, with a length of about 591.6 km, a width of 20-90 km from north to south, and an area of 23,575 km 2 , are fanned out in an arc. The administrative departments belong to Alxa League, Bayannaoer, Erdos, Baotou, and Hohhot in the Inner Mongolia Autonomous Region, and the geomorphology of the Hetao Plain is mainly alluvial Plain, the average slope is 1.8º ( Figure S1). The elevation of the study area is high in the west and low in the east. The Hetao Plain is located in an arid and semi-arid region. It is hot and humid in the summer, rainy and hot in the same period, and cold and dry in the winter, and there is a large temperature difference between the day and night. From 1989 to 2019, the annual average temperature is 8.9 °C, and the temperature change decreases gradually from west to east. The annual cumulative precipitation is 272 mm, and the trend increases from west to east. The lowest annual cumulative precipitation was 173 mm in 1991, and the highest annual cumulative precipitation value of 454 mm occurred in 2012. Due to a large amount of evaporation and lack of precipitation, the agricultural water consumption in the Hetao Plain mainly comes from the Yellow River, which is a typical irrigated agricultural region. The main crops are spring wheat, corn, and sunflower.

Landsat Data
Based on the surface reflection data set of Landsat TM, Landsat ETM+, and Landsat OLI on the GEE platform (https://earthengine.google.org/), a total of 7534 images of the study area were obtained (1 January 1989 to 31 December 2019) (Figure 2a). According to the Earth Resources Satellite Global Reference System (WRS-2), a total of nine images can cover the Hetao Plain. The image data covered by each scene are relatively balanced (Figure 2c), and the distribution of images in each month is also relatively balanced (Figure 2d). Landsat TM and Landsat ETM+ surface reflectance data are corrected by the Landsat Ecosystem Interference Adaptive Processing System (LEDAPS) algorithm. The Landsat OLI surface reflectance data use the Landsat Surface Reflectance Coding (LASRC) algorithm for atmospheric correction [57,58]. At the same time, these data sets also use the cloud, shadow, water, and snow masks generated by the code based on the Function of Mask algorithm (CFMASK) in the Surface Reflectance (SR) collection on GEE as bad observations and masked before the water detection [59,60]   Based on the surface reflection data set of Landsat TM, Landsat ETM+, and Landsat OLI on the GEE platform (https://earthengine.google.org/), a total of 7534 images of the study area were obtained (January 1, 1989 to December 31, 2019) ( Figure 2a). According to the Earth Resources Satellite Global Reference System (WRS-2), a total of nine images can cover the Hetao Plain. The image data covered by each scene are relatively balanced (Figure 2c), and the distribution of images in each month is also relatively balanced (Figure 2d). Landsat TM and Landsat ETM+ surface reflectance data are corrected by the Landsat Ecosystem Interference Adaptive Processing System (LEDAPS) algorithm. The Landsat OLI surface reflectance data use the Landsat Surface Reflectance Coding (LASRC) algorithm for atmospheric correction [57,58]. At the same time, these data sets also use the cloud, shadow, water, and snow masks generated by the code based on the Function of Mask algorithm (CFMASK) in the Surface Reflectance (SR) collection on GEE as bad observations and masked before the water detection [59,60]

Climate Data
The annual mean temperature and annual cumulative precipitation were calculated based on the 3-hour Global Land Data Assimilation System (GLDAS), for GLDAS fusion satellite and groundbased observation data products; using advanced surface modeling and data assimilation technology, the best surface state and flux field can be generated [61]. GLDAS-2 includes the two datasets of GLDAS-2.0 (1989-1999 years) and GLDAS-2.1 (2000-2019 years).

Sentinel-2 Multispectral Instruments (MSI) Images Data
Sentinel-2 MSI is a wide-scan, high-resolution, multi-spectral earth observation sensor for carrying out Copernicus land monitoring research, it has been successfully applied for mapping and monitoring vegetation, soil, and water cover [26,27]. In the study, Sentinel-2 MSI images were collected from July 1, 2019 to July 31, 2019. The accuracy of the maximum surface water body area extracted based on Landsat images in the same period was verified. In total, 140 Sentinel-2 MSI images were required to cover this study area.

Globeland 30 Data
The accuracy of the maximum surface water area extracted from Landsat images was verified, ensuring a uniform distribution of verification points in water and non-water bodies. We used the 2010 Global 30-meter land cover remote sensing data product (GlobeLand30) (http://www.globallandcover.com) for the classification of data. The surface cover was divided into water bodies and non-water bodies, 1500 random points of water bodies and non-water bodies were

Climate Data
The annual mean temperature and annual cumulative precipitation were calculated based on the 3-hour Global Land Data Assimilation System (GLDAS), for GLDAS fusion satellite and ground-based observation data products; using advanced surface modeling and data assimilation technology, the best surface state and flux field can be generated [61]. GLDAS-2 includes the two datasets of GLDAS-2.0 (1989-1999 years) and GLDAS-2.1 (2000-2019 years).

Sentinel-2 Multispectral Instruments (MSI) Images Data
Sentinel-2 MSI is a wide-scan, high-resolution, multi-spectral earth observation sensor for carrying out Copernicus land monitoring research, it has been successfully applied for mapping and monitoring vegetation, soil, and water cover [26,27]. In the study, Sentinel-2 MSI images were collected from 1 July 2019 to 31 July 2019. The accuracy of the maximum surface water body area extracted based on Landsat images in the same period was verified. In total, 140 Sentinel-2 MSI images were required to cover this study area.

Globeland 30 Data
The accuracy of the maximum surface water area extracted from Landsat images was verified, ensuring a uniform distribution of verification points in water and non-water bodies. We used the 2010 Global 30-meter land cover remote sensing data product (GlobeLand30) (http://www.globallandcover. com) for the classification of data. The surface cover was divided into water bodies and non-water bodies, 1500 random points of water bodies and non-water bodies were generated in the Hetao Plain, and a total of 3000 verification points were uniformly distributed in the study area [62,63] (Figure 3).
Water 2020, 12, x FOR PEER REVIEW 6 of 22 generated in the Hetao Plain, and a total of 3000 verification points were uniformly distributed in the study area [62,63] (Figure 3).

Surface Water Extraction Algorithm
The commonly used water extraction algorithms can be divided into two categories：machine learning algorithms (SVM [38], RF [39], DT [40], DL [41], etc.) and traditional algorithms (single-band method and multi-band method (the inter-spectral relation method and water index method)) [43,44]. Previous studies have shown that there is a significant difference in the response of the different reflectivity of remote sensing images to water bodies and non-water bodies, and water bodies can be effectively extracted by using the index of different spectral combinations and their thresholds [64][65][66][67]. In this paper, the water body index is employed (modified normalized difference water index (mNDWI, Equation (1)), normalized differential vegetation index (NDVI, Equation (2)), and enhanced vegetation index (EVI, Equation (3)). The existing water algorithm is selected and verified, when the spectral index of the pixel satisfies mNDWI-NDVI > 0 and EVI < 0.1 or mNDWI-EVI > 0 and EVI < 0.1, the pixels are water body; otherwise, they are non-water body (Equation (4)

Surface Water Extraction Algorithm
The commonly used water extraction algorithms can be divided into two categories: machine learning algorithms (SVM [38], RF [39], DT [40], DL [41], etc.) and traditional algorithms (single-band method and multi-band method (the inter-spectral relation method and water index method)) [43,44]. Previous studies have shown that there is a significant difference in the response of the different reflectivity of remote sensing images to water bodies and non-water bodies, and water bodies can be effectively extracted by using the index of different spectral combinations and their thresholds [64][65][66][67]. In this paper, the water body index is employed (modified normalized difference water index (mNDWI, Equation (1)), normalized differential vegetation index (NDVI, Equation (2)), and enhanced vegetation index (EVI, Equation (3)). The existing water algorithm is selected and verified, when the spectral index of the pixel satisfies mNDWI-NDVI > 0 and EVI < 0.1 or mNDWI-EVI > 0 and EVI < 0.1, the pixels are water body; otherwise, they are non-water body (Equation (4)) [23,55,56]. The verification results of the algorithm in the study area are shown in Figure 4.

Calculation Method of Maximum, Seasonal, Permanent, and Annual Average Water Body
According to the water frequency (WF) (Equation (5)), the surface water extent was divided into four indicators: (1) the maximum extent of the water body in a year is called the maximum water body (WF ≥ 0.25); (2) the seasonal extent of the water body, that is, the difference between the maximum and permanent water body, is called the seasonal water body (0.25 ≤ WF < 0.75); (3) a persistent yearly water body extent is called a permanent water body (WF ≥ 0.75); and (4) the annual average extent of the water body is called the annual average water body [23,55,56].
Here, WF is the frequency of the water body, Water N represents the number of times that all images were identified as water bodies during the year, and Good N represents the number of highquality images during the year. (e) pixels that satisfy the conditions of ((mNDWI > NDVI) vs. (EVI < 0.1)) are water bodies, marked with red scattered points; (f) the red part of the image corresponds to the red scatter in the (e) diagram, indicating the water body; (g) pixels scatter diagram between (mNDWI > EVI) and EVI; (h) pixels that satisfy the conditions of ((mNDWI > EVI) vs. (EVI < 0.1)) are water bodies, marked with red scattered points; (i) the red part of the image corresponds to the red scatter in the (h) diagram, indicating the water body.

Calculation Method of Maximum, Seasonal, Permanent, and Annual Average Water Body
According to the water frequency (WF) (Equation (5)), the surface water extent was divided into four indicators: (1) the maximum extent of the water body in a year is called the maximum water body (WF ≥ 0.25); (2) the seasonal extent of the water body, that is, the difference between the maximum and permanent water body, is called the seasonal water body (0.25 ≤ WF < 0.75); (3) a persistent yearly water body extent is called a permanent water body (WF ≥ 0.75); and (4) the annual average extent of the water body is called the annual average water body [23,55,56].
Here, WF is the frequency of the water body, N Water represents the number of times that all images were identified as water bodies during the year, and N Good represents the number of high-quality images during the year.

Surface Water Area
The maximum, seasonal, and permanent water body area refers to the number of such pixels multiplied by the area of a single pixel (Equation (6)) [49]. The annual average water area is the product of the effective water pixels and the water frequency in one year.
Here, S is the area of the maximum, seasonal, and permanent water body of different types (km 2 ), N is the total number of pixels, A is the number of pixels, and pixelsize is the pixel area (m 2 ).

Dynamic Degree of Surface Water
The dynamic degree of surface water represents the change speed of the surface water body area and water body number in unit time (Equation (7) Here, K is the dynamic degree of surface water; S a and S b represent the area and number of the surface water body at the beginning and end of the study, respectively; and T is the research period, in years. The larger the K value, the faster the change rate of surface water in unit time.

Accuracy Evaluation Method
In this study, Sentinel-2 MSI images with a 10 m spatial resolution were used to verify the surface water extraction results with a 30-m spatial resolution (Landsat). Firstly, to maintain a balanced distribution of water and non-water sample points, 1500 random points of water and non-water were generated based on GlobeLand30 water and non-water types. Secondly, 3000 verification points and the water extracted in this study were superimposed on Sentinel-2 MSI images. Visual discrimination of the accuracy of water extraction results and the accuracy evaluation was based on the producer accuracy, user accuracy, overall accuracy, and Kappa coefficient of the confusion matrix [35] (Figure 3).

Accuracy of Water Body Mapping
Since the maximum, seasonal, permanent, and average water body area were obtained from the same method and images, the classification accuracies of these maps are comparable. Therefore, this study selects the maximum water body area as the verification of the accuracy of our results. As shown in Table 1, the confusion matrix is based on the accurate evaluation of the Sentinel-2 MSI images and the largest water system map in the Hetao Plain in 2019. The overall accuracy was 95.9%, the producer accuracy was 93.24%, and the Kappa coefficient was 0.91, which indicates that the water extraction results have high accuracy.

Spatial Distribution of Surface Water
There are 1,430,000 and 830,000 water pixels in the frequency map of the water body in 2019 and the frequency map of the cumulative water body from 1989 to 2019, respectively. The pixels area with WF ≥ 0.25 accounted for 3.17% and 5.46% of the total area of the study area, respectively (Figure 5a,b). There is a great difference in the spatial distribution of water pixels with different frequencies between 2019 and 1989-2019, where the number of water pixels with WF ≥ 0.75 accounted for 51.97% and 32.99% of the total number of pixels, respectively (Figure 5c,d). These water pixels represent the main rivers, large lakes, and reservoirs that make up the study area and can maintain water throughout the year, which are the main sources of surface water. The number of water pixels with 0.25 ≤ WF < 0.75 accounts for 48.03% and 67.01% of the total, respectively. The water pixels are the edge of the stream, pond, and large water body, which forms the wet season and the dry season with the seasonal change. For example, the frequency of the water body in the central area of the Wuliangsuhai Lake is close to 1, which means that the water level is relatively deep and the water body exists during the high and low water periods. There are more water bodies in the edge part and the frequency is lower than 0.5, which indicates that the water level is shallow and the water body exists in the wet season and may dry up in the dry season. Figure

Spatial Distribution of Surface Water
There are 1,430,000 and 830,000 water pixels in the frequency map of the water body in 2019 and the frequency map of the cumulative water body from 1989 to 2019, respectively. The pixels area with WF ≥ 0.25 accounted for 3.17% and 5.46% of the total area of the study area, respectively (Figure 5a,  b). There is a great difference in the spatial distribution of water pixels with different frequencies between 2019 and 1989-2019, where the number of water pixels with WF ≥ 0.75 accounted for 51.97% and 32.99% of the total number of pixels, respectively (Figure 5c,d). These water pixels represent the main rivers, large lakes, and reservoirs that make up the study area and can maintain water throughout the year, which are the main sources of surface water. The number of water pixels with 0.25 ≤ WF < 0.75 accounts for 48.03% and 67.01% of the total, respectively. The water pixels are the edge of the stream, pond, and large water body, which forms the wet season and the dry season with the seasonal change. For example, the frequency of the water body in the central area of the Wuliangsuhai Lake is close to 1, which means that the water level is relatively deep and the water body exists during the high and low water periods. There are more water bodies in the edge part and the frequency is lower than 0.5, which indicates that the water level is shallow and the water body exists in the wet season and may dry up in the dry season. Figure

Temporal Distribution of Surface Water
From 1989 to 2019, the maximum, seasonal, and annual average water body area exhibited a downward trend, while the permanent water body area displayed an upward trend, as shown in Figure 6. The maximum water body area is 816-1827 km 2 , which is −26% to 65% of the average value (1107 km 2 ). The permanent water body area is 346-757 km 2 , which is −33% to 46% of the average value (519 km 2 ), after 2001, the permanent water body area showed an upward trend and was

Temporal Distribution of Surface Water
From 1989 to 2019, the maximum, seasonal, and annual average water body area exhibited a downward trend, while the permanent water body area displayed an upward trend, as shown in Figure 6. The maximum water body area is 816-1827 km 2 , which is −26% to 65% of the average value (1107 km 2 ). The permanent water body area is 346-757 km 2 , which is −33% to 46% of the average value (519 km 2 ), after 2001, the permanent water body area showed an upward trend and was positively correlated with annual cumulative precipitation ( Figure S2, Table S1). The seasonal water body area is 298-1161 km 2 , which is −49% to 97% of the average value (588 km 2 ). The annual average water body area is 598-1225 km 2 , which is −28% to 47% of the average value (831 km 2 ). The surface water body area in the Hetao Plain displays an overall downward trend, with the average annual water body area decreasing by 1.9 km 2 every year.
Water 2020, 12, x FOR PEER REVIEW 12 of 22 positively correlated with annual cumulative precipitation ( Figure S2, Table S1). The seasonal water body area is 298-1161 km 2 , which is −49% to 97% of the average value (588 km 2 ). The annual average water body area is 598-1225km 2 , which is −28% to 47% of the average value (831 km 2 ). The surface water body area in the Hetao Plain displays an overall downward trend, with the average annual water body area decreasing by 1.9 km 2 every year. The amount of surface water was divided into the maximum, seasonal, and permanent water bodies, as shown in Figure 7. From 1989 to 2019, the number of maximum water bodies was 232,000-635,000, which is −40% to 63% of the average value (389,000). The number of seasonal water bodies was 173,000-495,000, which is −39% to 74% of the average value (284,000). The number of permanent water bodies was 59,000-162,000, which is −45% to 50% of the average value (108,000). The number of maximum water bodies (R 2 = 0.3763), seasonal water bodies (R 2 = 0.2188), and permanent water bodies (R 2 = 0.7058) showed an overall upward trend over the 31 years. The number (a) and area (b) distribution of the maximum surface water at different spatial scales from 1989 to 2019 are shown in Figure 8. According to the size of the maximum water body area, the water body area was divided into 10 categories. The number of water bodies larger than 0.5 km 2 was 85, accounting for 99.95% of the total water body area and 0.022% of the total number of water bodies. The number of water bodies less than 0.005 km 2 was 377,803, accounting for 0.00008% of the total water body area and 96.37% of the total water body number. Therefore, the change of water body area in the Hetao Plain is mainly affected by large water bodies, while the change in the number of water bodies is mainly affected by small water bodies. The amount of surface water was divided into the maximum, seasonal, and permanent water bodies, as shown in Figure 7. From 1989 to 2019, the number of maximum water bodies was 232,000-635,000, which is −40% to 63% of the average value (389,000). The number of seasonal water bodies was 173,000-495,000, which is −39% to 74% of the average value (284,000). The number of permanent water bodies was 59,000-162,000, which is −45% to 50% of the average value (108,000). The number of maximum water bodies (R 2 = 0.3763), seasonal water bodies (R 2 = 0.2188), and permanent water bodies (R 2 = 0.7058) showed an overall upward trend over the 31 years. positively correlated with annual cumulative precipitation ( Figure S2, Table S1). The seasonal water body area is 298-1161 km 2 , which is −49% to 97% of the average value (588 km 2 ). The annual average water body area is 598-1225km 2 , which is −28% to 47% of the average value (831 km 2 ). The surface water body area in the Hetao Plain displays an overall downward trend, with the average annual water body area decreasing by 1.9 km 2 every year. The amount of surface water was divided into the maximum, seasonal, and permanent water bodies, as shown in Figure 7. From 1989 to 2019, the number of maximum water bodies was 232,000-635,000, which is −40% to 63% of the average value (389,000). The number of seasonal water bodies was 173,000-495,000, which is −39% to 74% of the average value (284,000). The number of permanent water bodies was 59,000-162,000, which is −45% to 50% of the average value (108,000). The number of maximum water bodies (R 2 = 0.3763), seasonal water bodies (R 2 = 0.2188), and permanent water bodies (R 2 = 0.7058) showed an overall upward trend over the 31 years. The number (a) and area (b) distribution of the maximum surface water at different spatial scales from 1989 to 2019 are shown in Figure 8. According to the size of the maximum water body area, the water body area was divided into 10 categories. The number of water bodies larger than 0.5 km 2 was 85, accounting for 99.95% of the total water body area and 0.022% of the total number of water bodies. The number of water bodies less than 0.005 km 2 was 377,803, accounting for 0.00008% of the total water body area and 96.37% of the total water body number. Therefore, the change of water body area in the Hetao Plain is mainly affected by large water bodies, while the change in the number of water bodies is mainly affected by small water bodies. The number (a) and area (b) distribution of the maximum surface water at different spatial scales from 1989 to 2019 are shown in Figure 8. According to the size of the maximum water body area, the water body area was divided into 10 categories. The number of water bodies larger than 0.5 km 2 was 85, accounting for 99.95% of the total water body area and 0.022% of the total number of water bodies. The number of water bodies less than 0.005 km 2 was 377,803, accounting for 0.00008% of the total water body area and 96.37% of the total water body number. Therefore, the change of water body area in the Hetao Plain is mainly affected by large water bodies, while the change in the number of water bodies is mainly affected by small water bodies. In this study, the dynamic attitudes of the maximum number and area of surface water in six different periods of 1989-1994, 1995-1999, 2000-2004, 2005-2009   In this study, the dynamic attitudes of the maximum number and area of surface water in six different periods of 1989-1994, 1995-1999, 2000-2004, 2005-2009, 2010-2014, and 2015-2019 were compared (Tables 2 and 3). The maximum number of water bodies appeared in 2005-2009, with a value of 464,241, and the minimum value appeared in 1989-1994, with a value of 300,169. The maximum water body area appeared in 1989-1994, and its value was 1246.471 km 2 , whilst the minimum value appeared in 2000-2004, and its value was 935.649 km 2 . From 1989 to 2019, the maximum number of water bodies exhibited an upward trend, with an overall dynamic change of 45.75%. However, the maximum water body area showed a downward trend, and the overall dynamic change was −12.09%.

Influence of Drought and Rainy Years on Surface Water Change
Precipitation is one of the important factors affecting the change in the surface water body area and water body number. In 2011 and 2012, the annual cumulative precipitation in the Hetao Plain was 246 and 454 mm, respectively. Compared with the 31-year annual cumulative precipitation of 272 mm, 2011 had a lower value than the annual cumulative precipitation, which is called a dry year, and 2012 had a higher value than the annual cumulative precipitation, which is called a rainy year. On the whole, the area and water body number of surface water in 2012 was larger than those in 2011 (Figure 9). The number of water bodies in 2012 was 900,000, which was 340,000 more than in 2011 (560,000). The number of water bodies less than 0.005 km 2 accounted for 85.3% of the total. Precipitation is an important reason for the formation of small water bodies, and the amount of precipitation directly affects the number of water bodies. The water body area in 2012 was 1827 km 2 , which was 1011 km 2 more than that in 2011 (816 km 2 ). Water body areas of more than 50 km 2 account for 86.6% of the total water body area, indicating that the change in the surface water body area is mainly affected by large surface water bodies.
Water 2020, 12, x FOR PEER REVIEW 14 of 22

Influence of Drought and Rainy Years on Surface Water Change
Precipitation is one of the important factors affecting the change in the surface water body area and water body number. In 2011 and 2012, the annual cumulative precipitation in the Hetao Plain was 246 and 454 mm, respectively. Compared with the 31-year annual cumulative precipitation of 272 mm, 2011 had a lower value than the annual cumulative precipitation, which is called a dry year, and 2012 had a higher value than the annual cumulative precipitation, which is called a rainy year. On the whole, the area and water body number of surface water in 2012 was larger than those in 2011 ( Figure 9). The number of water bodies in 2012 was 900,000, which was 340,000 more than in 2011 (560,000). The number of water bodies less than 0.005 km 2 accounted for 85.3% of the total. Precipitation is an important reason for the formation of small water bodies, and the amount of precipitation directly affects the number of water bodies. The water body area in 2012 was 1827 km 2 , which was 1011 km 2 more than that in 2011 (816 km 2 ). Water body areas of more than 50 km 2 account for 86.6% of the total water body area, indicating that the change in the surface water body area is mainly affected by large surface water bodies.

Influence of Climate Change and Human Activities on Surface Water Change
The drivers of open surface water dynamics include both climate change and human activities [63,[69][70][71]; here, we used annual cumulative precipitation (AP) and annual average temperature (AT) as measures of climatic factors, and irrigation as measures of human activities factors. The results of multiple linear regression analysis are shown in Table 4. Precipitation had statistically significant positive effects on the maximum water body area and water body number, seasonal water body area and water body number, and permanent water body area and water body number. The more precipitation, the larger the water body area and the higher the water body number. The temperature had statistically significant negative effects on the maximum water body area and water body number, seasonal water body area and water body number, and permanent water body area and water body number. A higher temperature will reduce the air pressure, reduce the content of water vapor in the air, and increase evaporation, thus increasing the demand for agricultural water, resulting in a reduction of the surface water body area [23]. Figure 10 shows the overall upward trend of precipitation and temperature in the Hetao Plain in the past 31 years.

Influence of Climate Change and Human Activities on Surface Water Change
The drivers of open surface water dynamics include both climate change and human activities [63,[69][70][71]; here, we used annual cumulative precipitation (AP) and annual average temperature (AT) as measures of climatic factors, and irrigation as measures of human activities factors. The results of multiple linear regression analysis are shown in Table 4. Precipitation had statistically significant positive effects on the maximum water body area and water body number, seasonal water body area and water body number, and permanent water body area and water body number. The more precipitation, the larger the water body area and the higher the water body number. The temperature had statistically significant negative effects on the maximum water body area and water body number, seasonal water body area and water body number, and permanent water body area and water body number. A higher temperature will reduce the air pressure, reduce the content of water vapor in the air, and increase evaporation, thus increasing the demand for agricultural water, resulting in a reduction of the surface water body area [23]. Figure 10 shows the overall upward trend of precipitation and temperature in the Hetao Plain in the past 31 years. Table 4. Multiple linear regression analyses of the water body area and water body number with climatic and anthropogenic driving factors in the Hetao Plain. The six dependent variables are maximum water body area, maximum water body number, year-long water body area, year-long water body number, seasonal water body area, and seasonal water body number. AP and AT are the Hetao Plain annual cumulative precipitation and annual average temperature respectively, the three dependent variables are AP, AT, and Irrigation. R 2 is the proportion of variance in the dependent variable, which can be explained by the selected explanatory variables. SEE is the standard error of the estimate. F and Sig. are the F-statistic and the p-value associated with it. For each linear regression, the variables were statistically significant (P < 0.05).

Maximum
Seasonal Permanent  Irrigation directly consumed groundwater and surface water resources during irrigation transportation as part of the evaporation or leakage loss, so irrigation had statistically significant negative effects on the water body area and water body number. The irrigated area in the Hetao Plain increased from 1989 (5161km 2 ) to 2019 (12410.8 km 2 ), and it showed a continuous rise at a rate of 217.91 km 2 /year ( Figure 11). The sown area in the Hetao Plain increased from 1989 (8725.8km 2 ) to 2019 (20388.6 km 2 ), and showed a continuous rise at a rate of 415.09km 2 /year ( Figure S3). Figure S6, it can be seen that the irrigation area accounts for more than 50% of the sown area in Hetao Plain. The population and gross domestic product (GDP) of the Hetao Plain have increased rapidly in the past 31 years ( Figure S4). Regarding the development of industrialization and urbanization, some surface water body areas have become building land ( Figure S5). Due to the backward irrigation technology and insufficient awareness of saving water in the Hetao Plain irrigated area, the surface water has gradually decreased. The change of surface water in arid and semi-arid areas is a comprehensive reflection of the impact of climate change and human activities on regional water resources. Irrigation directly consumed groundwater and surface water resources during irrigation transportation as part of the evaporation or leakage loss, so irrigation had statistically significant negative effects on the water body area and water body number. The irrigated area in the Hetao Plain increased from 1989 (5161 km 2 ) to 2019 (12,410.8 km 2 ), and it showed a continuous rise at a rate of 217.91 km 2 /year ( Figure 11). The sown area in the Hetao Plain increased from 1989 (8725.8 km 2 ) to 2019 (20,388.6 km 2 ), and showed a continuous rise at a rate of 415.09 km 2 /year ( Figure S3). Figure S6, it can be seen that the irrigation area accounts for more than 50% of the sown area in Hetao Plain. The population and gross domestic product (GDP) of the Hetao Plain have increased rapidly in the past 31 years ( Figure S4). Regarding the development of industrialization and urbanization, some surface water body areas have become building land ( Figure S5). Due to the backward irrigation technology and insufficient awareness of saving water in the Hetao Plain irrigated area, the surface water has gradually decreased. The change of surface water in arid and semi-arid areas is a comprehensive reflection of the impact of climate change and human activities on regional water resources.
it can be seen that the irrigation area accounts for more than 50% of the sown area in Hetao Plain. The population and gross domestic product (GDP) of the Hetao Plain have increased rapidly in the past 31 years ( Figure S4). Regarding the development of industrialization and urbanization, some surface water body areas have become building land ( Figure S5). Due to the backward irrigation technology and insufficient awareness of saving water in the Hetao Plain irrigated area, the surface water has gradually decreased. The change of surface water in arid and semi-arid areas is a comprehensive reflection of the impact of climate change and human activities on regional water resources.  Table 4. Multiple linear regression analyses of the water body area and water body number with climatic and anthropogenic driving factors in the Hetao Plain. The six dependent variables are maximum water body area, maximum water body number, year-long water body area, year-long water body number, seasonal water body area, and seasonal water body number. AP and AT are the Hetao Plain annual cumulative precipitation and annual average temperature respectively, the three dependent variables are AP, AT, and Irrigation. R 2 is the proportion of variance in the dependent variable, which can be explained by the selected explanatory variables. SEE is the standard error of the estimate. F and Sig. are the F-statistic and the p-value associated with it. For each linear regression, the variables were statistically significant (P < 0.05).

Comparison with the JRC Dataset
This research uses the Yearly Water Classification History (V1.1) dataset from the Joint Research Center (JRC) of the European Commission [21]. It is composed of 3,865,618 Landsat TM, Landsat ETM+, and Landsat OLI images captured on March 16, 1984 and December 31, 2018. The data set uses an expert system to classify each pixel to determine whether it is a water body, and the surface water is finally divided into a seasonal water body and permanent water body. In this study, the frequency of water bodies was divided into 11 intervals (>0, ≥0.05, ≥0.1, ≥0.15, ≥0.2, ≥0.25, ≥0.35, ≥0.45, ≥0.55, ≥0.65, and ≥0.75) (Figure 12). Different frequency thresholds correspond to the interannual variation of surface water. The number of permanent water pixels selected from 1989 to 2019 is consistent with the number of permanent water pixels in the JRC dataset.

Research Uncertainty
This study used Landsat images obtained from three different sensors from 1989 to 2019. Due to the small difference in surface water extraction when the spectral indices of different sensors (Landsat TM/ETM/OLI) are used, this study did not compare the differences [72,73]. Although the frequency threshold of the water body can remove most of the noise when extracting the water body, it will cause part of the water body information to be deleted and the water body area to be underestimated. The mixed pixels at the edge of the water body are the main cause of the pixels classification error. The next step can be combined with higher resolution images (e.g., Sentinel-1 or 2 A/B, Gaofen images) for water body research. This study was located in arid and semi-arid areas. The interannual and intra-annual variation of surface water in the Hetao Plain is significant, but only the interannual variation of surface water is considered in this study, and the intra-annual variation of surface water is not described. It is of great significance to understand the law of continuous surface water change in arid and semi-arid areas (taking the Hetao Plain as an example) for the study of global change, which is worthy of attention in the future.

Research Uncertainty
This study used Landsat images obtained from three different sensors from 1989 to 2019. Due to the small difference in surface water extraction when the spectral indices of different sensors (Landsat TM/ ETM/ OLI) are used, this study did not compare the differences [72,73]. Although the frequency threshold of the water body can remove most of the noise when extracting the water body, it will cause part of the water body information to be deleted and the water body area to be underestimated. The mixed pixels at the edge of the water body are the main cause of the pixels classification error. The next step can be combined with higher resolution images (e.g., Sentinel-1 or 2 A/B, Gaofen images) for water body research. This study was located in arid and semi-arid areas. The interannual and intra-annual variation of surface water in the Hetao Plain is significant, but only the interannual variation of surface water is considered in this study, and the intra-annual variation of surface water is not described. It is of great significance to understand the law of continuous surface water change in arid and semi-arid areas (taking the Hetao Plain as an example) for the study of global change, which is worthy of attention in the future.

Conclusions
This study investigated the spatial and temporal changes of surface water in the Hetao Plain from 1989 to 2019 based on all available Landsat TM, ETM+, and OLI images in the GEE platform. Using these data, the area and number of different water body range indexes were analyzed, as we investigated the variability of water bodies in the past 32 years and determined the changing trend. The results showed that the maximum, seasonal, and annual average water body area exhibited a downward trend on the whole, while the maximum, seasonal, and permanent water body number displayed an upward trend overall. The change of surface water in arid and semi-arid areas is a comprehensive reflection of the impact of climate change and human activities on regional water resources. Precipitation had statistically significant positive effects on the water body area and water body number, and it had statistically significant negative effects with temperature and irrigation. The Hetao Plain is an important commodity grain base in China, while agriculture is the pillar industry of the region's economy. The dynamic change of surface water has an important impact on agricultural development. It is essential to understand the area, number, and influencing factors of surface water changes to formulate reasonable plans in the fields of water resources management, agricultural irrigation, and ecological protection to improve the sustainable use of water resources and achieve economic and social sustainability.
Supplementary Materials: The following are available online at http://www.mdpi.com/2073-4441/12/11/3010/s1, Figure S1: The slope value in the Hetao Plain. Figure S2: The time series of the permanent water body area and annual cumulative precipitation in the Hetao Plain from 1989 to 2019. Figure S3: The trend of sown area in the Hetao Plain from 1989 to 2018. Figure S4: The trend of population and GDP in the Hetao Plain from 1989 to 2018. Figure S5: Changes in the land-use area in the Hetao Plain from 1990 to 2015. Figure S6: The ratio of irrigated area to the sown area in the Hetao Plain from 1989 to 2018. Table S1: The trend and correlation in the permanent water body area and precipitation in the Hetao Plain.
Author Contributions: This research was carried out under the cooperation of all authors. H.X. provided the writing ideas for the research, R.W. completed data collection, analysis and wrote the paper, and H.X., Y.Q., W.N., L.P., R.L., X.Z., X.B. and P.F. all contributed to the discussion and revision of the paper. All authors have read and agreed to the published version of the manuscript.