A field investigation on rill development and flow hydrodynamics under different upslope inflow and slope gradient conditions

S

Slope gradient (sine of the slope angle), %

q

Inflow discharge per unit width of the plot, L min⁻¹m⁻¹

h_r

Rill flow depth, m

v_s

Surface rill flow velocity, m s⁻¹

v_r

Mean rill flow velocity, m s⁻¹

α

Correction factor in determining mean rill flow velocity

g

Acceleration due to gravity, m s⁻²

Re

Reynolds number

Fr

Froude number

f

Darcy–Weisbach resistance coefficient

τ

Shear stress, N m⁻²

W

Stream power, N m⁻¹ s⁻¹

φ

Unit stream power, m s⁻¹

SLR

Soil loss rate, g m⁻² min⁻¹

RER

Rill erosion rate, g m⁻² min⁻¹

TRL

Total rill length, m

RD

Rill density, m m⁻²

MRW

Mean rill width, cm

MRD

Mean rill depth, cm

R_WD

The rill width–depth ratio

Soil erosion is a critical concern for the sustainable development of agricultural regions worldwide (Bennett et al. 2015), and soil erosion could cause non-point source pollution (Wang et al. 2019a, 2019b), soil organic carbon loss (Fang et al. 2018), and biodiversity reduction (Li et al. 2020). Thus, soil erosion is a serious threat to the environment and agricultural security (Jiang et al. 2020). Rill erosion is one of the main patterns of soil erosion by water on sloping croplands and rangelands (Bewket & Sterk 2003; Kimaro et al. 2008; Wu et al. 2019).

Rill erosion primarily results from concentrated surface flow (Kimaro et al. 2008; He et al. 2017), which is the main source of sediment yield in hillslope erosion processes (Cerdan et al. 2002; Wagenbrenner et al. 2010; Mirzaee & Ghorbani 2018). Generally, rill erosion occurs when the erosivity of flowing water exceeds a certain threshold of soil resistance; as a result, the micro-terrain of hillslope surface changes and gradually forms a rill (Govers et al. 2007; Knapen et al. 2007; Wang et al. 2015). The main factors affecting soil detachment and sediment transport on hillslopes include slope gradients (Koulouri & Giourga 2007; Fu et al. 2011), overland flow velocity, and flow depth (Gimenez & Govers 2002). During rill erosion processes, the evolution of rill networks greatly affects the confluence path of runoff, soil loss, and the micromorphology of the slope surface (Consuelo et al. 2007; Auerswald et al. 2009; Fang et al. 2015; Tian et al. 2017). Especially on steep loess hillslopes with sparse vegetation cover in northern China (Han et al. 2018; Yinglan et al. 2019), rills develop rapidly, causing high soil loss (Lei et al. 2001; Zhang et al. 2017), and there is a lack of effective prevention methods to control rill erosion on bare steep slopes in this region (Kou et al. 2020). Therefore, a further investigation of rill erosion processes on steep hillslopes is urgently needed.

Rill erosion may be affected by hydrological parameters, such as rill flow stream power (Mahmoodabadi et al. 2014a), runoff rate in rills (Wirtz et al. 2012), hydraulic slope, and rainfall intensity (Shen et al. 2016), as well as topographic factors (Kinnell 2009; Mahmoodabadi et al. 2014b). The rill flow has a significant effect on the transportation of eroded soils, which will accordingly affect the soil erosion rate (Di Stefano et al. 2017). For rills formed at the hillslope, they quickly adapted their topographic characteristics in response to the changes of upslope inflow discharge and slope gradient (Nearing 1997; Nord & Esteves 2010; Stefano et al. 2013). Moreover, the erosivity of flowing water is closely related to upslope inflow discharge and slope gradient, thus, upslope inflow discharge and slope gradient are two important factors for rill erosion (Berger et al. 2010; Sun et al. 2013; Shen et al. 2015).

