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

Abstract

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 modified Normalized Difference Water Index (mNDWI), Enhanced Vegetation Index (EVI), and Normalized Difference 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 significant upward trend. (2) The variation of the surface water area in the Hetao Plain is mainly affected by the maximum water body area, while the variation of the water body number is mainly affected by the number of minimum water bodies. (3) Precipitation has statistically significant positive effects on the water body area and water body number, which has statistically significant negative effects with temperature and irrigation. The findings 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 ecological protection.

1. 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 area in China. Ecological protection and high-quality development of the Yellow River Basin have become national strategies [6]. The Hetao Plain, as one of the main areas of the Yellow River Basin, is an important commercial grain base in China. Irrigation directly consumed groundwater and surface water resources. In the process of the gradual development of agriculture, the contradiction between the supply and demand for water has become increasingly prominent. Monitoring the dynamic change of water resources in the Hetao Plain and improving the efficiency of water resources are very important to the economic development and agricultural production in the Hetao Plain. In the past hundred years, due to the intensification of human activities and global climate change, surface water resources have generally experienced problems such as a sharp decline in the water volume, deterioration of the water quality, and area shrinkage, which have changed the form of the surface water hydrological cycle and hydrological dynamics [2,3]. These problems have caused huge losses to the life and agricultural production of the residents in the Hetao Plain and greatly affected the development of the local economy. Problems such as the shrinkage of the surface water area, fragmentation of the water landscape, and deterioration of the ecological environment have aroused widespread concern from all walks of life [7].

The dynamic monitoring of surface water is a key issue in research on resources and environmental change [8,9]. The timely and accurate acquisition of surface water change information is an important foundation for water resource protection, utilization, and sustainable development. Remote sensing images data have the characteristics of wide coverage, high revisit frequency, rich information, and low cost; these with different time-space spectra are used in land use and land cover mapping [10]. The satellite data commonly used for water mapping and change monitoring include the Advanced Very High-Resolution Radiometer (AVHRR) [11], Moderate Resolution Imaging Spectroradiometer (MODIS) [12,13,14,15], Systeme Probatoire d’Observation de la Terre (SPOT) [16], Visible Infrared Imaging Radiometer Suite Visible Infrared Imaging Radiometer (VIIRS) [17], Landsat [4,8,15,18,19,20,21,22,23,24,25], and Sentinel [26,27], among others. Landsat images data are characterized by long time series (1972 to present), high spatial resolution, strong spectral consistency, and free access; Landsat images data represent the main data source for monitoring the dynamic changes of water bodies with long time series [21,28,29,30]. Previous studies have frequently selected a single remote sensing image in a year or in different years at the same time to represent the surface water area and interannual variation [8,9,31,32]. Due to the strong interannual and intra-annual variation characteristics of surface water [8,9], an annual or interannual single image, which was selected to characterize the area and time series of surface water, will lead to inaccuracies of the surface water body area and water body number [31,33,34,35,36]. Therefore, the use of all existing long-term sequence images is essential for accurately characterizing the interannual and intra-annual changes in surface water.

Surface water extraction algorithms based on remote sensing data can be divided into two categories: machine learning algorithms and traditional algorithms [37]. One type is machine learning algorithms, including Support Vector Machine (SVM) [38], Random Forest (RF) [39], Decision Tree (DT) [40], Deep Learning (DL) [41], etc. Machine learning algorithms require professional knowledge on sample selection and algorithm training, and it is difficult to quickly map large regional scales such as the world and a country [42]. Traditional algorithms based on a combination of bands, including the single-band method and multi-band method, can be divided into the inter-spectral relation method and water index method [43]. The multi-band method can make full use of the spectral difference between water and non-water, and it is better than the single-band method in water extraction [44]. The water index method is a commonly used multi-band extraction method, which is combined with other remote sensing indexes to select the threshold of significant difference between a water body and a non-water body. It has high accuracy in water extraction. At the same time, this method has the advantages of rapid extraction, high accuracy, simplicity, robustness, and good repeatability when extracting large areas of water [9,25,34,42,45,46,47,48].

Since 2008, Petabyte-Scale Landsat images have been shared free of charge, making it possible for researchers to use all available Landsat images to better capture surface water change information. 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.

2. Data Sets

2.1. 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², 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.

2.2. Data Sources

2.2.1. 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] (Figure 2b). The number of images obtained from different Landsat sensors in 1989–2019 (Figure 2e) varies greatly between years. The number of images obtained since the launch of Landsat ETM+ in 1999 and Landsat OLI in 2013 has increased, but only Landsat ETM+ was available in 2012.

2.2.2. 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).

2.2.3. 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.

2.2.4. 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).

3. Research Methods of Surface Water

3.1. 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.

$$mNDWI=\frac{\left({\rho }_{green}-{\rho }_{swir1}\right)}{\left({\rho }_{green}+{\rho }_{swir1}\right)}$$

(1)

$$NDVI=\frac{\left({\rho }_{nir}-{\rho }_{red}\right)}{\left({\rho }_{nir}+{\rho }_{red}\right)}$$

(2)

