Abstract

Few studies focus on the quantitative impact of upslope inflow rate and slope gradient on rill development and erosion processes. Field plot experiments under varying inflow rates (6–36 L min−1m−1) and slope gradients (26, 42 and 57%) were conducted to address this issue. The results showed soil loss rates significantly demonstrated temporal variability in relevance to the rill developing process. Rill erosion and its contribution to soil loss increased with increasing inflow rates and slope gradients by power functions. There was a threshold inflow discharge (12–24 L min−1m−1), under which, rill erosion became the dominant erosion pattern. At the initial stage, downcutting of rill bottom and headward erosion were obvious, whereas rill broadening was significant at the actively rill developing period. Rill density increased with slope gradient increasing from 26% to 42%, and then decreased. For the 57% slope under high inflow rates (24–36 L min−1m−1), gravity caused an increase in the collapse of rills. Mean rill width increased with increasing inflow rates but decreased as slope gradients increased, while mean rill depth increased with increasing inflow rates and slope gradients. Stream power and rill flow velocity were the best hydrodynamic parameter to simulate rill erosion and rill morphology, respectively.

HIGHLIGHTS

  • Rill erosion and its contribution to soil loss increased with inflow rate and slope gradient by a power function.

  • The threshold inflow rate that made rill erosion dominated decreased with increasing slope gradient.

  • There was a threshold slope gradient between 42.3% and 57.4% where the rill network development began to weaken.

  • Stream power was the best parameter simulating rill erosion with a good linear relationship.

SYMBOLS AND ABBREVIATIONS

     
  • S

    Slope gradient (sine of the slope angle), %

  •  
  • q

    Inflow discharge per unit width of the plot, L min−1m−1

  •  
  • hr

    Rill flow depth, m

  •  
  • vs

    Surface rill flow velocity, m s−1

  •  
  • vr

    Mean rill flow velocity, m s−1

  •  
  • α

    Correction factor in determining mean rill flow velocity

  •  
  • g

    Acceleration due to gravity, m s−2

  •  
  • Re

    Reynolds number

  •  
  • Fr

    Froude number

  •  
  • f

    Darcy–Weisbach resistance coefficient

  •  
  • τ

    Shear stress, N m−2

  •  
  • W

    Stream power, N m−1 s−1

  •  
  • φ

    Unit stream power, m s−1

  •  
  • SLR

    Soil loss rate, g m−2 min−1

  •  
  • RER

    Rill erosion rate, g m−2 min−1

  •  
  • TRL

    Total rill length, m

  •  
  • RD

    Rill density, m m−2

  •  
  • MRW

    Mean rill width, cm

  •  
  • MRD

    Mean rill depth, cm

  •  
  • RWD

    The rill width–depth ratio

INTRODUCTION

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.

MATERIALS AND METHODS

Study site

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

Figure 1

(a) The study area and (b) the location of the experimental hillslopes with field plots.

Figure 1

(a) The study area and (b) the location of the experimental hillslopes with field plots.

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 km2 (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).

Experimental plots

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.

Table 1

Particle size distribution for the topsoil soil (0–20 cm) of the sloping field

Particle size distribution (%)
2–1 mm1–0.5 mm0.5–0.25 mm0.25–0.1 mm0.1–0.05 mm0.05–0.02 mm0.02–0.002 mm < 0.002 mm
1.29 ± 0.07 1.21 ± 0.19 3.53 ± 0.66 7.42 ± 0.15 25.34 ± 0.35 30.40 ± 0.62 11.52 ± 0.24 19.50 ± 0.59 
Particle size distribution (%)
2–1 mm1–0.5 mm0.5–0.25 mm0.25–0.1 mm0.1–0.05 mm0.05–0.02 mm0.02–0.002 mm < 0.002 mm
1.29 ± 0.07 1.21 ± 0.19 3.53 ± 0.66 7.42 ± 0.15 25.34 ± 0.35 30.40 ± 0.62 11.52 ± 0.24 19.50 ± 0.59 

As an example, 1.29 ± 0.07 refers to the mean with standard deviations, the same in other cells.

Table 2

Basic physical and chemical properties for the topsoil (0–20 cm) of the sloping field