The slope gradient is an important topographic factor which has direct effects on flow erosivity and infiltration, thus indirectly influences rill erosion processes (Shen et al. 2016; Zhang et al. 2017). Giménez et al. (2004) indicated that slope roughness is enhanced with the slope gradient becoming steeper, whereas flow velocity decreases until a critical hydraulic condition (Froude number is 1) of the flow is reached. It has also been recognized that the critical flow hydraulic condition for rill initiation is closely correlated to the slope gradient (Ziadat & Taimeh 2013; Ran et al. 2018). In addition, many researchers found an increase in the slope gradient can lead to an increase in rill erosion amount (Berger et al. 2010; Fang et al. 2015). Ziadat & Taimeh (2013) reported that soil erosion in uncultivated land was mostly influenced by slope gradient; thus, the determination of slope gradient factor was an essential part of rill erosion prediction (Sun et al. 2013). However, Liu et al. (1994) observed that the relationship between slope gradient and soil erosion varied when slope gradient reached certain thresholds. For instance, Yair & Klein (1973) found the amount of soil erosion decreased with the increase of slope gradients. Koulouri & Giourga (2007) stated that when the slope reached 40%, the dominant factor for soil loss was no longer land use, but slope gradient. In addition, the study in a vineyard area by Pijl et al. (2020) indicated that rittochino (vertical cultivation) with a slope of 39.5% had the most severe soil loss.

Typically, power functions are commonly used to fit the correlation between the slope gradient and the amount of soil erosion (Di Stefano et al. 2018). In some experiments conducted on gentle slopes, significant relationships were found between rill erosion and slope gradient (Liu et al. 1994; Sirjani & Mahmoodabadi 2012), whereas in some other studies, no significant relationship was reported (e.g., Chaplot & Le Bissonnais 2003). Zhang et al. (2002) found the soil detachment rate could be well predicted by a power function of upslope inflow discharge and slope gradient, and detachment rates were more sensitive to upslope inflow discharge than to slope gradients. Therefore, there is an interaction effect between slope gradient and upslope inflow discharges on the amount of erosion and the critical hydraulic conditions of rill initiation. Thus, it is necessary that more attention should be focused on the effects of upslope inflow discharges and slope gradients on rill flow hydraulics and its dynamic mechanism at steep hillslopes.

Rill network and its morphological characteristics play an important role in determining the flow characteristics of surface runoff and sediment delivery in hillslope erosion processes (Zhang et al. 2016; Wu et al. 2018). In fact, the formation and development process of rill morphology timely responds to the rill erosion process. After the emergence of rills, with the continuous flowing water, rills develop and show the phenomenon of bifurcation, merging and connecting on the slope, and finally evolve into rill networks (He et al. 2014; Shen et al. 2015). The indicators describing rill morphology are usually rill length, rill width, rill depth, and rill networks (Bewket & Sterk 2003; Raff et al. 2004). After the formation of rill networks, overland flows concentrated into rills can promote the development of rills, broadening the rill width and deepening the rill length (Tian et al. 2017). Furthermore, the formation and development of rills contribute to the runoff connectivity and concentration along the rill networks (Romero et al. 2007; Heras et al. 2011). However, the influence mechanisms of inflow discharge and slope gradient on rill morphological characteristics and rill network formation need to be further investigated.

It is important to clarify the hydrodynamic mechanism of rill erosion since sediment detachment and transport by overland flow are energy-consuming processes. Rill initiation is determined by the established threshold flow hydraulic conditions, and further rill development is controlled by different thresholds (Wirtz et al. 2012). Rill erosion usually increases with increasing upslope inflow discharge due to the increase of flow hydrodynamics (Di Stefano et al. 2017), such as flow pattern (Foster et al. 1984), flow velocity (Govers et al. 2007), depth (An et al. 2012), shear stress (Nearing et al. 1991), and stream power (Berger et al. 2010). Hillslopes' soil loss significantly increases after rill formation due to increases in rill flow depth, rill flow velocity, and shear stress compared to sheet flow (Liu et al. 2011). Numerous studies have verified that the rill flow hydraulic parameters are related to upslope inflow discharge and slope gradient, whereas the relationships are different in the studies (Mirzaee & Ghorbani-Dashtaki 2018). For instance, Giménez & Govers (2001) indicated that rill flow velocity is independent of the slope gradient on the mobile bed, while it becomes slope-dependent on the fixed bed. The flow in rills has higher velocities, depth and sediment transport capacity than the flow in the inter-rill areas (Gatto 2000; Mirzaee & Ghorbani-Dashtaki 2018). Furthermore, the rill bed form varies with erosion processes, which in turn, changes the hydraulics and erosion capacity of rill flow (Zhang et al. 2017; Jiang et al. 2018).

