- Email: [email protected]
Unsteady soil erosion model, analytical solutions and comparison with experimental results
_• ELSEVIER
Journal of
Hydrology Journal of Hydrology 178 (1996) 351-367
Unsteady soil erosion model, analytical solutions and comparison with experimental results G.C. Sander a'*, P.B. Hairsine b, C.W. Rose c, D. Cassidy a, J.-Y. Parlange d, W.L. Hogarth c, I.G. L i s l e c aFaculty of Science and Technology, Griffith University, Nathan, Qld. 4111, Australia bCSIRO Division of Soils, P.O. Box 639, Canberra ACT, Australia CFaculty of Environmental Sciences, Griffith University, Nathan, Qld. 4111, Australia dDepartment of Agricultural and Biological Engineering, Cornell University, Ithaca, N Y 14853, USA Received 9 March 1995; accepted 31 March 1995
Abstract Hairsine and Rose developed a soil erosion model which described the erosion transport of the multiparticle sizes in sediment for rain-impacted flows in the absence of entrainment in overland flow. In this paper we extend their steady-state solutions to account for the time variation of suspended sediment concentration during an erosion event. A very simple approximate analytical solution is found which agrees extremely well with experimental data obtained from nine experiments. We are able to reproduce the rapid initial increase to a peak in the total sediment concentration, which occurs about 3-5 min after the commencement of rainfall, as well as the subsequent declining exponential tail towards steady-state conditions. We are also able to show that the fraction of shielding of the original soil bed resulting from depositing sediment reaches its equilibrium value on about the same time-scale as the total peak suspended sediment concentration. Interestingly, we find that the masses of the individual particles which form this deposited layer are far from equilibrium, and that there is a great deal of continuous reworking and sorting of this material during the erosion event. Finally, our solution shows that the initial peak in the total sediment concentration is due to the enrichment of this sediment by the finer size classes and that as the event continues their percentage contribution diminishes.
1. Introduction I n a recent p a p e r b y H a i r s i n e a n d R o s e (1991) a soil e r o s i o n m o d e l h a s been d e v e l o p e d to d e s c r i b e the e r o s i o n a n d t r a n s p o r t o f s e d i m e n t in r a i n - i m p a c t e d flows * Corresponding author. 0022-1694/96/$15.00 © 1996 - Elsevier Science B.V. All rights reserved SSDI 0022-1694(95)02810-2
352
G.C. Sander et al. / Journal of Hydrology 178 (1996) 351-367
in the absence of entrainment by shallow overland flow. Soils are naturally composed of a range of particle sizes and aggregates, and in the past, soil erosion models have tended to lump or average over all size classes, resulting in a single conservation equation for the total suspended sediment concentration, c, as a function of distance x and time t (Govindaraju and Kavvas, 1991; Laguna and Giraldez, 1993). In contrast, the approach of Hairsine and Rose (1991) is to consider the contributions of the individual size classes to the total suspended sediment concentration separately. This leads to separate coupled conservation equations for the suspended sediment concentration ci(x, t) for each particle size class i. While this does increase the number of equations and complicates their solution, there are sound physical grounds for justifying this approach. Once particles have become suspended, they will naturally begin to settle back towards the soil surface owing to their immersed weight. The fine particles, having relatively small settling velocities, can be carried significant distances in the direction of the surface water flow before depositing back to the soil surface. The larger particles however, will return very quickly to the soil surface from suspension and would make their way to the end of the plane by saltation. Therefore, at the beginning of an erosion event, most of the contribution to the total suspended sediment concentration c will be coming from the finer particles, but, as time increases and most of the finer material has been removed, the larger particles will begin to increase their contribution to c. Prediction of size distributions of eroded particles then enables estimates to be made of quantities of sorbed pollutants, (e.g. Palis et al., 1990). An understanding of the dynamics of the suspended sediment concentrations ci of the individual particles, as well as the overall total concentration c, therefore has significant implications on the understanding of the supply of non-point source pollutants to waterways. In this paper we will only be considering erosion from situations where raindrop impact is the only erosive agent. Under these conditions the role of overland flow is simply to transport the sediment eroded under the action of raindrops. The fluid flow velocity is low enough so that no sediment is entrained into the fluid by the shear stresses acting between the soil surface and the fluid. This allows us to compare the predictions of the model with the experimental results of Proffitt et al. (1991) and considerably simplifies the mathematical solutions. As the development of the model is discussed in detail in Hairsine and Rose (1991) and Hairsine et al. (1995) we will only give a brief description of it here.
2. Theory Once sediment has become suspended in the overland flow, this sediment will begin to settle or deposit back towards the soil surface. The rate of deposition di is directly dependent on the size of the sediment and can vary by orders of magnitude from fine sediment to large aggregates. The deposition rate di of a given sediment size class i is
G.C. Sander et aL / Journal of Hydrology 178 (1996) 351-367
353
then given by d i = vic i
i = 1,2, ..., I
(1)
where vi is the settling velocity, ci is the sediment concentration in size class i and I is the total number of classes. The settling velocity classes vi are chosen so that there is an equal mass of soil in each class. Since deposition returns sediment to the soil bed, then at any given time some fraction H of the original soil surface is shielded by deposited sediment. Since this deposited layer will be weaker than the original soil, we must distinguish between this layer and the original soil when determining the amount of sediment being detached (or lifted) from the soil surface through raindrop impact. Assuming that the rate of rainfall detachment per unit area of the original soil to be non-selective with respect to size class i, then Rose (1960, 1961) shows that e i = a P / I where a is the rainfall detachability of bare soil and P is the rainfall rate. Since the deposited layer shields a portion H of the surface from the action of raindrops, then to determine the rate of detachment of the original soil surface, e; must be weighted by a factor of (1 - H). Thus during the erosion event, ei is given by ei = - aP i - [1 - H ( x , t)]
i = 1,2, ..., I
(2)
It is assumed that at the commencement of rainfall (t = 0) no deposited layer has formed thus H ( x , t = 0) = 0 and ei = a P / I . The distribution of particles that form the deposited layer is size class dependent owing to Eq. (1), thus let Mdi be the mass per unit area in this layer of sediment in size class i. The total mass Mdt of sediment in this layer is therefore given by I
Mdt(X, t) = E
Mdi(X, t)
(3)
i=1
Then assuming that the redetachment rate of size class i, edi , by rainfall of sediment in the deposited layer will be proportional to the fraction of particles of that size in the layer, and the fractional shielding of that layer, we have ~gd/
.
edi = a d r Md---~~-I i = 1,2,..,, I
(4)
The deposited layer is expected to have different strength characteristics compared to the original soil surface so the detachability of the deposited layer ad in Eq. (4) wil! be greater than a in Eq. (2). The change in space and time of the fraction of shielding H, must depend on the total mass of sediment per unit area Mdt in the deposited layer. Thus if we define M~t as the mass per unit area of sediment necessary for complete effective shielding (i.e. H = 1), H is given as (Hairsine et al., 1995) H ( x , t) = Mdt(x, t)
(5)
For a plane land element of uniform slope under a spatially uniform excess rainfall
G.C. Sander et al. / Journal of Hydrology 178 (1996) 351-367
354
R(t), conservation o f mass of sediment of size class i requires that O(Dci) O(qci) ~ ----ei+edi--di i=1,2,...,I (6) Ot ax where D(x, t) and q(x, t) are the depth o f flow and the volumetric flux per unit width o f plane, respectively. In the case of shallow kinematic overland flow, D and q are related by
01)0t ~--~xOq= R( t)
(7)
q = KD m
(8)
and
where K and m are constants depending on the flow regime and R(t) is the rainfall excess. In particular K = sl/2/n where s is the surface slope, n is a measure of the surface roughness and m is approximately 5/3 for turbulent flow and 3 for laminar flow. With D and q known Eqs. (6) and (7) can also be combined to give
Oci q Oci 1 Ot + D Ox - 1) (ei + edi - di - Rci)
i = 1,2, ..., I
(9)
Finally, we need a conservation equation which governs the variation o f Mdi in the deposited layer. The only processes which effect the deposited layer are di and edi, thus the change in mass Mdi in the deposited layer from the commencement of rainfall detachment is simply the difference between di and edi or OMdi
--=di-edi Ot
i---- 1,2,...,I
(10)
Therefore Eqs. (9) and (10) form a system of N ( N -- 21) partial differential equations for determining the N unknowns ci(x, t) and Mdi(X, t).
3. Approximate analytical solution In this section we seek an approximate analytical solution to Eqs. (9) and (10) with a view to reproducing the experimental data ofProffitt et al. (1991). Their experiments were performed for constant rainfall rates P, low slopes ( ~< 1%), three shallow mean flow depths o f 2, 5 and 10 mm, a n d no losses o f rainfall by infiltration, hence R(t) = P = constant. For Proffitt et al. (1991) to be able to maintain these mean flow depths for t/> 0 throughout their experiments " . . . a barrier was temporarily installed at the end of the flume and clear water introduced gently so that there was no turbulence. Once the correct mean depth o f ponding had been achieved, rainfall was started and the barrier removed simultaneously. Transient unsteady [water] flow occurring prior to achieving a steady-state appeared to be negligible after approximately one minute." If we are prepared to ignore the initial transient hydrological effects on ci, we can therefore simplify Eq. (9) by replacing D and q with their averaged quantities which are independent of x and t, without introducing any
G.C. Sander et al. / Journal o f Hydrology 178 (1996) 351-367
355
significant errors. So, for the experimental conditions of Proffitt et al. (1991) the following dimensionless variables can be defined
Pt
T = -=D
(lla)
x
z = -
(1 lb)
L
ci Ci = - -
(llc)
a
mdi -~
Mdi
(1 ld)
aD
I
mdr
-- gdt
aye)
--
Emdi
(lle)
i=1
and dimensionless parameters vi
ui = ~
(1 lf)
e = --
(llg)
PL
o~ --
adZ)
- -
(1
lh)
and = -ada
(1 li)
where L is length of the flow domain. Combining Eqs. (1), (2), (4), (5), (9), (10) and (11) gives
OCi~_eOCi 1 ( a ) Or Oz = I 1 - ~ m d ~ +amdi--(l+ui)Ci
i=1,2,...,I
(12)
and Omdi O'F ~- v i C i -- a m d i
i=
1,2,...,I
(13)
Even though Eqs. (12) and (13) form a system of linear partial differential equations there is no exact analytical solution available. Consequently we seek a simplified approximate analytical solution which will still capture the underlying physics of the erosion processes. Experimental data (field or laboratory) are typically collected at the end of the field slope e.g. Hudson (1981), or at the end of the t i m e , e.g. Proffitt et al. (1991). So in order to compare with this data we only need solutions of Eqs. (12) and (13) at z = L. For low slopes and low water flow velocities it is reasonable to
G.C. Sander et aL / Journal of Hydrology 178 (1996) 351-367
356
assume that t h e sediment concentration d o e s not vary greatly with z near z = L, certainly for early times in an erosion event the greatest changes in Ci near z = L are time dependent. Now if we assume that this is true for all time i.e. OCi/t~'F > > E(OCi/OZ)z=l., then the solution for Ci(L , 7-) can be approximated through the much simpler set of linear autonomous ordinary differential equations dr =
1--~md.r
+amdi--(l+vi)Ci
i=1,2,...,I
(14)
and dmdi
d
= l/iCi - amdi
i = 1,2,...,I
(15)
subject to the initial conditions 7- = 0,
C i -~ 0,
mdi = 0
i = 1,2, ...,I
(16)
Eqs. (14) and (15) are now easily solved by standard techniques. First write Eqs. (14) and (15) in vector notation as dw dr
(17)
Aw + b
where W = ( C l , C 2 , ..., e l , m d l , md2, ..., r o d / ) T
(18)
b=
(19)
(}ll
,j,...?,o,o,...,ao
)T
and A is a block (N x N) square matrix A=
[ DIS ] D2 D3
(20)
The matrices D1, DE and D3 are diagonal (I x I) matrices and Sis a symmetric (I x I) matrix. D1 has diagonal elements - ( 1 + vi), i = 1,2, ...,I, 1)2 has diagonal elements vi, i = 1,2, ..., I, and D3 has diagonal elements a, i = 1,2, ..., I. The matrix S is full, with diagonal elements [a - a / ( t ~ I ) ] , i = 1,2,..., I and all the off-diagonal elements are given by - a / ( t ~ I ) . The solution of Eq. (17) is then given by W = ~l~.gl exp(--Alr) + ~2~Lg2exp(--A2r) + ... + ~N~IN exp(--ANr) + Ws
(21)
where Aj and t~j,j = 1,2, ..., N, are the eigenvalues and eigenvectors of A respectively. The r/j,j = 1,2, ..., N are the constants of integration found from satisfying the initial condition r = 0, w = 0, Ci = 0 and mai = 0, or from solving En = -ws
(22)
where E i s the matrix of eigenvectors t~j with the eigenveetor t~j forming colunmj, and Ws is the steady-state solution of Eq. (17).
G.C. Sander et aL /
Journal of Hydrology 178 (1996) 351-367
357
4. Steady-state or equilibrium solution The steady-state or equilibrium solution of Eq. (17) is given when both dCi/d'r and
dmdi/d'r are zero or li ( 1 - ~ md~ ) + a m d i -- (1 + v i ) C i = 0
(23)
viC i - otmdi = 0
(24)
the solution of which is given by Hairsine and Rose (1991) as Ci =
1+
and
md~ = -
i = 1,2, ...,I
Vi[j.=~l( ~)]--1 1+
(25)
i = 1,2, . . . , I
(26)
We note at equilibrium Ci is independent of i, that is, vary with i.
C 1 =
C2,
... , =
CI
but mdi does
5. Comparison with the experimental data of Proflitt et al. (1991) The experiments performed by Proffitt et al. (1991) were on two bare soils of
•
40.00
E 30.00
20.00
d
cO (..) 10.00 O o
0.0o o
.o
6 ........
' .........
lo.oo
' .................
20.00
30.00
' .........
40.00
'
50.00
Time (mins) Fig. 1. Sediment concentration for the Solonchak as a function of time when P = 100 m m h -1 and D = 2 nun. Parameter values are a = 1000, ~ = 20 and a = 1233. Experimental data points given by the stars.
358
G.C. Sander et al. /Journal of Hydrology 178 (1996) 351-367
E
•.x- 3 0 . 0 0
E E "~'~20.00
6
c- 1 0 . 0 0 0 LD ~D -4.--'
0 F-
0.00
i 1 , , 1 , , 1
o.oo
i 1 , 1 1 , , , , 1
lo.oo
ii,
i i i i i i i i i i i i i i i i
20.00
30.00
i
i , , , i i i i I
40.00
50.00
Time (mins) Fig. 2. Sediment concentration for the Solonchak as a function of time when P = 100 mm h -1 and D = 5 nun. Parameter values are a = 1200, ~ = 20 and a = 718. Experimental data points given by the stars. different type (a vertisol referred to as Black E a r t h a n d a n aridisol referred to as S o l o n c h a k ) at slopes between 0.1 a n d 1.0%. Both sediment c o n c e n t r a t i o n a n d settling velocity characteristics were m e a s u r e d t h r o u g h time u n d e r two different c o n s t a n t rainfall rates o f 56 a n d 100 m m h -1 . T h e experiments were also p e r f o r m e d for three different average c o n s t a n t depths o f water o f 2, 5 a n d 10 ram. Figs. 1 - 3 give the e x p e r i m e n t a l c(t) (stars) for the S o l o n c h a k for depths o f 2 m m ,
E
.
10.00
E E t2n 5.o0
d
tO
0 b-
Fig. 3.
0.00
''''''I;'i''',li,,,l,1,,,l,,,l,,,,,,,,,l,,a,,,,,i
o.oo
lo.oo
20.00
30.00
Time (mins)
Sediment concentration for the Solonchak as
a
40.00
function of time when P
l
50.00
=
100 nun h-1 and D
=
I0 m m . P a r a m e t e r values are c~ = 1600, ~ = 20 and a = 412. E x p e r i m e n t a l d a t a points given by the stars.
359
G.C. Sander et al. / Journal of Hydrology 178 (1996) 351-367
5 m m and 10 mm, respectively, when P = 100 mm h -1. The solid lines in Figs. 1-3 are the predicted solutions obtained from Eqs. (21) and (22) for ten size classes, i.e. I = 10. The corresponding fall velocities vi, i = 1,2, ..., 10 used, were obtained from the soil curve in Fig. 2 of Proffitt et al. (1991 ) labelled original soil (wetted plus rainfall at 100 m m h -l for 40 min). The ten velocities were then simply read off at equal intervals along their y-axis. The finest sediment class size is given by i = 1 (smallest velocity) with the largest class size corresponding to i = 10 (largest velocity). The parameter values o f c~, ~; and a are obtained by matching the predicted concentration curve to the experimental data points. The level o f agreement between the experimental data and our predicted curves is excellent, well within the experimental error of +15%. The most interesting feature of the solutions is the very rapid increase in the concentration from 0 at t = 0, to its peak concentration which occurs around 1 rain, followed by a continuous decline towards its steady-state or equilibrium concentration. The peak concentration predicted by the theory is not absolute in that by making slight variations to the values of a, n and a, we can either increase or decrease the peak predicted concentration while still maintaining the excellent agreement with the experimental data. This is discussed in more detail in a later section. What is needed is to obtain more measured sediment concentrations at times earlier than 1.5 min in order to fix the peak concentration. It was also found that when P = 56 m m h -] exactly the same level of agreement was obtained between theory and experimental data as that displayed for P = 100 mm h -]. Note from now on the bar has been dropped from the depth D though references to D will still imply averaged depths. Fig. 4 shows the development of the deposited layer during the erosion event. While the sediment concentration takes quite a while to obtain equilibrium, the level of shielding obtained by the deposited layer attains equilibrium very rapidly. The equilibrium shielding value o f approximately 98% seems to be achieved on the 1.00
0.75
~- 0.50
0.25
0.00
0.00
IIJlllllllSlll
10.00
Jllll
Il'
j'Ll~llll*~hjl'ljl'ljlllll
20.00 30.00 TIME (mlns)
40.00
50.00
Fig. 4. Developmentof the fractionalcoverageH as a functionof time for P = 100mm h-I and D
= 5
ram.
G.C. Sander et al. / Journal of Hydrology 178 (1996) 351-367
360
0.030
i=
10
i=
8
E ,w. E .2£
"-" 0.020 (/3
5 0.010 ~J3 .<
0.000 0.00
10.00
20.00 30.00 TIME (mlns)
4O.uu
~u.O0
Fig. 5. Variation of the masses of the individual particle size classes making up the deposited layer as a function oftime, P = 100 ram h -I and D = 5 ram. same time-scale as the predicted p e a k concentrations shown in Fig. 2. While the level o f shielding appears to be at equilibrium, Fig. 5 shows that mass o f sediment in the various class sizes which f o r m the deposited layer is still varying significantly with time. Even at 50 min it is obvious that the majority o f size classes have not achieved their equilibrium values, Thus, as the erosion event proceeds there is a continuous change in the distribution o f the individual particle sizes which m a k e up the deposited layer. As time goes on, Fig. 5 shows that the deposited layer becomes m o r e and m o r e d o m i n a t e d by the coarser material since the finer material is transported further following each detachment. This is in agreement with prediction o f the steady-state solution for mdi given by Eq. (26). True equilibrium is achieved on a time-scale o f 5.00
~, 2.00 U
.N_ "6 1.00 u 0 £.)
"~.1o 0.00
,,l,,,,,,l '''''''''l'''l't=''l
0.00
10.00
'''''''''I'''''''''I
20.00 30.00 Time (mins)
40.00
50.00
Fig. 6. Variation of the suspended sediment concentration of the individual size classes with time for P = 100 mm h -I and D = 5 mm.
361
G.C. Sander et aL / Journal of Hydrology 178 (1996) 351-367
approximately 1500-2000 min being of the order of one day. In Fig. 6 we see the time-varying behaviour of the suspended sediment concentration for the various class sizes corresponding to Fig. 2. This demonstrates quite clearly that early on, it is the finer sediment which forms the major contribution to the predicted concentration peak in Fig. 2. This confirms the well-known experimental results (Meyer et al., 1975; Alberts et al., 1980, 1983) that sediment leaving an eroded slope is finer than the original soil. It has been noted by Rose and Dalai (1988) that this is a key mechanism in the enrichment of soil-sorbed nutrients leaving an eroded slope. The results in Fig. 6 can also be used to calculate settling velocity distributions of the eroded sediment and then as input to the method of Palis et al. (1990) for predicting losses of sorbed nutrients. Since Eq. (13) shows that at equilibrium all the ci's are equal, then we see from Fig. 6 that steady state has not been reached at 50 min but it is obvious that all the ci's are converging to a single value. In total, there were 12 different experiments performed (three different depths by two rainfall rates by two soil types). In three experiments for the Black Earth soil there was inflow of fresh water from the top of the flume. The parameter values for these runs were P = 56 mm h -l, D = 5 and 10 ram, and P = 100 mm h -l with D = 10 mm. No comparisons of this experimental data with the approximate analytic solution was performed as it was felt that the fresh water inflow at x = 0 would violate the assumptions under which the solution was developed. However, in all remaining nine experiments, agreement between measured and predicted c(t) was obtained to the same level as that demonstrated in Figs. 1-3. In fact, the type of behaviour displayed throughout Figs. 4-6 is very similar in all experiments. The detachabilities for the original soil a and the deposited layer a d for all nine experiments are summarised in Table 1. On physical grounds, as the depth of water covering the soil surface increases, then this water provides greater protection to the soil surface from rainfall. Hence, as D increases the rainfall detachability of both the original soil and the deposited layer should decrease. This behaviour is clearly visible in Table 1. It is also clear that for a given soil type, rainfall rate and depth of flow, ad >> a, confirming the greater detachability of the deposited layer over the original Ta~el Listofparametervalues ~ v i n g b e s t a ~ m e n t w i t h t h e e x p e f i m e n t ~ d a t a Solonchak P (mmh -i)
Black Earth D (mm)
a (kgm -3)
ad (kgrn -3)
M~t (kgm -2)
56
2 5 10
1113 358 319
22260 7160 6380
0.074 0.18 0.071
100
2 5 10
1233 718 412
24660 14360 8240
0.05 0.06 0.051
P (mmh -1)
D (mm)
a (kgrn -3)
ad M~t (kgm -3) (kgm -E)
56
2
3910
7429
0.15
100
2 5
3738 1950
7102 3705
0.18 0.23
G.C. Sander et al. / Journal of Hydrology 178 (1996) 351-367
362
soil. Being aware of this behaviour helps restrict the range of possible parameter values a and t~ which yield good agreement with the experimental c(t) data. We note that the values of the detachability a found in Table 1 are much higher than those reported experimentally in Proffitt et al. (1991). This difference arises because the experimental values o f a were based on an equilibrium value of H = 0.9 obtained by visual observation. In contrast Fig. 4 shows that H should be much closer to one, or to be precise 0.98. Since the detachability increases very rapidly as a function of H for H near one (see Fig. 4 o f Proffitt et al., 199 I), then small errors in estimating H near one are magnified, and cause very large errors in determining a. Thus it is not really appropriate to compare the experimental and analytical values of a except to note that the sign of the difference is correct. On the other hand, ad depends only mildly on H for H near one and we find good agreement between the experimental and analytical values. Fig. 4 of Profitt et al. (1991) gives experimental values ofad for D = 2, 5 and 10 mm of the order of 23 000 kg m -3, 12 000 kg m -a and 7500 kg m -3, respectively, which all agree with the values reported in Table 1. Hairsine and Rose (1991) developed theory for the general case where D is a function o f distance x and the detachabilities are functions o f depth. They considered expressions for a and ad of the form a = ao
(27a)
D < Do
a = a o ( D / D o ) -b
D >.Do
(27b)
and a d = ado
(28a)
D < Do
as = ado(Z)/no) -b
D >1Do
(28b)
In Eqs. (27) and (28) Do is usually taken to be o f the order o f two to three raindrop diameters being about 2 mm for the experiments of Proffitt et al. (1991). Using the data available in Table 1 we can estimate the parameters in Eq. (27) by performing a least squares curve fit. The parameter ado in Eq. (28) then follows automatically from Eqs. (1 li), (27) and (28) for D >~Do a / a d = ao/ado = n
(29)
For the purposes of the least squares fit Eq. (27) can be simplified to a = w D -b where w = aoDbo and results in a = 1730D -°'8
P = 56 mmh -1
(30a)
a = 2017D -°'7
P = 100 mmh -1
(30b)
and
Taking Do as approximately 2 m m gives ao = 994 kg m -3, ado = 19880 kg m -3 for P = 56 mm h -1 and ao = 1242 kg m -3, ado = 24840 kg m -3 for P = 100 mm h -1 . The agreement between Eq. (5.1) and the data in Table 1, is shown in Fig. 7. The results for P = 100 are excellent, however those for P = 56 are not as good. This is due mainly to
363
G.C. San~retal./JournalofHy~ology178(1996) 351-367 1400.0
1200.0 A
1000.0
800.0
25
Q
600.0 Q)
400.0 A
=
200.0
0.O. 00 ........ 2.001......
q
q
i
T
~
i
t
i
i
~
i
~
L
~
4.00 6.00 Depth (mm)
i
i
i
i
~
i
i
i
i
i
8.00
,
i
i
i
i
i
i
J
i
10.00
Fig. 7. Detachability as a function of depth for different rainfall rates. Solid lines are from Eq. (27) based on a least squares regression fit to the data points shown in Table 1.
the data point at D = 5 for P = 56 being too low. We also note that the values of b derived from the least squares procedure are in general agreement with those found by Proffitt et al. (1991). Table 1 also gives the variation of M~t with D for various rainfall rates. These are the least favourable results, since M~t represents the mass per unit area (kg m -2) of the deposited layer which provides complete shielding, these values appear to be lower than physically expected. There also does not appear to be any uniform relation between M~t and depth. This could mean that the definition of H(t) given by Eq. (7) is not quite correct and needs further development. However all other comparisons between theory and experiment (i.e. Figs. 1-6) are very encouraging and lead us to believe that considering the individual class sizes as opposed to just a total concentration per se is more useful in explaining soil erosion behaviour. Certainly it is not possible to predict the sudden peak concentration of Figs. 1-3 by using less than three particle sizes (I = 3). Even then these three size classes must cover several orders of magnitude in their settling velocities vi or the peak will still not appear.
5.1. Curve fitting to experimental c(t) data To curve fit the experimental data and determine the parameter values for a, a d and M~tt we used the solution in dimensionless form as given by Eqs. (21) and (22) since
364
G.C. Sander et al. / Journal of Hydrology 178 (1996) 351-367
only two unknown parameters a and n appear rather than three. The settling velocity distribution vi is obtained directly from settling tube experimentals (as in Fig. 2 o f Proffitt et al., 1991) once the number of size classes/has been set. For a given rainfall rate experiment the vi's are then determined from Eq. (1 If) as vi = vi/P, i = 1,2, ..., I, with Ci(r) and mdi(7-) then calculated from Eqs. (21) and (22) for a preset a and ~. The total dimensionless concentration C(r) is then found by summing the Ci(r) for i = 1,2,...,I. To find the detachability a we match the experimental sediment concentration at t = 40 min, ce(40), with the theoretical sediment concentration C at 7- = 4 0 P / D and then use Eq. (1 lc) to give a = ce(40)/C(4OP/D). The time of 40 min was chosen since this was the largest time for which experimental c data were available and as such this data point would be the most reliable owing to a minimum influence of transient flow effects. A consequence of this method for finding the detachability is that the theoretical solution will always go through the final experimental data point. Having a, the entire theoretical c(t) solution can be plotted and compared with all the experimental data. If the level o f agreement is unsatisfactory then the whole process is repeated with new values of a and n until the 'best' fit is found. With the final values o f a, ~; and a, Eq. (1 li) gives a d and Eq. (1 lh) gives M~t. The effect o f the parameters a and ~ on the curve fitting o f the solution is best seen through their influence on the magnitude and timing of the peak sediment concentration. Define Cmaxas the peak or maximum concentration reached and tmax as the time at which this peak occurs. If we rerun the model several times during which we keep fixed while a is decreased for each additional run, then we find that both tmax and Cmax increase, that is the peak concentration gets larger and is delayed, while the associated values o f a and ad decrease with M~t increasing. This occurs because if the detachabilities a and ad decrease, then to maintain the supply of soil particles in order to match the measured experimental sediment concentrations, the fraction H of shielding must decrease, In practice, this occurs through increasing nit in Eq. (5) so that the product a(1 - H), being the source of detached material, remains constant. In terms of the distribution o f the size classes in the eroded sediment, a fixed ~; with a decreasing a will decrease the contribution from the larger particles. Fixing the value of a and letting ~; decrease for each run, results in an increased peak concentration which occurs at an earlier time. However M~t and ad will remained fixed while the detachability a increases. Furthermore, the total mass o f suspended sediment increases, but there is a tendency towards a reduction in the percentage o f smaller particles contributing to this sediment. Taking the two parameters in combination, then Cmax can be reduced but kept at the same time tmax by increasing both a and ~;, or conversely, the peak concentration can be held constant and be made to occur earlier by increasing a and decreasing ~;, or to occur later by decreasing a and increasing n. We also note that when n is fixed and we vary a, the detachability only varies slightly. This occurs because we chose to find the detachability by matching the theoretical C to ce(40 ) which is near to the steady-state concentration. If we were able to use the true steady-state concentration to find a by a = Ce(t ~ o o ) / C ( r ~ c~ ) then a would remain independent of a as Eq. (25) shows that
G.C. Sander et aL / Journal of Hydrology 178 (1996) 351-367
365
C(r ~ oo ) is independent of a. Thus the small variations we see in calculating the detachability for fixed ~ simply is due to the difference between ce(40) and ce(t ~ oo ). Conversely when we hold a fixed and vary ~ during the curve fitting, then large variations in the calculated detachabilities will be found which are proportional to the changes in ~. 5.2. Sensitivity to number o f size classes It was discussed earlier that it was possible to obtain a very good fit with the experimental data with a range o f values for a and ~;, and therefore a, ad and M~t. In this section we are interested in how these ranges are affected by the number of particle size classes. Specifically, we are looking for an optimal number of particle class sizes lp above which the suitable ranges o f values for t~ and ~ which reproduce the experimental data are essentially independent of L This is particularly important when seeking solutions for the spatial dependence of the concentration through time from Eqs. (9) and (10). Any reduction in the number of particle size classes leads to a reduction of twice that number of partial differential equations, and would therefore have a significant impact on the computation time required to solve numerically for c(x, t) from Eqs. (9) and (10). F o r the purposes o f examining the sensitivity of the model t o / , the experiment on the Black Earth soil for a 5 m m depth and a 100 mm h -l rainfall rate is used. The results for I = 5, 10, 15 and 20 are summarised in Table 2. The parameter ranges given in Table 2 are based on good agreement between the theoretical c(t) and the experimental concentrations, and also on whether the settling velocity distribution curves are acceptable. For all particle size classes it was found that ~cmust always lie in the range 1.9 ~
Table 2 Effect o f n u m b e r o f particle size I on parameter values I
a
~-
1.9
5 10 15 20
160~a~220 6 0 ~ a ~ 120 40 ~ a ~ 80 40 ~
2547 ~
~ 2735 ~ 2097 ~ 1892 ~< 1809
~ = 10.0
M,~t
4 8 5 ~ a ~ 500 352 ~ a ~ 4 0 0 318 ~ a ~ 361 795 ~
0.12 ~
366
G.C. Sander et al. / Journal of Hydrology 178 (1996) 351-367
6. Conclusion I n conclusion we have presented a time-dependent solution for the soil erosion model o f Hairsine and Rose (1991) and shown that we can reproduce the experimental data o f Proffitt et al (1991). This simple analytical solution is based on the assumption that the spatial dependence o f the sediment concentration at x = L is small c o m p a r e d with its time dependence. D a t a f r o m nine experiments were used to verify the reliability o f the theoretical model with excellent agreement f o u n d for all nine. We were able to reproduce the initial rapid increase to a peak in sediment concentration which occurred approximately 5 rain after the c o m m e n c e m e n t o f rainfall, as well as the tailing exponential decline f r o m this peak t o steady-state conditions. W e show that while the total sediment concentration and the fraction o f shielding resulting f r o m the deposited layer appear to reach equilibrium fairly quickly, the distribution o f particle size classes which comprise this sediment is far f r o m equilibrium. True steady state for the Solonchak experiments was achieved between 1500 and 2000 min being on the order o f one day. Sensitivity to the n u m b e r o f particle size classes was investigated and it was decided that ten classes were sufficient to reproduce the structure and physics underlying the experimental data. As it is only possible to solve for c(x, t) numerically, and that any reduction in the n u m b e r o f size classes leads to a reduction o f twice as m a n y partial differential equations, then k n o w i n g Ip = 10 has significant practical implications for determining the spatial as well as the temporal dependence o f the sediment concentration. This problem is covered in future publications.
References Alberts, E.E., Moldenhauer, W.C. and Foster, G.R., 1980. Soil aggregates and primary particles transported in rill and interill flow. Soil Sci. Soc. Am. J., 44: 590-595. Alberts, E.E., Wendt, R.C. and Priest, R.F., 1983. Physical and chemical properties of eroded soil aggregates. Trans. ASAE, 26: 465-471. Govindaraju, R.S. and Kavvas, M.L., 1991.Modelling the erosion process over steep slopes: approximate analytical solutions. J. Hydrol., 127: 279-305. Hairsine, P.B. and Rose, C.W., 1991. Rainfall detachment and deposition: sediment transport in the absence of flow-driven processes. Soil Sci. Sot:. Am. J., 55: 320-324. Hairsine, P.B., Sander, G.C., Rose, C.W., Parlange, J.-Y. and Hogarth, W.L., 1995. Unsteady soil erosion due to rainfall impact: I. Theoretical considerations and practical applications. Water Resour. Res., submitted. Hudson, N., 1981. Soil Conservation, 2nd edn. Cornell University, Ithaca, NY. Laguna, A. and Giraldez, J.V., 1993. The description of soil erosion through a kinematic wave model. J. Hydrol., 145: 65-82. Meyer, L.D., Foster, G.R. and Romkens, M.J.M., 1975. Source of soil eroded by water from upland slopes. In: Present and Prospective Technology for Predicting Sediment Yields and Sources, Proc. Sediment Yield Workshop, Oxford, MS, 28-30 November 1972. USDA-ARS-40, pp. 177-189. Palis, R.G., Okwach, G., Rose, C.W. and Saffigna, P.G., 1990. Soil erosion processes and nutrient loss. I. The interpretation of enrichment ratio and nitrogen loss in runoff sediment. Aust. J. Soil Res., 28: 623-639.
G.C. Sander et al. / Journal of Hydrology 178 (1996) 351-367
367
Proffitt, A.P.B., Rose, C.W., and Hairsine, P.B., 1991. Rainfall detachment and deposition: experiments with low slopes and significant water depths. Soil Sci. Soc. Am. J., 55: 325-332. Rose, C.W., 1960. Soil detachment caused by rainfall. Soil Sci., 89: 28-35. Rose, C.W., 1961. Rainfall and soil structure. Soil Sci., 91: 49-54. Rose, C.W. and Dalal, R.C., 1988. Erosion and runoff of nitrogen. In: J.R. Wilson (Editors), Advances in Nitrogen Cycling in Agricultural Ecosystems. CAB International, Wallingford, UK, pp. 212-235.
Journal of
Hydrology Journal of Hydrology 178 (1996) 351-367
Unsteady soil erosion model, analytical solutions and comparison with experimental results G.C. Sander a'*, P.B. Hairsine b, C.W. Rose c, D. Cassidy a, J.-Y. Parlange d, W.L. Hogarth c, I.G. L i s l e c aFaculty of Science and Technology, Griffith University, Nathan, Qld. 4111, Australia bCSIRO Division of Soils, P.O. Box 639, Canberra ACT, Australia CFaculty of Environmental Sciences, Griffith University, Nathan, Qld. 4111, Australia dDepartment of Agricultural and Biological Engineering, Cornell University, Ithaca, N Y 14853, USA Received 9 March 1995; accepted 31 March 1995
Abstract Hairsine and Rose developed a soil erosion model which described the erosion transport of the multiparticle sizes in sediment for rain-impacted flows in the absence of entrainment in overland flow. In this paper we extend their steady-state solutions to account for the time variation of suspended sediment concentration during an erosion event. A very simple approximate analytical solution is found which agrees extremely well with experimental data obtained from nine experiments. We are able to reproduce the rapid initial increase to a peak in the total sediment concentration, which occurs about 3-5 min after the commencement of rainfall, as well as the subsequent declining exponential tail towards steady-state conditions. We are also able to show that the fraction of shielding of the original soil bed resulting from depositing sediment reaches its equilibrium value on about the same time-scale as the total peak suspended sediment concentration. Interestingly, we find that the masses of the individual particles which form this deposited layer are far from equilibrium, and that there is a great deal of continuous reworking and sorting of this material during the erosion event. Finally, our solution shows that the initial peak in the total sediment concentration is due to the enrichment of this sediment by the finer size classes and that as the event continues their percentage contribution diminishes.
1. Introduction I n a recent p a p e r b y H a i r s i n e a n d R o s e (1991) a soil e r o s i o n m o d e l h a s been d e v e l o p e d to d e s c r i b e the e r o s i o n a n d t r a n s p o r t o f s e d i m e n t in r a i n - i m p a c t e d flows * Corresponding author. 0022-1694/96/$15.00 © 1996 - Elsevier Science B.V. All rights reserved SSDI 0022-1694(95)02810-2
352
G.C. Sander et al. / Journal of Hydrology 178 (1996) 351-367
in the absence of entrainment by shallow overland flow. Soils are naturally composed of a range of particle sizes and aggregates, and in the past, soil erosion models have tended to lump or average over all size classes, resulting in a single conservation equation for the total suspended sediment concentration, c, as a function of distance x and time t (Govindaraju and Kavvas, 1991; Laguna and Giraldez, 1993). In contrast, the approach of Hairsine and Rose (1991) is to consider the contributions of the individual size classes to the total suspended sediment concentration separately. This leads to separate coupled conservation equations for the suspended sediment concentration ci(x, t) for each particle size class i. While this does increase the number of equations and complicates their solution, there are sound physical grounds for justifying this approach. Once particles have become suspended, they will naturally begin to settle back towards the soil surface owing to their immersed weight. The fine particles, having relatively small settling velocities, can be carried significant distances in the direction of the surface water flow before depositing back to the soil surface. The larger particles however, will return very quickly to the soil surface from suspension and would make their way to the end of the plane by saltation. Therefore, at the beginning of an erosion event, most of the contribution to the total suspended sediment concentration c will be coming from the finer particles, but, as time increases and most of the finer material has been removed, the larger particles will begin to increase their contribution to c. Prediction of size distributions of eroded particles then enables estimates to be made of quantities of sorbed pollutants, (e.g. Palis et al., 1990). An understanding of the dynamics of the suspended sediment concentrations ci of the individual particles, as well as the overall total concentration c, therefore has significant implications on the understanding of the supply of non-point source pollutants to waterways. In this paper we will only be considering erosion from situations where raindrop impact is the only erosive agent. Under these conditions the role of overland flow is simply to transport the sediment eroded under the action of raindrops. The fluid flow velocity is low enough so that no sediment is entrained into the fluid by the shear stresses acting between the soil surface and the fluid. This allows us to compare the predictions of the model with the experimental results of Proffitt et al. (1991) and considerably simplifies the mathematical solutions. As the development of the model is discussed in detail in Hairsine and Rose (1991) and Hairsine et al. (1995) we will only give a brief description of it here.
2. Theory Once sediment has become suspended in the overland flow, this sediment will begin to settle or deposit back towards the soil surface. The rate of deposition di is directly dependent on the size of the sediment and can vary by orders of magnitude from fine sediment to large aggregates. The deposition rate di of a given sediment size class i is
G.C. Sander et aL / Journal of Hydrology 178 (1996) 351-367
353
then given by d i = vic i
i = 1,2, ..., I
(1)
where vi is the settling velocity, ci is the sediment concentration in size class i and I is the total number of classes. The settling velocity classes vi are chosen so that there is an equal mass of soil in each class. Since deposition returns sediment to the soil bed, then at any given time some fraction H of the original soil surface is shielded by deposited sediment. Since this deposited layer will be weaker than the original soil, we must distinguish between this layer and the original soil when determining the amount of sediment being detached (or lifted) from the soil surface through raindrop impact. Assuming that the rate of rainfall detachment per unit area of the original soil to be non-selective with respect to size class i, then Rose (1960, 1961) shows that e i = a P / I where a is the rainfall detachability of bare soil and P is the rainfall rate. Since the deposited layer shields a portion H of the surface from the action of raindrops, then to determine the rate of detachment of the original soil surface, e; must be weighted by a factor of (1 - H). Thus during the erosion event, ei is given by ei = - aP i - [1 - H ( x , t)]
i = 1,2, ..., I
(2)
It is assumed that at the commencement of rainfall (t = 0) no deposited layer has formed thus H ( x , t = 0) = 0 and ei = a P / I . The distribution of particles that form the deposited layer is size class dependent owing to Eq. (1), thus let Mdi be the mass per unit area in this layer of sediment in size class i. The total mass Mdt of sediment in this layer is therefore given by I
Mdt(X, t) = E
Mdi(X, t)
(3)
i=1
Then assuming that the redetachment rate of size class i, edi , by rainfall of sediment in the deposited layer will be proportional to the fraction of particles of that size in the layer, and the fractional shielding of that layer, we have ~gd/
.
edi = a d r Md---~~-I i = 1,2,..,, I
(4)
The deposited layer is expected to have different strength characteristics compared to the original soil surface so the detachability of the deposited layer ad in Eq. (4) wil! be greater than a in Eq. (2). The change in space and time of the fraction of shielding H, must depend on the total mass of sediment per unit area Mdt in the deposited layer. Thus if we define M~t as the mass per unit area of sediment necessary for complete effective shielding (i.e. H = 1), H is given as (Hairsine et al., 1995) H ( x , t) = Mdt(x, t)
(5)
For a plane land element of uniform slope under a spatially uniform excess rainfall
G.C. Sander et al. / Journal of Hydrology 178 (1996) 351-367
354
R(t), conservation o f mass of sediment of size class i requires that O(Dci) O(qci) ~ ----ei+edi--di i=1,2,...,I (6) Ot ax where D(x, t) and q(x, t) are the depth o f flow and the volumetric flux per unit width o f plane, respectively. In the case of shallow kinematic overland flow, D and q are related by
01)0t ~--~xOq= R( t)
(7)
q = KD m
(8)
and
where K and m are constants depending on the flow regime and R(t) is the rainfall excess. In particular K = sl/2/n where s is the surface slope, n is a measure of the surface roughness and m is approximately 5/3 for turbulent flow and 3 for laminar flow. With D and q known Eqs. (6) and (7) can also be combined to give
Oci q Oci 1 Ot + D Ox - 1) (ei + edi - di - Rci)
i = 1,2, ..., I
(9)
Finally, we need a conservation equation which governs the variation o f Mdi in the deposited layer. The only processes which effect the deposited layer are di and edi, thus the change in mass Mdi in the deposited layer from the commencement of rainfall detachment is simply the difference between di and edi or OMdi
--=di-edi Ot
i---- 1,2,...,I
(10)
Therefore Eqs. (9) and (10) form a system of N ( N -- 21) partial differential equations for determining the N unknowns ci(x, t) and Mdi(X, t).
3. Approximate analytical solution In this section we seek an approximate analytical solution to Eqs. (9) and (10) with a view to reproducing the experimental data ofProffitt et al. (1991). Their experiments were performed for constant rainfall rates P, low slopes ( ~< 1%), three shallow mean flow depths o f 2, 5 and 10 mm, a n d no losses o f rainfall by infiltration, hence R(t) = P = constant. For Proffitt et al. (1991) to be able to maintain these mean flow depths for t/> 0 throughout their experiments " . . . a barrier was temporarily installed at the end of the flume and clear water introduced gently so that there was no turbulence. Once the correct mean depth o f ponding had been achieved, rainfall was started and the barrier removed simultaneously. Transient unsteady [water] flow occurring prior to achieving a steady-state appeared to be negligible after approximately one minute." If we are prepared to ignore the initial transient hydrological effects on ci, we can therefore simplify Eq. (9) by replacing D and q with their averaged quantities which are independent of x and t, without introducing any
G.C. Sander et al. / Journal o f Hydrology 178 (1996) 351-367
355
significant errors. So, for the experimental conditions of Proffitt et al. (1991) the following dimensionless variables can be defined
Pt
T = -=D
(lla)
x
z = -
(1 lb)
L
ci Ci = - -
(llc)
a
mdi -~
Mdi
(1 ld)
aD
I
mdr
-- gdt
aye)
--
Emdi
(lle)
i=1
and dimensionless parameters vi
ui = ~
(1 lf)
e = --
(llg)
PL
o~ --
adZ)
- -
(1
lh)
and = -ada
(1 li)
where L is length of the flow domain. Combining Eqs. (1), (2), (4), (5), (9), (10) and (11) gives
OCi~_eOCi 1 ( a ) Or Oz = I 1 - ~ m d ~ +amdi--(l+ui)Ci
i=1,2,...,I
(12)
and Omdi O'F ~- v i C i -- a m d i
i=
1,2,...,I
(13)
Even though Eqs. (12) and (13) form a system of linear partial differential equations there is no exact analytical solution available. Consequently we seek a simplified approximate analytical solution which will still capture the underlying physics of the erosion processes. Experimental data (field or laboratory) are typically collected at the end of the field slope e.g. Hudson (1981), or at the end of the t i m e , e.g. Proffitt et al. (1991). So in order to compare with this data we only need solutions of Eqs. (12) and (13) at z = L. For low slopes and low water flow velocities it is reasonable to
G.C. Sander et aL / Journal of Hydrology 178 (1996) 351-367
356
assume that t h e sediment concentration d o e s not vary greatly with z near z = L, certainly for early times in an erosion event the greatest changes in Ci near z = L are time dependent. Now if we assume that this is true for all time i.e. OCi/t~'F > > E(OCi/OZ)z=l., then the solution for Ci(L , 7-) can be approximated through the much simpler set of linear autonomous ordinary differential equations dr =
1--~md.r
+amdi--(l+vi)Ci
i=1,2,...,I
(14)
and dmdi
d
= l/iCi - amdi
i = 1,2,...,I
(15)
subject to the initial conditions 7- = 0,
C i -~ 0,
mdi = 0
i = 1,2, ...,I
(16)
Eqs. (14) and (15) are now easily solved by standard techniques. First write Eqs. (14) and (15) in vector notation as dw dr
(17)
Aw + b
where W = ( C l , C 2 , ..., e l , m d l , md2, ..., r o d / ) T
(18)
b=
(19)
(}ll
,j,...?,o,o,...,ao
)T
and A is a block (N x N) square matrix A=
[ DIS ] D2 D3
(20)
The matrices D1, DE and D3 are diagonal (I x I) matrices and Sis a symmetric (I x I) matrix. D1 has diagonal elements - ( 1 + vi), i = 1,2, ...,I, 1)2 has diagonal elements vi, i = 1,2, ..., I, and D3 has diagonal elements a, i = 1,2, ..., I. The matrix S is full, with diagonal elements [a - a / ( t ~ I ) ] , i = 1,2,..., I and all the off-diagonal elements are given by - a / ( t ~ I ) . The solution of Eq. (17) is then given by W = ~l~.gl exp(--Alr) + ~2~Lg2exp(--A2r) + ... + ~N~IN exp(--ANr) + Ws
(21)
where Aj and t~j,j = 1,2, ..., N, are the eigenvalues and eigenvectors of A respectively. The r/j,j = 1,2, ..., N are the constants of integration found from satisfying the initial condition r = 0, w = 0, Ci = 0 and mai = 0, or from solving En = -ws
(22)
where E i s the matrix of eigenvectors t~j with the eigenveetor t~j forming colunmj, and Ws is the steady-state solution of Eq. (17).
G.C. Sander et aL /
Journal of Hydrology 178 (1996) 351-367
357
4. Steady-state or equilibrium solution The steady-state or equilibrium solution of Eq. (17) is given when both dCi/d'r and
dmdi/d'r are zero or li ( 1 - ~ md~ ) + a m d i -- (1 + v i ) C i = 0
(23)
viC i - otmdi = 0
(24)
the solution of which is given by Hairsine and Rose (1991) as Ci =
1+
and
md~ = -
i = 1,2, ...,I
Vi[j.=~l( ~)]--1 1+
(25)
i = 1,2, . . . , I
(26)
We note at equilibrium Ci is independent of i, that is, vary with i.
C 1 =
C2,
... , =
CI
but mdi does
5. Comparison with the experimental data of Proflitt et al. (1991) The experiments performed by Proffitt et al. (1991) were on two bare soils of
•
40.00
E 30.00
20.00
d
cO (..) 10.00 O o
0.0o o
.o
6 ........
' .........
lo.oo
' .................
20.00
30.00
' .........
40.00
'
50.00
Time (mins) Fig. 1. Sediment concentration for the Solonchak as a function of time when P = 100 m m h -1 and D = 2 nun. Parameter values are a = 1000, ~ = 20 and a = 1233. Experimental data points given by the stars.
358
G.C. Sander et al. /Journal of Hydrology 178 (1996) 351-367
E
•.x- 3 0 . 0 0
E E "~'~20.00
6
c- 1 0 . 0 0 0 LD ~D -4.--'
0 F-
0.00
i 1 , , 1 , , 1
o.oo
i 1 , 1 1 , , , , 1
lo.oo
ii,
i i i i i i i i i i i i i i i i
20.00
30.00
i
i , , , i i i i I
40.00
50.00
Time (mins) Fig. 2. Sediment concentration for the Solonchak as a function of time when P = 100 mm h -1 and D = 5 nun. Parameter values are a = 1200, ~ = 20 and a = 718. Experimental data points given by the stars. different type (a vertisol referred to as Black E a r t h a n d a n aridisol referred to as S o l o n c h a k ) at slopes between 0.1 a n d 1.0%. Both sediment c o n c e n t r a t i o n a n d settling velocity characteristics were m e a s u r e d t h r o u g h time u n d e r two different c o n s t a n t rainfall rates o f 56 a n d 100 m m h -1 . T h e experiments were also p e r f o r m e d for three different average c o n s t a n t depths o f water o f 2, 5 a n d 10 ram. Figs. 1 - 3 give the e x p e r i m e n t a l c(t) (stars) for the S o l o n c h a k for depths o f 2 m m ,
E
.
10.00
E E t2n 5.o0
d
tO
0 b-
Fig. 3.
0.00
''''''I;'i''',li,,,l,1,,,l,,,l,,,,,,,,,l,,a,,,,,i
o.oo
lo.oo
20.00
30.00
Time (mins)
Sediment concentration for the Solonchak as
a
40.00
function of time when P
l
50.00
=
100 nun h-1 and D
=
I0 m m . P a r a m e t e r values are c~ = 1600, ~ = 20 and a = 412. E x p e r i m e n t a l d a t a points given by the stars.
359
G.C. Sander et al. / Journal of Hydrology 178 (1996) 351-367
5 m m and 10 mm, respectively, when P = 100 mm h -1. The solid lines in Figs. 1-3 are the predicted solutions obtained from Eqs. (21) and (22) for ten size classes, i.e. I = 10. The corresponding fall velocities vi, i = 1,2, ..., 10 used, were obtained from the soil curve in Fig. 2 of Proffitt et al. (1991 ) labelled original soil (wetted plus rainfall at 100 m m h -l for 40 min). The ten velocities were then simply read off at equal intervals along their y-axis. The finest sediment class size is given by i = 1 (smallest velocity) with the largest class size corresponding to i = 10 (largest velocity). The parameter values o f c~, ~; and a are obtained by matching the predicted concentration curve to the experimental data points. The level o f agreement between the experimental data and our predicted curves is excellent, well within the experimental error of +15%. The most interesting feature of the solutions is the very rapid increase in the concentration from 0 at t = 0, to its peak concentration which occurs around 1 rain, followed by a continuous decline towards its steady-state or equilibrium concentration. The peak concentration predicted by the theory is not absolute in that by making slight variations to the values of a, n and a, we can either increase or decrease the peak predicted concentration while still maintaining the excellent agreement with the experimental data. This is discussed in more detail in a later section. What is needed is to obtain more measured sediment concentrations at times earlier than 1.5 min in order to fix the peak concentration. It was also found that when P = 56 m m h -] exactly the same level of agreement was obtained between theory and experimental data as that displayed for P = 100 mm h -]. Note from now on the bar has been dropped from the depth D though references to D will still imply averaged depths. Fig. 4 shows the development of the deposited layer during the erosion event. While the sediment concentration takes quite a while to obtain equilibrium, the level of shielding obtained by the deposited layer attains equilibrium very rapidly. The equilibrium shielding value o f approximately 98% seems to be achieved on the 1.00
0.75
~- 0.50
0.25
0.00
0.00
IIJlllllllSlll
10.00
Jllll
Il'
j'Ll~llll*~hjl'ljl'ljlllll
20.00 30.00 TIME (mlns)
40.00
50.00
Fig. 4. Developmentof the fractionalcoverageH as a functionof time for P = 100mm h-I and D
= 5
ram.
G.C. Sander et al. / Journal of Hydrology 178 (1996) 351-367
360
0.030
i=
10
i=
8
E ,w. E .2£
"-" 0.020 (/3
5 0.010 ~J3 .<
0.000 0.00
10.00
20.00 30.00 TIME (mlns)
4O.uu
~u.O0
Fig. 5. Variation of the masses of the individual particle size classes making up the deposited layer as a function oftime, P = 100 ram h -I and D = 5 ram. same time-scale as the predicted p e a k concentrations shown in Fig. 2. While the level o f shielding appears to be at equilibrium, Fig. 5 shows that mass o f sediment in the various class sizes which f o r m the deposited layer is still varying significantly with time. Even at 50 min it is obvious that the majority o f size classes have not achieved their equilibrium values, Thus, as the erosion event proceeds there is a continuous change in the distribution o f the individual particle sizes which m a k e up the deposited layer. As time goes on, Fig. 5 shows that the deposited layer becomes m o r e and m o r e d o m i n a t e d by the coarser material since the finer material is transported further following each detachment. This is in agreement with prediction o f the steady-state solution for mdi given by Eq. (26). True equilibrium is achieved on a time-scale o f 5.00
~, 2.00 U
.N_ "6 1.00 u 0 £.)
"~.1o 0.00
,,l,,,,,,l '''''''''l'''l't=''l
0.00
10.00
'''''''''I'''''''''I
20.00 30.00 Time (mins)
40.00
50.00
Fig. 6. Variation of the suspended sediment concentration of the individual size classes with time for P = 100 mm h -I and D = 5 mm.
361
G.C. Sander et aL / Journal of Hydrology 178 (1996) 351-367
approximately 1500-2000 min being of the order of one day. In Fig. 6 we see the time-varying behaviour of the suspended sediment concentration for the various class sizes corresponding to Fig. 2. This demonstrates quite clearly that early on, it is the finer sediment which forms the major contribution to the predicted concentration peak in Fig. 2. This confirms the well-known experimental results (Meyer et al., 1975; Alberts et al., 1980, 1983) that sediment leaving an eroded slope is finer than the original soil. It has been noted by Rose and Dalai (1988) that this is a key mechanism in the enrichment of soil-sorbed nutrients leaving an eroded slope. The results in Fig. 6 can also be used to calculate settling velocity distributions of the eroded sediment and then as input to the method of Palis et al. (1990) for predicting losses of sorbed nutrients. Since Eq. (13) shows that at equilibrium all the ci's are equal, then we see from Fig. 6 that steady state has not been reached at 50 min but it is obvious that all the ci's are converging to a single value. In total, there were 12 different experiments performed (three different depths by two rainfall rates by two soil types). In three experiments for the Black Earth soil there was inflow of fresh water from the top of the flume. The parameter values for these runs were P = 56 mm h -l, D = 5 and 10 ram, and P = 100 mm h -l with D = 10 mm. No comparisons of this experimental data with the approximate analytic solution was performed as it was felt that the fresh water inflow at x = 0 would violate the assumptions under which the solution was developed. However, in all remaining nine experiments, agreement between measured and predicted c(t) was obtained to the same level as that demonstrated in Figs. 1-3. In fact, the type of behaviour displayed throughout Figs. 4-6 is very similar in all experiments. The detachabilities for the original soil a and the deposited layer a d for all nine experiments are summarised in Table 1. On physical grounds, as the depth of water covering the soil surface increases, then this water provides greater protection to the soil surface from rainfall. Hence, as D increases the rainfall detachability of both the original soil and the deposited layer should decrease. This behaviour is clearly visible in Table 1. It is also clear that for a given soil type, rainfall rate and depth of flow, ad >> a, confirming the greater detachability of the deposited layer over the original Ta~el Listofparametervalues ~ v i n g b e s t a ~ m e n t w i t h t h e e x p e f i m e n t ~ d a t a Solonchak P (mmh -i)
Black Earth D (mm)
a (kgm -3)
ad (kgrn -3)
M~t (kgm -2)
56
2 5 10
1113 358 319
22260 7160 6380
0.074 0.18 0.071
100
2 5 10
1233 718 412
24660 14360 8240
0.05 0.06 0.051
P (mmh -1)
D (mm)
a (kgrn -3)
ad M~t (kgm -3) (kgm -E)
56
2
3910
7429
0.15
100
2 5
3738 1950
7102 3705
0.18 0.23
G.C. Sander et al. / Journal of Hydrology 178 (1996) 351-367
362
soil. Being aware of this behaviour helps restrict the range of possible parameter values a and t~ which yield good agreement with the experimental c(t) data. We note that the values of the detachability a found in Table 1 are much higher than those reported experimentally in Proffitt et al. (1991). This difference arises because the experimental values o f a were based on an equilibrium value of H = 0.9 obtained by visual observation. In contrast Fig. 4 shows that H should be much closer to one, or to be precise 0.98. Since the detachability increases very rapidly as a function of H for H near one (see Fig. 4 o f Proffitt et al., 199 I), then small errors in estimating H near one are magnified, and cause very large errors in determining a. Thus it is not really appropriate to compare the experimental and analytical values of a except to note that the sign of the difference is correct. On the other hand, ad depends only mildly on H for H near one and we find good agreement between the experimental and analytical values. Fig. 4 of Profitt et al. (1991) gives experimental values ofad for D = 2, 5 and 10 mm of the order of 23 000 kg m -3, 12 000 kg m -a and 7500 kg m -3, respectively, which all agree with the values reported in Table 1. Hairsine and Rose (1991) developed theory for the general case where D is a function o f distance x and the detachabilities are functions o f depth. They considered expressions for a and ad of the form a = ao
(27a)
D < Do
a = a o ( D / D o ) -b
D >.Do
(27b)
and a d = ado
(28a)
D < Do
as = ado(Z)/no) -b
D >1Do
(28b)
In Eqs. (27) and (28) Do is usually taken to be o f the order o f two to three raindrop diameters being about 2 mm for the experiments of Proffitt et al. (1991). Using the data available in Table 1 we can estimate the parameters in Eq. (27) by performing a least squares curve fit. The parameter ado in Eq. (28) then follows automatically from Eqs. (1 li), (27) and (28) for D >~Do a / a d = ao/ado = n
(29)
For the purposes of the least squares fit Eq. (27) can be simplified to a = w D -b where w = aoDbo and results in a = 1730D -°'8
P = 56 mmh -1
(30a)
a = 2017D -°'7
P = 100 mmh -1
(30b)
and
Taking Do as approximately 2 m m gives ao = 994 kg m -3, ado = 19880 kg m -3 for P = 56 mm h -1 and ao = 1242 kg m -3, ado = 24840 kg m -3 for P = 100 mm h -1 . The agreement between Eq. (5.1) and the data in Table 1, is shown in Fig. 7. The results for P = 100 are excellent, however those for P = 56 are not as good. This is due mainly to
363
G.C. San~retal./JournalofHy~ology178(1996) 351-367 1400.0
1200.0 A
1000.0
800.0
25
Q
600.0 Q)
400.0 A
=
200.0
0.O. 00 ........ 2.001......
q
q
i
T
~
i
t
i
i
~
i
~
L
~
4.00 6.00 Depth (mm)
i
i
i
i
~
i
i
i
i
i
8.00
,
i
i
i
i
i
i
J
i
10.00
Fig. 7. Detachability as a function of depth for different rainfall rates. Solid lines are from Eq. (27) based on a least squares regression fit to the data points shown in Table 1.
the data point at D = 5 for P = 56 being too low. We also note that the values of b derived from the least squares procedure are in general agreement with those found by Proffitt et al. (1991). Table 1 also gives the variation of M~t with D for various rainfall rates. These are the least favourable results, since M~t represents the mass per unit area (kg m -2) of the deposited layer which provides complete shielding, these values appear to be lower than physically expected. There also does not appear to be any uniform relation between M~t and depth. This could mean that the definition of H(t) given by Eq. (7) is not quite correct and needs further development. However all other comparisons between theory and experiment (i.e. Figs. 1-6) are very encouraging and lead us to believe that considering the individual class sizes as opposed to just a total concentration per se is more useful in explaining soil erosion behaviour. Certainly it is not possible to predict the sudden peak concentration of Figs. 1-3 by using less than three particle sizes (I = 3). Even then these three size classes must cover several orders of magnitude in their settling velocities vi or the peak will still not appear.
5.1. Curve fitting to experimental c(t) data To curve fit the experimental data and determine the parameter values for a, a d and M~tt we used the solution in dimensionless form as given by Eqs. (21) and (22) since
364
G.C. Sander et al. / Journal of Hydrology 178 (1996) 351-367
only two unknown parameters a and n appear rather than three. The settling velocity distribution vi is obtained directly from settling tube experimentals (as in Fig. 2 o f Proffitt et al., 1991) once the number of size classes/has been set. For a given rainfall rate experiment the vi's are then determined from Eq. (1 If) as vi = vi/P, i = 1,2, ..., I, with Ci(r) and mdi(7-) then calculated from Eqs. (21) and (22) for a preset a and ~. The total dimensionless concentration C(r) is then found by summing the Ci(r) for i = 1,2,...,I. To find the detachability a we match the experimental sediment concentration at t = 40 min, ce(40), with the theoretical sediment concentration C at 7- = 4 0 P / D and then use Eq. (1 lc) to give a = ce(40)/C(4OP/D). The time of 40 min was chosen since this was the largest time for which experimental c data were available and as such this data point would be the most reliable owing to a minimum influence of transient flow effects. A consequence of this method for finding the detachability is that the theoretical solution will always go through the final experimental data point. Having a, the entire theoretical c(t) solution can be plotted and compared with all the experimental data. If the level o f agreement is unsatisfactory then the whole process is repeated with new values of a and n until the 'best' fit is found. With the final values o f a, ~; and a, Eq. (1 li) gives a d and Eq. (1 lh) gives M~t. The effect o f the parameters a and ~ on the curve fitting o f the solution is best seen through their influence on the magnitude and timing of the peak sediment concentration. Define Cmaxas the peak or maximum concentration reached and tmax as the time at which this peak occurs. If we rerun the model several times during which we keep fixed while a is decreased for each additional run, then we find that both tmax and Cmax increase, that is the peak concentration gets larger and is delayed, while the associated values o f a and ad decrease with M~t increasing. This occurs because if the detachabilities a and ad decrease, then to maintain the supply of soil particles in order to match the measured experimental sediment concentrations, the fraction H of shielding must decrease, In practice, this occurs through increasing nit in Eq. (5) so that the product a(1 - H), being the source of detached material, remains constant. In terms of the distribution o f the size classes in the eroded sediment, a fixed ~; with a decreasing a will decrease the contribution from the larger particles. Fixing the value of a and letting ~; decrease for each run, results in an increased peak concentration which occurs at an earlier time. However M~t and ad will remained fixed while the detachability a increases. Furthermore, the total mass o f suspended sediment increases, but there is a tendency towards a reduction in the percentage o f smaller particles contributing to this sediment. Taking the two parameters in combination, then Cmax can be reduced but kept at the same time tmax by increasing both a and ~;, or conversely, the peak concentration can be held constant and be made to occur earlier by increasing a and decreasing ~;, or to occur later by decreasing a and increasing n. We also note that when n is fixed and we vary a, the detachability only varies slightly. This occurs because we chose to find the detachability by matching the theoretical C to ce(40 ) which is near to the steady-state concentration. If we were able to use the true steady-state concentration to find a by a = Ce(t ~ o o ) / C ( r ~ c~ ) then a would remain independent of a as Eq. (25) shows that
G.C. Sander et aL / Journal of Hydrology 178 (1996) 351-367
365
C(r ~ oo ) is independent of a. Thus the small variations we see in calculating the detachability for fixed ~ simply is due to the difference between ce(40) and ce(t ~ oo ). Conversely when we hold a fixed and vary ~ during the curve fitting, then large variations in the calculated detachabilities will be found which are proportional to the changes in ~. 5.2. Sensitivity to number o f size classes It was discussed earlier that it was possible to obtain a very good fit with the experimental data with a range o f values for a and ~;, and therefore a, ad and M~t. In this section we are interested in how these ranges are affected by the number of particle size classes. Specifically, we are looking for an optimal number of particle class sizes lp above which the suitable ranges o f values for t~ and ~ which reproduce the experimental data are essentially independent of L This is particularly important when seeking solutions for the spatial dependence of the concentration through time from Eqs. (9) and (10). Any reduction in the number of particle size classes leads to a reduction of twice that number of partial differential equations, and would therefore have a significant impact on the computation time required to solve numerically for c(x, t) from Eqs. (9) and (10). F o r the purposes o f examining the sensitivity of the model t o / , the experiment on the Black Earth soil for a 5 m m depth and a 100 mm h -l rainfall rate is used. The results for I = 5, 10, 15 and 20 are summarised in Table 2. The parameter ranges given in Table 2 are based on good agreement between the theoretical c(t) and the experimental concentrations, and also on whether the settling velocity distribution curves are acceptable. For all particle size classes it was found that ~cmust always lie in the range 1.9 ~
Table 2 Effect o f n u m b e r o f particle size I on parameter values I
a
~-
1.9
5 10 15 20
160~a~220 6 0 ~ a ~ 120 40 ~ a ~ 80 40 ~
2547 ~
~ 2735 ~ 2097 ~ 1892 ~< 1809
~ = 10.0
M,~t
4 8 5 ~ a ~ 500 352 ~ a ~ 4 0 0 318 ~ a ~ 361 795 ~
0.12 ~
366
G.C. Sander et al. / Journal of Hydrology 178 (1996) 351-367
6. Conclusion I n conclusion we have presented a time-dependent solution for the soil erosion model o f Hairsine and Rose (1991) and shown that we can reproduce the experimental data o f Proffitt et al (1991). This simple analytical solution is based on the assumption that the spatial dependence o f the sediment concentration at x = L is small c o m p a r e d with its time dependence. D a t a f r o m nine experiments were used to verify the reliability o f the theoretical model with excellent agreement f o u n d for all nine. We were able to reproduce the initial rapid increase to a peak in sediment concentration which occurred approximately 5 rain after the c o m m e n c e m e n t o f rainfall, as well as the tailing exponential decline f r o m this peak t o steady-state conditions. W e show that while the total sediment concentration and the fraction o f shielding resulting f r o m the deposited layer appear to reach equilibrium fairly quickly, the distribution o f particle size classes which comprise this sediment is far f r o m equilibrium. True steady state for the Solonchak experiments was achieved between 1500 and 2000 min being on the order o f one day. Sensitivity to the n u m b e r o f particle size classes was investigated and it was decided that ten classes were sufficient to reproduce the structure and physics underlying the experimental data. As it is only possible to solve for c(x, t) numerically, and that any reduction in the n u m b e r o f size classes leads to a reduction o f twice as m a n y partial differential equations, then k n o w i n g Ip = 10 has significant practical implications for determining the spatial as well as the temporal dependence o f the sediment concentration. This problem is covered in future publications.
References Alberts, E.E., Moldenhauer, W.C. and Foster, G.R., 1980. Soil aggregates and primary particles transported in rill and interill flow. Soil Sci. Soc. Am. J., 44: 590-595. Alberts, E.E., Wendt, R.C. and Priest, R.F., 1983. Physical and chemical properties of eroded soil aggregates. Trans. ASAE, 26: 465-471. Govindaraju, R.S. and Kavvas, M.L., 1991.Modelling the erosion process over steep slopes: approximate analytical solutions. J. Hydrol., 127: 279-305. Hairsine, P.B. and Rose, C.W., 1991. Rainfall detachment and deposition: sediment transport in the absence of flow-driven processes. Soil Sci. Sot:. Am. J., 55: 320-324. Hairsine, P.B., Sander, G.C., Rose, C.W., Parlange, J.-Y. and Hogarth, W.L., 1995. Unsteady soil erosion due to rainfall impact: I. Theoretical considerations and practical applications. Water Resour. Res., submitted. Hudson, N., 1981. Soil Conservation, 2nd edn. Cornell University, Ithaca, NY. Laguna, A. and Giraldez, J.V., 1993. The description of soil erosion through a kinematic wave model. J. Hydrol., 145: 65-82. Meyer, L.D., Foster, G.R. and Romkens, M.J.M., 1975. Source of soil eroded by water from upland slopes. In: Present and Prospective Technology for Predicting Sediment Yields and Sources, Proc. Sediment Yield Workshop, Oxford, MS, 28-30 November 1972. USDA-ARS-40, pp. 177-189. Palis, R.G., Okwach, G., Rose, C.W. and Saffigna, P.G., 1990. Soil erosion processes and nutrient loss. I. The interpretation of enrichment ratio and nitrogen loss in runoff sediment. Aust. J. Soil Res., 28: 623-639.
G.C. Sander et al. / Journal of Hydrology 178 (1996) 351-367
367
Proffitt, A.P.B., Rose, C.W., and Hairsine, P.B., 1991. Rainfall detachment and deposition: experiments with low slopes and significant water depths. Soil Sci. Soc. Am. J., 55: 325-332. Rose, C.W., 1960. Soil detachment caused by rainfall. Soil Sci., 89: 28-35. Rose, C.W., 1961. Rainfall and soil structure. Soil Sci., 91: 49-54. Rose, C.W. and Dalal, R.C., 1988. Erosion and runoff of nitrogen. In: J.R. Wilson (Editors), Advances in Nitrogen Cycling in Agricultural Ecosystems. CAB International, Wallingford, UK, pp. 212-235.