Surface treatmentSoil bulk density (g cm−3)Soil water content (%)Soil porosity (%)Total N (mg g−1)Organic C (mg g−1)
Bare slope 1.35 ± 0.50 13.0 ± 2.0 44.57 ± 2.32 0.56 ± 0.18 422 ± 5.49 
Surface treatmentSoil bulk density (g cm−3)Soil water content (%)Soil porosity (%)Total N (mg g−1)Organic C (mg g−1)
Bare slope 1.35 ± 0.50 13.0 ± 2.0 44.57 ± 2.32 0.56 ± 0.18 422 ± 5.49 

As an example, 1.35 ± 0.50 refers to the mean with standard deviations, the same in other cells.

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.

Figure 2

Field runoff plot system at the study site.

Figure 2

Field runoff plot system at the study site.

Experimental design

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−1m−1 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.

Experimental procedures and measurements

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−3, 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 (vs) were measured by the dye (KMnO4) 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 (hr) 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.

Data calculation and analysis

The slope gradient (S, %) in this study is the sine of the slope angle (θ, °) for field plots: 
formula
(1)
The surface runoff rate at the plot outlet (qout) is calculated as follows: 
formula
(2)
where qout is surface runoff rate at the plot outlet, L min−1m−1; M is the total mass of a runoff sediment sample, g; m is the dry weight of sediment in a runoff sediment sample, g; ρ is the density of water used in the experiment, g L−1; t is the duration of collecting a runoff sediment sample, min; b is the width of the plot, m; 1,000 is the dimensional adjustment coefficient.
The runoff coefficient (RC) corresponding to each runoff sample is obtained by the following equation (Tian et al. 2017): 
formula
(3)
where qin is upslope inflow discharge (L min−1 m−1).
Mean rill flow velocity (vr) is estimated by the measured surface rill flow velocity (vs) multiplied by a correlation factor α as the equation: 
formula
(4)
where a is 0.67 for laminar flow, 0.7 for transitional flow, and 0.8 for turbulent flow (Horton et al. 1934; Li et al. 1996).
Reynolds number (Re) and Froude number (Fr) can be used to reflect the rill flow regime, which are calculated as: 
formula
(5)
where hr is the mean rill flow depth (m); νm is kinematical viscosity (m2 s−1); and g is acceleration of gravity (m s−2).
The resistance to rill flow is estimated by Darcy–Weisbach friction coefficient (f), which is calculated as: 
formula
(6)
where g is gravitational acceleration (m s−2); J is surface slope (m m−1), which is equal to S in value; hr is the depth of rill flow (mm).
Shear stress (Foster et al. 1984), stream power (Bagnold 1966), and unit stream power (Yang 1996) are used to reflect the hydrodynamic parameters of rill flows, which are calculated as: 
formula
(7)
where τ is flow shear stress (N m−2); W is flow stream power (N m−1s−1); φ is unit stream power (m s−1); ρ is the density of water (kg m−3).
Rill density (RD) is the total rill lengths (TRL) divided by the whole hillslope area (m m−2) (Bewket & Sterk 2003). The RD value is suitable for roughly describing the morphological characteristics of rills developed at hillslopes. The formulas of TRL and RD are: 
formula
(8)
 
formula
(9)
where TRL is the total lengths of all rills on the studied hillslope, m; RD is the rill density, m m−2; TRLi is the total length of a rill and its bifurcations (m); and i= 1, …, n represents the number of the rills; A0 is the plane area of the studied hillslope (m2).
In addition, for a certain cross section of a rill, the rill width–depth ratio (RWD) can reflect the shape of the rill cross section. For a certain cross section of a rill, RWD is calculated as follows: 
formula
(10)
where MRW is the mean rill width, cm; MRD is the mean rill depth, cm; RWD is the rill width–depth ratio.

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.

RESULTS AND DISCUSSION

Temporal changes in surface runoff and soil loss

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.

Figure 3

Temporal changes in surface runoff rates under three slope gradients and six inflow discharges. S stands for the slope gradient. (a) S = 26%, (b) S = 42% (c) S = 57%.

Figure 3

Temporal changes in surface runoff rates under three slope gradients and six inflow discharges. S stands for the slope gradient. (a) S = 26%, (b) S = 42% (c) S = 57%.

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.

Figure 4

Mean runoff coefficients for different slope gradients and inflow discharges.

Figure 4