In summary, it is still not thoroughly understood what mechanisms are dominating the inflow discharge and slope gradient on the rill erosion process and its hydrodynamic mechanisms, and the rill development and its morphological characteristics, especially for steep hillslopes. Therefore, in this study, we performed field plot experiments under combinations of different upslope inflow discharges and slope gradients at steep loess hillslopes. The objectives of this study are to: (1) quantify the effects of slope gradient and inflow discharge on rill erosion processes and rill morphological characteristics; (2) evaluate rill flow hydraulic characteristics and the hydrodynamic mechanism of rill erosion.

The study site is located in the hilly loess area of northern China, where the Guanting Reservoir is located (Figure 1).

This region belongs to a transition zone between the northwest plateau of Beijing and the North China Plain, which is one of the key areas implementing soil and water conservation in China. The Guanting Reservoir was built in the early 1950s and long served as an important drinking water source for Beijing City. The hilly areas account for approximately 37% of the total area of Guanting Reservoir watershed (38°15′–41°14.2′ N, 112°8.3′–116°20.6′ E), and most of the hilly areas are highly erodible loess, covering a total area of 500 km² (Hu & Wang 2004). The main causes leading to severe soil erosion here include erodible loessial soils, relatively sparse vegetation cover, and short-duration rainstorms (Tian et al. 2017). The climate is semi-arid continental, with a mean annual temperature of 9 °C (Li & Pan 2018). Mean annual precipitation, most of which occurs in the summer months as short-duration and high-intensity rainstorms, varies from 370 to 480 mm. Field runoff plot experiments were performed at the Huailai Field Experimental Station (40°15′32″ N, 115°37′ 01″ E), Beijing Normal University, which is located along the east side of the Guanting Reservoir (Figure 1).

A fallow loess hillslope (35 m length and 25 m width) was selected as the experimental field (Figure 1). The particle size distribution and basic physical and chemical properties of the soils in the 0–20 cm surface layer are presented in Tables 1 and 2, respectively.

The runoff plot experiments were conducted on the selected hillslope. Lei et al. (2001) found sediment transport capacity would not increase beyond certain values of rill length and slope. Furthermore, considering the characteristics of small erosion plots, we designed the runoff plot of 5 m in length and 2 m in width, which is bounded by concrete walls or aluminum plastic plates. A specifically designed PVC pipe with a row of equally spaced (1.0 cm) grooves (0.5 cm width, 1.5 cm arc length) was fixed at the top of the runoff plot and set along the longitudinal direction of the PVC pipe to make sure the upslope inflow water flowed uniformly over the whole runoff plot. Meanwhile, a PVC pipe was set at the bottom of the runoff plot to collect runoff sediment samples during the experiment. The schematic diagram of the runoff plot and its affiliated facilities is shown in Figure 2.

In this study, we monitor the comparison results by setting the slope gradient at three levels, and inflow discharge at six levels. For each slope gradient, two runoff plots with the same size and treatment were constructed at the selected hillslope. In total, six runoff plots with the same size of three different slope gradients were used in our experiments. All the experiments were repeated twice to ensure the reproducibility of results. Given the fact that this experiment was performed in the same area as that by Tian et al. (2017), six values of the inflow discharge per unit width (1 m) of the plot ranging from 6 to 36 L min⁻¹m⁻¹ were used. Each inflow discharge was introduced at the head of the plot. Considering the characteristics of hillslope erosion driven by short-duration rainstorms in the local area, the simulated scouring test lasts 60 minutes for each experimental treatment. Rill erosion was the most significant on bare loess slopes with slope gradients of 17.6–57.7% (Shen et al. 2015). The rill is prone to intersect on loess hillslopes with a slope gradient over 15° (26%), which is generally considered as the lower limit of steep slopes in the hilly loess area of China. Although 25° (42%) is the upper limit of the slope gradient for cultivation in China set by the government, the phenomenon of cultivation in steeper slopes is still common in the hilly loess area (Li & Pan 2018). Therefore, to investigate rill development and erosion processes driven by upslope inflow on a wide range of steep hillslopes, slope gradient (sine of the slope angle) of 26, 42, and 57% were selected. According to the characteristics of short-term rainfall–runoff erosion in this area, the duration of each trial lasted 60 minutes after the initiation of surface runoff at the plot outlet.

