Tracking Changing Evidence of Water Erosion in Ordos Plateau, China Using the Google Earth Engine

Abstract

Ordos Plateau is one of the primary sources of sediment in the Yellow River, and changes in regional soil erosion directly affect the ecological status of the lower reaches of the Yellow River. Many recent studies have been published using remote sensing (RS) and geographic information systems (GIS) to evaluate soil erosion. In this study, much satellite remote sensing data in the Google Earth Engine (GEE) can better track soil erosion protection, which is significant in guiding the ecological protection and restoration of the Ordos Plateau and the Yellow River basin. In this study, we used GEE to observe the changes in soil erosion in the Ordos area from 2013 to 2021. The Theil–Sen procedure and Mann–Kendall significance test methods were used to evaluate the trend of land erosion in the Ordos area from 2013 to 2021. Based on GEE, the RUSLE is applied to evaluate soil erosion and analyze the changing trend. As a result, (1) we found that the annual change of soil and water loss in the Ordos Plateau showed three stages: 2013–2015, 2016–2018, and 2018–2021. After 2018, soil loss decreased from 14 × 10¹⁷ Mg in 2018 to 4 × 10¹⁷ Mg in 2021, which indicates that the environmental restoration project started in 2018 has achieved encouraging results. (2) The results showed that 40.9% of the regional soil erosion trend showed a significant decline, and 45.7% of the regional soil erosion trend showed a slight decline. Only 13.3% of the regional soil erosion is on the rise. (3) The test results of different land use types show that 87.3% of soil erosion occurs in bare and cultivated land. Because the terrain of Ordos is relatively flat, 95.39–96.17% of soil erosion occurs in areas with a slope of 0 to 5. (4) The reliability of the RUSLE model based on the GEE platform is proved by regression model verification of observation data and model prediction results. (5) GEE’s cloud-based features can provide data and scripts to users in developing countries which lack sufficient observation data or the necessary computing resources to develop these data. The results show that GEE has robust analysis and processing ability, can analyze a large amount of data, and can provide efficient digital products for soil erosion research.

1. Introduction

Soil and water loss is a widespread environmental problem in the world, which seriously endangers the survival and sustainable development of human beings [1]. The Yellow River basin is China’s region with the most serious soil erosion. Because of the severe soil erosion, the Yellow River has become the river with the largest sediment concentration in the world, and even the suspended river phenomenon has been formed in the middle and lower reaches. The Ordos Plateau is one of the primary sediment sources in the Yellow River basin. Early mineral exploitation and cultivated land reclamation caused a large amount of sediment loss.

At present, the commonly used hydrological models are HSPF model, AGNPS model, L-THIA model, SWAT model, and so on. The SWAT model is one of the most widely used distributed hydrological models in water resource analysis and evaluation. When predicting soil erosion modulus, the Water Erosion Prediction Project (WEPP) is often used [2,3,4,5]. It is a standard research method to quantitatively estimate soil erosion using the soil erosion model [6,7,8]. The United States established the revised universal soil erosion equation (RUSLE) based on the universal soil erosion equation (USLE) [9,10,11]. RUSLE is one of the widely used soil erosion equations and has the advantages of simple structure, readily available parameters, and more straightforward calculation [12,13]. However, using RUSLE depends on image processing software, and the data collection and processing workload is enormous. Therefore, it is difficult to use this model for long-term soil and water loss monitoring. With the development of technology, the appearance of Google Earth Engine [14], a cloud geographic information processing platform, makes it possible to monitor soil erosion over a long time.

