Applicability of Different Hydraulic Parameters to Describe Soil Detachment in Eroding Rills

This study presents the comparison of experimental results with assumptions used in numerical models. The aim of the field experiments is to test the linear relationship between different hydraulic parameters and soil detachment. For example correlations between shear stress, unit length shear force, stream power, unit stream power and effective stream power and the detachment rate does not reveal a single parameter which consistently displays the best correlation. More importantly, the best fit does not only vary from one experiment to another, but even between distinct measurement points. Different processes in rill erosion are responsible for the changing correlations. However, not all these procedures are considered in soil erosion models. Hence, hydraulic parameters alone are not sufficient to predict detachment rates. They predict the fluvial incising in the rill's bottom, but the main sediment sources are not considered sufficiently in its equations. The results of this study show that there is still a lack of understanding of the physical processes underlying soil erosion. Exerted forces, soil stability and its expression, the abstraction of the detachment and transport processes in shallow flowing water remain still subject of unclear description and dependence.

In most cases, a linear equation describes the relation between the hydraulic parameters mentioned above and the detachment rate. By exceeding a certain threshold, erosion by concentrated flow begins and detachment rate increases. This threshold has a positive x-axis intercept, which means that there is no detachment below this point.
Another option is to consider concentrated flow erosion as a nonlinear threshold phenomenon or as a two-part linear threshold phenomenon: below the threshold soil detachment takes place (first linear relationship) but after exceeding the threshold, detachment rate increases much faster (second linear relationship) [14]. But it is unclear if this linear relationship is really suitable.
Knapen et al. [14] calculated the correlation between shear stress, unit length shear force, stream power and Reynolds number and the detachment rate from several WEPP datasets. The best average correlation was determined for stream power with R 2 = 0.59. The WEPP-used shear stress is a variable that reaches only low R 2 values for all of the tested data sets. Knapen et al. [14] describes the shear stress as follows (p. 80 f.): ''Although the use of flow shear stress as soil detachment predictor can be contested, critical shear stress (t cr ) and concentrated flow erodibility KC (…) have been selected as the most universal parameters to describe soil erosion resistance to concentrated flow.'' The correlations between these factors and the soil detachment rate show very varying results. There is not a single parameter that always reveals the best correlation. These considerations lead to two main questions: 1. Are soil erosion, detachment and transport, directly dependent on water flow characteristics? 2. Are these concepts, as implemented in soil erosion models, suitable to describe rill erosion?
Due to difficulties to measure certain parameters, models have to be calibrated. During this process, the phenomenon of equifinality can appear: different parameter sets show the same result. Another weakness of rill erosion models is that the model parameters are often adapted from river hydrodynamics equa-tions. Govers and his colleagues [13,44] showed that these equations are not suitable for rill erosion processes. Therefore, there is often a mismatch between model results and observed or measured ''reality'' [43]. Additionally, models only project the concepts of the designer, not necessarily the reality.
In laboratory experiments, the initial and boundary conditions are well controlled. Soil parameters are well known and rill forms and slope can be adapted to the specific question. Thus, physical laws can be tested in a well-defined environment. However, Giménez and Govers [5] showed that parameters determined under laboratory conditions are not easily transformable to natural environments. One disadvantage of former laboratory experiments or field observations is the fact that in most cases only total runoff and sediment output are measured while the relative contribution of the individual processes is not considered [45].
Field data currently reflect the reality as close as possible. Nevertheless, observations as well as experiments show certain disadvantages: (1) Measurement techniques may disturb the observed processes, (2) time scale of human observations is shorter than that of the process under study, (3) some processes cannot be measured directly or indirectly and (4) some processes are chaotic and the spatial and temporal variations are difficult to specify [43].
The relationship between soil detachment and hydraulic parameters used in soil erosion models is in most cases deduced from laboratory experiments but the transferability of these results to natural rills is not generally given. Our setup in natural rills enables to measure the input parameters for calculating hydraulic parameters combining the advantages of laboratory experiments with the advantages of testing natural rills.
The main purpose of the field experiments was to quantify in a detailed temporal and spatial resolution the soil erosion dynamics in natural rills under concentrated flow for comparison of the measured sediment dynamics with those calculated by means of the most common detachment and transport equations.
Specifically, this study's objectives are: 1. elucidating the relationship between hydraulic parameters such as shear stress, unit length shear force, unit stream power, stream power, effective stream power and the Reynolds number and soil detachment in natural rills, 2. providing an explanation why physically-based soil erosion models do not capture rill erosion processes and 3. addressing the question whether current modelling approaches are generally suited to describe rill erosion processes.
The overall aim of this study is to have a critical view on concepts for modelling rill erosion based on experiments performed in naturally developed rills.