Mean runoff coefficients for different slope gradients and inflow discharges.

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.

Figure 5

Temporal changes in soil loss rates under different slope gradients and inflow discharges. (a) Slope gradient = 26%, (b) slope gradient = 42%, (c) slope gradient = 57%.

Figure 5

Temporal changes in soil loss rates under different slope gradients and inflow discharges. (a) Slope gradient = 26%, (b) slope gradient = 42%, (c) slope gradient = 57%.

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−1m−1), 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.

Contributions of slope gradient and inflow discharge to rill erosion

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−1m−1, and then decreased as inflow discharges continued to increase.

Table 3

Rill erosion, sheet erosion, hillslope soil loss, and the contribution of rill erosion to hillslope soil loss for different slope gradients and inflow discharges

Slope gradient (%)Inflow discharge (L min−1 m−1)Rill erosion (kg)Sheet erosion (kg)Hillslope soil loss (kg)Contribution of rill erosion to hillslope soil loss (%)
26 1.08 3.83 4.92 22.0 
12 3.16 6.42 9.58 33.0 
18 9.58 13.78 23.36 41.0 
24 18.36 13.85 32.21 57.0 
30 33.97 18.29 52.27 65.0 
36 43.34 19.47 62.81 69.0 
42 3.28 11.76 15.04 21.8 
12 12.40 19.39 31.78 39.0 
18 27.73 17.07 44.80 61.9 
24 39.86 18.30 58.16 68.5 
30 64.89 21.35 86.24 75.2 
36 81.47 17.87 99.33 82.0 
57 5.19 8.83 14.01 37.0 
12 14.61 12.96 27.57 53.0 
18 37.56 16.87 54.43 69.0 
24 56.75 18.92 75.66 75.0 
30 73.10 19.43 92.53 79.0 
36 90.46 15.96 106.43 85.0 
Slope gradient (%)Inflow discharge (L min−1 m−1)Rill erosion (kg)Sheet erosion (kg)Hillslope soil loss (kg)Contribution of rill erosion to hillslope soil loss (%)
26 1.08 3.83 4.92 22.0 
12 3.16 6.42 9.58 33.0 
18 9.58 13.78 23.36 41.0 
24 18.36 13.85 32.21 57.0 
30 33.97 18.29 52.27 65.0 
36 43.34 19.47 62.81 69.0 
42 3.28 11.76 15.04 21.8 
12 12.40 19.39 31.78 39.0 
18 27.73 17.07 44.80 61.9 
24 39.86 18.30 58.16 68.5 
30 64.89 21.35 86.24 75.2 
36 81.47 17.87 99.33 82.0 
57 5.19 8.83 14.01 37.0 
12 14.61 12.96 27.57 53.0 
18 37.56 16.87 54.43 69.0 
24 56.75 18.92 75.66 75.0 
30 73.10 19.43 92.53 79.0 
36 90.46 15.96 106.43 85.0 

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−1m−1. This result showed that rill erosion became the dominant soil erosion pattern instead of sheet erosion when inflow discharge increased to 24 L min−1m−1. Therefore, we can conclude that there was a threshold of inflow discharge ranging from 18 and 24 L min−1m−1, 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−1m−1, which indicated that the critical inflow discharge controlling the main soil erosion pattern was between 12 and 18 L min−1m−1. However, with slope gradient increasing to 57%, rill erosion became the dominant soil erosion pattern when inflow discharge was 12 L min−1m−1, which demonstrated that the critical inflow discharge controlling the main soil erosion pattern was between 6 and 12 L min−1m−1. 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).

Soil loss rate (SLR) and rill erosion rate (RER) can be expressed by inflow discharges (q) and slope gradients (S) as dual power functions, respectively, as follows: 
formula
(11)
 
formula
(12)
where SLR is soil loss rate, g m−2 min−1; RER is rill erosion, g m−2 min−1; q is inflow discharge per unit width, L min−1m−1; and S is slope gradient, %.

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).

Rill networks and morphology

Rill networks

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 

Figure 6

Rill networks at the end of the experiment for different slope gradients and inflow discharges. (a) Slope gradient = 26%, (b) slope gradient = 42%, (c) slope gradient = 57%.

Figure 6

Rill networks at the end of the experiment for different slope gradients and inflow discharges. (a) Slope gradient = 26%, (b) slope gradient = 42%, (c) slope gradient = 57%.

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−1m−1, 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−1m−1. 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.

