In this study, the soil loss driven by rill and inter-rill erosion was estimated using the Revised Universal Soil Loss Equation (RUSLE) model, which enables the spatial pattern of soil erosion to be estimated (Renard et al., 1991, 1997). This model has been widely used in watersheds throughout the world (Thomas et al., 2018; Batista et al., 2019). Liu et al. (1994) adapted the model to be more applicable to the Loess Plateau region by considering erosion on a steep slope. Researchers have also improved the accuracy of the model under extreme rainfall events by introducing precipitation concentration degree (Xu et al., 2022). To summarize, the total erosion of hillslopes mainly comprises rill, inter-rill, and gully erosion, and the erosion rate can be calculated using the following equation: T E m = R I m * K * L S * C m * P * α , (1) where TE is the total erosion amount per unit area in a given timestep (t·ha⁻¹); m is the calculation month; RIm is the rainfall erosivity factor introducing the rainfall concentration (MJ·mm·ha⁻¹ h⁻¹); K is the soil erodibility factor (t·h·MJ⁻¹ mm⁻¹); LS is the slope length and gradient factor (dimensionless); and C m and P are the vegetation cover factor and soil conservation factor, respectively (dimensionless). Cai et al. (2020) determined that multi-year average gully erosion accounted for 49% of the total erosion in the Wuding River basin using a model with a physical mechanism, and α is the amplification coefficient used to represent the gully erosion amount, with a value, in this study, of 1.96. The specific calculations for each factor can be referenced in the Xu et al. (2022).

Borselli et al. (2008) proposed an index of connectivity (IC) describing the hydrological linkage between sediment sources and sinks. The IC consists of the following two components: an upslope component ( D u p ) that represents the potential for down-routing at a given location and a downslope component ( D d n ) that accounts for potential flow sinks between that location and the stream network (Vigiak et al., 2016). The calculation formula is as follows: I C = ⁡ log 10 D u p D d n = ⁡ log 10 W ¯ * S ¯ * A ∑ i = 1 n d i / W i * S i . (2) The range of the IC is [−∞, +∞], and a higher value indicates a higher degree of connectivity; W ¯ is the average weight factor of the contributing area, which reflects the surface roughness to measure the resistance of runoff passing through a location; S ¯ is the average slope of the upslope contributing area (m/m); A is the contributing area (m²); d i is the flow length to the downstream main channel (m); W i is the weight factor of the calculated cell; S i is the slope in the calculated cell (m/m); i refers to the calculated cell; and n is the total cell number from the point to the stream network along the flow path.

The structural SDR (HSDR) is calculated using the IC (in Eq. 2) with the following function, which is now included in the InVEST model (Borselli et al., 2008; Vigiak et al., 2012; Cavalli et al., 2013): H S D R = S D R max 1 + ⁡ exp I C 0 − I C i / k , (3) where S D R max is the maximum theoretical structural sediment delivery ratio, which is taken as 1 in this study; I C i is the connectivity index value of each calculated cell; and I C 0 and k are the calibration parameters.

In the main sediment yield area of the Yellow River basin, where soil and water conservation measures have been implemented, soil erosion occurred but cannot transport the eroded sediment to the channels. Zhang (2017) pointed out that the eroded sediment can enter the channels when the daily rainfall reaches 25 mm/day; otherwise, it is trapped in the hillslope system. In addition, the benefits of soil and water conservation measures may be reduced and the deposited sediment be transported again when daily rainfall reaches 50 mm/day (Zhang, 2017). Consequently, we developed a transport threshold factor (Er) according to the rainfall thresholds of 25 and 50 mm/day. The definition formulae are as follows: E r j = 0   P j < 25   m m / d a y , E r j = 1   25   m m / d a y ≤ P j < 50   m m / d a y , E r j = P j / 50 β   P j ≥ 50   m m / d a y , (4) where P j is the daily rainfall; j is the calculation day; and β is the adjustment coefficient. Then, the Er factor is calculated as the ratio of the sum of the daily rainfall erosivity factor multiplied by the E r j during the calculation month to the monthly rainfall erosivity. The Er factor raster is obtained by interpolation, using the inverse distance weighting method.

