3. Results and discussion

3.1. Soil erosion rate

Under a rainfall intensity of 90 mm/h, the soil erosion variation range rate for the four soil types was 1.8-35.4 g min-1m-2 (SL), 27.1-135.4 g min-1m-2 (LL), 16.5-76.0 g min-1m-2 (ML) and 105.2-467.1 g min-1m-2 (HL), having the order of SL < LL ≈ ML < HL (Fig. 2). When rainfall intensity was 120 mm/h, the rate of soil erosion for the four soil types increased by 11.9 times (SL), 1.6 times (LL), 2.3 times (ML) and 1.4 times (HL). Erosion rate at SL was the most sensitive to an increase in rainfall intensity, and the overall order of soil erosion rate was ML ≈ LL < SL < HL. These results indicate that rainfall intensity plays a major role in soil erosion. Soil erosion rate with slope recorded an initial increase for SL and HL before decreasing as the angle of slope increased; maximum erosion rates were recorded at 15° and 20°, respectively. The soil erosion rates for LL and ML had no regularity with an increase in slope angle.
In order to further analyze causes for different types of soil erosion, differences in erosion rate before rills were formed were analyzed. Results indicate that, although no significant differences in SL, LL or ML were recorded, a significant difference between SL, LL, ML and HL was evident (Fig. 3). As shown in Fig. 2, significant differences among the four soils during the whole soil erosion process were also recorded. Comparison of these results indicate that differences of SL, LL and ML soil erosion were mainly caused by differences in rill erosion. These results also indicated that an increase of soil erosion rate caused by different types of soil rill erosion was not obvious when rainfall intensity was 90 mm/h. When rainfall intensity was 120 mm/h, SL rill erosion caused a sharp increase in the soil erosion rate; LL rill erosion caused a strong increase in the soil erosion rate; and ML rill erosion caused a slight increase in the soil erosion rate. An increase in the soil erosion rate caused by HL rill erosion was not obvious.

3.2. Morphological development and occurrence of rills

3.2.1. Occurrence of rills

Runoff characteristics and rill generation time (Table 2) indicate that, when rainfall intensity was 90 mm/h, rills were evident in all slopes under the HL soil type. For SL, LL and ML, rills were only recorded with slope angles of 15°, 20° and 25°, respectively. The occurrence time of rills gradually decreased with a decrease in particle size, a finding that is consistent with the variation law of runoff generation time (Table 2). Under a rainfall intensity of 120 mm/h, rills appeared in HL and SL at all four slopes; LL recorded an absence of rills at 15° and 25°, and ML had no rills at 10°. Rill frequency results indicate that HL was the most prone to the formation of rills, and rills occurred more frequently at SL, LL and ML under heavy rain, with SL being the most sensitive to heavy rain. Under different rainfall intensities and slopes, a decrease in loess particle size resulted in a gradual decrease in runoff generation time, and average runoff rates before rill formation gradually increased (Table 1 and Table 2). However, the uniformity of soil particle size composition and the size of aggregates affects topsoil stability (Neyshabouri et al., 2011; Madenoglu et al., 2020). In general, the stability of topsoil decreases with a decrease in uniformity of soil particle size composition and an increase in aggregate size (Neyshabouri et al., 2011). Our results indicated that particle size composition of ML and HL was uniform (Table 1), and that particle size of aggregates was larger. Although these characteristics were conducive to runoff stripping and scouring, organic matter content in ML was higher than in the other soil types, thus enhancing surface soil stability. Before rill formation, HL recorded rates of runoff and soil erosion. A large volume of sediment was stripped and transported, resulting in the formation of an uneven slope, thereby creating conditions for concentrated flow. HL was also recorded to have the highest clay content of the soil types, thus being more readily stripped and scoured, thus being the main reason why this soil type was prone to rill formation.
As previously recorded (He et al. 2017), rills are mainly caused by heavy rainfall intensity; consistent with our results, rills typically increase with an increase in rainfall intensity (Berger et al., 2010). Rainfall intensity affects soil erosion through raindrop strike and runoff (Vaezi et al., 2018). Large rainfall intensity increases instantaneous runoff, and an increase in raindrop strike increases runoff disorder, which increases runoff erosivity, thus increasing the probability of rill generation. By comparing runoff rate and runoff velocity under different types of soil before rill formation indicated that rill formation is more sensitive to an increase in runoff rate caused by an increase in rainfall intensity. Under a rainfall intensity of 90 mm/h, SL runoff was low, and runoff was difficult to converge, therefore not being conducive to the generation of rills. Under a rainfall intensity of 120 mm/h, the increase in raindrop splash erosion and runoff promoted the removal and transportation of sediment, resulting in the formation of an uneven slope. As SL contained the largest sand content, a loose structure and the lowest content of organic matter, this soil type had poor structural stability. Once uneven points are present in a soil, it will quickly collapse, promoting the generation of rills. Heavy rainfall intensity also increased stripping and transportation of sediment under LL and SM soil types. However, due to LL and ML having higher clay and organic matter contents, soil structural stability of these soil types increased, resulting in a reduction of soil collapse. At the same time, because silt and sand contents were still large, erosion and scouring were not easily undertaken, thus rills did not readily form under LL and ML soils. Even when rills did form under these soil types, formation took a long time to occur.
Although SL, LL and ML rills did not record an obvious change with slope gradient, rills still mainly occurred on slopes with a greater gradient, indicating that slope gradient was an important factor affecting rill development (Mancilla et al., 2005). Slope affects rill generation by affecting rainfall, the runoff dynamic gradient and soil stability (Bagarello et al., 2010; Langhans et al., 2014; Jerzy et al., 2016). In general, the hydraulic gradient of runoff from a steep slope is large. This increases the velocity of slope flow and increases the shear force of slope runoff; at the same time, the instability of a steep soil slope increases, and the component force along the slope surface increases. Sediment is therefore removed by runoff, thus being conducive to the occurrence of rills (Komatsu et al, 2011; Ries et al., 2014; Zhang et al., 2016). However, under the influence of reduced the area of land experiencing rainfall, the impact of slope on runoff has a critical value. The sensitivity of different types of soil to a reduction of rainfall differs, and the occurrence time and probability of rills vary with slope angle.