Before the experiment, upslope inflow discharge was calibrated to the designed values manually. In addition, the topsoils (0–20 cm) were sampled by cutting rings on the top, middle, and bottom of the slope surface to measure the soil water content and dry bulk density. For each treatment, the antecedent soil water content and soil bulk density were controlled as 13 ± 2% and 1.3 ± 0.1 g cm⁻³, respectively, before the start of the experiment. To reduce the usage of backfill for the surface soil treatment, experiments with inflow discharges ranging from low to high were conducted.

Runoff sediment samples were collected from the outlet of the plot at intervals of 4 minutes, and the duration of collecting each sample was recorded by a stopwatch. These samples were used to calculate soil loss rates and surface runoff rates during the experimental process. Surface flow velocities in rills (v_s) were measured by the dye (KMnO₄) tracing method with the slope segment of 0.5–1.5 m, 1.5–2.5 m, 2.5–3.5 m, and 3.5–4.5 m, respectively. The water depth in rills (h_r) was measured perpendicularly to the rill bed using a thin steel ruler with 0.1 mm precision and the measurement sections were the same as the velocity measurements. The measurements of rill morphology (length, width, depth) were performed along each rill at a 5 cm interval to compute the rills' volume. Then, the rill erosion amount can be obtained by multiplying the rills' volume with soil bulk density. In addition, photographs were taken during the experimental process. The maps of rill networks at the end of the experiment were plotted, which included three steps. First, as soon as the scouring test was over, we used a high-definition digital camera to take photos of the entire plot. Second, according to the distribution of rills formed on the plot, we drew a sketch of the rill network on white paper according to a certain scale. Third, for each rill, we measured its length and the width every 3–5 cm along the length of the rill. Last, the maps of rill networks were made by Adobe Photoshop, combining the photos from the above three steps.

The amount of rill erosion was estimated by volumetric measurements of rills on each soil erosion plot. In our experiments, hillslope soil loss is the sum of the amount of sheet erosion and rill erosion. Thus, sheet erosion was obtained by subtracting the rill erosion from the hillslope soil loss. The statistical tests include variance analysis, correlation analysis, equation simulation, and were conducted using the IBM SPSS Statistics for Windows Version 19.0. The data graphs were drawn using ORIGIN 9.0.

Figure 3 shows the temporal variations of surface runoff rates at the runoff plot outlet during the experimental process under different slope gradients and inflow discharges.

According to Figure 3, there were no significant differences in the trend of surface runoff rate over time under different treatments. For each treatment, surface runoff rates increased with time in the initial stage of the experiment, and then gradually became steady. With the same slope gradient, the increasing speed of surface runoff rate was higher under higher upslope inflow discharge conditions. This indicated the temporal variation increased with the increase of inflow discharges. With the same inflow discharge, the surface runoff rates were higher under higher slope gradient conditions. This indicated that, to a certain extent, soil infiltration rates decreased with increasing slope gradients. Overall, under the same inflow discharge condition, there were no significant differences in surface runoff processes among different treatments.

To quantify the overall effect of slope gradients and inflow discharges on surface runoff production, the mean runoff coefficients during the whole experimental process for different treatments are presented in Figure 4.

Under the same slope gradient condition, the mean runoff coefficients slightly increased with the increasing inflow discharges (Figure 4), indicating soil infiltration rates decreased as inflow discharges increased. Under the same inflow discharge condition, the mean runoff coefficients increased with the increase of slope gradients, suggesting the soil infiltration rates decreased as slope gradients increased. However, the rate of increment in mean runoff coefficients with slope gradients gradually decreased as inflow discharges increased. This phenomenon implied that the decreasing effect of slope gradients on soil infiltration rates became weakened with the increases of inflow discharges.

To explore the temporal variations of soil loss during the experiment, the changes in soil loss rates with duration for different inflow discharges and slope gradients are illustrated in Figure 5.

According to Figure 5, the soil loss rates between different treatments were of obvious differences, especially for different inflow discharges. Under different conditions, the starting time for rill initiation varied. Specifically, the starting time of rill formation became earlier with the increasing slope gradient and inflow discharge, which indicated that rills formed more easily under steeper slopes and high inflow discharges. In addition, the duration of rill initiation reduced as the slope gradient increased. For most treatments, the soil loss rate curves were high in the early minutes during the experimental process, which was mainly because the loose soil particles were relatively abundant on the hillslope surface at the beginning of the experiment, leading to numerous sediments being detached and transported by surface sheet flow, especially under high inflow discharges. This phenomenon corresponds to the experimental findings of Parsons & Stone (2006).

