Effects of Vegetation Change on Soil Erosion by Water in Major Basins, Central Asia

Abstract

The uncertainties in soil erosion (SE) are further intensified by various factors, such as global warming, regional warming and humidification, and vegetation cover changes. Moreover, quantitative evaluations of SE in major basins of Central Asia (CA) under changing environments have rarely been conducted. This study conducted quantitative evaluation of SE in four major basins (Syr Darya Basin (SDB), Amu Darya Basin (ADB), Ili River Basin (IRB) and Tarim River Basin (TRB) using the Revised Universal Soil Loss Equation (RUSLE) and analyzed the main driving factors. SE quantities in the basins presented relatively consistent upward fluctuating trends from 1982 to 2017. Vegetation cover variation fluctuated significantly from 1982 to 2017. Specifically, vegetation cover decreased continuously in SDB, ADB, and IRB, but increased gradually in TRB. Pixels with positive spatial variation of vegetation mainly occurred around lakes and oases near rivers. The Normalized Difference Vegetation Index (NDVI) showed higher correlation with precipitation (80.5%) than with temperature (48.3%). During the study period, the area of arable land (AL) exhibited the largest change among all land use types in CA. Under long-term human activities, the proportion of NDVI of other land types converting to AL was the highest. In the structural equation model (SEM), precipitation, temperature, Shannon Diversity Index (SHDI), and NDVI strongly influenced SE. Overall, the major basins in CA were jointly affected by climate, human activities, and vegetation. Specifically, climatic factors exerted the strongest influence, followed by SHDI (human activities). SE was found to be relatively serious in ADB, SDB, and IRB, with SE in SDB even approaching that in the Loess Plateau. Under the background of global changes, appropriate water and land resource management and optimization configurations should be implemented in CA with reference to TRB in order to relieve local SE problems.

1. Introduction

The fifth climate change evaluation report of the Intergovernmental Panel on Climate Change (IPCC) pointed out that the global surface average temperature increased by 0.89 °C from 1901 to 2012 and the rate of increase in the 20th century was the highest during 1980–2010 [1]. Moreover, the rate of temperature increase has been increasing since 1980, with the highest rate in mid-latitude regions of the Northern Hemisphere [2,3]. Hence, climatic changes dominated by global warming have become one of important environmental problems in the world [4]. In global environmental change, vegetation, which is a key component of land ecosystems, is extremely sensitive to climatic changes, and dynamic changes of vegetation are often used as an important biological indicator of climatic change on Earth [5,6,7]. A previous study pointed out that the global greening of vegetated areas increased by 5% from 2000 to 2017, 1/3 of which was attributed to China and India. The growth of greening of vegetated areas in China accounted for 25% of global growth [8]. Although the greening of vegetated areas was increased significantly, vegetation areas in China and India account for only 9% of the global total vegetation area [9]. The analysis shows that afforestation and intensive agriculture contributed 42% and 32% respectively to China’s greening. Intensive agriculture accounts for 82% of India’s greening [10]. According to a previous study, the increase in vegetation cover and extension of the growing season in mid and high-latitude areas in the Northern Hemisphere was mainly attributed to climatic warming [11]. Furthermore, some studies have found a positive correlation between Normalized Difference Vegetation Index (NDVI) and temperature in high latitude areas of the Northern Hemisphere in the growing season [12,13]. Piao et al. [14] concluded that the increase in NDVI was attributable to temperature increase at the national scale, while it was related to precipitation at the regional scale. The increasing frequency of human activities since the industrial revolution has also introduced some ecological degradation problems, such as forest reduction, biodiversity reduction, and soil erosion (SE) [15,16,17]. Specifically, SE is the most serious problem [18].

SE refers to damages and loss of water-land resources and productivity of lands caused by various exogenic forces [19]. It not only causes land resource depletion, soil fertility reduction, and ecological imbalance, but also influences production development activities of various industries [20,21,22]. Climatic changes and human activities all influence the development rate of SE, and these influences are mainly reflected by the promotion of SE by global temperature rise [23,24]. As the near-surface temperature increases, the wind speed increases, which also accelerates the melting of ice and snow, and increases the flow of rivers, which in turn leads to increased SE [25]. Land-use and land-cover change (LUCC) caused by human activities and climatic changes alter the original surface vegetation types and coverage, runoff conditions, and physicochemical properties of soil, thus influencing the dynamics and anti-erosion resistance system of SE [26]. Eybergen et al. [27] studied the key process of climate sensitivity in 1989 and believed that SE was influenced by climatic changes. Kirkby et al. [28] performed model predictions and found that an increase in global temperature 2–3 °C will cause land vegetation changes at the local scale, thus causing serious SE. Furthermore, vegetation, which is the hub for material circulation and energy exchange of ecosystems, can control water and soil loss fundamentally. It is an important factor for controlling SE and the most positive and effective measure for soil and water conservation [29,30]. The 8th (2009–2013) national forest resource survey of China reported that the national forest area reached 208 million hm² and the forest cover was 21.6%; the construction of key protection forest systems through the “Three North” Shelter Forest Program and in Yangtze River Basin decreased the water and soil loss area by 110,000 km², and 40% of the water-land loss area in the engineering region could be controlled effectively [31]. Since the implementation of the policy of returning arable land (AL) into forest land (FL) and grass land (GL) in 1999, vegetation coverage in the Loess Plateau increased by 28% and water and soil loss has been inhibited effectively [32].

Central Asia (CA) is located in Eurasia’s hinterland, and it belongs to a typical temperate desert and temperate continental climate. In response to global warming, CA is experiencing a significant temperature rise [33]. Chen et al. [34] studied the characteristics of temperature change in arid regions of CA in the last 100 years (1901–2003) using empirical orthogonal functions and pointed out that the regional average rate of temperature rise was 0.18 °C per decade, reaching 0.21 °C per decade in winter. This was far higher than the global rate of temperature rise. The annual precipitation in CA has been increasing since the 20th century, but with significant spatial differences [35]. Based on the Climatic Research Unit (CRU) data, Chen et al. [36] pointed out that annual precipitation in arid regions of CA has generally been increasing in the past 80 years (1930–2009). Nevertheless, the arid regions of CA are different from other global arid regions because it features high mountain-basin structures with relatively abundant precipitation [37]. Zhang et al. [38] discussed vegetation changes in CA from 1982–2012 and found that vegetation activities were improving, with prominent improvements in mountainous areas [39]. In arid regions, SE mainly occurred on basin geomorphic units according to mountain–oasis–lake regions [40]. Influenced by climatic changes and human activities, some rivers in CA are cut off at lower reaches and lakes at lower reaches dry up gradually [41], thus intensifying local SE and exposing river basin ecosystems to more risks [42]. Although CA is generally warm and humid, such trends still cannot change the degree of droughts [43]. Therefore, SE continues to pose serious threats to CA [44]. However, water erosion by water in CA under the background of rapid global changes has rarely been investigated [45,46].

To address these issues, this study investigated SE in four major basins of CA and analyzed the major driving factors. The main objectives were to: (1) perform quantitative assessment of spatial-temporal changes in water erosion by water in major basins of CA; (2) explore environmental changes (global warming, global greening, regional warming, humidification, and LUCC) in CA and their impacts on water erosion by water in major river basins; (3) analyze the main controlling factors and regulation mechanism of water erosion by water in major basins of CA.

2. Materials and Methods

2.1. Study Area

CA (46°24′–97°06′E, 33°48′–55°48′N) covers five countries (Kazakhstan, Uzbekistan, Turkmenistan, Kyrgyzstan, and Tajikistan) and Xinjiang, China (Figure 1). It is located in Eurasia’s hinterland and belongs to the typical temperate continental climate. CA is an arid region with the largest transverse span in the Northern Hemisphere [47,48]. The annual average precipitation in CA is 17–800 mm. Extremely arid regions are located in the hinterland of Taklimakan Desert in the Tarim Basin, and regions with high precipitation are distributed in the Tianshan Mountains and Ili River Valley [49,50]. The main source of water in arid regions of CA is ice and snow melting in mountainous regions. The Tianshan Mountains are the “water tower” of CA and they serve as the major water source of CA [51].

Major river basins in CA include the Amu Darya Basin (ADB) (4.65 × 10⁵ km²), Syr Darya Basin (SDB) (2.19 × 10⁵ km²), Ili River Basin (IRB) (4.16 × 10⁵ km²), and Tarim River Basin (TRB) (8.46 × 10⁵ km²). Among them, TRB has the largest mobile desert in the world, which is the Taklimakan Desert [52]. Amu Darya and Syr Darya both originate from the Tianshan Mountains, and they are the two major rivers flowing into the Aral Sea [53,54]. The Ili River is the major inflowing river of Balkhash Lake. Global warming has accelerated ice and snow melting in mountainous areas [55]. Although this increases river runoff to some extent, the inflow volume to Balkhash Lake in the lower reaches is decreasing every year because of large-scale water diversion for irrigation in the upper reaches [56,57]. As a result, lake, and water areas in the lower reaches of CA, including the Aral Sea, Taitma Lake, and Balkhash Lake, are decreasing every year [58,59]. In arid regions, river basins support most human and social economic activities [60]. Therefore, problems such as the degradation of the ecological environment and water shortage in CA require urgent attention.

2.2. Data Sources

2.2.1. Meteorological and Soil Data

CA has extremely few meteorological stations, which are mainly distributed in oases and plains [61]. Moreover, most stations have been abandoned since the collapse of the Soviet Union. Consequently, it is very difficult to obtain long-term continuous meteorological data [62]. In this study, the CRU dataset of the University of East Anglia of East Anglia, UK (http://www.cru.uea.ac.uk/data (accessed on 14 April 2021)) was used for meteorological data. The CRU dataset has good adaptability over global areas and provides high accuracy in the characterization of drought and moisture [63,64]. Therefore, the temperature and precipitation data of CA from 1982 to 2017 were extracted from CRU TS v.4.04 (CRU Time series version 4.04), with a resolution of 0.5 × 0.5°. We resampled all data to a spatial resolution of 1 km by a statistical downscaling method in order to ensure consistency across data products.

Soil data used in this study were obtained from the world soil database (http://www.fao.org/soils-portal/en/ (accessed on 4 January 2021)). The Harmonized World Soil Database (HWSD) is the global soil database jointly developed by the Food and Agriculture Organization of the United Nations (FAO) and the International Institute for Applied Systems (IIASA). Based on the simplified American standards, HWSD divides soils into four types according to grain size: clay particles (<0.002 mm), silt (0.002–0.05 mm), sand (0.05–2 mm), and gravel (>2 mm). The spatial resolution of data is 1 km.

2.2.2. Vegetation and Land Use Data

As an important component of ecosystems, vegetation is the link between the atmosphere and pedosphere [65]. With the continuous development of remote sensing technology, many vegetation indexes have become available for determining surface vegetation cover and its spatial heterogeneity, with support for global and regional studies [66,67]. NDVI, which is acquired by the AVHRR sensor of the National Oceanic and Atmospheric Administration and MODIS satellite sensor developed by the National Aeronautics and Space Administration (NASA), is characterized by wide coverage, long time series, and better accuracy. Thus far, it has been extensively used to study dynamic changes of global or regional vegetation [68,69]. In this study, two sets of global databases were applied, namely GIMMS3g NDVI and MODIS NDVI. The spatial resolution, interval, and time range of GIMMS3g NDVI (Global Inventory Modelling and Mapping Studies) data were 0.0833°, 15 d and 1982–2015, respectively. The MODIS NDVI data were collected from the Land Processes Distributed Active Archive Center (LP DAAC/NASA) of NASA MODIS. The MODIS NDVI database used in this study was the MOD13A2 V6 product, provided by the MODIS sensor onboard the Terra Satellite. The spatial resolution and time resolution were 1 km and 16 d, respectively. The data were collected during 2000–2017. The above databases have been widely used to study long-term changes in vegetation coverage in many regions worldwide [70,71,72]. Hence, GIMMS 3g NDVI was calibrated with the variance matching method using MODIS NDVI data as the benchmark, and the spatial resolution was unified to 1 km. Finally, an NDVI dataset of the major basins in CA from 1982 to 2017 was obtained.

LUCC data were obtained from ESA CCI-LC (European Space Agency, Climate Change Initiative, Land Cover, https://www.esa-landcover-cci.org (accessed on 18 June 2021)), with a spatial resolution of 300 m. ESA CCI-LC has been extensively applied in studies on CA [73,74]. In this study, ESA CCI-LC products in 1992 and 2017 were used as the basic data and TM images of CA in 1982 were selected as the background maps. Artificial visual interpretation was carried out with reference to ESA CCI-LC in 1992. According to the interpretation results, CA was divided into 8 land use/cover types: arable land (AL), grass land (GL), forest land (FL), potential desert land (PDL), desert land (DL), water body (WB), construction land (CL), and others. High-resolution Google Earth images were used for verification and the accuracy could reach more than 86% [75,76]. This indicates that the land use/cover data after artificial visual interpretation has good applicability to CA. In this study, we resample all raster data to 1-km.

2.3. Methods

The Revised Universal Soil Loss Equation (RUSLE) model is an improved SE analysis model based on the USLE model of the Agricultural Research Service (ARS) of the United States Department of Agriculture (USDA). With its strong practicability and comprehensive abilities, the RUSLE model has been extensively used in Chinese and foreign studies [77,78]. Based on the geographic information system (GIS), SE in CA was quantified using remote sensing images as the input data. The basic formula of the RUSLE model is:

$$SE=R\times K\times LS\times C\times P$$

(1)

where SE is the annual average soil erosion quantity (t km² a⁻¹) (Figure 2 and Figure 3). R is the rainfall erosivity factor (MJ mm km⁻² h⁻¹ a⁻¹). K is the soil erodibility factor (h hm² km⁻² MJ⁻¹ mm⁻¹). LS is the slope length (L) and steepness (S) factor. C is the surface vegetation cover factor and P is a dimensionless conservation practice factor.

2.3.1. Rainfall Erosivity Factor (R)

Rainfall erosivity factor (R) is an important driving force of SE and it is an important factor for predicting SE [79,80]. This study uses the empirical formula derived by Wischmeier and Smith [81]. The formula, which fully considers the characteristics of rainfall, is one of the most scientific methods for calculating R and has been widely used in different regions.

$$R={{\displaystyle \sum }}_{i=12}^{12}1.375\times {10}^{\left[1.5\times {\log }_{10}\left(\frac{P_i{}^2}{P}\right)-0.08188\right]}$$

(2)

where, i represents the month in a year; P_i is the rainfall in month i (mm); P_a represents the annual rainfall (mm). R is the annual average rainfall erosivity (MJ mm/[ha·h·a]). We compared the R-factor results for Central Asia with the published R-factor datasets for the pan-third pole region (20 countries containing key regions, http://data.tpdc.ac.cn/zh-hans/ (accessed on 23 December 2021)), and the fitting results were good (R² = 0.853) (Figure S2), indicating that our results are accurate enough to be used in this study.

2.3.2. Soil Erodibility Factor (K)

Soil is the body on which SE occurs and soil erodibility significantly influences the occurrence probability of SE [82]. In this study, the extensively used EPIC model proposed by Williams et al. [83] was selected to calculate the soil erodibility factor (K value). The formula of the EPIC model is as follows:

$$\begin{array}{cc}{K}_{EPIC}=& \left\{0.2+0.3exp\left[-0.0256S_a\left(1-\frac{S_i}{100}\right)\right]\right\}\times {\left(\frac{S_i}{C_l+S_i}\right)}^{0.3}\\ & \times \left[1-\frac{0.25C}{C+exp\left(3.72-2.95C\right)}\right]\times \left[1-\frac{0.7S_n}{S_n+exp\left(-5.51+22.9S_n\right)}\right]\end{array}$$

(3)

$$S_n=1-\frac{S_a}{100}$$

(4)

where S_a is the sand content (2–0.05 mm), %. S_i is the silt content (0.05–0.002 mm), %. C_l is the clay content (<0.002), and %. C is the organic carbon content, %. For the K factor, we compared it with the published K factor data of the pan-polar third pole 20 countries (http://data.tpdc.ac.cn/zh-hans/ (accessed on 23 December 2021)), which was calculated based on the soil attribute data of the International Soil Reference and Information Centre. The two fitted well (R² = 0.768), indicating that the results of our K factor were more accurate and could be applied to Central Asia (Figure S3).

2.3.3. Slope Length and Steepness Factor (LS)

The LS factor was derived based on DEM (digital elevation model, 30-m http://hydro.iis.u-tokyo.ac.jp/~yamadai/MERIT_DEM (accessed on 25 October 2021)). The LS factor algorithm in this study was based on the studies of McCool et al. [84] as follows:

$$L={\left(\frac{\lambda }{22.1}\right)}^m$$

(5)

$$S=\{\begin{array}{c}10.8sin\left(\theta \right)+0.03 \theta \le 9{}^{°}\\ 16.8sin\left(\theta \right)-0.50 \theta \ge 10{}^{°}\end{array}$$

(6)

where λ is the horizontal projection slope length (m); m is the variable slope length index; and θ is the steepness of slope (°).

2.3.4. Cover Management (C)

Vegetation cover and management factor (C) reflects the influences of vegetation cover and crop management measures on SE quantity [85]. As vegetation cover is highly correlated with NDVI, C was calculated using the method proposed by Cai et al. [86]:

$$f=\frac{\left(NDVI-NDV{I}_{min}\right)}{\left(NDV{I}_{max}-NDV{I}_{min}\right)}$$

(7)

$$C=\{\begin{array}{c}1,f=0\\ 0.6508-0.3436\times {\log }_{10}f,0\le f\le 78.3\%\\ 0,f>78.3\%\end{array}$$

(8)

where f is the vegetation cover. NDVI_max and NDVI_min are maximum and minimum NDVI, respectively.

2.3.5. Support Practice Factor (P)

Support practice factor (P) reflects the inhibition effect of water and soil conservation measures on the occurrences of SE. It can be assigned according to the degree of water and soil conservation, which is reflected by different land use types [87,88]. P values range between 0 and 1. p = 1 indicates no conservation measure. Smaller p values indicate more effective water and soil conservation measures for inhibiting SE [89,90]. The P-factors in this study refer to the p-values of AL, GL, FL and PDL land use types in the relevant literature [32,45,88] to generate the final P-factor maps for the four major basins in CA. Finally, a diagram of P factor was generated (Figure 3).

2.4. Data Analysis

2.4.1. Trend Analysis

The unitary linear recursive analysis can simulate the variation trend of a single pixel. It uses the least square method to fit the slope of vegetation cover changes during a certain period [91]. The calculation formula is as follows:

$${\theta }_{slope}=\frac{n\times {{\displaystyle \sum }}_{i=1}^n\left(i\times A_i\right)-{{\displaystyle \sum }}_{i=1}^{n}i{{\displaystyle \sum }}_{i=1}^n{A}_i}{n\times {{\displaystyle \sum }}_{i=1}^n{i}^2-{\left({{\displaystyle \sum }}_{i=1}^n{A}_i\right)}^2}$$

(9)

where slope represents the slope of the variation trend; A_i is the numerical value of samples (SE, NDVI) in year i; and n is the research time. When θ_slope > 0 or <0 or =0, the sample value increases gradually or decreases continuously or remains constant, respectively, during this period.

2.4.2. Correlation Analysis

The correlation between vegetation cover and climatic factor was analyzed using the spatial analysis method based on the single pixel approach. Similar to Jiang et al. [92], the correlation coefficients of vegetation cover with temperature and precipitation were calculated using the following formula:

$$r_{xy}=\frac{{{\displaystyle \sum }}_{i=1}^n\left(x_i-\overline{x}\right)\left(y_i-\overline{y}\right)}{\sqrt{{{\displaystyle \sum }}_{i=1}^n{\left(x_i-\overline{x}\right)}^2{{\displaystyle \sum }}_{i=1}^n{\left(y_i-\overline{y}\right)}^2}}$$

(10)

where r_xy is the correlation coefficient between vegetation cover and climatic factors; x_i refers to the temperature or precipitation value in year i; x is the annual average temperature or precipitation; y_i is vegetation coverage in year i; y is the annual average value of vegetation cover; and i is the number of study years.

2.4.3. Landscape Pattern Analysis

The Shannon Diversity Index (SHDI) is a comprehensive index. It can reflect not only the number of vegetation types, but also the evenness of the spatial distribution of different vegetation types [93]. Based on the LUCC data in CA, the landscape diversity index from 1982 to 2017 was calculated by using the Fragstats 4.2 software (University of Massachusetts, Amherst, MA, USA). The calculation formula is as follows:

$$SHDI=-{\displaystyle \sum }_{i=1}^n{p}_i\ln p_i$$

(11)

where P_i is the proportion of type i in the whole landscape and n refers to the total number of plaque types in the landscape. SHDI represents the overall complexity degree of the landscape. SHDI = 0 indicates that the whole landscape is only composed of one plaque. With the increase in SHDI (SHDI ≥ 0), the plaque types increase or the distribution of different plaque types in the landscape becomes increasingly more uniform. The variations of SHDI in the study area from 1982 to 2017 are shown in Figure S1. In general, SHDI gradually increased in CA (+0.0003/year).

2.5. Model Validation

In the process of global changes, increases in atmospheric CO₂ concentration gradually led to global warming [94]. SE can cause loss of soil organic carbon, thus altering global carbon reserves and carbon cycles [95] and accelerating the global warming process [96]. Berhe et al. [97] discussed the role of SE in the biogeochemical cycles of essential elements such as carbon, nitrogen, and phosphorus and found that the area of land with SE increases with soil carbon density, and it was influenced significantly by the distribution of soil carbon density. Therefore, SE has significant correlation with soil carbon accumulation. On this basis, the soil carbon cumulative data in CA was downloaded from FAO (http://faostat.fao.org/site/630/default.aspx (accessed on 26 September 2021)) for fitting with the SE quantity of the major river basins in CA. The results showed good correlation (R² = 0.79, RMSE = 830.12, MAE = 830.11). To further verify the accuracy of the simulation results, annual sediment delivery data from 2006 to 2008 were collected from three hydrometric stations in TRB. The data of the water transmission is derived from the Annual Hydrological Report, China. The annual sediment delivery data is calculated by accumulating the daily sediment delivery data. A linear fitting analysis was carried out between the sediment delivery data in TRB and annual erosion quantity at hydrometric stations, and the results showed good degree of fitting (R² = 0.82, RMSE = 189.39, MAE = 186.64). Although we used two different data sources to validate the simulation results and obtained larger RMSE values (Figure 4), we found that the RMSE values in our validation results differed less from the MAE values. Therefore, our simulation accuracy is better for SE by water in the main watersheds of CA by means of parameter localization and other means.

3. Results

3.1. Spatial and Temporal Variation in Water Erosion by Water

As shown in Figure 5, there is a certain pattern in the spatial variation of the slope of SE in each basin. The grids with positive SE changes were mainly distributed near the rivers. According to the statistics of the area of pixels with slope of SE change greater than 0, the area of pixels with an increasing trend accounted for 59.3% of the total area of the four basins, indicating an overall increasing trend of SE. Compared with the other three basins, TRB has the smallest SE variation and fewer pixel points with high slope of SE variation are distributed. SDB shows the most prominent increasing trend, with pixels showing an increasing trend accounting for 64.2% of the total basin area. Among the four basins, SE was the most severe in SDB, which had the highest SE increasing trend, followed by ADB (60.2%), IRB (55.7%) and TRB (57.4%).

From 1982 to 2017, the average SE quantity of all four basins presented relatively consistent rising trend (Figure 6). Since 2000, the SE quantity in ADB, IRB, and SDB began to increase quickly and peaked in around 2003, but it dropped sharply in the following 3–5 years. The SE quantity in TRB did not show sudden changes in the end of 2000 but increased slowly. Moreover, the extremum did not occur within this period. According to the linear fitting analysis of SE changes in the four basins from 1982 to 2017, SE quantity changed the most rapidly in SDB and the change slope was 0.9, with the maximum reaching as high as 2728.68 t km⁻². The variation rate of SE quantity was the lowest (only 0.05) in TRB and the maximum erosion quantity was 279.58 t km⁻². The highest growth rate of SE quantity was observed in SDB (0.9), followed by ADB (0.88), IRB (0.33), and TRB (0.05).

3.2. Spatial and Temporal Variation of Vegetation

The spatial variation trend of NDVI was explored through unitary linear recursive analysis and the features of vegetation variation were simulated pixel by pixel. The statistics on the percentage of pixel area with NDVI growth in the total basin area were also obtained (Figure 7). According to the results, NDVI in most areas declined gradually during this period and the proportion of pixels with NDVI growth was extremely low. The proportion of areas with NDVI growth was the highest in ADB (15.5%), but it was lower than 10.0% in SDB, IRB, and TRB. As shown in Figure 7, the distribution of regions with NDVI growth followed some spatial laws. Areas with increasing vegetation in ADB and SDB were mainly around the Aral Sea and the Tianshan Mountains in the east. NDVI in IRB was mainly distributed near Balkhash Lake. Areas with increasing NDVI in TRB were mainly concentrated in oasis regions near the Tarim River. In general, pixels with decreasing NDVI in CA were far higher than those with increasing NDVI from 1982 to 2017, indicating that vegetation activities were generally weakening during this period.

Statistics on annual average values of NDVI in the four basins were obtained. Furthermore, the variation characteristics of vegetation over time were analyzed. As shown in Figure 8, the interannual variation of NDVI was wide in all basins and NDVI fluctuated within a broad range. Although the NDVI fluctuated widely, the trends of NDVI in ADB, IRB, and SDB were relatively consistent, rising from 1982 to 2000 and then declining after 2000. The variation characteristics of NDVI in TRB were different from those in the other three basins. Specifically, NDVI in TRB increased with continuous fluctuation from 1982 to 2017. To explore overall variation trends of NDVI in different basins, a linear fitting analysis was performed on the average NDVI in different years. NDVI was found to increase in TRB (+0.0002/year), but it generally declined continuously in the other three basins. Overall, only vegetation activity in TRB was strengthened from 1982 to 2017. NDVI in generally presented a declining trend, indicating that local vegetation growth was declining.

4. Discussion

4.1. Effects of Climatic Factors on Vegetation Changes

Under the background of current climatic changes, the growth condition and spatial distribution of vegetation around the world are changing significantly [98,99,100]. Vegetation also exhibits significantly different responses to different climatic changes [101,102,103]. In the response of vegetation growth to climate, temperature and precipitation are two key factors [104,105]. Climatic warming driven by global temperature rise might change the vegetation phenology and primary productivity of ecosystems, thus influencing ecosystem services [106,107]. Meanwhile, global warming may increase the growth period of vegetation in spring and autumn, thus extending the growing season [108] and increasing vegetation cover [109]. As the main factor controlling changes in vegetation growth, precipitation can influence vegetation growth directly by changing soil water content [110,111]. In CA, NDVI exhibited significantly different responses to climate. With respect to correlations of temperature and precipitation with NDVI, precipitation showed distinctly higher sensitivity to climatic changes than temperature. In this study (Figure 9), areas with positive correlation between NDVI and temperature accounted for 48.3% of total area. Areas with positive correlation between NDVI and precipitation reached as high as 80.4%, indicating that precipitation is the major influencing factor of vegetation conditions in CA. Yin et al. [112] investigated the correlation between spatial-temporal changes of vegetation and climatic changes in CA from 1982 to 2012. They concluded that NDVI had a weak negative correlation with annual temperature changes, but a positive correlation with precipitation, which is consistent with the results of this study.

In terms of the spatial distribution of correlations between climate and NDVI, vegetation activities were affected more prominently by precipitation, especially in plains in northwest regions of CA (Figure 9a). The influences of temperature on vegetation activities were mainly observed in eastern regions with relatively high altitudes (Figure 9b). This agrees with the research results of Zhang et al. [38] on the responses of vegetation changes to climatic changes in CA from 1982 to 2012. At the basin level, correlations between NDVI and precipitation in ADB, SDB, and IRB all conformed to above characteristics, but the correlation was opposite in TRB. In TRB, although temperature and precipitation showed distinct correlations with NDVI, the area with correlation coefficients between NDVI and temperature higher than zero (69.3%) was significantly larger than that of precipitation (51.9%). Specifically, regions where NDVI showed good responses to precipitation were mainly concentrated in mountainous areas of TRB and oases near the mainstream of the Tarim River. Investigating vegetation dynamics in TRB and their responses to hydrologic climatic factors, Wang et al. [113] concluded that vegetation changes generally presented significantly positive correlations with annual average temperature and local precipitation in most mountainous areas. This verified the spatial differences in the correlations of NDVI with temperature and precipitation in TRB.

4.2. Effects of Human Activities on Vegetation Changes

LUCC is not only an important content in studying global environmental changes, but also a clear reflection of the close relationship between human activities and the natural environment [114,115]. Analyzing the dynamic change process of land uses can reflect the influences of human activities on vegetation cover directly or indirectly [116,117]. The economics and policies of countries in CA have changed greatly since the collapse of the Soviet Union, thus resulting in significant changes in land use [118,119]. The Aral Sea problem has been the key concern of people [120,121]. During the existence of the Soviet Union, irrigation systems and water conservancy facilities were well established in the Aral Sea Basin [122]. Specifically, the Aral Sea Basin provided irrigation water to 6.60 × 10⁶ km² of farmland and cotton land in the 1960s, but the inflow volume and water level of Aral Sea dropped dramatically due to construction of water conservancy facilities such as reservoirs and canals [123,124]. After the collapse of the Soviet Union, edges of oases were reclaimed on a large scale over trans-boundary river basins such as ADB and SDB due to the lack of overall management at the country level as well as doubled populations [125,126]. Field irrigation intercepted inflow volume, thus continuously accelerating the shrinkage of the Aral Sea and worsening vegetation growth around it [127,128]. Regional climatic changes caused by global climatic changes are one of the causes of the Aral Sea crisis, but human activities have had a stronger influence on the Aral Sea [129,130]. In arid regions of China, large amounts of the water of the Tarim River have been used for irrigation due to disordered reclamation into farmlands [131]. Moreover, new reservoirs were constructed in the upper reaches. Consequently, the lower reaches of the Tarim River were in a cutoff state [132,133]. Taitma Lake is a terminal lake that has been associated with the large-scale death of Populus euphratica and shrubs along riverbanks and serious degradation of the ecological environment [134,135].

Human activities can change NDVI by altering the land use/cover mode [136,137]. An unreasonable mode of land use can destroy the soil structure, leading to SE and ultimately decreasing soil productivity and surface coverage [138,139,140]. In this study, LUCC data of CA in 1982 and 2017 were chosen for comparative analysis (Figure 10). Results showed that AL, GL, FL, DL, and WB underwent significant changes (>100 km²). Among them, the largest changes were observed in AL, GL, and FL, followed by WB and DL. The total area of GL, FL, and DL converting to AL amounted to 1195.74 km² and the total area of GL and DL converting to FL was 1298.46 km², which are far greater than that of areas converting to GL (867.32 km²) and DL (460.27 km²) (

Table 1. Land use transfer matrix from 1982 to 2017 (km²).

	2017
	Land-Use	AL	GL	FL	PDL	CL	DL	WB	Others
	AL	8399.36	65.85	62.93	0.10	60.32	10.25	3.08	0
	GL	268.03	17,410.85	598.23	1.51	21.87	103.58	5.83	0
	FL	682.83	442.85	12,043.65	1.24	10.04	48.57	7.80	0
1982	PDL	0.02	1.61	3.27	4206.14	0.31	5.63	0.70	0
	CL	0.03	0	0	0	38.6	0	0	0
	DL	244.88	424.47	700.26	4.90	5.14	22,256.18	21.32	0
	WB	6.12	30.04	19.91	25.43	0.13	356.69	1412.82	0
	Others	0	0	0	0	0	0	0	527.73

Note: AL represents arable land, GL represents grassland, FL represents forestland, PDL represents potential desert land, CL represents construction land, DL represents desert land, WB represents water body, and Others represents other land types.

). Generally, AL experienced the largest changes from 1982 to 2017, with an increase of 11.90%.

In the process of land use and land cover changes, AL changes are the primary influencing factor of regional NDVI [141,142]. An overlay analysis was performed between regions with large land use changes and the spatial variation trends of NDVI, and the results are shown in Figure 11. In the process of land use transfer, the regional area transforming into AL accounted for the highest proportion of NDVI, whereas the area of other land use types transforming into FL and GL accounted for significantly lower proportions. Li et al. [143] reported that the greening of vegetation in the arid region of northwest China from 1990 to 2010 due to anthropogenic activities such as tillage expansion has increased the regional NDVI by 26.7%. The expansion of AL area became the main driver influencing the increase in regional NDVI. This agrees with the research results of this study.

4.3. Feedback Relationship between Changing Environment and Water Erosion by Water

This study performed a comprehensive analysis on climate, vegetation, and human activity. The results (Figure 12b) showed that precipitation, temperature, Shannon diversity index, and vegetation cover influenced SE significantly. Among them, precipitation showed the most significant direct effect on SE (Figure 12a) and the path coefficient was 0.42. Moreover, the influences of precipitation on NDVI were found to be the most prominent with interactions among different factors, indicating that precipitation can influence SE indirectly by influencing vegetation growth. This might be because vegetation can relieve the spattering effect of raindrops on soils and the scouring effect of runoff on soils through the canopy and root systems, thus enabling the prevention of SE [144,145]. Zhou et al. [146] investigated SE in the lower reaches of the Minjiang River Basin and found that SE did not significantly increase in regions with high SE risk and vegetation coverage lower than 40% when vegetation recovered to 60%. If vegetation coverage was recovered to 78.3%, the C factor would increase by at least 0.001 and SE risk would decrease by more than 50%. On the whole, the correlation between precipitation and SE was mainly determined by land use and vegetation cover. In a similar study in America, Nearing et al. [147] found that given a fixed land cover, SE would increase by 25–55% in a year if precipitation increases by 10% during the erosion period. According to this study, temperature is the second primary influencing factor of SE after precipitation. Its influences on SE are slightly higher than those of SHDI and NDVI. The global temperature rise changed water circulation between the oceans and land, thus making SE more uncertain [148,149]. Ma et al. [95] carried out a meta-analysis on global SE and found that SE was intensified by 2.1% from 1982 to 2019 due to global warming. SE was the most intense in mid and high-altitude regions. Moreover, global warming can also cause changes of soil moisture content by influencing regional precipitation, which will directly influence the K factor in the RUSLE model [150,151]. In general, SE in CA is jointly influenced by climate, human activities, and vegetation. Among them, climatic factors are the primary influencing factors, followed by SHDI dominated by human activities.

The state of SE in the four basins was compared with that in the surrounding areas (

Table 2. Comparative analysis of SE with other surrounding areas.

Study Areas	SE Ranges (t/km2·a)	Periods	References
Houzhai River Basin	104–223	1973–2013	Li et al. [158]
Loess Plateau	2399–2957	2000–2008	Fu et al. [159]
Xinjiang Province	1350–2084	1985–2011	Zhang et al. [160]
Three-River Headwaters Region	38–54	2005–2015	He et al. [161]
European Union	2220–2460	2010	Panagos et al. [162]
Taohe River Basin	1099–1424	2000–2018	Wang and Zhao [163]
ADB	1490–1525	1982–2017	This study
SDB	2677–2729
IRB	1306–1321
TRB	273–297

). SE in ADB was found to be similar to that in Xinjiang. Xinjiang is a region with relatively serious drought, salinization, and desertification as well as water and soil loss [152]. SE in IRB was close to that in the Taohe River Basin. The Taohe River is the second tributary of the Yellow River and precipitation changes are the major cause of local SE changes in the corresponding basin [153]. IRB generally receives more precipitation than the other three basins in this study, which may explain the similarity in SE between IRB and the Taohe River Basin. The Loess Plateau represents the area with the most prominent SE problems in China [154], and SE in SDB was slightly lower than that in the Loess Plateau. This might be the collaborative consequences of unreasonable use of water resources and soil in SDB [155]. SE in TRB was very similar to that of the Houzhai River Basin in a karst region in South China, which was only 273–297 t hm⁻² a⁻¹. This may be attributable to the implementation of ecological water conveyance projects and artificial protection measures in TRB since 2000, which improved the ecological environment and thereby relieved land degradation in the basin [156,157]. In general, SE in the major basins of CA was relatively serious and the ecological environment in the basins require significant improvements.

4.4. Limitations and Prospects

Although the GIMMS NDVI3g database used in this study has relatively long time series, the spatial resolution is relatively low. Consequently, the quantitative study on vegetation cover changes based on remote sensing data might have some uncertainties [164,165]. This study mainly referred to parameters in previous studies for evaluating SE using the RUSLE model and it is an experimental formula, soil erosion amount between the simulation and the observation differs, which may introduce some errors in the measurement of SE in different basins [45,166]. Regarding the influences of human activities on SE, although the influences of LUCC-induced vegetation cover changes and spatial distribution characteristic changes on SE were analyzed, several other human activities, such as over grazing and agricultural management technology, that may influence vegetation cover were excluded [167,168]. Although we used higher resolution remote sensing data, there are large differences in the data products themselves (e.g., DEM), which can also cause some uncertainty in the assessment results. This will also be reflected in the factor calculation of the RUSLE model. For example, the calculation of K factor and R factor involves a variety of data parameters, and their results are correspondingly uncertain, while the calculation process of other factors is relatively simple, and fewer parameters are involved in the calculation, and their impact will be less than that of R factor and K factor [169,170]. In addition, the validation data used in this study are regional soil organic carbon and sediment delivery data of rivers. Although it can reflect the SE condition well, there are some errors. In the future study, we will conduct field tests of soil isotope tracing for this area to verify our simulation results more accurately.

Vegetation coverage is closely related to soil loss, and it has extremely important significance to the study of vegetation changes. The monitoring of vegetation changes is a long-term dynamic process [171]. Therefore, it is necessary to study vegetation cover changes in a longer time series. This study only discussed the relationship between vegetation cover in major basins of CA and regional climate factors preliminarily on a yearly scale, but annual vegetation changes as well as the relationship between interannual vegetation changes and climate in different seasons were not taken into account. Future studies should explore the responses of vegetation cover changes to temperature and precipitation at different spatial-temporal scales, aiming to further deepen our understanding on the relationship between vegetation and climate. As a preliminary study, this study reveals the different collaborative effects of climatic changes and LUCC on SE. However, the influencing mechanisms of climatic changes and LUCC on SE remain unclear. Therefore, the mechanisms also need to be further analyzed in depth at the earliest.

5. Conclusions

This study performed a quantitative evaluation of SE in major basins in CA from 1982 to 2017 using the RUSLE model and analyzed the major driving factors of SE changes in different basins under changing environments. Spatially, areas with increasing SE accounted for 58.9% of the total area of the four basins and pixels with high variation rates of SE were mainly concentrated in regions with high relief. Temporally, the SE quantity of all four basins presented relatively consistent rising trend with fluctuations. Specifically, the largest increase in SE quantity was observed in SDB (64.2%), followed by ADB (59.4%), IRB (56.5%), and TRB (55.7%). Regarding vegetation changes, pixels with positive spatial variation trends of NDVI were mainly around tail-end-lakes and rivers. In all four basins, areas with decreasing NDVI accounted for more than 90.0%. Except for TRB, vegetation cover exhibited a decreasing trend in all three basins. In terms of the correlation between climate and vegetation, the correlation between NDVI and precipitation in CA (80.5%) was significantly higher than that of temperature (48.3%) and precipitation was the major climatic factor affecting vegetation growth. However, temperature had a stronger influence on vegetation (69.3%) than precipitation (51.9%) in TRB. Among land use types in CA, AL underwent the largest change from 1982 to 2017, increasing by 999.38 km². In the process of human activities influencing vegetation in CA, the NDVI proportion of other land use types converting to AL was the largest.

According to a comprehensive analysis of climate, vegetation, and human activities, SE in major basins in CA is jointly influenced by climate, human activities, and vegetation. Among them, climatic factors exert the strongest influence, followed by SHDI (human activities). Comparing the results of this study and other studies on surrounding areas, we found that SE is more serious in ADB, SDB, and IRB, with SE quantity in SDB even nearing that in the Loess Plateau, China. Under the background of global changes, appropriate water and soil resource management and optimization configuration should be implemented in CA with references to that in TRB. In this manner, SE problems in CA may be overcome.