3.2.2. Morphological development of rills

Different indicators are currently used to quantify the development process of rills from the aspects of rill cross-section morphology, rill network characteristics and connectivity. In this study, cumulative rill length (CDL), average rill width (WA), average rill depth (Ha), rill density (DS), rill merging node (J) and rill number (N) were selected from these aspects to quantify rill morphology variation with rainfall. These variables were also used to clarify differences in rill development processes for the different types of soil, such as headward erosion, rill wall collapse and rill bottom undercutting.
Results for change process of CDL with rainfall (Fig. 4) indicated that SL, LL and ML had an approximate linear increase with an increase in rainfall; the increment of increase, however, was small. CDL of HL recorded an initial rapid increase as rainfall increased; after this the rate of increase declined. The increment for HL was large. At the end of the rainfall experiment, CDL range for SL, LL and ML was 1.5-4.8 m, and there was no obvious sequence relationship among these soils. CDL range for HL was 8.2-41.7 m. CDL for SL increased with an increase in rainfall intensity, and CDL changes for LL, ML and HL were almost not affected by rainfall intensity. CDL for SL initially increased before decreasing with an increase in slope, reaching a maximum when the angle of slope was 15°. CDL for LL and ML did not record a change with an increase in slope angle; CDL results for HL recorded an increase as the angle of slope increased.
Results for Wa change process with rainfall (Fig. 5) recorded an increase with an increase in rainfall duration. Results for SL and LL recorded an initial rapid increase before the rate of increase declined; ML recorded a rapid increase at 20° and 25° under a rainfall intensity of 120 mm/h, having an initial rapid increase and then increasing slowly under other conditions. At the end of the rainfall experiment, Wa for the different soil types was: HL < ML ≈ LL < SL, in which Wa of SL ranged from 12.6 - 24.2 cm, ML and LL ranged from 3.4 - 15.0 cm, and HL ranged from 1.6 - 2.5 cm. The influence of rainfall intensity on Wa of SL and HL was not obvious; Wa of LL and ML increased with an increase in rainfall intensity. Wa recorded a decrease for SL with an increase slope angle; LL and ML recorded an increase in Wa with an increase in slope angle; and HL recorded no change.
Ha recorded a general trend of increase with rainfall (Fig. 6). Ha results for SL, LL and HL all recorded a rapid initial increase before the rate of increase declined; Ha for ML recorded an initial rapid increase under a rainfall intensity of 120 mm/h. At the end of the rainfall experiment, no obvious sequence relationship was identified between Ha under all soil types in the range of 2.5-10.7 cm. The influence of rainfall intensity on Ha under SL, ML and HL is therefore not obvious. Ha under LL recorded an increase with rainfall intensity. Ha under SL recorded an initial decrease before increasing as the angle of slope increased, with the maximum value occurring when the angle of slope was 10°. Ha of LL and HL increased with an increase in slope angle. Under a rainfall intensity of 120 mm/h, Ha under ML decreased with an increase in slope angle.
The change process of DS with rainfall (Fig. 7) recorded an increase as rainfall duration and rainfall intensity increased for all soil types. However, under a rainfall intensity of 120 mm/h, the increasing trend of DS under different types of soil gradually declined. At the end of the rainfall experiment, DS for the different soil types was in the order of: ML < SL < LL < HL. DS ranged from 0.2 - 2.2 m/m2 for SL, LL and ML, and from 3.2 - 7.6 m/m2 for HL. DS of SL initially increased before decreasing as the angle of slope increased, reaching a maximum at 20°. DS under LL and HL increased as the angle of slope increased, and DS under ML was irregular.
Results for the change process of J with rainfall (Fig. 8) indicated that J under SL easily occurred under a rainfall intensity of 120 mm/h and a slope angle ranging from 10 ° to 20°. Results indicated that J did not easily form under LL and ML soil types; J also easily formed under HL. When rainfall intensity was 90 mm/h, J under HL rapidly increased with an increase in rainfall duration. When rainfall intensity was 120 mm/h, J recorded a rapid initial increase before increasing slowly with rainfall duration under SL and HL soil types. At the end of the rainfall experiment, the variation range of SL and HL was 2-6 and 5-24, respectively. Although J recorded an increase with an increase in slope intensity under SL and HL, it was not readily apparent at 25° under SL.
Change processes of N under all soil types recorded an increase with rainfall duration (Fig. 9), being characterized by an initial rapid increase which then declined. HL recorded the most N, ranging from 9 to 40. There was no obvious order relationship among SL, LL and ML, with N ranging from 1 to 7. Results for rainfall intensity recorded N to increase under all soil types. With an increase in slope angle, N under SL initially increased before decreasing, peaking at 20o. N under LL decreased as slope angle increased. N under ML was predominantly not affected by slope, and N under HL recorded an increase as the angle of slope increased.
Our results indicate that the development process of rill morphological parameters under different soil types are similar, recording rapid increases under a light rain intensity and a gentle slope, however the level of increment was small. Under a heavy rain intensity and a steep slope, rill morphological parameters initially rapidly increased before the rate of increase declined; the level of increment of rapid increase was larger. In the early stage of rill development, the converging effect of rill development on runoff gradually increased, and converging runoff further promoted the rapid development of rills. With the continuous development of rills, the converging effect of rills on runoff gradually weakened, decreasing rill development. Our results are consistent to those of Wu et al. (2018) who examined feedback coupling effects between rill morphology and rill erosion. Our results also indicated that morphological characteristics of rill development in different soil types differ. The development of rills under SL was mainly due to an increase of Wa and J; rill development under HL was mainly due to an increase of CDL, J and N. The index parameters of LL and ML rill development were between SL and HL. At the end of the rainfall experiment, Ha had no obvious order relationship, however the speed of increase differed during the rainfall events. The main reason for this phenomenon is that rill widening is mainly caused by collapse of the rill wall and merging of adjacent rills under gravity (Mancilla et al., 2005; Shen et al., 2015). An increase in rill length is mainly caused by an increase in the number of rills and rill head advance under traceable erosion (Schneider et al., 2013; Wu et al 2018). An increase in rill depth is mainly caused by further retrogressive erosion of rills due to undercutting, caused by runoff shear dispersion and the reemergence of rill head undercutting in rills (Gómez et al., 2003; Zhang et al., 2017; Zhao et al., 2018). A higher level of soil cohesion, due to higher levels of clay and organic matter content, make it easier for soil to form a mass structure, especially as clay content in a soil can significantly enhance the anti-dispersion ability of the wet soil layer. Soil structure stability is the important factor here. Therefore, Wa of HL slopes was smaller, being less than 5 cm. Under the condition that rills easily occur, stripping erosion due to runoff further intensifies the headward erosion of rills. SL soil only had a clay content of 12.1% and the lowest organic matter content of the four soils, resulting in this soil to have poor structural stability. In addition, the sand content of SL peaked at 68.5%, resulting in the soil to be loose and porous. Once a rill formed in SL, the rill wall and gully head readily collapsed under runoff effects, thus aggravating slope erosion. Due to differences in rill development mode between SL and HL, Ha in SL was mainly developed by further retrogressive erosion of the rill head after rill collapse, characterized by a phased increase; Ha in HL was mainly developed by rill bottom cutting whilst the rill length increased, characterized by a gradual increase. At the same time, due to differences in rill development mode between SL and HL, the negative effect of SL rill development on runoff convergence slowly weakened, and the negative effect of HL rill development on runoff convergence weakened faster. Therefore, after rill formation, SL erosion is more severe and lasts longer.
Results from this analysis indicate that rainfall intensity has a positive effect on rill morphological parameters for all soils, and sensitivity of rill parameters to rainfall intensity was: N > DS > J > Wa > CDL ≈ Ha. The angle of slope had both positive and negative effects on SL, LL and ML rill morphological parameters, with only positive effects on HL. The sensitivity of rill parameters to slope response had the order of: Ha > Wa ≈ J > DS ≈ N > CDL.