Overall, for most experimental treatments, two peak values of soil loss rates existed during the whole experimental process. At the initial stage of the experiment, there was much infiltration on the 57% slope and the erosion capacity of the surface runoff is limited. In addition, the soil loss rate was relatively low at the beginning of the experiment, as shown in Figure 5(c). After 8 minutes, there appeared a first peak value of soil loss rate, which was because the runoff increased and there were abundant erodible soil particles on the hillslope surface. Thus, the soil loss rate also increased rapidly afterwards. As the experiment continued, the soil loss rate decreased with the reduction of erodible soil particles. During 20 to 35 minutes, with the active evolution of rill erosion, the soil loss rate increased to the second peak value. After 35 minutes, the rill development entered a relatively stable stage, and the soil loss rate slowly decreased.

Under the lowest inflow discharge (6 L min⁻¹m⁻¹), the curves of soil loss rate generally declined with slight fluctuation during the entire experiment. For relatively high inflow discharges, the curves of soil loss rate sharply decreased to the valley at approximately 10 to 20 minutes, which was mainly due to decreases in the detachment rate of soil particles and the number of loose particles on the slope surface. Then, the soil loss rate curves reached another crest at approximately 25 to 45 minutes after the run, which was mainly due to the active development of rills, contributing to a sharp increment in the soil loss rate. Thereafter, the soil loss rate curve descended and tended to be an equilibrium since the development of rills gradually entered into a relatively steady stage.

For each slope gradient treatment, rill erosion amount, hillslope soil loss, and the contribution of rill erosion to hillslope soil loss gradually increased with the increase of inflow discharges (Table 3). This phenomenon indicated that the erosion capacity of flowing water greatly increased as inflow discharges increased. For the lowest slope gradient condition (26%), sheet erosion increased with increasing inflow discharges. However, with the slope becoming steeper, sheet erosion increased as inflow discharges increased up to 30 L min⁻¹m⁻¹, and then decreased as inflow discharges continued to increase.

For the lowest slope gradient condition, the contribution of rill erosion to hillslope soil loss was 41.0% and 57.0%, respectively, under inflow discharges of 18 and 24 L min⁻¹m⁻¹. This result showed that rill erosion became the dominant soil erosion pattern instead of sheet erosion when inflow discharge increased to 24 L min⁻¹m⁻¹. Therefore, we can conclude that there was a threshold of inflow discharge ranging from 18 and 24 L min⁻¹m⁻¹, under which, the dominant soil erosion pattern started to change from sheet erosion to rill erosion. With the slope gradient increasing to 42%, rill erosion became the dominant erosion pattern when inflow discharge was 18 L min⁻¹m⁻¹, which indicated that the critical inflow discharge controlling the main soil erosion pattern was between 12 and 18 L min⁻¹m⁻¹. However, with slope gradient increasing to 57%, rill erosion became the dominant soil erosion pattern when inflow discharge was 12 L min⁻¹m⁻¹, which demonstrated that the critical inflow discharge controlling the main soil erosion pattern was between 6 and 12 L min⁻¹m⁻¹. In summary, the threshold inflow discharge that controlled the dominant soil erosion pattern was related to slope gradients. Overall, the critical inflow discharge that determined the dominant soil erosion pattern decreased as slope gradients increased.

For the same inflow discharge, with increasing slope gradients, both rill erosion and the contribution of rill erosion to hillslope soil loss generally increased, whereas the variations of sheet erosion showed no definite trend with slope gradients. However, hillslope soil loss significantly increased with slope gradients varying from 26% to 42%, and then slightly varied as slope gradients ranged from 42% to 57%. The variations in the contribution of rill erosion to hillslope soil loss and hillslope soil loss with slope gradients showed that the increase of slope gradients enhanced rill erosion, although it implied that there was probably a threshold of slope gradient where soil erosion began to weaken. For example, rill erosion increased from 3.28–81.47 kg to 5.19–90.46 kg with the increase of slope gradient from 42% to 57%, whereas hillslope soil loss changed very slightly (15.04–99.33 kg and 14.01–106.43 kg at slope gradients of 42 and 57%, respectively).

