Modeling the erosion process over steep slopes: approximate analytical solutions

Journal o f Hydrology, 127 (1991) 279-305

279

Elsevier Science Publishers B.V., A m s t e r d a m

[1]

Modeling the erosion process over steep slopes: approximate analytical solutions Rao S. Govindaraju and M. Levent Kavvas Department (~['Civil Engineering, University of California, Davis. CA 95616. USA (Received 8 October 199(/; revised and accepted 21 November 19901

ABSTRACT Govindaraju, R.S. and Kavvas, M.L., 1991. Modeling the erosion process over steep slopes: approximate analytical solutions. J. Hydrol., 127: 279-305. Analytical expressions are developed for the rainfall-runoff-erosion process on steep hillslopes subjected to time-varying rainfall events. The erosion equation is essentially represented as a lirst-order reaction with the reaction rate being represented by the soil erodibility. The analytical transient solutions are based upon the assumption that the flow and sediment discharge have a constant relationship as during steady-state conditions. The analytical solution for the sediment discharge performs well when compared with numerical and experimental results. The approximate analytical solution for the concentration profile is the asymptotic limit of the transient numerical solutions. An error analysis shows that the analytical solutions improve with increasing slope length and that the solution model presented here is applicable to a wide range of physical situations.

INTRODUCTION

This study analyzes the problem of rainfall induced surface erosion on steep slopes. Modern highways are often built to very demanding standards, resulting in very steep slopes adjacent to roads. For example, the freshly cut decomposed granite slopes along the new alignment of Route 299 in Shasta county, northern California, present a serious threat in terms of erosion of the hillslope surfaces. Decomposed granite is a cohesionless soil and the bare hillslopes are very prone to erosion when subjected to rainfall. The eroded soil fills up the ditches and sedimentation basins along the highways, plugs the culvert drains and then causes the washing of the road by surface flows. Cleaning up these draining facilities after heavy rainfall events by maintenance crews is very expensive and not a lasting solution. To formulate effective abatement strategies, an understanding of the physical phenomena of soil erosion on such surfaces is essential. The steep nature of the hillslope and the cohesionless soils on the exposed surface make this problem a unique one. Soil erosion by water involves the

0022-1694/91/$03.50

© 1991 - - Elsevier Science Publishers B.V. All rights reserved

280

RS. GOVINDARAJU AND M.L. KAVVAS

detachment of particles and their subsequent transportation by surface water. The two important agents in sheet erosion are rainfall and overland flow. Experimental studies of Moss et al. (1979) indicated that on mild slopes, sediment detachment is primarily the result of the impact action of rainfall. On slopes less than 3.0%, rainfall is the dominating detaching agent and the surface overland flow merely acts as a transporting agent leading to what is called 'rainflow transportation' (see Moss et al., 1979). The sediment is splashed into the air in a random fashion, but as the slopes are mild, there is negligible preferential downstream splash (Walker et al., 1978). Once the overland flow develops sufficiently, the impact of raindrops is cushioned and thus soil removal as a result of impact of rainfall is considerably reduced (Palmer, 1964). As the slope increases, the ability of the overland flow itself to entrain and transport sediment increases. The surface flow exerts a tractive force on the bed, causing the sediment particles on the hillslope to be dislodged and become part of the suspended sediment. For slopes greater than 5.0%, entrainment by overland flow is the dominant eroding mechanism (Moss et al., 1979) used the term 'rheologic flow' for sediment transported by this procedure). For very steep slopes, the influence of rainfall impact in sediment detachment is practically negligible. Moreover, when the soil surface is made up of cohesionless material (such as decomposed granite), the influence of rainfall impact in soil detachment and transportation can be ignored (Kilinc and Richardson, 1973). There are distinct advantages in representing the significant features of a physical process such as erosion over land surfaces in mathematical terms. With the use of digital computers, we are able to simulate many scenarios and then evaluate the resulting consequences or formulate strategies and preventive measures. The first type of such models are the black-box type models. A common example is the use of gross soil erosion equations such as the universal soil loss equation (Wischmeier, 1960; Agricultural Research Service, 1961). Some researchers tried to formulate predictive equations based on watershed parameters (Flaxman, 1972) or regression equations (Kilinc and Richardson, 1973; Leaf, 1974; Megahan, 1974). Bennett and Sabol (1973) related sediment transport to water discharge, and, until recently, such methods have been utilized as predictive models by many public and private agencies. However, as our understanding of the physical process of erosion improves, we try to incorporate this new knowledge into the analysis, leading to a more physics-based modeling effort. Such models provide spatial distribution of unknown quantities and are applicable to ungauged experimental sites. This study falls into the category of such conceptually based mathematical models. Examples of such modeling efforts are the works of Meyer and Wischmeier

MODELING THE EROSION PROCESS OVER STEEP SLOPES

281

(1969) and Rowlinson and Martin (1971), who split up the erosion process into component processes of detachment and transportation. Foster and Meyer (1972) distinguished between sediment transport occurring through rills and through the more shallow overland flow (inter-rill areas). Most of these models have a continuity equation for the conservation of the sediment mass and another equation that relates the sediment load to the flow detachment (called the erosion equation in this paper). Similar ideas, with some modifications, were used by David and Beer (1975) and Negev (1967). Smith (1976) used a dynamic simulation model that involved the kinematic wave approach for surface flow modeling coupled to a sediment transport model. He also used a sophisticated infiltration component. In this study, we use a physics based hydrologic model for surface overland flows (see Govindaraju et al., 1990), which provides analytical solutions for temporal variations in rainfall and infiltration. This is coupled to the erosion model of Foster and Meyer (1972), and analytical solutions for the time-space varying sediment discharges over steep hillslopes are developed. These analytical solutions are tested against numerical and experimental results. A recent review of physics-based models for predicting erosion from overland flows and channels was provided by Foster (1982). MATHEMATICAL F O R M U L A T I O N

There are two distinct physical phenomena involved in the study of erosion and deposition of sediment particles in hillslopes. First we need to be able to analyze the flow dynamics and obtain analytical solutions for the time-space distribution of the flow depth and discharge. These quantities are then used in the prediction of the sediment concentrations and sediment discharges. This approach explicitly assumes that the sediment concentrations in the overland flow regime are sufficiently small so that the flow dynamics are not affected by the suspended sediment and may therefore be treated independently. This assumption uncouples the water flow from the sediment transport process. Of course, if this assumption is not valid (for example, as a result of very high sediment concentrations, as in mud flows), then both the flow and sediment dynamics need to be solved simultaneously. In this study, we shall assume that the flow process is not affected by the suspended sediment. Further assumptions have been made in this study with regard to the sediment transport component. The complicated processes involved in sediment sorting are not considered. The changes in bed resistance and transient effects have been ignored. Gully and rill formation is an important feature on steep slopes. The rill geometry on soil surfaces is a dynamic quantity and its influence is represented only in terms of an average soil erodibility. Better

282

R.S GOVINDARAJU AND M I

KAVVAS

quantitative representation of the rilled surface structure is important for modeling purposes and is currently being studied by the authors. Surface.flow dynamics

The shallow water equations (or the Saint Venant equations) are frequently used as the mathematical representation of flows on planar land segments. These are comprised of the continuity and momentum equations as ~73,

6(V>,)

~3t +

~.,<

Sc =

SO

-

R ( x , t) -- f i x ,

t)