Rill morphological characteristics

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

Table 4

Rill morphological characteristics for different slope gradients and inflow discharges

Slope gradient (S, %)Inflow discharge (q, L min−1m−1)TRL (m)RD (m m−2)MRW (cm)MRD (cm)RWD
26 12.24 1.22 1.8 0.8 2.31 
12 17.85 1.79 2.9 1.1 2.59 
18 22.52 2.25 4.7 1.7 2.86 
24 27.74 2.77 5.5 2.4 2.27 
30 32.68 3.27 6.2 3.1 2.01 
36 36.36 3.64 6.5 3.4 1.91 
42 15.30 1.53 1.6 1.0 1.58 
12 22.04 2.20 2.5 1.4 1.75 
18 26.64 2.66 3.9 2.2 1.78 
24 33.72 3.37 4.7 2.9 1.63 
30 38.46 3.85 5.6 3.7 1.50 
36 44.28 4.43 6.0 4.0 1.50 
57 17.24 1.72 1.4 1.4 1.05 
12 20.10 2.01 2.0 1.9 1.08 
18 24.64 2.46 3.6 2.6 1.37 
24 28.36 2.84 4.5 3.3 1.38 
30 32.22 3.22 5.3 4.4 1.20 
36 33.15 3.32 5.8 4.7 1.23 
Slope gradient (S, %)Inflow discharge (q, L min−1m−1)TRL (m)RD (m m−2)MRW (cm)MRD (cm)RWD
26 12.24 1.22 1.8 0.8 2.31 
12 17.85 1.79 2.9 1.1 2.59 
18 22.52 2.25 4.7 1.7 2.86 
24 27.74 2.77 5.5 2.4 2.27 
30 32.68 3.27 6.2 3.1 2.01 
36 36.36 3.64 6.5 3.4 1.91 
42 15.30 1.53 1.6 1.0 1.58 
12 22.04 2.20 2.5 1.4 1.75 
18 26.64 2.66 3.9 2.2 1.78 
24 33.72 3.37 4.7 2.9 1.63 
30 38.46 3.85 5.6 3.7 1.50 
36 44.28 4.43 6.0 4.0 1.50 
57 17.24 1.72 1.4 1.4 1.05 
12 20.10 2.01 2.0 1.9 1.08 
18 24.64 2.46 3.6 2.6 1.37 
24 28.36 2.84 4.5 3.3 1.38 
30 32.22 3.22 5.3 4.4 1.20 
36 33.15 3.32 5.8 4.7 1.23 

TRL: total rill lengths of the whole plot, m; RD: rill density, m m−2; MRW: mean rill width, cm; MRD: mean rill depth, cm; RWD: the rill width-depth ratio.

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.

Table 5

The fitted equation of rill morphological indicators, inflow discharges and slope gradients

Regression equationR2Sig.(p)n
RD= 0.346S0.125q0.533 0.92 <0.01 18 
MRW=0.027S+ 0.157q + 1.947 0.96 <0.01 18 
MRD= 0.033S0.519q0.798 0.96 <0.01 18 
RWD = 32.734 S−0.798q−0.016 0.85 <0.01 18 
Regression equationR2Sig.(p)n
RD= 0.346S0.125q0.533 0.92 <0.01 18 
MRW=0.027S+ 0.157q + 1.947 0.96 <0.01 18 
MRD= 0.033S0.519q0.798 0.96 <0.01 18 
RWD = 32.734 S−0.798q−0.016 0.85 <0.01 18 

RD: rill density, (m m−2); MRW: mean rill width (cm); MRD: mean rill depth (cm); RWD: the rill width-depth ratio; S: slope gradient (%); q: inflow discharge (L min−1 m−1).

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.

Hydrodynamic mechanisms of rill erosion and morphology

Rill flow hydraulics

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

Table 6

Mean rill flow hydraulic parameters for different slope gradients and inflow discharges

