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

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


INTRODUCTION
Soil erosion is a critical concern for the sustainable development of agricultural regions worldwide (Bennett et al. ), and soil erosion could cause non-point source pollution (Wang et al. a, b), soil organic carbon loss (Fang et al. ), and biodiversity reduction (Li et al. ).
Thus, soil erosion is a serious threat to the environment and agricultural security (Jiang et al. ). Rill erosion is one of the main patterns of soil erosion by water on sloping croplands and rangelands (Bewket &  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. ; Zhang et al. ). Giménez et al. () 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 ; Ran et al. ). In addition, many researchers found an increase in the slope gradient can lead to an increase in rill erosion amount (Berger et al. ; Fang et al. ). Ziadat & Taimeh () 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. ). However, Liu et al. () observed that the relationship between slope gradient and soil erosion varied when slope gradient reached certain thresholds. For instance, Yair & Klein () found the amount of soil erosion decreased with the increase of slope gradients. Koulouri & Giourga () 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 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.  (Mirzaee & Ghorbani-Dashtaki ). For instance, Giménez & Govers () 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 ; Mirzaee & Ghorbani-Dashtaki ). Furthermore, the rill bed form varies with erosion processes, which in turn, changes the hydraulics and erosion capacity of rill flow (Zhang et al.

; Jiang et al. ).
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.

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

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. (), six values of the inflow discharge per unit width (1 m) of the plot ranging from 6 to 36 L min À1 m À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. ). 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 ). 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 (v s ) were measured by the dye (KMnO 4 ) 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 highdefinition 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: (1) As an example, 1.29 ± 0.07 refers to the mean with standard deviations, the same in other cells. As an example, 1.35 ± 0.50 refers to the mean with standard deviations, the same in other cells.
Figure 2 | Field runoff plot system at the study site.
The surface runoff rate at the plot outlet (q out ) is calculated as follows: where q out is surface runoff rate at the plot outlet, L min À1 m À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. ): where q in is upslope inflow discharge (L min À1 m À1 ).
Mean rill flow velocity (v r ) is estimated by the measured surface rill flow velocity (v s ) multiplied by a correlation factor α as the equation: where a is 0.67 for laminar flow, 0.7 for transitional flow, Reynolds number (Re) and Froude number (Fr) can be used to reflect the rill flow regime, which are calculated as: where h r is the mean rill flow depth (m); ν m is kinematical viscosity (m 2 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: where g is gravitational acceleration (m s À2 ); J is surface slope (m m À1 ), which is equal to S in value; h r is the depth of rill flow (mm).
Shear stress (Foster et al. ), stream power (Bagnold ), and unit stream power (Yang ) are used to reflect the hydrodynamic parameters of rill flows, which are calculated as: where τ is flow shear stress (N m À2 ); W is flow stream power (N m À1 s À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 RD value is suitable for roughly describing the morphological characteristics of rills developed at hillslopes. The formulas of TRL and RD are: where TRL is the total lengths of all rills on the studied hillslope, m; RD is the rill density, m m À2 ; TRL i is the total length of a rill and its bifurcations (m); and i ¼ 1, …, n represents the number of the rills; A 0 is the plane area of the studied hillslope (m 2 ).
In addition, for a certain cross section of a rill, the rill width-depth ratio (R WD ) can reflect the shape of the rill cross section. For a certain cross section of a rill, R WD is calculated as follows: where MRW is the mean rill width, cm; MRD is the mean rill depth, cm; R WD is the rill width-depth ratio.

RESULTS AND DISCUSSION
Temporal changes in surface runoff and soil loss 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  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 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: SLR ¼ 0:0199q 1:239 S 1:075 (R 2 ¼ 0:954, p < 0:01, n ¼ 54) RER ¼ 0:0004q 1:868 S 1:524 (R 2 ¼ 0:977, p < 0:01, n ¼ 54) 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 À1 m À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
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)  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  L min À1 m À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.
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.
(), 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, 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 concen- 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 (), who conducted flume experiments under upslope inflow conditions. RD: rill density, (m m À2 ); MRW: mean rill width (cm); MRD: mean rill depth (cm); R WD : the rill width-depth ratio; S: slope gradient (%); q: inflow discharge (L min À1 m À1 ).

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.
Overall, mean rill flow velocity (v r ) increased with increases in either slope gradient or inflow discharge ( to 57%, respectively. For all slope gradients, the Re value began to exceed 2,000 when the inflow discharge increased to 24 L min À1 m À1 , which implied that rill flows became tur- Rill erosion rates significantly increased with an increase of τ, W, and φ. The relationships between rill erosion rates and where RER is rill erosion rate (g m À2 min À1 ), τ is shear stress (N m À2 ),  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.

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