Ethics Statement
No specific permits were required for the described field studies. The mayors of the towns next to the study sites or the owners of the fields were informed about the intended activities and were asked for permission. The test sites Freila, Negratin and Salada are abandoned fields which are sporadically used as pasture for goats or sheep and in Belerda the experiment was accomplished on an almond field. The locations Freila, Negratin end Salada are not privately-owned and permission was granted from the owner of the study site Belerda. None of the study sites are protected in any way and the field studies did not involve endangered or protected species.

Study areas
The four study areas in Andalusia are located at Negratin, Freila, Salada and Belerda. UTM coordinates of the tested rills are given in Table 1.
Negratin and Freila. The areas are located within the Hoya de Baza sedimentary basin and composed of marls, in which calcareous Regosols have developed. The climate is semi-arid and vegetation is dominated by low shrubs and Stipa tenacissima grass tussocks. The land cover at the south side of the Negratin-dam is dominated by abandoned cereal fields, which are extensively grazed by sheep and agricultural land comprised mainly of cereal dry-farming and almond grooves [46].
Salada. Located at the SE-margin of the Betic range (SE-Spain), inside the penibetic complex. The area is composed of conglomerates with a clayey to loamy matrix, in which Regosols as well as to fairly developed (Calcic) Cambisols have developed. Vegetation is similar to that found in the Freila and Negratin-area. The climate is semi-arid too, but less accentuated than in the previously mentioned area [46]. Here the land use consists of rain fed agricultural areas (where cereals, olives, and almonds are cultivated), and abandoned or uncultivated areas.
Belerda. This test area is located in the Guadix basin. The parent material consists of tertiary and quaternary conglomerates, sands, silts and clays. The soil texture class following the FAO [47] is a silty clay loam. The land use is separated into cultivated areas, with almond and olive groves, and abandoned agricultural fields [48]. The climate is, though still semi-arid, characterised by higher average annual temperatures and precipitations in comparison with the other test zones.
The climatic parameters of the test fields are summarized in Table 1.

Tested rills
The main descriptors of the rills are summarized in Table 2. In this table, grain size class limits are from [49], texture class is determined following [47]. Photographies of the rills are presented in Figure 1.
The tested rills in Freila have developed on a sandy loam with high gravel content. Sand content is 57% with a relatively homogeneous contribution between coarse, medium, fine and very fine sand. The same is true in the silt fraction, the 34% are homogeneously contributed in the complete silt fraction between 63 and 2 mm. The rills show all a dense rock fragment cover and the highest vegetation cover of the four test sides.
In Negratin, the soil material is nearly gravel free, coarse, medium and fine sand also show low amounts, most of the fine material is in the grain size class ,20 mm. The rock fragment cover in the rill is higher than the gravel content of the soil material thus it is possible that residual rock fragment accumulation has occurred.
In Salada the grain size distribution is similar to Negratin. The highest account of the fine soil material is in the class ,63 mm. The residual rock fragment accumulation is formed even more clearly as in Negratin; the vegetation cover is relatively high compared to the other test sites.
The rill in Salada is the only rill that has developed in a field being used for agriculture. The soil material is composed by a mixture of all particle size classes from gravel to clay. The rock fragment cover is high compared to the other test sites and the vegetation cover comparatively low. This test site shows the highest dry bulk density which can be declared by the actual agricultural use.

Rill experiment (RE)
The rill experiments consist of two runs: first the rill is tested under field conditions (run a); in a second run (run b), approximately 15 minutes later, the same rill is tested under almost saturated soil conditions. A constant discharge of 250 L (or 330 L, respectively) is maintained during 4 minutes (or 3 minutes, respectively), using a motor-driven pump, resulting in a total water inflow of 1000 L. Mobilisation of material at the inflow has been avoided.
The flow velocity within the rill is characterized by the travel time of the waterfront and of two colour tracers (started at 1 and 2 minutes of the experiment), measured for every meter using a chronograph. By means of this procedure, three velocity curves are recorded and changes in flow dynamics can be detected. As colour tracers, food colourings (E 124 (red) and E 13 (blue)) are used for reasons of safety.
The rill's slope is characterized by measuring with a spring bow of 1 m range and a digital spirit level. It must be considered that slope measuring provides only average slopes for 1 meter. A step or a knick-point in the rill is not accounted, but its position and height are recorded.
Four water samples are taken at three different measuring points (MP1-MP3). The first sample is taken as soon as the waterfront has reaching the sampling point, the second 30 seconds later, the third 90 second later, and the fourth 150 seconds later.
The (suspended) sediment concentration SSC is determined by filtration of the samples in laboratory [50].
At each measuring point, rill cross section is measured. With a laser rangefinder, the distance between sensor and rill bottom is measured in 0.002 m steps. This allows an accurate calculation of the rills cross section area and an estimation of the rills volume.  Table 2. doi:10.1371/journal.pone.0064861.g001 Water level is continuously measured by ultrasonic sensors at each measuring point.