Slope gradient (%)Inflow discharge (L min−1 m−1)Rill flow velocity (m s−1)Rill flow depth (×10−3 m)Reynolds numberFroude numberFriction coefficient
26 0.284 2.41 544 1.9 1.28 
12 0.368 3.50 1,021 2.0 1.11 
18 0.455 4.57 1,653 2.2 0.94 
24 0.527 5.97 2,498 2.2 0.92 
30 0.607 7.21 3,477 2.3 0.84 
36 0.659 8.86 4,640 2.2 0.87 
42 0.329 2.13 556 2.3 1.46 
12 0.448 3.28 1,166 2.5 1.22 
18 0.564 4.30 1,927 2.7 1.01 
24 0.640 5.51 2,798 2.8 1.00 
30 0.699 6.64 3,690 2.7 1.01 
36 0.733 8.12 4,728 2.6 1.13 
57 0.362 1.99 571 2.6 1.70 
12 0.464 3.13 1,152 2.6 1.63 
18 0.574 3.88 1,768 2.9 1.32 
24 0.629 4.63 2,311 3.0 1.31 
30 0.711 5.40 3,050 3.1 1.20 
36 0.748 6.21 3,689 3.0 1.24 
Slope gradient (%)Inflow discharge (L min−1 m−1)Rill flow velocity (m s−1)Rill flow depth (×10−3 m)Reynolds numberFroude numberFriction coefficient
26 0.284 2.41 544 1.9 1.28 
12 0.368 3.50 1,021 2.0 1.11 
18 0.455 4.57 1,653 2.2 0.94 
24 0.527 5.97 2,498 2.2 0.92 
30 0.607 7.21 3,477 2.3 0.84 
36 0.659 8.86 4,640 2.2 0.87 
42 0.329 2.13 556 2.3 1.46 
12 0.448 3.28 1,166 2.5 1.22 
18 0.564 4.30 1,927 2.7 1.01 
24 0.640 5.51 2,798 2.8 1.00 
30 0.699 6.64 3,690 2.7 1.01 
36 0.733 8.12 4,728 2.6 1.13 
57 0.362 1.99 571 2.6 1.70 
12 0.464 3.13 1,152 2.6 1.63 
18 0.574 3.88 1,768 2.9 1.32 
24 0.629 4.63 2,311 3.0 1.31 
30 0.711 5.40 3,050 3.1 1.20 
36 0.748 6.21 3,689 3.0 1.24 

Overall, mean rill flow velocity (vr) increased with increases in either slope gradient or inflow discharge (Table 6); however, there were significant differences in vr 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 (hr) 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−1m−1, 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.

Correlations of rill erosion and rill flow hydrodynamics

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

Figure 7

Rill erosion rates versus rill flow shear stress, stream power, and unit stream power, respectively, for all treatments.

Figure 7

Rill erosion rates versus rill flow shear stress, stream power, and unit stream power, respectively, for all treatments.

Average rill flow shear stress (τ), stream power (W), and unit stream power (φ) varied from 0.645 to 4.345 N m−2, 0.183 to 3.252 N m−1 s−1, and 0.061 to 0.419 m s−1, respectively, and the values were increasing with the increases of either upslope inflow discharge or slope gradient (Figure 7). Rill erosion rates significantly increased with an increase of τ, W, and φ. The relationships between rill erosion rates and the above three hydrodynamic parameters were described well by the linear functions, which are listed as follows: 
formula
(13)
where RER is rill erosion rate (g m−2 min−1), τ is shear stress (N m−2), 
formula
(14)
W is stream power (N m−1 s−1), 
formula
(15)
φ is unit stream power (m s−1).
The critical hydrodynamic parameter could be obtained by the above functions. When no rill erosion occurred, i.e., RER = 0, then the critical value was determined. All the critical values increased as slope gradient increased. The increase in slope gradient will result in higher values of hydrodynamic parameters, and further enhance the detachment power and transport capacity of rill flows (Wang et al. 2016). As a result, rill erosion rates increased with the increases of slope gradients. Moreover, as the slope becomes steeper, the gravity influenced the flow geometry in rills and increased the rill bed roughness, which led to an increase in the critical shear stress and the critical stream power for rill erosion (Nearing 1997). At the same time, for all treatments, the functional relations between rill erosion rates and those three hydrodynamic parameters are analyzed, which are expressed as: 
formula
(16)
where RER is rill erosion rate (g m−2 min−1), τ is rill flow shear stress (N m−2), W is stream power (N m−1 s−1), φ is unit stream power (m s−1). The critical shear stress, the critical stream power, and the critical unit stream power computed by Equation (16) were 0.916 N m−2, 0.255 N m−1 s−1, and 0.051 m s−1, respectively. In contrast, the best hydrodynamic parameter to fit the rill erosion rate was unit stream power with a good linear relationship.