$$EVI=2.5\times \frac{{\rho }_{nir}-{\rho }_{red}}{{\rho }_{nir}+6\times {\rho }_{red}-7.5\times {\rho }_{blue}+1}$$

(3)

$$Water=\left\{\begin{array}{c}1,EVI<0.1\mathit{and}(\mathit{mNDWI}>\mathit{EVI}\mathrm{or}\mathit{mNDWI}>\mathit{NDVI})\\ 0,\mathit{Other}\mathit{conditions}\end{array}\right\}$$

(4)

Here, ${\rho }_{\mathit{red}}$, ${\rho }_{green}$, ${\rho }_{\mathit{blue}}$, ${\rho }_{nir}$, and ${\rho }_{swir1}$ are the reflectance of the red band (0.63–0.69 μm), green band (0.52–0.6 μm), blue band (0.45–0.52 μm), near-infrared band 1 (0.77–0.9 μm), and shortwave infrared band 1 (1.55–1.75 μm), respectively.

3.2. 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].

$$WF=\frac{N_{Water}}{N_{\mathit{Good}}}$$

(5)

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_{\mathit{Good}}$ represents the number of high-quality images during the year.

3.3. Surface Water Area and Dynamic Degree

3.3.1. 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.

$$S={\displaystyle \sum _{\mathit{i}=1}^{\mathit{N}}\mathit{A}}\times \mathit{pixelsize}\times {10}^{\mathit{-}6}$$

(6)

Here, S is the area of the maximum, seasonal, and permanent water body of different types (km²), N is the total number of pixels, A is the number of pixels, and pixelsize is the pixel area (m²).

3.3.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)) [68]

$$K=\frac{s_a-s_b}{s_a}\times \frac{1}{T}\times 100\mathit{\%}$$

(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.

3.4. 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).

4. Result

4.1. 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 was employed for an accurate evaluation of this study.

	Sentinel-2 MSI
Water Body Map (2019)	Water	No Water	Total	User Accuracy (%)
Water	952	55	1007	94.54%
No-Water	69	1924	1993	96.54%
Total	1021	1979	3000	Overall Accuracy = 95.9%
Producer Accuracy (%)	93.24%	97.22%		Kappa Coefficient = 0.91

, 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.

4.2. 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 5e shows the number of water pixels observed in eight intervals of the water body frequency from 1989 to 2019. From 1989 to 2019, the number of water pixels with a frequency range of 0.95–1 is the largest. The number of water pixels in the frequency range of 0.25–0.35 is more in 2008 and 2012 than in other years, and the water frequency in other regions is roughly the same.

4.3. 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², which is −26% to 65% of the average value (1107 km²). The permanent water body area is 346–757 km², which is −33% to 46% of the average value (519 km²), 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², which is −49% to 97% of the average value (588 km²). The annual average water body area is 598–1225 km², which is −28% to 47% of the average value (831 km²). 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² 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² = 0.3763), seasonal water bodies (R² = 0.2188), and permanent water bodies (R² = 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² 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² 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, 2010–2014, and 2015–2019 were compared (

Table 2. The dynamic attitude of the maximum water body numbers from 1989 to 2019.

Year	Number	Number Change	Dynamic Index	Overall Dynamic Change
1989–1994	300,169	-	-	45.75%
1995–1999	328,336	28,167	9.38%
2000–2004	403,831	75,495	22.99%
2005–2009	464,241	6041	14.96%
2010–2014	420,587	−43,654	−9.4%
2015–2019	453,496	32,909	7.82%

and

Table 3. The dynamic attitude of the maximum water body area from 1989 to 2019.

Year	Area(km2)	Area Change	Dynamic Index	Overall Dynamic Change
1989–1994	1246.471	-	-	−12.09%
1995–1999	1122.563	−123.908	−9.94%
2000–2004	935.649	−186.914	−16.65%
2005–2009	1070.912	135.263	14.46%
2010–2014	1178.781	107.869	10.07%
2015–2019	1060.578	−118.203	−10.03%

). 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², whilst the minimum value appeared in 2000–2004, and its value was 935.649 km². 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%.

4.4. 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² 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², which was 1011 km² more than that in 2011 (816 km²). Water body areas of more than 50 km² 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.

4.5. 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. 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² 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
Variables	Area	Number	Area	Number	Area	Number
	Coef.	Coef.	Coef.	Coef.	Coef.	Coef.
AP	2.899	0.976	2.016	0.752	0.785	0.225
AT
Irrigation	−8.155		−6.269		−1.864
Constant	1164.863	126.219	674.052	79.172	489.671	46.479
Model summary
R2	0.51	0.391	0.463	0.392	0.245	0.23
SEE	175	76	142	59	81	26
F	14.065	18.627	11.64	18.665	4.384	8.661
Sig.	0.000	0.000	0.000	0.000	0.022	0.006

. 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.

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²) to 2019 (12,410.8 km²), and it showed a continuous rise at a rate of 217.91 km²/year (Figure 11). The sown area in the Hetao Plain increased from 1989 (8725.8 km²) to 2019 (20,388.6 km²), and showed a continuous rise at a rate of 415.09 km²/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.

5. Discussion

5.1. 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.

5.2. 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.

6. 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.