According to Equations (11) and (12), both soil loss and rill erosion rates increased by a power function with an increase in either inflow discharge or slope gradient. Compared with the exponents of inflow discharge, the slope gradients were lower, which indicated that the inflow discharge played a more important role in both soil loss and rill erosion rates than the slope gradient. In addition, both the exponents of inflow discharge and slope gradient in Equation (11) were lower than those in Equation (12).

To explore the characteristics of rill networks developed at different slopes, the rill networks at the end of the experiment for different inflow discharges and slope gradients are illustrated in Figure 6 

Based on the rills generated on the slope surface at the end of the experiment (Figure 6), for each slope gradient treatment, both the uniformity and density of rill network distribution were enhanced with the increases of inflow discharges.. In addition, the cross-merging phenomenon between the rills became more obvious under higher inflow discharges, which indicated that the rill network became more complex. Under the same inflow discharge condition, the density of rill network increased with slope gradient increasing from 26% to 42%, and then decreased with the slope becoming steeper. This phenomenon was more significant under higher inflow discharges. As the slope gradient increased from 42% to 57%, the total length of rills decreased by 8.6–25.1% under inflow rates of 12–36 L min⁻¹m⁻¹, while the average rill depth increased by 17.5–35.7%. Although the rill density decreased, the average depth of rills increased. This phenomenon was probably different from previous studies, and was mainly related to our experimental condition. In this study, rills were formed only by the upslope inflow without raindrop impacts. Under this circumstance, when the slope gradient increased to a certain level (57%), the undercutting erosivity of runoff increased significantly with the influence of gravity, making the rills continue to deepen. As a result, more runoff flowed into the deep rills and the number of rills formed was reduced compared to the 42% slope condition. Therefore, the total length of rills formed under the 57% slope condition was less than that under the 42% slope, which was reflected in the decrease of rill density. This may also imply that there was a critical slope gradient between 42 and 57%, under which, the development of rill networks began to weaken. For the largest slope (57%) condition, rills on the upper slope were significantly denser than those on the lower part of the slope, especially under relatively low inflow discharges. This phenomenon was probably due to the strong erosion capacity of flowing water with low sediment concentration in the upper part of the slope on the steep hillslopes.

As reported in a previous study (Zhang et al. 2017), the formation and evolution of drop pits on the slope surface was an essential procedure during rill development processes. Their finding has agreed with our experiments, in which the drop pit phenomenon was clearly observed at the initial stage of rill formation during the experimental process for all treatments. For the 26 and 42% slope, all rills were generally continuous at the end of the experiment. However, on the steeper slope (57%), there were some drop pits and discontinuous rills at the lower-middle part of the runoff plot under relatively low inflow discharges at the end of the experiment. A possible explanation for this phenomenon was related to the deep rills developed at the upper part of the slope. In this case, most of the surface runoff converged into the deep rills, and intensively scoured the deep soil layer. During this process, some of the surface runoff might turn into the interflow, which decreased the surface upslope inflow discharge over the lower part of the runoff plot. Thus, the drop pits and discontinuous rills at the lower part of the slope were unable to develop into continuous rills due to the weakened energy of flowing water.

Furthermore, based on our observation, the gravitational erosion obviously occurred within the plot of 57% slope under upslope inflow rates of 24, 30, and 36 L min⁻¹m⁻¹. Under these circumstances, we can conclude that gravity plays an important role in slope erosion and rill evolution, which causes an increase in the collapse of the rills. In contrast, under other slope gradients and inflow rate conditions, the main driving force of soil erosion and rill development is the upslope inflow. For these treatments, the role of gravity on erosion and rill evolution was not obvious.

Different statistical indicators of rill morphology for different slope gradients and inflow discharges are presented in Table 4.