3.3. Correlation between influencing factors and relevant parameters of erosion processes

In order to examine the influence of influencing factors on slope erosion, Pearson correlation analysis was used on influencing factors, runoff and sediment parameters, and rill shape parameters (Table 3 and Table 4). Results indicate that clay, silt and MWD were significantly positively correlated with the amount of erosion and rill erosion, and sand and OM was significantly negatively correlated with the amount of erosion and rill erosion (Table 3). The correlation of soil characteristic factors with rill morphological parameters indicated that clay and silt are significantly positive correlated with CDL, DS, J and N. Significant negative correlations were recorded between clay and Wa, and a significant negative correlation was recorded between silt and Wa, Ha; a significant positive correlation between sand and Wa, Ha was recorded, and a significant negative correlation was recorded between sand and CDL, DS, J and N. A significant positive correlation was found between OM and MWD and CDL, DS, J and N. finally, a significant negative correlation existed between OM, MWD and Wa, Ha. These findings correlate the rationality of our previous analysis.
As soil characteristics also affect slope roughness after erosion, they also affect flow velocity between interrills and rill, thereby having an important indirect impact on the development of erosion and rill morphology. As shown in Table 4, clay, silt and MWD were significantly positively correlated with Q and V, and significantly negatively correlated with T1 and VL. Sand was significantly positively correlated with T1 and VL and significantly negatively correlated with Q and V. OM was positively correlated with T1 and V, and negatively correlated with Q. Wa and Ha were negatively correlated with V, CDL, DS and J; N was negatively correlated with VL and positively correlated with Wa.

