Dynamic Changes of Soil Erosion in the Taohe River Basin Using the RUSLE Model and Google Earth Engine

The Taohe River Basin is the largest tributary and an important water conservation area in the upper reaches of the Yellow River. In order to investigate the status of soil erosion in this region, we conducted a research of soil erosion. In our study, several parameters of the revised universal soil loss equation (RUSLE) model are extracted by using Google Earth Engine. The soil erosion modulus of the Taohe River Basin was calculated based on multi-source data, and the spatio-temporal variation characteristics of the soil erosion intensity were analyzed. The results showed the following: (1) the average soil erosion modulus of the Taohe River Basin in 2000, 2005, 2010, 2015 and 2018 were 1424, 1195, 1129, 1099 and 1124 t·ha−1·year−1, respectively, and the overall downward trend was obvious. (2) The ranges of soil erosion in the Taohe River Basin in 2000, 2005, 2010, 2015 and 2018 are basically the same—mainly with slight erosion—and the soil erosion in the middle and lower reaches was more serious. (3) When dealing with the vegetation cover factor and conservation practice factor in the RUSLE model, Google Earth Engine provided a new approach for soil erosion investigation and monitoring over a large area.


Introduction
Soil erosion can lead to riverbed siltation, trigger mountain torrents, destroy surface soil structure and reduce soil productivity, which has become the focus of environmental science, agriculture, soil and water conservation and other disciplines [1]. Serious soil erosion threatens the production and the life of human beings, which has become an environmental problem faced by human beings all over the world [2]. The problem of soil erosion in China is extremely serious. According to relevant data, the total area of soil erosion in China reached 2.73 × 10 6 km 2 in 2018, accounting for 28.80% of the total area [3]. The amount of soil erosion is much greater than the amount of soil loss tolerance [4]. Therefore, it is important to clarify the dynamic changes of soil erosion for the comprehensive management of soil and water conservation, the related construction of water conservancy projects and ecological civilization construction. The soil erosion model is a common method for the quantitative estimation of soil erosion [5]. In 1986, the United States established the revised universal soil loss equation (RUSLE) based on the universal soil loss equation (USLE) model [6,7]. Compared with other soil erosion models, the RUSLE model has the advantages of having a simple formula, few parameter requirements and high estimation accuracy, and has become a widely used quantitative estimation model of soil erosion globally [8,9]. The extraction of the RUSLE model parameters relies heavily on image processing software such as the

Study Area
The Taohe River originates from the eastern part of Xiqing Mountain in Qinghai Province, China. Most of the area is located in the southwest of Gansu Province, with the range of 101.60 •~1 04.33 • E and 34.10 •~3 6.02 • N. The Taohe River is the second largest tributary of the Yellow River, with a total length of 673 km and a River Basin area of 2.55 × 10 4 km 2 . Xizhai village and Haidianxia divide the Taohe River into upper, middle and lower reaches. The Taohe River Basin is divided into two geomorphic units: the Gannan plateau and the Longxi loess hill. The valley in the Gannan plateau area is narrow and the terrain on both sides is high and steep, while the loess hilly area is low and flat, and the elevation of the whole basin is between 1002-4866 m ( Figure 1). The Taohe River Basin is located in the hinterland of the mainland. It has obvious plateau continental climate characteristics. The Gannan plateau is dominated by natural alpine grassland meadows and mountain forests. The loess hilly areas are dominated by semi-arid sparse woodland and shrub grassland. The soil types in the upper, middle and lower reaches of the river are respectively alpine meadow soil, mountain brown soil and loess [23].
Water 2019, 11, x FOR PEER REVIEW 2 of 16 software such as the Environment for Visualizing Images (ENVI). In 2010, Google, Carnegie Mellon University and the United States Geological Survey jointly developed and established a cloud platform for processing satellite images and Earth observation data (Google Earth Engine). Compared with traditional image processing tools, Google Earth Engine has the function of mass data storage, which can quickly invoke and batch process remote sensing data [10][11][12]. Currently, it has become the most advanced cloud geographic information processing platform in the world [13]. Users can access the database online through the application programming language interface (API) provided by Python and a JavaScript platform based on web interactive development environment (IDE). It can conduct cloud computing and the visualization of Earth data over a large range and a long timescale in a short time, which can provide strong technical support for remote sensing monitoring in a large area [10]. Relevant scholars have achieved good results in vegetation monitoring [14][15][16], land use [17,18], crop classification [19,20] and water quality inversion [21,22] using Google Earth Engine. At present, there are few reports on the combination of soil erosion and Google Earth Engine. Therefore, this paper took the Taohe River Basin as the research area, tried to combine the RUSLE model with Google Earth Engine, extracted some parameters of the RUSLE model with Google Earth Engine, calculated the soil erosion modulus in 2000, 2005, 2010, 2015 and 2018 in the Taohe River Basin, and analyzed the spatial and temporal characteristics of soil erosion.