(73'

V ?V

1 ?V

?x

g ?x

g ?t

-~

i(x,t) iV gy

0 < x < L,t

> 0

(1) (2)

Figure 1 is a definition sketch of overland flow on a typical hillslope. Here y(x, t) is the overland flow depth, V(x, t) is the depth-averaged flow velocity, L is the length of the hillslope whose slope is So, i (x, t) is the net lateral inflow into the overland flow section (and equals the rainfall R ( x , t ) minus the infiltration f ( x , t)) and g is the gravitational acceleration. The friction slope (or the slope of the total energy line) Sf is related to the flow depth and velocity through a friction relationship. One of the common relationships used in such instances is the Chezy law: V(x, t) =

C,,/y(x, t) x Sf(x, t)

(3)

where C is a coefficient expressing the roughness of the soil surface. There are no known analytical solutions to the full Saint Venant eqns. ( 1) and (2), and thus researchers have considered simplifications whenever the physical conditions justified them. One such simplification is obtained by neglecting

Rainfall R(x,t) Rainfall Causing Soil Splash

Reference (Datum)

Fig. I. S c h e m a t i c sketch o f o v e r l a n d flow o v e r a hillslope.

283

MODELING THE EROSION PROCESS OVER STEEP SLOPES

all but the first two terms on the right-hand side ofeqn. (2) (see Morris, 1979; Morris and Woolhiser, 1980; Govindaraju et al., 1988), leading to Sf =

0y 0x

So

(4)

Combining eqns. (1), (3) and (4) leads to the diffusion wave approximation

#Y ~ l Cy3!2 (So - -~x #Y )~/2] = i(x,t) #t _ O-Jx

(5)

Govindaraju et al. (1990) used a sine transformation for the spatial part of the solution to obtain analytical expressions for the overland flow and discharges. The choice of this spatial representation of the flow profile is not arbitrary. It was observed by Govindaraju et al. (1988) that the sine functions are the eigenfunctions of the corresponding homogeneous problem while utilizing the diffusion wave approximation for overland flows. We shall use their analysis in this study. The depth of the overland flow profile may be approximately expressed as

y(x,t)

= h(t) sin ~

(6)

The flow discharge at any space and time point is the product of the flow depth and depth-averaged velocity, and is given as

q(x,t)

= Ch 3/2sin 3/2 ~-~

x

SO - h ~ c o s

2L)]

(7)

The initial and boundary conditions which are applicable to steep slopes were discussed by Govindaraju et al. (1988) and are used in this study. These end conditions are

y(x,O) = 0 y(0, t)

Oy(L, t) #x

=

0 <~ x <~ L

0 -

(8a) (Sb)

0

(8c)

It should be noted that the solution model suggested in eqn. (6) satisfies the boundary conditions in eqns. (8b) and (8c). Equation (6) constitutes an approximation to the true solution and, in general, is different from the true solution. Substituting eqn. (6) into the diffusion wave model of eqn. (5) and integrating over the spatial domain, results in the following ordinary differential

284

R.S. G O V I N D A R A J U A N D M.L. KAVVAS

equation for the flow depth at the downstream end (which equals h(t)):

dh 7r CS(i/2h3, 2 d~ + ~-~

g(t)

=

0

(9)

where g(t) is the spatially averaged net lateral inflow given by L

g(t)

-

2Lrr f i ( x , t ) d x

(t0),

0

It is reasonable to assume that there is no spatial variation of rainfall over the scale of a single hillslope. Equation (9) thus uses an average spatial infiltration. If the soil surface is fairly impervious, then there is no approximation in using the representation in eqn. (10). If the soil is pervious, one first needs to calculate infiltration to estimate i (x, t) as a function of rainfall and infiltration; g(t) is then estimated by eqn. (10). However, g(t) is a time-dependent function. This function may be approximated to any degree of accuracy by discretizing it into a series of piecewise constants. Let us say that this value is g, over the time interval (z, r + Az). The solution to eqn. (9) is then obtained as (Govindaraju et al., 1990)

1 log 6a

au + a:

1

(u + a) 2 7~

+ a-~3 tan ' \

a,]3 /

k

CvS(So)(r ..... t)

4L

(11)

where k and a are given by the following expressions: k -

6a

~

a -

rcC\/(So)]

+ a) 2

+ ~ t,a n

" u = u(t) = h , ~

'\

ax/3

ot

(12a)

