The Assessment of the Spatiotemporal Characteristics of Net Water Erosion and Its Driving Factors in the Yellow River Basin

Abstract

The Yellow River Basin (YRB) is an important grain production base, and exploring the spatiotemporal heterogeneity and driving factors of soil erosion in the YRB is of great significance to the ecological environment and sustainable agricultural development. In this study, we employed the Revised Universal Soil Loss Equation (RUSLE) in conjunction with Transport-Limited Sediment Delivery (TLSD) to explore a modified RUSLE-TLSD for use assessing net water erosion. This modification was performed using sediment data, and the explanatory power of driving factors was assessed utilizing an optimal parameters-based geographical detector (OPGD). The results demonstrated that the modified RUSLE-TLSD can accurately simulate the spatiotemporal distribution of net water erosion (NSE = 0.5766; R² = 0.6708). From 2000 to 2020, the net water erosion modulus in the YRB ranged between 1.62 and 5.33 t/(ha·a). Specifically, the net water erosion modulus decreased in the YRB and the middle reaches of the YRB (MYRB), but it increased in the upper reaches of the YRB (UYRB). The erosion occurred mainly in the Loess Plateau region, while the deposition occurred mainly in the Hetao Plain and Guanzhong Plain. The Normalized Difference Vegetation Index (NDVI) and slope emerged as significant driving factors, and their interaction explained 31.36% of YRB net water erosion. In addition, the redistribution of precipitation by vegetation and the slope weakened the impact of precipitation on the spatial pattern of net water erosion. This study provides a reference, offering insights to aid in the development of soil erosion control strategies within the YRB.

1. Introduction

Rapid climate changes and frequent human activities have led to significant changes in soil water erosion and sediment delivery [1,2]. These changes have had further implications for agricultural productivity, food security, and soil fertility. As an important agricultural production area in China, the YRB has also attracted increasing attention from scholars due to the changes in water erosion intensity and spatiotemporal patterns. Some scholars [3,4,5] have analyzed the potential water erosion intensity (i.e., water erosion modulus) based on the Revised Universal Soil Loss Equation (RUSLE) [6]. The RUSLE model integrates a multitude of factors, such as rainfall, topography, vegetation, soil, etc., and has been widely used worldwide because of the limited difficulty it exhibits in acquiring data and its ease of operation.

During the research process, relevant scholars [5,7,8] have revised the RUSLE factor calculation formula to improve the reliability of the simulation results regarding potential water erosion based on prior experience and slope runoff tests. However, the RUSLE model faces challenges in assessing simulation accuracy when combined with measurement data, and it cannot simulate the sediment transport and deposition. As a result, recent studies have increasingly adopted net water erosion (i.e., sediment yield) analysis models such as SWAT [9], EPIC [10], WaTEM/SEDEM [11], the RUSLE-coupled index of connectivity (RUSLE-IC) [12], RUSLE-coupled Transport-Limited Sediment Delivery (RUSLE-TLSD) [13], etc. These models can present the whole simulation process of soil separation–transportation–deposition seen in water erosion. However, it is noteworthy that current research predominantly focuses on sub-regions within the YRB. For example, Liu and Fu (2016) [14], Yan et al. (2022) [15], and Zhao et al. (2020) [16] carried out the analysis of net water erosion intensity in the Yanhe River Basin based on one of the aforementioned models. Based on the WaTEM/SEDEM, Pal et al. (2018) [17] and Zhang et al. (2023) [18] analyzed the characteristics of net water erosion in the Shejiagou Catchment and the Taohe River Basin separately. Furthermore, utilizing the RUSLE-TLSD model, Li et al. (2022) [19] uncovered the spatiotemporal patterns of potential erosion and net erosion in the Upper Yellow River Basin (UYRB).

Previous research by other scholars successfully calibrated the model using hydrological station data, achieving reliable accuracy. However, given the vast extent of the YRB, obtaining detailed sediment measurement data proved challenging. Additionally, the effectiveness of erosion control measures exhibited significant scale effects [1]. Consequently, the calibration model derived from sub-regions cannot be directly applied to the entire basin, presenting a challenge for model modification and verification in large-scale basins. To address this challenge, our study drew on the research conducted by Li et al. (2022) [19]. The RUSLE-TLSD model was a semi-distributed and semi-empirical soil water erosion and sediment delivery model. This model can be calibrated in combination with hydrological station data to describe the spatial pattern of net water erosion and sediment transport. In comparison to complex hydrological models such as SWAT, the RUSLE-TLSD model offers simpler and more efficient operation. Although it is currently applied in only a few small watersheds in India [13,20,21,22], the Hexi Corridor in China [23], and the UYRB [19], the model’s principles are mature and rigorous. The RUSLE-TLSD model considers that sediments are detached from the soil surface and transported downhill along the water flow, with the ability of rivers to transport sediment being variable and limited [13]. At some point along the downslope, the decrease in slope leads to a weakening of streamflow transport capacity and causes sediment deposition [24].

In addition to assessing water erosion and sediment transport, it is also important to investigate the spatiotemporal patterns of driving factors and propose control measures. Currently, the widely utilized methods for analyzing soil erosion-driving factors encompass correlation analysis, the optimal parameters-based geographical detector (OPGD), geographically weighted regression (GWR), double-mass curves, attribution theory, etc. Among these methods, the OPGD model stands out for its ability to unveil rich geographic features and information. It excels in accurately and efficiently analyzing spatial heterogeneity and facilitates the assessment of various factors’ impact on water erosion [25,26]. And it is widely used in work on the driving factors of soil erosion because of its non-linear assumptions and clear physical meaning.