Study Area
The Taohe River originates from the eastern part of Xiqing Mountain in Qinghai Province, China. Most of the area is located in the southwest of Gansu Province, with the range of 101.60°~104.33° E and 34.10°~36.02° N. The Taohe River is the second largest tributary of the Yellow River, with a total length of 673 km and a River Basin area of 2.55 × 10 4 km 2 . Xizhai village and Haidianxia divide the Taohe River into upper, middle and lower reaches. The Taohe River Basin is divided into two geomorphic units: the Gannan plateau and the Longxi loess hill. The valley in the Gannan plateau area is narrow and the terrain on both sides is high and steep, while the loess hilly area is low and flat, and the elevation of the whole basin is between 1002-4866 m ( Figure 1). The Taohe River Basin is located in the hinterland of the mainland. It has obvious plateau continental climate characteristics. The Gannan plateau is dominated by natural alpine grassland meadows and mountain forests. The loess hilly areas are dominated by semi-arid sparse woodland and shrub grassland. The soil types in the upper, middle and lower reaches of the river are respectively alpine meadow soil, mountain brown soil and loess [23].

Data Sources
The data used in this study include (1) soil texture at 1:100,000 in the Taohe River Basin, which was derived from the Resource and Environment Data Cloud Platform (http://www.resdc.cn/); (2) soil organic matter at 1:100,000 in the Taohe River Basin, which came from Soil Science Data Center

RUSLE Model
The RUSLE model was used to estimate soil erosion in the Taohe River Basin. The calculation is as follows [7]: where A refers to the amount of soil loss in per unit time and per unit area (t·ha −1 ·year −1 ). R refers to the rainfall erosivity factor (MJ·mm·ha −1 ·h −1 ·year −1 ). K refers to the soil erodibility factor, which is the soil loss rate of specific soil rainfall erosivity per unit measured in a standard plot (t·ha·h· ha −1 MJ −1 ·mm −1 ). L and S refer to the topographic factor (dimensionless). C refers to the vegetation cover factor (dimensionless). P refers to the conservation measure factor, which includes engineering and tillage measure factors (dimensionless) [24].

Rainfall Erosivity Factor (R)
Precipitation is one of the important external forces that causes soil erosion, and it is also a dynamic index to objectively evaluate the separation and transportation of soil caused by rainfall. A method for estimating rainfall erosivity factor using annual and monthly rainfall proposed by Wischmeier is studied [6]. The calculation is as follows: where P i is monthly precipitation (mm) and P is annual precipitation (mm). According to the annual and monthly average precipitation data of 14 meteorological stations in and around the Taohe River Basin in 2000, 2005, 2010, 2015 and 2018, each station's rainfall erosivity factor was calculated. Kriging interpolation was programmed in Google Earth Engine to obtain the spatial distribution of the rainfall erosivity factor at a 30 m spatial resolution in the Taohe River Basin ( Figure 2).

Soil Erodibility Factor (K)
Soil erosion is affected by a combination of the physical and chemical properties of the soil. Soil erodibility is an important indicator of soil sensitivity to rainfall and runoff erosion. In this study, soil texture and soil organic carbon were used to estimate soil erodibility factors. The calculation is as follows [25]:

Soil Erodibility Factor (K)
Soil erosion is affected by a combination of the physical and chemical properties of the soil. Soil erodibility is an important indicator of soil sensitivity to rainfall and runoff erosion. In this study, soil texture and soil organic carbon were used to estimate soil erodibility factors. The calculation is as follows [25]: where Sand, Silt and Clay represent the sand, silt and clay content in the soil percentage (%), and C is the percentage of organic carbon content in soil (%); Sn1 = 1 − Sand/100. The greater the K value of the Water 2020, 12, 1293 5 of 15 soil erodibility factor, the more susceptible the soil to erosion; conversely, the soil is less susceptible to erosion. According to the relationship between organic matter and organic carbon (soil organic carbon content = soil organic matter content × 0.58) [26], the spatial distribution data of 1:100,000 soil organic matter in the Taohe River Basin were transformed into 1:100,000 soil organic carbon spatial distribution data. According to Formula (3), the spatial distribution of the soil erodibility factor with 30 m spatial resolution in the Taohe River Basin was calculated by programming on Google Earth Engine ( Figure 3A).
where S is the slope factor and θ is the value of the slope (°); the value of the slope can be extracted from DEM data.
The slope length factor of the Taohe River Basin was calculated by using overland flow according to Meth [29]. The calculation is as follows: where L is the slope length factor, Aout is the cell exit confluence area (m 2 ), Ain is the cell inlet confluence area (m 2 ), x is the grid size (m), NCSL is the cell non-cumulative slope length (m), m is the slope length index, and the value of m is related to the slope. If the value of the slope is less than 0.5°, m is 0.2; if the value of the slope is greater than 0.5° and less than 1.5°, m is 0.3; if the value of slope is greater than 1.5° and less than 3°, m takes a value of 0.4; if the value of slope is greater than 3°, m takes a value of 0.5. In this study, with the aid of the topographic factor calculation tool, the LS value of the topographic factor with a 30 m spatial resolution in the Taohe River Basin was obtained ( Figure 3B).

Topographic Factor (LS)
Topographic factors are called slope and slope length factors, which determine the movement state and direction of surface runoff. At a small scale, the slope length factor is generally measured by field data, but at a large scale, it is mainly obtained by DEM data.
The slope factor is calculated as follows [27,28]: where S is the slope factor and θ is the value of the slope ( • ); the value of the slope can be extracted from DEM data. The slope length factor of the Taohe River Basin was calculated by using overland flow according to Meth [29]. The calculation is as follows: where L is the slope length factor, A out is the cell exit confluence area (m 2 ), A in is the cell inlet confluence area (m 2 ), x is the grid size (m), NCSL is the cell non-cumulative slope length (m), m is the slope length index, and the value of m is related to the slope. If the value of the slope is less than 0.5 • , m is 0.2; if the value of the slope is greater than 0.5 • and less than 1.5 • , m is 0.3; if the value of slope is greater than 1.5 • and less than 3 • , m takes a value of 0.4; if the value of slope is greater than 3 • , m takes a value of 0.5. In this study, with the aid of the topographic factor calculation tool, the LS value of the topographic factor with a 30 m spatial resolution in the Taohe River Basin was obtained ( Figure 3B).

Vegetation Cover Factor (C)
Vegetation can protect the surface soil and slow down the rate of soil erosion. The normalized difference vegetation index (NDVI) is the most commonly used data source for calculating vegetation Water 2020, 12, 1293 6 of 15 cover factors. The vegetation cover factor in the Taohe River Basin was calculated by NDVI data. The calculation is as follows [30]: where a and b are the parameters that determine the NDVI-C relationship curve. The Van der Knijff's experiment shows that when a = 2 and b = 1, the C value is the most appropriate value [30]. The higher the C factor, the worse the vegetation growth. In this study, the code was compiled in Google Earth Engine to realize the mosaic of the Landsat satellite images, NDVI calculation and annual maximum synthesis of NDVI in 2000, 2005, 2010, 2015 and 2018 at the Taohe River Basin (Figure 4).
Water 2019, 11, x FOR PEER REVIEW 6 of 16 Figure 3. Spatial distribution of soil erodibility factor and topographic factor in the Taohe River Basin.

Vegetation Cover Factor (C)
Vegetation can protect the surface soil and slow down the rate of soil erosion. The normalized difference vegetation index (NDVI) is the most commonly used data source for calculating vegetation cover factors. The vegetation cover factor in the Taohe River Basin was calculated by NDVI data. The calculation is as follows [30]: where a and b are the parameters that determine the NDVI-C relationship curve. The Van der Knijff's experiment shows that when a = 2 and b = 1, the C value is the most appropriate value [30]. The higher the C factor, the worse the vegetation growth. In this study, the code was compiled in Google Earth Engine to realize the mosaic of the Landsat satellite images, NDVI calculation and annual maximum synthesis of NDVI in 2000, 2005, 2010, 2015 and 2018 at the Taohe River Basin (Figure 4). The conservation measure factor refers to the percentage of soil loss to planting down the slope after adopting soil and water conservation measures. The p value is between 0-1. If the value is 0, then the area is not affected by soil erosion; if the value is 1, the area has not been subjected to any soil or water conservation measures [24]. The field investigation of soil and water conservation measures is time-consuming and laborious. Therefore, the estimation of the p value in this study is realized by assigning values to land use types [31]. This study uses the LANDSAT/LT05/C01/T1_TOA and LANDSAT/LC08/C01/T1_RT datasets to screen out the Landsat images with the smallest cloud cover in the Taohe Table 1 [24], the spatial distribution of the conservation measure factor at a 30 m spatial resolution of the Taohe River Basin is obtained ( Figure 5).

Results
According to field data from the Soil and Water Conservation Bureau of the Gansu Provincial Department of Water Resources, the total amounts of soil erosion were 2.96 × 10 8

Results
According to field data from the Soil and Water Conservation Bureau of the Gansu Provincial Department of Water Resources, the total amounts of soil erosion were 2.96 × 10 8 t in the TaoheRiver Basin. This study concluded that the total amounts of soil erosion in the Taohe (Figure 6), and the areas and percentage of each erosion intensity level over the whole study area was statistically calculated (Table 3; Table 4). Table 2. Classification standard of soil erosion intensity grades [32].

Soil Erosion Intensity Grades Soil Erosion Modulus t/(km 2 ·a)
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 spatial distributions of soil erosion intensity grades in 2000, 2010 and 2018 were analyzed to clarify the characteristics of soil erosion intensity transfer (Figure 7). From 2000 to 2010, the soil erosion intensity grade in the Taohe River Basin remained unchanged at 23,153.41 km 2 , accounting for 90.75% of the total area. The area with a rise in the soil erosion intensity grade was 36.99 km 2 , accounting for 0.14% of the total area. The area with a decline in the soil erosion intensity grade was 2322.73 km 2 , accounting for 9.10% of the total area, which is mainly distributed in the middle and downstream areas. From 2010 to 2018, the soil erosion intensity grade in the Taohe River Basin remained unchanged at 24,475.79 km 2 , accounting for 95.93% of the total area. The area with a rise in the soil erosion intensity grade was 509.43 km 2 , accounting for 2.00% of the total area. It is mainly  (Table 3; Table 4).
The spatial distributions of soil erosion intensity grades in 2000, 2010 and 2018 were analyzed to clarify the characteristics of soil erosion intensity transfer (Figure 7). From 2000 to 2010, the soil erosion intensity grade in the Taohe River Basin remained unchanged at 23,153.41 km 2 , accounting for 90.75% of the total area. The area with a rise in the soil erosion intensity grade was 36.99 km 2 , accounting for 0.14% of the total area. The area with a decline in the soil erosion intensity grade was 2322.73 km 2 , accounting for 9.10% of the total area, which is mainly distributed in the middle and downstream areas. From 2010 to 2018, the soil erosion intensity grade in the Taohe River Basin remained unchanged at 24,475.79 km 2 , accounting for 95.93% of the total area. The area with a rise in the soil erosion intensity grade was 509.43 km 2 , accounting for 2.00% of the total area. It is mainly distributed in the lower reaches (East bank of the lower reaches of the Taohe River), but also sporadically distributed in the upper reaches. The area with a decline in soil erosion intensity grade was 527.91 km 2 , accounting for 2.07% of the total area. It is mainly distributed in the lower reaches (west bank of the lower reaches of the Taohe River) and the middle reaches.

Time Variation of Soil Erosion Intensity
The average soil erosion modulus of the Taohe

Time Variation of Soil Erosion Intensity
The average soil erosion modulus of the Taohe (Figure 8C). The higher the precipitation erosion factor was, the more serious the soil erosion was. The precipitation erosivity factor showed a downward trend from 2000 to 2015, and it increased significantly in 2015-2018 (the northwestern region had an abnormally abundant year in 2018 [33]), which is consistent with the change of soil erosion modulus in the Taohe River Basin. Therefore, the change of precipitation is the direct cause of the change of soil erosion in the Taohe River Basin. The areas of different erosion intensity grades were counted, and a transfer matrix of the soil erosion intensity grade was made. From 2000 to 2005, the stability rates of soil micro erosion, mild erosion, moderate erosion, strong erosion, pole strong erosion and violent erosion in the Taohe River Basin were 99.92%, 81.43%, 73.24%, 58.74%, 62.41% and 67.27%, respectively (Table 5). From 2005 to 2010, the stability rates of soil micro erosion, mild erosion, moderate erosion, strong erosion, pole strong erosion and violent erosion in the Taohe River Basin were 99.34%, 91.20%, 86.26%, 78.17%, 80.11% and 81.18%, respectively (Table 6). From 2010 to 2015, the stability rates of soil micro erosion, mild erosion, moderate erosion, strong erosion, pole strong erosion and violent erosion in the Taohe River Basin were 99.16%, 90.31%, 85.90%, 79.34%, 82.18% and 86.37%, respectively (Table 7). From 2015 to 2018, the stability rates of soil micro erosion, mild erosion, moderate erosion, strong erosion, pole strong erosion and violent erosion in the Taohe River Basin were 98.63%, 91.63%, 87.32%, 79.45%, 82.68% and 87.80%, respectively (Table 8). From 2000 to 2018, the stability rates of soil micro erosion, mild erosion, moderate erosion, strong erosion, pole strong erosion and violent erosion in the Taohe River Basin were 95.57%, 80.81%, 73.70%, 60.02%, 69.47% and 99.50%, respectively (Table 9). On the whole, the area change of violent soil erosion was the least in the Taohe River Basin. The area of strong soil erosion changed was largest in the Taohe River Basin.   (Table 9). On the whole, the area change of violent soil erosion was the least in the Taohe River Basin. The area of strong soil erosion changed was largest in the Taohe River Basin. In the Taohe River Basin from 2000 to 2018, the proportions of the area transferred from micro erosion to mild erosion, moderate erosion, strong erosion, pole strong erosion and violent erosion were 0.14%, 0.00%, 0.00%, 0.00% and 0.00%, respectively. Only an area of 24.48 km 2 changed from micro erosion to mild erosion. The proportions of the area which transferred from mild erosion to micro erosion, moderate erosion, strong erosion, pole strong erosion and violent erosion were 21.45%, 0.23%, 0.01%, 0.00% and 0.00%, respectively. An area of 827.42 km 2 changed from mild erosion to micro erosion, and an area of 8.89 km 2 changed from mild erosion to moderate erosion. The proportions of the area which transferred from moderate erosion to micro erosion, mild erosion, strong erosion, pole strong erosion and violent erosion were 0.36%, 34.70%, 0.13%, 0.00% and 0.01%, respectively. An area of 673.17 km 2 changed from moderate erosion to mild erosion. The proportions of the area which transferred from strong erosion to micro erosion, mild erosion, moderate erosion, pole strong erosion and violent erosion were 0.08%, 2.31%, 50.09%, 0.14% and 0.02%, respectively. An area of 394.90 km 2 changed from strong erosion to moderate erosion, and an area of 18.24 km 2 changed from strong erosion to mild erosion. The proportions of the area which transferred from pole strong erosion to micro erosion, mild erosion, moderate erosion, strong erosion and violent erosion were 0.07%, 0.22%, 7.81%, 42.38% and 0.07%, respectively. An area of 234.92 km 2 changed from pole strong erosion to strong erosion, and an area of 43.31 km 2 changed from pole strong erosion to moderate erosion. The proportions of the area which transferred from violent erosion to micro erosion, mild erosion, moderate erosion, strong erosion and pole strong erosion were 0.40%, 0.07%, 0.25%, 3.70%, and 40.77%, respectively. An area of 119.06 km 2 changed from violent erosion to pole strong erosion. From 2000 to 2018, the unchanged area of soil erosion intensity grade in the Taohe River Basin was 23,140.47 km 2 , accounting for 90.70%; the weakened area of soil erosion intensity grade was 2333.26 km 2 , accounting for 9.15%; the intensified area of soil erosion intensity grade was 39.40 km 2 , accounting for 0.15%; and the overall soil erosion intensity grade was slightly weakened (Table 9).

Conclusions
In this study, Google Earth Engine and the RUSLE models were used to calculate the soil erosion modulus of the Taohe  and lower reaches of the Taohe River Basin belong to the loess soil, and the vegetation growth is weaker than that of the upper reaches, so the soil erosion in the middle and lower reaches is more serious. (3) From 2000 to 2018, the stability rates of micro erosion, mild erosion, moderate erosion, strong erosion, pole strong erosion and violent erosion in the Taohe River Basin were 95.57%, 80.81%, 73.70%, 60.02%, 69.47% and 99.50%, respectively. The area change of violent erosion in the Taohe River Basin was the smallest and the area change of strong erosion was the largest. (4) From 2000 to 2018, the area with an unchanged soil erosion intensity grade accounts for 90.70%; the area with a weakened soil erosion intensity grade accounts for 9.15%; and the area with an intensified soil erosion intensity grade accounts for 0.15%. On the whole, the soil erosion intensity grade is slightly weakened. (5) Google Earth Engine has a strong advantage in dealing with the vegetation cover factor and Conservation measure factor of the RUSLE models. This will provide strong technical support for the analysis of dynamic changes in soil erosion in large areas in the future.
In this paper, Google Earth Engine and the RUSLE model are combined for the first time to study the spatial and temporal changes of soil erosion in 2000, 2005, 2010, 2015 and 2018 and to find a new method for the data processing and acquisition required by the RUSLE model.