Descriptors for soil detachment
Soil detachment can be described by shear stress t, unit length shear force C, stream power v, unit stream power v U and effective stream power v eff . Reynolds number describes the balance between the inertial flow forces represented by the product in the numerator and the viscous forces as described by the dynamic viscosity in the denominator. It is a criterion for stability of a flowing medium. When Reynolds number is small, viscous forces dominate the motion and inertial ones can be ignored whereas at high Reynolds numbers inertial forces dominate and it is often possible to ignore viscosity [51]. Reynolds Number Re is calculated as follows: Liquid density is calculated using sediment concentration and grain density. The use of water's density is not practicable due to sediment concentrations of more than 400 g L 21 . Grain density was measured by a capillary pycnometer following DIN Table 2. Rill parameters: Grain size class limits are from [49], texture class is determined following [47].  Land use rangeland rangeland rangeland rangeland rangeland cropland K t is a transport coefficient, which has been adopted from the WEPP dataset. The WEPP-location is given. Measured values are starting soil moisture, maximum width, maximum depth, grain density, dry bulk density, org. material; parameters estimated in the field are vegetation cover and rock fragment cover; critical shear stress is calculated following WEPP. doi:10.1371/journal.pone.0064861.t002 between three measured velocities (arrival of the waterfront and arrival of the two colour tracers). Hydraulic radius and wetted cross section area can be calculated by measuring water level and the rill profile.
The viscosity of the sediment suspensions was measured with a shear rate controlled rheometer (Haake MARS from Thermo Fisher Scientific, Karlsruhe, Germany) and a cone-plate geometry with an angle of 2u and a diameter of 60 mm [53]. The shear rate c is defined as: with v = fluid velocity and y = the gap between the cone and base plate. The rheomter controls the shear rate and measures the shear stress t, from which the viscosity g is calculated via The sample volume is always 2.0 ml and the cell is tempered to 20uC+/20.01uC. Data points are taken at shear rates between 150 s 21 and 1500 s 21 . The viscosity does not depend on the shear rate. This is according to theoretical considerations. For a suspension of monodisperse particles one expects a linear relation [54,55] For quantification of the different processes in the rill, the transport rate T R [kg s 21 ] and the transport capacity T C [kg s 21 ] are calculated: Kt [s 2 m 0.5 kg 20.5 ] is a transport coefficient depending on soil substrate. The Kt value of the WEPP substrate which was most similar to the given test site conditions is used.

Quantification of different erosion processes
Following shear stress based model concepts, the transport rate cannot exceed the transport capacity [56]. Shear stress of the flowing water controls also the detachment. Therefore the transport rate up to the transport capacity is considered here as shear stress dependent uptake. The transport rate exceeding the transport capacity is considered as shear stress independent erosion caused by processes such as bank failure and headcut retreat. The resulting quantities are set into relation and given in percent of total transport rate.

Dynamic viscosity
The dynamic viscosity of the liquid shows a clear positive correlation with sediment concentration, i.e. dynamic viscosity increases with sediment concentration (see Figure 2). However, clear deviations from the trend line were observed for samples with low sediment concentrations, which were often rich in transported organic material. The small branchlets with low weight imply a low sediment concentration, but in rheometer measurements, they tilt and a high shear stress is erroneously measured. The trend line equation has been calculated for samples from different test sites, the R 2 -value of 0.92 indicates that this equation can be used for further experiments.

