- Email: [email protected]
Effects of transport capacity and erodibility on rill erosion processes: A model study using the Finite Element method
Geoderma 146 (2008) 114–120
Contents lists available at ScienceDirect
Geoderma j o u r n a l h o m e p a g e : w w w. e l s ev i e r. c o m / l o c a t e / g e o d e r m a
Effects of transport capacity and erodibility on rill erosion processes: A model study using the Finite Element method L.J. Yan a, X.X. Yu b, T.W. Lei a,d,⁎, Q.W. Zhang c, L.Q. Qu a a
Key Laboratory of Modern Precision Agriculture Integration Research, College of Hydraulic and Civil Engineering, China Agricultural University, Beijing 100083, PR China Key Laboratory of Soil and Water Conservation and Desertification Combat (Beijing Forest University), Ministry of Education, Beijing, 100083, PR China Institute of Environment and Sustainable Development in Agriculture, CAAS, Beijing, 100081, PR China d State Key Laboratory of Soil Erosion and Dryland Farming on the Loess Plateau, Institute of Soil and Water Conservation, CAS & MWR, Yangling, Shanxi, 712100, PR China b c
A R T I C L E
I N F O
A B S T R A C T
Article history: Received 16 March 2007 Received in revised form 6 May 2008 Accepted 18 May 2008 Available online 24 June 2008
Soil erosion prediction models are of great significance for soil and water conservation management. Rill erosion is the most important component of hillslope soil erosion processes. Therefore, predicting hillslope erosion requires that rill erosion is well understood and predictable. In this study, a physical process-based model for rill erosion was developed. A series of mathematical models were advanced to simulate rill hydrodynamics, soil detachment, and transport capacity for well-defined rill channels. The Finite Element method and a Visual C++ program were used to numerically solve the hydrodynamic and sediment continuity model equations. Model parameters, estimated by a different method, were used in simulation case studies to determine the impacts of their values on rill sedimentation processes, as compared with the experimental data. The comparison of experimental results with simulation outputs of the model validated the model equations and the computer program. The model is capable of simulating the dynamics of rill erosion processes. Comparisons of the simulated results with different model parameters indicated that: (i) the developed model can well simulate the rill erosion processes, provided that the model parameters were correctly estimated; (ii) the transport capacity limits the sedimentation process and the maximum possible sediment concentration at the rill outlet; (iii) rill erodibility determines the rate at which sediment concentration approaches its allowed value; and (iv) rill erodibility and transport capacity combined determine the distribution of soil erosion along a rill. Therefore, this study not only supplies a tool for further development of dynamic soil erosion prediction models but also verifies the importance of correctly estimating model parameters. © 2008 Elsevier B.V. All rights reserved.
Keywords: Finite Element method Simulation Rill erosion Dynamics Well-defined rill
Symbols list τ γ C Dc
shear stress of flow water specific gravity of water Chezy coefficient potential detachment
θ τc c DH
Dr
net detachment rate
f
g h
gravitational acceleration depth of water flow measured in the vertical direction the initial water depth along the rill rill erodibility parameter suggested by Lei rill length unit flow discharge
G H0
slope in degrees critical shear stress of soil sediment concentration hydrodynamic diffusion coefficient Darcy-Weisbach hydraulic frictional coefficient sediment load the water depth at the inlet
Kr Krw Q S
rill erodibility parameter rill erodibility parameter by WEPP flow rate slope grade
Hi Krl L q
⁎ Corresponding author. China Agricultural University, Qinghua Donglu, Beijing, 100083, PR China. Fax: +86 10 6273 6367. E-mail address: [email protected] (T.W. Lei). 0016-7061/$ – see front matter © 2008 Elsevier B.V. All rights reserved. doi:10.1016/j.geoderma.2008.05.009
Sfx frictional slope Sx rill slope steepness u depth- and width-averaged velocity in the x direction Vi the initial velocity along the rill x Cartesian coordinate ρ water density ω stream power function
Sss Tc V0
sediment source/sink term sediment transport capacity the velocity at the inlet
w β σ
rill width empirical depositional coefficient rainfall intensity
1. Introduction Rill erosion is the detachment by scouring and transport of sediments by a concentrated flow of water in a narrow, erodible channel (Bagnold, 1966). It is the main source of sediments and mechanism for sediment transportation in hillslope erosion processes. On the Loess Plateau of China, for example, rill erosion accounts for approximately 70% of the total erosion from upland areas (Zheng and Tang, 1997). Since the process of rill erosion is very complex, it is difficult to quantitatively identify how runoff and sediment load affect erosion under different rainfall, runoff, and surface conditions. It has been a
L.J. Yan et al. / Geoderma 146 (2008) 114–120
great challenge for erosion scientists to quantify processes contributing to rill erosion prediction. In recent years, process-based erosion prediction models have received increasing attention for various theoretical and practical reasons. Models such as WEPP (Water Erosion Prediction Project) (Nearing et al., 1989), GUEST (Hairsine and Rose 1992), and EuroSEM (Morgan et al., 1992; Morgan, 1994) as well as several others show great potential for application. These models have the capability to describe temporal and spatial distributions of soil loss. Furthermore, since the model is process-based it can be extrapolated to cover a wide range of situations (Flanagan et al., 1995). With the increase in computer capabilities, the Finite Element technique has been successfully applied to solve partial differential equations. The advantages of the method are its flexibility, generality, and consistency when compared to other numerical schemes (Sharda and Nearing, 1999). Some scientists have applied it in simulating runoff and soil erosion processes. For example, Sharda and Singh's Finite Element model simulated runoff and soil erosion from agricultural lands reasonably well (Sharda et al., 1994). Many significant studies in modeling and numerical schemes have been developed. For example, Lei et al. (1998) developed a model which simulated the physics of rill evolution. The model used a full form of the hydrodynamic routing equation, a probabilistically based soil detachment model, and a new sediment transport equation for rill erosion. Finite Element techniques were used to solve the erosion equations in this model. Models, such as WEPP, include transport capacity and rill erodibility parameters that are of great importance in rill erosion simulation. However, the means by which to quantify the impacts of such model parameters on erosion processes and to determine the model parameters themselves still remains a problem. Lei et al. (1998), in their simulation work, found that in order to achieve reasonably good modeling results the sediment transport capacity had to be calibrated before simulation. With the existing methodology, even though a model may accurately predict the sediment at rill/runoff plot outlets, very poor prediction of the rill development and sedimentation processes along a rill is a common problem. The present study had four main objectives: 1) to develop a numerical procedure based on the Finite Element method to simulate rill erosion processes; 2) to compare the simulated data with experimental results to validate the model and the numeric procedures, as well as the model parameters; 3) to test the influence of changes in model parameters, such as erodibility and transport capacity, on rill erosion processes; and 4) to compare experimentally observed rill sedimentation processes with simulated results that adopt parameters estimated by different methods to illustrate their influences on rill erosion processes and outlet sediments. 2. Mathematical models for rill erosion
115
depth- and width-averaged velocity in the x direction, σ(x, t) (m s− 1) is the rainfall intensity; g (=9.82 m s− 2) is gravitational acceleration; Sx (m m− 1) is rill slope steepness; it is equal to tan(θ),θ (°) is slope in degrees; and Sfx (m m− 1) is frictional slope given by (Chow, 1959): Sfx ¼
u2 1 C
ð3Þ
h
where C (m0.5 s− 1) is the Chezy coefficient given by (Chow, 1959): sffiffiffiffiffiffi 8g C¼ f
ð4Þ
f is the Darcy-Weisbach hydraulic frictional coefficient. Sediment continuity is given in the following equation: h
Ac Ac A Ac þ hu ¼ hDH þ ðSss −cσ Þ At Ax Ax Ax
ð5Þ
where c(x, t) (kg m− 3) is sediment concentration; DH (m2 s− 1) is the hydrodynamic diffusion coefficient for the sediments; and Sss(x, t) (kg m− 2 s− 1)) is the sediment source/sink term that, in detachmenttransport processes, is given by: cq Sss ¼ Dc 1− Tc
ð6Þ
where Dc (kg s− 1 m− 2) is potential detachment; Tc (kg m− 1 s− 1) is sediment transport capacity; and q (m2 s− 1) is unit flow discharge, while in depositional processes Sss is computed by: Sss ¼ −βðcq−Tc Þ
ð7Þ
−1
where β (1 m ) is an empirical depositional coefficient. 2.2. Model parameters 2.2.1. Sediment transport capacity Sediment transport capacity is defined as the maximum sediment load that can be carried along by a given water flow (Mancilla and Chen, 2005). Many attempts have been made to determine the sediment transport capacity based on data from various sources. Using laboratory generated experimental data for rill erosion, Nearing et al. (1997) related erosion to stream power, advanced by Bagnold (1966) as: ω ¼ ρgSx q
ð8Þ −3
−3
where ω (kg s ) is the stream power function; ρ (kg m ) is water density. Based on this concept, Nearing et al. (1997) suggested an equation to relate sediment transport capacity with stream power as: Be½aþblogðωÞ 1 þ e½aþblogðωÞ
2.1. Dynamic equations for hydrodynamic and sedimentation processes
log10 ðTc Þ ¼ A þ
Mass balance of water flow under rainfall and concentrated flow in the rill has been described in many studies (Milne-Thomson, 1960; Freeze, 1978; Ziekiewica and Taylor, 1991; Shames, 1992). Based on the research of Lei et al. (1998), the mass and momentum conservation/ balance equations are given by the following models:
where a = 0.845, b = 0.412, A = − 34.47, B = 38.61 were empirically derived coefficients. A laboratory flume experimental method for measuring transport capacity was advanced by Lei et al. (2001) where the transport capacity was linearly related to the slope and the flow rate as:
Ah AðuhÞ þ ¼σ At Ax
ð1Þ
Tc ¼ a þ bθ þ cQ
! AðuhÞ A u2 h g A h2 Sx þ þ ¼ gh pffiffiffiffiffiffiffiffiffiffiffiffiffiffi þ Sfx At 2 Ax Ax 1 þ S2x
ð2Þ
where x (m) is a Cartesian coordinate; t (s) is time; h(x, t) (m) is depth of water flow measured in the vertical direction; u(x, t) (m s− 1) is
ð9Þ
ð10Þ
where a = −0.3109, b = 0.01718, c = 0.1203 were regression coefficients and Q (L min− 1) was flow rate. 2.2.2. Rill erodibility The WEPP model provided a new erosion prediction tool based on water erosion physical processes. Rill erodibility is a key parameter for
116
L.J. Yan et al. / Geoderma 146 (2008) 114–120
rill erosion prediction in the WEPP model, as shown in the following equation: cq Dr ¼ Kr ðτ−τ c Þ 1− Tc −1
−2
Dc ¼ ð11Þ −2
where Dr (kg s m ) is the net detachment rate; τ and τc (N m ) are the shear stress of flowing water and the critical shear stress of soil, respectively; Kr (kg N− 1 s− 1) is the rill erodibility parameter; cq (kg m− 1 s− 1) is the sediment yield in the rill. The effects of soil erodibility (Kr) play a great role in soil erosion dynamic processes. The higher the Kr value, the more readily the soil particles are detached and, hence, more soil material is entrained in the water flowing from upstream rill sections. Consequently, there is a faster increase in the amount of sediment load carried in the flowing water and the transport capacity is approached more rapidly. When the sediment load tends to reach the transport capacity, the erosion intensity is reduced downstream. Eq. (11) shows that the net detachment rate has its maximum value, the potential detachment rate, when water is clear, i.e., when cq = 0, which is given by: Drmax ¼ Dc ¼ Kr ðτ−τc Þ:
ð12Þ
Values for Kr and τc were estimated by linear regression of a set of detachment rates and the corresponding hydraulic shears (Foster and Meyer, 1972; Nearing et al., 1989; Ghebreiyessus et al., 1994). A rational method for determining the soil erodibility and critical shear stress of rill erosion under concentrated flow was advanced by Lei et al. (2008). Using their method, soil erodibility of a loess soil was determined to be 0.3211 ± 0.0005 kg N− 1 s− 1, and the critical shear stresses τc were found to increase with slope gradient, i.e., τc was 3.191 N m− 2, 3.960 N m− 2, 4.130 N m− 2, 4.385 N m− 2, 4.565 N m− 2 for slopes of 5°, 10°, 15°, 20°, and 25° respectively. Gilley et al. (1993) reported a method for rill erodibility determination for WEPP. In their experiments to determine the parameter, simulated rainstorms were applied to field runoff plots, 0.46 by 9 m, to generate rills. Subsequently, concentrated flows of different flow rates were introduced into the rills to simulate rill erosion. The sediment yields thus obtained from the runoff plots were used to estimate detachment rates, for the Dc values in Eq. (12), by: Dc ¼
cq L
ð13Þ
where L (m) is the rill length, L = 9 m. The estimated detachment rates and corresponding hydraulic shear stresses under different hydraulic conditions were then used to determine rill erodibility and critical shear stress by regressing the data using Eq. (12). Zhang et al. (2001) experimentally confirmed the sedimentation process in a well-defined rill to be described in the following form: c ¼ A 1−e−Bx
ð14Þ
where A (kg m− 3) is a regression coefficient, representing the highest possible sediment concentration in the rill, under a given hydraulic regime; B is a decay coefficient, (1 m− 1). This B in Eq. (14) indicates how fast the sediment concentration reaches the value at the transport capacity, c = A. As indicated in Eq. (17), B is the pure soil parameter to determine the potential detachment rate, while Tc determines the hydraulic condition for potential detachment rate. Lei et al. (2001) also verified the following relationship: Aq ¼ Tc
Combining Eqs. (14) and (13) yields:
ð15Þ
where q (m2 s− 1) is unit width flow rate of the given flow; A is same as Eq. (14).
Tc 1−e−BL cq Aq 1−e−BL ¼ ¼ L L L
ð16Þ
Lei et al. (2002) also suggested a method to estimate the potential detachment rate: dA 1−e−Bx q ¼ ABq ¼ BTc dx xY0
Dc ¼ lim Dr ¼ lim xY0
ð17Þ
Eqs. (16) and (17) give two different methods by which to estimate rill erodibility. One, used by WEPP, is based upon Eq. (16), and is given by: Krw ¼
Tc 1−e−BL Dr : ¼ τ−τc ðτ−τc ÞL
ð18Þ
The other, suggested by Lei et al. (2002) and based upon Eq. (17), takes the following form: Krl ¼
Dr Tc B ¼ : τ−τc τ−τ c
ð19Þ
In order to verify the difference between the two erodibility parameters estimated by these two methods, Eq. (18) can be divided by Eq. (19) to give: T B
c Krl BL c ¼ T ðτ−τ −BL Þ ¼ c 1−e 1−e−BL Krw
ð20Þ
ðτ−τc ÞL
Krw ¼
Krl 1−e−BL BL
ð21Þ
where Krl is rill erodibility suggested by Lei and Krw is rill erodibility used by WEPP. This functional relationship (Eq. (21)) suggests that the estimate for soil erodibility used by WEPP depends upon the value of B, which is a complex coefficient affected by factors such as water flow, rill slope, and soil characteristics. Different B values lead to different soil erodibilities and in turn affect the simulation of the erosion process. In order to quantitatively understand the influences of B and Krw values on rill erosion processes, different Krw values and those computed from various B values for a 25° rill slope are listed in Table 1 to be used in this study. Experimental data, from a study by Zhang (2002) on a silty-clay loess soil, were used to evaluate the impacts of model parameters on rill erosion processes. The experiments were conducted in a flume, 8 m long and 1 m wide, which was sub-divided into strips, 8 m by 0.1 m, to imitate rills. The experiments involved five slopes (5°, 10°, 15°, 20°, and 25°), 9 slope lengths (0.5, 1, 2, 3, 4, 5, 6, 7, 8 m), and three flow rates (8 and 12 L min− 1 for 5° slopes; 2, 4, and 8 L min− 1 for steeper slopes). Numerical simulations were made for the same conditions as those applied in the experiments for the model verification. Simulated sediment concentrations along the rills were also compared with those from the experiments. Table 1 Krl and Krw under different slope gradients; different Krw values computed from different B values at 25° Slope
Krl (kg N− 1 s− 1)
Krw (kg N− 1 s− 1)
5° (B = 0.3) 10° (B = 0.2) 15° (B = 0.3) 20° (B = 0.4) 25° (B = 0.3)
0.3209 0.3211 0.3214 0.3214 0.3216
0.1109 0.1489 0.1110 0.0868 0.1111 (B = 0.3) 0.0707 (B = 0.5) 0.0446 (B = 0.8) 0.0357 (B = 1.0)
L.J. Yan et al. / Geoderma 146 (2008) 114–120
Fig. 1. a–m) Comparisons of simulated rill erosion processes with experimental data. 117
118
L.J. Yan et al. / Geoderma 146 (2008) 114–120
2.3. Boundary and initial conditions With reference to the experimental setup and the simulation of the hydraulic and erosion processes in a well-defined rill, the boundary conditions are listed as follows: 8 Vx ¼ V0 > > > 0 > > > > > > Hx0 ¼ H0 > > > > > > > > dc > > ¼ Dc > > > dx x¼x0 > > > > > : τð0; t Þ ¼ γSH0
ð22Þ
Fig. 2. Comparison of sediment concentration for simulations and experiments.
where γ is specific gravity of water, 9800 (N m− 3); S is hydraulic slope, sin(θ); V0 (m s− 1), and H0 (m), the water velocity and depth at the inlet respectively, are interrelated and determined by the flow in-flow rate, q (m3 s− 1), as: q ¼ V0 H0 w
ð23Þ
where w (m) is rill width. The initial conditions for the governing equations are given as: 8 < Vxi ðx; t ¼ 0Þ ¼ Vi ð24Þ H ðx; t ¼ 0Þ ¼ Hi : i ci ðx; t ¼ 0Þ ¼ 0 where Vi (m s− 1) and Hi (m), the initial water velocity and depth values along the rill, can be assigned arbitrarily with any values, except negative and zero Hi values. 3. Results and discussion Using the numerical formulations discussed above, a program was written in Visual C++ to simulate the rill erosion processes using the Finite Element method. Various model parameter values, for the transport capacity, rill erodibility and critical shear stress were used in the simulation with the aim of comparing the impacts of these different parameters, as estimated by different methods, on rill erosion processes. 3.1. Comparison of simulated rill erosion processes with experimental data The simulated sediment concentrations along rills, under different slope and flow rate conditions with parameter sets including transport capacity, rill erodibility and critical shear stress estimated by Lei et al. (2001, 2008), were compared with the experimental data collected under the same conditions, as shown by the Rational Tc and Kr line in Fig. 1a–m. The sediment concentration hygrographs, for both the simulated and the experimental data, were in close agreement with the relationship between sediment concentration and rill length described by Lei et al. (1998) and Huang et al. (1996). Fig. 1a–m for all the combinations of slope gradients and flow rates indicate very good agreement between the simulated results and those from the experiments. Sediment concentration increases with rill length. In the upper part of the rill, i.e., when the down-slope distance is short (Fig. 1a to m), the greatest increases in sediment concentration were observed indicating higher detachment rates. The rate of increase in sediment concentration becomes lower with increasing rill length, and sediment concentration tends to approach a stable maximum value that is limited by the transport capacity of the water flow. Using all the data sets consisting of 117 data points from the 13 experimental runs, a comparison was made against simulated data, as
shown in Fig. 2. Linear regression indicated that the simulated sediment concentrations were closely correlated with the measured data, with a slope coefficient of 0.989 and a coefficient of determination of 0.90. Therefore, there exists a very good one to one relationship between the two sets of data. These results imply that the simulation model, including the mathematical formulations and the numerical algorithms of the Finite Element method, as well as the computational procedure codes, were correct. Moreover, the model parameters set, including the transport capacity, rill erodibility, and the critical shear stress estimated by Lei et al. (2001, 2008), must likewise be correct in order to produce such good predictions of sediment yields not only at the rill outlets but also at the various locations along the rills. 3.2. Impacts of transport capacity on rill erosion processes Simulations were also made, using the transport capacity given in Eq. (9) reported by Lei et al., (1998) and the rill erodibility and critical shear stress estimated by Lei et al. (2008), to check the impacts of this parameter on the rill sedimentation processes. The results are indicated by the Tc — Limited curves in Fig. 1a–m. Comparisons were made between the experimentally measured data and the simulated data sets. The comparisons show that generally the simulated sediment concentrations in the water flow in the upper parts of the rill (b1 m down-slope distance) were in agreement with the experimental data. The sediment concentration increases rapidly along this section of the rill. As sediment concentrations approach the limiting transport capacity of the water flow, the simulated sedimentation loads gradually and then rapidly deviate from the experimental data to considerably lower amounts. At the end of the rills, the simulated sediment yields were dramatically reduced as compared with the measured data. These simulated results indicated how the transport capacity parameter limited the sedimentation processes behavior. This might, to some extent, have explained the reason why calibration of the transport capacity was needed before a reasonably good prediction of sediment concentrations at rill outlets was obtained by the rill erosion model of Lei et al. (1998). These results indicate that the transport capacity parameter is of great importance in rill erosion process modeling. 3.3. Impacts of rill erodibility on rill erosion processes To check the impact of the erodibility parameter on rill sedimentation processes, simulations were also made with rill erodibility values estimated by the method used in the WEPP model, as given by Eq. (21), with a B value of 0.2, 0.3, and 0.4 for various slopes. The simulated results are depicted by the Kr — Influence curves in Fig. 1a– m. The results show that the simulated sediment concentration curves consistently underestimated rill erosion considerably for all rill sections. Therefore, rill erodibility estimated by the WEPP method could be too low to make a good prediction of rill sedimentation processes.
L.J. Yan et al. / Geoderma 146 (2008) 114–120
Fig. 3. Simulated results between different Krw values computed from different B values for a 25o slope.
Nevertheless, Fig. 1a–m also showed that, generally, the sediment concentration curves tended to increase gradually down the length of the rills, approaching the level of the potential sediment concentrations as determined by the transport capacity. Potentially, as rill length further increases, the sediment concentration may reach the levels defined by the transport capacity, provided the rills are long enough. These results imply that, although the rill erodibility value might be too low to estimate the sedimentation processes along a short rill, the sediment yields at rill outlets could be reasonably close to the actual values provided the rill is long enough and the transport capacity is correctly defined. The overall error of the sediment concentrations at the rill outlets between the simulated results and the experimental data for the rills of 8 m long was 30%. These simulation results indicate that, although the lowered values of estimated rill erodibility could give reasonably good predictions of sediment yields at the outlets of a long rill, if the transport capacity is correctly defined, a good estimation of sediment yields at rill outlets alone cannot guarantee that there was good simulation of erosion dynamics along the rills, unless comparisons are made along the length of the rill, as carried out in the previous section. Thus, simulated results estimated by WEPP might need reappraisal with sedimentation data collected along the rills instead of using the sediment yields at outlets alone. A comparison was also made between the experimentally measured rill sedimentation processes, at a slope gradient of 25° and a flow rate of 4 L min− 1, and the simulated results predicted when using rill erodibilities estimated by Eq. (21) for different B values (0.3, 0.5, 0.8 and 1), as used by the WEPP model, as well as the simulated result for rill erodibility estimated by Lei et al. Eq. (19). The transport capacity was given by Eq. (10). The results are shown in Fig. 3 and indicate that, with an increase in the B value, rill erodibility as determined by Eq. (21) decreases, and that the estimated sedimentation loads were considerably reduced. The sediment yield at the outlet of an 8 m long rill is reduced by 70% when B = 1 compared with the results produced using a rill erodibility value suggested by Lei et al. (2008). That is, the sediment yield at the rill outlet predicted using the rill erodibility value of Lei et al. (2008) was more than twice that calculated with the erodibility estimated by the WEPP model method. 3.4. Combined impacts of transport capacity and rill erodibility on erosion processes Further comparisons were made using various transport capacity and rill erodibility values to see the combined impacts of these two important parameters on rill erosion dynamics. The results are illustrated in Fig. 1a–m. The results indicate that the simulated sedi-
119
mentation processes are in good agreement with the experimental data when used with the correctly correctly-estimated model parameters of Lei et al. (2008). When the simulation is limited by transport capacity alone (i.e., with the rill erodibility correctly estimated), the sedimentation processes in the upper sections of the rill relate well to the experimental data. However, the simulation derived curves gradually diverge from the experimental data as down-slope distance increases. The predicted sediment loads finally reach the limit defined by the imposed transport capacity value. This would inevitably cause low estimations of the sediment yields at the rill ends as well as for the sediment detachment rates along the rills. In the cases where the rill erodibility value influenced the model, the simulated sedimentation processes differed from the experimental curves for all sections of the rills. However, the sediment load generally continued to increase with down-slope distance towards a level that corresponded to the transport capacity, which potentially could be attained if the rills were long enough. Thus, the estimated sediment yields at the rill outlets could be very close to the actual value for long rills. However, the estimated sediment loads from specific rill locations would differ from those values derived used a correctlyestimated rill erodibility value. 4. Conclusions Mathematic models were formulated based on hydrodynamics and physical sedimentation concepts, outlined by Nearing et al. (1997) and advanced by Lei et al. (1998b). Numerical simulations were made using the Finite Element method and a Visual C++ program to simulate rill erosion dynamics under different rill slope and hydraulic conditions. Comparisons made between the experimental data and the simulated results showed that the models, numerical algorithms, and the program all work well. Rill sedimentation processes were accurately simulated by the model (as given by Eq. (19)) and rational Kr and Tc. Sediment concentration data from the simulations were closely correlated with the experimental data and these very good predictions were made without calibration of the model parameters, which indicated that the model parameters estimated were correct and emphasized the importance of model parameter estimation. Reducing the transport capacity values produced transport capacity limited sedimentation processes and lowered sediment yield estimations below those measured. Simulation using lower rill erodibility values could make reasonably good predictions of sediment yields at the outlets of long rills given good estimated transport capacity values. However, the sediment sources might not be the same as in real situations. Checking the simulation results of both sediment yields at the rill outlets and the sedimentation processes along the rills is important. This study not only supplies a tool for further development of dynamic soil erosion prediction models but also verifies the importance of correct estimation of model parameters. Acknowledgments This paper is based on the work supported by the Key Laboratory of Soil and Water Conservation & Desertification Combat (Beijing Forest University), Ministry of Education (2005k-01-01), by the Natural Science Foundation of China (40635027), Program for Changjiang Scholars and Innovative Research Team in University, as well as supported by the CAS/SAFEA International Partnership Program for Creative Research Teams, 2007-380. References Bagnold, R.A., 1966. An approach to the sediment transport problem from general physics. US Geological Survey Paper 422-∣, Washington.
120
L.J. Yan et al. / Geoderma 146 (2008) 114–120
Chow, V.T., 1959. Open-channel Hydraulics. McGraw-Hill, New York, N Y. Flanagan, D.C., Ascough, J.C., Nicks, A.D., Nearing, M.A., Laflen, J.M., 1995. Overview of the WEPP erosion prediction model, technical documentation. National Soil Erosion Research Laboratory, West Lafayette.USDA-Water Erosion Prediction Project. Foster, G.R., Meyer, L.D., 1972. Transport of soil particles by shallow flow. Transactions of the ASAE 15, 99–102. Freeze, R.A., 1978. Mathematical models of hillslope hydrology. In: Kirkby, M.J. (Ed.), Hillslope Hydrology. Wiley Interscience, New York, pp. 177–225. Ghebreiyessus, Y.T., Gantzer, C.J., Alberts, E.E., Lentz, R.W., 1994. Soil erosion by concentrated flow: shear stress and bulk density. Transactions of the ASAE 37, 1791–1797. Gilley, J.E., Elliot, W.J., Laflen, J.M., Simanton, J.R., 1993. Critical shear stress and critical flow rates for initiation of rilling. Journal of Hydrology 142, 251–271. Hairsine, P.B., Rose, C.W., 1992. Modeling water erosion due to overland flow using physical principles 2: rill flow. Water Resources Research 28, 245–250. Huang, C.H., Laflen, J.M., Bradford, J.M., 1996. Evaluation of the detachment-transport coupling concept in the WEPP rill erosion equation. Soil Sci. Soc. Am. J. 60, 734–739. Lei, T.W., Nearing, M.A., Haghighi, K., Bralts, V.F., 1998. Rill erosion and morphological evolution: a simulation model. Water Resources Research 34, 3157–3168. Lei, T.W., Zhang, Q.W., Yan, L.J., Zhao, J., Pan, Y.H., 2008. A rational method for estimating erodibility and critical shear stress of an eroding rill. Geoderma 244 (3–4), 628–633. Lei, T.W., Zhang, Q.W., Zhao, J., Nearing, M.A., 2002. Analytic method for determination of detachment rate of concentrated flow in erosion rills. Acta Pedologica Sinica 39, 788–793 (in Chinese with English abstract). Lei, T.W., Zhang, Q.W., Zhao, J., Tang, Z.J., 2001. Laboratory study on sediment transport capacity in the dynamic process of rill erosion. Transactions of the ASAE 44, 1537–1542. Mancilla, G.A., Chen, S.L., 2005. A field study on sediment transport capacity on rills. ASAE Annual International Meeting Sponsored by ASAE. Tampa, Florida.Paper No. 052051. ASAE, St Joseph, Mich. Milne-Thomson, L.M., 1960. Theoretical Hydrodynamics. Macmillan, New York. Morgan, R.P.C., 1994. The European soil erosion model: an update on its structure and research base. In: Rickson, R.J. (Ed.), Conserving Soil Resources: European
Perspectives. Selected Papers from the First International Congress of the European Society for Soil Conservation. Wallingford. UK, pp. 286–299. Morgan, R.P.C., Quinton, J.N., Rickson, R.J., 1992. EUROSEM: Documentation Manual. Silsoe College, Silsoe, UK. Nearing, M.A., Foster, G.R., Lane, L.J., Finkner, S.C., 1989. A process based soil erosion model for USDA. Water Erosion Prediction Project Technology. Transaction of the ASAE 32, 1587–1593. Nearing, M.A., Norton, L.D., Bulgakov, D.A., Larionov, G.A., West, L.T., Dontsova, K.M., 1997. Hydraulics and erosion in eroding rills. Water Resources Research. Res 33, 865–876. Shames, I.H., 1992. Mechanics of Fluids, 3rd Edn. McGraw-Hill Book Company, New York. Sharda, V.N., Nearing, M.A., 1999. Finite element modeling of erosion on agricultural lands. In: Stott, D.E., Mohtar, R.H., Steinhardt, G.C. (Eds.), Sustaining the Global Farm. Selected Papers from the 10th International Soil Conservation Organization Meeting. Purdue University and the USDA-ARS National Soil Erosion Research Laboratory, West Lafayette, IN, pp. 1008–1017. Sharda, V.N., Singh, S.R., Sastry, G., Dhruvanarayana, V.V., 1994. A finite element model for simulating runoff and soil erosion from mechanically treated agricultural lands 2. Field validation and applications. Water Resources Research 30, 2299–2310. Zhang, Q.W., Lei, T.W., Pan, Y.H., Gao, P.L., 2002. Rational computational method of soil erodibility and critical shear stress from experiment data. Journal of the Graduate School of the Chinese Academy of Sciences 21, 468–475 (in Chinese with English abstract). Zhang, Q.W., Lei, T.W., Zhao, J., 2001. Coupling relationship between detachment rate with sediment load and rill length. Journal of soil and water conservation 15, 92–95 (in Chinese with English abstract). Zheng, F.L., Tang, K.L., 1997. Rill erosion process of steep slope land of the Loess Plateau. International Journal of S. Ziekiewica, C., Taylor, R.L., 1991. 4th Edn. The Finite Element Method, vol. 2. McGrawHill Book Company, New York.
Contents lists available at ScienceDirect
Geoderma j o u r n a l h o m e p a g e : w w w. e l s ev i e r. c o m / l o c a t e / g e o d e r m a
Effects of transport capacity and erodibility on rill erosion processes: A model study using the Finite Element method L.J. Yan a, X.X. Yu b, T.W. Lei a,d,⁎, Q.W. Zhang c, L.Q. Qu a a
Key Laboratory of Modern Precision Agriculture Integration Research, College of Hydraulic and Civil Engineering, China Agricultural University, Beijing 100083, PR China Key Laboratory of Soil and Water Conservation and Desertification Combat (Beijing Forest University), Ministry of Education, Beijing, 100083, PR China Institute of Environment and Sustainable Development in Agriculture, CAAS, Beijing, 100081, PR China d State Key Laboratory of Soil Erosion and Dryland Farming on the Loess Plateau, Institute of Soil and Water Conservation, CAS & MWR, Yangling, Shanxi, 712100, PR China b c
A R T I C L E
I N F O
A B S T R A C T
Article history: Received 16 March 2007 Received in revised form 6 May 2008 Accepted 18 May 2008 Available online 24 June 2008
Soil erosion prediction models are of great significance for soil and water conservation management. Rill erosion is the most important component of hillslope soil erosion processes. Therefore, predicting hillslope erosion requires that rill erosion is well understood and predictable. In this study, a physical process-based model for rill erosion was developed. A series of mathematical models were advanced to simulate rill hydrodynamics, soil detachment, and transport capacity for well-defined rill channels. The Finite Element method and a Visual C++ program were used to numerically solve the hydrodynamic and sediment continuity model equations. Model parameters, estimated by a different method, were used in simulation case studies to determine the impacts of their values on rill sedimentation processes, as compared with the experimental data. The comparison of experimental results with simulation outputs of the model validated the model equations and the computer program. The model is capable of simulating the dynamics of rill erosion processes. Comparisons of the simulated results with different model parameters indicated that: (i) the developed model can well simulate the rill erosion processes, provided that the model parameters were correctly estimated; (ii) the transport capacity limits the sedimentation process and the maximum possible sediment concentration at the rill outlet; (iii) rill erodibility determines the rate at which sediment concentration approaches its allowed value; and (iv) rill erodibility and transport capacity combined determine the distribution of soil erosion along a rill. Therefore, this study not only supplies a tool for further development of dynamic soil erosion prediction models but also verifies the importance of correctly estimating model parameters. © 2008 Elsevier B.V. All rights reserved.
Keywords: Finite Element method Simulation Rill erosion Dynamics Well-defined rill
Symbols list τ γ C Dc
shear stress of flow water specific gravity of water Chezy coefficient potential detachment
θ τc c DH
Dr
net detachment rate
f
g h
gravitational acceleration depth of water flow measured in the vertical direction the initial water depth along the rill rill erodibility parameter suggested by Lei rill length unit flow discharge
G H0
slope in degrees critical shear stress of soil sediment concentration hydrodynamic diffusion coefficient Darcy-Weisbach hydraulic frictional coefficient sediment load the water depth at the inlet
Kr Krw Q S
rill erodibility parameter rill erodibility parameter by WEPP flow rate slope grade
Hi Krl L q
⁎ Corresponding author. China Agricultural University, Qinghua Donglu, Beijing, 100083, PR China. Fax: +86 10 6273 6367. E-mail address: [email protected] (T.W. Lei). 0016-7061/$ – see front matter © 2008 Elsevier B.V. All rights reserved. doi:10.1016/j.geoderma.2008.05.009
Sfx frictional slope Sx rill slope steepness u depth- and width-averaged velocity in the x direction Vi the initial velocity along the rill x Cartesian coordinate ρ water density ω stream power function
Sss Tc V0
sediment source/sink term sediment transport capacity the velocity at the inlet
w β σ
rill width empirical depositional coefficient rainfall intensity
1. Introduction Rill erosion is the detachment by scouring and transport of sediments by a concentrated flow of water in a narrow, erodible channel (Bagnold, 1966). It is the main source of sediments and mechanism for sediment transportation in hillslope erosion processes. On the Loess Plateau of China, for example, rill erosion accounts for approximately 70% of the total erosion from upland areas (Zheng and Tang, 1997). Since the process of rill erosion is very complex, it is difficult to quantitatively identify how runoff and sediment load affect erosion under different rainfall, runoff, and surface conditions. It has been a
L.J. Yan et al. / Geoderma 146 (2008) 114–120
great challenge for erosion scientists to quantify processes contributing to rill erosion prediction. In recent years, process-based erosion prediction models have received increasing attention for various theoretical and practical reasons. Models such as WEPP (Water Erosion Prediction Project) (Nearing et al., 1989), GUEST (Hairsine and Rose 1992), and EuroSEM (Morgan et al., 1992; Morgan, 1994) as well as several others show great potential for application. These models have the capability to describe temporal and spatial distributions of soil loss. Furthermore, since the model is process-based it can be extrapolated to cover a wide range of situations (Flanagan et al., 1995). With the increase in computer capabilities, the Finite Element technique has been successfully applied to solve partial differential equations. The advantages of the method are its flexibility, generality, and consistency when compared to other numerical schemes (Sharda and Nearing, 1999). Some scientists have applied it in simulating runoff and soil erosion processes. For example, Sharda and Singh's Finite Element model simulated runoff and soil erosion from agricultural lands reasonably well (Sharda et al., 1994). Many significant studies in modeling and numerical schemes have been developed. For example, Lei et al. (1998) developed a model which simulated the physics of rill evolution. The model used a full form of the hydrodynamic routing equation, a probabilistically based soil detachment model, and a new sediment transport equation for rill erosion. Finite Element techniques were used to solve the erosion equations in this model. Models, such as WEPP, include transport capacity and rill erodibility parameters that are of great importance in rill erosion simulation. However, the means by which to quantify the impacts of such model parameters on erosion processes and to determine the model parameters themselves still remains a problem. Lei et al. (1998), in their simulation work, found that in order to achieve reasonably good modeling results the sediment transport capacity had to be calibrated before simulation. With the existing methodology, even though a model may accurately predict the sediment at rill/runoff plot outlets, very poor prediction of the rill development and sedimentation processes along a rill is a common problem. The present study had four main objectives: 1) to develop a numerical procedure based on the Finite Element method to simulate rill erosion processes; 2) to compare the simulated data with experimental results to validate the model and the numeric procedures, as well as the model parameters; 3) to test the influence of changes in model parameters, such as erodibility and transport capacity, on rill erosion processes; and 4) to compare experimentally observed rill sedimentation processes with simulated results that adopt parameters estimated by different methods to illustrate their influences on rill erosion processes and outlet sediments. 2. Mathematical models for rill erosion
115
depth- and width-averaged velocity in the x direction, σ(x, t) (m s− 1) is the rainfall intensity; g (=9.82 m s− 2) is gravitational acceleration; Sx (m m− 1) is rill slope steepness; it is equal to tan(θ),θ (°) is slope in degrees; and Sfx (m m− 1) is frictional slope given by (Chow, 1959): Sfx ¼
u2 1 C
ð3Þ
h
where C (m0.5 s− 1) is the Chezy coefficient given by (Chow, 1959): sffiffiffiffiffiffi 8g C¼ f
ð4Þ
f is the Darcy-Weisbach hydraulic frictional coefficient. Sediment continuity is given in the following equation: h
Ac Ac A Ac þ hu ¼ hDH þ ðSss −cσ Þ At Ax Ax Ax
ð5Þ
where c(x, t) (kg m− 3) is sediment concentration; DH (m2 s− 1) is the hydrodynamic diffusion coefficient for the sediments; and Sss(x, t) (kg m− 2 s− 1)) is the sediment source/sink term that, in detachmenttransport processes, is given by: cq Sss ¼ Dc 1− Tc
ð6Þ
where Dc (kg s− 1 m− 2) is potential detachment; Tc (kg m− 1 s− 1) is sediment transport capacity; and q (m2 s− 1) is unit flow discharge, while in depositional processes Sss is computed by: Sss ¼ −βðcq−Tc Þ
ð7Þ
−1
where β (1 m ) is an empirical depositional coefficient. 2.2. Model parameters 2.2.1. Sediment transport capacity Sediment transport capacity is defined as the maximum sediment load that can be carried along by a given water flow (Mancilla and Chen, 2005). Many attempts have been made to determine the sediment transport capacity based on data from various sources. Using laboratory generated experimental data for rill erosion, Nearing et al. (1997) related erosion to stream power, advanced by Bagnold (1966) as: ω ¼ ρgSx q
ð8Þ −3
−3
where ω (kg s ) is the stream power function; ρ (kg m ) is water density. Based on this concept, Nearing et al. (1997) suggested an equation to relate sediment transport capacity with stream power as: Be½aþblogðωÞ 1 þ e½aþblogðωÞ
2.1. Dynamic equations for hydrodynamic and sedimentation processes
log10 ðTc Þ ¼ A þ
Mass balance of water flow under rainfall and concentrated flow in the rill has been described in many studies (Milne-Thomson, 1960; Freeze, 1978; Ziekiewica and Taylor, 1991; Shames, 1992). Based on the research of Lei et al. (1998), the mass and momentum conservation/ balance equations are given by the following models:
where a = 0.845, b = 0.412, A = − 34.47, B = 38.61 were empirically derived coefficients. A laboratory flume experimental method for measuring transport capacity was advanced by Lei et al. (2001) where the transport capacity was linearly related to the slope and the flow rate as:
Ah AðuhÞ þ ¼σ At Ax
ð1Þ
Tc ¼ a þ bθ þ cQ
! AðuhÞ A u2 h g A h2 Sx þ þ ¼ gh pffiffiffiffiffiffiffiffiffiffiffiffiffiffi þ Sfx At 2 Ax Ax 1 þ S2x
ð2Þ
where x (m) is a Cartesian coordinate; t (s) is time; h(x, t) (m) is depth of water flow measured in the vertical direction; u(x, t) (m s− 1) is
ð9Þ
ð10Þ
where a = −0.3109, b = 0.01718, c = 0.1203 were regression coefficients and Q (L min− 1) was flow rate. 2.2.2. Rill erodibility The WEPP model provided a new erosion prediction tool based on water erosion physical processes. Rill erodibility is a key parameter for
116
L.J. Yan et al. / Geoderma 146 (2008) 114–120
rill erosion prediction in the WEPP model, as shown in the following equation: cq Dr ¼ Kr ðτ−τ c Þ 1− Tc −1
−2
Dc ¼ ð11Þ −2
where Dr (kg s m ) is the net detachment rate; τ and τc (N m ) are the shear stress of flowing water and the critical shear stress of soil, respectively; Kr (kg N− 1 s− 1) is the rill erodibility parameter; cq (kg m− 1 s− 1) is the sediment yield in the rill. The effects of soil erodibility (Kr) play a great role in soil erosion dynamic processes. The higher the Kr value, the more readily the soil particles are detached and, hence, more soil material is entrained in the water flowing from upstream rill sections. Consequently, there is a faster increase in the amount of sediment load carried in the flowing water and the transport capacity is approached more rapidly. When the sediment load tends to reach the transport capacity, the erosion intensity is reduced downstream. Eq. (11) shows that the net detachment rate has its maximum value, the potential detachment rate, when water is clear, i.e., when cq = 0, which is given by: Drmax ¼ Dc ¼ Kr ðτ−τc Þ:
ð12Þ
Values for Kr and τc were estimated by linear regression of a set of detachment rates and the corresponding hydraulic shears (Foster and Meyer, 1972; Nearing et al., 1989; Ghebreiyessus et al., 1994). A rational method for determining the soil erodibility and critical shear stress of rill erosion under concentrated flow was advanced by Lei et al. (2008). Using their method, soil erodibility of a loess soil was determined to be 0.3211 ± 0.0005 kg N− 1 s− 1, and the critical shear stresses τc were found to increase with slope gradient, i.e., τc was 3.191 N m− 2, 3.960 N m− 2, 4.130 N m− 2, 4.385 N m− 2, 4.565 N m− 2 for slopes of 5°, 10°, 15°, 20°, and 25° respectively. Gilley et al. (1993) reported a method for rill erodibility determination for WEPP. In their experiments to determine the parameter, simulated rainstorms were applied to field runoff plots, 0.46 by 9 m, to generate rills. Subsequently, concentrated flows of different flow rates were introduced into the rills to simulate rill erosion. The sediment yields thus obtained from the runoff plots were used to estimate detachment rates, for the Dc values in Eq. (12), by: Dc ¼
cq L
ð13Þ
where L (m) is the rill length, L = 9 m. The estimated detachment rates and corresponding hydraulic shear stresses under different hydraulic conditions were then used to determine rill erodibility and critical shear stress by regressing the data using Eq. (12). Zhang et al. (2001) experimentally confirmed the sedimentation process in a well-defined rill to be described in the following form: c ¼ A 1−e−Bx
ð14Þ
where A (kg m− 3) is a regression coefficient, representing the highest possible sediment concentration in the rill, under a given hydraulic regime; B is a decay coefficient, (1 m− 1). This B in Eq. (14) indicates how fast the sediment concentration reaches the value at the transport capacity, c = A. As indicated in Eq. (17), B is the pure soil parameter to determine the potential detachment rate, while Tc determines the hydraulic condition for potential detachment rate. Lei et al. (2001) also verified the following relationship: Aq ¼ Tc
Combining Eqs. (14) and (13) yields:
ð15Þ
where q (m2 s− 1) is unit width flow rate of the given flow; A is same as Eq. (14).
Tc 1−e−BL cq Aq 1−e−BL ¼ ¼ L L L
ð16Þ
Lei et al. (2002) also suggested a method to estimate the potential detachment rate: dA 1−e−Bx q ¼ ABq ¼ BTc dx xY0
Dc ¼ lim Dr ¼ lim xY0
ð17Þ
Eqs. (16) and (17) give two different methods by which to estimate rill erodibility. One, used by WEPP, is based upon Eq. (16), and is given by: Krw ¼
Tc 1−e−BL Dr : ¼ τ−τc ðτ−τc ÞL
ð18Þ
The other, suggested by Lei et al. (2002) and based upon Eq. (17), takes the following form: Krl ¼
Dr Tc B ¼ : τ−τc τ−τ c
ð19Þ
In order to verify the difference between the two erodibility parameters estimated by these two methods, Eq. (18) can be divided by Eq. (19) to give: T B
c Krl BL c ¼ T ðτ−τ −BL Þ ¼ c 1−e 1−e−BL Krw
ð20Þ
ðτ−τc ÞL
Krw ¼
Krl 1−e−BL BL
ð21Þ
where Krl is rill erodibility suggested by Lei and Krw is rill erodibility used by WEPP. This functional relationship (Eq. (21)) suggests that the estimate for soil erodibility used by WEPP depends upon the value of B, which is a complex coefficient affected by factors such as water flow, rill slope, and soil characteristics. Different B values lead to different soil erodibilities and in turn affect the simulation of the erosion process. In order to quantitatively understand the influences of B and Krw values on rill erosion processes, different Krw values and those computed from various B values for a 25° rill slope are listed in Table 1 to be used in this study. Experimental data, from a study by Zhang (2002) on a silty-clay loess soil, were used to evaluate the impacts of model parameters on rill erosion processes. The experiments were conducted in a flume, 8 m long and 1 m wide, which was sub-divided into strips, 8 m by 0.1 m, to imitate rills. The experiments involved five slopes (5°, 10°, 15°, 20°, and 25°), 9 slope lengths (0.5, 1, 2, 3, 4, 5, 6, 7, 8 m), and three flow rates (8 and 12 L min− 1 for 5° slopes; 2, 4, and 8 L min− 1 for steeper slopes). Numerical simulations were made for the same conditions as those applied in the experiments for the model verification. Simulated sediment concentrations along the rills were also compared with those from the experiments. Table 1 Krl and Krw under different slope gradients; different Krw values computed from different B values at 25° Slope
Krl (kg N− 1 s− 1)
Krw (kg N− 1 s− 1)
5° (B = 0.3) 10° (B = 0.2) 15° (B = 0.3) 20° (B = 0.4) 25° (B = 0.3)
0.3209 0.3211 0.3214 0.3214 0.3216
0.1109 0.1489 0.1110 0.0868 0.1111 (B = 0.3) 0.0707 (B = 0.5) 0.0446 (B = 0.8) 0.0357 (B = 1.0)
L.J. Yan et al. / Geoderma 146 (2008) 114–120
Fig. 1. a–m) Comparisons of simulated rill erosion processes with experimental data. 117
118
L.J. Yan et al. / Geoderma 146 (2008) 114–120
2.3. Boundary and initial conditions With reference to the experimental setup and the simulation of the hydraulic and erosion processes in a well-defined rill, the boundary conditions are listed as follows: 8 Vx ¼ V0 > > > 0 > > > > > > Hx0 ¼ H0 > > > > >
ð22Þ
Fig. 2. Comparison of sediment concentration for simulations and experiments.
where γ is specific gravity of water, 9800 (N m− 3); S is hydraulic slope, sin(θ); V0 (m s− 1), and H0 (m), the water velocity and depth at the inlet respectively, are interrelated and determined by the flow in-flow rate, q (m3 s− 1), as: q ¼ V0 H0 w
ð23Þ
where w (m) is rill width. The initial conditions for the governing equations are given as: 8 < Vxi ðx; t ¼ 0Þ ¼ Vi ð24Þ H ðx; t ¼ 0Þ ¼ Hi : i ci ðx; t ¼ 0Þ ¼ 0 where Vi (m s− 1) and Hi (m), the initial water velocity and depth values along the rill, can be assigned arbitrarily with any values, except negative and zero Hi values. 3. Results and discussion Using the numerical formulations discussed above, a program was written in Visual C++ to simulate the rill erosion processes using the Finite Element method. Various model parameter values, for the transport capacity, rill erodibility and critical shear stress were used in the simulation with the aim of comparing the impacts of these different parameters, as estimated by different methods, on rill erosion processes. 3.1. Comparison of simulated rill erosion processes with experimental data The simulated sediment concentrations along rills, under different slope and flow rate conditions with parameter sets including transport capacity, rill erodibility and critical shear stress estimated by Lei et al. (2001, 2008), were compared with the experimental data collected under the same conditions, as shown by the Rational Tc and Kr line in Fig. 1a–m. The sediment concentration hygrographs, for both the simulated and the experimental data, were in close agreement with the relationship between sediment concentration and rill length described by Lei et al. (1998) and Huang et al. (1996). Fig. 1a–m for all the combinations of slope gradients and flow rates indicate very good agreement between the simulated results and those from the experiments. Sediment concentration increases with rill length. In the upper part of the rill, i.e., when the down-slope distance is short (Fig. 1a to m), the greatest increases in sediment concentration were observed indicating higher detachment rates. The rate of increase in sediment concentration becomes lower with increasing rill length, and sediment concentration tends to approach a stable maximum value that is limited by the transport capacity of the water flow. Using all the data sets consisting of 117 data points from the 13 experimental runs, a comparison was made against simulated data, as
shown in Fig. 2. Linear regression indicated that the simulated sediment concentrations were closely correlated with the measured data, with a slope coefficient of 0.989 and a coefficient of determination of 0.90. Therefore, there exists a very good one to one relationship between the two sets of data. These results imply that the simulation model, including the mathematical formulations and the numerical algorithms of the Finite Element method, as well as the computational procedure codes, were correct. Moreover, the model parameters set, including the transport capacity, rill erodibility, and the critical shear stress estimated by Lei et al. (2001, 2008), must likewise be correct in order to produce such good predictions of sediment yields not only at the rill outlets but also at the various locations along the rills. 3.2. Impacts of transport capacity on rill erosion processes Simulations were also made, using the transport capacity given in Eq. (9) reported by Lei et al., (1998) and the rill erodibility and critical shear stress estimated by Lei et al. (2008), to check the impacts of this parameter on the rill sedimentation processes. The results are indicated by the Tc — Limited curves in Fig. 1a–m. Comparisons were made between the experimentally measured data and the simulated data sets. The comparisons show that generally the simulated sediment concentrations in the water flow in the upper parts of the rill (b1 m down-slope distance) were in agreement with the experimental data. The sediment concentration increases rapidly along this section of the rill. As sediment concentrations approach the limiting transport capacity of the water flow, the simulated sedimentation loads gradually and then rapidly deviate from the experimental data to considerably lower amounts. At the end of the rills, the simulated sediment yields were dramatically reduced as compared with the measured data. These simulated results indicated how the transport capacity parameter limited the sedimentation processes behavior. This might, to some extent, have explained the reason why calibration of the transport capacity was needed before a reasonably good prediction of sediment concentrations at rill outlets was obtained by the rill erosion model of Lei et al. (1998). These results indicate that the transport capacity parameter is of great importance in rill erosion process modeling. 3.3. Impacts of rill erodibility on rill erosion processes To check the impact of the erodibility parameter on rill sedimentation processes, simulations were also made with rill erodibility values estimated by the method used in the WEPP model, as given by Eq. (21), with a B value of 0.2, 0.3, and 0.4 for various slopes. The simulated results are depicted by the Kr — Influence curves in Fig. 1a– m. The results show that the simulated sediment concentration curves consistently underestimated rill erosion considerably for all rill sections. Therefore, rill erodibility estimated by the WEPP method could be too low to make a good prediction of rill sedimentation processes.
L.J. Yan et al. / Geoderma 146 (2008) 114–120
Fig. 3. Simulated results between different Krw values computed from different B values for a 25o slope.
Nevertheless, Fig. 1a–m also showed that, generally, the sediment concentration curves tended to increase gradually down the length of the rills, approaching the level of the potential sediment concentrations as determined by the transport capacity. Potentially, as rill length further increases, the sediment concentration may reach the levels defined by the transport capacity, provided the rills are long enough. These results imply that, although the rill erodibility value might be too low to estimate the sedimentation processes along a short rill, the sediment yields at rill outlets could be reasonably close to the actual values provided the rill is long enough and the transport capacity is correctly defined. The overall error of the sediment concentrations at the rill outlets between the simulated results and the experimental data for the rills of 8 m long was 30%. These simulation results indicate that, although the lowered values of estimated rill erodibility could give reasonably good predictions of sediment yields at the outlets of a long rill, if the transport capacity is correctly defined, a good estimation of sediment yields at rill outlets alone cannot guarantee that there was good simulation of erosion dynamics along the rills, unless comparisons are made along the length of the rill, as carried out in the previous section. Thus, simulated results estimated by WEPP might need reappraisal with sedimentation data collected along the rills instead of using the sediment yields at outlets alone. A comparison was also made between the experimentally measured rill sedimentation processes, at a slope gradient of 25° and a flow rate of 4 L min− 1, and the simulated results predicted when using rill erodibilities estimated by Eq. (21) for different B values (0.3, 0.5, 0.8 and 1), as used by the WEPP model, as well as the simulated result for rill erodibility estimated by Lei et al. Eq. (19). The transport capacity was given by Eq. (10). The results are shown in Fig. 3 and indicate that, with an increase in the B value, rill erodibility as determined by Eq. (21) decreases, and that the estimated sedimentation loads were considerably reduced. The sediment yield at the outlet of an 8 m long rill is reduced by 70% when B = 1 compared with the results produced using a rill erodibility value suggested by Lei et al. (2008). That is, the sediment yield at the rill outlet predicted using the rill erodibility value of Lei et al. (2008) was more than twice that calculated with the erodibility estimated by the WEPP model method. 3.4. Combined impacts of transport capacity and rill erodibility on erosion processes Further comparisons were made using various transport capacity and rill erodibility values to see the combined impacts of these two important parameters on rill erosion dynamics. The results are illustrated in Fig. 1a–m. The results indicate that the simulated sedi-
119
mentation processes are in good agreement with the experimental data when used with the correctly correctly-estimated model parameters of Lei et al. (2008). When the simulation is limited by transport capacity alone (i.e., with the rill erodibility correctly estimated), the sedimentation processes in the upper sections of the rill relate well to the experimental data. However, the simulation derived curves gradually diverge from the experimental data as down-slope distance increases. The predicted sediment loads finally reach the limit defined by the imposed transport capacity value. This would inevitably cause low estimations of the sediment yields at the rill ends as well as for the sediment detachment rates along the rills. In the cases where the rill erodibility value influenced the model, the simulated sedimentation processes differed from the experimental curves for all sections of the rills. However, the sediment load generally continued to increase with down-slope distance towards a level that corresponded to the transport capacity, which potentially could be attained if the rills were long enough. Thus, the estimated sediment yields at the rill outlets could be very close to the actual value for long rills. However, the estimated sediment loads from specific rill locations would differ from those values derived used a correctlyestimated rill erodibility value. 4. Conclusions Mathematic models were formulated based on hydrodynamics and physical sedimentation concepts, outlined by Nearing et al. (1997) and advanced by Lei et al. (1998b). Numerical simulations were made using the Finite Element method and a Visual C++ program to simulate rill erosion dynamics under different rill slope and hydraulic conditions. Comparisons made between the experimental data and the simulated results showed that the models, numerical algorithms, and the program all work well. Rill sedimentation processes were accurately simulated by the model (as given by Eq. (19)) and rational Kr and Tc. Sediment concentration data from the simulations were closely correlated with the experimental data and these very good predictions were made without calibration of the model parameters, which indicated that the model parameters estimated were correct and emphasized the importance of model parameter estimation. Reducing the transport capacity values produced transport capacity limited sedimentation processes and lowered sediment yield estimations below those measured. Simulation using lower rill erodibility values could make reasonably good predictions of sediment yields at the outlets of long rills given good estimated transport capacity values. However, the sediment sources might not be the same as in real situations. Checking the simulation results of both sediment yields at the rill outlets and the sedimentation processes along the rills is important. This study not only supplies a tool for further development of dynamic soil erosion prediction models but also verifies the importance of correct estimation of model parameters. Acknowledgments This paper is based on the work supported by the Key Laboratory of Soil and Water Conservation & Desertification Combat (Beijing Forest University), Ministry of Education (2005k-01-01), by the Natural Science Foundation of China (40635027), Program for Changjiang Scholars and Innovative Research Team in University, as well as supported by the CAS/SAFEA International Partnership Program for Creative Research Teams, 2007-380. References Bagnold, R.A., 1966. An approach to the sediment transport problem from general physics. US Geological Survey Paper 422-∣, Washington.
120
L.J. Yan et al. / Geoderma 146 (2008) 114–120
Chow, V.T., 1959. Open-channel Hydraulics. McGraw-Hill, New York, N Y. Flanagan, D.C., Ascough, J.C., Nicks, A.D., Nearing, M.A., Laflen, J.M., 1995. Overview of the WEPP erosion prediction model, technical documentation. National Soil Erosion Research Laboratory, West Lafayette.USDA-Water Erosion Prediction Project. Foster, G.R., Meyer, L.D., 1972. Transport of soil particles by shallow flow. Transactions of the ASAE 15, 99–102. Freeze, R.A., 1978. Mathematical models of hillslope hydrology. In: Kirkby, M.J. (Ed.), Hillslope Hydrology. Wiley Interscience, New York, pp. 177–225. Ghebreiyessus, Y.T., Gantzer, C.J., Alberts, E.E., Lentz, R.W., 1994. Soil erosion by concentrated flow: shear stress and bulk density. Transactions of the ASAE 37, 1791–1797. Gilley, J.E., Elliot, W.J., Laflen, J.M., Simanton, J.R., 1993. Critical shear stress and critical flow rates for initiation of rilling. Journal of Hydrology 142, 251–271. Hairsine, P.B., Rose, C.W., 1992. Modeling water erosion due to overland flow using physical principles 2: rill flow. Water Resources Research 28, 245–250. Huang, C.H., Laflen, J.M., Bradford, J.M., 1996. Evaluation of the detachment-transport coupling concept in the WEPP rill erosion equation. Soil Sci. Soc. Am. J. 60, 734–739. Lei, T.W., Nearing, M.A., Haghighi, K., Bralts, V.F., 1998. Rill erosion and morphological evolution: a simulation model. Water Resources Research 34, 3157–3168. Lei, T.W., Zhang, Q.W., Yan, L.J., Zhao, J., Pan, Y.H., 2008. A rational method for estimating erodibility and critical shear stress of an eroding rill. Geoderma 244 (3–4), 628–633. Lei, T.W., Zhang, Q.W., Zhao, J., Nearing, M.A., 2002. Analytic method for determination of detachment rate of concentrated flow in erosion rills. Acta Pedologica Sinica 39, 788–793 (in Chinese with English abstract). Lei, T.W., Zhang, Q.W., Zhao, J., Tang, Z.J., 2001. Laboratory study on sediment transport capacity in the dynamic process of rill erosion. Transactions of the ASAE 44, 1537–1542. Mancilla, G.A., Chen, S.L., 2005. A field study on sediment transport capacity on rills. ASAE Annual International Meeting Sponsored by ASAE. Tampa, Florida.Paper No. 052051. ASAE, St Joseph, Mich. Milne-Thomson, L.M., 1960. Theoretical Hydrodynamics. Macmillan, New York. Morgan, R.P.C., 1994. The European soil erosion model: an update on its structure and research base. In: Rickson, R.J. (Ed.), Conserving Soil Resources: European
Perspectives. Selected Papers from the First International Congress of the European Society for Soil Conservation. Wallingford. UK, pp. 286–299. Morgan, R.P.C., Quinton, J.N., Rickson, R.J., 1992. EUROSEM: Documentation Manual. Silsoe College, Silsoe, UK. Nearing, M.A., Foster, G.R., Lane, L.J., Finkner, S.C., 1989. A process based soil erosion model for USDA. Water Erosion Prediction Project Technology. Transaction of the ASAE 32, 1587–1593. Nearing, M.A., Norton, L.D., Bulgakov, D.A., Larionov, G.A., West, L.T., Dontsova, K.M., 1997. Hydraulics and erosion in eroding rills. Water Resources Research. Res 33, 865–876. Shames, I.H., 1992. Mechanics of Fluids, 3rd Edn. McGraw-Hill Book Company, New York. Sharda, V.N., Nearing, M.A., 1999. Finite element modeling of erosion on agricultural lands. In: Stott, D.E., Mohtar, R.H., Steinhardt, G.C. (Eds.), Sustaining the Global Farm. Selected Papers from the 10th International Soil Conservation Organization Meeting. Purdue University and the USDA-ARS National Soil Erosion Research Laboratory, West Lafayette, IN, pp. 1008–1017. Sharda, V.N., Singh, S.R., Sastry, G., Dhruvanarayana, V.V., 1994. A finite element model for simulating runoff and soil erosion from mechanically treated agricultural lands 2. Field validation and applications. Water Resources Research 30, 2299–2310. Zhang, Q.W., Lei, T.W., Pan, Y.H., Gao, P.L., 2002. Rational computational method of soil erodibility and critical shear stress from experiment data. Journal of the Graduate School of the Chinese Academy of Sciences 21, 468–475 (in Chinese with English abstract). Zhang, Q.W., Lei, T.W., Zhao, J., 2001. Coupling relationship between detachment rate with sediment load and rill length. Journal of soil and water conservation 15, 92–95 (in Chinese with English abstract). Zheng, F.L., Tang, K.L., 1997. Rill erosion process of steep slope land of the Loess Plateau. International Journal of S. Ziekiewica, C., Taylor, R.L., 1991. 4th Edn. The Finite Element Method, vol. 2. McGrawHill Book Company, New York.