The movement of sediment from the erosion source to channels is influenced by structural factors, such as topography and soil conservation measures, and dynamic factors, such as rainfall (Zhang, 2017). Therefore, we propose a dynamic SDR ( S D R s l p ) model by considering the influence of the structural characteristics of the hillslopes and transport threshold on the sediment. The equation for this is represented as follows: S D R s l p = H S D R * E r , (5) where HSDR is the structural sediment delivery ratio, reflecting the potential transport capacity of hillslope sediment into the channel system; and Er (in Eq. 4) is the transport threshold factor, reflecting the transport power of rainfall events on eroded soil.

The watershed can be divided into hillslopes and channels. In this study, the hillslope sediment yield ( S Y s l p ) was defined as the amount of sediment entering the channel system from the hillslope system over a specific period of time. The hillslope SDR is the ratio of S Y s l p to the total erosion amount (TE) on hillslopes. To analyze the applicability of the dynamic SDR model, we coupled the hillslope SDRs (HSDR and S D R s l p ) and soil erosion model to simulate the hillslope sediment yield. The potential sediment yield on the hillslope ( P S Y s l p ) was defined based on TE (in Eq. 1) and HSDR (in Eq. 3), represented as follows: P S Y s l p = ∑ i = 1 N T E * H S D R i . (6)

S Y s l p was calculated based on the TE and S D R s l p (in Eq. 5) using the following equation: S Y s l p = ∑ i = 1 N T E * S D R s l p i , (7) where i is the calculated raster and N is the total number of the raster.

Calibration of hillslope sediment yield models remains a challenge, owing to a lack of long-term measured hillslope data (Wen and Deng, 2020). The measured sediment yield of a watershed ( S Y w s d ), which is the total amount of sediment transported through the observed cross-section in the channel, provides an alternative for calibration and validation. However, deposition and erosion may still occur after the sediment-laden flow in the hillslope enters the channel system. In this study, the values of S Y w s d were selected for model calibration and validation in the Mahuyu and Heimutouchuan watersheds when only channel erosion occurred based on the sediment-carrying capacity.

As for the sediment-carrying capacity of the channel, many scholars have pointed out that a linear relationship exists between the flow and sediment concentration for saturated flows in hilly loess areas (Zhang et al., 2018; Zheng, 2018). In this study, the observed monthly maximum daily flow (Q) and monthly S Y w s d were used for curve fitting analysis. The points on the upper side of the fitted curve were then selected for the next fitting analysis until a well-fitted linear curve was obtained, which was regarded as a saturated flow status. The linear relationship was then considered the estimation equation for the sediment-carrying capacity of the channel.

The hillslope sediment yield model based on the S D R s l p was calibrated using the observed S Y w s d from 2006 to 2018 at the Mahuyu hydrological station. For validation, the hillslope sediment yield model was applied to the simulation in the Heimutouchuan watershed (Figure 1D). The S Y w s d from 2006 to 2018 at the Dianshi hydrological station was used as the simulated S Y s l p comparison. Model performance was evaluated using the coefficient of determination (R ²), percentage deviation (PBIAS), and Nash–Sutcliffe efficiency (NSE) coefficients (Nash and Sutcliffe, 1970; Gupta et al., 1999; Yesuf et al., 2015).

The relationship between the monthly maximum daily flow and sediment-carrying capacity of the channel was evaluated by three-time curve fitting. The fitting process and results are shown in Figures 3A‒C.

The power relationship between the flow and sediment yield in unsaturated flows is shown in Figures 3A, B, while the linear relationship between the flow and sediment yield in saturated flows is shown in Figure 3C. It can be seen from Figure 3C that 16 points were selected, which is seen as the saturated flow status, and there was a linear relationship between the monthly maximum daily flow and monthly sediment yield of the watershed. The equation was S Y w s d = 0.527 *Q − 0.498, and the determination coefficient (R ²) was 0.995. Therefore, this fitted relationship can be used as a tool for estimating the sediment-carrying capacity.

The months with the sediment-carrying capacity greater than simulated S Y s l p in the Mahuyu watershed were selected as calibration months. The model parameters were then calibrated by comparing the observed S Y w s d in the calibration months during the period 2006–2018. The results are shown in Figure 4A. The main parameters involved in the simulation of the HSDR, IC ₀ and k, were −4 and 4, respectively, which are determined by referring to the results of the HSDR in the watershed near the Mahuyu watershed (Zhao et al., 2020), and then, the parameter involved in the Er factor, β, was 1.2.

