Resources and Ecology in the Qinghai-Tibet Plateau

Changes of Soil Erosion and Possible Impacts from Ecosystem Recovery in the Three-River Headwaters Region, Qinghai, China from 2000 to 2015

  • WANG Zhao 1, 2 ,
  • WANG Junbang , 1, *
  • 1.Key Laboratory of Ecosystem Network Observation and Modeling, Institute of Geographic Sciences and Natural Resources Research, Chinese Academy of Sciences, Beijing 100101, China
  • 2.School of Earth Science and Resources, China University of Geosciences (Beijing), Beijing 100083, China
* FU Gang, E-mail:

Received date: 2019-04-30

  Accepted date: 2019-07-20

  Online published: 2019-10-11

Supported by

National Key Research and Development Program of China(2016YFC0500203)

Science and Technology Program of Qinghai Province(2018-ZJ-T09)

Science and Technology Program of Qinghai Province(2017-SF-A6)


Copyright reserved © 2019


Soil erosion poses a great threat to the sustainability of the ecological environment and the harmonious development of human well-being. The revised universal soil loss equation (RUSLE) was used to quantify soil erosion in the Three-River Headwaters region (TRH), Qinghai, China from 2000 to 2015. The possible effects of an ecosystem restoration project on soil erosion were explored against the background of climatic changes in the study area. The model was validated with on-ground observations and showed a satisfactory performance, with a multiple correlation coefficient of 0.62 from the linear regression between the estimations and observations. The soil erosion modulus in 2010-2015 increased 6.2%, but decreased 1.2% compared with those in the periods of 2000-2005 and 2005-2010, respectively. Based on the method of overlay analysis, the interannual change of the estimated soil erosion was dominated by climate (about 64%), specifically by precipitation, rather than by vegetation coverage (about 34%). Despite some uncertainties in the model and data, this study quantified the relative contribution of ecological restoration under global climatic change; meanwhile the complexity, labor-intensiveness and long-range character of ecological restoration projects have to be recognized. On-ground observations over the long-term, further parameterization, and data inputs with higher quality are necessary and essential for decreasing the uncertainties in the estimations.

Cite this article

WANG Zhao , WANG Junbang . Changes of Soil Erosion and Possible Impacts from Ecosystem Recovery in the Three-River Headwaters Region, Qinghai, China from 2000 to 2015[J]. Journal of Resources and Ecology, 2019 , 10(5) : 461 -471 . DOI: 10.5814/j.issn.1674-764X.2019.05.001

1 Introduction

Soil erosion refers to the processes of destruction, erosion, handling and deposition of soil and its parent materials under the external forces of hydraulic pressure, wind, freezing and thawing, and gravity (Shi, 1998). Soil erosion by water or wind is a major threat to humanity, with significant effects on land degradation, agricultural production, ecosystem services and direct social economic losses (García-Ruiz et al., 2017; Zerihun et al., 2018; Zhou et al., 2018). Increasing soil erosion poses a great threat to soil structure, agricultural production, water quality and the environment, and it has become one of the major global environmental problems (Qiao et al., 2002). At present, research on soil erosion has drawn a great deal of attention in a variety of disciplines, such as agronomy, soil science, hydrology, and environmental science (Vrieling et al., 2010).
Soil erosion has been studied through methods that can be mainly summarized as modelling, rainfall simulation and erosion plots, though at least 17 different methods have been used to assess and quantify soil erosion (Rodrigo-Comino, 2018). The revised universal soil loss equation (RUSLE) has become a widely used quantitative model for soil erosion globally due to its simple calculation formula and low data requirements (Jabbar, 2003; Qin et al., 2009; Chen, 2011; Zeng et al., 2011; Li et al., 2012; Yi et al., 2015). The model is generally combined with GIS and RS to overcome its limitations in terms of reliability, to allow for spatial coverage of large areas, and also to improve accuracy (Zerihun et al., 2018). However, the RUSLE should be parameterized according to local soil erosivity (Zhang et al., 2003).
Three-River Headwaters (TRH) of the Yangtze River, the Yellow River and the Lancang River is located in the hinterland of the Qinghai-Tibet Plateau. It is considered the water tower for China, and even for Asia (Immerzeel et al., 2010). Being sensitive and vulnerable to climate change and human activities, it is considered as the ecological security barrier; and the status and trends of its ecosystem have attracted a great deal of attention (Liu et al., 2018). Much research has examined soil erosion and its driving forces, including climate change, ecosystem protection and recovery programs (Liu et al., 2005; Qi et al., 2007; Wu et al., 2009; Wu et al., 2010; Cao et al., 2018).
However, research on the impact of ecological engineering on soil erosion changes in TRH is surprisingly lacking. In this study, the RUSLE model was selected, and through a combination of model simulation and GIS spatial analysis, the soil erosion was quantitatively simulated in order to analyze the soil erosion that occurred before and after the implementation of the ecological project. In addition, we analyzed the characteristics and causes of soil erosion changes, with the objective of evaluating the effects of ecological engineering project implementation. This analysis can provide a scientific basis for the rolling implementation and long-term planning of ecological engineering projects, as well as the prevention and control of soil erosion in TRH.

2 Materials and methods

2.1 Site description