The Yellow River previously held the record for the highest sediment load [27], with over 130 Gt of soil eroded between 1950 and 2010 [28]. Over the past few decades, the sediment load in the YRB has changed dramatically. Sediment deposition in the main channel of the UYRB has intensified, while the sediment transport ratio in the midstream area has decreased significantly. In the 1970s, the Tongguan Hydrological Station transported approximately 1.6 billion tons of sediment annually [29]. However, the sediment load decreased to just 240 million tons in 2020. Before the year 2000, engineering measures such as check dams were the main reason for weakening water erosion in the YRB [30]. After 2000, the implementation of various comprehensive management measures in small watersheds became widespread. Amid the influence of complex climate changes and human various activities, the spatiotemporal patterns and driving mechanisms of water erosion and sediment transport in the YRB have undergone profound transformations [31,32].

Therefore, this study focused on the YRB to explore the trend of water erosion changes from 2000 to 2020 while analyzing the spatiotemporal characteristics of driving factors. We attempted to (1) modify the RUSLE-TLSD model based on the sediment measurement data; (2) simulate and analyze the spatiotemporal pattern of net water erosion and sediment deposition in the YRB; and (3) utilize the OPGD model to unveil the driving factor variations across different spatiotemporal scales, with year-scale, month-scale, large-basin-scale, and sub-region-scale analyses. The research objective is to provide a theoretical reference for soil and water resource conservation in the YRB by conducting soil water erosion and driver factor assessments, and to provide a model reference for the modeling and analysis of soil water erosion in other large basins.

2. Materials and Methods

2.1. Study Area and Flowchart