Correlations of rill morphology and rill flow hydrodynamics

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

Table 7

Correlation coefficients between rill morphological indicators and rill flow hydrodynamic parameters for all treatments

vrhrReFrfτWφ
RD 0.939** 0.919** 0.965** 0.438 −0.548* 0.785** 0.818** 0.478* 
MRW 0.859** 0.954** 0.934** 0.233 −0.764** 0.633* 0.689* 0.298 
MRD 0.971** 0.782** 0.880** 0.661* −0.326 0.944** 0.963** 0.753** 
RWD -0.376 0.871 0.579 −0.785** −0.611* −0.609* −0.533* −0.780** 
vrhrReFrfτWφ
RD 0.939** 0.919** 0.965** 0.438 −0.548* 0.785** 0.818** 0.478* 
MRW 0.859** 0.954** 0.934** 0.233 −0.764** 0.633* 0.689* 0.298 
MRD 0.971** 0.782** 0.880** 0.661* −0.326 0.944** 0.963** 0.753** 
RWD -0.376 0.871 0.579 −0.785** −0.611* −0.609* −0.533* −0.780** 

*: significant at p= 0.05; **: significant at p= 0.01; vr: rill flow velocity; hr: rill flow depth; Re: Reynolds number; Fr: Froude number; f: Darcy–Weisbach friction coefficient; τ: shear stress; W: stream power; φ: unit stream power; MRD: mean rill depth; MRW: mean rill width; RD: rill density; RWD: the rill width-depth ratio.

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, RWD 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 RWD 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.

In addition, there was a good linear function relationship between mean rill depth (MRD) and rill flow velocity (vr), which was expressed as: 
formula
(17)

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−1. 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).

CONCLUSIONS

Plot experiments under six upslope inflow discharges (6–36 L min−1m−1) 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−1m−1), 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−1s−1. 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.

ACKNOWLEDGEMENTS

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.

DATA AVAILABILITY STATEMENT

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

REFERENCES