Correlations between detachment rate and hydraulic parameters
The R 2 values of the correlations between the detachment rate and different hydraulic parameters show the complete possible range from R 2 = 0 up to R 2 = 0.99 (see Table S19). Trend lines are increasing, decreasing and almost constant and thus it is not possible to find any clear dependency. Notably, only 40 of 252 correlations (about 16%) show an increasing trend line with an R 2 value$0.7. Table 4 shows that the highest average R 2 -value is calculated for the (t-t cr ) -detachment rate -relationship if all  Figure 3 shows the relationships between the measured transport rates and the predicted transport capacities. From 144 samples, in 82 cases the transport rate exceeds the capacity, corresponding to approximately 57% of all cases. Tables S20 and S21 present the differences between transport rates and transport capacities (S20) and the percentage of transport rate exceeding the capacity (S21) and hence the percentage of processes which are not controlled by the influence of shear stress. The percentage of material which is transported by processes independent of shear stress is on average 41.5% (see Table 5). Remarkably, the distribution is uneven, i.e. in the three Freila-experiments, the mean is 24.3% while in Negratin, Salada and Belerda, the average value is as high as 58.7% (see Table 5). The second group shows clearly higher sediment concentrations, meaning that the processes independent of shear stress provide higher sediment concentrations than the shear stress-based processes. This indicates that the influence of hydraulic parameters is higher for low sediment concentrations, or, in other words that high sediment concentrations are not caused by hydraulic parameters.