It can be seen from Figure 4A that simulated S Y s l p had a good agreement with the observed S Y w s d in the calibration months. The evaluation index R ² was 0.663, NSE was 0.589, and PBIAS was 34.2%, indicating good performance in hillslope sediment yield prediction. For validation of model applicability using the same parameters, the months with the sediment-carrying capacity greater than simulated S Y s l p at the Dianshi station of the Heimutouchuan watershed were selected. The validation results are shown in Figure 4B. The evaluation indices R ², NSE, and PBIAS were 0.575, 0.537, and 26.6%, respectively. The values of evaluation indices met the standard, indicating that the parameters were reasonable, and the dynamic hillslope sediment yield model was applicable in this watershed.

The watershed average TE, P S Y s l p (in Eq. 6), S Y s l p (in Eq. 7), HSDR, and S D R s l p were calculated at a monthly scale from 2006 to 2018. To analyze the influence of the Er factor on the hillslope sediment yield, 28 months with erosive rainfall but no transport rainfall events and 9 months with monthly average rainfall >50 mm were selected. The variations in the sediment yield and delivery ratio were plotted, as shown in Figures 5 and 6, respectively.

As shown in Figure 5, the average HSDR of the watershed was 0.535, and this essentially remained unchanged. P S Y s l p was >0, and the maximum value was 2.3 t/ha, while S D R s l p and S Y s l p were 0, which indicates that there was eroded sediment on the hillslope of the watershed in these months, but not enough transport power. According to the S Y w s d data obtained from the Mahuyu station, S Y w s d was 0 in 86% of the months covered by Figure 5, and the maximum value was only 0.06 t/ha, which indicates that the calculated S D R s l p was rational. The transport rainfall threshold of 25 mm in the Er factor was shown to be acceptable and able to reflect the situation where erosion occurred but no eroded sediment entered the channels.

As shown in Figure 6, the variation in S D R s l p was striking when the monthly average transport rainfall was >50 mm. S D R s l p was larger than HSDR, and S Y s l p was larger than P S Y s l p . This was especially true for July 2017, where P S Y s l p was 55.9 t/ha, while S Y s l p reached 115.5 t/ha. According to the observed rainfall data, the maximum average daily rainfall of the seven rainfall stations in the Mahuyu watershed was 98 mm in July 2017, and the maximum daily rainfall at the Guoxingzhuang rainfall station was 115.6 mm. Many scholars have obtained hillslope sediment yield via the inversion of dam sediment retention in regions with similar underlying surfaces near the Mahuyu watershed. For example, the hillslope sediment yield in the Chabagou watershed ranged from 107 to 490 t/ha under the heavy rainfall events of 26 July 2017 (Bai et al., 2020), and the S Y s l p in the Wangmaogou watershed reached 253 t/ha under the 148-mm rainfall of the same day (Gao et al., 2018). Comparing the above findings with our simulated results, it can be seen that the simulated S Y s l p based on S D R s l p is closer to the actual hillslope sediment yield than is the P S Y s l p that is calculated by the HSDR. Thus, S D R s l p can more accurately reflect the sediment delivery characteristics under key rainfall events in this watershed.

On the basis of the monthly output raster results, the annual TE, P S Y s l p , S Y s l p , and their accumulated values were tallied, and then the related annual HSDR and S D R s l p were calculated. The results are shown in Figures 7A, B. In terms of annual variation, S Y s l p was significantly larger than P S Y s l p in 2017, while they were more similar in other years. S D R s l p varied markedly on the annual scale during the study period, ranging from 0.231 to 0.878. As can be seen from Figure 7B, S D R s l p had no notable change, indicating that hillslope sediment delivery was relatively stable on the multi-year scale.

To analyze the influence of the Er factor on the spatial variation in sediment delivery characteristics, we chose July 2006 as the representative month due to its uneven spatial distribution of rainfall. The spatial distributions of TE, P S Y s l p , and S Y s l p in July 2006 were plotted, as shown in Figure 8. It can be seen that the amount of soil erosion and sediment yield in the Mahuyu watershed varied dramatically in spatial distribution. The average P S Y s l p and S Y s l p of the watershed were 23.7 and 25.5 t/ha, respectively; and the difference between them was not significant. However, compared with P S Y s l p , the percentage of S Y s l p area <10 t/ha and >150 t/ha increased by 4.6% and 1.5%, respectively, indicating that the S D R s l p , which considers the Er factor, can better reflect the key area of sediment yield.