The Three-River Headwaters region (TRH) (89°24ʹ-102°23ʹ E and 31°39ʹ-36°16ʹ N) is located in the hinterland of the Qinghai-Tibet Plateau, and the southern part of Qinghai Province, China (Fig. 1). The total land area of TRH is about 3.63×105 km2. Grassland is the main ecosystem type in TRH, covering 65.4% (Xu et al., 2008) to 69.53% of the total area. The main grassland vegetation types include alpine steppe (26%), alpine meadow (23%) and alpine desert steppe (20%). The terrain from southeast to northwest gradually rises from an altitude of 3500 m to 4200 m above sea level. It has a typical plateau continental climate, characterized by alternating hot and cold seasons and distinct wet and dry seasons (Chen et al., 2003). According to the 50-year averages of 18 meteorological stations in the region, the annual mean temperature is below 0° and the annual total precipitation is below 400 mm.

Fig. 1 The location of the Three-River Headwaters region and its ecosystem types

To mitigate grassland degradation in the region, the government implemented the ecological protection and recovery plan. It has finished the first stage from 2005 to 2012, and the second stage has been under implementation since 2013. One of its major planning goals is to reduce soil erosion by 1.13948×107 m3. To compare the temporal and spatial changes of soil erosion, the study period was divided into the three phases; that is, Phase I: Before the implementation of the project (2000-2005), Phase II: Pre-project implementation (2005-2010) and Phase III: Late-project implementation (2010-2015).

2.2 Model