The satellite observation data accumulated by GEE can help us better understand soil erosion changes. This paper selected the typical Ordos Plateau, which is strongly disturbed by human activities, as the research area. Based on the platform of Google Earth Engine, long-term satellite observation is used to review the changes in soil erosion and monitor the effects of ecological restoration and protection. Google Earth Engine has a massive data storage function compared to traditional image processing tools, which can quickly call and process remote sensing data in batches. GEE can realize fast, multi-method and batch processing of images by programming in JavaScript and Python. In recent years, many scholars have applied GEE to water surface dynamic change monitoring, agricultural yield estimation, urban development, and so on [14]. Relevant scholars have made exemplary achievements in vegetation monitoring [15,16] using Google Earth Engine for land use [17,18], crop classification [18,19,20] , and water quality inversion [19,20], and its application for soil erosion monitoring is still in the preliminary exploration stage.

In this paper, we show the results of 9 years of long-term observation of soil and water loss from 2013 to 2021 and discuss the changing trend and pattern of soil and water loss in the Ordos Plateau. In this study, a long-term observation case study using Google Earth Engine was put forward to provide strategic information for soil and water loss restoration. The primary purpose of this study is to explore the deployment of end-to-end applications on GEE and take the Ordos as an example to conduct a fully automatic and practical soil loss assessment.

2. Study Area

The Ordos area is located in the southwest of the Inner Mongolia Autonomous Region and the hinterland of the Ordos Plateau. It is about 400 km long from east to west and 340 km wide from north to south, with a total area of 87,000 square kilometers. Its geographical location is shown in Figure 1.

The Ordos area has an average altitude of more than 1400 m, and its terrain is low in the southeast and high in the northwest. The area has a complex terrain connected to the Loess Plateau in the south and is surrounded by the Yellow River in the east, west, and north. There are five landforms: plains, plateaus, desert, hills, and sandy land.

The Ordos is a semi-arid continental climate zone in the north temperate zone, with drought and little rain. The precipitation is mainly concentrated from June to August every year. The Yellow River is 590 km long in the Ordos, and there is water all year round. The soil types are mainly aeolian sandy soil, swamp soil, saline soil, and chestnut soil. The land is mainly divided into four categories: alluvial plain, hilly and gully, desert plateau, and wavy plateau.

3. Materials and Methods

3.1. Satellite Data

This study was conducted in Ordos area, and all the data from 1 January 2013 to 31 December 2021 were analyzed to estimate the soil loss in the Ordos area. The data used in this paper are shown in

Table 1. Data list, sources, datasets, and original resolution.

	Variables	Source	Dataset	Resolution
RUSLE	Daily rainfall data	UCSB/CHG	CHIRPS	5566 m
	Soil organic carbon	Beijing Normal University	GSDE	10 km
	Sand content	Beijing Normal University	GSDE	10 km
	Silty content	Beijing Normal University	GSDE	10 km
	Clay content	Beijing Normal University	GSDE	10 km
	Digital elevation model (DEM)	University of Tokyo	MERIT	3 arc-second
	NDVI	USGS	Landsat 8	30 m
	Land cover	European Space Agency	WorldCover (2020)	10 m

.