Based on the dynamic changes of total rill length, mean rill width, and the rill width–depth ratio with increasing slope gradients (Table 4), the vertical erosion in rills was exaggerated with the slope becoming steeper since the gravitational erosion occurred. This result is similar to He et al. (2016), in which the authors investigated the dynamic changes of total rill length under different simulated rainfall conditions. For each slope gradient treatment, total rill length, rill density, and mean rill depth gradually increased as inflow discharges increased, while the increment rate declined with the increases in inflow discharges. However, under the same inflow discharge condition, mean rill width generally showed a decreasing trend with the increases of slope gradients. Furthermore, the rill width–depth ratio varied slightly with the increase of inflow discharges, demonstrating no obvious variation trend. This result is similar to the laboratory rill experiments conducted by Shen et al. (2016) and He et al. (2016). According to a least significant difference (LSD) test, mean rill depth under the lowest slope gradient was significantly different from that under the highest slope gradient condition, whereas mean rill width and total rill length among three different slope gradient treatments were not significantly different.

Under the same inflow discharge condition, rill density first increased with the increase of slope gradients up to 42% or 57%, and then decreased as the slope gradient continued to increase. For each inflow discharge treatment, although both mean rill width and rill depth increased with increasing slope gradients, the rill width–depth ratio showed a decreasing trend as slope gradients increased. This phenomenon indicated that with the increase of slope gradients, deepening rills were more significant than the widening rills because of the higher shear stress of concentrated rill flow. The findings were consistent with the observations in rill experiments under field conditions by Mirzaee & Ghorbani-Dashtaki (2018). For the same inflow discharge condition, mean rill width and rill depth under different slope gradients were not significantly different according to an LSD test. However, the rill density was significantly different among six different inflow discharges according to an LSD test. For the highest slope gradient, the side wall collapse of rills was more significant than that of lower slope gradients, mainly due to the obvious phenomenon of gravity erosion under steep slopes.

During the experimental process, one of the main sediment delivery processes was observed to be headwall retreat erosion at small knickpoints or steps, which contributed to the formation of continuous rills. Furthermore, for steeper slopes (57%), bank collapse processes of rills due to the interactions of flowing water and gravity were another significant erosion pattern resulting in soil loss. These phenomena were consistent with the observations for field rill experiments under inflow conditions by Wirtz et al. (2012) and Tian et al. (2017). Therefore, the limitation of soil detachment or sediment transport equations was that important sediment preparing processes frequently occurred at steep hillslopes, such as headcut retreat, step-pool effects, and bank collapse, which were not adequately considered (Dong et al. 2019; Zegeye et al. 2020). Through the observations during the experiment, sediment yield at the runoff plot outlet primarily resulted from knickpoints, chutes, meanders, and sidewall sloughing. This was consistent with the findings by Stefanovic & Bryan (2009), who investigated rill development processes on loamy sand and sandy loam under laboratory conditions.

According to the analysis above, the rill morphology was influenced by the interaction effects of inflow discharges and slope gradients. At the same time, the regression analysis also indicated that rill morphological indicators could not be well expressed by only inflow discharges or slope gradients. The best fitting functions for rill morphological indicators, inflow discharges (q), and slope gradients (S) are presented in Table 5.

Rill density, mean rill depth, and the rill width–depth ratio can be described by the power function of inflow discharges and slope gradients, whereas mean rill width was well expressed by the linear function relationship of inflow discharges and slope gradients (Table 5). All the fitted equations were highly significant (p < 0.01). In addition, the coefficient of the inflow discharge was greater than that of the slope gradient for all equations, which indicated that the inflow discharge factor played a more important role in the rill morphological indicator than the slope gradient factor. The function form of mean rill width, inflow discharge, and slope gradient was consistent with the result by Lei & Nearing (2000), who conducted flume experiments under upslope inflow conditions.

Rill flow hydraulic parameters over the whole course of the experiment for different inflow discharges and slope gradients are presented in Table 6.

Overall, mean rill flow velocity (v_r) increased with increases in either slope gradient or inflow discharge (Table 6); however, there were significant differences in v_r among six different inflow discharges, but slight differences among three different slope gradients. This result indicated that flow velocities in rills were mainly dependent on the magnitude of inflow discharges, which is similar to the observations by Nearing (1997), and Li et al. (2016), who found that flow velocities were largely independent of slope gradients.