3.4. A multi-type soil erosion prediction model based on rill development

Correlation analysis shows that soil characteristic factors seriously affect slope erosion, thus the establishment of a quantitative index of soil characteristics in this study took erosion amount as the objective function. Soil characteristic factors with a significant correlation with erosion amount were selected to establish a regression relationship, thereby forming a comprehensive quantitative parameter (C). By using multiple linear regression in SPSS, the silt factor was eliminated due to collinearity, and the equation was established as:
C=4.220 clay+0.392 sand+49.986 MWD-10.55 OM R2 =0.814 (1)
This analysis indicates that the change of erosion amount has a synchronous effect with rill morphology evolution. In this study, we adopted the same method as that used by Zhang et al. (2019), constructing the rill morphology expression (G) as:
G=6.271 CAL+1.108 Wa+2.100 DS+0.862 N R2 =0.926 (2)
Based on the constructed soil characteristic factor C and rill form G, combined with the main factors of rainfall, rainfall intensity and slope that affect soil erosion, the prediction model of slope erosion in the loess area was established as:
R2=0.919 (3)
where, D is the erosion module of individual rainfall under a bare slope (kg m-2); P is rainfall (mm); I is rainfall intensity (mm/h); and S is slope (°). The model was validated using 78 independent data sets (Fig. 10), and R2 between the measured value and the simulated value was 0.919, indicating that the equation could suitably predict slope erosion.