Among them, (1) soil data come from global soil data GSDE, including soil sand, silt, clay, and organic matter content (http://poles.tpdc.ac.cn/zh-hans/data/2e46eb77-3ca2-4b90-9a42-fd49f10630d4/ accessed on 6 December 2020) [21]. (2) Topographic data come from MERIT DEM data provided by Yamazaki D et al., including processed flow direction and discharge data (https://developers.google.com/earth-engine/datasets/catalog/MERIT_Hydro_v1_0_1 accessed on 6 December 2020) [22]. (3) Rainfall data come from the CHIRPS precipitation data set provided by Funk et al., which combines 0.05 resolution satellite images with in situ station data to create grid rainfall time series (https://developers.google.com/earth-engine/datasets/catalog/UCSB-CHG_CHIRPS_DAILY accessed on 6 December 2020) [23]. (4) Land Use data come from ESA’s global 10 m resolution land use images (https://developers.google.com/earth-engine/datasets/catalog/ESA_WorldCover_v100 accessed on 6 December 2020) [24]. (5) USGS provides annual median synthetic data of Landsat8 images from 2013 to 2021 (https://developers.google.com/earth-engine/datasets/catalog/landsat accessed on 6 December 2020).

3.2. Method

3.2.1. Calculation Based on Goggle Earth Engine Platform

GEE, the world’s advanced geographic data scientific analysis platform [15], can realize the fast calculation of large-area remote sensing index in the cloud. The approach to processing and analysis is shown in Figure 2. More details can be found in our GUI application: (https://yqx0903.users.earthengine.app/view/erdssoil accessed on 15 December 2022).

According to the calculation formula of RUSLE in reference [25], the various factors were calculated using the GEE platform’s powerful cloud computing capability, and RUSLE calculation from 2013 to 2021 was completed.

A = R · K · LS · C · P

(1)

where $A$ refers to the amount of soil loss per unit time and area. R represents the rainfall erosivity factor. $K$ is the soil erodibility factor. $L$ and $S$ represent the topographic factor (dimensionless). $C$ is the vegetation cover factor (dimensionless). $P$ refers to the conservation measure factor (dimensionless). K factor and LS factor are regarded as constants. R, K, and P factor are variables.

The soil erodibility (K) factor is calculated as follows:

$$K=\left\{0.2+0.3\exp \left[-0.0256Sand\left(1-\frac{Silt}{100}\right)\right]\right\}\times {\left[\frac{Silt}{Clay+Silt}\right]}^{0.3}\times \left[1-\frac{0.25C}{C+\exp \left(3.72-2.95C\right)}\right]\times \left[1-\frac{0.7Sn1}{Sn1+\exp \left(-5.51+22.9Sn1\right)}\right].$$

(2)

where Sand, Silt and Clay represent the percentage content of sediment and sedimt and clay in the soil, respectively, and C represents the organic carbon content in the soil; Sn1 = 1-sand/100. Soil date comes from global soil data GSDE, which is uploaded to GEE storage for the next calculation.

LS factor is composed of slope factor (S) and slope length factor (L). The LS factor can be calculated in GEE by using MERIT DEM data [22,26,27,28].

To fully use the image information of the study area and overcome the influence of cloudiness, this study adopts the minimum cloud cover image synthesis method. This paper selects 9 years of Landsat 8 data from 2013 to 2021. The Landsat cloud mask algorithm is used in the GEE platform (https://code.earthengine.google.com/?Script-path=examples:datasets/Landsat_LT05_C01_T1_SR, accessed on 15 December 2022) to calculate the input Landsat-8 data set, which accords with the time and space range, removes the cloud pixels and reconstructs the composite image with the minimum cloud cover of the target year for the cloudless pixels. Furthermore, the processed Landsat data set is used to calculate the annual maximum NDVI in Google Earth Engine. The formula of C-factor is as follows:

$$C=\exp \left[-a\times \frac{NDVI}{b-NDVI}\right]$$

(3)

where a = 2, b = 1.

To obtain the annual land use images, this study used the random forest method in GEE to classify the annual LandSat8 image data in the Ordos area. The specific process is as follows: the global land use data of 10 m resolution in 2020 provided by ESA is used as the label, and using stratified sampling strategy, 1000 training points are randomly selected every year. After manual visual interpretation, these sampling points are used as the final training labels. Then, the median synthetic data of Landsat are calculated. Finally, the random forest classifier is used to train the annual images to obtain the annual land use images. In the study, we use the GEE random forest function (ee.Classifier.smileRandomForest) to conduct the land cover and land use mapping. The number of trees and min leaf population were 100 and 1, respectively. The latter parameter means that when the number of nodes becomes 1, the tree stops growing. According in

Table 2. P-values of different land use types.

Land Use Types	Trees	Shrubland	Grassland	Cropland	Built-Up	Barren	Water
P	1	1	1	0.24	0	1	0

[29], land use categories were reclassified to the P-factor.

3.2.2. Trend Analysis

In this study, the Theil–Sen procedure and Mann–Kendall significance test was used to analyze the trend of soil erosion in the Ordos from 2013 to 2021. According to the value of the Theil–Sen procedure and Mann–Kendall significance test, the changing trend of soil erosion can be divided into five categories: significant increase, slight increase, insignificant change, slight decrease, and significant decrease.

$$S_D=Median\left(\frac{D_j-D_i}{j-i}\right),$$

(4)

$$S={{\displaystyle \sum }}_{j-i}^{n-1}{\displaystyle \sum }_{j-i+1}^{n}sgn\left(D_i-D_j\right),$$

(5)

$$sgn\left(D_j-D_i\right)=\left\{\begin{array}{c}1, D_i-D_j>0\\ 0, D_i-D_j=0\\ -1, D_i-D_j<0\end{array}\right.,$$

(6)

$$s\left(S\right)=\frac{n\left(n-1\right)\left(2n+5\right)}{18},$$

(7)

$$Z=\left\{\begin{array}{c}\frac{S-1}{\sqrt{s\left(S\right)}}, S>0\\ 0, S=0\\ \frac{S+1}{\sqrt{s\left(S\right)}}, S<0\end{array}\right.,$$

(8)

where S_D indicates the trend of soil erosion from 2013 to 2021, S_D > 0 indicates the increase in soil erosion from 2013 to 2021, and S_D < 0 indicates the decrease in soil erosion from 2013 to 2021. $D_i$ and $D_j$ are the observations at $i$ and $j$. n is the length of the period. The trends are shown in

Table 3. Trend features.

SD	Z	Trend Feature
SD > 0	|Z| > 1.65	significant increase
	|Z| < 1.65	slight increase
SD = 0	Z	insignificant change
SD < 0	|Z| > 1.65	slight decrease
	|Z| < 1.65	significant decrease

.

3.2.3. Model Validation

In this study, we use the linear regression model to evaluate the accuracy of the RUSLE model [2]. The specific methods are as follows: (1) The sediment transport moduli in the Yangtze River and Yellow River basins of China are selected as the observation data. (2) The observation data and prediction data are linearly fitted. (3) The effectiveness of the soil erosion model is evaluated by using the fitting situation.

4. Results

4.1. Model Validation

We performed a regression analysis on the soil erosion intensity of the RUSLE model and used the sediment transport modulus as the observed value of soil loss [2,30,31]; the observed data are from Wang et al. [1]. The regression results are shown in Figure 3. The results show that it is reliable to estimate the soil erosion in the Ordos area by using RUSLE and its parameters (R² = 0.5812).

4.2. Soil Erosion Changes in Loess Plateau

The intensity of soil erosion can be divided into six grades from slight erosion to severe erosion, as shown in

Table 4. Classification standard of soil erosion intensity grades [32].

Soil Erosion Intensity Grades	Soil Erosion Modulus (Mg/(km2 yr))
Micro-erosion	<1000
Mild erosion	1000~2500
Moderate erosion	2500~5000
Strong erosion	5000~8000
Pole strong erosion	8000~15,000
Violent erosion	>15,000

.

The analysis of soil erosion intensity shown in Figure 4 shows that the soil erosion intensity of 32.1% of the study area decreased from 2016 to 2021, and the light erosion gradually spread to the western region (Figure 4a–d). Total soil erosion decreased from ~12 × 10¹⁷ Mg in 2016 to 4 × 10¹⁷ Mg in 2021. In different soil erosion grades, the area of violent erosion decreased significantly from 352 × 10² km² in 2016 to 49 × 10² km² in 2021. The area proportion and total amount of erosion at all levels are counted yearly to obtain Figure 4f. Figure 4f shows that in each year, the proportions of slight and mild erosion are the highest, the proportions of moderate and intensive erosion are lower, and the proportion of severe erosion fluctuates significantly. It can be seen from Figure 4e that the changing trend of erosion area for each grade is consistent with the total amount of erosion every year. Before 2016, the total amount of erosion in the region showed an overall upward trend, reaching a peak in 2016 and then began to decline, which is inseparable from the measures taken by the country in the field of soil and water conservation. The total amount of erosion in 2021 has decreased by 66.7% compared with the peak in 2016, showing remarkable achievements in controlling soil and water loss in recent years.

4.3. Staged Characteristics of Soil Erosion Changes in the Ordos

Figure 5 shows the spatial distribution of soil erosion types in different periods. From 2013 to 2021, the spatial distribution of soil and water loss in the Ordos is relatively consistent, with the south and north being mainly micro erosion and mild erosion, and the areas with severe erosion being concentrated in the northeast. From 2013 to 2016, the intensity of soil erosion in the Ordos increased yearly, from 577.9 × 10¹⁵ Mg to 1350.22 × 10¹⁵ Mg. The soil erosion of pole strong erosion and violent erosion increased by 305.7847× 10² km², accounting for 27.16% of the total area. They are mainly distributed in the western region. From 2016 to 2018, this trend remained unchanged. After 2018, the erosion intensity in the western Ordos decreased yearly. In 2018, the erosion rates of pole strong erosion and violent erosion were 227.82 × 10³ km² and 343.196 × 10³ km², respectively. However, in 2021, the erosion area of these two types of soil was only 125.32 × 10³ km² and 48.35 × 10³ km², respectively.

4.4. Change Tendency of Soil Erosion

The change trend chart shows that the central and western regions have the most significant changes in the study area (Figure 6). Therefore, regarding significant changes, the downward trend accounts for the main component, while the proportion of significant and non-significant upward trends is tiny. According to the conclusion drawn from Figure 4, the erosion intensity in the southwest has gradually declined since 2016, with 40.94% of the total pixels showing a significant decline and 45.74% of the total pixels with a slight downward trend. Therefore, the decline of soil erosion is the primary trend in the Ordos.

5. Discussion

5.1. Soil Erosion Characteristics of Different Land Use Types

Land use type is the most intuitive manifestation of human influence on soil erosion, and different land use types have different abilities to maintain the soil [33]. Therefore, there are significant differences in the possibility of soil erosion among different land types. According to the data in 2020, the main land use types in the Ordos include cultivated land, woodland, grassland, water bodies, construction land, and bare soil. Among them, the area of bare soil is the largest, accounting for 51.53%, followed by grassland and cultivated land, accounting for 33.92% and 8.23%, respectively. These three land use types account for 97% of the total land use area and are the three main land use types in the Ordos. Because there are too few samples of the other types, the central statistics are soil erosion on the three main types in bare soil, grassland, and cultivated land.

Figure 7 shows that, among the different land use types, the most intense soil erosion occurs in bare soil, followed by cultivated land. From 2013 to 2018, the erosion modulus of bare soil gradually increased and gradually returned to the 2013 level from 2018 to 2021. The erosion amount decreased from 20.37 × 10¹¹ Mg·km⁻²·yr⁻¹ to 13.36 × 10¹¹ Mg·km⁻²·yr⁻¹. The erosion modulus of grassland remained unchanged, and the changing trend of the erosion modulus of cultivated land was similar to that of bare soil, reaching its peak in 2018 and then decreasing.

5.2. Soil Erosion Influenced by Climate and Human Activities

In Section 4.4, we show the changing trend of soil erosion in the Ordos, and we can find that the trend of soil erosion in most areas of Ordos is declining to different degrees. The main reason is related to the increasing vegetation coverage year by year. Increasing vegetation coverage can reduce the soil erosion caused by rainfall to a certain extent and enhance the service capacity of soil and water conservation.

The widespread distribution of barren land and farmland in the western Ordos leads to the decline of soil erosion. Mechanically speaking, the barren land is not covered by vegetation, which leads to the direct exposure of the soil to the wind. At the same time, the farmland is seriously disturbed by human activities and production activities. Therefore, poor land and farmland are easily affected by wind, which leads to severe soil and nutrient loss, thus reducing soil erosion capacity.

Soil erosion is influenced by land use types, rainfall, and soil properties. The northwest is mainly plains and plateaus, with flat terrain and mostly woodland and grassland, so the soil conservation function in this area is mainly affected by vegetation cover.

There is a large amount of soil erosion in the central desert area, mainly due to the scarcity of vegetation and poor soil and water conservation measures, which leads to low soil conservation intensity. Therefore, the local soil and water conservation measures significantly improve the soil conservation function in areas with severe soil erosion and a bad ecological environment.

The northeast of Ordos has the highest degree of soil erosion, mainly because of the excellent climate variability. The precipitation and temperature change significantly in the northeast of Ordos. Consistent with previous studies, farmland soil erosion is high. In addition, farmland is artificially managed, which makes it challenging to control soil erosion in the farmland ecosystem.

5.3. Advantages of the Google Earth Engine Platform

In this study, the GEE platform was used to calculate the soil erosion in the Ordos from 2013 to 2020. The GEE platform can obtain a large amount of data quickly and provide powerful data processing capabilities. Compared with traditional remote sensing image processing tools (such as ENVI and ArcGIS), the GEE platform includes many encapsulation functions, such as cloud removal, cropping, mosaic, operation, and classification. Through code programming, data acquisition, processing, analysis, and result output can be carried out quickly and conveniently, significantly improving research efficiency.

The calculation process for soil erosion needs a lot of intermediate processes, and rainfall, DEM, soil, NDVI, and land use can all be collected and processed in the cloud platform. Conventional data collection and processing work is heavy, and it is challenging to calculate regional soil erosion quickly due to the influence of clouds and rain. The GEE platform can quickly collect and display data and judge the critical period of data quality selection. All the methods proposed in this study are based on the GEE remote sensing cloud platform, which eliminates the limitations of data storage and computing capacity and has great application potential.

Google Earth Engine’s online data acquisition and powerful cloud computing capabilities can effectively solve the problems of data collection difficulties, complex data preprocessing, and inefficiency of the RUSLE model. Compared with the traditional method of using GIS software to calculate RUSLE, to extract soil erosion factors, the GEE cloud computing platform can only process and calculate large-scale remote sensing data in batches by writing codes and designing user interfaces. Under the same order of magnitude of data processing, the efficiency has been dramatically improved, and the real-time visualization module can make faster decisions on the experimental results.

6. Conclusions

In this study, the temporal and spatial changes of soil erosion in the Ordos Plateau from 2013 to 2021 were calculated by GEE. The main results are as follows:

(1)

From 2013 to 2021, the trend of soil and water loss increased at first, then decreased after treatment, which was mainly divided into three stages: 2013–2015, 2106–2018, and 2018–2021. Soil erosion in 2021 was slightly lower than that in 2013, and it showed an improving trend compared with 2016 under the restoration and protection projects implemented by the government.

(2)

The results of testing the trend of land erosion in the Ordos area from 2013 to 2021 showed that 40.9% of the region experienced a significant decrease in soil erosion, while 45.7% experienced a slight decline. Only 13.3% of the region had an increasing trend in soil erosion.

(3)

Compared with the measured data, it shows that the result of soil erosion detection using Google Earth Engine is reliable. Using Google Earth Engine, we can better understand the current situation and changes in soil and water loss and identify the original or less manufactured landscapes by tracking the evidence of historical changes to provide support for environmental restoration and protection management. The long-term changes in soil erosion can be used as one of the target indicators of sustainable development, which can be used to evaluate the sustainable management of land resources and provide a basis for evidence-based decision-making.

One limitation of the validation of the model and the calculations with RUSLE is that they do not accurately represent the relationship between erosion on slopes and sediment transport in rivers. These are two separate processes, and in some cases, only a small portion of material removed from slopes ends up in river channels. The catchment area can also impact this relationship [32]. Future research work will focus on further expanding the research scale, entirely using the cloud computing power of GEE, and using various soil erosion models to comprehensively verify the study results to cope with environmental changes in different regions.