Mean rill flow depth (h_r) increased with increasing inflow discharges but decreased as slope gradients increased. Reynolds number (Re) increased as inflow discharges increased, whereas Re increased and then decreased with increases of slope gradients, in which the slope gradients increased from 26% to 42% and from 42% to 57%, respectively. For all slope gradients, the Re value began to exceed 2,000 when the inflow discharge increased to 24 L min⁻¹m⁻¹, which implied that rill flows became turbulent flow according to the criterion of open channel flow (Tian et al. 2017). Fr was greater than 1.0 for all treatments, indicating that rill flows belonged to supercritical flow conditions. This observation is similar to the experimental results for steep calcareous soils (Mirzaee & Ghorbani 2018). In addition, Darcy–Weisbach resistance coefficient (f) decreased with the increases of inflow discharges but increased as the slope became steeper.

The variations of rill erosion rates with three hydrodynamic parameters of the rill flow, respectively, under different slope gradients are illustrated in Figure 7.

The Pearson correlation coefficients between rill morphological indicators and rill flow hydrodynamic parameters for all treatments are shown in Table 7.

There was a negative correlation between the flow resistance coefficient and all the rill morphological indicators (Table 7). However, RD, MRW, and MRD were positively correlated with the other hydrodynamic parameters. Specifically, R_WD was positively related to mean rill flow depth and Reynolds number but negatively correlated with other parameters. The most sensitive rill flow hydrodynamic parameter to RD, MRW, MRD, and R_WD was Reynolds number, mean rill flow depth, mean rill flow velocity, and Froude number, respectively. Overall, the common sensitive parameter to various rill morphological indicators was shear stress and stream power.

According to Equation (17), when no rills formed, i.e., MRD= 0, then the critical rill flow velocity was determined. Thus, in our experiments, the critical flow velocity for rill occurrence could be estimated as 0.224 m s⁻¹. Furthermore, this critical value was close to the measured flow velocity at the initial stage of rill development during the experiment (He et al. 2018).

Plot experiments under six upslope inflow discharges (6–36 L min⁻¹m⁻¹) and three slope gradients (26–57%) were conducted to explore the quantitative effects of inflow discharge and slope gradient on rill erosion and rill morphology, as well as the hydrodynamic mechanism of rill erosion. The results showed soil loss rates significantly related to the fluctuations during the process of rill development. Both rill erosion and its contribution to hillslope soil loss increased with the increases in inflow discharge and slope gradient, which can be described by a power function. The inflow discharge played a more important role than the slope gradient. There was a threshold of inflow discharge (12–24 L min⁻¹m⁻¹), under which, rill erosion instead of sheet erosion became the dominant erosion pattern, and the critical inflow discharge reduced as slope gradient increased. The uniformity of rill network distribution enhanced with increasing inflow discharges. Rill density increased (17.7–25.4%) with slope gradient increasing from 26% to 42%, and then decreased (7.5–25.1%) with the slope becoming steeper. Mean rill width increased with the increases of inflow discharge but decreased as the slope gradient increased. However, mean rill depth increased with increases in inflow discharge and slope gradient. The rill width–depth ratio decreased as slope gradients increased, which indicated that deepening rills were more significant than the widening rills with the slope becoming steeper. In addition, stream power was the best hydrodynamic parameter correlated with rill erosion with a good linear relationship, and the critical stream power of rill erosion was 0.255 N m⁻¹s⁻¹. The most sensitive hydrodynamic parameter to rill density, mean rill width, rill depth, and the rill width–depth ratio was Reynolds number, rill flow depth, rill flow velocity, and Froude number, respectively. Although the most sensitive hydrodynamic parameters to various rill morphological indicators were different, the common sensitive parameters were rill flow stream power and shear stress. The findings will help to understand the impacts of upslope inflow rate on morphological variations of rills, rill network development, and hydrodynamic mechanism of rill erosion on steep hillslopes. Weakening the erosion power of upslope runoff, such as taking measures to cover the soil surface, is useful for preventing rill formation and development, as well as the downslope soil loss.

The authors wish to thank all of the technicians involved in field and laboratory work. This research is jointly funded by the Natural Science Foundation of China Project (41907061), the Research Center on Mountain Torrent & Geologic Disaster Prevention of the Ministry of Water Resources, Changjiang River Scientific Research Institute (CKWV2019761/KY), and the Research Program from the State Key Laboratory of Soil Erosion and Dryland Farming on the Loess Plateau (A314021402-2005). We declare that we have no financial and personal relationships with other people or organizations that can inappropriately influence our work. The authors declare no conflicts of interest.

All relevant data are included in the paper or its Supplementary Information.