The YRB covers an area of approximately 797,000 km², comprising an internal drainage area (IDA) spanning 45,000 km²; the UYRB, occupying 371,000 km²; the midstream area (MYRB), totaling 356,000 km²; and the downstream area (DYRB), occupying 25,000 km² (Figure 1). The topography within the basin varies significantly, resulting in notable differences in terrain. The Yellow River flows through diverse regions. Affected by atmospheric circulation and monsoon circulation, the climate in the YRB is complex, including a plateau climate, an arid climate, a semi-humid climate, etc. Moreover, the YRB is the main area of grain production, and the YRB’s grain output accounts for about 35% of the total grain output in China. The Yellow River is well known for its high sediment content, which has the spatial characteristics of “different sources of water and sediment”. The UYRB is the major source region of the Yellow River’s water volume, with sediment concentrations primarily resulting from torrential rain and floods in the Loess Plateau area (http://www.yrcc.gov.cn/gzfw/nsgb/, accessed on 21 May 2022). Under natural conditions, in the past, the UYRB featured an alluvial channel with a gradual rise in the riverbed. Conversely, the MYRB had no obvious sedimentation, and the sediment delivery ratio was close to 1. Sediment deposition predominantly occurred within the DYRB’s channel.

The overall flow chart of this study is shown in Figure 2. In order to further analyze the changes in net water erosion characteristics and driving factors in the YRB, the study firstly explored a modified RUSLE-TLSD model for water erosion simulation in the YRB by using meteorological data, topographic data, vegetation indices, soil data, and observed sediment data from hydrological stations. Combined with the sediment data, the accuracy of the model before and after the modification was compared using Nash and Sutcliffe efficiency (NSE). This study analyzed the spatial and temporal characteristics of water erosion at the basin scale and sub-basin scale, and further analyzed the water erosion and deposition area in the upstream and middle reaches by combining the terrain niche index. After that, the driving factors of soil water erosion in the YRB were quantitatively analyzed based on the OPGD model.

2.2. Data

The sources and purpose of the datasets used in this study are shown in

Table 1. The dataset used in the research and its source, purpose.

Dataset	Resolution	Source	Purpose
Basic boundary dataset of soil and water conservation in the Yellow River Basin	-	National Cryosphere Desert Data Center (http://www.ncdc.ac.cn, accessed on 26 August 2022)	Define the study’s scope.
Meteorological data	-	China Meteorological Information Center (https://data.cma.cn/, accessed on 10 August 2021)	Calculating R Factor, etc.
Hydrological stations sediment observation data	-	Ministry of Water Resources of the People’s Republic of China (http://www.mwr.gov.cn/sj/#tjgb, accessed on 30 October 2022)	Model calibration and validation.
ASTER GDEM V3.0	30 m	Geospatial Data Cloud (https://www.gscloud.cn/, accessed on 18 January 2022)	Calculation of the LS factor, As, etc.
MOD13Q1 V006	250 m	NASA Earth Science Data Systems (https://www.earthdata.nasa.gov/, accessed on 6 October 2022)	Calculation of the C factor, etc.
Grid Data on Soil Erodibility in China	30 m	Loess Plateau Data Center, National Earth System Science Data Center (http://loess.geodata.cn, accessed on 12 December 2022)	Extracting the K factor.
Soil Map-Based Harmonized World Soil Database (v1.2)	1 km	National Tibetan Plateau/Third Pole Environment Data Center (https://data.tpdc.ac.cn/, accessed on 17 April 2022)	Exploring the effect of soil type on erosion
China land cover dataset	30 m	Zenode (https://doi.org/10.5281/zenodo.4417810, accessed on 9 April 2023)	Calculation of P factor, etc.
Soil thickness map	90 m	Soil Sub Center, National Earth System Science Data Center (http://soil.geodata.cn, accessed on 12 May 2023)	Calculation of T factor, etc.

. Among them, the LUCC data from the China Land Cover Dataset, which is a Landsat-derived product created by Yang and Huang (2021) [33] based on the Google Earth Engine (GEE), have an accuracy of 79.31%. The soil thickness data from Liu et al. (2022) [34] are listed. These data were generated based on machine learning and approximately 5000 representative soil profiles. This was more detailed and accurate than previous research methods. Due to limitations imposed by NDVI and personal computer computing power, all data were resampled to 250 m. Detailed descriptions of data applications can be found in the relevant chapters.

2.3. Modified RUSLE-TLSD Model

2.3.1. Modified RUSLE Model

In the previous study [35], preliminary analysis was carried out on the temporal–spatial features of potential YRB soil water erosion. The RUSLE is calculated as follows:

A = R × K × LS × C × P,

(1)

where A (t·ha⁻¹·a⁻¹) represents the potential water erosion rate; R (MJ·mm·ha⁻¹·h⁻¹·a⁻¹) is the rainfall erosivity factor; K (t·h·MJ⁻¹·mm⁻¹) is the soil erodibility factor; LS is the slope factor and the slope length factor; C is the vegetation cover and management factor; and P is the conservation practices factor.

Based on the formula of Zhang et al. (2022) [36], this study calculated the R factors of 163 meteorological stations in the YRB and its vicinity using daily precipitation data and then interpolated them using ANUSPLIN 4.37 software to obtain the R factor dataset. The K factor raster dataset was extracted by mask from the Grid Data on Soil Erodibility in China. The LS factor was calculated by applying the LS Tool 2.2.1 from the work of Fu et al. (2015) [37]. The calculation equation of K factor and LS factor is given in the Supplementary Materials. In addition, based on the method proposed by Cai et al. (2000) [38], this study estimated the vegetation coverage of the YRB using the MOD13Q1 data and then, based on this, evaluated the C factors. Referring to related research [3,5,16,19,23,39,40], the detailed values of the P factor are shown in

Table 2. P factor values in the YRB.

LUCC Type	Cropland	Forest	Shrub	Grassland	Water	Sonw/Ice	Barren	Impervious	Wetland
p Value	0.2 + 0.03 × s	1	1	1	0	0	0	0	1

, where s represents the slope (%).

To modify RUSLE model, this research added the M factor (dimensionless). This factor is used as a modified layer for model calibration on a large basin scale. The modified model is as follows:

A = R × K × LS × C × P × M,

(2)

As a direct object of water erosion, the soil’s characteristics will have a degree of influence on the water erosion intensity. Soil thickness, as an important indicator of soil, reflects one of the physical properties of soil, plays a role in controlling the hydrological processes of hillslopes, and is able to determine hydrological and erosive behavior to a certain extent [41]. Under the same erosion conditions, a thicker soil layer provides more material, implying that thicker soil may suffer more erosion and produce more sediment. Therefore, this study explores the M factor using the soil thickness dataset (Figure S1), and the formula is as follows:

$$M=a\times ST+b$$

(3)

where ST is the soil thickness; a and b are calibration coefficients, and their values are corrected by the sediment data; and M ≥ 0.

2.3.2. Transport-Limited Sediment Delivery

The net water erosion in the catchment area is limited by the supply of material and the transport capacity; that is, the net water erosion per grid depends on the potential soil erosion that may occur within the grid and the upstream sediment inflow [42]. As shown in Equations (4) and (5) [43], if the sediment transport capacity (TC, t ha⁻¹ a⁻¹) is greater than the sum of potential water erosion and upstream sediment influx, then the summation will be the sediment flux generated in the grid. That is, the amount of sediment engendered by the grid unit itself and the sediment imported from the upstream will be transferred to the downstream grid. Otherwise, the sediment flux value of the grid is equivalent to TC, and sediment deposition occurs in the grid.

$$T_{out}=\min \left(A+\sum T_{in},TC\right),$$

(4)

$$D=A+\sum T_{in}-T_{out},$$

(5)

where $T_{out}$(t·ha⁻¹·a⁻¹) is the unit sediment flux of outflow; A is the potential soil water erosion; $T_{in}$ (t·ha⁻¹·a⁻¹) is the amount of sediment flowing into the grid upstream; and D (t·ha⁻¹·a⁻¹) is the modulus of sediment deposition. TC is affected by vegetation cover, rainfall runoff, terrain, etc. Equation (6) is proposed by Verstraeten et al. (2007) [44].

$$TC=K_{TC}\times R\times K\times A_s^{1.4}\times S^{1.4},$$

(6)

where $K_{TC}$ is the correction parameter of TC (m); $A_s$ is the upslope confluence area per width of contour (m²/m); and S is the slope gradient (m/m). $K_{TC}$reflects the vegetation composition within the TC [13,44], and its calculation formula is shown in Equation (7).

$$K_{TC}=\beta \times \exp \left({\displaystyle \frac{-NDVI}{1-NDVI}}\right),$$

(7)

where $\beta$ is the model correction coefficient, which is used to reduce the error between the simulated value of the modified RUSLE-TLSD model and the measured value of the hydrological station.

In the study, the flow direction was assigned using the D-Infinity Flow Directions [42], and net water erosion (t·ha⁻¹·a⁻¹) was acquired by subtracting the D from the A. And the net water erosion intensity was divided into seven levels according to the erosion modulus (t·ha⁻¹·a⁻¹): deposition (<0), tolerable (0–10), slight (10–25), moderate (25–50), severe (50–80), very severe (80–150), and destructive (>150). The batch processing code of the modified RUSLE-TLSD model was written based on the TauDEM (a suite of DEM tools for the extraction and analysis of hydrologic information from topography, as represented by a DEM, found at https://hydrology.usu.edu/taudem/taudem5/, accessed on 1 September 2022) and the ArcPy (https://pro.arcgis.com/en/pro-app/latest/arcpy/get-started/what-is-arcpy-.htm, accessed on 12 September 2022).

2.4. Terrain Niche Index

TNI is a comprehensive index reflecting information on elevation and slope, and its level represents the degree of spatial variation in topographic conditions, or the complexity of the terrain in the study area. The higher the elevation and the greater the slope, the higher the TNI level of the area. Combining slope and elevation together to study the effect of changes in topographic gradient on water erosion intensity can accurately reflect the distribution trend of soil erosion intensity in relation to topographic changes [45,46]. This is evaluated using Equation (8):

$$TNI=\log \left[\left({\displaystyle \frac{E}{\overline{E}}}+1\right)\times \left({\displaystyle \frac{SL}{\overline{SL}}}+1\right)\right],$$

(8)

where $E$ and $\overline{E}$ represent the altitude value of any place and the average altitude of the region where the place is located; and $SL$ and $\overline{SL}$ represent the slope of any place and the slope mean of the area where the place is located, respectively. To further analyze the spatiotemporal pattern of water erosion on a different TNI, the natural break method was used to divide the TNI into eight levels (as shown in Figure S2).

2.5. Optimal Parameters-Based Geographical Detector

The Geodetector model is the foundation of the OPGD model used to analyze spatial variability and reveal influence ability [25]. In this study, Net water erosion was used as the Y variable, and LUCC, soil type, landform type, TNI, slope, elevation, precipitation, NDVI, and soil thickness were used as the X variables. The q-value in OPGD model is an important statistic in assessing spatial dissimilarity as well as the importance of factors explaining spatial variability. A higher q-value indicates that the factor explains the spatial variability of water erosion more strongly. The q value ∈ [0–1], but the grading method and number of intervals of continuous driving factors affect the q value. Therefore, this study used the “GD” package (https://cran.r-project.org/web/packages/GD/vignettes/GD.html, accessed on 16 November 2021) to achieve continuous factor optimal discretization and geographic detection batch processing.

2.6. Model Accuracy Evaluation

In this research, the RUSLE-TLSD model was modified and verified based on the sediment data of hydrological stations (Figure 1), and the Nash and Sutcliffe efficiency (NSE) [47] was used to evaluate the accuracy.

$$NSE=1-{\displaystyle \frac{\sum _{i=1}^n{\left(O_i-{MS}_i\right)}^2}{\sum _{i=1}^n{\left(O_i-\overline{O}\right)}^2}},$$

(9)

where n is the count of observations, $O_i$ is the measured value, ${MS}_i$ is the model simulated value, and $\overline{O}$is the average value of the observed value. NSE ∈ [−∞, 1].

3. Results

3.1. Model Modification and Verification

During this study, a total of 358 sediment observation data from the China River Sediment Bulletin (http://www.mwr.gov.cn/sj/#tjgb, accessed on 30 October 2022) were used to modify and verify the RUSLE-TLSD model. Among them, 192 data were employed for model calibration and came from Lanzhou Station, Shizuishan Station, Tongguan Station, Sanmenxia Station, Hongqi Station, Baijiachuan Station, Xianyang Station, Zhuangtou Station, Wuzhi Station, and Hejin Station. The remaining 166 data were used for model validation and came from Tangnaihai Station, Toudaoguai Station, Longmen Station, Xiaolangdi Station, Gangu Station, Zhangjiashan Station, and Huaxian Station. In Figure 3, we present a comparison of the NSE changes at different β values before and after modifying the RUSLE-TLSD model and the linear regression results when β = 1 at the calibration station. Before the modification, when β = 1, the NSE was only 0.0078 (Figure 3a). The linear regression results were also unsatisfactory (R² = 0.0713, Figure 3b). Upon modification, the reliability of the modified RUSLE-TLSD model significantly improved when applied at larger scales. At β = 1, with T factor parameters a = 0.0085 and b = -0.4, the calibration station’s NSE increased to 0.5766 (Figure 3a), and the linear regression effect was relatively better (R² = 0.6708, Figure 3b). Under these conditions, the NSE at the verification station was 0.5369. Therefore, when with β set to 1, a set to 0.0085, and b set to -0.4, the modified RUSLE-TLSD model is suitable for simulating and analyzing net water erosion in the YRB for the period from 2000 to 2020.

To further verify the reliability of the modified models, we compared the results of the current modified RUSLE model with those used in other studies (Table S1). The results showed a strong correlation (R² = 0.87) between the current research results and other related research results. At the same time, we combined the online historical images provided by ArcGIS Pro and selected areas A, B, and C (Figure 4) to compare the simulation details of net water erosion patterns in 2020 with the remote sensing images from July 2020 and analyze them. The modified RUSLE-TLSD model effectively simulated the spatial distribution details of water erosion and sediment deposition in the YRB to a certain extent. The sediment deposition in the three regions mainly occurred in the ravines, and the erosion pattern was closely related to the distribution of ravines. Although the details of the model simulation results are limited by data resolution, sediment deposition in area A shows a northwest–southeast trend similar to that seen in gullies. While vegetation coverage in area B is better than in area C, the net water erosion intensity of B is significantly lower than that in area C. The net water erosion patterns in areas B and C are highly consistent with vegetation coverage, river network systems, etc., indicating the reliability of the model.

3.2. Model Factor Characteristics

The soil erosion process in the YRB is full of uncertainty and complexity, being under the joint influence of various factors such as rainfall, vegetation, topography, and human activities. As shown in the Figure 5, the average value of the multi-year R factor in the YRB is 1353 (MJ·mm)/(ha·h·a), with a spatial trend of increasing from northwest to southeast. The K factor is in the range of 0–0.0258 (t·h)/(MJ·mm), with an average value of 0.0108 (t·h)/(MJ·mm), and the areas with high values are mainly concentrated in the Loess Plateau area. The average value of the LS factor is 12.662, and the maximum value is 243.74, which fully reflects the role of topography in determining the LS factor in space, and generally shows the distribution characteristics of high LS factor values in highland and hilly areas and low LS factor values in plain areas. The average value of the multi-year C factor is 0.1309, and the areas with higher C factor values are mainly distributed in the Hetao plain sides. The mean value of the multi-year P factor is 0.7883, and the areas with low P factor values are mainly located in Ningxia Yellow Diversion Irrigation Area, Hetao Plain, Guanzhong Plain, etc., which are dominated by cultivated land and built-up areas. The mean value of M factor is 0.75, and spatially the areas with low values of M factor (0–1.25) are located mainly in the western, northern, and eastern parts of the YRB, while those with high-value (1.25–2.49) areas are concentrated in the Loess Plateau area in the central part of the basin.

3.3. Spatiotemporal Distribution of Net Water Erosion

3.3.1. Annual Changes of Net Water Erosion Spatiotemporal Pattern

During the research period, the multi-year average net water erosion rate in the YRB was 3.12 t/(ha·a), with a multi-year average sediment load of 248 million tons. Within the YRB, the multi-year average net water erosion modulus in the UYRB was about 1.66 t/(ha·a), and the multi-year average sediment load was 61 million tons, constituting approximately 24.60% of the total YRB sediment load. The multi-year average net water erosion modulus in the MYRB was about 5.27 t/(ha·a), with a multi-year average sediment load of about 187 million tons, representing about 75.40% of the YRB sediment load (Table S2). Spatially, net water erosion was relatively concentrated in the Loess Plateau. The sediment from water erosion followed the river courses through the valley or plain, with the river’s sediment-carrying capacity diminishing as the slope decreases. Therefore, sediment produced from UYRB erosion was accumulated to a relatively significant degree in the Hetao Plain, while sediment from the MYRB settled in significant quantities in the Guanzhong Plain and the North China Plain. From 2000 to 2020, the overall change in net erosion in the Yellow River Basin was not significant, the regions with the significant enhancement of net erosion intensity were scattered, and the regions with significant weakening of net erosion were mainly located in the Loess Plateau region (Figure 6).

The YRB net water erosion modulus from 2000 to 2020 ranged from 1.62 to 5.33 t/(ha·a), indicating a trend of decreasing fluctuations (slope: −0.063, Figure 7). The net water erosion in the MYRB was relatively strong, and erosion rate ranged from 2.46 to 10.61 t/(ha·a), but it showed decreasing fluctuations (slope: −0.1688, Figure 7), similar to the entire basin, and the risk of soil erosion faced by the YRB and MYRB will gradually decrease in the future. Although the UYRB net water erosion was relatively weak, with an erosion modulus between 1.09 and 3.61 t/(ha·a), it generally showed a slow increasing trend (slope: 0.0269, Figure 7). Regarding this trend, we need to pay more attention to the risk of the UYRB soil erosion in the future.

Working with different TNI grades, the net water erosion exhibited partial similarities but also demonstrated partial heterogeneity (Figure 8). Within TNI levels 1–6, both the MYRB and UYRB displayed a pattern of initially decreasing and then increasing in terms of net water erosion modulus. The largest sediment deposition rate occurred at TNI level 3, while the most intense erosion was observed at TNI level 6. However, where the TNI level exceeded 6, different net water erosion change trends emerged between the UYRB and MYRB. In the UYRB, net water erosion exhibited insignificant changes, whereas in the MYRB, the net water erosion modulus experienced a significant decrease. Specifically, the erosion was relatively intense in the TNI 5–8 level area of the UYRB (UYRB_TNI) and the TNI 5–6 level area of the MYRB (MYRB_TNI), and the multi-year average net water erosion modulus values were 4.52 t/(ha·a) and 21.21 t/(ha·a), respectively.

3.3.2. Monthly Changes in Net Water Erosion

During the winter months (December to February), most areas within the YRB have less precipitation. Consequently, these months were excluded from our study. In general, monthly net water erosion, influenced by variations in rainfall, exhibited a temporal pattern characterized by initial intensification followed by attenuation over time (Figure 9). The spatial distribution changes in water erosion remained relatively stable during this period. But the gravity center (Figure S3) shifted roughly in an east–west direction. From October to March, net water erosion was generated sporadically in the MYRB and UYRB, and the water erosion intensity was relatively weak, standing at less than the YRB multi-monthly average of 0.7799 t/(ha·a) (Figure 9). In spring (March to May) and autumn (September to November), net water erosion displayed distinct trends of gradual intensification and subsequent weakening, respectively. Notably, the gravity center of net water erosion in the spring shifted westward (Figure S3). In the autumn, the distribution range decreased for net erosion, with the gravity center of erosion moving eastward (Figure S3). In summer (June to August), the hazard of net water erosion was relatively high, especially in the MYRB, where substantial changes were observed (Figure 9). In the MYRB, the net water erosion rates from June to August were 1.8475 t/(ha·a), 2.6082 t/(ha·a), and 2.0667 t/(ha·a), respectively, markedly surpassing those seen in the YRB and UYRB. The severely eroded areas manifested as northeast–southwest strips, and were mainly concentrated on the Loess Plateau. Unlike the MYRB, where the net water erosion was strongest in July, the UYRB saw the strongest net water erosion in August, which may be related to the delayed rainfall belt and the increase in ice and snow melting caused by rising temperatures.

3.4. Driving Forces of Net Water Erosion

At the basin scale, the explanatory ability of each driving factor was not high (Figure 10). This was closely related to the complexity of soil water erosion processes, which are to some extent caused in large basins by multiple interrelated factors [48]. In addition, factor detection revealed that the q-value of driving forces in the YRB displayed an overall decreasing trend, along with fluctuations. This trend signified the increasingly intricate interplay of factors influencing net water erosion in the YRB, with individual factors exerting a more limited influence. In addition, compared with Figure 7, it can be seen that in years with strong erosion, the q-value increased significantly, while in years with weak erosion, the q-value decreased significantly.

The multi-year average q-value indicated that, in the YRB, the slope had the strongest explanatory capability (10.48%), followed by soil thickness (10.01%), TNI (8.11%), NDVI (6.88%), soil type (5.36%), elevation (4.35%), rainfall (4.02%), LUCC (2.85%), and landform type (2.77%) (Figure 11). At the sub-region scale, significant variations in the explanatory capacity of driver factors were observed between the UYRB and the MYRB. In the UYRB, the driving factors with strong explanatory power were soil thickness (7.71%), the slope (7.57%), and elevation (7.51%). Except for elevation, which displayed a higher q-value than that in the YRB, the explanatory capabilities of other factors were lower than those in the YRB, implying that there was a more complex driving mechanism behind net water erosion in the UYRB. In the MYRB, the explanatory capability of NDVI and LUCC increased significantly, being capable of explaining 13.62% and 11.45% of phenomena, respectively. This suggested that human-related factors such as LUCC and vegetation cover had a more pronounced influence in the MYRB. In addition, this study conducted an analysis of the driving forces in areas with strong net water erosion. In the UYRB_TNI, topographic factors (slope and TNI) exerted a greater impact on soil water erosion, with the q-value of the slope exceeding 15%. In the MYRB_TNI region, NDVI explained 34.89% of net water erosion, while LUCC explained 16.22% of net water erosion. And the explanatory power values of the slope and TNI decreased to 6.18% and 5.51%, respectively.

The interaction detection evinced that the strongest explanatory power was observed when NDVI interacted with slope. And in the YRB, UYRB, MYRB, UYRB_TNI, and MYRB_TNI, the explanatory power values for this interaction were 31.36%, 27.09%, 47.70%, 22.71%, and 52.83%, respectively (Figure 12). Across regions of different scales, elevation and slope exhibited strong explanatory capabilities when interacting with other factors. In the MYRB and the MYRB_TNI regions, NDVI and LUCC also demonstrated high explanatory power when interacting with other variables. These findings further proved the factor detection results, indicating that in the UYRB, net water erosion was primarily influenced by topographical factors, whereas in the MYRB, it was influenced jointly by topographical and human-related factors. Due to the complex topography, although elevation and slope had high explanatory power, the explanatory capability of TNI when interacting with other variables was relatively weak.

The explanatory power of influence factors varied significantly across different months, with topographical factors (slope, TNI) emerging as the primary controlling factors overall. During the concentrated rainfall season (June–August), the explanatory capacity of multiple driving factors demonstrated significant augmentation (Figure 13). In the MYRB and MYRB_TNI regions especially, there was a significant increase in the explanatory power of precipitation and NDVI in July. In the MYRB_TNI region, the explanatory power of NDVI and rainfall reached 0.2233 and 0.1257, respectively, in July. But in the UYRB, the explanatory power of precipitation displayed a trend of decreasing first and then increasing from June to August, which could likely be attributed to glacial snowmelt water and the time-lag effect of vegetation on precipitation changes [49]. The spatiotemporal heterogeneity of overall and local changes in erosion-driving factors presented some challenges in terms of understanding and analyzing water erosion mechanisms at large scales.

4. Discussion

4.1. Model Availability and Error Sources Analysis

It became evident that the unmodified RUSLE model significantly underestimated the erosion rates in several tributaries as compared to the hydrological station data (Table S3). This discrepancy was, in part, attributable to the RUSLE factor calculation formula employed in this study. Specifically, the calculation of R values utilized daily rainfall data. Recent research by Yue et al. (2022) [50] has demonstrated that results based on daily rainfall data [51] tend to overestimate erosion compared to those derived from high-density meteorological stations, providing hourly rainfall data. Furthermore, in the calculation of the C factor, our study did not account for the influence of seasonal changes in vegetation growth [52] and relied on MODIS products. Nevertheless, Yang et al. (2020) [53] discovered that the C factor calculated using MODIS was lower than that determined using Landsat products, particularly in the Loess Plateau. Setting different cumulative percentages of NDVI when estimating vegetation coverage also impacted the results of C factor calculation. In addition, during the calibration process of the unmodified RUSLE-TLSD model, this study observed a complex relationship between the accuracy of simulating net water erosion in the control areas of certain hydrological stations and the parameter β. Specifically, the accuracy was found to be positively correlated with β in some control areas, while in others, it displayed a negative correlation (Table S4). On the whole, the unmodified RUSLE-TLSD model tends to overestimate net water erosion in the UYRB and underestimate it in the MYRB.

In this situation, this study modified the RUSLE model and coupled it with TLSD model to conduct a comprehensive assessment of net water erosion on a large basin scale. The results of accuracy evaluation revealed that, although the accuracy of the modified RUSLE-TLSD model in the YRB was lower than that achieved in small-scale net water erosion studies, the overall simulation outcomes were reasonably reliable. Calibration stations achieved an NSE of 0.5766, while validation stations achieved an NSE of 0.5369. In assessing net water erosion within the control areas of each hydrological station, it was noted that 8 stations displayed NSE values exceeding 0.3 (Figure S4), indicating the general effectiveness of the modified model.

4.2. In-Depth Analysis of Net Water Erosion Pattern and Driving Factors

From 2000 to 2020, net water erosion in different regions of the YRB showed different spatiotemporal change characteristics. This was closely related to precipitation, topography, vegetation cover, and soil texture [54,55]. In the UYRB, despite the relatively low vegetation coverage (Table S5), the main soil types were Alpine soil, arid soils, and Caliche soils (Figure S5). This region experienced relatively weak erosion due to the scarcity of erosive rainfall. In contrast, the MYRB featured Skeletol primitive soils (Figure S5), characterized by surface vertical and horizontal gullies, rugged terrain, and concentrated summer rainfall. Although the “Grain for Green” program has achieved its set objectives, brought ecological benefits, and significantly increased the vegetation cover of the region (Table S5), the water erosion intensity in this region remained higher than that observed in the UYRB.

Furthermore, the natural ecology of the UYRB, notably the Qinghai–Tibet Plateau, exhibited fragility and a limited capacity to withstand anthropogenic influences [56]. The increase in vegetation coverage in this region was gradual (slope = 0.0018, Table S5). However, the UYRB experienced an increasing trend in rainfall (Slope = 5.4503, Table S5), resulting in an increase in erosion intensity. In the MYRB, particularly within the Loess Plateau, fractional vegetation cover exhibited a significant upturn. The augmentation in fractional vegetation cover played a crucial role in reducing the kinetic energy of raindrops upon reaching the surface [57]. This weakened the impact of surface runoff on sediment transport [58], thus serving as the predominant factor behind the decline in net water erosion in the MYRB.

During the evaluation process, we considered the influence of rainfall, topography, NDVI, and LUCC, etc., on net water erosion. The OPGD results showed that slope and NDVI emerged as the primary controlling factors. In previous studies [59], it was explicitly established that vegetation coverage played a central role in governing potential water erosion in the YRB. This was due to the fact that runoff, driven by gravity, accumulates from higher elevations to lower ones. In regions characterized by elevated terrain and significant topographical variations, longer slope lengths and steeper slopes inherently generate greater kinetic energy, facilitating sediment transport in flowing water [19]. In the process of sediment transport, the use of erosion control measures could reduce the river sediment load [1]. Within the MYRB_TNI, the explanatory capacity of NDVI reached 34.89%, whereas the explanatory capacity of NDVI∩elevation and NDVI∩slope for net water erosion reached 52.83% and 48.52%, respectively. NDVI∩slope exhibited a more robust explanatory power for net water erosion compared to other factor combinations.

Despite the precipitation being a direct factor leading to water erosion, in this study, its explanatory power for net water erosion patterns was not outstanding. To further analyze the reasons for this situation, this study analyzed the relationship between net water erosion patterns and the multi-year average precipitation, multi-year average NDVI, and slope (Figure 14). The results indicated that in the UYRB, the net water erosion intensity was higher in areas with NDVI values less than 0.4, a precipitation of 200–400 mm or greater than 800 mm, and slopes greater than 10°. In the MYRB, the net water erosion intensity was higher in areas with NDVI values of 0.2–0.6, slopes of 5–30°, and precipitation of 400–600 mm. The increase in precipitation did not directly lead to an increase in net water erosion intensity, especially where NDVI was greater than 0.6 or slope was less than 5°, indicating that vegetation and terrain largely limit the impact of rainfall on net water erosion patterns. Moreover, areas with high rainfall often had higher vegetation coverage. When using LUCC data, it was found that areas with lower vegetation coverage were mainly located in the Gobi Desert region, where there was less precipitation and the basic conditions for strong water erosion were lacking. Overall, there were some mutual influences between the driving factors that further affect the water erosion pattern. The redistribution of precipitation by NDVI (vegetation) and slope greatly limited the impact of precipitation on the net water erosion space pattern in the YRB.

Erosion-driving factors exhibited variations, not only across different spatial patterns (Figure 11 and Figure 12) but also in different timescales (Figure 10 and Figure 13). Regarding spatial patterns, it was evident that as the scale of analysis diminished, the geographical characteristics and anthropogenic influences within a region became more distinct, thereby elucidating the mechanisms driving erosion with greater clarity [60]. On a monthly scale, alterations in erosion-driving factors can be attributed to seasonal fluctuations in rainfall and vegetation. On an annual scale, the change in vegetation coverage and net water erosion significantly affects the explanatory power of driving factors (Table S6). It is important to note that as net water erosion intensity decreases, the effectiveness of a single erosion control measure becomes more limited. Consequently, focusing erosion control efforts on regions with high erosion potential or on areas prone to erosion can help to achieve water and soil conservation goals with fewer control measures.

4.3. Policy Implications

This study found complex interactions between net water erosion intensity and driving forces within the YRB. It became evident that each driving factor exerted influence on the distribution of net water erosion and that change in net water erosion engendered variations in the efficacy of each driving factor. Furthermore, it is crucial to recognize that soil conservation measures not only directly curtail water erosion and sediment delivery but are also intricately linked to regional climate, grain production, the carbon cycle, etc. [61,62]. It is imperative to consider the geographical and ecological attributes unique to different YRB regions in future endeavors related to ecological protection and land use planning. These efforts should aim to bolster the water conservation capacity in the UYRB and address erosion issues in the MYRB. Moreover, when examining the interaction between erosion and other elements, it is vital to adopt a comprehensive perspective to analyze the impact of relevant erosion control projects or policies on basin-wide coordination. The analysis of influences beyond the basin potentially requires a broader temporal and spatial framework and research perspective.

4.4. Limitations and Research Prospects

In this paper, the modification of the RUSLE-TLSD model was realized by incorporating the M factor, initially meeting our accuracy requirements. Soil is the direct object of hydraulic erosion, and its own characteristics will have a certain impact on the intensity of water erosion. Soil undergoes the process of being transported and deposited after water erosion, and the soil thickness affects the amount of material involved in the whole process. Therefore, this study mainly utilized soil thickness data to construct the M factor. Although it is commonly believed that thicker soil layers tend to have better erosion resistance, the results of this study show that in the Yellow River Basin, the thicker the soil layer, the higher the erosion modulus. Therefore, the suitability of the modified model in other basins needs to be further experimentally verified and explored by us and by related scholars. However, this corrective idea undoubtedly provides an effective reference for other studies.

Due to the difficulty of assessing the distribution of check dams and the model’s insufficiency, this study did not include conservation projects in the analysis and did not consider the impact of freeze–thaw erosion and wind erosion on sediment data. In addition, this study unifies the data to a resolution of 250 m, resulting in certain size of error. In future, the distribution dataset of erosion and sediment control projects should be further collected, and high-spatial-resolution and large-scale water erosion research will be carried out based on GEE, etc.

5. Conclusions

The main aim of this study was to explore a modified RUSLE-TLSD model suitable for the YRB based on measured sediment data from hydrological stations. It was carried out in order to simulate and analyze the spatial pattern of soil net water erosion from 2000 to 2020, and to analyze the driving factors affecting the net water erosion pattern at different spatiotemporal scales based on the OPGD model. The main conclusions are as follows:

(1)

The M factor is explored for model correction based on soil thickness data and sediment observation data, and the modified RUSLE-TLSD model achieved sufficient accuracy in the YRB. The NSE values at calibration station and verification station increase to 0.5766 and 0.5369, respectively.

(2)

The net water erosion rate within the YRB fluctuates between 1.62 t/(ha·a) and 5.33 t/(ha·a) from 2000 to 2020. Furthermore, the multi-year average net water erosion rates are 3.12 t/(ha·a) for the YRB, 1.66 t/(ha·a) for the UYRB, and 5.27 t/(ha·a) for the MYRB. On an annual scale, net water erosion rates display an overall decrease within the YRB and UYRB, while within the UYRB, they experience an overall increase trend. The period from June to August shows the highest intensity of water erosion. Spatially, regions with pronounced erosion intensity are primarily concentrated on the Loess Plateau. The Hetao Palin and Guanzhong Plain, as important areas for grain cultivation, are mainly characterized by sediment deposition.

(3)

Over the period from 2000 to 2020, the explanatory power of the driving forces behind net water erosion displayed a trend of decreasing fluctuations. The multi-year average value of the factor detection results indicates that, in the YRB, UYRB, MYRB, UYRB_TNI, and MYRB_TNI, slope (10.84%), soil thickness (7.71%), NDVI (13.62%), slope (15.46%), and NDVI (34.89%), respectively, exert the strongest explanatory influence on net water erosion. In the UYRB, terrain-related factors dominate in terms of shaping the spatial framework of net water erosion, whereas in the MYRB, human activities play a more significant role. The interactions between NDVI and slope have a dominant impact on net water erosion and can explain 22.71–52.83% of net water erosion. In addition, the redistribution of precipitation by vegetation and terrain greatly limits the impact of precipitation on the net water erosion space pattern in the YRB.