Discussion
A comparison with results of other research groups shows that the measured values are in a realistic range. Ghebreiyessus [3]     Therefore, these outcomes indicate that linear models may generally not be sufficient in order to describe the complex processes in natural rills. Four possible improvements may help to improve this important concept which has been studied already for over thirty years, (1) including a clear description of the employed parameters, (2) including the turbulence, (3) considering the impact of processes that do not depend on the shear stress and likewise (4) consider the high spatial and temporal variability observed in natural rills. These potential improvements will be discussed in more detail below.
For instance, the flow shear stress, a hydraulic parameter, and the critical shear stress, a soil parameter (similar to soil strength), must be differentiated. In particular, the flow shear stress must exceed the critical shear stress for erosion to occur. A number of hydraulic parameters, such as the flow velocity or the fluid density, water depth or width and roughness are used for the computation of the flow shear stress. The actual version of the shear stress equation calculates the average shear stress by depth averaging of momentum equation for steady uniform flow per area and time. Some factors used in shear stress calculation have been developed from empirical studies [15][16][17][18][19][20][21][22][23][24][25][26]. In most cases, the theoretical basis of the equations is however not clear. The formula applied by Chisci et al. [27] is derived from Landau and Lifchitz [57]. Other versions of the Landau-Lifchitz equation can be found in the literature [2][3][4][5].
The critical shear stress is the force needed to detach a soil particle. So it corresponds to a soil parameter and therefore, input  for calculations should also depend on soil characteristics. However, this is the case for the WEPP model [34] only, where the critical shear stress is calculated using soil parameters such as texture, organic matter content and dry bulk density. In other cases, both hydraulic and soil parameters are used [34]. The discrepancies in the methods of computation of the shear stress may be due to the conditions under which the equations are deduced, as these equations are based on empirical observations. The empirical nature of the development of the different expressions is clearly highlighted in previous work [30][31][32]. That means the equations are not deduced from physical laws but from empirical studies. In many studies [12,35,[37][38][39][40][41], neither critical shear stress nor shear stress are used for the calculation of the transport capacity at all. In other studies shear stress is used to calculate transport capacity and detachment capacity [36] or transport rate [42] and critical shear stress to calculate the detachment capacity [36]. In both cases it is clear that shear stress and critical shear stress operate against each other, the important parameter is the difference between these two variables.
A summary of these equations can be found in Reid and Dunne [58], on the EPA-homepage [59] and in Hessel and Jetten [60].
The second reason for the low R 2 -values in the correlations between hydraulic parameter and soil detachment can be the lack of turbulence parameters in the equations.
In the study of Knapen et al. [14] the Reynolds number shows very different correlations to the detachment rate, and this holds as well for the results of this study. The reason could be that the turbulence, described by the Reynolds number, does not directly operate on substrate, it influences the acting shear stress, that means the calculated shear stress is much lower than the operating shear stress, a relation which has been confirmed in several studies: Nearing et al. [61] found that turbulence can increase the active shear stress by a factor of several thousands. They measured flow shear stresses ranging from 0.5 to 2 Pa, and tensile strengths ranging from 1 to 2 kPa. Despite the fact that the tensile strengths are 1000 times larger than the flow shear stresses, the authors also measured detachment rates in the order of 300 g m 22 s 21 . Such large detachment rates were attributed to turbulent burst events. Another study about the influence of turbulence on detachment rates was published by Nearing & Parker [62]. They showed that under turbulent flow conditions the same shear stress value caused a clearly higher detachment rate. In their flume experiments the difference between detachment rate caused by turbulent and laminar flow increased with increasing shear stress value, i.e., if given hydraulic conditions lead to a high shear stress value, the influence of turbulence on soil erosion is higher than in low shear stress value ranges.
The shear stress equation, as well as the equations describing other hydraulic parameters, assumes that drag forces are dominant for controlling erosion. But rill erosion is the result of the combination of different processes including headcut erosion, sidewall sloughing, tunnelling, micro-piping, slaking piping and sapping [14,45,[63][64][65][66][67]. This is the third possible improvement for the problems of the model equations. The percentage of headcutting in the different studies ranges between ''four times higher than the contribution of bed scours'' [67] to ''60% of total rill erosion'' [68]. Stefanovic and Bryan [69] showed that concentrated flow causes sediment production primarily from knickpoints, chutes, meanders and bank failure. Govers [45] distinguished between hydraulic erosion, mass wasting processes on rill sidewalls, gullying and piping. Hydraulic rill erosion mostly occurred during three extreme runoff events. Mass wasting processes caused 37% of total erosion in rills. Gullying, the retreat erosion at knickpoints and headcuts caused about 12% of rill erosion rates. In the experiments presented here, the main mechanisms causing rill erosion were mass wasting and gullying processes, hence the correlations between hydraulic parameters and detachment rate are generally low. However, the hydraulic rill erosion only occurs in extreme runoff events, in most cases, the runoff values are too low to cause hydraulic rill erosion. The percentage of material which is transported independent of shear stress is very high on the water front samples. Here the transport of loose material is probably more important than in the other samples meaning that this process is mainly independent of shear stress. In these cases of transport rate vs. transport capacity ,1 the independence of shear stress cannot be excluded, in the other cases the processes controlled by shear stress can occur. Thus, it can be deduced that, in the case of T R .T C , not only shear stress controlled processes provide the material; at least the difference between T R and T C is caused by processes independent of shear stress.
The experiments presented here show that the correlation between hydraulic parameters and detachment rate does neither change from one experiment to another, nor from one run to another, but from one measuring point and run to another. Thus, sediment producing processes have a high spatial and temporal variability. This is the fourth possible improvement for models. It is very difficult to propose a single factor that always describes the soil detachment satisfactory. The high variability of erosion processes, even under controlled experimental conditions, has been highlighted in different studies. Measured variability shows a wide range between 3.4% and 173.2% [70][71][72][73][74][75]. This is partially the result of non-homogeneous parameters concerning soil characteristics and rainfall. On experimental plots, infiltration rates and soil aggregate stability can be highly variable [76] and rainfall also shows a high spatial and temporal variability [77]. Therefore, the input parameters to the different measurements reflected in the mentioned studies were not really comparable. Nevertheless, the results also make clear that modelling soil erosion has to include uncertainty in model input, as well as in the data used for model calibration and validation.
In field experiments, the spatial and temporal variability of soil conditions cannot be avoided, and is, furthermore, part of the investigations. Thus, additional input parameters as rainfall or flow should be maintained constant in the experiments to generate reproducible data. The high variability in soil erosion processes cannot be represented by a single factor like shear stress. The results show that there is not a simple linear correlation between a certain hydraulic parameter and soil detachment rate. Depending on model purpose and scale, the factors can be used to predict the magnitude of rill detachment but they are not applicable for the simulation of rill erosion with high-resolution spatial and temporal change in processes.
A newer approach is the use of probability density functions to predict soil detachment [78,79]. Sidorchuk gives two sources of stochasticity in erosion modelling: (1) the necessity of spatial and temporal averaging when determining deterministic equations, which describe concentrated flow erosion and (2) the fact that the main erosion factors, if these can be determined anyway, can only be measured with limited accuracy. This is not the first attempt to model erosion by relating the probability of soil detachment with the excess of erosion driving forces over soil erosion resistance forces, other articles using a stochastic approach to describe soil erosion were published by Nearing [80], Wilson [81] and Sidorchuk [82][83][84][85][86][87]. Notably, one of the earliest articles about stochastic in erosion processes has been published by Einstein [88]. These stochastic models reduce the number of empirical components. Applying these models to the experiments presented here is beyond the scope of the current study.

Conclusions
The results show that a linear correlation between hydraulic parameter and soil detachment is not sufficient to describe processes in natural rills. The reason for this behaviour is the combination of various processes that can cause different amounts of soil erosion. The shear stress, for instance, only describes one process, while the results clearly show that there is not one fixed parameter that always predicts soil detachment best. Applicability of one certain hydraulic parameter to predict the sediment concentration changes at a certain point in time within a few minutes, because the temporal and spatial distribution of the different erosion processes is highly randomly determined. Therefore, it might be more useful to formulate results in probabilistic terms, an approach which has already been implemented by previous researchers, but is beyond the current work.