- Email: [email protected]
The description of soil erosion through a kinematic wave model
Journal of Hydrology, 145 (1993) 65-82
65
Elsevier Science Publishers B.V., A m s t e r d a m
[3]
The description of soil erosion through a kinematic wave model A. Laguna" and J.V. Girgldezb aDepartment of Applied Physics, University of C6rdoba, Apdo. 3048, 14080 C6rdoba, Spain bDepartment of Agronomy, University of C6rdoba, Apdo. 3048, 14080 C6rdoba, Spain (Received 21 May 1992; revision accepted 20 October 1992)
ABSTRACT Laguna, A. and Gir~ildez, J.V., 1993. The description of soil erosion through a kinematic wave model. J. Hydrol., 145: 65-82. Kinematic wave models offer a powerful tool for interpreting erosion experiments and for forecasting soil and water losses under different management systems. One important problem is the identification of the model parameters. To estimate the influence of error in the estimation of parameters a sensitivity analysis was performed. A previous survey of the soil erosion literature provided data to establish the expected ranges for each parameter. Of the three main parameters, the interrill and rill erosion coefficients, and the rill water depth coefficient, the last is the most important since sediment concentration is most sensitive to variations in this parameter, although the range of variation in the model is affected by the reference values of the other parameters. A simulated rainfall experiment in rectangular plots of 5 × 15 m 2 on a 20% slope provided data which were fitted to the model. The agreement between the model and the data was reasonably good. Sediment yield was related linearly to total runoff volume in individual simulation events.
INTRODUCTION
Soil erosion poses a serious threat to agriculture and the environment. Erosion removes the most fertile fraction of the soil, thus reducing crop productivity. Sediments plug canals and culverts, fill up reservoirs, and carry away chemicals adsorbed in soil particles, thereby polluting the environment. In Mediterranean climates, water erosion is especially severe, since after a warm dry summer there is only sparse vegetation cover to protect the soil surface from early autumn rains. Agricultural practices exacerbate the situation by enhancing the erosive effect of water. Farmers frequently plough the soil to control weeds. The topsoil becomes a loose layer on top of a layer C o r r e s p o n d e n c e to: A. Laguna, D e p a r t m e n t of Applied Physics, University o f C 6 r d o b a , Apdo. 3048, 14080 C 6 r d o b a , Spain.
0022-1694/93/$06.00
© 1993 - - Elsevier Science Publishers B.V. All rights reserved
66
A. L A G U N A A N D ,I.V. G I R A L D E Z
NOTATION a
b B C
C
El & g h i ia
k K
Ke Kr L m n
q, qe qp r
S
s, SL !
tc le /r 0 lsc
T
L Tr Vc
v~ x X,
exponent exponent interrill erosion coefficient [M L -m-2 T ~- D] sediment concentration [M L -3] cover factor from the USLE interrill erosion [M L 2 T ~] rill erosion[M L 2 T 1] acceleration of gravity [L T 2] flow depth [L] infiltration rate [L T ~] average water abstraction [L T -t] rill erosion coefficient [M L i b T i] soil erodibility factor from the USLE coefficient [T 2] coefficient [L ~ T] slope length [L] exponent Manning's coefficient water discharge rate per unit width [L 2 T -~] equilibrium runoff flow rate per unit width [L 2 T ~] rainfall intensity [L T i] total sediment yield [M] sensitivity of c to P~ parameter plane slope time [T] water concentration time [T] runoff equilibrium time [T] runoff initiation time [T] sediment concentration time [T] rainfall duration [T] transport capacity [M L- t T - t ] runoff duration [T] settling velocity [L T 1] total runoff volume [L 3] distance along the flow direction [L] threshold of slope erosion [L]
Greek letters
depth-discharge coefficient [L 2 a T q ] rill erosion coefficient [L i] O. Oo initial and final soil moisture [M M ~] coefficient [M L -2 T -~] water density [M L -3] P critical shear stress of the flow [M L ~ T 2] •c shear stress of the flow [M L 1 T-2] f~ stream power [M T -3] threshold stream power [M T 3] f~o Ot
7
K I N E M A T I C WAVE M O D E L O F SOIL EROSION
67
of hard plough sole penetrable only with great difficulty by water and plant roots. When water filters into the soil, buoyancy forces raise the topsoil above the plough sole and carry it downslope. The resulting gullies are thenceforth U-shaped. Nevertheless, agriculture is not the only participant responsible for soil erosion. Forestry practices require the establishment of roads and fire breaks along mountain slopes, practices which require clearing the vegetation. Civil works such as the opening of new roads or simply urbanizing vegetated areas remove plant protection from the soils, thus accelerating natural erosion. In spite of soil erosion being a menace for agriculture and the environment, the number of studies in many countries is surprisingly low. Blaikie (1985, section 2.2) remarked on this deficiency, especially in less-developed countries, which are severely affected by this problem. There is general concern with respect to the effects of erosion, or to its related problem of land degradation and desert encroachment or even desertification, but in most of the cases the real extent of the problem is not known precisely. The use of an erosion model might be very helpful for predicting the consequences of soil management practices, and land-use changes. The complexity of the phenomena makes the formulation of a simple model difficult, even when a single factor like rain splash of water runoff is the only agent of erosion. Many researchers prefer to use empirical relations between observed soil loss rates or sediment gains with selected soil or environmental parameters, via a statistical equation. Although this approach is convenient for a rough estimate, its use is not recommended in conditions differing from those considered originally. An appropriate way of studying soil erosion processes is through the formulation of the fundamental transport equations of water and sediment (e.g. Bennett, 1974, or Rose, 1985). Once posed, these equations can be solved either by numerical or analytical methods. The kinematic wave simplification for surface water flow on a plane yields an analytical solution for the cases of constant rainfall and infiltration rates, as Singh and Prasad (1982), among others, have demonstrated. Other authors, such as Govindaraju and Kavvas (1991), have worked on an approximate solution of a diffusion wave model previously derived by Govindaraju et al. (1990). This is a more complete hydrological model since it allows for a varying time infiltration rate in discrete time steps, and assumes a constant sediment concentration for the whole erosion processs. In spite of being very useful from the hydrological point of view, the assumption of constancy of sediment concentration is a shortcoming of this model. The kinematic wave model may be applied to the interpretation of erosion experiments after their parameters have been identified. The purpose of this report is to explore the fit of the kinematic wave
68
A LAGUNA AND J.V GIRALDEZ
model to experimental results under controlled conditions and to analyse its sensitivity to the most relevant parameters. T H E K I N E M A T I C SOIL E R O S I O N M O D E L
The model based on the kinematic wave approach to the Saint-Venant equations consists of three expressions. Expression (1) is the mass conservation expression for water 0h 0q a--t + 0x
r -- i
(1)
where h is the water depth [L]; q is the flow rate per unit width [L 2 T - ~]; r and i are the rainfall and infiltration rates respectively [L T - ~]; t refers to time [T], and x to space [L]. The relationship between depth and flow rate constitutes the m o m e n t u m conservation equation (expression 2): q =
~" h ~
(2)
where ~ is a coefficient related to Manning's n or Chezy's C, [L 2-" T ~], and a is an exponent. The mass conservation equation for sediment furnishes the last expression
(3): a(hc)
c~----~+
a(qc)
8x
-
El + ER
(3a)
where c is the sediment concentration in water flow [M L-3]. El and ER respectively represent the sediment production rate due to interrill and rill erosion processes [M L -2 T-~]. There are several proposals for the El and E R terms. Singh and Prasad (1982) selected a potential expression o f El with respect to the rainfall rate: E,
=
(3b)
B . rm
where B is a coefficient and m an exponent. Rill erosion depends on the water depth in the rill, and on the flow rate which will induce deposition if it is not too high: ER
---
?(k" h b -
c " q)
(3c)
where 7 is another parameter and b an exponent which according to Shirley and Lane's (1978) findings is close to a. Therefore it would be assumed that b=a.
Other approaches to the expression of sediment generation by the interaction of erosion-deposition processes have been proposed by Rose et al. (1983)
69
KINEMATIC WAVE MODEL OF SOIL EROSION
and Croley and Foster (1984). Croley and Foster (1984) proposed a similar equation to (3c) which was later modified by Nearing et al. (1989) in the U S D A - W E P P framework. Their expression is ER =
3" ( 1
(3d)
qTc c )
where Tc is the transport capacity of rill flow and ~ is a coefficient whose value is either Kr" (rr - r~) when the term in parentheses in (3d) is positive, or (v c • T~)/q when it is negative. The settling velocity is Vc; rf represents the shear stress of the flow, and its critical value is r~. Transport capacity is related to the flow's shear stress, -of, through a potential equation. Alternatively Rose et al. (1983) separated the distinct contributions of the rainfall detachment rate, essentially (3b); they use the sediment entrainment rate, and sediment deposition rate, Vc • c. Bagnold stream power is defined as the shear stress by the mean flow velocity, 92 = p • g . SL " qL " x / L (e.g. Bagnold, 1977). In that expression p is water density, qL is the flow rate at the end of the plane, and SL its slope. This equation was adopted as the parameter controlling the net bed load transport rate, over a certain threshold value 920, defined as the stream power for the distance at which there was no sediment entrainment, x , . The final expression for ER is ER
=
p " g " SL " K~ "
+
qL
1 -- - -
X
+
(hc)
(3e)
The inclusion of the partial derivative of the product of water depth by sediment concentration with respect to time allows the elimination of such a term in eqn. (3a). This converts a partial differential equation into an ordinary differential equation of simpler solution. Time-dependent parameters are evaluated in steps. Different size sediments are introduced once their mass fraction is known. For the purposes of the present work the expression (3c) is preferred since its integration is straightforward. The initial conditions and boundary conditions are as follows, where T is rainfall duration, and L is the length of the slope: =
O,
0 <<. x <<. L
(4a)
h(0, t) =
0,
t >~ 0
(4b)
h(x,O)
The rainfall is assumed to be homogeneously distributed over the whole surface, a plane, with constant rate: r(x,t)
=
r > O,
r(x, t)
=
O,
0 < t < T t >1 T
(5a) (5b)
70
A. LAGUNA AND J.V. GIR,~LDEZ
Water infiltrates the soil at a constant rate, whenever the water depth over the soil is positive: i(x, t) =
i > O,
i(x, t) =
O,
(6a)
h(x, t) > 0 h(x, t) =
0
(6b)
Note that r must be greater than i. Solutions o f the equations
The solution of eqns. (1)-(3) subjected to conditions (4)-(6) may be found through the m e t h o d of characteristics. The solutions are found for the pair of values (x, t), whose paths are the characteristics. The times at which the characteristics, starting at the origin of coordinates in the t - x plane, reach the downstream boundary are k n o w n as the concentration times, to, and tsc, for water and sediments respectively. Singh and Prasad (1982) developed solutions for the different domains b o u n d e d by the limiting characteristics of water and sediment, the recession curve and the lines t = T, t = 0, x = 0, and x = L. There are a few misprints; in their eqn. (60), the characteristic for the sediments in domain D3b, is corrected to t -
T
=
r(a -
]) + i
"{[(r - i)xx] '/" - (r'x'o -
i ' x ) '/a}
(7)
and the value of the auxiliary function #(xT, T) used in their eqn. (64) for the sediment concentration in domain D3a should be /~(xx, T)
=
7~(r-
i) Ca 1)[ ~ ( r -
i)¢' ~)la(xw/~)'l~-- T ] a -- 7
(8)
M o d e l parameters
The solution of eqns. (1)-(3) will characterize completely the soil erosion processes for any point in the plane at all times. The main problem in applying the model is the identification of the parameters, since they are not immediately measured at the experimental site. In this section the pertinent literature will be reviewed, looking for the reported range of variation as a preliminary step in the sensitivity analysis of the model. There are six parameters in all: the ~, B, 7 and k coefficients, and the a and m exponents. Parameters o f the dynamic equation Parameters ~ and a determine surface water movement. Under a laminar
K I N E M A T I C WAVE M O D E L O F SOIL EROSION
71
TABLE 1 Selected values of the coefficient ~ for some slopes and surfaces Surface
Manning's n
Bare clay-loam soil (eroded)
0.012-0.033
Sparse vegetation cover
0.053-0.130
Short grass prairie
0.100-0.200
Bluegrass Sod
0.170-0.480
Slope
(%)
(m~/:s ' )
l0 20 10 20 10 20 10 20
26.35-9.58 37.27-13.55 5.97-2.43 8.44-3.44 3.16-1.58 4.47-2.24 1.86-0.66 2.63-0.93
is computed from eqn. (9), with Manning's n values from Woolhiser (1975).
flow regime, a assumes a value of 3 in the Darcy-Weisbach formula, whereas under a turbulent regime, values of 5/3 in Manning's equation, and 1.5 in Ch6zy's equation have been proposed. However, Croley and Foster (1984) estimated values of a close to 0.25 and argued the convenience of the parameter being smaller than unity. These values were computed from the estimates of the time of concentration and the average velocity of water flow in the rills of experimental erosion plots. Croley and Foster justified their estimates on the basis of the irregularity of canal geometry. Nevertheless, values for a that are smaller than unity imply multiple values of water depth at a given time in the recession stage. Huggins and Burney (1982) supported Manning's exponent. Singh and Regl (1983) proposed 1.5 on the grounds of its commonly accepted usage. In this study, a 1.5 value was used for a. Coefficient ~ has a wider range than a. This coefficient depends basically on the slope and on surface roughness. Manning's equation is
(&),/2 -
n
(9)
where SL is the slope and n is Manning's coefficient. Some possible values for ~, for different slopes and surfaces are gathered in Table 1. They have been computed from eqn. (9) with Manning's n values from Woolhiser (1975). Reasonable values for ~ are between 0.5 and 40 m ~/2s -1 . Figure 1 shows the concentration of sediments in the model for several ~ values within the given interval, keeping the other parameters constant. The selected values for rainfall and infiltration rate were somewhat high, since one of the objectives of the experiment was the production of high erosion events. The main effect of the variation of ~ is in time of concentration. The higher ~ is (this implies
A. L A G U N A AND J.V. G I R A L D E Z
72 50
75 B = 5.05e05 k9 s-'m 7 =0.1 m -~
•
= 10 m i ~ s-' 7 = 0.1 m - ' k = 50 kg m -=,a S-'
\
k = 50 kg m-~I~ s - '
'~kt)
B (kg s
5O
E ~25
E
m-4 )
7.525e5 5.05e5
v
I .Oe4
u 25
40
.
.
.
3~0
.
.
Time
50j |
a =
.
.
6~0
.
,
,
,
(s)
,
310~
50 |
10 m ~r~ s -I
kB = 5 . 0 5 e 0 5 kg s - ' en-° 50 kg m-=/z s - '
(m-') .~.
0.055 0.1
~25
_~ 25
0.3 0.5
v
t
'
,
,
6O0
cr = 10 rn '/= s -~ B = 5 . 0 5 e 0 5 kg s - ' rn ~ 7 = 0.1 m - '
k (kg m "fC*s-')
{il,o2~o 75 50
(5) 0
~
(s)
Time
I[4)r ~ '
'
'
'
0
I 300
'
'
r
0
, 600
,
,
0
Time (s)
Fig.
I. S e d i m e n t
r =
1.376e-5ms -t , L =
concentration
,
t
I 300
Time
for s e v e r a l c(, B, 7 a n d k v a l u e s (n =
'
f
, ,"'T'--~ ' 600
(s)
1.5, m = 2, i = 9 . 6 5 e - 6 m s
~,
1 0 0 m , T = 5 0 0 s ) . T i m e f r o m the start o f rainfall a n d runoff.
steep slope and smooth surface), the shorter is the time of concentration. For the lowest c( value the concentration time is larger than rainfall duration; in that case the steady state is not reached. Sediment concentrations in the steady state are surprisingly high for small values. Since the peak flow rate is independent of this parameter, the sediment load [M T - 1 L- 1] will also be high in the case of small ~. One can expect that for small ~ values, surface flow velocity is low and so, consequently, is the erosive capacity of the flow. A possible explanation for this paradox is that the k coefficient in the transport capacity of the flow (eqn. (3c)) must be linked to ~. Therefore if ~ changes, k should change as well, to a higher rate. After a certain threshold value, there are no great variations among the curves. These differences are greater the smaller 0( is. If the value of ~ is very small the curves increase with time.
Parameters of the interrill erosion equation Parameters B and m influence interrill or rain splash erosion. Singh and Regl (1983) and Blau et al. (1988) proposed a value o f m = 1. Croley (1982)
KINEMATIC
WAVE MODEL
73
OF SOIL EROSION
referred to Foster's (1982) preference for m = 2 based on a more realistic interpretation of splash experiments. The latter author suggested the experimental expression E l ( k g m - 2 h i) =
O.O138.Kl(kghN-im-2).C.r2(mmh
i)
(10)
where Kj is the soil erodibility factor for raindrop impact. In the absence of K, values, the universal soil loss equation K factor may be adopted, even though this factor assumes both interrill and intrarill erosion processes. The C factor is the cover factor from the universal soil loss equation. The parameter B may be evaluated experimentally or estimated through expressions similar to eqn. (10), for further refinement after calibration experiments. Croley and Foster (1984) obtained values for B between 1 x 106 and 1.5 x 106kg s m -2 after examining the results of an experimental plot under simulated rainfall. From the results of the model, for the time zero c =
B.r m : r
--
(11)
l
This equation implies that the initial concentration of sediments is proportional to B. After eqn. (10) a reasonable range for the values of B is (1 x 104 to 1 x 106kgs Im-4), which corresponds to the pair of values for K a n d C of (0.02, 0.01) and (0.04, 0.6) respectively.
Parameters of the intrarill erosion equation Parameters Y and k are associated with the intrarill erosion process. According to eqn. (3c) the sediment yield is proportional to the difference between transport capacity and sediment load. Parameter y [L ~] assumes values in the range 3-33 m - i in the work of Foster and Huggins, and David as quoted by Foster (1982). For y values higher than unity most of the resulting sediment concentrations are negligible. Singh and Regl (1983) suggested a reasonable value of 0.03 m -I . Parameter k, [M L -(a+l) T-l], is related to rill erodibility. Singh and Regl (1983) proposed a value of 1.5 kg s ' m -25. I f k is smaller than unity there are no differences in the results for different values of the parameter. In the rest of the work, values of a = 1.5 and m = 2 will be assumed. As was stated previously in eqn. (11) and as Fig. 1 shows, B is proportional to the concentration of sediments at time zero. In the steady state, holding the rest of the parameters constant, both sediment concentration and the value of parameter B follow a straight line with positive slope. The same is true of parameter k. Nevertheless y is inversely related to sediment concentration in the steady state. If this parameter assumes small values, both net soil loss and
74
A. LAGUNA AND J.V. GIRALDEZ
sediment deposition are relatively small. In such a case, changes from the initial conditions take place slowly (Blau et al., 1988). Sediment concentration increases monotonically with time from initial conditions until the steady state if B'~
k
r~ - -
<
r-i
1
(12)
since the whole phenomenon is controlled by rill erosion, decreasing in the negative case. ANALYSIS OF THE MODEL SENSITIVITY
The determination of the model parameters is difficult because most of them are not directly measurable. Therefore these parameters need to be estimated through indirect measurements in erosion experiments. In these cases a sensitivity analysis may ascertain the consequences to the final result of a small change in the value of the parameter. Since the final result in the present work is the sediment concentration c = f(PI, P2,...,Pn)
(13)
where PI, P2,. • • P, are the parameters and input variables for the model. The sensitivity of c to any one of the P~ is (Beven, 1979) ~C
Pi
. . c Si .- . 6p
(14)
or the equivalent S~ =
Ac _Pi AP" c
(15)
representing the fractional change in P transmitted to c, which affords a comparison of sensitivity with parameters whose range is of different order of magnitude. One shortcoming of this method is the assumption of parameter independence. Among the five model parameters, three were selected for the sensitivity analysis since it was assumed a fixed value of 1.5 for the exponent and the coefficient ~ is easily obtained from eqn. (9). Sensitivity to one parameter may change as the reference values of the other parameters change. Nearing et al. (1990) computed the sensitivity of the model to the parameters by applying the extreme values of the respective ranges. Therefore the sensitivity of the sediment concentration with respect to parameters B, 7 and k was computed for the maximum, minimum, and intermediate values found in the literature and indicated above as a function of time, as shown in Fig. 2. The
KINEMATIC WAVE MODEL OF SOIL EROSION
75
1.5 7=0.1 k = 50
rn-, kg rn-vz
s-' • B
~.0
• ~
03
•
®
=
le06
k~q
a B
5.05e0~
0 (3 =
le04
s -' m-*
o
>,
O
o
>
o
e
®
o
® o
05
o
m o
wl
a
0 0
O0
:
i
o
e
e
e
®
e
e
Q
o
o
o
~
o
o
o
® o
*°°ooo,ooooo,~ i
i
I
i
]
]
T
600
300 Time
(s)
].5-
15 B k
= =
505e05 kg 50 kg m-$/,
s-' s-'
B=§OSeO$kgs-,m 7 = 0 1 m - *
rn-,
-10
10
z>
• 7
=
0.5
a 7
=
0 "y
:
01 001
r n -~
•
k
o ¢
k= k =
=
100
kg
m-s/~
~@ • ® ® o o ® m e m ® o
8o00°°°**0oo o
e
o e
00
e
®
o
o
o
o
o
o
o
e
0
0.0
i
o
,
~
,
[
300 Time
(s)
i
0
o o o o o o o
@ o a
O.S u3
®
@
a
o
o ~^00
°
i
o o
°
o e e e e e o ® e m
i
o
o °
a
s-~
50 I ®o
-0.5
~
=l
0.0
'
'
l
0000000000 o° l
l
r
l
600 Time
l
l
l
I
300 (s)
Fig. 2. Sensitivity coefficients of sediment concentration to B, ;, a n d k p a r a m e t e r s = 10mW2s ~ , i = 9 . 6 5 e - 6 m s - ] ~ r = 1 . 3 7 6 e - 5 m s i L = 100m).
600
(n =
1.5, m = 2,
values of rainfall and infiltration rate adopted were, as in Fig. 1, obtained in the experiments described below in the erosion plots. In the curves shown in Fig. 2 the variation over time in the sensitivity coefficient to any parameter can be seen. In the rising stage, erosion processes are controlled by the detachment of soil particles by raindrop impact. At this stage there is a maximum value in sensitivity to the interrill erosion parameter B, although the sensitivity to the other parameters is low. However, in the recession stage, when soil loss is primarily due to shear by water flow, the sensitivity with respect to B is smaller than during other stages, and the sensitivity to k acquires its maximum value. The highest values for the sensitivity coefficient are found for the parameter values that imply greater sediment concentrations. The sensitivity coefficient for parameter 7 is negative since its value increases when the value of c diminishes. These results agree with the sensitivity analysis of Blau et al. (1988): the higher the value of parameter 7, the lower the sensitivity of the model. In the steady state where the greater fraction of soil loss occurs, the most important parameter in terms
76
A. LAGUNA AND J.V. GIR,~LDEZ
of model sensitivity appears to be k, except for the range of values close to or lower than unity. D E S C R I P TI ON OF E X P E R I M E N T S
The trials were carried out in a rectangular plot 5 m wide and 15 m long, on a homogeneous 20% slope. The soil formed on a fluvial alluvium on the boundary between the fourth terrace and the present day alluvial plain of the Guadalquivir River near the town of C6rdoba. The soil has a loamy texture with 15% cobbles of about 5cm in diameter on the surface. The presence of the cobbles affects the soil loss processes but it does not change the interpretation of the results. In stoneless soils the parameters ~, B and 7 would be somewhat different. The plot was protected against external runoff at the borders by a plastic strip 20 cm high above the soil surface. A properly lined ditch at the lower border, similar to the runoff collector of Kinnell (1987), conveyed water and sediment to a tipping bucket (e.g. Barfield and Hirschi, 1986) where flow rates were measured. Sediments were sampled periodically in fixed volume bottles for subsequent laboratory determination of sediment concentrations. A sediment trap at the outlet acted as a final check of the total sediment yield. Rainfall was simulated by eight interconnected sectorial sprinklers fed by an external supply at 981 kPa pressure with self-regulators. The design was similar to that of Hart (1984). The uniformity of application was evaluated by 15 rain gauges located regularly on the plot. The Christiansen uniformity coefficient, UC (unity minus the average absolute deviation from the mean depth divided by the mean depth), was used to find the uniformity in a number. A simpler, equivalent parameter is the complementary parameter with respect to the unity of the coefficient of variation (1 - CV). In 60% of the plot the soil was covered by weeds and the stubble of herbicide-treated plants from the first trials, in June and July 1988. The second set of trials, in September 1988, were made on freshly tilled soil where all the vegetation was removed.
Results The amount of water applied to the plot in the different trials appears in Table 2 with the uniformity of distribution. In all trials but the first, the water was received with a Christiansen uniformity coefficient over 75%. The coefficient of variation remained below 0.31 in all the trials but the first. For these treatments a normal distribution could be fitted to the data as suggested by Warrick et al. (1989).
KINEMATICWAVEMODELOF SOILEROSION
77
TABLE 2 Applied rainfall characteristics Trial date
Number of rain gauges
Average rain depth (ram)
Standard deviation (mm)
UC (%)
I - CV (%)
23 June 1988 27 July 1988 28 July 1988 29 July 1988 7 September 1988 8 September 1988 9 September 1988 12 September 1988 13 September 1988
8 lI 15 15 11 15 15 15 15
16.2 23,6 7,7 8.6 94.4 12.7 14.1 14.8 14.0
7.82 6.68 2.28 2.53 18.46 2.41 2.37 4.31 4.08
57 76 77 75 84 86 89 80 79
0.48 0.70 0.69 0.71/ 0.80 0.80 0.81 0.70 0.70
Two of the runoff hydrographs are shown in Figs. 3 and 4, where the zero time refers to the beginning of runoff flow. The hydrograph of Fig. 3 corresponds to a constant rainfall rate with a concave downward rising limb, a plateau starting at 410 s and ending with a conspicuous peak when the rain stops at 595 s. Chow (1959, section 18.8) indicated the presence of similar peaks in constant rainfall experiments, when there is a sudden increase in flow after the rain stops impacting on the surface. The second hydrograph of 9 September, Fig. 4, quickly reaches the equilibrium state as a consequence of the soil's wetness after the previous rain. A summary of hydrographic data is presented in Table 3. Total runoff volume, Vr, is the result of numerical integration of hydrograph areas. The equilibrium runoff flow rate, qp, w a s taken as the average value within the steady-state time range. The average water abstraction i~, including rain water evaporation, surface retention, weeds and stubble interception, and infiltration into the soil, was derived from
a
=
2.80 rnSf2s-~ °
~
°
~°
~° o
J 5E-005
o
E
E
E
Oo
o
°
oO~
o B = 1 . 5 0 e 5 kg s ' m ~ 7 = 0.1 r n - ' k = 7 k
OE+O00
--
r
'
'
i
i
i
440
200
Time
(s)
600
i
200
Time
Fig. 3. Experimental data and model equation results for the 9 September trial.
440
i
(s)
78
A.
LAGUNA AND J.V. GIR,~LDEZ
10-
a
=
1.,33 m ' / = s
o
' o
o
oo
"~ 5 E - O O 5 -
E S
E
o o
o
° ° o o o Oo o o
oo
o v
E B = 1 5 1 e 5 kg s-* m - " y 0 1 m -~ k = 1.65 kg m - ~ / ~ s - '
OE+O00
J
i
I
i
L
~
,
i
300 Time
i
i
i
3~0
600
Time
(s)
o
'
(s)
Fig. 4. Experimental data and model equation results for the 8 September trial.
the equilibrium formula qp
=
L(r
-
(16)
ia)
where r is the average rainfall rate. Sediment concentration varied in every trial in the pattern exhibited in Figs. 3 and 4. Sediment concentration decreases little with time except for a sudden increase at the end of the rainfall, coincident with the observed runoff peak. This concentration peak might be attributable to the same origin as the flow rate peak. There is no large change in sediment concentration over the whole time range, so that a constant value may be accepted, as an approximation. This fact supports the hypothesis of Govindaraju and Kavvas (1991) that there is a constant relationship between flow and sediment discharge. To explore the relationship between sediment yields and some parameters, the correlation matrix of Table 4 was computed. Concentration time was estimated from each hydrograph. Sediment yield, S (kg), is more closely correlated with runoff volume, V~(m3), and rainfall rate, r (m s-~), respectively, than with other parameters. The linear regression between these two parameters is S =
- 0 . 9 2 5 + 18.6Vr
r 2 = 0.901
(17)
Figure 5 shows how close the fit is. This result is roughly similar to the regression expression of Singh and Chen (1982), since in their formula, the exponent of the volume was near unity. Nevertheless, their equation was derived in a larger catchment, which induced a poorer linear relation and a lower coefficient, owing to the smaller sediment yield. In erosion events with constant rainfall rate, the sediment sampling intervals could be larger without loss of precision. There is no clear relationship between sediment yield and peak runoff rate. The sediment graphs show that sediment yield increases with rain duration:
79
K I N E M A T I C WAVE M O D E L O F SOIL EROSION
I
7
I
l
q
l
~
80
A. LAGUNA
TABLE
AND
J.V. GIR,~LDEZ
4
Correlation matrix of soil loss parameters T
T
T~
-
t~
qp
1.00
.
T,
0.89
1.00
. .
.
.
t~
0.89
0.80
1.00
qp
- 0.80
- 0.52
- 0.61
1.00
r
-0.80
-0.83
-0.67
r
~
S
.
.
.
. -
-
-
-
0.40
1.00
V~
0.34
0.70
0.26
0.15
- 0.62
1.00
S
0.51
0.82
0.47
0.02
-0.73
0.96
1.00
the constancy of sediment concentration explains the linear relation between runoff volume and sediment yield. DETERMINATION
OF
EXPERIMENTAL
MODEL
PARAMETERS
FROM
THE
DATA
The shape of the experiment hydrographs is very similar to those derived from the kinematic wave model. The rising limb shows some departure from the theoretical hydrograph, which may be attributable to the influence of the neglected terms in the Saint-Venant equation for the conservation of momentum. Overton and Meadows (1976, section 4-4) discussed possible effects of different flow regimes in this region on the hydrograph shape in the second half of the rise. A similar departure is seen in the sediment graph. The constant infiltration rate assumed in the model is an approximation similar to Wooding's (1966) • index use. This approximation might not be a
ob e el e d warsoi
0.3
•
.EE 0.2
m
o
~0.1 o E o
0.0
I
0
I
I
I
I
I
I
I
1
I
2
I
I 3
I
I
I
4
Sediment yield ( k g ) Fig. 5.
Relationship between sediment yield and total runoff volume for individual simulation events.
KINEMATIC WAVE MODEL OF SOIL EROSION
81
reasonable hypothesis when the soil is initially dry. The constant rate was obtained from the rainfall rate and peak flowrate measured through eqn. (16). Parameter ~ was determined in the hydrographs according to eqn. (2), since the exponent a is assumed to be 1.5, as stated above. The concentration time, tc, was measured experimentally and inserted in the defining equation for ~. The parameters related to sediment movement were determined from the concentration of solids measured in the runoff water during the steady state and in the initial stage. Parameter B can be found in eqn. (11) from the initial sediment concentration; m is assumed to be 2 as stated above. Similarly, steady-state data allow the characterization of the relationship between 1' and k, from the expression for the sediment concentration in that period. Figures 3 and 4 show the fit for two different trials; the second trial was made on wet soil. Consequently, coefficient ~ acquired a higher value than in the first fit. A similar result obtained with parameter k, as noted above, whereas B has very similar values in both cases. All parameters fall within the ranges found in the literature and used on the sensitivity analysis. CONCLUSIONS
The experimental results demonstrate the linear relationship between the mass of soil lost and the volume of runoff. Sediment yield is more closely related to runoff volume than to other parameters such as rain intensity or concentration time. The kinematic wave model allows a close fit to experimental data; this makes it a versatile tool for the analysis of erosion trials, at least in simple cases such as those presented in this study. All the parameters were found to be within the range reported in the literature reviewed. The sensitivity of the model to the changes in the values of the three parameters studied is dependent on the values assumed by them, but in general is greater for those values which produce the highest sediment concentrations. The model is most sensitive to k, the rill erosion coefficient. REFERENCES Bagnold, R.A., 1977. Bed load transport by natural rivers. Water Resour. Res., 13: 303-311. Barfield, B.J. and Hirschi, M.C., 1986. Tipping bucket flow measurements on erosion plots. Trans. ASAE, 29: 1600-1604. Bennett, J.P., 1974. Concepts of mathematical modeling of sediment yield. Water Resour. Res., 10: 485-492. Beven, K., 1979. A sensitivity analysis of the Penman-Montheith actual evapotranspiration estimates. J. Hydrol., 44: 169-190. Blaikie, P., 1985. The Political Economy of Soil Erosion in Developing Countries. Longman, London.
82
A. LAGUNA AND J.V. GIR,~LDEZ
Blau, J.B., Woolhiser, D.A. and Lane, L.J., 1988. Identification of erosion model parameters. Trans. ASAE, 31: 859-854, 854. Chow, V.T., 1959. Open-Channel Hydraulics. McGraw-Hill, New York. Croley, II, T.E., 1982. Unsteady overland sedimentation. J. Hydrol., 56: 325-346. Croley, II, T.E. and Foster, G.R., 1984. Unsteady sedimentation in nonuniform rills. J. Hydrol., 70: 101-122. Foster, G.R., 1982. Modelling the erosion process. In: C.T. Haan, H.P. Johnson and D.L. Brakensiek (Editors), Hydrologic Modeling of Small Watersheds. ASAE Monograph, St. Joseph, MI. Govindaraju, R.S. and Kavvas, M.L., 1991. Modeling the erosion process over steep slopes: approximate analytical solutions. J. Hydrol., 127: 279-305. Govindaraju, R.S., Kavvas, M.L. and Jones, S.E., 1990. Approximate analytical solutions for overland flows. Water Resour. Res., 26: 2903-2912. Hart, G.E., 1984. Erosion from simulated rainfall on mountain range-land in Utah. J. Soil Water Cons., 39: 330-334. Huggins, L.F. and Burney, J.R., 1982. Surface runoff, storage and routing. In: C.T. Haan, H.P. Johnson and D.L. Brakensiek (Editor), Hydrologic Modeling of Small Watersheds. ASAE Monograph, St. Joseph, MI. Kinnell, P.I.A., 1987. A runoff collector and flume for use on bare fallow plots. J. Agric. Eng. Res., 38: 99-104. Nearing, M.A., Foster, G.R., Lane, L.J. and Finkner, S.C., 1989. A process-based soil erosion model for USDA-Water erosion prediction project technology. Trans. ASAE, 32: 15871593. Nearing, M.A., Deer-Ascough, L. and Laflen, J.M., 1990. Sensitivity analysis of the WEPP hillslope profile erosion model. Trans. ASAE, 33: 839-849. Overton, D.E. and Meadows, M.E., 1976. Stormwater Modeling. Academic Press, New York. Rose, C.W., 1985. Developments in soil erosion and deposition models. Adv. Soil Sci., 2: 1-63. Rose, C.W., Williams, J.R., Sander, G.C. and Barry, D.A., 1983. A mathematical model of soil erosion and deposition processes: I and II. Soil Sci. Soc. Am. J., 47: 991-1000. Shirley, E.D. and Lane, L.J., 1978. A sediment yield equation from an erosion simulation model. Hydrol. Water Resour. Ariz. Southwest, 8: 90-96. Singh, V.P. and Chen, V.J., 1982. On the relation between sediment yield and runoff volume. In: V.P. Singh (Editor), Modeling Components of Hydrologic Cycle. Water Resources, Littleton, CO. Singh, V.P. and Prasad, S.N., 1982. Explicit solutions to kinematic equations for erosion on an infiltrating plane. In: V.P. Singh (Editor), Modeling Components of Hydrologic Cycle. Water Resources, Littleton, CO, pp. 515-538. Singh, V.P. and Regl, R.R., 1983. Analytical solutions of kinematic equations for erosion on a plane. I: rainfall of indefinite duration. Adv. Water Resour., 6: 2-10. Warrick, A.W., Hart, W.E. and Yitayew, M., 1989. Calculation of distribution and efficiency for nonuniform irrigation. J. Irrig. Drain. Eng., 115: 674-686. Wooding, R.A., 1966. A hydraulic model for the catchment stream problem: III. Comparison with runoff observations. J. Hydrol., 4: 21-37. Woolhiser, D.A., 1975. Simulation of unsteady overland flow. In: K. Mahmood and V. Yevjevich (Editors), Unsteady Flow in Open Channels, Water Resources, Fort Collins, CO.
65
Elsevier Science Publishers B.V., A m s t e r d a m
[3]
The description of soil erosion through a kinematic wave model A. Laguna" and J.V. Girgldezb aDepartment of Applied Physics, University of C6rdoba, Apdo. 3048, 14080 C6rdoba, Spain bDepartment of Agronomy, University of C6rdoba, Apdo. 3048, 14080 C6rdoba, Spain (Received 21 May 1992; revision accepted 20 October 1992)
ABSTRACT Laguna, A. and Gir~ildez, J.V., 1993. The description of soil erosion through a kinematic wave model. J. Hydrol., 145: 65-82. Kinematic wave models offer a powerful tool for interpreting erosion experiments and for forecasting soil and water losses under different management systems. One important problem is the identification of the model parameters. To estimate the influence of error in the estimation of parameters a sensitivity analysis was performed. A previous survey of the soil erosion literature provided data to establish the expected ranges for each parameter. Of the three main parameters, the interrill and rill erosion coefficients, and the rill water depth coefficient, the last is the most important since sediment concentration is most sensitive to variations in this parameter, although the range of variation in the model is affected by the reference values of the other parameters. A simulated rainfall experiment in rectangular plots of 5 × 15 m 2 on a 20% slope provided data which were fitted to the model. The agreement between the model and the data was reasonably good. Sediment yield was related linearly to total runoff volume in individual simulation events.
INTRODUCTION
Soil erosion poses a serious threat to agriculture and the environment. Erosion removes the most fertile fraction of the soil, thus reducing crop productivity. Sediments plug canals and culverts, fill up reservoirs, and carry away chemicals adsorbed in soil particles, thereby polluting the environment. In Mediterranean climates, water erosion is especially severe, since after a warm dry summer there is only sparse vegetation cover to protect the soil surface from early autumn rains. Agricultural practices exacerbate the situation by enhancing the erosive effect of water. Farmers frequently plough the soil to control weeds. The topsoil becomes a loose layer on top of a layer C o r r e s p o n d e n c e to: A. Laguna, D e p a r t m e n t of Applied Physics, University o f C 6 r d o b a , Apdo. 3048, 14080 C 6 r d o b a , Spain.
0022-1694/93/$06.00
© 1993 - - Elsevier Science Publishers B.V. All rights reserved
66
A. L A G U N A A N D ,I.V. G I R A L D E Z
NOTATION a
b B C
C
El & g h i ia
k K
Ke Kr L m n
q, qe qp r
S
s, SL !
tc le /r 0 lsc
T
L Tr Vc
v~ x X,
exponent exponent interrill erosion coefficient [M L -m-2 T ~- D] sediment concentration [M L -3] cover factor from the USLE interrill erosion [M L 2 T ~] rill erosion[M L 2 T 1] acceleration of gravity [L T 2] flow depth [L] infiltration rate [L T ~] average water abstraction [L T -t] rill erosion coefficient [M L i b T i] soil erodibility factor from the USLE coefficient [T 2] coefficient [L ~ T] slope length [L] exponent Manning's coefficient water discharge rate per unit width [L 2 T -~] equilibrium runoff flow rate per unit width [L 2 T ~] rainfall intensity [L T i] total sediment yield [M] sensitivity of c to P~ parameter plane slope time [T] water concentration time [T] runoff equilibrium time [T] runoff initiation time [T] sediment concentration time [T] rainfall duration [T] transport capacity [M L- t T - t ] runoff duration [T] settling velocity [L T 1] total runoff volume [L 3] distance along the flow direction [L] threshold of slope erosion [L]
Greek letters
depth-discharge coefficient [L 2 a T q ] rill erosion coefficient [L i] O. Oo initial and final soil moisture [M M ~] coefficient [M L -2 T -~] water density [M L -3] P critical shear stress of the flow [M L ~ T 2] •c shear stress of the flow [M L 1 T-2] f~ stream power [M T -3] threshold stream power [M T 3] f~o Ot
7
K I N E M A T I C WAVE M O D E L O F SOIL EROSION
67
of hard plough sole penetrable only with great difficulty by water and plant roots. When water filters into the soil, buoyancy forces raise the topsoil above the plough sole and carry it downslope. The resulting gullies are thenceforth U-shaped. Nevertheless, agriculture is not the only participant responsible for soil erosion. Forestry practices require the establishment of roads and fire breaks along mountain slopes, practices which require clearing the vegetation. Civil works such as the opening of new roads or simply urbanizing vegetated areas remove plant protection from the soils, thus accelerating natural erosion. In spite of soil erosion being a menace for agriculture and the environment, the number of studies in many countries is surprisingly low. Blaikie (1985, section 2.2) remarked on this deficiency, especially in less-developed countries, which are severely affected by this problem. There is general concern with respect to the effects of erosion, or to its related problem of land degradation and desert encroachment or even desertification, but in most of the cases the real extent of the problem is not known precisely. The use of an erosion model might be very helpful for predicting the consequences of soil management practices, and land-use changes. The complexity of the phenomena makes the formulation of a simple model difficult, even when a single factor like rain splash of water runoff is the only agent of erosion. Many researchers prefer to use empirical relations between observed soil loss rates or sediment gains with selected soil or environmental parameters, via a statistical equation. Although this approach is convenient for a rough estimate, its use is not recommended in conditions differing from those considered originally. An appropriate way of studying soil erosion processes is through the formulation of the fundamental transport equations of water and sediment (e.g. Bennett, 1974, or Rose, 1985). Once posed, these equations can be solved either by numerical or analytical methods. The kinematic wave simplification for surface water flow on a plane yields an analytical solution for the cases of constant rainfall and infiltration rates, as Singh and Prasad (1982), among others, have demonstrated. Other authors, such as Govindaraju and Kavvas (1991), have worked on an approximate solution of a diffusion wave model previously derived by Govindaraju et al. (1990). This is a more complete hydrological model since it allows for a varying time infiltration rate in discrete time steps, and assumes a constant sediment concentration for the whole erosion processs. In spite of being very useful from the hydrological point of view, the assumption of constancy of sediment concentration is a shortcoming of this model. The kinematic wave model may be applied to the interpretation of erosion experiments after their parameters have been identified. The purpose of this report is to explore the fit of the kinematic wave
68
A LAGUNA AND J.V GIRALDEZ
model to experimental results under controlled conditions and to analyse its sensitivity to the most relevant parameters. T H E K I N E M A T I C SOIL E R O S I O N M O D E L
The model based on the kinematic wave approach to the Saint-Venant equations consists of three expressions. Expression (1) is the mass conservation expression for water 0h 0q a--t + 0x
r -- i
(1)
where h is the water depth [L]; q is the flow rate per unit width [L 2 T - ~]; r and i are the rainfall and infiltration rates respectively [L T - ~]; t refers to time [T], and x to space [L]. The relationship between depth and flow rate constitutes the m o m e n t u m conservation equation (expression 2): q =
~" h ~
(2)
where ~ is a coefficient related to Manning's n or Chezy's C, [L 2-" T ~], and a is an exponent. The mass conservation equation for sediment furnishes the last expression
(3): a(hc)
c~----~+
a(qc)
8x
-
El + ER
(3a)
where c is the sediment concentration in water flow [M L-3]. El and ER respectively represent the sediment production rate due to interrill and rill erosion processes [M L -2 T-~]. There are several proposals for the El and E R terms. Singh and Prasad (1982) selected a potential expression o f El with respect to the rainfall rate: E,
=
(3b)
B . rm
where B is a coefficient and m an exponent. Rill erosion depends on the water depth in the rill, and on the flow rate which will induce deposition if it is not too high: ER
---
?(k" h b -
c " q)
(3c)
where 7 is another parameter and b an exponent which according to Shirley and Lane's (1978) findings is close to a. Therefore it would be assumed that b=a.
Other approaches to the expression of sediment generation by the interaction of erosion-deposition processes have been proposed by Rose et al. (1983)
69
KINEMATIC WAVE MODEL OF SOIL EROSION
and Croley and Foster (1984). Croley and Foster (1984) proposed a similar equation to (3c) which was later modified by Nearing et al. (1989) in the U S D A - W E P P framework. Their expression is ER =
3" ( 1
(3d)
qTc c )
where Tc is the transport capacity of rill flow and ~ is a coefficient whose value is either Kr" (rr - r~) when the term in parentheses in (3d) is positive, or (v c • T~)/q when it is negative. The settling velocity is Vc; rf represents the shear stress of the flow, and its critical value is r~. Transport capacity is related to the flow's shear stress, -of, through a potential equation. Alternatively Rose et al. (1983) separated the distinct contributions of the rainfall detachment rate, essentially (3b); they use the sediment entrainment rate, and sediment deposition rate, Vc • c. Bagnold stream power is defined as the shear stress by the mean flow velocity, 92 = p • g . SL " qL " x / L (e.g. Bagnold, 1977). In that expression p is water density, qL is the flow rate at the end of the plane, and SL its slope. This equation was adopted as the parameter controlling the net bed load transport rate, over a certain threshold value 920, defined as the stream power for the distance at which there was no sediment entrainment, x , . The final expression for ER is ER
=
p " g " SL " K~ "
+
qL
1 -- - -
X
+
(hc)
(3e)
The inclusion of the partial derivative of the product of water depth by sediment concentration with respect to time allows the elimination of such a term in eqn. (3a). This converts a partial differential equation into an ordinary differential equation of simpler solution. Time-dependent parameters are evaluated in steps. Different size sediments are introduced once their mass fraction is known. For the purposes of the present work the expression (3c) is preferred since its integration is straightforward. The initial conditions and boundary conditions are as follows, where T is rainfall duration, and L is the length of the slope: =
O,
0 <<. x <<. L
(4a)
h(0, t) =
0,
t >~ 0
(4b)
h(x,O)
The rainfall is assumed to be homogeneously distributed over the whole surface, a plane, with constant rate: r(x,t)
=
r > O,
r(x, t)
=
O,
0 < t < T t >1 T
(5a) (5b)
70
A. LAGUNA AND J.V. GIR,~LDEZ
Water infiltrates the soil at a constant rate, whenever the water depth over the soil is positive: i(x, t) =
i > O,
i(x, t) =
O,
(6a)
h(x, t) > 0 h(x, t) =
0
(6b)
Note that r must be greater than i. Solutions o f the equations
The solution of eqns. (1)-(3) subjected to conditions (4)-(6) may be found through the m e t h o d of characteristics. The solutions are found for the pair of values (x, t), whose paths are the characteristics. The times at which the characteristics, starting at the origin of coordinates in the t - x plane, reach the downstream boundary are k n o w n as the concentration times, to, and tsc, for water and sediments respectively. Singh and Prasad (1982) developed solutions for the different domains b o u n d e d by the limiting characteristics of water and sediment, the recession curve and the lines t = T, t = 0, x = 0, and x = L. There are a few misprints; in their eqn. (60), the characteristic for the sediments in domain D3b, is corrected to t -
T
=
r(a -
]) + i
"{[(r - i)xx] '/" - (r'x'o -
i ' x ) '/a}
(7)
and the value of the auxiliary function #(xT, T) used in their eqn. (64) for the sediment concentration in domain D3a should be /~(xx, T)
=
7~(r-
i) Ca 1)[ ~ ( r -
i)¢' ~)la(xw/~)'l~-- T ] a -- 7
(8)
M o d e l parameters
The solution of eqns. (1)-(3) will characterize completely the soil erosion processes for any point in the plane at all times. The main problem in applying the model is the identification of the parameters, since they are not immediately measured at the experimental site. In this section the pertinent literature will be reviewed, looking for the reported range of variation as a preliminary step in the sensitivity analysis of the model. There are six parameters in all: the ~, B, 7 and k coefficients, and the a and m exponents. Parameters o f the dynamic equation Parameters ~ and a determine surface water movement. Under a laminar
K I N E M A T I C WAVE M O D E L O F SOIL EROSION
71
TABLE 1 Selected values of the coefficient ~ for some slopes and surfaces Surface
Manning's n
Bare clay-loam soil (eroded)
0.012-0.033
Sparse vegetation cover
0.053-0.130
Short grass prairie
0.100-0.200
Bluegrass Sod
0.170-0.480
Slope
(%)
(m~/:s ' )
l0 20 10 20 10 20 10 20
26.35-9.58 37.27-13.55 5.97-2.43 8.44-3.44 3.16-1.58 4.47-2.24 1.86-0.66 2.63-0.93
is computed from eqn. (9), with Manning's n values from Woolhiser (1975).
flow regime, a assumes a value of 3 in the Darcy-Weisbach formula, whereas under a turbulent regime, values of 5/3 in Manning's equation, and 1.5 in Ch6zy's equation have been proposed. However, Croley and Foster (1984) estimated values of a close to 0.25 and argued the convenience of the parameter being smaller than unity. These values were computed from the estimates of the time of concentration and the average velocity of water flow in the rills of experimental erosion plots. Croley and Foster justified their estimates on the basis of the irregularity of canal geometry. Nevertheless, values for a that are smaller than unity imply multiple values of water depth at a given time in the recession stage. Huggins and Burney (1982) supported Manning's exponent. Singh and Regl (1983) proposed 1.5 on the grounds of its commonly accepted usage. In this study, a 1.5 value was used for a. Coefficient ~ has a wider range than a. This coefficient depends basically on the slope and on surface roughness. Manning's equation is
(&),/2 -
n
(9)
where SL is the slope and n is Manning's coefficient. Some possible values for ~, for different slopes and surfaces are gathered in Table 1. They have been computed from eqn. (9) with Manning's n values from Woolhiser (1975). Reasonable values for ~ are between 0.5 and 40 m ~/2s -1 . Figure 1 shows the concentration of sediments in the model for several ~ values within the given interval, keeping the other parameters constant. The selected values for rainfall and infiltration rate were somewhat high, since one of the objectives of the experiment was the production of high erosion events. The main effect of the variation of ~ is in time of concentration. The higher ~ is (this implies
A. L A G U N A AND J.V. G I R A L D E Z
72 50
75 B = 5.05e05 k9 s-'m 7 =0.1 m -~
•
= 10 m i ~ s-' 7 = 0.1 m - ' k = 50 kg m -=,a S-'
\
k = 50 kg m-~I~ s - '
'~kt)
B (kg s
5O
E ~25
E
m-4 )
7.525e5 5.05e5
v
I .Oe4
u 25
40
.
.
.
3~0
.
.
Time
50j |
a =
.
.
6~0
.
,
,
,
(s)
,
310~
50 |
10 m ~r~ s -I
kB = 5 . 0 5 e 0 5 kg s - ' en-° 50 kg m-=/z s - '
(m-') .~.
0.055 0.1
~25
_~ 25
0.3 0.5
v
t
'
,
,
6O0
cr = 10 rn '/= s -~ B = 5 . 0 5 e 0 5 kg s - ' rn ~ 7 = 0.1 m - '
k (kg m "fC*s-')
{il,o2~o 75 50
(5) 0
~
(s)
Time
I[4)r ~ '
'
'
'
0
I 300
'
'
r
0
, 600
,
,
0
Time (s)
Fig.
I. S e d i m e n t
r =
1.376e-5ms -t , L =
concentration
,
t
I 300
Time
for s e v e r a l c(, B, 7 a n d k v a l u e s (n =
'
f
, ,"'T'--~ ' 600
(s)
1.5, m = 2, i = 9 . 6 5 e - 6 m s
~,
1 0 0 m , T = 5 0 0 s ) . T i m e f r o m the start o f rainfall a n d runoff.
steep slope and smooth surface), the shorter is the time of concentration. For the lowest c( value the concentration time is larger than rainfall duration; in that case the steady state is not reached. Sediment concentrations in the steady state are surprisingly high for small values. Since the peak flow rate is independent of this parameter, the sediment load [M T - 1 L- 1] will also be high in the case of small ~. One can expect that for small ~ values, surface flow velocity is low and so, consequently, is the erosive capacity of the flow. A possible explanation for this paradox is that the k coefficient in the transport capacity of the flow (eqn. (3c)) must be linked to ~. Therefore if ~ changes, k should change as well, to a higher rate. After a certain threshold value, there are no great variations among the curves. These differences are greater the smaller 0( is. If the value of ~ is very small the curves increase with time.
Parameters of the interrill erosion equation Parameters B and m influence interrill or rain splash erosion. Singh and Regl (1983) and Blau et al. (1988) proposed a value o f m = 1. Croley (1982)
KINEMATIC
WAVE MODEL
73
OF SOIL EROSION
referred to Foster's (1982) preference for m = 2 based on a more realistic interpretation of splash experiments. The latter author suggested the experimental expression E l ( k g m - 2 h i) =
O.O138.Kl(kghN-im-2).C.r2(mmh
i)
(10)
where Kj is the soil erodibility factor for raindrop impact. In the absence of K, values, the universal soil loss equation K factor may be adopted, even though this factor assumes both interrill and intrarill erosion processes. The C factor is the cover factor from the universal soil loss equation. The parameter B may be evaluated experimentally or estimated through expressions similar to eqn. (10), for further refinement after calibration experiments. Croley and Foster (1984) obtained values for B between 1 x 106 and 1.5 x 106kg s m -2 after examining the results of an experimental plot under simulated rainfall. From the results of the model, for the time zero c =
B.r m : r
--
(11)
l
This equation implies that the initial concentration of sediments is proportional to B. After eqn. (10) a reasonable range for the values of B is (1 x 104 to 1 x 106kgs Im-4), which corresponds to the pair of values for K a n d C of (0.02, 0.01) and (0.04, 0.6) respectively.
Parameters of the intrarill erosion equation Parameters Y and k are associated with the intrarill erosion process. According to eqn. (3c) the sediment yield is proportional to the difference between transport capacity and sediment load. Parameter y [L ~] assumes values in the range 3-33 m - i in the work of Foster and Huggins, and David as quoted by Foster (1982). For y values higher than unity most of the resulting sediment concentrations are negligible. Singh and Regl (1983) suggested a reasonable value of 0.03 m -I . Parameter k, [M L -(a+l) T-l], is related to rill erodibility. Singh and Regl (1983) proposed a value of 1.5 kg s ' m -25. I f k is smaller than unity there are no differences in the results for different values of the parameter. In the rest of the work, values of a = 1.5 and m = 2 will be assumed. As was stated previously in eqn. (11) and as Fig. 1 shows, B is proportional to the concentration of sediments at time zero. In the steady state, holding the rest of the parameters constant, both sediment concentration and the value of parameter B follow a straight line with positive slope. The same is true of parameter k. Nevertheless y is inversely related to sediment concentration in the steady state. If this parameter assumes small values, both net soil loss and
74
A. LAGUNA AND J.V. GIRALDEZ
sediment deposition are relatively small. In such a case, changes from the initial conditions take place slowly (Blau et al., 1988). Sediment concentration increases monotonically with time from initial conditions until the steady state if B'~
k
r~ - -
<
r-i
1
(12)
since the whole phenomenon is controlled by rill erosion, decreasing in the negative case. ANALYSIS OF THE MODEL SENSITIVITY
The determination of the model parameters is difficult because most of them are not directly measurable. Therefore these parameters need to be estimated through indirect measurements in erosion experiments. In these cases a sensitivity analysis may ascertain the consequences to the final result of a small change in the value of the parameter. Since the final result in the present work is the sediment concentration c = f(PI, P2,...,Pn)
(13)
where PI, P2,. • • P, are the parameters and input variables for the model. The sensitivity of c to any one of the P~ is (Beven, 1979) ~C
Pi
. . c Si .- . 6p
(14)
or the equivalent S~ =
Ac _Pi AP" c
(15)
representing the fractional change in P transmitted to c, which affords a comparison of sensitivity with parameters whose range is of different order of magnitude. One shortcoming of this method is the assumption of parameter independence. Among the five model parameters, three were selected for the sensitivity analysis since it was assumed a fixed value of 1.5 for the exponent and the coefficient ~ is easily obtained from eqn. (9). Sensitivity to one parameter may change as the reference values of the other parameters change. Nearing et al. (1990) computed the sensitivity of the model to the parameters by applying the extreme values of the respective ranges. Therefore the sensitivity of the sediment concentration with respect to parameters B, 7 and k was computed for the maximum, minimum, and intermediate values found in the literature and indicated above as a function of time, as shown in Fig. 2. The
KINEMATIC WAVE MODEL OF SOIL EROSION
75
1.5 7=0.1 k = 50
rn-, kg rn-vz
s-' • B
~.0
• ~
03
•
®
=
le06
k~q
a B
5.05e0~
0 (3 =
le04
s -' m-*
o
>,
O
o
>
o
e
®
o
® o
05
o
m o
wl
a
0 0
O0
:
i
o
e
e
e
®
e
e
Q
o
o
o
~
o
o
o
® o
*°°ooo,ooooo,~ i
i
I
i
]
]
T
600
300 Time
(s)
].5-
15 B k
= =
505e05 kg 50 kg m-$/,
s-' s-'
B=§OSeO$kgs-,m 7 = 0 1 m - *
rn-,
-10
10
z>
• 7
=
0.5
a 7
=
0 "y
:
01 001
r n -~
•
k
o ¢
k= k =
=
100
kg
m-s/~
~@ • ® ® o o ® m e m ® o
8o00°°°**0oo o
e
o e
00
e
®
o
o
o
o
o
o
o
e
0
0.0
i
o
,
~
,
[
300 Time
(s)
i
0
o o o o o o o
@ o a
O.S u3
®
@
a
o
o ~^00
°
i
o o
°
o e e e e e o ® e m
i
o
o °
a
s-~
50 I ®o
-0.5
~
=l
0.0
'
'
l
0000000000 o° l
l
r
l
600 Time
l
l
l
I
300 (s)
Fig. 2. Sensitivity coefficients of sediment concentration to B, ;, a n d k p a r a m e t e r s = 10mW2s ~ , i = 9 . 6 5 e - 6 m s - ] ~ r = 1 . 3 7 6 e - 5 m s i L = 100m).
600
(n =
1.5, m = 2,
values of rainfall and infiltration rate adopted were, as in Fig. 1, obtained in the experiments described below in the erosion plots. In the curves shown in Fig. 2 the variation over time in the sensitivity coefficient to any parameter can be seen. In the rising stage, erosion processes are controlled by the detachment of soil particles by raindrop impact. At this stage there is a maximum value in sensitivity to the interrill erosion parameter B, although the sensitivity to the other parameters is low. However, in the recession stage, when soil loss is primarily due to shear by water flow, the sensitivity with respect to B is smaller than during other stages, and the sensitivity to k acquires its maximum value. The highest values for the sensitivity coefficient are found for the parameter values that imply greater sediment concentrations. The sensitivity coefficient for parameter 7 is negative since its value increases when the value of c diminishes. These results agree with the sensitivity analysis of Blau et al. (1988): the higher the value of parameter 7, the lower the sensitivity of the model. In the steady state where the greater fraction of soil loss occurs, the most important parameter in terms
76
A. LAGUNA AND J.V. GIR,~LDEZ
of model sensitivity appears to be k, except for the range of values close to or lower than unity. D E S C R I P TI ON OF E X P E R I M E N T S
The trials were carried out in a rectangular plot 5 m wide and 15 m long, on a homogeneous 20% slope. The soil formed on a fluvial alluvium on the boundary between the fourth terrace and the present day alluvial plain of the Guadalquivir River near the town of C6rdoba. The soil has a loamy texture with 15% cobbles of about 5cm in diameter on the surface. The presence of the cobbles affects the soil loss processes but it does not change the interpretation of the results. In stoneless soils the parameters ~, B and 7 would be somewhat different. The plot was protected against external runoff at the borders by a plastic strip 20 cm high above the soil surface. A properly lined ditch at the lower border, similar to the runoff collector of Kinnell (1987), conveyed water and sediment to a tipping bucket (e.g. Barfield and Hirschi, 1986) where flow rates were measured. Sediments were sampled periodically in fixed volume bottles for subsequent laboratory determination of sediment concentrations. A sediment trap at the outlet acted as a final check of the total sediment yield. Rainfall was simulated by eight interconnected sectorial sprinklers fed by an external supply at 981 kPa pressure with self-regulators. The design was similar to that of Hart (1984). The uniformity of application was evaluated by 15 rain gauges located regularly on the plot. The Christiansen uniformity coefficient, UC (unity minus the average absolute deviation from the mean depth divided by the mean depth), was used to find the uniformity in a number. A simpler, equivalent parameter is the complementary parameter with respect to the unity of the coefficient of variation (1 - CV). In 60% of the plot the soil was covered by weeds and the stubble of herbicide-treated plants from the first trials, in June and July 1988. The second set of trials, in September 1988, were made on freshly tilled soil where all the vegetation was removed.
Results The amount of water applied to the plot in the different trials appears in Table 2 with the uniformity of distribution. In all trials but the first, the water was received with a Christiansen uniformity coefficient over 75%. The coefficient of variation remained below 0.31 in all the trials but the first. For these treatments a normal distribution could be fitted to the data as suggested by Warrick et al. (1989).
KINEMATICWAVEMODELOF SOILEROSION
77
TABLE 2 Applied rainfall characteristics Trial date
Number of rain gauges
Average rain depth (ram)
Standard deviation (mm)
UC (%)
I - CV (%)
23 June 1988 27 July 1988 28 July 1988 29 July 1988 7 September 1988 8 September 1988 9 September 1988 12 September 1988 13 September 1988
8 lI 15 15 11 15 15 15 15
16.2 23,6 7,7 8.6 94.4 12.7 14.1 14.8 14.0
7.82 6.68 2.28 2.53 18.46 2.41 2.37 4.31 4.08
57 76 77 75 84 86 89 80 79
0.48 0.70 0.69 0.71/ 0.80 0.80 0.81 0.70 0.70
Two of the runoff hydrographs are shown in Figs. 3 and 4, where the zero time refers to the beginning of runoff flow. The hydrograph of Fig. 3 corresponds to a constant rainfall rate with a concave downward rising limb, a plateau starting at 410 s and ending with a conspicuous peak when the rain stops at 595 s. Chow (1959, section 18.8) indicated the presence of similar peaks in constant rainfall experiments, when there is a sudden increase in flow after the rain stops impacting on the surface. The second hydrograph of 9 September, Fig. 4, quickly reaches the equilibrium state as a consequence of the soil's wetness after the previous rain. A summary of hydrographic data is presented in Table 3. Total runoff volume, Vr, is the result of numerical integration of hydrograph areas. The equilibrium runoff flow rate, qp, w a s taken as the average value within the steady-state time range. The average water abstraction i~, including rain water evaporation, surface retention, weeds and stubble interception, and infiltration into the soil, was derived from
a
=
2.80 rnSf2s-~ °
~
°
~°
~° o
J 5E-005
o
E
E
E
Oo
o
°
oO~
o B = 1 . 5 0 e 5 kg s ' m ~ 7 = 0.1 r n - ' k = 7 k
OE+O00
--
r
'
'
i
i
i
440
200
Time
(s)
600
i
200
Time
Fig. 3. Experimental data and model equation results for the 9 September trial.
440
i
(s)
78
A.
LAGUNA AND J.V. GIR,~LDEZ
10-
a
=
1.,33 m ' / = s
o
' o
o
oo
"~ 5 E - O O 5 -
E S
E
o o
o
° ° o o o Oo o o
oo
o v
E B = 1 5 1 e 5 kg s-* m - " y 0 1 m -~ k = 1.65 kg m - ~ / ~ s - '
OE+O00
J
i
I
i
L
~
,
i
300 Time
i
i
i
3~0
600
Time
(s)
o
'
(s)
Fig. 4. Experimental data and model equation results for the 8 September trial.
the equilibrium formula qp
=
L(r
-
(16)
ia)
where r is the average rainfall rate. Sediment concentration varied in every trial in the pattern exhibited in Figs. 3 and 4. Sediment concentration decreases little with time except for a sudden increase at the end of the rainfall, coincident with the observed runoff peak. This concentration peak might be attributable to the same origin as the flow rate peak. There is no large change in sediment concentration over the whole time range, so that a constant value may be accepted, as an approximation. This fact supports the hypothesis of Govindaraju and Kavvas (1991) that there is a constant relationship between flow and sediment discharge. To explore the relationship between sediment yields and some parameters, the correlation matrix of Table 4 was computed. Concentration time was estimated from each hydrograph. Sediment yield, S (kg), is more closely correlated with runoff volume, V~(m3), and rainfall rate, r (m s-~), respectively, than with other parameters. The linear regression between these two parameters is S =
- 0 . 9 2 5 + 18.6Vr
r 2 = 0.901
(17)
Figure 5 shows how close the fit is. This result is roughly similar to the regression expression of Singh and Chen (1982), since in their formula, the exponent of the volume was near unity. Nevertheless, their equation was derived in a larger catchment, which induced a poorer linear relation and a lower coefficient, owing to the smaller sediment yield. In erosion events with constant rainfall rate, the sediment sampling intervals could be larger without loss of precision. There is no clear relationship between sediment yield and peak runoff rate. The sediment graphs show that sediment yield increases with rain duration:
79
K I N E M A T I C WAVE M O D E L O F SOIL EROSION
I
7
I
l
q
l
~
80
A. LAGUNA
TABLE
AND
J.V. GIR,~LDEZ
4
Correlation matrix of soil loss parameters T
T
T~
-
t~
qp
1.00
.
T,
0.89
1.00
. .
.
.
t~
0.89
0.80
1.00
qp
- 0.80
- 0.52
- 0.61
1.00
r
-0.80
-0.83
-0.67
r
~
S
.
.
.
. -
-
-
-
0.40
1.00
V~
0.34
0.70
0.26
0.15
- 0.62
1.00
S
0.51
0.82
0.47
0.02
-0.73
0.96
1.00
the constancy of sediment concentration explains the linear relation between runoff volume and sediment yield. DETERMINATION
OF
EXPERIMENTAL
MODEL
PARAMETERS
FROM
THE
DATA
The shape of the experiment hydrographs is very similar to those derived from the kinematic wave model. The rising limb shows some departure from the theoretical hydrograph, which may be attributable to the influence of the neglected terms in the Saint-Venant equation for the conservation of momentum. Overton and Meadows (1976, section 4-4) discussed possible effects of different flow regimes in this region on the hydrograph shape in the second half of the rise. A similar departure is seen in the sediment graph. The constant infiltration rate assumed in the model is an approximation similar to Wooding's (1966) • index use. This approximation might not be a
ob e el e d warsoi
0.3
•
.EE 0.2
m
o
~0.1 o E o
0.0
I
0
I
I
I
I
I
I
I
1
I
2
I
I 3
I
I
I
4
Sediment yield ( k g ) Fig. 5.
Relationship between sediment yield and total runoff volume for individual simulation events.
KINEMATIC WAVE MODEL OF SOIL EROSION
81
reasonable hypothesis when the soil is initially dry. The constant rate was obtained from the rainfall rate and peak flowrate measured through eqn. (16). Parameter ~ was determined in the hydrographs according to eqn. (2), since the exponent a is assumed to be 1.5, as stated above. The concentration time, tc, was measured experimentally and inserted in the defining equation for ~. The parameters related to sediment movement were determined from the concentration of solids measured in the runoff water during the steady state and in the initial stage. Parameter B can be found in eqn. (11) from the initial sediment concentration; m is assumed to be 2 as stated above. Similarly, steady-state data allow the characterization of the relationship between 1' and k, from the expression for the sediment concentration in that period. Figures 3 and 4 show the fit for two different trials; the second trial was made on wet soil. Consequently, coefficient ~ acquired a higher value than in the first fit. A similar result obtained with parameter k, as noted above, whereas B has very similar values in both cases. All parameters fall within the ranges found in the literature and used on the sensitivity analysis. CONCLUSIONS
The experimental results demonstrate the linear relationship between the mass of soil lost and the volume of runoff. Sediment yield is more closely related to runoff volume than to other parameters such as rain intensity or concentration time. The kinematic wave model allows a close fit to experimental data; this makes it a versatile tool for the analysis of erosion trials, at least in simple cases such as those presented in this study. All the parameters were found to be within the range reported in the literature reviewed. The sensitivity of the model to the changes in the values of the three parameters studied is dependent on the values assumed by them, but in general is greater for those values which produce the highest sediment concentrations. The model is most sensitive to k, the rill erosion coefficient. REFERENCES Bagnold, R.A., 1977. Bed load transport by natural rivers. Water Resour. Res., 13: 303-311. Barfield, B.J. and Hirschi, M.C., 1986. Tipping bucket flow measurements on erosion plots. Trans. ASAE, 29: 1600-1604. Bennett, J.P., 1974. Concepts of mathematical modeling of sediment yield. Water Resour. Res., 10: 485-492. Beven, K., 1979. A sensitivity analysis of the Penman-Montheith actual evapotranspiration estimates. J. Hydrol., 44: 169-190. Blaikie, P., 1985. The Political Economy of Soil Erosion in Developing Countries. Longman, London.
82
A. LAGUNA AND J.V. GIR,~LDEZ
Blau, J.B., Woolhiser, D.A. and Lane, L.J., 1988. Identification of erosion model parameters. Trans. ASAE, 31: 859-854, 854. Chow, V.T., 1959. Open-Channel Hydraulics. McGraw-Hill, New York. Croley, II, T.E., 1982. Unsteady overland sedimentation. J. Hydrol., 56: 325-346. Croley, II, T.E. and Foster, G.R., 1984. Unsteady sedimentation in nonuniform rills. J. Hydrol., 70: 101-122. Foster, G.R., 1982. Modelling the erosion process. In: C.T. Haan, H.P. Johnson and D.L. Brakensiek (Editors), Hydrologic Modeling of Small Watersheds. ASAE Monograph, St. Joseph, MI. Govindaraju, R.S. and Kavvas, M.L., 1991. Modeling the erosion process over steep slopes: approximate analytical solutions. J. Hydrol., 127: 279-305. Govindaraju, R.S., Kavvas, M.L. and Jones, S.E., 1990. Approximate analytical solutions for overland flows. Water Resour. Res., 26: 2903-2912. Hart, G.E., 1984. Erosion from simulated rainfall on mountain range-land in Utah. J. Soil Water Cons., 39: 330-334. Huggins, L.F. and Burney, J.R., 1982. Surface runoff, storage and routing. In: C.T. Haan, H.P. Johnson and D.L. Brakensiek (Editor), Hydrologic Modeling of Small Watersheds. ASAE Monograph, St. Joseph, MI. Kinnell, P.I.A., 1987. A runoff collector and flume for use on bare fallow plots. J. Agric. Eng. Res., 38: 99-104. Nearing, M.A., Foster, G.R., Lane, L.J. and Finkner, S.C., 1989. A process-based soil erosion model for USDA-Water erosion prediction project technology. Trans. ASAE, 32: 15871593. Nearing, M.A., Deer-Ascough, L. and Laflen, J.M., 1990. Sensitivity analysis of the WEPP hillslope profile erosion model. Trans. ASAE, 33: 839-849. Overton, D.E. and Meadows, M.E., 1976. Stormwater Modeling. Academic Press, New York. Rose, C.W., 1985. Developments in soil erosion and deposition models. Adv. Soil Sci., 2: 1-63. Rose, C.W., Williams, J.R., Sander, G.C. and Barry, D.A., 1983. A mathematical model of soil erosion and deposition processes: I and II. Soil Sci. Soc. Am. J., 47: 991-1000. Shirley, E.D. and Lane, L.J., 1978. A sediment yield equation from an erosion simulation model. Hydrol. Water Resour. Ariz. Southwest, 8: 90-96. Singh, V.P. and Chen, V.J., 1982. On the relation between sediment yield and runoff volume. In: V.P. Singh (Editor), Modeling Components of Hydrologic Cycle. Water Resources, Littleton, CO. Singh, V.P. and Prasad, S.N., 1982. Explicit solutions to kinematic equations for erosion on an infiltrating plane. In: V.P. Singh (Editor), Modeling Components of Hydrologic Cycle. Water Resources, Littleton, CO, pp. 515-538. Singh, V.P. and Regl, R.R., 1983. Analytical solutions of kinematic equations for erosion on a plane. I: rainfall of indefinite duration. Adv. Water Resour., 6: 2-10. Warrick, A.W., Hart, W.E. and Yitayew, M., 1989. Calculation of distribution and efficiency for nonuniform irrigation. J. Irrig. Drain. Eng., 115: 674-686. Wooding, R.A., 1966. A hydraulic model for the catchment stream problem: III. Comparison with runoff observations. J. Hydrol., 4: 21-37. Woolhiser, D.A., 1975. Simulation of unsteady overland flow. In: K. Mahmood and V. Yevjevich (Editors), Unsteady Flow in Open Channels, Water Resources, Fort Collins, CO.