(12bi

The solutions in eqns. (11) and (12) are applicable over the interval (r. r + At). After the rain stops (and assuming that there is no infiltration), there is no net lateral inflow contribution to overland flow. The solution to eqn. (9) for the receding part of the hydrograph is then given as (Govindaraju et al., 1990)

h(t)

=

I

4Lu~ 12 4L - rcurCx/(So)(t r - t)

(13)

where tr is the duration of rainfall and u, = ~ . The value of ur is therefore determined from eqns. (11) and (12) for the case when the surface is pervious

MODELING THE EROSION PROCESS OVER STEEP SLOPES

285

and there is infiltration. The applicability of the solutions in eqns. (11)-(13) were examined by Govindaraju et al. (1990). Sediment transport dynamics

The detachment and subsequent transport of sediment is strongly dependent of the flow dynamics (even though the reverse is not true for small sediment concentrations). Continuity of sediment mass in the surface flow regime leads to the following flow equation (Foster and Meyer, 1972; Bennett, 1974; Li, 1979): 0

~t (cy) + ~x (cq)

=

(Dr + D,)/ps

(14)

where c(x, t) is the concentration of the sediment by volume, p~ is the mass density of the particles, Dr is the rill erosion rate (mass per unit area per unit time) and D1 is the sediment delivery rate from inter-rill areas (mass per unit area per unit time). The spreading of the sediment owing to dispersion effects are neglected in eqn. (14). The inter-rill erosion rate is frequently expressed as Di

=

(15)

clR":

where c l and c2 are constants and R is the rainfall as defined in eqn. (1). Foster et al. (1977) used a value of c2 of 1.0. Equation (15) suggests that the inter-rill erosion depends on the rainfall intensity. Thus, this equation may also be conceived as incorporating the effects of soil splash as a result of impact action of raindrops. The detachment of particles inside rills depends on the specific flow conditions and is the dominant process of sediment entrainment and transportation in steep hillslopes. It is frequently expressed as (see Foster and Meyer, 1972; Foster, 1982) Dr =

psa(Tc/ps-

cq)

(16)

where Tc is the flow transport capacity (mass per unit area per unit time). This equation suggests that the erosion rate is proportional to the difference between the sediment transport capacity and the sediment load in the flow. Thus the flow has a maximum eroding capacity when it is free of suspended sediment. When the sediment load is greater than the transport capacity of the flow, deposition occurs. Most soils have some resistance to surface erosion, and therefore the overland flow has to develop enough tractive force to overcome this critical shear stress value. The transport capacity is commonly expressed as (Li, 1979; Foster, 1982) Tc =

Ct(TyS o -

zc~)p

(17)

286

R.S. GOV1NDARAJU AN[) M t . KAVVA%

where C t is a coefficient determining the erodibility of the soil as a result of sheet erosion and 7 is the density of water. The term 7ySo is a measure of the tractive force exerted by the surface flow on the soil particles on the bed. The exponent p varies between 1.0 and 2.5. Foster and Meyer (1972) chose p = 1.5 for mathematical convenience and the same value was later used by Li (1979). The critical shear stress, ~cr, in eqn. (17) is very small for cohesionless soils and is often neglected. Foster (1972) suggested a value of ('~ of 0.6 in eqn. (17). The parameter ~ in eqn. (16) behaves like a rate constant and may vary over a wide range depending on the soil type. Substituting eqns. ( 15)-(17) into eqn~ (14) yields the following:

8(cq)

~(cy) Ti -+

--

[ ('~ (7)'5,,-- ~ c , ) " ...... c q ] + (clR'~-)/p~ i,--,

(18)

Equation (18) is the governing differential equation for sediment transport in surface flows. It requires two boundary conditions and one initial condition for the problem to be well posed. The initial condition for the sediment transport must be one of zero concentration, dictated by the dry initial condition in eqn. (8a). Similarly, we choose a zero concentration value at the upstream end of the flow domain to coincide with the zero depth condition of eqn. (8b). The downstream boundary condition for the concentration is taken to be similar to the flow condition at that end. Thus we have the following end conditions for the flow concentration:

c(x,O)

=

0

(19~)

c(0, t) =

0

(19b)

~c( L, t)

-

0

(19c)

For the sake of simplicity, we start our analysis with the steady-state solution where there are no temporal variations. F r o m eqns. (8) and (9), the steadystate value of the flow depth at the downstream end of the flow section is

h(t = ~ )

= 11,~ =

_

iL -:C,,/(S,,)..

(20)

The steady-state flow discharge over the flow plane using the diffusion wave approximation in eqn. (4) follows from eqn. (7) as q~(.v) =

C

C~-(-S0)

sin

~

C~/(S,~).

_~L cos \ ~-~/~ (21)

287

MODELING THE EROSION PROCESS OVER STEEP SLOPES

It should be noted that at x = L, the outflow discharge from the whole plane is iL, as would be expected from continuity considerations. We let Y represent the sediment volume discharge (i.e. Y = cq). The steady-state equation for sediment discharge is then obtained as --dx + o-Y

-

7Soh~s sin

+ oP~

rex

-

rcr

(22)

P~

Using the integration factor method, the analytical steady-state solution is obtained as Yr(x)

-

c l R ~2

(1 - e ..... ) + e ....

O-P,;

o-C t

| e ~ [7S0h~s sin ( ~ ) Ps :o

- r~,~]"dx (23)

This is the analytical solution for the steady-state sediment discharge over the whole plane. The subscript p of the left-hand side indicates that the form of the solution, Y ( x ) , depends on this value. The integral on the right-hand side cannot be presented analytically for fractional values ofp. It may be evaluated numerically to any desired degree of accuracy by using any quadrature formula (see, for example, Abramowitz and Stegun, 1968). However, implicit analytical solutions are available for integral values ofp. Thus, for p = 1, the solution to eqn. (22) is expressed as Y~(x)

_

c_ l_R " : (1 -

o-p,;

+

e ~") +

Ctrcr(e

p,;

~' - 1 )

o-C~ ~'S0&,; [ (rrx) = p,~ o-2 + (rc/2k)2 o-sin ~ - ~cos

(rt.v) ,,. rt ] Z + e (24)

This equation clearly shows the influence of raindrop impact, slope and length of overland flow plane, the critical shear stress and o- on the steady-state sediment discharge. The influence of the surface roughness is introduced through the steady-state depth at the outflow section, hs,;, as given in eqn. (20). In many physical situations, the value of L is large (note that L is a macroscopic value). Many hillslopes have values of L greater than 40 m. Equation (24) may be considerably simplified as the value of L tends to infinity asymptotically. Thus Yl(x) -

_clR"2 _ (1 -- e "~) + Ct-ccr . (e ~_ -- 1) + -C- t ;'Soh~ sin ( g x )

(25)

28g

R.S. GOVINDARAJU AND M . [ KAVVAS

where Yr(x) is the asymptotic value of Yp(x) as L tends to infinity. The third term on the right-hand side of eqn. (25) is the most important term when considering slopes composed of cohesionless soils such as decomposed granite. This term incorporates the influence of the tractive forces exerted by the overland flow in entraining and subsequent transportation of sediment particles. The second term on the right-hand side shows the influence ol" critical shear stress and may be neglected for cohesionless soils (Foster and Meyer, 1972). The first term on the right-hand side of eqn. (25) includes the effect of rain splash and inter-rill erosion and is also often neglected. When the exponent p = 2 in eqn. (22), the explicit solution after simplification is obtained as --

c,p,

(1

....

~')+-

e

_2,<

,s0t,,

+ p~[o-2+ 0z/Lf]

+~

~

.

(I

p~

-

e

[

cos

sin ( ~ - / ) [ a s i n

....

- ~ cos

+

,

,e0, o-

The expression for steady sediment discharge in eqn. (26) for p - 2 is much more complicated than eqn. (24) for p = 1. Again, as L is likely to be fairly large, the asymptotic expression for very large L is presented: -

(1

-

- 2L.,./~-

e

~') + -

7Soh. sin ~

(I

- e

+

P~

(27)

There are certain similarities in eqns. (25) and (27). The fourth term on the right-hand side of eqn. (27) is the d o m i n a n t term for hillslopes composed of cohesionless soils, and incorporates the effect of the tractive force exerted by the overland flow. The terms including the critical shear stress and rain intensity are relatively less important. Foster and Meyer (1972) used only the effect of overland flow in soil detachment in their analysis. In eqns. (25) and (27), the increase of critical shear stress leads to a reduction in the steady sediment discharge (as would be expected). It is not possible to write explicit analytical solutions (like eqns. (24)-(27)) when p is not

289

M O D E L I N G THE EROSION PROCESS OVER STEEP SLOPES

an integer. The form of eqns. (25) and (27) suggest the following relationship for steady-state sediment discharge for large L:

Y (x) -

c l R '2

--

Ct -~- - -

Ctr 7S0hss sin ( ~ x ] l ~

(1 - e -"x) + e ..... )

P(1

~'cr

L

-

-

\ 2L Jl

C¢

-

-

P'rc,- - - 7Sohss sin

/\[~x]

(28)

Equation (28) is not exact and has been obtained through inspection of eqns. (25) and (27). It reduces to eqn. (27) for p = 2. It should be noted that the approximation in eqn. (28) is restricted to terms involving the critical shear stress. As it has been argued that the critical shear stress has a small influence, eqn. (28) is likely to be a good alternative analytical expression for the steady-state sediment discharge. The usefulness of this model will be evaluated later in this study. The depth-averaged steady-state concentration profile can easily be obtained with any of eqns. (24)-(28) by using the following relationship: %(x)

=

Y(x)/qss(X )

(29)

where qss(X) is obtained from eqn. (21). We shall be dealing with depthaveraged concentrations in this study. In general, it is possible for the concentration to be dependent on the depth. However, overland flows are very shallow by nature and concentrations in the vertical direction are not important. It is very difficult to sample depth-dependent concentrations because of the thin nature of the flow. Almost all researchers, working in the laboratory or in the field, report depth-averaged concentrations. Moreover, rain splash would cause further mixing in the vertical direction, lending more credence to the depth-averaged approach. In deeper flows (such as channels and estuaries), it is conceivable that concentrations can vary very significantly with depth (Foster (1982) and Li (1979) dealt briefly with channel problems). The analytical expressions for the transient solutions are based on the premise that the water and sediment discharge have the same relationship at all times as they do at steady state. This implies that the sediment discharge depends on the flow discharge in a similar manner at all times. The sediment discharge is therefore assumed to rise and fall with the water discharge. The outflow hydrographs and sedimentographs have their peaks at the same time with no time lag. This assumption is supported by observed water and sediment discharge profiles (see Lane and Shirley, 1982). Therefore, the asymptotic time-space dependent sediment discharge for large L is

290

R.S. GOVINDARAJU AND M.L. KAVVAS

given by V (x, f)

-

clR':

(l

--e

C;"f)

~") + P~ 75oh(t) sin ~,2L./J

C~ ('t (rex) --- Y,ir(l .... e ~;') -- pr<,.- 7Soh(t) sin

(30)

where the value of h(t) in the above equations is obtained from either eqn. ( 11 or (13), depending on the nature of the rainfall function. The depth-averaged time-space dependent concentration can then be determined as

c(x,t) = Y(x,t)/q(x,t)

(31)

Thus we are able to obtain analytical solutions for the flow and sediment transport over steep hillslopes. In the following sections, we examine the validity and applicability of this model by comparison with numerical and experimental results. DISCUSSION OF RESULTS

Before we consider the usefulness of the analytical solutions developed in the previous section, we first examine the ability of these solutions to replicate physical situations. We first compare the performance of these solutions with those obtained from experiments conducted by Singer and Walker (1983). Their experimental set-up consisted of a laboratory flume 0.55 m by 3.0 m in which 200 kg of moist fresh soil (bulk density 1200 k g m 3) was packed to 80 mm thickness for each erosion event. Water was applied as rainfall (through the use of a rainfall simulator), overland flow or a combination of both. During their experiments, their flume was placed at 9% slope and the duration of each run was 30.0min. This allowed sufficient time for steady state to develop without causing major changes to the bed surface, and also gave the bed a chance to react to the rainfall and overland flow treatments (see Singer and Walker, 1983). Figure 2(a) shows the observed and analytically computed outflow hydrographs for the two different cases of rainfall intensity. In these experiments, approximately 20% of the simulated rainwater seems to have infiltrated. The agreement between the experimental results and the analytical solutions presented in Fig. 2(a) is reasonably good. Figure 2(b) shows the corresponding sedimentographs for the hydrographs shown in Fig. 2(a). The sediment discharge curves are obtained from the experimentally observed values of Singer and Walker (1983) and numerical solution of eqns. (18) and (19). A centered implicit finite difference scheme has been used for all the numerical solutions presented in this study. Figure 2(b) shows that eqns. (18)

MODELING

100,0

THE

EROSION

PROCESS

OVER

STEEP

291

SLOPES

-

80,0'

_-

--_-_

_-

-

~, - -

Numerical (100)

60,0' o O

40.0'

•

Observed (100)

•

Numerical (50)

•

Observed (50)

r¢ 20.0

0,0 0.0

5.0

10.0

15.0

20.0

25.0

30.0

T i m e (mins)

(a) 60.0-

50.0

rV

40.0

........ Numerical (100)

3O.O o t~

20.0 _

----_~

--

~.

~.

_----

~.~.~..~.

•

Observed (1001

-

Numerical (50)

•

Observed (50)

- -

10.0 bets indicate rain intensity in mm/hr) 0.0

,,i--"

0.0

(b)

, ,

!

5.0

.

.

.

.

i

.

10.0

.

.

.

n

.

15.0

.

.

.

i

.

20.0

.

.

.

u

. . . .

25.0

l

30.0

Time (rains)

Fig. 2. C o m p a r i s o n of (a) n u m e r i c a l outflow d i s c h a r g e h y d r o g r a p h s and (b) n u m e r i c a l s e d i m e n t o g r a p h s with the e x p e r i m e n t a l l y observed results of Singer a n d W a l k e r (1983).

and (19) are reasonable representations of the erosion process. In Figs. 2(a) and 2(b), the length of the flow plane is only 3 m. The analytical results are applicable for longer flow sections and are not considered for the cases in Fig. 2. Figs 3(a)-3(d) show the numerical and analytical results of this study and the observed values of Kilinc and Richardson (1973), who made 24 experimental runs and measured water and sediment discharges under simulated

292

R.S. GOVINDARAJU AND M.L. KAVVA.," Sediment discharge vs time for rain intensity = 3.65 in/hr

~. 20.0 16.0' Numerical (15)

12.o o ,.m

8.0'

~3

¢"

Analytical (15)

•

(Yoserved ( I 5 )

•

Numea'ical (20)

•

Analytical (20)

4.0"

Ot~,wexl (20)

0.0 . . . . 0.0

(Numbers

in tnackets indmat~ slope in percent)

| ....

i ....

10.0

20.0

t ....

i ....

i ....

30.0 40.0 Time (rain)

50.0

n

60.0

(a)

Sediment discharge vs time for rain intensity = 3.65 in/hr

-~ 5 0 " 0 ~ , - - - - - ~ -

-

=-:

= :

777=

3~

~

I

a

. f

0.01 0,0

(Numbcls

Ntmmiscal(30) ,~ytieat (30) Otm'ved(30) Nm',n'imd(40) Analytica(40) Otmtved(40)

in brackets indicatc slope in pctcem)

10.0

20.0

30.0 40.0 Time (rain)

50.0

60.0

(b)

Sediment discharge vs time for rain intensity = 4.60 in/hr ~. 40.0

* 30.0 ~

~

. . . . . . . . . . . . . .

!

Nunmical(20)

IO.O

• (Numbers

0.0n

o.o

Numerical(15)

, , , . m ....

10.0

in Mackcts indicate sl n ....

20.0

I ....

o ....

30.0 40.0 Time (rain)

analytical(20)

---o--- Observed(20)

in percent) I i t I I I

50.0

60.0

(c)

Fig. 3. N u m e r i c a l , a n a l y t i c a l a n d e x p e r i m e n t a l l y o b s e r v e d ( K i l i n c a n d R i c h a r d s o n , 1973) results. (a) R a i n intensity, 3 . 6 5 i n h -I . Slopes, 15 a n d 2 0 % . (b) R a i n intensity, 3 . 6 5 i n h - I . Slopes, 30 a n d 4 0 % . (c) R a i n intensity, 4 . 6 0 in h-~'. Slopes, 15 a n d 2 0 % . (d) R a i n intensity, 4.60 in h L . Slopes, 30 a n d 4 0 % .

MODELING

THE

EROSION

PROCESS

OVER

STEEP

293

SLOPES

Sediment discharge vs time for rain intensity = 4.60 in/hr 80.0 S .

70.0 ~

. . . . . . . . . . . . . . . .

60.0

~5o.o

Numerical (30) 4

"~ ~ 40.0

~

~

.

~

_~

~e-

~q

-

.

--

.;

_

-

.a

a 7.

• ~.

30.0

---o---

10.0

Observed (30) Numerical (40)

•

~ 20.0

Analytical (30)

Analytical (40) Observed (40)

(Numbers in brackets indicate slope in percent)

0.0

.... 0.0

, .... 10.0

(d)

, .... 20.0

, .... 30.0

, .... 40.0

, .... 50.0

, 60.0

Time (min)

Fig. 3. (continued).

rainfall. The sediment density used in that study was 93.61bft 3. The laboratory flume was 4 ft high, 5 ft wide and 16 ft long, with adjustable slope. Six bed slopes were investigated (5.7, 10, 15, 20, 30 and 40%). Four different rainfall intensities were simulated (1.25, 2.25, 3.65 and 4.60 in h ~). A constant infiltration rate was assumed for each run. Each experiment was carried out for a total duration of 1 h, with sediment concentration data being collected at 5-10 min intervals. Kilinc and Richardson (1973) used a Darcy-Weisbach friction relationship, which is similar to the Chezy relationship used in this study. F r o m their data, Chezy's C was estimated as C = 15.0ftl/2s 1. Li et al. (1977) estimated the value o f p as 2.5 and the Ct value as 65.2 for the results of Kilinc and Richardson (1973). In this study, we have used p = 2.5 and C t = 40.0. The value of o- was chosen as 10.0 ft -~ , which is within the range of values suggested by Foster (1982). In Fig. 3(a) comparisons are made for cases of 15 and 20% slopes with simulated rainfall intensity o f 3.65 in h -1 . The sediment discharge values are presented for three different methods. The numerical results are obtained through the numerical solution of eqns. (18) and (19). The analytical solutions are obtained through eqn. (30). The observed results are evaluated from the data of Kilinc and Richardson (1973). There seems to be some spread in the observed values but the numerical and analytical results are able to simulate observed values well. Figure 3(b) shows results for slopes of 30 and 40% and rain intensity o f 3.65inh -~. Figure 3(c) and 3(d) show similar results for rainfall intensity o f 4.60 in h -j . Figure 3(c) considers instances of 15 and 20%

294

R.S. G O V I N D A R A J U

A N D M.L. KAVVAS

Steady state sediment discharges over the hiUslope Length of slope = 20m, Chezy's C = 15.0, Slope = 10%, Rain intensity = 50 m m / h r 100.0

-

80.0 "~ ~

60.0

8

40.0 E

•

20.0

Exact Lo = 1.0) Approx. (p = 1,0) Exact (p = 1,5) Approx. (p = 1.5) Exact (p = 2,0) Approx. (1a = 2.0)

0.0

-20.0 0.0

4,0

8.0

12.0

16.0

20.0

Distance along the slope (m) Fig. 4. Numerical and analytical steady-state sediment discharge profiles for p = 1.0, 1.5 anti 2.(].

slopes, and the slopes are 30 and 40% in Fig. 3(d). The results in these figures show that the model developed in this study may be used to predict erosion caused by overland flows. A n o t h e r feature of interest in Figs. 3(a)-3(d) is that the analytical solutions approximate the numerical solution very well. The analytical solution is supposed to perform well for long slopes, and the study of Kilinc and Richardson used a 16 ft flume. This suggests that the analytical solution presented in eqn. (30) is a good alternative to the numerical solution of eqns. (18) and (19). Equations (24) and (26) are the analytical steady-state solutions fbr p = 1 and p = 2 respectively. Figure 4 shows the shape of these two curves for a particular slope of length 2 0 m whose surface roughness is represented by C = 15 m 1/2 s i, and the slope is subjected to rain intensity of 50 m m h ~. The soil erodibility is set to o- = 10m ~ ~, Ct = 0.6, and the slope of the soil surface is 10%. The exact solutions for p = 1.0 and 2.0 in Fig. 4 are obtained from eqns. (24) and (26) respectively. The approximate solutions in Fig. 4 for p = 1.0 and 2.0 are obtained from the asymptotic solutions for large slope lengths presented in eqns. (25) and (27) respectively. The exact solutions for p = 1.0 and 2.0 are very well matched by the approximate solutions for p = 1.0 and 2.0. Analytical solutions are available only for integral values of

295

MODELING THE EROSION PROCESS OVER STEEP SLOPES

p. The exact solution f o r p = 1.5 in Fig. 4 represents the steady-state solution obtained from the numerical solution of eqn. (22) subjected to the boundary conditions of eqns. (19b) and (19c). The approximate solution for p = 1.5 is obtained from eqn. (28). The analytical approximation of eqn. (28) is a very good alternative to the numerical solution for p = 1.5 and these two solutions are practically identical in Fig. 4. Therefore the approximate Numerical sediment discharge profiles Length of slope = 20m, Chezy's C = 15.0, Slope = 10 %, Rain Intensity = 50 mm/hr 40.0 ]

35.° 1

1

Time = 0.25 min Time = 1.0 min • Time = 2.0 min Time = 3.0 min • Time = 4.0 min --------o--- Time = 5.00 min

25.0 "I

1

"

15.0 -4 .~

10.0 5.0 0.0 0.0

4.0

8.0

12.0

16.0

20.0

Distance along the slope (m)

(a)

Analytical sediment discharge profiles Length of slope = 20m, Chezy's C = 15.0, Slope = 10%, Rain intensity = 50 mm/hr 40.0 35.0 30.0 Time = 0.25 min Time = 1.0 min Time = 2.0 min • Time = 3.0 min • Time = 4.0 min .------n--- Time = 5.00 min

25.0 20.0 15.0 10.0 5.0 0.0 0.0

4.0

8.0 12.0 Distance along slope (m)

16.0

20.0

(b) Fig. 5. Transient (a) numerical and (b) analytical sediment discharge profiles for the rising phase.

296

R.S. GOVINDARAJIJ AND M.L. KAVVAS

analytical expressions for the steady state developed in this study may be used for slopes longer than 20 m. Figure 4 shows the remarkable variation in the steady-state profiles for different values of p. There is, however, no physical way of estimating p (which is a function of soil type, and frequently is a calibration parameter). Foster and Meyer (1972) used p = 1.5 for mathematical simplicity. Figure 4 shows the usefulness of the analytical steady-state solution in predicting surface erosion. We now investigate the ability of the approximate transient analytical solution in predicting the sediment discharge profiles at various times. Figure 5(a) shows the rising sediment discharge profiles over a 20m slope. The surface roughness is expressed by C = 15mJ/2s ~ p = 1.5, the critical shear stress was chosen as 0.00t 5 kg m 2 and the other parameters are the same as in Fig. 4. The results in Fig. 5(a) were obtained by the numerical solution of eqns. (18) and (19). This figure shows that the steady state is reached in about 5.0 min for the physical situation considered in this example. (It should be noted that there is hardly any difference between the profiles at times 4.0 and 5.0 min.) Figure 5(b) shows the sediment discharge profiles obtained from the analytical solution described in eqn. (30). The times for these sediment profiles coincides with those presented in Fig. 5(a),

Sediment concentrationprofiles Length of slope = 20m, C h e z y ' s C = 15.0, Slope = 10%, Rain intensity = 50 m m / h r

120.0 :

. . . . . .

~ 100.0 90.0

8°°t

• 8

/

.ollf! o011/

70.0t/t["/ / 50.0

•

30;0

A

Time = 0.25 min Time = 0.50 rain Time = 0.75 min Time = 1.0 rain Time = 3.00 min Time = 5.00 min Theoretical

20.o|1 10.0

-f

0.0 i [ 0.0

.

.

. 4.0

.

. 8.0

,

. 12.0

16.0

20.0

Distance along slope (m) Fig. 6. T r a n s i e n t s e d i m e n t c o n c e n t r a t i o n profiles for the rising phase and the a n a l y t i c a l sediment concentration at steady state.

297

M O D E L I N G THE EROSION PROCESS OVER STEEP SLOPES

showing the numerically obtained results. Figures 5(a) and 5(b) are very similar, implying that the analytical solution presented for the transient sediment discharge is a possible alternative to numerical solutions when large slope lengths are under investigation. Figure 6 shows the numerically obtained concentration profiles at various times for the case considered in Figs. 5(a) and 5(b). It may be observed that the sediment concentration reaches its steady-state value within half a minute at the downstream end. With increasing time, the concentration values at the downstream end do not change while the concentration approaches the uniform steady-state value over increasing portions of the upstream sections. As time progresses, the concentration over the entire flow plane seems to be nearing the steady-state value. The theoretical solution shown in Fig. 6 is obtained from eqns. (28) and (29). The transient concentration predicted by eqn. (31) remains practically unchanged over time and space and therefore only a single theoretical profile (that of the steady-state soution in eqns. (28) and (29)) is shown in Fig. 6. The numerical solution asymptotically reaches the analytical solution. Figure 6 suggests that the analytical models presented in this paper are better at predicting the sediment discharge than at predicting the sediment concentrations when compared with numerical results. The steep nature of the concentration profiles at x = 0 is the result of the imposed boundary condition of zero concentration at the upstream end (see eqn. (19b)). From the nature of the concentration profiles in Fig. 6, it may be concluded that the rise of the sediment discharge profiles in Figs. 5(a) and 5(b) is the result of the water discharge while the concentration remains essentially unchanged.

Outflow hydrographs and sedimentographs 20m, Slope = 10%, Rain imensity = 50ram/hr. Rain duration - 10.0 rain

Length of slope 400

"

~ 35.0 .

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

30/)

?-' ~ 25.0-

=~2o.o~ 15.0" :~ I0.0i ~: u~ 5.0" O 0

.

O0

.

.

.

.

,

3.0

.

.

.

.

.

,

.

.

.

.

,

q [)

60 Time

.

.

.

.

.

,

120

.

.

.

.

.

,

15.11

(rain,

Fig. 7. R i s i n g a n d r e c e d i n g l i m b s o f t h e f l o w a n d s e d i m e n t d i s c h a r g e s at t h e o u t f l o w s e c t i o n .

298

R.S. GOVINDARAJU AND M I . KAVVAS;

An important descriptor of the erosion process occurring over a planar land section is the sedimentograph. This graph determines the total amount of sediment lost over the whole hillslope section as a function of time. It also determines the amount of sediment released to downstream sections or receiving waters. For the analytical solution developed in this study to be a useful modeling tool, it should be able to match this feature of the solution process. Some evidence of the analytical model's capability has already been indicated during model validation when comparing with the laboratory results of Kilinc and Richardson (1973) in Figs. 3(a)-3(d). Figure 7 shows the flow and sediment discharge at the outlet of a 20 m hillslope. This figure shows the outflow hydrographs and sedimentographs for the case considered in Figs. 5 and 6. The rain intensity is 50 m m h - ~and the rain duration is 10.0 min. The rising and recession limbs of both the water and sediment flow curves are shown in Fig. 7. The water discharge rises initially and reaches a steady-state value by 5.0min, and once the rain stops, at 10.0min, the water flow recedes very rapidly. The numerical solution in Fig. 7 refers to the numerically obtained sediment outflow discharge from eqns. (18) and (19). From eqn. (30), the approximate analytical solution for the outflow sedimentograph is obtained as Yv(L,t) -

c_ l _R '2 ap~

-

C't (1 - e ~L) + __[TSoh(t)]p p~

pT~r - -

P~

7Soh(t)

+

-Ct- r ~ ( 1 - e

c,L )

p~

(32)

This equation is the analytical sediment discharge shown in Fig. 7. This figure shows that the numerical and analytical solutions are practically identical, thus establishing the validity and usefulness of the approximate analytical solution (eqn. (32)) in prediction of surface erosion. It has been argued that the critical shear stress and rainfall impact arc negligible when considering steep slopes composed of cohesionless soils. From eqn. (32), the sedimentograph is proportional to he, where h is obtained from eqn. (11) or (13). The sedimentograph rises at the start of the simulation, achieves steady state in 5.0 rain, and once the rain stops, at 10.0 min, falls very rapidly. This behavior is very similar to the outflow hydrograph. The rapid drop in the sedimentograph at the cessation of rain is owing to the fast recession of the hydrograph. Neglecting the critical shear stress and rainfall impact in eqn. (32), a simple calculation yields a value of 35.164gs ~m for the steady-state sediment discharge. This value is in good agreement with the steady-state value achieved by both the analytical and numerical sedimentographs in Fig. 7.

MODELING

10000'

THE EROSION

299

PROCESS OVER STEEP SLOPES

Error in steady state analytical solution at outflow section Slope = 10%, Chezy's C = 15.0, Rain intensity = 50 mm/hr

1000" 1oo.

~

tO"

~

1

Sigma = 00.1 Sigma = 00.5 Sigma = 01.0 Sigma = 10.0

.1 .01

Sigma = 40.0 . . . . . . .

!

. . . . . . .

10 Length of slope (m)

(a)

100

i

100

Error in steady state analytical solution at outflow section Slope = 10%, Chezy's C = 15.0, Rain intensity = 50 mm/hr

I0"

+ .1

,o •

s

.01

p=l.0 p = 1.5 p=2.0 p=2.5

. . . . . . .

= ~

i

~

. . . . . .

l0 Length of slope (m)

100

(b) Error in steady state analytical solution at outflow section 100

Slope = 10%, Chezy's C = 15.0, Rain intensity = 50 mm/hr

10

+ g.

.1

s

Ct = 00.60

•

Ct = 06.00

Ct = 60.00

.01 (c)

. . . . . . .

, 10

. . . . . . .

,

100

Length of slope (m)

Fig. 8. E r r o r v a r i a t i o n as a function of slope length for v a r i o u s (a) a (sigma) values, (b) p values+ and (c) C~ values.

300

R.S. G O V I N D A R A J U AND M.L,. KAVVAS

After comparing the analytical solution with numerical and observed results, we now address the issue of its range of applicability. As the development of the analytical solution was based on some simplifying assumptions, some error is introduced when using the analytical solution. For the sake of comparison, we assume that the numerical solution is 'exact' in all cases and compute the percentage error as a measure of the deviation of the analytical solution from the numerical solution. We use the steady-state value at the outlet section as an indicator of the error in the following analysis. Figure 8(a) shows the error percentage as a function of slope length for different values of the soil erodibility ~ (written in the figure as sigma). The other parameters are as in Figs. 5 and 6. Figure 8(a) shows that the log of the error decreases with increasing o- and that the log of the percentage error decreases linearly with increasing log of hillslope length. This decrease would be expected as the analytical solutions are developed on the premise of asymptotically large slope length. Figure 8(b) shows the decrease in error with increasing hillslope length for various p vlaues. The four p values considered are 1.0, 1.5, 2.0 and 2.5. The sensitivity of the results to p value has already been indicated in Fig. 4. F r o m Fig. 8(b), it is not possible to find a consistent relationship between the error and p. For small hillslope lengths, the error is smallest when p = 2.5, but, when the hillslope increases beyond 25 m the error is smallest for p = 1.0. Figure 8(c) shows the decrease in error percentage for increasing hillslope length for various Ct values. From this figure, it appears that the error is not particularly sensitive to changes in C,. Figures 8(a)-8(c) may be used to determine the level of error when using the analytical solutions for various physical parameters. Generally, an error of 5% or less would be acceptable in most physical situations. The three parameters chosen in Figs. 8(a)-8(c) are not directly measurable in the field or the laboratory, and are essentially treated as calibration parameters. The results in Fig. 8 indicate that the analytical solutions are valid for a wide range of parameters with acceptable error, especially when the slope length is greater than 20 m. SUMMARY AND CONCLUSIONS In this study, a physics based surface overland flow model has been coupled with an erosion, transportation and deposition model to study the phenomena of rainfall-runoff-erosion over steep hillslopes which are primarily composed of cohesionless soils. Analytical solutions to overland flow have been presented after Govindaraju et al. (1990). These solutions yield the time-space behavior of the water depth and discharge over the overland flow section under temporal variation in rainfall. The lateral inflow function is averaged

MODELING THE EROSION PROCESS OVER STEEP SLOPES

301

over the hillslope. However, at the scale of 20-100m, the rainfall function is likely to be spatially constant and temporal variations are more important. The rainfall-runoff model is assumed to be essentially unaffected by the presence of sediment and is treated as being independent of the erosion process. The erosion model, however, needs the time-space behavior of the flow, as the overland flow not only exerts the tractive force on the surface, causing the soil particles to dislodge, but is also the medium of sediment transport. The erosion process is conceived as occurring in three stages - - dislodging of soil particles from the bed, entrainment and transportation of the soil, and subsequent deposition of the excess sediment load. All three stages are very dependent on the surface flow status and therefore the erosion process cannot be treated independently of the flow process. There is no perfect erosion equation as the process is not yet completely understood at the microscopic scale, but the first-order reaction model in eqn. (16) has some physical basis and has been used successfully by earlier researchers (see, for example, Foster and Meyer, 1972). The model used in this study for the erosion process assumes that as soon as the transporting capacity of the flow is less than the sediment load, this excess is deposited on the flow bed. It would be useful to know where deposition starts in physical situations. If we ignore the effects of critical shear stress and raindrop impact, then for deposition to occur, Dr = 0 in eqn. (16), which implies that d Y/dx = 0 in eqn. (22). Using the approximate analytical solution in eqn. (28) and solving for d Y/dx = 0 yields the solution x = L, i.e. the downstream end of the flow section. This implies that there is no deposition anywhere along the slope as long as the slope is straight (not concave). With increasing distance along the slope, the sediment load in the flow increases and there is a corresponding increase in the volume of water and the carrying capacity of the flow (and also the stream power or tractive force) moving down as runoff. For the deposition of sediment to occur, the flow must slow down and lose some of its transporting capacity. This can occur either because of some obstruction in the flow path, or more often, because of a change in the slope. In case of concave hillslopes, the slope of the surface becomes milder with increasing distance from the upstream end, causing the flow to slow down and deposit some of its load. This phenomenon has been noted in studies of other researchers (Foster and Huggins, 1977). Many hillslopes adjacent to roads in northern California are very steep and straight. The hillslope ends abruptly at the road, which is comparatively flat. Thus, there is a drastic reduction in the slope at the foot of the hillslope. The flow slows down appreciably when it reaches the road and it has practically no transporting capacity at this point. Thus almost all the sediment load from the hillslope is deposited on the road (this phenomenon was observed by the

302

R,S GOVINDARAJ[) AND M I

KAVV~S

authors on many of the hillslopes in the mountain regions in Shasta County in northern California). This causes a serious hazard to traffic and results in maintenance costs of millions of dollars each year. This analysis suggests that, .just as on straight slopes, no deposition would occur on convex slopes, where the slope is always increasing, which results in an increase in the transporting capacity. In such situations, the carrying capacity of the flow is not likely to be a limiting factor in the erosion process. The analytical solutions presented in this paper are applicable tbr long slopes (L > 20 m). The steady-state analytical solution is found to be a good alternative to the numerical solution. The approximate transient analytical solution is derived on the basis that the relation between the sediment flow and water discharge remains the same at all times as during steady state. The analytical steady-state solution matches the numerical solution very well (see Fig. 4). The sedimentograph is also well represented by the transient analytical solution as in Fig. 7. The sediment discharge profiles from the numerical and analytical solutions are very similar in Figs. 5(a) and 5(b). The concentration profiles are, however, not very well represented by the analytical model towards the upstream end of the flow section (see Fig. 6), The analytical solution predicts an almost constant concentration value over the whole section (Fig. 6). It is the limit of the numerical transient solutions. The primary reason for this discrepancy in the concentration profiles is the upstream boundary condition, which prescribes a zero concentration (because of zero depth) at the upstream end, resulting in a steep jump in both the numerical and analytical solutions. The analytical solution has a wide range of applicability, as seen in Fig. 8. Figures 8(a)-8(c) show that the discrepancy in the analytical and numerical solutions reduces with increasing slope length. The analytical solutions are based on large slope length. Figure 8(a) shows that the error is rather large for o < 0.5 and the slope length would have to be much greater than 100m for the error to be within tolerable limits. As o increases, the error reduces very fast for a given slope length. The variations in the error percentage for different p and C~ values are shown in Figs. 8(b) and 8(c). The variables o , C, and p influence the erosion process only and not the flow process. In this study, all the suspended sediment is lumped together. In some studies (see, for example, Li, 1979), the total sediment is divided into several size classes, and a continuity equation similar to eqn. (14) is written for each of the size classes. The solution process then follows along similar lines, and the total sediment concentration is obtained by summing the concentrations of each of the individual component classes. However, for many practical applications, the total sediment load is the desired quantity. Also, subdividing

M O D E L I N G THE EROSION PROCESS OVER STEEP SLOPES

303

the sediment into size classes increases the number of parameters that are not measurable and need to be obtained through some form of calibration procedure. The analytical model presented here has the limitation that the surface slope remains constant throughout the simulation. In actuality, the slope constantly changes in time and space as the sediment material is eroded or deposited over the flow plane. However, at the hillslope scale it may be assumed that the macroscopic slope of the entire hillslope section does not change appreciably. This assumption was used in this study. We have also neglected the influence of canopy and vegetation cover in this study (see Singer and Walker, 1983) and have restricted our attention to erosion on bare slopes. In conclusion, we have presented a simple analytical model for prediction of sediment discharge from hillslopes subjected to time-varying rainfall events. The analytical solutions are applicable to a wide range of physical parameters. It is hoped that such analytical solutions provide more insight into our understanding of the physical process of erosion over hillslopes and help us evaluate the relative significance of the parameters that influence the rainfall-runoff-erosion process.

ACKNOWLEDGMENTS

This work was supported by California Department of Transportation under Agreement RTA-65K207. This support is gratefully acknowledged.

REFERENCES Abramowitz, M. and Stegun, I., 1968. Handbook of Mathematical Functions. Dover, New York, pp. 875-924. Agricultural Research Service, 1961. A universal equation for predicting rainfall-erosion losses. US Dept. Agric., Washington, DC, Rep. ARS 22-66, 11 pp. Bennett, J.P., 1974. Concepts of mathematical modeling of sediment yield. Water Resour. Res., 10(3) 485-492. Bennett, J.P. and Sabol, G.V., 1973. Investigation of sediment transport curves constructed using periodic and aperiodic samples. Proc. Int. Association for Hydraulic Research Int. Symp. on River Mechanics. Asian Institute of Technology, Bangkok, Vol. 2, pp. 49-61. David, W.P. and Beer, C.E., 1975. Simulation of soil erosion. Part I. Development of a mathematical erosion model. Trans. ASAE, 18(1): 125-133. Flaxman, E.M., 1972. Predicting sediment yield in western United States, J. Hydraul. Div. ASCE, 98(HY!2): 2073-2085. Foster, G.R., 1982. Modeling the erosion process. In: C.T. Haan, H.P. Johnson and D.L. Brakensiek (Editors), Hydrologic Modeling of Small Watersheds. ASAE 295-380.

304

R.S. GOVINDARAJU AND M.L. KAVVAS

Foster, G.R. and Huggins, L.F., 1977. Deposition of sediment by overland flow on concave slopes. In: Soil Erosion: Prediction and Control. Special Publ. Soil Conserv. Soc. Am., 21: 167-180. Foster, G.R. and Meyer, L.D., 1972. A closed-form soil erosion equation for upland areas. In: H.W. Shen (Editor), Sedimentation Symp. to Honor Prof. H.A. Einstein, pp. 12.1 12.9, Fort Collins, CO. Foster, G.R., Meyer, L.D. and Onstad, C.A., 1977. An equation derived from basic erosion principles. Trans. Am. Soc. Agric. Eng., 20(4): 678-682. Govindaraju, R.S., Jones, S.E. and Kavvas, M.L., 1988. On the diffusion wave modeling of overland flow 1. Solution for steep slopes. Water Resour. Res., 25(5): 734-744. Govindaraju, R.S., Kavvas, M.L. and Jones, S.E., 1990. Approximate analytical solutions for overland flows. Water Resour. Res., 26(12): 2903-2912. Kilinc, M. and Richardson, E.V., 1973. Mechanics of soil erosion from overland flow generated by simulated rainfall. Hydrology Papers, Colorado State University, Fort Collins, Paper 63, 54 pp. Lane, L.J. and Shirley, E.D., 1982. Modeling erosion in overland flow. Proc. Workshop on Estimating Erosion and Sediment Yields on Rangelands, Arizona. US Dept. Agric., Agric. Res. Serv., Agric. Rev. and Manuals, ARM-W-26: 120-28. Leaf, C., 1974. A model for predicting erosion and sediment yield from secondary forest road construction. Rocky Mountain Forest and Range Exp. Sta., US Dept. Agric., Ft. Collins, CO, Forest Service Res. Note RM-274. Li, R.M., 1979. Water and sediment routing from watersheds. In: H.W. Shen (Editor), Modeling of Rivers, Wiley, New York, pp. 9.1-9.88. Li, R.M., Simons, D.B. and Carder, D.R., 1977. Mathematical modeling of soil erosion by overland flow. In: Soil Erosion: Prediction and Control. Special Publ. Soil Conserv. Soc. Am. 21: 210-216. Megahan, W.F., 1974. Erosion over time on severely distributed granitic soils: a model. Intermountain Forest and Range Exp. Sta., US Dept. Agric., Ogden, UT. Foresl Service Res. Paper INT- 156, 14 pp. Meyer, L.D. and Wischmeier, W.H.. 1969. Mathematical simulation of the process of soil erosion by water. Trans. ASAE, 12(6): 754-762. Morris, E.M., 1979. The effect of the small-slope approximation and lower boundar~ conditions on solutions of the Saint-Venant equations. J. Hydrol., 40: 31-47. Morris, E.M. and Woolhiser, D.A., 1980. Unsteady one-dimensional flow over a plane: partial equilibrium and recession hydrographs. Wate Resour. Res., 16(2) 355-360. Moss, A.J., Walker, P.H. and Hutka, J., 1979. Raindrop-stimulated transportation in shallow water flows: an experimental study. Sediment Geol., 22: 165-184. Negev, N., 1967. A sediment model on a digital computer. Tech. Rep. Stanford Univ., Stanford, CA, Rept. No. 76. 109 pp. Palmer, R.S., 1964. The influence of a thin water layer on water drop impact forces, lnt. Assoc. Sci. Hydrot. Publ., 65: 141-148. Rowlinson, D.L. and Martin, G.L., 197t. Rational model describing slope erosion. J. lrrig. Drainage Div. ASCE, 97(1): 39-50. Singer, M.J. and Walker, P.H., 1983. Rainfall-runoff in soil erosion with simulated rainfall, overland flow and cover. Aust. J. Soil Res., 21: 109-t22. Smith, R.E., 1976. Simulating erosion dynamics with a deterministic distributed watershed model. Proc. Third Federal Interagency Sedimentation Conf., Water Resources Council, Washington, DC, Vol. 1, pp. 163-173.

MODELING THE EROSION PROCESS OVER STEEP SLOPES

305

Walker, P.H., Kinnell, P.I.A. and Green, P., 1978. Transport of a noncohesive sandy mixture in rainfall and runoff experiments. Soil Sci. Am. J., 42: 793-801. Wischmeier, W.H., 1960. Cropping-management factor evaluations for a universal soil loss equation. Soil Sci. Soc. Am. Proc., 24(4): 176-193.