2.2.1 RUSLE model
The RUSLE quantifies soil erosion by water from hill slopes through a linear equation (Renard, 2008):
$A=100\times R\times K\times LS\times C\times P$
where A is the annual average soil erosion modulus (t km-2 yr-1); R is the rainfall erosivity factor (MJ mm (km2 h yr)-1); K is the soil erodibility factor (t km2 h (km2 MJ mm)-1); LS is the slope length gradient factor (dimensionless); C is the surface vegetation cover and management factor (dimensionless); and P is the soil and water conservation measure factor (dimensionless). Multiplication by 100 is used to convert the units from t ha-1 yr-1 to t km-2 yr-1.
(1) Rainfall erosivity factor R
The rainfall erosivity factor (R) reflects the potential erosion due to rainfall factors on soil, and it is the main driving force of soil erosion. We use the rainfall observation data from meteorological stations distributed in the TRH to determine the average monthly rainfall and annual average rainfall of the study area. The monthly scale calculation formula is (Wischmeier et al., 1965):
$R=17.02\times \underset{i=1}{\overset{12}{\mathop \sum }}\,1.735\times {{10}^{(1.5\lg \frac{P_{i}^{2}}{P}-0.8188)}}$
where, P is the average annual rainfall (mm); and Pi is the average monthly rainfall (mm).
(2) Soil erodibility factor K
Soil is the main component of soil erosion, and the soil erodibility factor is an indicator of the sensitivity of erosion to soil properties. Different soil types have different K values, and there are many methods for its estimation (Wischmeier et al., 1971). In this study, according to Williams et al. (Williams et al., 1990), the soil erodibility factor was calculated from the silt, sand and clay percentages of the soil texture and the soil organic carbon content:
$\begin{matrix} & \text{ }\!\!~\!\!\text{ }K=0.1317\times \left\{ 0.2+0.3\exp \left[ -0.0256\times SAN\times \left( 1-\frac{SIL}{100} \right) \right] \right\}\times \\ & \ \ \ {{(\ \frac{SIL}{CLA-SIL})}^{0.3}}\times (1-\frac{0.25C}{C+\exp 3.72-2.95C})\times \\ & \ \ \ \ \left( 1-\frac{0.7SN1}{SN1+\text{exp}(-5.51+22.9SN1)} \right)\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \\ \end{matrix}$
where, 0.1317 is the conversion factor which converted the US unit system to the international units; SAN, SIL, CLA, and C are the percentages of sand (0.050-2.000 mm), silt (0.002-0.050 mm), clay (<0.002 mm), and organic carbon content (%); and SN1=1- SAN/100.
(3) Slope length factor LS
Topography is the basic natural geographical factor affecting soil erosion. It affects the formation and development of soil and vegetation, restricts the redistribution of material and energy on the surface, and determines the state and direction of movement of surface runoff. In the RUSLE model, the slope and slope length are indicators for measuring the influence of topography on soil erosion. The greater the slope, the greater the runoff energy, and the stronger the erosion ability on the slope. Likewise, the larger the slope length, the larger the runoff, and the stronger the erosion. In the evaluation of soil erosion on the slope scale, the slope length is generally obtained by field measurements, but on the basin scale, the slope length can be extracted by DEM.
The slope length factor L is calculated as follows:
$L={{(\lambda /22.13)}^{\alpha }}$
$\alpha =\beta /(\beta +1)$
$\beta ={\frac{\sin \theta }{0.0896}}/{(3.0\sin {{\theta }^{0.8}}+0.56)}\;$
$\lambda =cellsize\times \sqrt{flowdir}\times \cos \theta$
Among the variables, $\lambda $ is the horizontal slope length (m), $\alpha $ is the slope length index, $\beta $ is the ratio of rill erosion and surface erosion, 22.13 is the slope length (m) of the standard plot, $\theta $ is the slope extracted with DEM data, $cellsize$ is the resolution of the unit raster (m), and flow direction is the water flow direction code, which uses 1-8 to represent eight different directions.
The DEM data was used to fill in the depression and to calculate the flow direction and water accumulation of the non-depression. The runoff path was modified according to the slope ratio (the ratio of the elevation difference and distance between a grid in the flow direction and its adjacent grid), and $\lambda $ was extracted according to the revised runoff path. Slope factor S adopts the following formula established in the 1980s (Mccool et al., 1987):
$S=\left\{ \begin{matrix} & 10.8\times \sin \theta +0.03\ \ \ \ \ \ \ \ \ \ (\theta <5{}^\circ ) \\ & 16.8\times \sin \theta -0.50\ \ \begin{matrix} {} & {} \\ \end{matrix}(5{}^\circ \le \theta <10{}^\circ ) \\ \end{matrix} \right. $
Under steep slope conditions (θ>10o), the following formula was applied to calculate the S factor as proposed by Liu et al. (Liu et al., 2000):
$S=20.91\times \sin \theta -0.9\begin{matrix} {} & {} \\ \end{matrix}(\theta \ge 10{}^\circ )$
(4) Vegetation cover factor C
The vegetation cover factor (C, dimensionless) is defined as the ratio of the amount of soil lost in a certain vegetation type to the amount of soil lost in a reference plot under the same conditions of soil, slope and rainfall. At the regional scale, remotely sensed NDVI is used to estimate factor C (Asis et al., 2007; Zhou et al., 2008; Kouli et al., 2009; Markose et al., 2016) according the following method (Cai et al., 2000):
$C=\text{exp}\left[ -\alpha \times \frac{NDVI}{(\beta -NDVI)} \right]$
where, NDVI is the normalized vegetation index averaged from the 8-day data in the growing season (defined from June to August). The α and β coefficients are given values of 2.5 and 1.0, respectively (Cao et al., 2018).
(5) Soil and water conservation factor P
P reflects the maintenance of soil erosion by man-made engineering measures. Its value ranges from 0-1. The smaller the value, the better the effect of the manual measures. The larger the value, the worse the effect of the manual measures. When P=0, there is no erosion; and when P=1, there are no measures to control soil erosion. In studies of large-scale watershed soil erosion, P is usually determined by combining field survey data to assign different land use types. Combining the literature values with the land use types of the study area, the P value of agricultural land is 0.15; water, wetland, bare rock, and ice and snow have no erosion, so their P values are 0; and the remaining land use types basically do not have any water conservation measures, so their P values are 1 (Zhao et al., 2007; Fu et al., 2008; Li and Zheng, 2012; Sun et al., 2013).
(6) The total amount of soil erosion
The factors are stacked and multiplied by a unit conversion coefficient of 100. The obtained A is the soil erosion modulus (t km-2 yr-1), and the total amount of soil erosion (QSE, t yr-1) is given by the formula:
$QSE=\underset{i=1}{\overset{n}{\mathop \sum }}\,A\times {{M}_{i}}$
In the formula, A is the amount of soil erosion per unit area (t km-2 yr-1), Mi is the area (km2) of the i pixel. The size of the pixel in this study is 250 m × 250 m.
2.2.2 Classification of soil strength
According to the Standards for classification and gradation of soil erosion (SL 190-96) in China, implemented in 1997, the erosion intensity of the study area is divided into 6 grades (Table 1).
Table 1 The standards for classification and gradation of soil erosion (SL 190-96)
Class Soil erosion modulus (t km-2 yr-1)
Micro erosion <500
Mild erosion 500-2500
Moderate erosion 2500-5000
Serious erosion 5000-8000
Most serious erosion 8000-15000
Intense erosion >15000
2.2.3 Trend analysis
The trend analysis was applied to explore the impacts from vegetation and climate changes on soil erosion. The trend was calculated from the slope of a linear regression through a least squares algorithm:
$\beta =\frac{n\sum\limits_{i=1}^{n}{i\times Y}-\sum\limits_{i=1}^{n}{i}\sum\limits_{i=1}^{n}{{{Y}_{i}}}}{n\times \sum\limits_{i=1}^{n}{{{i}^{2}}}-{{\left( \sum\limits_{i=1}^{n}{i} \right)}^{2}}}$
where, β is the interannual trend of the dependent variable Yi as soil erosion rate, NDVI, or precipitation in the research period of n years, and i is the year from 2000 to 2015.

2.3 Data

In this study, the satellite based NDVI, the interpolated meteorological grid, and the soil data were used as inputs, and their sources are summarized in Table 2.
Table 2 Data and data sources
Data Dataset name Reference/website
DEM SRTM of 90 m resolution
Meteorological Data Daily observations of meteorological stations China meteorological science data sharing service network
Soil texture Harmonized World Soil Database (HWSD) v1.2 Harmonized World Soil Database(
Vegetation index MODIS data products MOD09Q1 V006
Soil classification National soil database of China at 1:100,000 scale Resource and Environment Data Cloud Platform (
Land cover and use China Cover Resource and Environment Data Cloud Platform (
Vegetation classification National vegetation map of China at 1:100,000 scale Resource and Environment Data Cloud Platform (
2.3.1 NDVI data from the MODIS product
In this study, the NDVI was calculated from the red and near-infrared band of MODIS 8-day 250 m Surface Reflectance products (MOD09) for the THR from 2000 to 2015. The adaptive Savitzky-Golay (SG) method in TIMESAT software (Jönsson et al., 2004) was applied to remove the noise in the NDVI timeseries, then the NDVI data were used to calculate the vegetation coverage factor C in equation 10.
2.3.2 Meteorological data
In this study, the meteorological data of 250 m resolution were interpolated with the observations from the 836 stations in China and the 345 stations in surrounding countries The meteorological variables include daily total precipitation (mm), daily minimum (Tmin, °C), maximum (Tmax, °C) and average (Tavg, °C) temperatures, daily relative humidity
(RHU, %), and sunshine hours (SH, hour). The software of ANUSPLIN version 4.2, an algorithm based on the multi-variate thin plate smoothing splines (Hutchinson et al., 2000; Gill et al., 2002), was used to produce the meteorological data of 250-m spatial resolution and 8-day time step. The stations in the surrounding countries were used to produce the interpolated area, then the research area was clipped as an area that is smaller than the interpolated area in order to remove boundary effects, especially in the THR area where there are fewer stations. Compared to the observations on the flux towers of ChinaFLUX and AsiaFLUX, the interpolated data of 1 km resolution showed much higher consistency with the averaged multiple correlation coefficient (R2) above 0.95 for the temperature data, and 0.70 for the precipitation data (Wang et al., 2017). Using software Matlab R2014b, we calculated average annual rainfall and monthly rainfall, which were used to calculate the rainfall erosivity factor R in equation (2).
The digital elevation model (DEM) data of 250 m resolution was resampled from the 90 m Shuttle Radar Topographic Mission (SRTM) database v4.1 (http://srtm.csi.cgiar. org/SELECTION/inputCoord.asp) and used as input for the spatial interpolation of the meteorological data. The DEM data were used to fill in the depression and calculate the flow direction and water accumulation of the non-depression area. Then, the runoff path was modified according to the slope ratio (the ratio of elevation difference and distance between a grid in the flow direction and its adjacent grid), and the slope length factor LS in equation (4)-(9) was calculated.
2.3.3 Soil data
The soil texture and nutrient contents were from the Harmonized World Soil Database (HWSD v1.2), which includes the percent of sand (0.050-2.000 mm), silt (0.002- 0.050 mm) and clay (<0.002 mm), and organic carbon content. In the database, soil maps of China have been compiled from the Second National Soil Survey of China at a scale of 1:1 million (FAO et al., 2012).
With the support of ArcGIS 10.2 and the soil map of Qinghai Province, soil data were extracted from the national soil classification map of China at 1:100000 scale, which was downloaded from the Chinese Academy of Sciences data center. Then, sand, silt, and clay percentages in the topsoil were extracted from the HWSD database, resampled to a resolution of 250 m, and used to calculate the soil erodibility factor K in equation (3).
2.3.4 Validation data
The soil erosion data sampled in 2010 were collected from the 24 stations distributed over the whole THR area. Erosion pegs were applied in the measurements taken in September or October in each year and the observed soil erosion was calculated by the following formula (Shao et al., 2018):
${{A}_{obs}}=Z\times S\times \text{cos}\theta /1000$
${{M}_{obs}}=1000\times A\times \gamma /Sq$
Where, Aobs and Mobs are the field measured soil erosion (m3) and soil erosion rate (t km-2 yr-1), respectively, Z is erosion thickness (mm), θ is hill slope, γ is soil bulk density (t m-3), and Sq is horizontal projected area (m2). The details of the sampling methods and station information can be found in Shao et al. (Shao et al., 2018).

3 Results and analysis

3.1 Validation of estimated soil erosion

Estimated data were compared to the observed data from the 24 stations monitoring soil and water conservation in TRH in 2010 (Shao et al., 2016). The results showed that the estimation was significantly linearly correlated with the observations and the correlation coefficient is 0.617 (significance level, P = 0.002, sample number, n = 24), which means that the estimations can explain about the 62% of the spatial variance in the observations (Fig. 2).

Fig. 2 The verification results of the soil erosion modulus

Compared to other published results (Table 3), the approach used for estimation in this study was more rational.
Table 3 Estimated soil erosion in the published literature for the TRH area
Citation Method Average Value
(t km-2 yr-1)
Standard Deviation
(t km-2 yr-1)
(t km-2 yr-1)
(t km-2 yr-1)
Total amount of soil
erosion (×108 t yr-1)
Cao et al. (2018) Used RUSLE Model to Simulate 1140 3.2
Lin et al. (2017) Used RUSLE Model to Simulate 8946 31
Shao et al. (2011) 137Cs tracing method 574.4 487.974 1614 4
Li et al. (2009) 137Cs and 210Pbex tracing methods 1152 527.102 1623 453
Liu et al. (2005) Used RUSLE Model to Simulate 313 0.823
This study Used RUSLE Model to Simulate 2018.716 3213.14 77305.95 0 7.21
In previous studies, the validation data was very limited and generally only used data from a few stations. In our study, the 24 stations used were distributed over the whole THR, and so the validation is more representative of the spatial pattern.

3.2 Spatial and temporal patterns of soil erosion

Mild erosion was dominant, with the soil erosion modulus of 2018.716 t km-2 yr-1 and the total amount of erosion of 7.21 × 108 t yr-1 over the whole TRH in the period from 2000 to 2015 (Fig. 3). According to the Standards for classification and gradation of soil erosion (SL 190-96), the areas of mild, middle and strong gradation accounted for 74.8%, 12.36% and 11.53% of the total area, and the severe erosion zone was about 1.16%.

Fig. 3 The annual anomaly of soil erosion modulus (a), annual NDVI (b) and annual precipitation (c) from 2000 to 2015.

Compared with that of the pre-project (the period of 2000-2005), the soil erosion modulus (intensity) increased to 144.1 t km-2 yr-1 (7.4%) and 120.2 t km-2 yr-1 (6.2%) in the early stage (the period of 2005-2010) and the latter stage (the period of 2010-2015), respectively. The areas of the different erosion gradations changed less with mild erosion decreasing from 23.6% in 2000-2005 to 21.6% in 2005-2010, and then increasing to 22.7% in 2010-2015, and the severely eroded areas increased from 0.92% in 2000-2005 to 1.06% in 2005-2010 and then to 1.16% in 2010-2015 (Fig. 4).

Fig. 4 Spatial patterns of soil erosion intensity classification in the periods of 2000-2005 (a), and 2010-2015 (b), and the changes of soil erosion modulus in 2010-2015 compared to 2000-2005 (c) and 2005-2010 (d) in Three-River Headwaters.

3.3 Soil erosion in different ecosystems

Grasslands are the principal contributors of soil erosion, with their higher soil erosion modulus and wider coverage (Table 4). An increase of NDVI will inhibit the increase of the soil erosion modulus, and the amount of precipitation is positively correlated with the amount of the soil erosion modulus. The shrub ecosystem had the highest precipitation, reaching 602 mm, while the precipitation in the alpine desert steppe was lower, at 361 mm. The NDVI in the forest and shrub ecosystems are higher, reaching 0.458 and 0.391, respectively, but their soil erosion moduli are still high, reaching 3160 t km-2 yr-1 and 4007 t km-2 yr-1, respectively. These levels might be due to the precipitation, which reached 562 mm and 602 mm, respectively. Compared with other types, the alpine steppe has a lower NDVI of 0.216, and its precipitation is 414 mm, which is slightly higher than the average precipitation in TRH. Its soil erosion modulus was the same as the soil erosion modulus of the whole region, and its area is about 28% of the total area, so the total quantity of soil erosion is relatively large, accounting for 26% of the whole TRH.
Table 4 Soil erosion in different ecosystems
Ecosystem NDVI Precipitation (mm) Soil erosion modulus (t km-2 yr-1) Quantity of soil erosion (104 t yr-1) Area (km²)
Alpine meadow 0.34 549.03 3390.65 28150.68 83024.38
Alpine steppe 0.22 414.70 2195.42 20461.01 93198.50
Forest 0.46 562.40 3160.82 258.50 817.81
Shrub 0.39 602.58 4007.36 6067.26 15140.31
Farmland 0.30 444.96 136.19 7.49 550.63
Alpine desert steppe 0.16 361.98 2069.57 14934.26 72161.31
Wetland 0.30 490.16 0 0 21420.13
Waterland 0 312.55 0 0 11279.13
Desert 0.11 306.91 1097.43 2218.23 20212.88
Others 0.10 441.78 5.52 21.76 39447.56
Total 0.26 529.19 2018.72 72119.20 357252.63
The soil erosion in grassland was 87% of the total erosion averaged across 2000-2015, which included contributions of 39%, 28% and 21% by Alpine meadow, Alpine steppe and Alpine desert steppe, respectively. Except Alpine meadow, all the other ecosystems played roles in decreasing soil erosion. For example, in alpine steppe, the soil erosion modulus decreased from 10% in 2005-2010 to 6% in 2010-2015 (Fig. 5). In Alpine meadow, the problem appeared less optimistic, since its soil erosion modulus increased from 5% in 2005-2010 to 7% in 2010-2015.

Fig. 5 Three stages of soil erosion modulus and total soil erosion in different ecosystems

Note: AM: Alpine meadow ecosystem; AS: Alpine steppe ecosystem; FO: Forest ecosystem; SH: Shrub ecosystem; FA: Farmland ecosystem; SV: Alpine desert steppe; WE: Wetland ecosystem; WA: Waterland ecosystem; DE: Desert ecosystem; OT: Other ecosystem.

The implementation of the “Phase I Project” has slowed the trend of increasing soil erosion in various ecosystems to some extent (Fig. 5). Even in phase III, the erosion amounts of alpine grassland, forest and sparse vegetation ecosystems decreased significantly compared with phase II. In the future, soil erosion protection measures should still focus on the alpine meadow ecosystem.

3.4 Impacts of the changes in vegetation and climate

The Fig. 6 shows the interannual trends of annual total precipitation, NDVI and the soil erosion modulus from 2000 to 2015 in the study area (Fig. 6). In the overall spatial pattern, these three trends showed insignificantly increasing fluctuations. Overall, the NDVI increased insignificantly at the rate of 0.0005 yr-1 (R2=0.1001, P=0.316) from 2000-2015 in TRH. The NDVI increased in the northeast and central areas, and decreased in the central and western TRH. The precipitation increased insignificantly at a rate of 4.4296 mm yr-1 (R2=0.1427, P=0.377), with a contrasting spatial pattern that was increasing over the southeast and central areas and decreasing over the southwest area. For the rate of change of the erosion modulus, it increased significantly at the rate of 14.72 t km2 yr-1 (R2=0.0774, P=0.278), and the amount of erosion in the central and eastern regions has gradually increased, while the erosion in parts of the south and the northeast have decreased.

Fig. 6 The interannual trends of NDVI (a), precipitation (b) and the soil erosion modulus (c) from 2000 to 2015.

In order to further explore the impacts of vegetation and precipitation changes on soil erosion, the interannual trends were overlaid for the soil erosion modulus, precipitation and NDVI; and there are eight combinations of spatial patterns (Fig. 7) and Table 5 shows the overlaid trends and area percentages.

Fig. 7 Overlay of the trends of the soil erosion modulus, precipitation, and NDVI in different regions.

Note: 1. Increasing precipitation, increasing NDVI and increasing erosion;

2. Increasing precipitation, decreasing NDVI, and increasing erosion;

3. Increasing precipitation, increasing NDVI and decreasing erosion;

4. Decreasing precipitation, increasing NDVI and decreasing erosion;

5. Decreasing precipitation, increasing NDVI and increasing erosion;

6. Decreasing precipitation, decreasing NDVI and increasing erosion;

7. Decreasing precipitation, decreasing NDVI and decreasing erosion;

8. Increasing precipitation, decreasing NDVI and decreasing erosion.

Table 5 The overlaid trends and their area percentages of the NDVI, precipitation and soil erosion in the Three-River Headwaters for 2000-2015
Types Change rate
Change rate of
Change rate of
soil erosion
1 0.0010 4.5620 15.0537 64.380
2 -0.0009 4.3570 39.4090 25.151
3 0.0026 4.2974 -12.1285 9.442
4 0.0011 -0.4644 -5.2669 0.470
5 0.0003 -0.3066 1.2424 0.252
6 -0.0008 -0.4505 3.5909 0.255
7 -0.0001 -0.7922 -0.4779 0.004
8 -0.0002 1.1022 -5.9287 0.046

Note: The trends numbered 1 to 8 have the same definitions given in the legend of Fig. 7.

Theoretically, the changes in precipitation will result in a positive effect and vegetation coverage change will have a negative effect on the soil erosion modulus. As shown in Table 4, for the area with the increasing soil erosion (90.04%), it was driven predominantly by the increasing precipitation (64.38%) and decreasing vegetation coverage (25.15%). For the areas with the decreasing soil erosion (9.96%), it should reflect the negative effects of the increasing vegetation coverage, rather than the positive effects of increasing precipitation.
Overall, the increase in the soil erosion modulus in the TRH is mainly due to the significant increase in rainfall erosivity (about 64% area), followed by decreasing vegetation coverage (25%), while the decline in the soil erosion modulus benefits from a significant increase in vegetation coverage in the region (about 9% area).

4 Discussion

4.1 Uncertainty

This study applied a RUSLE model to estimate soil erosion in the Three-River Headwater. The major uncertainties may come from both the model itself and the input data. On the one hand, the RUSLE model was originally used for site- specific soil erosion estimation, and in the model localization process, so using a RUSLE model to simulate spatial patterns of soil erosion may have some uncertainties. On the other hand, in order to facilitate the calculation, we determined the P factor by combining field survey data in order to assign different land use types. The P value of agricultural land is 0.15; water, wetland, bare rock and ice and snow have no erosion, so the P value is 0; the remaining land use types basically do not employ any water conservation measures, so their P values are 1. This method may not fully reflect the actual soil and water conservation measures of each grid, so the results may have some uncertainties.
However, the estimation can explain about 62% of the spatial variability in the observations (Fig. 2), and was very consistent with the values from the literature (Table 3), which means the estimation was satisfactory despite some uncertainties. In the future, parameterization would be applied to improve the model based on more observations.
For input data, the interpolated precipitation data was applied. Due to the limited number of meteorological observation stations that were available, and their uneven spatial distribution, the data is expected to introduce some uncertainties into the estimation. But the validation showed that the interpolated precipitation was significantly correlated with the observations from the ecological stations with the linear correlation equation of y=0.99x-2.06 (R2=0.78, P< 0.001). The interpolated mean air temperature and total precipitation could explain 99% and 78% of variability in the observations, respectively (Wang et al., 2017). In the future, by the fusion of multiple sources of data, such as satellite-based potential atmospheric precipitation (or TRMM), the precipitation data could be improved and the uncertainty could be decreased in the estimation of the soil erosion.

4.2 Effects of ecological projects

The soil erosion modulus of the study area showed an increasing trend under global climatic changes, though the ecological projects still play an important role in the prevention and control of soil erosion in the study area (Shao et al., 2016; Cao et al., 2018). One key problem for understanding the further implications of the project is quantifying the contributions of the project and climate changes to soil erosion. As a fundamental step, however, the contributions from precipitation and vegetation coverage should first be quantified.
Based on our method of overlay analysis, the contribution from precipitation played a dominant role compared with the changes in vegetation coverage. According to the area impacted, increasing precipitation results in more serious soil erosion in 64% of the area, while the changes of soil erosion in over 34% of the area would be attributed to vegetation changes (with 34% being the sum of the areas of 25% decreasing and 9% increasing vegetation coverage). This conclusion differs from that reached by the analysis of controlling factors, in which the changes of precipitation contributed a positive effect of 180% and an effect of -80% from the increasing vegetation coverage (Cao et al., 2018).
However, this method ignores the effect of changing precipitation on vegetation coverage. Generally, increasing precipitation induces an increase of vegetation coverage. The changes of precipitation will result in a positive effect if only this single factor is considered, although it is more complex if both precipitation and temperature was considered (Wang et al., 2003; Wang et al., 2004; Wang et al., 2009). Therefore, it is a difficult challenge to assess the effects of ecological restoration under climate change, and more appropriate methodologies will be necessary and essential for improving the results.

4.3 Soil erosion in different vegetation types

The intensity of overall soil erosion is determined by the two aspects of vegetation and precipitation. In alpine desert steppe, there is less annual precipitation but vegetation coverage is still lower. Therefore, due to the sparse vegetation, the relatively weak rainfall can easily lead to relatively strong soil erosion in the alpine desert steppe. In this case, significantly reduced soil erosion could be attributed to the implementation of the ecosystem protection and recovery project.
In contrast, the erosion intensity in forestland ecosystems is lower than in grasslands due to the roots of forest being stronger and deeper than those of grass (Sun et al., 2013). In addition, the alpine steppe and alpine meadow are considered to have strong stability and self-recovery ability thanks to their higher coverage and dense root systems (Zhang et al., 2016). Due to the greater annual precipitation, however, there are higher soil erosion intensities in those areas (Wang et al., 2003, Wang et al., 2004). From the viewpoint of soil erosion, there is a great challenge for the implementation of ecosystem protection and recovery projects, and they point to the need for further soil conservation efforts in alpine steppe and meadow in the future.

5 Conclusions

The revised universal soil loss equation (RUSLE) was used to estimate soil erosion in the Three-River Headwaters region (TRH), Qinghai, China, from 2000 to 2015. We tried to quantify the effects of an ecosystem protection and restoration project on soil erosion under the background of climate changes in the study area. According to the validation, the estimations could explain about 62% of the spatial variance in the observations. From 2000 to 2015, the average soil erosion modulus of the TRH was 2018.716 t km-2 yr-1 with a total amount of erosion of 7.21 × 108 t yr-1. The soil erosion showed a fluctuating increasing trend mainly due to the significant increase in rainfall erosivity (about 64% of the area), followed by decreasing vegetation coverage (25%), while the declining soil erosion modulus may contribute to a significantly increasing vegetation coverage in the region (about 9% area). The implementation of the ecological restoration project may mitigate the increasing trend of increasing precipitation to some extent. The ecosystem restoration remains a great challenge under global climatic changes. On the methodology, more on-ground observations are necessary and essential for parameterizing the model and decreasing uncertainty in the estimation.
Asis A M D, Omasa K . 2007. Estimation of vegetation parameter for modeling soil erosion using linear Spectral Mixture Analysis of Landsat ETM data. ISPRS Journal of Photogrammetry & Remote Sensing, 62(4):309-324.

Cai C F, Ding S W, Shi Z H , et al. 2000. Study of applying USLE and geographical information system IDRISI to predict soil erosion in small watershed. Journal of Soil and Water Conservation, 14(2):19-24. (in Chinese)

Cao W, Liu L L, Dan W U . 2018. Soil erosion changes and driving factors in the Three-River Headwaters region. Acta Prataculturae Sinica, 27(6):10-22. (in Chinese)

Chen G C, Lu X F, Peng M , et al. 2003. Basic characteristics and protection of ecosystem in Sanjiangyuan region of Qinghai province. Qinghai Science and Technology, 10(4):14-17. (in Chinese)

Chen Z F . 2011. Research on Soil Loss Equation of ChongQing City based on Revised Universal Soil Loss Equation (RUSLE). Southwest University. (in Chinese)

FAO, IIASA, ISRIC, et al. 2012. Harmonized World Soil Database (version 1.2). FAO and IIASA.

Fu S F, Zha X . 2008. Study on predicting soil erosion in Dongzhen Watershed based on GIS and USLE. Journal of Geo-Information Science, 10(3):390-395. (in Chinese)

García-Ruiz J M, Beguería S, Lana-Renault N , et al. 2017. Ongoing and emerging questions in water erosion studies. Land Degradation & Development, 28(1):5-21.

Gill R, Kelly R, Parton W , et al. 2002. Using simple environmental variables to estimate below-ground productivity in grasslands. Global Ecology and Biogeography, 11(1):79-86.

Hutchinson G, McIntosh P . 2000. A case study of integrated risk assessment mapping in the Southland Region of New Zealand. Environmental Toxicology and Chemistry, 19(4):1143-1147.

Immerzeel W W, Van Beek L P H, Bierkens M F P . 2010. Climate change will affect the Asian Water Towers. Science, 328(5984):1382-1385.

Jabbar M T . 2003. Application of GIS to Estimate Soil Erosion Using RUSLE. Geo-Spatial Information Science, 6(1):34-37.

Jönsson P, Eklundh L . 2004. TIMESAT—a program for analyzing time- series of satellite sensor data. Computers & Geosciences, 30(8):833-845.

Kouli M, Soupios P, Vallianatos F . 2009. Soil erosion prediction using the Revised Universal Soil Loss Equation (RUSLE) in a GIS framework, Chania, Northwestern Crete, Greece. Environmental Geology, 57(3):483-497.

Li T H, Zheng L N . 2012. Soil erosion changes in the Yanhe Watershed from 2001 to 2010 based on RUSLE model. Journal of Natural Resources, 27(7):1164-1175. (in Chinese)

Liu B Y, Nearing M A, Shi P J , et al. 2000. Slope length effects on soil loss for steep slopes. Soil Science Society of America Journal, 64(5):1759-1763.

Liu M C, Li D Q, Wen Y M , et al. 2005. The spatial analysis of soil retention function in Sanjiangyuan region and its value evaluation. China Environmental Science, 25(5):627-631. (in Chinese)

Liu S, Zamanian K, Schleuss P M , et al. 2018. Degradation of Tibetan grasslands: Consequences for carbon and nutrient cycles. Agriculture Ecosystems & Environment, 252(5):93-104. (in Chinese)

Markose V J, Jayappa K S . 2016. Soil loss estimation and prioritization of sub-watersheds of Kali River basin, Karnataka, India, using RUSLE and GIS. Environmental Monitoring & Assessment, 188(4):1-16.

Mccool D K, Brown L C, Foster G R , et al. 1987. Revised slope steepness factor for the universal soil loss equation. Transactions of the ASAE - American Society of Agricultural Engineers (USA), 30(5):1387-1396.

Qi Y G, Zhang W, Zhang X Z . 2007. Practice and exploration of soil and water conservation ecological monitoring in the three-rivers’ headstream region. Soil and Water Conservation in China, (11):23-25. (in Chinese)

Qiao Y, Yun Q . 2002. Fast soil erosion investigation and dynamic analysis in the loess plateau of China by using information composite technique. Advances in Space Research, 29(1):85-88. (in Chinese)

Qin W, Zhu Q, Zhang Y . 2009. Soil erosion assessment of small watershed in Loess Plateau based on GIS and RUSLE. Transactions of the Chinese Society of Agricultural Engineering, 25(8):157-163. (in Chinese)

Renard A L . 2008. Poverty and charity in Saudi Arabia: the royal family, the private sector and the welfare state. International Critic, 41(4):137-156. (in French)

Rodrigo-Comino J . 2018. Five decades of soil erosion research in “terroir”. The State-of-the-Art. Earth-Science Reviews, 179:436-447.

Shao Q Q, Fan J W, Liu J Y. 2018. Monitoring and assessment on the effects of the ecological conservation and restoration project in the Three-River Source Region. Science Press. (in Chinese)

Shao Q Q, Fan J W, Liu J Y , et al. 2016. Assessment on the effects of the first-stage ecological conservation and restoration project in Sanjiangyuan region. Acta Geographica Sinica, 71(1):3-20. (in Chinese)

Shi D M . 1998. Discussion on some terminology of soil and water conservation. Journal of Soil and Water Conservation, 4(4):90. (in Chinese)

Sun J, Cheng G, Li W , et al. 2013. On the variation of NDVI with the principal climatic elements in the Tibetan Plateau. Remote Sensing, 5(4):1894-1911. (in Chinese)

Vrieling A, Sterk G, Vigiak O . 2010. Spatial evaluation of soil erosion risk in the West Usambara Mountains, Tanzania. Land Degradation & Development, 17(3):301-319.

Wang G, Ding Y, Wang J , et al. 2004. Land ecological changes and evolutional patterns in the source regions of the Yangtze and Yellow Rivers in recent 15 years. Acta Geographica Sinica, 59(2):163-173. (in Chinese)

Wang G X, Sheng Y P, Qian J , et al. 2003. Study on the influence of vegetation change on soil moisture cycle in Alpine Meadow. Journal of Glaciology and Geocryology, 25(6):653-659. (in Chinese)

Wang J B, Liu J Y, Shao Q Q , et al. 2009. Spatial-temporal patterns of net primary productivity for 1988-2004 based on glopem-cevsa model in the “three-river headwaters” region of Qinghai Province, China. Chinese Journal of Plant Ecology, 33(2):254-269. (in Chinese)

Wang J B, Wang J W, Ye H , et al. 2017. An interpolated temperature and precipitation dataset at 1-km grid resolution in China (2000 - 2012). China Scientific Data, 2(1):72-80.

Williams J R, Sharpley A N . 1990. EPIC-erosion/productivity impact calculator: 2. User manual. Technical Bulletin - United States Department of Agriculture, 4(4):206-207.

Wischmeier W H, Johnson C B, Cross B V . 1971. Soil erodibility nomograph for farmland and construction sites. Journal of Soil & Water Conservation, 26(5):189-193.

Wischmeier W H, Smith D D . 1965. Predicting rainfall-erosion losses from cropland east of the Rocky Mountains. Agricultural Handbook, 282.

Wu W, Liu F, Chen Q , et al. 2009. Study on soil erosion types in Three-river Headstream Region. Journal of Earth Sciences & Environment, 31(4):423-426. (in Chinese)

Wu W Z, Liu F G . 2010. Dynamic analysis and distribution characteristics of freeze-thaw erosion in the Three-rivers Headstream Region. Journal of Qinghai Normal University, 26(1):57-61. (in Chinese)

Xu X L, Liu J Y, Shao Q Q , et al. 2008. The dynamic changes of ecosystem spatial pattern and structure in the Three-River Headwaters region in Qinghai Province during recent 30 years. Geographical Research, 27(4):829-838. (in Chinese)

Yi K, Wang S Y, Wang X , et al. 2015. The characteristics of spatial-temporal differentiation of soil erosion based on RUSLE model: A case study of Chaoyang City, Liaoning Province. Scientia Geographica Sinica, 35(3):365-372. (in Chinese)

Zeng L Y, Wang M H, Li C M . 2011. Study on soil erosion and its spatio-temporal change at Hongfeng Lake watershed based on RUSLE model. Hydrogeology & Engineering Geology, 38(2):113-118. (in Chinese)

Zerihun M, Mohammedyasin M S, Sewnet D , et al. 2018. Assessment of soil erosion using RUSLE, GIS and remote sensing in NW Ethiopia. Geoderma Regional, 12:83-90.

Zhang H Y, Fan J W, Shao Q Q , et al. 2016. Ecosystem dynamics in the Returning Rangeland to Grassland Program. Acta Prataculturae Sinica, 25(4):1-15. (in Chinese)

Zhang W, Xie Y, Liu B . 2003. Spatial distribution of rainfall erosivity in China. Journal of Mountain Science, 21(1):33-40. (in Chinese)

Zhao L, Yuan G L, Zhang Y , et al. 2007. The amount of soil erosion in Baoxiang Watershed of Dianchi Lake based on GIS and USLE. Bulletin of Soil and Water Conservation, 27(4):42-46. (in Chinese)

Zhou B B, Chen X P, Wang Q J , et al. 2018. Effects of nano carbon on soil erosion and nutrient loss in a semi-arid loess region of Northwestern China. International Journal of Agricultural and Biological Engineering, 11(1):138-145.

Zhou P, Luukkanen O, Tokola T , et al. 2008. Effect of vegetation cover on soil erosion in a mountainous watershed. Catena, 75(3):319-325.