REFERENCES
Auerswald
K.
Fiener
P.
Dikau
R.
2009
Rates of sheet and rill erosion in Germany – a meta-analysis
.
Geomorphology
111
(
3
),
182
193
.
Bagnold
R. A.
1966
An Approach to the Sediment Transport Problem From General Physics
.
US Geological Survey Professional Paper 442-I
.
USGS
,
Reston, VA
,
USA
.
Bennett
S. J.
Gordon
L. M.
Neroni
V.
Wells
R. R.
2015
Emergence, persistence, and organization of rill network
.
Natural Hazards
79
,
S7
S24
.
doi: 10.1007/s11069-015-1599-8
.
Berger
C.
Schulze
M.
Rieke-Zapp
D.
Schlunegger
F.
2010
Rill development and soil erosion: a laboratory study of slope and rainfall intensity
.
Earth Surface Processes and Landforms
35
(
12
),
1456
1467
.
Cerdan
O.
Le Bissonnais
Y.
Couturier
A.
Bourennane
H.
Souchère
V.
2002
Rill erosion on cultivated hillslopes during two extreme rainfall events in Normandy, France
.
Soil and Tillage Research
67
(
1
),
99
108
.
Consuelo
C. R.
Stroosnijder
L.
Guillermo
A.
2007
Sheet and rill erodibility in the northern Andean Highlands
.
Catena
70
(
2
),
105
113
.
Di Stefano
C.
Ferro
V.
Palmeri
V.
Pampalone
V.
2017
Flow resistance equation for rills
.
Hydrological Processes
31
,
2793
2801
.
Di Stefano
C.
Ferro
V.
Palmeri
V.
Pampalone
V.
2018
Testing slope effect on flow resistance equation for mobile bed rills
.
Hydrological Processes
32
(
5
),
664
671
.
Foster
G. R.
Huggins
L. F.
Meyer
L. D.
1984
A laboratory study of rill hydraulics: II: Shear stress relationships
.
Transactions of the ASAE
27
,
797
804
.
Giménez
R.
Govers
G.
2001
Interaction between bed roughness and flow hydraulics in eroding rills
.
Water Resources Research
37
(
3
),
791
799
.
Gimenez
R.
Govers
G.
2002
Flow detachment by concentrated flow on smooth and irregular beds
.
Soil Sci Soc Am J
66
(
5
),
1475
1483
.
Giménez
R.
Planchon
O.
Silvera
N.
Govers
G.
2004
Longitudinal velocity patterns and bed morphology interaction in a rill
.
Earth Surface Processes and Landforms
29
(
1
),
105
114
.
Govers
G.
Giménez
R.
Van
O. K.
2007
Rill erosion: exploring the relationship between experiments, modelling and field observations
.
Earth-Science Reviews
84
(
3–4
),
87
102
.
Han
D.
Wang
G.
Liu
T.
Xue
B.
Kuczera
G.
Xu
X.
2018
Hydroclimatic response of evapotranspiration partitioning to prolonged droughts in semiarid grassland
.
Journal of Hydrology
563
,
766
777
.
Heras
M. M.-D.
Espigares
T.
Merino-Martín
L.
Nicolau
J. M.
2011
Water-related ecological impacts of rill erosion processes in Mediterranean-dry reclaimed slopes
.
Catena
84
,
114
124
.
Horton
R. E.
Leach
H. R.
Vliet
R. V.
1934
Laminar sheet-flow
.
Transactions American Geophysical Union
15
(
2
),
393
404
.
Hu
C.
Wang
Y.
2004
Study on water-sediment optimum allocation in upstream basin and comprehensive measures of sediment control in Guanting reservoir. Reservoir sedimentation and general plan of water and sediment control
.
Journal of Sediment Research
2
,
19
26
(in Chinese with English abstract)
.
Knapen
A.
Poesen
J.
Govers
G.
Gyssels
G.
Nachtergaele
J.
2007
Resistance of soils to concentrated flow erosion: a review
.
Earth Science Reviews
80
(
1–2
),
75
109
.
Kou
P.
Xu
Q.
Yunus
A.
Dong
X.
Pu
C.
Zhang
X.
Jin
Z.
2020
Micro-topographic assessment of rill morphology highlights the shortcomings of current protective measures in loess landscapes
.
Science of the Total Environment
737
,
139721
.
Lei
T.
Nearing
M. A.
2000
Flume experiments for determining rill hydraulic characteristic erosion and rill patterns
.
Journal of Hydraulic Engineering
l31
(
11
),
49
54
(in Chinese with English abstract)
.
Lei
T.
Zhang
Q.
Zhao
J.
Tang
Z.
2001
A laboratory study of sediment transport capacity in the dynamic process of rill erosion
.
Transaction of American Society and Agriculture Engineers
44
,
1537
1542
(in Chinese with English abstract)
.
Li
G.
Abrahams
A. D.
Atkinson
J. F.
1996
Correction factors in the determination of mean velocity of overland flow
.
Earth Surface Processes and Landforms
21
(
6
),
509
515
.
Li
G.
Zheng
F.
Lu
J.
Xu
X.
Hu
W.
Han
Y.
2016
Inflow discharge impact on hillslope erosion processes and flow hydrodynamics
.
Soil Science Society of America Journal
80
(
3
),
711
719
.
Li
Q.
Wang
G.
Wang
H.
Shrestha
S.
Xue
B.
Sun
W.
Yu
J.
2020
Macrozoobenthos variations in shallow connected lakes under the influence of intense hydrologic pulse changes
.
Journal of Hydrology
584
,
124755
.
Liu
B.
Nearing
M.
Risse
L.
1994
Slope gradient effects on soil loss for steep slopes
.
Transactions of the ASAE
37
(
6
),
1835
1840
.
Liu
G.
Yang
M.
Warrington
D.
Liu
P.
Tian
J.
2011
Using beryllium-7 to monitor the relative proportions of sheet and rill erosion from loessal soil slopes in a single rainfall event
.
Earth Surface Processes and Landforms
36
(
4
),
439
448
.
Mahmoodabadi
M.
Ghadiri
H.
Yu
B.
Rose
C.
2014a
Morpho-dynamic quantification of flow-driven rill erosion parameters based on physical principles
.
Journal of Hydrology
514
,
328
336
.
Mahmoodabadi
M.
Ghadiri
H.
Rose
C.
Yu
B.
Rafahi
H.
Rouhipour
H.
2014b
Evaluation of GUEST and WEPP with a new approach for the determination of sediment transport capacity
.
Journal of Hydrology
513
,
413
421
.
Nearing
M.
1997
Hydraulics and erosion in eroding rills
.
Water Resources Research
33
(
4
),
865
876
.
Nearing
M.
Bradford
M.
Parker
C.
1991
Soil detachment by shallow flow at low slope
.
Soil Science Society of America Journal
55
(
2
),
339
344
.
Pijl
A.
Reuter
L. E. H.
Quarella
E.
Vogel
T. A.
Tarolli
P.
2020
GIS-based soil erosion modelling under various steep-slope vineyard practices
.
Catena
193
,
104604
.
Raff
D. A.
Ramírez
J. A.
Smith
J. L.
2004
Hillslope drainage development with time: a physical experiment
.
Geomorphology
62
,
169
180
.
Romero
C. C.
Stroosnijder
L.
Baigorria
G. A.
2007
Sheet and rill erodibility in the northern Andean Highlands
.
Catena
70
,
105
113
.
Shen
H. O.
Zheng
F. L.
Wen
L. L.
Lu
J.
Jiang
Y. L.
2015
An experimental study of rill erosion and morphology
.
Geomorphology
231
,
193
201
.
Shen
H. O.
Zheng
F. L.
Wen
L. L.
Han
Y.
Hu
W.
2016
Impacts of rainfall intensity and slope gradient on rill erosion processes at loessial hillslope
.
Soil & Tillage Research
155
,
429
436
.
Sirjani
E.
Mahmoodabadi
M.
2012
Effects of sheet flow rate and slope gradient on sediment load
.
Arabian Journal of Geosciences
7
,
203
210
.
Stefanovic
J. R.
Bryan
R. B.
2009
Flow energy and channel adjustments in rills developed in loamy and sandy loam soils
.
Earth Surface Processes and Landforms
34
,
133
144
.
Sun
L.
Fang
H.
Qi
D.
Li
J.
Cai
Q.
2013
A review on rill erosion process and its influencing factors
.
Chinese Geographical Science
23
(
4
),
389
402
.
Tian
P.
Xu
X.
Pan
C.
Hsu
K.
Yang
T.
2017
Impacts of rainfall and inflow on rill formation and erosion processes on steep hillslopes
.
Journal of Hydrology
548
(
5
),
24
39
.
Wagenbrenner
J. W.
Robichaud
P. R.
Elliot
W. J.
2010
Rill erosion in natural and disturbed forests: 2. Modeling approaches
.
Water Resources Research
46
,
1
12
.
Wang
Z.
Yang
X.
Liu
J.
Yuan
Y.
2015
Sediment transport capacity and its response to hydraulic parameters in experimental rill flow on steep slope
.
Journal of Soil and Water Conservation
70
(
1
),
36
44
.
Wang
D.
Wang
Z.
Shen
N.
Chen
H.
2016
Modeling soil detachment capacity by rill flow using hydraulic parameters
.
Journal of Hydrology
535
,
473
479
.
Wu
S.
Chen
L.
Wang
N.
Yu
M.
Assouline
S.
2018
Modeling rainfall-runoff and soil erosion processes on hillslopes with complex rill network planform
.
Water Resources Research
54
(
12
),
10117
10133
.
Yang
C. T.
1996
Sediment Transport: Theory and Practice
.
McGraw Hill
,
New York, USA
.
Zegeye
A. D.
Langendoen
E. J.
Steenhuis
T. S.
Wolde
M.
Tilahun Seifu
A.
2020
Bank stability and toe erosion model as a decision tool for gully bank stabilization in sub humid Ethiopian highlands
.
Ecohydrology and Hydrobiology
20
,
301
311
.
Zhang
G. H.
Liu
B. Y.
Nearing
M. A.
Huang
C. H.
Zhang
K. L.
2002
Soil detachment by shallow flow
.
Transactions of the ASAE
45
,
351
357
.
Zhang
P.
Yao
W.
Tang
H.
Wei
G.
Wang
L.
2017
Laboratory investigations of rill dynamics on soils of the loess plateau of China
.
Geomorphology
293
,
201
210
.