USDA
*
8ig
fl) ~
-.
as e
o ££
E S3
££«>
•s"° S
030
c o «»
o£ ®
Q W DC
Unites States
Department of
Agriculture
Northeast Watershed
Center
University Park PA 1 6802
United States
Environmental Protection
Agency
Office of Environmental
Processes and Effects Research
Washington DC 20460
EPA-600/7-84-034
March 1984
Research and Development
Rill-lnterrill Erosion and
Deposition Model of
Stripmine Hydrology
Interagency
Energy/Environment
R&D Program
Report
L
-------
RILL-INTERRILL EROSION AND DEPOSITION MODEL
OF STRIPMINE HYDROLOGY
by
R. M. Khanbilvardi, A. S. Rogowski, and A. C. Miller
U.S. Department of Agriculture, ARS
Northeast Watershed Research Center
University Park, Pennsylvania 16802
EPA-IAG-D5-E763
Project Officer
Clinton W. Hall
Office of Energy, Minerals and Industry
Washington, D.C. 20460
Office of Research and Development
U.S. Environmental Protection Agency
Washington, D.C. 20460
U.S. Environmental Protection Agency
Region 5, Library (PL-12J)
17 West Jackson Boulevard, 12th Floor
Chicago. IL 60604-3590
-------
DISCLAIMER
This report has been reviewed by the Office of Energy, Minesoils
and Industry, U.S. Environmental Protection Agency, and approved for
publication. Approval does not signify that the contents necessarily
reflect the views and policies of the U.S. Environmental Protection
Agency, nor does mention'of trade names or commercial products
constitute endorsement or recommendation for use.
11
-------
FOREWORD
The Federal Water Pollution Control Act Amendments of 1972, in
part, stress the control of nonpoint source pollution. Sections 102
(C-l), 208 (b-2,F) and 304(e) authorize basin scale development of
water quality control plans and provide for area-wide waste treatment
management. The act and the amendments include, when warranted,
waters from agriculturally and silviculturally related nonpoint
sources, and requires the issuance of guidelines for both identifying
and evaluating the nature and extent of nonpoint source pollutants
and the methods to control these sources. Research at the Northeast
Watershed Research Center contributes to the aforementioned goals.
The major objectives of the Center are to:
• Study the major hydrologic and water-quality associated
problems of the Northeastern U.S. and
• Develop hydrologic and water quality simulation
capability useful for land-use planning. Initial
emphasis is on the hydrologically most severe
land uses of the Northeast.
Within the context of the Center's objectives, stripmining for
coal ranks as a major and hydrologically severe land use. In
addition, once the site is reclaimed and the conditions of the mining
permit are met, stripmined areas revert legally from point to nonpoint
sources. As a result, the hydrologic, physical, and chemical behavior
111
-------
of the reclaimed land needs to be understood directly and in terms
of control practices before the goals of Sections 102, 208 and 304
can be fully met.
Signed:
n
Harry B. Pionke
Director
Northeast Watershed
Research Center
IV
-------
ABSTRACT
An erosion-sediment yield model (labeled KEM) was developed
from the continuity consideration for sediment transport and from
equations describing rill and interrill erosion. This computerized
model is based on dividing the upland areas into a grid containing
rill and interrill zones and on the Universal Soil Loss Equation
(USLE). The USLE is used to estimate the sediment contribution from
the interrill areas. Prediction of soil loss from the interrill areas
is based on the premise that both raindrop impact and overland flow
energy can create soil erosion. The rill flow carries the interrill
erosion along with the rill scour. Rill transport capacity governs
the amount of removed soil from the site. If the flow transport
capacity is less than the available eroded soil, net erosion equals
the transport capacity and the excess sediment is deposited in the
flow paths. Otherwise, all eroded soil will move downslope and out
of the watershed.
The model was tested by simulating actual events on a small water-
shed in Central Pennsylvania for summer storms during 1981. Applying
the model to this stripmined and reclaimed area created a set of infor-
mation about the location and amount of watershed erosion and
deposition. The areal distribution of erosion and deposition was
compared with measured data. The model performed satisfactorily in
predicting soil loss from the site.
-------
CONCLUSIONS
An erosion sediment yield model (KEM) based upon partial area
hydrology has been developed and tested. The model was based on the
fundamental mechanics of erosion, rill transport, and deposition.
This computer model was formulated from the continuity consideration
for sediment transport and from equations describing rill and inter-
rill erosion. The model considers the watershed to be divided into
subwatersheds so that the error in selecting the parameters for the
Universal Soil Loss Equation can be reduced.
Prediction of soil loss was based on the premise that both
raindrop impact and overland flow energy can create soil erosion
anywhere on the watershed. Depending upon sediment transport capacity
of the flow, sediment may or may not move downslope.
The erosion process has been separated into rill and interrill
components. It was assumed that interrill erosion resulted from the
rainfall impact and the detached sediment was transported into micro-
channels (rills) by sheet flow occurring in the interrill areas.
Erosion in the rill was then considered to be a result of soil detach-
ment in the rill and the transport capacity of the rill flow.
The model is able to generate rill sources and rill patterns as
well as contributing interrill areas. Contributing interrill areas
VI
-------
are portions of interrill areas that contribute sediment to rill
microchannels. The model determines the eroded soil available for
transport at all points of a watershed and net amounts of soil eroded
for any area. The model can also predict sediment loads and timing
of stream inputs. The model predicts areal distribution of erosion
and deposition. This type of information could be used to plot iso-
erosion and deposition maps.
To test the model, it was applied to a stripmined and reclaimed
area in Central Pennsylvania. The model was executed for five storm
events which occurred in Summer 1981 (June 1 to September 1). The
rills and contributing interrill areas were generated for each storm
event. Then, the location and amount of watershed erosion and
deposition were determined. The outputs of all the storms were added
to each other. The end result was a map showing erosion or deposi-
tion in each subwatershed.
The model outputs were compared with soil loss measured using
erosion pins. The predicted results of erosion and deposition in
terms of being one or the other, were in good agreement with measured
data, the only discrepancy was in order of magnitude. Although the
simulated results overpredicted by 25% the measured values, it is
believed that the model is accurate enough to be useful on mined and
reclaimed areas.
Although this model requires many types of parameters, most of them
can be determined from published information or measured directly from
the maps or in the field. This model could increase the accuracy of
sediment yield predictions by allowing determination of zone and
vii
-------
subwatersheds contributions to the total sediment yield. The model
could be useful in evaluating the environmental impact of land use
practices. It also may serve as basis for reservoir and channel design
and land use planning.
Many problems have occurred while attempting to develop and apply
this model. It is believed that the following suggestions will be of
value for future studies involving erosion modeling from upland areas:
(1) Further studies are needed to enhance this model for unsteady
conditions.
(2) Many sediment transport equations are available to calculate
the sediment transport capacity. However, most of those
equations were derived for application in rivers where water
is deep. There is a need to evaluate the applicability of
these equations in a shallow flow condition.
(3) Model formulation did not include the ground water or
subsurface flow. In humid regions where subsurface or ground
water is dominant, there is a need to incorporate this effect
into the model.
(4) Further studies are needed to determine how land use activ-
ities such as road building, construction, and logging will
change parameters used in Philip's Equation and other
parameters related to sediment yield.
viii
-------
(5) Further studies are helpful to identify the difference
between the size distribution of the soil at source area
and the sediment transported in the watershed.
(6) This model similates the water and sediment yield for
individual storm. Further studies are suggested to incor-
porate the water and soil loss processes continuous
simulation.
(7) Further studies are needed to evaluate the influence of
rainfall energy on the transport capacity over land surfaces
and on turbulent mixing in the rill water.
(8) Further studies are needed to extend the application of
this model to a region where ground is frozen or covered
with ice.
ix
-------
CONTENTS
Page
Foreword iii
Abstract v
Conclusions vi
List of Figures xii
List of Tables xiv
1. Introduction 1
1.1 Need for Study 1
1.2 Erosion and its Problem 4
1.3 Objectives 5
1.4 Approach 5
2. Literature Review 7
2.1 Soil Erosion 7
2.1.1 Climate Condition 9
2.1.2 Topography 11
2.1.2.a Slope Steepness 11
2.1.2.b Slope Length 12
2.1.3 Soil 13
2.1.4 Vegetative Cover 14
2.1.5 Human Factors 15
2.2 Overland Flow 15
2.3 Soil Loss Determination 19
3. Model Formulation 41
3.1 Introduction 41
3.2 Physical Processes 42
3.2.1 Hydrology Components 48
3.2.2 Infiltration 48
3.3 Model Description and its Components. ... 49
3.3.1 Rill Sources 50
3.3.1.a Small Subwatersheds .... 51
3.3.1.b Large Subwatersheds .... 54
3.3.2 Rill Development 56
3.3.3 Contributing Areas-Partial
Area Concept 63
3.3.4 Governing Equation 64
3.3.4.a Interrill Erosion 66
3.3.4.b Contributions to a Rill . . 71
3.3.5 Governing Equations in Rills .... 71
3.3.5.a Rill Cross Section 72
3.3.5.b Rill Transport Capacity . . 76
3.3.5-c Rill Detachment Capacity. . 82
3.4 Data Requirements 84
-------
Page
4. Application of the Model and Results ........ 86
4.1 Introduction 86
4.2 The Study Area 86
4.3 Input Data ......... 89
4.3.1 Topography. ...... 89
4.3.2 Infiltration Parameters . 89
4.3.3 Precipitation 95
4.3.4 Soil Erodibility 99
4.3.5 Watershed Cover and Management
Factors 99
4.3.6 Watershed General Parameters. . . . 102
4.4 Subwatersheds Grouping .... 102
4.5 Results-and Discussion ..... 105
References .."............ 118
Appendices
A. Erosion Sediment Yield Model "KEM" .... 123
B. Sample Problem of Input and Output for KEM 165
XI
-------
LIST OF FIGURES
Figure Page
2.1. Simulation of soil erosion by water 28
3.1. Overall processes of erosion and sediment yield ... 44
3.2. Interrill areas and rill details (triangular ,_
approximation for rill cross section)
3.3. Factors influencing soil erosion in upland areas . . 46
3.4. Structure of the model 47
3.5. Variation of water viscosity with temperatures ... 53
3.6. Influence of runoff/infiltration ratio and soil
erodibility on rill source 55
3.7. A rill path on watershed 57
3.8. The possible directions of rill movement 58
3.9. Five possible directions for rill movement at
point A 60
3.10. Slope steepness in the five possible rill
directions 61
3.11. Shields diagram 79
3.12. Modified shields diagram 81
4.1. Flow diagram for the input-output files 87
4.2. The location of the study area 88
4.3. Location of area ABDC in the selected watershed ... 90
4.4. Topography of the selected watershed n-.
XI1
-------
Figure Page
4.5. Soil sorptivity (in/secl/2) of the selected 93
watershed, in three and two dimensions . . .....
4.6. A-values (in/sec) for the selected watershed, in
three and two dimensions .............. • ^
4.7. Distribution of site storms and SCS storm types ... 97
4.8. Soil er edibility of the selected watershed ..... 100
4.9. Vegetation density (lbs/100 sq ft) of the selected
watershed, in three and two dimensions ....... 101
4.10. Watershed ABDC and its grid system .... ..... 104
4.11. Rill patterns and contributing areas (shadow areas)
generated by model for storm 1 ..... . ..... 106
4.12. Eroded soil (Ibs) available for transport at the
end of storm 1 for section pqsr of area ABDC .... 107
4.13. Net amount (Ibs) of soil eroded (-) or deposited
(+) at the end of runoff period after
storm 1 for section pqsr of area ABDC ........ 108
4.14. Areal distribution of erosion and deposition on
each subwatershed .................. HQ
4.15. Details of erosion pin ............... H2
4.16. Erosion pins (•) distribution over the watershed
ABDC .................. ...... 113
4.17. Erosion and deposition details of section MM .... 114
4.18. Erosion and deposition details of section GG* .... 115
4.19. Comparison of the predicted and measured erosion
and deposition for watershed ABDC ..... .... 117
xiii
-------
LIST OF TABLES
Table Page
2.1. Classification of sediment yield models - 39
3.1. Values of a_ and b_ in Equation (3.6) for each
'SCS storm type 69
4.1. Precipitation information for Summer 1981 95
4.2. Summarized output for storm events of Summer 1981 . . 93
4.3. Assigned input information for the entire
watershed 103
xiv
-------
-------
SECTION 1
INTRODUCTION
1.1. Need for Study
Recently a sincere and dedicated concern has been shown by
the public for preserving our natural resources. Two of the most
important resources are soil and water. Continued sediment erosion
and pollution of our waterways have led to serious problems currently
being studied by engineers and the scientific community.
Among the existing problems facing the world, is the increase
in population. The area of land is fixed, but population continues
to grow. This increased population demands more homes, highways, and
industries which in turn reduce the amount of land available for food
production. Consequently, soil erosion and sedimentation are of par-
ticular importance in the evaluation of water and soil resources.
Soil erosion could result in significant deterioration of the land's
long term productive capability.
The erosion process and the resulting sediment are the most impor-
tant problems, because they not only affect the water quality in the
stream but also may restrict future land use. Also, most pesticides
and nutrients attach themselves to the clay particles, and thus are
-------
carried away when erosion occurs. These eroded soils may enter the
streams and pollute the natural water resources.
Soil erosion is the primary source of sedimentation, because
eroded soil particles are transported to rivers and deposited. These
sediments may result in the premature filling of lakes and reservoirs
and reduce the capacity of rivers to carry flood flows. In irrigation
systems, sediment may greatly reduce the carrying capacity of con-
structed artificial channels.
How can this situation be resolved to satisfy the needs and desires
of an ever more demanding society with a limited supply of land resources?
The problems associated with land use can be solved using the knowledge
of land characteristics from which various land use opportunities can
be interpreted. Mathematical modeling and system analysis techniques
offer a unique opportunity to analyze various control technologies
and provide methods to select those which are most efficient.
In upland areas, especially mountainous watersheds where soil
erosion is likely to be a more serious problem, erosion models are
needed to develop procedures for reservoir design and to assess the
impact of changes including irrigation, urbanization, mining, and
rural constructions. These models are highly beneficial as a design
tool in developing plans to control soil loss and erosion damage on
upland areas. Present models vary considerably with respect to
generality and applicability. Most of them depend on a set of
historical data to calibrate the unknown parameters. Thus, these
models must be used with caution if significant changes in the basin
-------
land use have occurred or the events differ greatly from those in the
calibration data.
Many methods are available to estimate soil erosion and sediment
yield and one may ask, why should another model be developed? The
answer to this has three parts:
(1) Although several equations have been developed to predict
soil erosion, it has been pointed out by many investigators
that the applicability of these soil loss equations may be
exceeded when they are applied to a watershed because they
were developed from research data on small plots for a
particular land use and there is no assurance that a tech-
nique good on a small scale in one land resource area will
apply to another.
(2) There is a new emphasis on evaluating the impact of changes
in water quality and quantity with changes in land use. Thus,
there is a need to develop an erosion model from basic prin-
ciples to any conditions. Erosion is a complex phenomenon
and should not be restricted to estimating only the amount
of eroded soil from a watershed.
(3) Most models developed to predict soil erosion from a water-
shed usually consider the entire drainage area as contributing
to erosion. Therefore, they reflect gradient and soil type
along the entire slope length. This implies that the entire
-------
watershed produces runoff. This type of approach cannot
be correct because realistically only portions of a
watershed can produce runoff.
1.2. Erosion and Its Problem
Erosion is the removal of soil particles from the land surface.
This process is the result of energy developed by the water as it
falls toward the earth and as it flows over the surface of the land.
Problems associated with soil erosion are many. Erosion on crop-
land degrades the productivity of the soil resource necessary for crop
and food production. As a result of soil erosion, plant nutrients and
fine particles are selectively removed causing poor soil tilth and
increasing runoff because of poor infiltration.
Eroded soil not only reduces soil productivity, but also is a
pollutant itself. It may carry soil-adsorbed chemicals, which are
themselves pollutants. Sediment pollutes by muddying the water,
clogging fish gills, covering nests and spawning grounds, and
increasing the dissolved oxygen demand. In those communities in
which surface water is a source of water supply, sediment removal
from water for public use can be costly, because, large investments in
the construction and operation of treatment facilities are generally
required. Similarly, sediment deposition in stream channels, rivers,
and lakes usually require costly removal.
-------
1.3. Objectives
The major objective of this study is to develop a mathematical
model having the following characteristics:
(1) Predict soil erosion from upland areas using the partial
area concept.
(2) Predict probable location of rills and the total amount of
sediment transported to a stream or sediment pond following
a specific design storm.
(3) Use readily available data that will be user oriented.
The model is to be based on basic principles of erosion mechanics
and experimental information available so that it can be used for
simulating soil loss and sediment yield due to changes in land use.
1.4. Approach
In order to accomplish the objectives of this study, a model that
describes the spatial changes in land use, geometry, and soil character-
istics is needed. The idea that both raindrop impact and overland flow
can create soil erosion and the partial area concept form the basis
for this model.
Rill development is simulated on a computer. After generating
the rill patterns, the contributing interrill areas are delineated.
-------
The eroded soil from rills and contributing interrill areas is calcu-
lated and routed along the rill patterns. When the rill transport
capacity exceeds the available detached soil, deposition does not
occur. If the transport capacity is less than the total soil avail-
able to be transported, the amount of deposition is assumed to be
the difference between the detached load and the rill transport
capacity.
-------
SECTION 2
LITERATURE REVIEW
Erosion prediction from upland areas is highly beneficial as
a design tool in developing plans to control sediment loss and
erosion damage. The relationship between rainfall, runoff (overland),
soil erosion, and sediment yield has been studied by many specialists
for many years. This review of literature assembles information con-
cerning runoff (overland flow), soil erosion, and sediment yields.
For clarity, this chapter is divided into three sections: soil
erosion, overland flow, and soil loss determination.
2.1. Soil Erosion
Erosion is defined as the abrading of the land. It is the removal
of soil materials by erosive agents and refers to two phases of the
process of detaching and transporting. Erosion could be caused by
wind or water, each has a different process and characteristics.
Since wind erosion is not part of this study, it will not be dis-
cussed here.
Brant et al. (1972) have listed the most important sources of
erosion and sediment yield:
-------
(1) Natural erosion occurs from the weathering of soils,
rocks, and uncultivated land, and geological abrasion.
(2) Agricultural erosion is a major sediment source due to
the large area involved and the land distrubing effects of
cultivation.
(3) Urban and rural erosion originates mainly from exposed
bare soil in areas under construction or from lands disturbed
by mining.
(4) Highway erosion is associated with the stripping of large
areas of their vegetative protection during road construc-
tion (i.e., landslides).
(5) Stream bank, channel, and shoreline erosion from concen-
trated water flows and wave action in channels and
floodplains.
Raindrops hitting the earth's surface initiate the process of soil
erosion by water. Rainfall detaches soil particles by drop impact and
transports them by splash. Detachment capacity of rainfall depends on
the diameter of the raindrop, distribution, velocity, and total mass
or kinetic energy at impact.
Raindrop impact and flowing water are the erosive agents involved
in the erosion process. However, the forces of gravity and cohesion
are working against the erosive agents. Therefore, erosion will take
-------
place whenever eroding forces exceed the forces of resistance. When
there is no more surface storage capacity and the rainfall intensity
exceeds the infiltration rate of the soil, there would be overland
flow which moves the eroded soil particles downslope.
The basic factors affecting soil erosion from watersheds are
climate, topography, soil, vegetative cover, and human activities. These
factors are briefly described here.
2.1.1. Climate Condition
Rainfall is obviously the most important factor among the other
major climatic factors of wind, temperature, and snow. Rainfall
erosivity depends on the rainfall intensity, duration, frequency, and
shear force exerted on the ground surface. Wischmeier and Smith
(1958) found that the rainstorm parameter most highly correlated
with soil loss from fallow ground was a product term: kinetic energy
of the storm times maximum 30-minute rainfall intensity. This
product is called the "rainfall erosion index" and explains 72 to 97
percent of the variation in individual storm erosion from tilled con-
tinuous fallow ground on each of six widely scattered soils (Smith
and Wischmeier, 1962).
Raindrop splash and overland flow are major agents in transporting
the eroded material. Overland flow is surface runoff which travels
over the ground surface to a channel. Overland flow tends to be
channelized and make rills and gullies, whereas splash erosion tends
to remove soil particles from the surface as a uniform sheet layer.
-------
Meyer et al. (1975a) studied effect of flow rate and canopy on
rill erosion and found that soil loss due to the rill flow rate can
be fitted by the equation:
DR(Q-QC) (2.1)
where;
E = total rill and interrill erosion per unit of rill length
(wt/time/length),
E = erosion rate on interrill areas per unit of rill length
(wt/time/length),
D = detachment rate in rills per unit of rill length per unit
of excess discharge (wt/time/length/discharge),
Q = rill discharge rate (discharge = vol/time),
Q = critical discharge below which rill erosion is negligible
(vol/time), and
E_ = ET when Q £ Q .
10
-------
2.1.2. Topography
The effect of topography can be explained in terms of slope steep-
ness and slope length.
2.1.2.a. Slope Steepness
Slope steepness is the main topographic feature that is highly
correlated with soil erosion. Zingg (1940) studied the effect of
degree of slope on soil loss and found that on slopes of less than
10 percent, erosion approximately doubled as slope, expressed as
percentage, increased two-fold.
Foster and Martin (1969) found that on slopes above 20 percent,
the erosion rate tended to level off for some conditions. Thus, they
found that depending on bulk density, erosion on short slopes from
35 to 100% reached a maximum and then decreased as the slope steep-
ened.
The influence of slope steepness on the transport capacity of
the flow is also very important. Most studies show that transport
capacity changes as some power greater than two of the energy grade
line which approximately is equal to the land slope. Thus, the
increase in land slope would increase transport capacity rapidly.
11
-------
2.1.2.b. Slope Length
Zingg (1940) also explained the influence of slope length on
the average erosion per unit area for a given slope. His expression
was in the form of:
A = Ln (2.2)
where;
A = Average erosion per unit area,
L = slope length,
n = a coefficient.
Wischmeier and Smith (1965) used a value of n near n = 0.5 for most
circumstances. However, Foster and Meyer (1972a) indicated that the
value of n depends on the relative susceptibility of different soils
to rilling and the resulting ratio of rill erosion to interrill
erosion. They indicated that where soil loss is primarily from rills,
n will approach one, but if the interrill erosion is dominant it will
approach zero. Young and Mutchler (1969) also showed that n increases
with increasing slope length because rill erosion increases faster
than interrill erosion.
12
-------
Smith et al^ (1945) used five plots with lengths up to 270 ft.
Erosion loss was measured at downslope intervals of 10 ft for five
years and the best fitted equation was reported as:
Y = 0.016 L°'57 (2.3)
where;
Y = average depth of soil loss (ft),
L = slope length (ft).
2.1.3. Soil
The major factors affecting soil erosion are texture, structure,
permeability, compactness, and infiltration capacity of the soil
profile. Soil texture determines the permeability and erodibility
of soil. Erodibility, detachability, and transportability of soil
directly influence the rate and amount of soil erosion. However,
under the same hydraulic, climatic, and vegetative cover different
types of soil might have different erodibility and soil losses.
Middleton's (1930) effort to determine the erodibility of soils
led to the suggestion that the "dispersion ratio" and "erosion ratio"
as erosion indices relate erosion to the physical properties of soil.
The "dispersion ratio" was obtained by dividing the amount of silt
and clay in a sediment sample by the total quantity of silt plus clay
13
-------
present in the soil, and the "erosion ratio" was obtained by dividing
the dispersion ratio by che colloid moisture equivalent ratio (which
is, the percentage of water retained by a sample of soil one centimeter
deep that has been saturated with water and drained under a centrifugal
force 1000 times gravity for 30 minutes). He concluded that the greater
its dispersion and erosion ratios, the greater the erosion of a soil.
Based on these criteria soils have been divided according to erod-
ibility and nonerodibility as follows:
(1) If "dispersion ratio" is greater than 10 and "erosion
ratio" is greater than 15, soil is erodible.
(2) If "dispersion ratio" is less than 10 and "erosion
ratio" is less than 15, the soil is nonerodible.
2.1.4. Vegetative Cover
Plant cover is one of the best protections against soil loss.
Plant cover affects both the infiltration rate and the susceptibility
of soil to erosion. It causes the absorption of raindrop impact and
the reduction of overland flow velocity and tractive force by
increasing the hydraulic roughness and decreasing the effective
slope (Baver, 1965).
Mulches and vegetation increase the hydraulic roughness, and
decrease the effective slope steepness, therefore, they reduce the
runoff velocity and erosion. The effect of crops and their
14
-------
management cannot be evaluated independently. A crop can be grown
continuously or in rotations. The sequences of crops within a
system can be varied, and therefore, different combinations of these
variables might have different effects on soil loss.
2.1.5. Human Factors
Most of the erosion control and conservation practices may be
considered as human factors. Human beings disturb the soil, manipu-
late the vegetation and change the natural sequence of evolution. The
legislative and administrative erosion control measures are the most
effective means of preventing soil erosion or reducing it.
Beasly (1974) explained major practices to prevent soil erosion.
Contour tillage, contour strip-cropping, and terracing with contour
farming are among them. He indicated that terracing with contour
farming is the most effective because with terracing, the sediment
deposits in the terrace channel and may equal up to 90 percent of all
the soil moved to the channel. Other practices such as diversion
waterways, ponds, reservoirs, check dams, and gully control struc-
tures also may be considered as technical and engineering measures.
2.2. Overland Flow
Overland flow formulation and solution has been of great interest
to engineering communities and the reported methods are divided into
15'
-------
three groups; regression models, frequency analysis, and physical
processes models (Linsley, 1971).
The use of hydraulic procedures for predicting overland flow is
associated with some problems. Overland flow depends on rainfall
supply which can be depleted by infiltration. Since both of these
elements vary with time and location, the overland flow could be both
unsteady and spatially varied. Depending on the rate of flow and
nature of land surface the flow could be laminar, turbulent, or
both. The impact of raindrop and formation of roll waves provide an
additional complication in overland flow (Robertson et al., 1964).
There are many factors affecting the volume of water obtained
from a rainfall or storm. These factors could be described as follows:
1. Rainfall characteristics
a. rainfall intensity
b. rainfall duration
c. time distribution
d. spatial distribution
2. Watershed characteristics
a. watershed size
b. watershed shape
c. slope of watershed
d. vegetative cover and its density
3. Soil characteristics
a. shape
16
-------
b. size distribution
c. unit weight
d. porosity
4. Climatic factor
a. temperature
b. wind direction
c. wind speed
d. antecedent moisture condition
5. Man-made factors
a. land use
b. construction.
As mentioned earlier, three types of overland flow formulation
exist: regression models, frequency analysis, and physical pro-
cesses. Regression models primarily attempt to establish a
mathematical expression to relate rainfall to runoff. Frequency
analysis uses statistical characteristics of the recorded rainfall
or runoff to generate or synthesize nonrecorded events. Physical
process models deal with the concept of water balance and divide
the processes of rainfall and runoff into components (or parameters)
Overland flow and runoff transport the detached soil materials.
Only when precipitation rate exceeds infiltration and all surface
depression storage is exhausted will runoff occur.
17
-------
Infiltration is defined as the rate at which water percolates
into the soil. The equations presented in the literature describe
the infiltration rate either by using the empirical concept (Horton,
1939; Holtan, 1961) or by equations based on the physical concept of
water entry into the soil (Green and Arapt, 1911; Philip, 1957).
Horton's equation describes the exponential decrease of the
infiltration rate. His equation is:
f = f + (f - f )e~kt (2.4)
c o c
where;
f = infiltration rate (in/hr),
f = the infiltration rate assumed similar to the saturation
c
permeability (in/hr),
f = the initial infiltration rate (in/hr),
o
t = time (hrs), and
k = a constant.
However, the difficulties with determining the parameters f and k
cause this equation to be of lesser importance to the rainfall-runoff
modeling of watershed.
18
-------
Philip's (1957) infiltration equation is based on soil physics
of water movement in porous media. The equation is:
1/2
I - S t + A t (2.5)
where;
I = cumulative infiltration (in),
S = sorptivity which depends upon moisture content and diffusivity
1 /?
of the soil (in/sec ' ),
A = a coefficient depending on conductivity at the wetting
front (in/sec), and
t = time (sec).
Infiltration can be partially controlled by engineering and agricul-
tural practices, such as tillage, raking of the surface, and compaction.
2.3. Soil Loss Determination
The soil erosion and consequently soil losses are among the
important problems that engineering communities have been faced with
during the last decades. The efforts of many researchers in predicting
soil losses from farm land and agricultural watersheds have been
19
-------
reported since the 1930's. Zingg (1940) expressed soil loss as a
function of slope length and steepness. The equation is:
= 0.026(Se)3-37(Le)1-60 (2.6)
where;
X = total soil loss (Ibs) ,
S = land slope in percentage, and
L horizontal length of land slope (ft)
Rainfall impact on the soil surface has long been considered
the initial phase of the water erosion process. Ellison (1944) pub-
lished an equation expressing soil erosion by splash as a function
of raindrop size, velocity, and intensity. His equation is of the
form:
(2.7)
where;
E = soil splashed in pounds during a 30-minute period,
S
20
-------
R = a constant,
V = raindrop velocity (ft/sec),
d = raindrop diameter (in), and
I = rainfall intensity (in/hr).
Musgrave (1947) found that the rate of sheet erosion was related
to a number of factors and that certain relationships existed between
these factors. He presented the following equation:
T = K (S/10)1'35(L/72.5)°-35(P_n/1.25)1'75 (2.8)
s g 30
where;
T = the probable soil loss (tons/acre/year),
S
K = a soil factor depending on the erodibility of soil and cover,
g
S = slope steepness (in percent),
L = slope length (in feet), and
P = the maximum 30-minute rainfall expected in a two year
period.
Sheet erosion is not easily observed in the field because the
irregularities on the soil surface cause minor rills to form. When
21
-------
these so-called micro channels form, they can be deepened and end
as gullies. The significant erosion by the surface water would be
within these channels. In scour erosion, rolling, lifting, and
abrading of soil particles are the main process of soil detachment.
The studies of many researchers led to the development of
Universal Soil Loss Equation (USLE). This equation (Wischmeier and
Smith, 1965) is currently the most comprehensive and popular regres-
sion model to estimate soil loss. The Universal Soil Loss Equation
is:
E=RKLSCP (2.9)
where;
E = average annual soil loss (tons/acre),
R = the rainfall factor (the number of erosion index units
in a normal year's rains),
K = soil erodibility factor (the erosion rate per unit of
erosion index for a specific soil in cultivated continuous fallow
on a 9 percent slope and 72.6 ft long),
L = slope length factor (the ratio of soil loss from a field
for a given slope length to that from a 72.6 ft slope length on the
same type and gradient),
22
-------
S = slope-gradient factor (the ratio of soil loss from the
field gradient to that from a 9 percent slope),
C = cropping-management factor (the ratio of soil loss from a
field with specific cropping and management to that from the fallow
condition on which the factor K is evaluated), and
P = erosion-control factor (the ratio of soil loss with con-
touring, strip-cropping, or terracing to that with straight; row
farming, up- and down-slope).
The Universal Soil Loss Equation was developed from more than 10,000
plot years of soil loss data. Usually, it is applied to farmlands
because of massive amounts of data available for such areas. The
USLE is not a complete sediment yield equation and does not include
transport and deposition phenomenon of sediment yield. It is limited
to annual soil loss for the entire watershed. However, the USLE has
been used by Soil Conservation Service for many years.
Williams (1974) showed that erosion could be related to the
transportation process more than the detachment process, and in order
to characterize the transportation process, he used the volume and rate
of runoff instead of rainfall factor. He proposed the Modified Uni-
versal Soil Loss Equation (MUSLE) in the form of:
E = 95(Q • q )°'56(K LS C P) (2.10)
23
-------
where;
E = sediment yield from an individual storm (tons),
Q = storm runoff volume (acre-feet),
q = peak runoff rate (cubic feet/sec),
K, L, 5, C, and P are the same factors previously defined in
the USLE.
In MUSLE a runoff factor composed of the runoff volume and peak
runoff rate has been used instead of rainfall factor (R) in USLE.
Williams (1974) found a close agreement between the predicted and
measured data while using the MUSLE for computing the sediment yield
for several different floods..
Foster et al. (1973) proposed a modified form of the USLE to
reflect both rainfall and runoff contribution. This equation is:
E = E K C P L S (2.11)
n
where;
E = 0.50 + 15 Q q 1/3,
Q = runoff volume (in),
q = runoff rate (in/hr),
P
24
-------
E, K. C, P, S, and L are the same as in the USLE.
One of the disadvantages of the regression type equations is
that they are limited to their own particular data base and therefore
restricted in application to those conditions and they lack a certain
flexibility and could not interpret the physical process involving the
soil erosion. Therefore, a need for improving the physical process
modeling was apparent since all the above mentioned soil loss equations
are typically based on regression analysis.
Meyer and Wischmeier (1969) subdivided the erosion process into
four sub-processes. These are soil detachment by rainfall, transport
by rainfall, detachment by runoff, and transport by runoff. They used
the following mathematical models to evaluate each subprocess.
(1) Soil detachment by rainfall, D :
(2) Transport by rainfall, T :
(2.12)
(2.13)
(3) Detachment by runoff, D :
F
2/1
25
-------
(4) Soil transport by runoff, T :
TF
where;
A. * the area of each increment (sq ft),
S = rainfall detachment factor,
ST = rainfall transport factor,
SDF = runoff detachment factor,
S = runoff transport factor,
IF
r = rainfall intensity (in/hr),
S = slope steepness, and
q = flow rate (cfs).
In fact, this was the first time that subprocesses of soil
erosion were tied together in a mathematical form. The interaction
between these four subprocesses determines the total erosion from
each part of the watershed. The four submodels could be combined
to route the soil downslope. The sediment load carried from each
increment is the lesser of the sediment lo.ads from the previous incre-
ment plus the detachment in that increment or the transport capacity
from that increment. Net erosion for each increment is the difference
26
-------
between the sediment loads entering and leaving it. This conceptual
model can be expressed in diagram form as shown in Figure 2.1.
Foster and Meyer (1972a) published the closed-form soil erosion
equation. They derived this equation from basic hydraulic and sediment
transport theory. According to their approach, the basic governing
equation of the erosion process is the continuity equation for sedi-
ment transport which can be shown by the following formula:
3G
where;
G_ = sediment load in flow (wt/time/unit of cross section width),
F
x = distance along the flow surface,
D.., = flow detachment (wt/unit area/time), and
= rainfall detachment rate (wt/unit area/time)
The sign convention is: when D > 0 there is detachment, and
when D_ < 0 there is deposition. An interrelationship between detach-
F
ment by runoff and sediment load carried by runoff was given by:
D G
JU,--!
27
-------
ro
oo
Detachment by
rainfall
*
i
'
Detachment by
runoff
F
_j
Detachment of
increment (DET)
Residual soil
I
Carried downslope
IF DET > TRANS
»
Deposition (DET-TRANS)
with residual, TRANS being
carried to next subarea
t
Compare t
with tran
f
Transport Transport
capacity of capacity of
rainfall Tn runoff T
K r
t *
otal detachment
sport capacity
*
Residual soil carried
t_
t
Transport capacity
of increment (TRANS)
IF TRANS > DET
t
No deposition occurring
and residual, DET being
carried to the next
subarea
f
downslope
Figure 2.1. Simulation of soil erosion by water.
-------
where;
D = the detachment capacity of the flow,
T = the transport capacity of the flow,
D and G were defined in the previous equation.
F F
The terms D , T , and R^ are independent variables defined by rain-
fall and runoff characteristics and by soil properties. Since rainfall
and runoff both provide the energy required to detach and transport
soil particles, the hydrologic process should be considered in the
erosion simulation. Therefore, the important characteristics of
rainfall and runoff were assumed either to be known or to be avail-
able from stochastic generation. Finally, the closed-form soil
erosion equations were expressed as (Foster and Meyer, 1972a) :
(2.18)
and
" (2-19)
where;
G /T = relative sediment load (dimensionless),
29
-------
x/L = relative distance (dimensionless) ,
a = L D /T = a flow erosion parameter (dimensionless),
o co co
6 = L RDT/T = a rainfall erosion parameter (dimensionless) ,
A = sediment yield (tons/ft/hr) ,
T = transport rate of flow at the end of a uniform slope
(wt/unit width/time),
L = length of slope (ft), and
D = detachment capacity of flow at the end of a uniform slope
(wt/unit area/time).
Foster and Wischmeier (1974) indicated that most equations were
derived from data which was obtained from a uniform slope and they do
not reflect the real phenomena of sediment movement. The effect of
slope irregularity on sediment load is not accurately simulated by
the overall average steepness. In relation to this fact and to
detachment and transport capacity of soil erosion, Onstad and Foster
(1975) have shown that the detachment capacity for any segment of
a complex slope composed of uniform slope segments can be represented
as follows:
W.(KCPS) ,
DCOJ = -(x- (2-20)
30
-------
where;
D . = detachment capacity of flow at the end of a uniform
coj f J
slope segment j (wt/unit area/time),
1/3
W. = an energy term for segment j (SI unit) =0.50R+15Q.q ,
K = soil erodibility factor of USLE (tons/acre),
C = crop management factor of USLE (dimensionless),
P = soil conservation practice factor of USLE (dimensionless),
S = slope steepness factor of USLE (dimensionless),
x. = distance from upper end of slope to lower end of segment j
(feet), "
x. = distance from upper end of slope to upper end of segment j
(feet),
R = storm rainfall factor of the USLE (Si unit),
S t
Q = storm runoff volume (inches), and
q = storm peak runoff rate (in/hr).
Equation (2.20) shows that each segment may have a unique set of
input parameters. Applying this concept to a watershed with complex
slope having n segments, the following equation would result:
31
-------
n
G = E D . (2.21)
• i
J-l
where;
G = the sediment yield of the slope.
They indicated that soil transport capacity is not the limiting factor,
However, they used the USLE to analyze the transport capacity. The
transport capacity (T ) at any point x in Ib/ft of width is:
T W K SCP x
c 185.58
where;
W = average rainfall energy term (SI unit),
K = an average USLE soil erodibility weighted on the basis of
contribution of each segment to sediment load.
After calculating the sediment yields and transport capacity at the
bottom of each slope segment of a complex profile, if transport
capacity exceeded the detached load of the segment plus any upstream
contribution, the sediment yield was the sum of the detached load and
the upstream contribution. Deposition did not occur on the segment.
32
-------
However, if the transport capacity was less than the total soil
available to be transported the sediment yield equaled the transport
capacity and the rest of the soil was deposited on the segment. Cal-
culations are carried out in this manner downslope to the channel
was reached.
Kuh et a. 1. (1976) used the aforementioned studies as input into
a two dimensional model which can predict both the total amount of
erosion from a watershed and areal distribution of erosion and sedi-
ment deposition. The end result was a map showing erosion or deposition
in each grid on the watershed. Subsequently, iso-erodent lines were
developed showing areas of erosion and deposition for one storm. How-
ever, they indicated that more research was needed for evaluating the
flexibility of their model for use under different conditions.
Bennet (1974) presented a conceptual model to estimate sediment
yield. The model divided the watershed into an upland and a lowland
channel area. The concept of water continuity, momentum, and sediment
continuity were used to formulate the erosion and sediment yield
processes of both areas. The transporting capacity in the overland
flow area and the channel systems, the variation of the channel bed,
and the meandering of the channel system are considered to affect the
rate of sediment transport and deposition.
Negev's (1967) model simulates generation and transport of soil
by raindrop impact and overland flow. He presented the following
formula for determining the production of fine soil particles by
raindrop splash:
33
-------
A(t) = (1 - FVC) x K x (P )AE (2.23)
n t
where;
A(t) = weight of soil particles produced during time interval t,
FVC = fraction of vegetative cover as a function of the rela-
tive value during the growing season,
K = the coefficient of soil properties (depending on the soil
erodibility),
P = precipitation during the time interval t, and
AE = an exponent.
Fine particles which are produced are available for transport by over-
land flow. Depending on the transport capacity of flow these soil
particles would be either removed or deposited. He modeled this
process as follows:
AN
SPT(t) = COT x RDS(t - 1) x OF(t) (2.24)
where;
SPT(t) = soil particles transported during time interval t
(weight),
34
-------
COT = coefficient of transport (depending on the transport
ability of overland flow),
RDS(t - 1) = reservoir of deposited soil particles existing
at the beginning of the time interval t,
OF(t) = the overland flow occurring during time interval t,
AN = an exponent.
David and Beer (1975) divided the erosion loss into rill and
interrill contributions. They believed that in the interrill zone
the soil erosion is mainly by the effect of raindrop splash. They
proposed the following equation to determine interrill soil erosion.
Ed=° (Scf)(Lsf)(Ir)" (2'25)
where;
E, = erosion caused by raindrop splash (weight)
S . = soil cover factor,
L f = land slope factor,
I = rainfall intensity,
d = overland flow depth,
2 and w = coefficients.
35
-------
They indicated that in the rill zone, soil detachment is mainly
by overland flow which takes place along the rill. The following
equation describes this part of erosion:
E = C'dU (2.26)
where;
E = amount of rill flow scour (weight),
C1 = a constant representing the soil characteristics and the
overland flow surface slope,
d = the overland flow depth, and
u = a constant.
They indicated that the detached material that is not transported
may be redeposited. These redeposited soil particles are left loosely
on the ground as detachment storage until the next overland flow
occurs. The following expression then describes the rate at which
the total detachment storage decreases (David and Beer, 1975):
—Rf
D = D e (2.27)
t o
36
-------
where;
D = total detachment storage at the end of the time interval,
D = total detachment storage at the beginning of the time
interval,
R = a /a , where a is a soil related factor and a a climatic
s c s c
factor, and
t = the time interval.
Simons e_t al. (1975) developed a model to simulate the processes
of erosion and sediment yield from a forested watershed. The mechanics
of water and sediment routing, the effect of particle size on erosion
rate and transporting capacity, and the processes of degradation and
aggradation in the channel system were considered in their model.
The governing equations were the water continuity, the momentum, and
the sediment continuity equations. To solve the flow on land surface
and in channels, kinematic wave approximation was used. Meyer-Peter,
Muller's bed equation (1935) and Einstein's (1950) suspension procedure
were used in computing the sediment transport capacity. Though this
model appears to be very accurate, it requires large amounts of data
and computer time to run. So, Simons e_t_ al_. (1977) developed a single
plane model which they call "physical process model." In that model
instead of routing flow over time and space using finite difference
formulations they average the physical processes over both time and
space to obtain a simple approximation of the complex model.
37
-------
The models or methods which have been discussed in this chapter
are useful for predicting the soil erosion and sediment yield from a
watershed. However, the advantages, disadvantages, limitations, and
the applicability of them were not mentioned. Without the same sets
of information to compute sediment yield for all of them, it is almost
impossible to make comparisons. However, Shiao (1978) tried to classify
them on the basis of their approach and theory. Table 2.1 shows this
classification scheme. According to this table, the sediment yield
models can be divided into three groups:
Group A - The Universal Soil Loss Equation (USLE) and its modifi-
cations. The models in this group were derived using
a regression technique. During the years, some of
these methods have been modified and their limitations
developed to a certain degree. These techniques have
been developed based on available data for specified
areas, thus, the use of them for other areas creates
some problems.
Group B - Those models use the concept of water continuity and
balance the ratio of detachment and transport. The
models are usually an improvement of the Group A type
equations, since physical processes of soil erosion
and sediment yield have been used in their formulations.
Most of these models were developed on small areas or
plots, and therefore their application is restricted.
38
-------
Table 2.1. Classification of sediment yield models.
Group
Models
Group A
USLE (Wischmeier and Smith, 1965),
Musgrave (1947), Williams (1975), Onstad and
Foster (1975)
Group B
Meyer and Wischmeier's (1969) conceptual model,
Foster and Meyer (1972), David and Beer (1975),
Negev (1967)
Group C
Bennet (1974), Simon et al. (1975, 1977)
39
-------
Group C - Those models with the complete processes of erosion
and sediment yield. The models in this group are
based on the physical and hydraulic processes of erosion
and are usually more complex than other groups. The
model described by Bennet (1974) is although only
conceptual, is a good example as well as Simons et a_l.
(1975) model which was tested with limited data.
40
-------
SECTION 3
MODEL FORMULATION
3.1. Introduction
The processes of erosion and sediment yield are a complex
physical phenomena. A mathematical model that attempts to predict
soil erosion caused by rainfall is generally based on the fundamental
factors involved in the erosion process.
In order to solve the complexity of the erosion process, the
overall process must first be subdivided into several subprocesses
that can be studied individually. In the development of this model,
not all the physical processes of water erosion and sediment yield are
included. However, an attempt has been made to develop a model that
is logically sound, scientifically correct, practical, and uses
easily obtainable data. This model has been called KEM and the
Fortran Program listing along with an example is given in Appendix
A and B. Some of the existing models are so complex that use of
them is difficult and at times impossible because of the lack of
available data.
Those processes which were considered in this model are presented
in the following sections. However, it should be noted that the
41
-------
drainage area in this model is represented by 10,000 square sub-
watersheds with homogeneous parameters, or by a square (100 x 100)
matrix in computer format.
3.2. Physical Processes
Soil erosion is the result of two principal physical phenomena:
the detachment of soil particles from the soil mass and the transport
of these particles. The erosion process is possible if erosive agents
are available. For water erosion, the area of interest of this
research, the erosive agents are rainfall and runoff.
Raindrop impact is the primary source of energy for detaching
soil from any land area. The impact breaks the soil aggregates into
particles. Rainfall reaches the ground and percolates into the soil.
If the subsoil is not saturated, the initial infiltration rate may be
higher than the rainfall intensity, and the water penetrates into the
soil. When the soil becomes saturated or the infiltration rate is
less than the rainfall intensity, then water accumulates on the ground
surface and begins to move by gravity as overland flow over the ground
surface. If rainfall continues, the depth of overland flow increases
and exerts a shear force, large enough to move the already loosened
soil particles or erode soil from the ground surface.
Erosion begins at a place where the soil is most susceptible to
erosion. Generally, if the binding forces of the soil particles are
small, then erosion is more probable.
42
-------
In general, runoff on erodible soil concentrates in many small
channels. These microchannels are called rills, and the areas between
the rills are then defined as interrill areas. The continuation of
erosion processes, increases the rill dimensions and ultimately a gully
will form. If the rainfall stops, or if it decreases, .the transport
capacity of the water is reduced significantly. Consequently, soil
particles being transported by the runoff will be deposited as the
runoff jnoves toward a stream or outlet point. Figure 3.1 shows the
overall processes of erosion and sediment yield, and Figure 3.2 shows
a typical rill cross section and its details.
Soil erosion from a watershed depends on many factors such as soil
moisture, natural topography, vegetative cover, and the forces created
by rainfall and the resulting runoff. Soils with a strong binding
force such as clay are less likely to erode and soils with large
particle sizes are more stable and therefore less likely to move.
Slope also is a major factor, the flatter the slope the less sus-
ceptible the particles are to transport. Ground surface with dense
vegetative cover not only reduce the detachment capacity of raindrop
impact, but may retard the flow rates which consequently reduces the
shear forces exerted by the flow. Figure 3.3 summarizes the factors
influencing the soil erosion and sediment yield processes.
The structure of the model is presented in Figure 3.4. This
model is formulated on the assumption that erosion process is divided
into rill and interrill erosion according to sources of eroded sedi-
ment.
43
-------
Precipitation
V
I
Rainfall
detachment
Sheet flow
Rill flow
K(s~ediment transport)
is
Infiltration .
Sediment
transport
in channel
_sr / ^~
Figure 3.1. Overall processes of erosion and sediment yield.
44
-------
Rill flows
Interrill
width
Rill
flow
width
OVERLAND FLOW IN RILLS
Figure 3.2. Interrill areas and rill details (triangular approxi-
mation for rill cross section).
45
-------
Soil erosion
ON
Climate
Soil
characteristics
Topography
Soil cover
condition
Management
activities
Rainfall intensity
and storm duration
Sorptivity and A-
value determine
Infiltration
Properties of soil
(i . I--. , specific
gravity)
Degree and length
of slope
Vegetation and
ground cover - C
factor in USLE
Ground surface
activities - p
factor in USLE
Figure 3.3. Factors influencing soil erosion in upland areas.
-------
Determine the infiltration
Determine the rainfall excess
Determine the rill sources
Determine the rill patterns
Determine the contributing areas
(effective interrill areas)
Determine the eroded soil from
contributing areas
Determine soil detachment in rills
Determine transport capacity of rill flow
Balancing the total eroded soil
available and transport capacity
Determine the net erosion or deposition
T
Figure 3.4. Structure of the model.
47
-------
3.2.1. Hydrology Components
Rainfall and runoff can exert enough kinetic energy necessary
to detach the soil particles from the ground surface. The amount of
water available for the processes of erosion, and surface runoff,
resulting from rainfall excess must first be determined. Rainfall
excess is that amount of rainfall available after losses such as
interception, depression storage, evapotranspiration, and infiltration
have been subtracted. In this model, the losses due to interception,
depression storage, and evapotranspiration are assumed to be negligible
during storm periods and assumed to be included in the infiltration.
The average rainfall intensity and steady state conditions are used
for all calculations.
3.2.2. Infiltration
Infiltration is the movement of water into the soil. The infil-
tration process depends on surface condition, soil permeability, and
soil moisture content. Several methods are available in the litera-
ture, to compute infiltration capacity (Horton, 1940; Green and Ampt,
1911; Philip, 1957). In this model, the cumulative infiltration is
calculated by Philip's Equation. Philip (1957) derived an infiltration
equation as follows:
I = St1/2 + At (3.1)
48
-------
where;
I = cumulative infiltration (in),
t = time (sec),
1/2
S = sorptivity of the soil (in/sec ),
A = a parameter depending on soil water content (in/sec).
This equation gives the infiltration at a point. Since by definition
each subwatershed is represented by a point, the equation is assumed
to be applicable to this model.
Sorptivity of the soil (S) is of greater importance at short
times in the beginning of infiltration, but A-value is of greater
importance at long times. Rogowski's (1980) conclusion following his
experimental evaluation of Philip's infiltration equation indicates
that there is a poor correlation between the S and A values. The
values of both S and A can however be measured experimentally (Tricker,
1978; Rogowski, 1980).
3.3. Model Description and Its Components
The impact of raindrops hitting the soil at a high velocity is
the first step in the erosion process. Falling raindrops detach and
transport soil particles. The transport capacity of rainfall depends
on the slope of the land surface and on sloping land more than half of
the detached soil by rainfall is moved downs lope as it falls back to
49
-------
the surface. As rainfall continues, the infiltration rate decreases
and consequently the rainfall excess begins. Thus, at first a thin
sheet of surface runoff will form which is called sheet flow. Sheet
flow can remove the lighter soil particles, organic matter, and
soluble nutrients from the land. Also, sheet flow occurs rather
uniformly over the surface and in some cases it may go unnoticed.
The flow regime in sheet flow is considered to be laminar flow. If
the surface were smooth and uniformly inclined, which is seldom the
case, it is possible to have sheet flow. However, when the accumu-
lated surface water moves downslope, it rarely moves as a uniform
sheet flow over the land surface. Because, the land surface is
usually irregular. Therefore, each portion of the surface runoff
takes the path of least resistance, concentrating in depressions and
gaining in velocity as the land slope and the runoff water depth
increases. The flow regime is no longer laminar. The velocity of
the runoff and its turbulence governs the erosiveness of overland
flow. As the surface runoff increases, it is more likely that water
will concentrate and sufficient soil may be removed to form small
but well defined channels which are called rills.
3.3.1. Rill Sources
The initiation of rills can be described under two categories
which depend on the subwatershed size.
50
-------
3.3.1.a. Small Subwatersheds
The development of the rills is in fact a result of flow detach-
ment. Concentration of overland flow (sheet flow) in rills and
development of the rill pattern have not been extensively studied.
Much of the published information concerning rill flow is directed
toward predicting the detachment or erosion in the rills (Foster and
Meyer, 1972). This is usually done by comparing the flow's average
shear stress to the critical shear stress of: the soil. However, it
is clear that when rills start to develop, there should be enough
energy associated with overland flow to develop a well defined path
in an erodible soil. Since sheet flow is assumed to be a laminar
flow regime and rill flow is a turbulent one, the distinction between
these two flow regimes can be approximated in a mathematical way.
Therefore, the initiation of rill depends on the following factors:
(1) The flow regime must be turbulent.
(2) The soil in that section must be erodible (soil erodibility
factor, K, greater than or equal to 0.10).
(3) The slope steepness must be great enough to convey the
water from the rill source.
The turbulence of the overland flow regime can be found from
Reynold's Number. According to Venard (1961), if Reynold's Number is
51
-------
high (greater than 500), the flow regime is turbulent, otherwise the
flow regime is laminar. The Reynold's Number is defined by Equation
(3.2).
Reynold's Number = Vd/v (3.2)
where;
V = velocity of surface flow (ft/sec),
d = depth of flow (ft), and
2
v = kinetmatic viscosity (ft /sec)
The variation in kinematic viscosity due to temperature variation is
shown in Figure 3.5. This is important because it may cause a 100
percent change in hydraulic conductivity (Shiao, 1978), which in turn
might affect the permeability of soil and consequently the infiltra-
tion.
The above technique for finding the rill sources may not be con-
sidered suitable for all conditions. More specifically, when the size
of subwatershed is small, this technique is more useful and more accurate
in predicting rill sources. However, the degree of accuracy will decrease
as the size of subwatershed (the horizontal and vertical distances between
grid points) increases. For large subwatersheds, not only the data
may not be very good, but the assumption that the subwatershed is
52
-------
xlO'
CO
Ul
O
O
CO
O
•H
4-1
cd
B
0)
C
J L
Kinematic viscosity
» ' ' L.
-I U
30 40 50 60 70 80
Temperature in °F
90
J L
100
110
Figure 3.5. Variation of water viscosity with temperatures.
-------
homogeneous and uniform may no longer be reasonable. The soil type
may not be the same over the subwatershed area and the ground surface
slope may not represent the details of slope steepness along the
slope length. Rogowski et al. (1974) indicated that by decreasing
length of slope length reading, more details of slope steepness for
a fixed slope length can be obtained. To overcome these sources of
errors, the following approach is considered to be more useful for
larger size subwatersheds.
3.3.1.b. Large Subwatersheds
The rill source depends not only on rainfall properties, but is
also influenced by infiltration. These two parameters are helpful
to find the excess rainfall which is the main factor in determining
rill source. The other important parameters are erodibility of soil
and available positive slope. However, the suggestion is that a rela-
tionship between runoff-infiltration ratio and soil erodibility for
developed rills can represent the conditions necessary in initiating
rill formation. Figure 3.6 shows the role and magnitude of runoff-
infiltration ratio and soil erodibility factor (K) in defining the
rill source in this model. The region that runoff starts to concen-
trate as well as the one where the rill does not form have been
indicated in Figure 3.6.
This type of curve, (Figure 3.6), is likely to be different for
different states or regions of the United States. However, a
54
-------
o
U-l
c
o
c
OS
1.0
0
Rill source region
No rill source region
0.1
0.5
Soil erodibility (K)
Figure 3.6. Influence of runoff/infiltration ratio and soil
erodibility on rill source.
55
-------
relationship between runoff-infiltration ration and soil erodibility
factor for initiating the rill can be obtained experimentally or
defined theoretically. In obtaining this type of graph, the
permeability of soil may be assumed to be the ultimate infiltration
rate. The necessary information for the graph can be obtained from
soil maps and/or field measurements. However, in those areas which
lack of available data, experimental techniques for measuring all
the parameters are suggested. Based on Rogowski's (1974) result, it
is recommended that if the width or length of subwatershed (grid point
distance) is less than 30 feet the Reynold's Number approach works
well, otherwise the second approach should be used.
3.3.2. Rill Development
The drainage area (watershed) is represented by 10,000 subwater-
sheds or a 100 by 100 matrix. The generation of rill patterns is
done by a series of moves between adjacent points in the matrix. A
rill consists of a series of connected, adjacent points in the matrix.
Figure 3.7 shows a completed rill path.
Each rill starts from its source and ends up at the bottom of the
watershed hillslope or joins an existing rill. Every point within the
interior of the watershed is surrounded by eight other matrix points.
These points represent possible directions of rill progression away
from the source or previous point of rill (Figure 3.8). By assuming
and judically picking the matrix orientation, the overall downhill
56
-------
Rill source
Figure 3.7. A rill path on watershed.
57
-------
Figure 3.8. The possible directions of rill movement.
58
-------
direction is from the top to the bottom of the matrix and thus there
are only five possible step choices (Figure 3.9).
Therefore, as the rill pattern is developed downslope, each
previously established point in the rill has to be connected to one
of the five possible directions. The rill can move straight downslope,
horizontally in either direction, or downslope at 45 degrees in either
direction. The final path of rill is a sequence of moves from point
to point in the matrix.
Since the rills are only affected by topography, the direction of
each move is chosen by the steepest slope of all the five possible
directions. The elevation of each point in the matrix and the hori-
zontal distance between them are known. Therefore, slope steepness
can be calculated as follows (see Figure 3.10):
- E(k,l + Dl/L- (3.3a)
- E(k + 1,1 + 1)]/L (3.3b)
- E(k + 1, 1)]/L (3.3c)
- E(k + 1,1 - 1)]/L (3.3d)
59
-------
• •
Higher elevation
Watershed general slope
Lower elevation
Figure 3.9. Five possible directions for rill movement at Point A.
60
-------
E(k,l+D
Figure 3.10. Slope steepness in the five possible rill directions.
61
-------
- E(k, 1 - 1)]L/1 (3.3e)
where;
E(k,l) = elevation of the point in the "k"th row and "l"th
column,
LI = the horizontal and vertical distance between two adjacent
points in the matrix,
L_ = the distance between two adjacent points at 45° direction.
After completion of the first rill by moving from its source to the
downslope boundary, the second rill path is developed. This process
is repeated until all the rills have been generated.
There are several special cases in which the rill patrern is
constrained. They are as follows:
(1) Rills are not allowed to exit at the side boundaries.
(2) If S, = S™, S™ would be selected.
(3) If S1 = S2 = S3, S3 would be selected.
(4) If S = S2 = S3 = S4 = S5, S3 would be selected.
62
-------
3.3.3. Contributing Areas-Partial Area Concept
Most of the soil loss equations for prediction of erosion from
irregular slopes were developed for estimating soil erosion for the
entire watershed. The partial-area concept is a different approach.
According to the concept of partial area hydrology, watershed areas
do not produce surface runoff with uniform frequency. Therefore, only
some portions of the watershed have high potential for surface runoff.
These areas of the watershed are likely to be close to the rivers,
gullies, or the small channels (rills), and those areas within a
watershed with low potential for surface runoff or far away from the
drainage system will be less likely to erode or if erosion occurs,
there will be little or no sediment contribution to the flow system.
As mentioned earlier, rills are considered the only flow system
capable of transporting the eroded soil. Therefore, according to the
concept of partial-area hydrology, runoff in the rills is produced
only from distinct portions of the watershed, in this case interrill
areas, which are adjacent to the rills. These contributing areas are
required to have surface runoff (sheet flow) to convey the eroded
soil to rills. Therefore, all subwatersheds one node away from a
rill could contribute to the rill. If the area contributes runoff
(sheet flow), it is considered as the contributing area, otherwise,
it is a non-contributing area.
63
-------
3.3.4. Governing Equation
The watershed areas are divided between interrill areas and rill
areas and different type of erosion mechanics are involved. Erosion
occurring in rills is defined as rill erosion, and erosion occurring
on the interrill areas is defined as interrill erosion.
Since soil erosion processes involve detachment and transport
of soil material by rainfall and runoff, the mechanics of soil erosion
is composed of four subprocesses described by Meyer and Wischmeier
(1969). These processes as described earlier in Chapter 2 are:
detachment by rainfall, transport by rainfall, detachment by runoff,
and transport by runoff. However, depending on the watershed character-
istics, some of these subprocesses may be negligible.
It has been assumed here that interrill erosion is mainly due to
detachment by rainfall and transport capacity of rainfall is only
high enough to carry the detached soil material from the interrill
areas to rills. According to this assumption, rills are the only
overland flow system responsible for carrying the interrill detached
soils to the downslope of watershed and ultimately the watershed out-
let.
Young and Wiersma (1973) determined the relative importance of
raindrop impact and overland flow on the erosion process. In their
laboratory study they used soils with three textures: loam, silt
loam, and sandy loam. They used simulated rainfall with preformed
rills in a 5 ft x 15 ft plot bed. They concluded that for all three
type soils tested; 80-85% of the soil loss originating in the interrill
64'
-------
areas was transported to rills before leaving the plot. Therefore,
although rainfall impact is primarily responsible for soil detach-
ment in interrill areas, the rill flow is the main transport mechanism
of the detached particles. However, the rill flow also contributes
soil erosion and this must be added to interrill material. Their
result supports the assumption for this model.
The total detached soil material in the rills is then compared to
the sediment transport capacity to determine if sediment will deposit
or if all of the material will be transported to the next node. The
routing procedure is done for all segments of the rills from upslope
to downslope of the watershed.
The basic governing equation of the erosion process in this model
is the continuity equation for sediment transport. This equation is
as follows (Foster and Meyer, 1972):
ax = Dr + Di
where:
G = sediment load (weight/time/unit of rill width),
X = distance downslope along the rill paths,
D = detachment (or deposition) rate of rill erosion (weight/
time/unit of total area, including rill and interrill areas), and
65
-------
D. = delivery rate of particles detached by interrill erosion
to rill flow (weight/time/unit of total area).
The above equation can be explained as a mass balance technique
for each rill path. The following sections describe different equa-
tions and techniques that are used to model the interrill areas and
the rill system.
3.3.4.a. Interrill Erosion
The dominant factor in the detachment of soil particles on
interrill areas is considered to be raindrop impact. Raindrop impact
and the sheet flow associated with the excess rainfall on the inter-
rill areas transport the detached soil particles from the effective
(contributing) interrill areas to the small flow channels (rills).
Because the depths and flow rates on interrill areas are small, the
shear stress due to the thin sheet flow is small and almost no detach-
ment occurs. Therefore, the detachment capacity of sheet flow in
interrill areas can be neglected. Thus, interrill erosion is pri-
marily due to soil detachment by raindrop impact and subsequent
transport of the detached particles by shallow interrill sheet flow.
After specifying the contributing interrill areas which contribute
erosion and runoff to rills, the amount of erosion from these areas
would be calculated. The Universal Soil Loss Equation (USLE) is the
most common estimator of annual potential soil loss caused by rainfall.
66
-------
The equation expresses annual soil loss per unit area due to erosion
by rainfall which is the only dominate erosive agent in interrill
areas. The USLE was formulated by Wischmeier and Smith (1965) as:
E = (R)(K)(LS)(C)(P) (3.5)
where;
E = computed annual soil loss (tons/acre),
R = rainfall energy factor,
K = soil erodibility factor,
LS = slope-length factor,
C = cropping management (vegetative cover) factor, and
P = erosion control practice factor.
This equation is generally for large watersheds. However, it was
applied successfully to small watersheds by Rogowski and Tamura
(1970) and should be applicable for each subwatershed in this model.
One of the important factors in the Universal Soil Loss Equation
is the rainfall energy factor (R). R-values (Wischmeier and Smith,
1978) normally used in the USLE for annual values of erosion do not
apply to individual storms. Cooley (1980) presented a more general
equation relating maximum 30-minute intensity for storms of any duration
67
-------
and volume of total precipitation. His equation provides R-values
for individual storm events of the four storm types defined by the
Soil Conservation Service (SCS Technical Notes, June, 1970). These
storm types represent typical rainfall distribution used in hydrologic
studies. Cooley's Equation is of the form:
Y
- (3.6)
where;
Y = 2.119 D°-°0086,
R = rainfall energy factor,
J = total storm rainfall in inches,
D = storm duration in hours,
a and b = constants depending on the storm type.
The values of a and b are obtained (Cooley, 1980) and Table 3.1
presents these values to be used in Equation (3.6).
The LS factor in the Universal Soil Loss Equation represents the
effects of slope length and slope steepness. The general magnitudes
of the LS factor can be estimated as follows (Wischmeier and Smith,
1965):
68
-------
Table 3.1. Values of a. and b in Equation (3.6) for each SCS storm
type1.
SCS Storm Type
Coefficients
IA
I
II
IIA
12.98
15.03.
17.90
21.51
0.7488
0.5780
0.4134
0.2811
"Cooley (1980).
69
-------
LS = (L)°'5[0.0076 + 0.0054(5) + 0.00076(S)2] (3.7)
where;
L = length in foot from the point of origin of the overland flow
to the point at which runoff enters a defined channel, and
S = the average slope over the runoff length (in percent).
The soil characteristics affecting raindrop impact detachment are
difficult to quantify. They vary with time and wetting of the soil
and formation of a surface seal affect them. Due to the difficulty
to obtain the tested parameter of these factors, the soil erodibility
factor in the Universal Soil Loss Equation may be the best available
indicator of the relation between raindrop impact and the soil pro-
perties controlling detachability.
According to Foster and Meyer (1975) , the transport capacity
of interrill flow is a function of several factors. Some of these
factors can be mentioned as runoff rate (sheet flow rate), slope
steepness, roughness of the ground surface, transportability of
detached soil particles, and effect of raindrop impact. However,
Podmore and Merva (1971) indicate the thin sheet flow on the inter-
rill areas without the raindrop impact is probably able to transport
only a small load. The raindrop impact significantly increases the
transport capacity on the interrill areas. In this model, it is
assumed that all the detached soil from the contributing interrill
70
-------
areas, based on the concept of partial-area hydrology, is transported
to the rills. The following section explains how this transport
does occur mathematically.
3.3.4.b. Contributions to a Rill
The soil particles eroded in the contributing interrill areas
would be transported to rill joints. The steepest slope of the con-
tributing interrill area to the adjacent rill joint governs where the
sediment enters the rill. Since rills are assumed to be the only flow
system that carries surface runoff, the excess rainfall (sheet flow)
from contributing interrill areas would also be transported to rill
joints in the same way as the interrill erosion. After completion
of these processes, there would be two numbers corresponding to each
rill: one representing the total runoff and other the total eroded
soil available. The next task will be to route the runoff and sediment
from rill source through rill path. However, the detachment capacity
and the transport capacity of rill flow would have to be taken into
consideration.
3.3.5. Governing Equations in Rills
The sediment routing procedure in this model is based on two
processes: Balancing the sediment supply rate and the transport
capacity of rill flow. The available eroded soil from the contributing
71
-------
areas is added to the soil detachment in rills and the total is
compared to the transport capacity of rills. If the total eroded
soil is greater than total transport capacity of rill flow, only that
portion of eroded soil that can be transported by the flow would be
transported to the next section of a rill. If the total eroded soil
is less than rill transport capacity, all the eroded soil would be
transported by rill flow.
3.3.5.a. Rill Cross Section
To route runoff and sediment, it is necessary to estimate the
rill cross section. However, to be able to achieve this goal, other
parameters need to be found first. One of these parameters which is
used to determine the flow rate in each rill section is the travel
time between rill joints.
The time that it takes for runoff to travel from one joint in a
rill to the next joint is called the travel time for that section of
rill. It depends on the rill length, land slope steepness (or rill
bed slope for pre-existing rills), and soil roughness. Depending
on the available information, there are several methods for finding
the surface flow time between two joints (points) of a rill. Since
one of the objectives of this study was to develop a simple model
with readily available data, the equation recommended by Federal
Aviation Agency (1970) was used to compute travel time. The equation
is:
72
-------
vs
(3>7)
where;
T = surface flow time or travel time (min),
L = length of flow path (ft),
S = slope of flow path (in percent),
y = a composite weighted factor (always less than one).
The value for y fact9r in the above equation is equal to the C factor
in Rational Formula (R = CIA) of hydrology. However, if this factor
is not known other approximations are applicable. In this model, the
value of the factor is approximated by the volume of runoff divided
by the volume of rainfall.
The rate of flow in the rill (Q) is found by dividing total
runoff at that section of rill by the travel time. Although rills
are small channels, they can be considered hydraulically as open chan-
nels. For uniform flow in open channel system the common Manning's
Formula is valid. This equation is:
(H)2/3(SS)1/2(F) - VF (3.8)
73
-------
where;
Q = flow rate in the rill (cubic feet per second),
F = flow cross sectional area (square feet),
n = Manning's roughness coefficient for rill,
H = hydraulic radius defined as flow cross sectional area
divided by the wetted perimeter (feet) = F/W
SS = slope of the energy gradient (under steady state, it equals
the slope of rill bed),
V = mean velocity of rill flow (ft/sec).
Manning's roughness coefficient, especially for shallow depth, depends
on flow depth. According to the available data (Beasley, 1974; Ree
and Palmer, 1949), the following equation is recommended for calculating
the value of n for rills which are considered to be shallow channel
systems in bare soil. The equation is:
n = 0.018 + 0.01(d) (3.9)
where;
n = Manning's roughness coefficient,
d = flow length in rills (ft).
74
-------
It was assumed that rill cross sections have a triangular shape
with the bed sides perpendicular to each other (Figure 3.2). There-
fore, hydraulic radius (H) and rill cross section area (F) can be
expressed in terms of flow depth (d).
F = d2 (3.10)
H = F/W = (d )/2/Td) - 0.3535 d " (3.11)
Substituting Equations (3.9) and (3.11) in Equation (3.8) results in
the following equation:
0.745(SS)1/2(d)8/3 - 0.01(Q)(d) - 0.018(Q) = 0 (3.12)
where;
SS = slope steepness,
d = rill flow depth (ft), and
Q = rill flow rate (cfs).
In order to solve Equation (3.12) for d, Newton's approach
(Carnahan, 1969) is used. Knowing d, the flow depth in the rill, the
velocity of flow and width of rill can be computed.
75
-------
3.3.5.b. Rill Transport Capacity
Meyer and Wischmeier (1969) illustrated their approach for
finding the runoff transport capacity by an equation of the form:
Tf = Stf(SS)5/3(Q)5/3 (3.13)
where;
T- = runoff transport capacity,
S - = a coefficient depending on the soil's transportability
properties,
SS = slope of the energy gradient of flow, and
Q = flow rate.
Since the coefficient (S f) is not sufficiently evaluated yet,
Yalin's Equation (1963) is recommended. Yalin's Equation is applicable
based on the assumption used for its derivation. The equation assumes
that sediment motion begins when the lift force of flow exceeds a
critical lift force. After a particle is lifted from the bed, it
would be in suspension until it settles back to bed by gravity force.
The concepts in the Yalin Equation seem most applicable to conditions
of concentrated shallow flow associated with upland erosion. The
76
-------
Yalin Equation requires only two common flow parameters, hydraulic
radius and slope of the energy gradient. The transportability of a
soil depends on the particle density, particle diameter, and the
critical lift force. The Yalin Equation is:
T = 0.635(SG)pwm(SV)g6[l - log(l +a)] (3.14)
where;
a = h • <5 (3.14a)
5 = (T - T )/T (when T < T ,5=0) (3.14b)
h = 2.45(SG)~°*4(T )°'5 (3.14c)
cr
T = (SV)2/(SG - 1.0)gm (3.14d)
SV = / gH(SS) (3.14e)
T- = transport capacity of rill flow (Ibs of sediment/second/
ft. of flow width),
SV = shear velocity of flow (ft/sec),
2
g = acceleration of gravity (ft/sec ),
p = mass density of the fluid,
w
77
-------
H = hydraulic radius (ft),
SS = slope steepness,
SG = particle density,
m = soil particle diameter (ft),
T = critical lift force, which is given in the Shields diagram
(Figure 3.11).
To obtain the critical lift force (T ), the size of soil particle
is required to use Shields diagram. Since the particle size in the
rill bed is not included in the input data of this model, it has to
be generated. The soil particle size can be found by relating it
to Manning's roughness coefficient. This relationship is (Shen, 1971):
n - 0.034(m)1/6 (3.15)
where;
n = Manning's roughness coefficient,
m = particle diameter (ft).
In order to find the critical lift force (T ), Shields diagram
was recommended by Yalin to be used in Equation (3.14b) and Equation
(3.14c). However, the Modified Shields diagram (Shen, 1971) is used
78
-------
c
l-t
(D
Shields parameter (dimensionless)
p-
(P
P-
W
t"
OQ
i-t
(U
-------
here since it is easier to use in a computer model. Figure 3.12
presents the Modified Shields diagram. This graph can be approximated
to fit a straight line with the following equation:
1 U°'92
13.2
(3.16)
where ;
Z = T
cr
1
2,
v 0
.2
's'V Pw
j
- 1/3
U
1/3
(3.16a)
(3.16b)
v = kinematic viscosity of flow (Figure 3.5),
•y and y = specific weight of rill flow and soil, respectively,
W S
(lb/ft),
T = critical shear stress (lb/ft ),
cr
m = soil particle diameter (ft),
p = mass density of the rill flow which is assumed to be
equal to that of water.
80
-------
10 i
x-s
CO
tn
-------
3.3.5.C. Rill Detachment Capacity
The detachment capacity of rill flow is defined as the rate per
unit of total area at which rill flow without a sediment load can
erode particles from the soil matrix at a location on the slope. The
detachment capacity of rill flow (D ) depends strongly on the average
bed shear. It is taken as proportional to a power of the difference
between actual shear stress and a critical shear stress. The critical
shear stress is believed to be the minimum requirement for initiation
of motion of sediment grains in the rill bed. A possible expression
for detachment capacity of rill flow (D ) is (DuBoy, 1879):
D - c (T - T )" (3.17)
c c cr
where;
D = detachment capacity of rill flow (weight/time),
2
T = critical shear stress (Ib/ft ),
cr
2
T = actual shear stress (Ib/ft ),
C and « coefficients depending on the soil and fluid properties.
Applying the DuBoy's type formula to narrow channel and neglecting the
channel side effects lead us to the following equation (Morris and
Wiggert, 1971):
82
-------
(3.18)
where;
g = detachment capacity (bed load in Ib per ft of width per
s
time),
B = rill width (ft),
y = specific weight of fluid (Ib/ft ),
G/s = total detachment (Ibs of bed load), and
-------
Knowing the particle size of sediment material, critical shear stress
(T ) can be found from Shields diagram or Modified Shields diagram
(represented by Equation 3.16).
All the rill segments would then be routed using equations and
techniques discussed above. The result would be a net amount of soil
accumulated at each joint of rill. By comparing it with the eroded
soil before the start of runoff, the erosion or deposition at each
joint would be obtained.
3.4. Data Requirements
Data needed for this model can be divided into two categories:
field data and model input data. Field data includes, contour map,
soil properties, vegetation distribution, ground cover characteristics
and rainfall distribution. Model input consists of data obtained from
maps that have been converted to conform to the format requirements
for the model.
Soil characteristic maps are difficult to obtain. However, SCS
maps usually contain most of the needed information. The vegetation
distribution map is generally difficult to find for a natural water-
shed, but a field trip may yield required information. For larger
watershed aerial surveys and land satellite information are useful.
Rainfall data is usually recorded near or at the site.
The model input includes three types of information: climatolog-
ical, watershed geometry, and watershed physical characteristics.
84
-------
Climatological data include rainfall intensity, duration, and average
temperature at the beginning of the storm. Watershed geometry data
consists of slopes and watershed dimension while watershed physical
characteristics data include soil characteristics, vegetation type,
and cover conditions.
85
-------
-------
SECTION 4
APPLICATION OF THE MODEL AND RESULTS
4.1. Introduction
To test the applicability of the model, it was applied to a
watershed and the results obtained were compared with experimental
data. The flow diagram of the input-output files is shown in
Figure 4.1. The information about the watershed and the procedure
used to obtain the input data and to choose the model parameters
are described in the following sections.
4.2. The Study Area
An area selected to test the model is a 10 acre watershed in
Central Pennsylvania described in detail by Pedersen et al. (1978).
This stripmined and reclaimed land is adjacent to undisturbed land.
The site is located less than a mile northeast of Kylertown in
Clearfield County, Pennsylvania. Figure 4.2 shows the general loca-
tion of the site. It lies within the Pittsburgh Plateau section of
the Appalachian Plateau Province. The geologic system is Pennsylvanian,
characterized by cyclic sequences of sandstone, shale, coal, and clay
(Pedersen et^ ad., 1978). The north most part of the area was reclaimed
86
-------
Weather information file
Watershed
parameters
Hydrologic file
Soil-hydrology
arameters
Land use and
conservation
practice parameter,
Erosion-sediment
yield file
X1
x'soil erodibility
^ parameters
Scour and deposition
output file
Figure 4.1. Flow diagram for the input-output files.
87
-------
CO
CO
PENNSYLVANIA
Selected site
SCALE = I: 2,000,000
Figure 4.2. The location of the study area.
-------
in 1969 and the south most reclaimed and topsoiled in 1974 (Rogowskis
1980) . On the periphery of the topsoiled and nontopsoiled areas are
small pockets of undisturbed soil (Rogowski, 1980).
To test the erosion model, a 460 ft x 320 ft portion ABDC (Figure
4.3) of the site was selected because it has all the parameters neces-
sary for numerical simulation of erosion. The required data were
measured from instruments located at the site for the past few years
(Pedersen et al., 1978; Rogowski, 1980).
4.3. Input Data
The important parameters and variables determined prior to appli-
cation of the model are summarized in the following sections.
4.3.1. Topography
Topography of the site was obtained from the field survey and is
shown in Figure 4.4. Although this figure presents four foot contour
intervals, two foot contour intervals has been used in this study
for better accuracy.
4.3.2. Infiltration Parameters
To determine the runoff available for transport of eroded sedi-
. ment estimates of infiltration were needed. Rogowski (1980) used the
89
-------
250
200 •
150
100
100
150
200
Figure 4.3. Location of area ABDC in the selected watershed. (Note:
The numbers are relative coordinates to be matched with
Figures 4.5, 4.6, 4.9, and 4.16).
90
-------
to?0
oi
-------
same site for an experimental study of infiltration on mine spoil.
He measured the infiltration rates and obtained the infiltration
parameters (S and A values for Philip's (1957) equation, Equation
3.1). These parameters (shown in Figures 4.5 and 4.6) were used to
compute infiltration over the area during the study period (June 1
to September 1, 1981).
To obtain these values, Rogowski (1980) used ring infiltrometers.
In this method, the measured infiltration rate may exceed the true
rate if the flow penetrates below the bottom of the ring. Depending
on the ring diameter, the measured infiltration rate can be off by
16 to 34 percent. Tricker (1978) in his study indicated that this
error can be appropriately corrected. His correction equation for
small ring (6 inches in diameter) which has the larger error can,
assuming uniform permeability and air-dry soil, be written as:
I = 0.65 I0'46 t-°'64 (4.1)
c m
where;
I = corrected infiltration rate (in/hr),
c
I = measured infiltration rate (in/hr), and
m
t = the cumulative time (hrs).
92
-------
eso
200 -
ISO -
100 -
SO 100 ISO 2,00
I/9
Figure 4.5. Soil sorptivity (in/sec ") of the selected watershed,
in three and two dimensions.
93
-------
SSOr—i—i—i—i—i—:—i—i—i
100 ISO EDO
Figure 4.6. A-values (in/sec) for Che selected watershed, three
and two dimensions.
94
-------
This correction seems to give better results for longer times where
largest errors occur (Rogowski, 1980). Therefore, this equation was
incorporated into the model to correct the experimental infiltration
rate.
A.3.3. Precipitation
Precipitation was measured by a universal rain gage. Table 4.1
presents the storm events which occurred in Summer 1981 (June 1 to
September 1) . According to Soil Conservation Service (SCS Technical
Note, June 1970), there are four types of storm events which represent
a range of rainfall distributions commonly seen in hydrologic studies.
Figure 4.7 shows these distributions and illustrates the storm types
and their distribution observed during the study period. To evaluate
the storm type dominant at the site, the distribution of these storm
events is plotted in Figure 4.7 and their characteristics are pre-
sented in Table 4.1. Because of the short duration (less than 0.50
hrs) the other storm events are not shown. The storms observed during
the study appear to follow closely the storm type IA. Rogowski (1981)
also concluded that storm type IA is the most prevalent type of storm
in Central Pennsylvania. Therefore, storm type parameters (Table 3.1)
were used to calculate the rainfall erosivity factor (R in USLE).
These values are shown in Table 4.2.
Storm durations reported in Table 4.1 were obtained by using time
compression approximation (Reeves and Miller, 1975). The "Time
95
-------
Table 4.1. Precipitation information for Summer 1981.
Magnitude Duration
Storms (inches) (hours)
I 1.50 3.60
2 0.35 0.50
3 1.35 4.00
4 0.85 0.40
5 1.20 3.50
96
-------
100
20 40 60
% Duration
80
100
Figure 4.7. Distribution of site storms and SCS storm types.
97
-------
Table 4.2. Summarized output for storm events of Summer 1981.
Storms
1
2
3
4
5
Duration
(hrs)
3.60
0.50
4.00
0.40
3.50
Average
Intensity
(in/sec.)
0.000116
0.00019
0.000094
0.00059
0.000095
a
Erosivity
11.86
2.39
8.75
18.32
7.51
Soil Loss
(Ibs)
228
0
93
53
61
a(T/A)/(1/100 ft-T/A in/hr).
Predicted soil leaving the area ABDC.
98
-------
Compression" is a method sometimes used for typical unsteady rainfall
events in watershed modeling. In the method the assumption is that
for a given soil the maximum infiltration rate depends on the cumula-
tive infiltration and not on the rainfall-time distribution. More
specifically, it was indicated (Reeves and Miller, 1975) that for
intermittent rainfall events, if the rainfall resumes after a period
of drizzle or no rain, the water profile will achieve a shape very
similar to that which would have developed after a continuous rainfall
event that had the same total infiltration. This method was used here
to reduce computationally the effect of interruptions on rainfall
intensity, which plays a major role in soil detachment.
4.3.4. Soil Erodibility
The erodibility factor (K) of each portion of the selected site
was reported (Rogowski, 1979) and is shown in Figure 4.8. This
figure also shows the topsoiled and nontopsoiled areas of the site.
Very low values of K on nontopsoiled spoil reflect extensive erosion
pavement development and extreme rockiness. Erodibility was assumed
to be constant during the study period.
4.3.5. Watershed Cover and Management Factors
The watershed vegetation density is presented in Figure 4.9. This
information based on site quadrant survey was used to estimate cover
99
-------
Location Map
Nontopsoiled
K. = 0.03
Figure 4.8. Soil erodibility of the selected watershed.
100
-------
2SO
£00 -
100-
SO 100 1SO 200
Figure 4.9. Vegetation density (lbs/100 sq ft) of the selected
watershed, three and two dimensions.
101
-------
factor (C) (Wischmeier and Smith, 1978) for the Universal Soil Loss
Equation. In this approach, the vegetation density of 0.06 Ib/sq ft
was taken to be equal to 50% ground cover (Wischmeier and Smith, 1978)
There is no conservation practice management in this watershed.
Therefore, the practice management factor (P in Equation 3.5) was set
equal to one for the entire area.
4.3.6. Watershed General Parameters
The other parameters used in the application of the model are
tabulated in Table 4.3. They were kept constant for all parts of the
watershed. Although their numerical values were assumed, they are
representative of the watershed conditions.
4.4. Subwatersheds Grouping
As it was discussed earlier, a portion (460 ft x 320 ft) of the
site was selected for model application. This portion (area ABDC in
Figure 4.3) was divided into 1472 square subwatersheds (subareas).
Figure 4.10 shows this grouping and also indicates the numbering of
the subwatersheds from upslope to downslope (in this case from south
to north). As it is indicated, there are 32 rows with 46 subwater-
sheds in each row. Each subwatershed is 10 by 10 sq ft.
102
-------
Table 4.3. Assigned input information for the entire watershed.
Parameter
Magnitude
Unit
Overland flow density
Soil particle density
Temperature
Overland flow viscosity
Gravity acceleration
62.5
165.5
60
0.000012
32.2
lb/ff
lb/ff
ft /sec
ft/sec2
103
-------
Location Mao
Figure 4.10. Watershed ABCD and its grid system.
104
-------
-------
4.5. Results and Discussion
The model was run for each storm event reported in Table 4.1.
After determining the infiltration at each subwatershed, the available
runoff for the start of a rill was checked. The rill patterns and
contributing interrill areas were then generated. Figure 4.11 shows
the rill patterns and the contributing interrill areas generated by
the model using the input data for storm 1. Figure 4.12 shows the
eroded-soil available for transport at each rill node after the storm 1
has ended, while Figure 4.13 indicates the net amount of soil removed
(negative) or deposited (positive), by rill flow after the runoff has
ended for a section of area ABDC. The runoff period was assumed to
equal the storm period. The deposition values in Figure 4.12 may be
thought of as the total amount of sediment recovered at the indicated
points and includes the components of both sheet erosion from con-
tributing areas and rill scour. Figure 4.12 indicates that although
a great amount of soil was removed from different locations of the
area ABDC, only 22 pounds of soil actually left the area. The rest
of the detached soil was deposited before reaching the outlet boundary.
The total amount of soil leaving the area ABDC for storm 1 was
determined to be 228 pounds. Using the same procedure, the same
type of results were obtained (not shown) for other four storm events
in Table 4.1.
Table 4.2 indicates the net amount of soil leaving the area ABDC
for each storm. Also presented in this table are the rainfall
105
-------
Figure 4.11. Rill patterns and contributing areas (shadow areas)
generated by model for storm 1.
106
-------
A
p S 9 10 11 12 13
c
15 16 9
r
1
1
1
1
1
1
:
i
i
i
i
i
i
i
i
• 2
]
1
1
1
1
1
1
1
1
1
1
2
2
2
2
1
I
1-
.1
1
1
-1
1
1
1
1
1
1
1
1
1
1
1
1
3
rr
\
^
^i_
" r 5 I 1
i
2
3
e,
S
6
7
8
9
10
11
12
13
14
15
15
17
15
19
20
21
22
23
24
25
26
27
2S
2S
30
31
32
JD
Figure 4.12.
Eroded soil (Ibs) available for transport at the end
of storm I for section pqr.s of area ABCD.
107
-------
N
P8 9 10 11 12 13 14 15 16 fl 4,
-1
-1
-1
41
-1
-1
-1
-1
-1
46
1
-1
-1
-1
-1
-I
-2
169
-1
-1
-1
1674
5
-1
-1
V
-1
-1
-1
-1
-1
-1
-2
-1
-1
-1
-1
-t
-1
-1
-1
-1
-1
-1
-1
-1
-2 I
144
-2
-2
891
2
-1
4
-1
-1
-1
101
-1
15
r
\
^1
^ r s \
*
X
B
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
D
Figure 4.13. Net amount (Ibs) of soil eroded (-) or deposited (+)
at the end of runoff period after storm 1 for section
pqrs of area ABDC.
108
-------
erosivity factor (R in USLE) which were determined using the model.
This table illustrates the roles of rainfall intensity, duration,
and erosivity on the magnitude of soil loss. Storm 4 was the most
erosive among all the storm events. But, due to its short duration,
it did not produce a great amount of soil loss compared to storm 1.
According to these results, storm 1 and 2 were the most and least
effective rainfall events with regards to the total soil loss from
the area ABDC, respectively. Even though their average intensities
were similar in magnitude both the duration and soil loss differed
considerably. The final result was obtained by summing up the total
erosion and deposition for each event. Figure 4.14 shows the com-
posite areas of erosion and deposition at the end of study period.
The erosion and deposition distribution on the watershed is one
of the benefits obtained from this model. Figure 4.14 illustrates
the utility of the model for locating areas that are susceptible
to sediment erosion or deposition on the area ABDC. This type of
information may be useful for evaluating the impact of non-point
sources of pollution in an area. Figure 4.14 could also be presented
as an iso-erosion and deposition map (contours of equal erosion and
deposition). However, because of the large range in magnitude of
erosion and deposition the results were not presented as an iso-
erosion and deposition map.
To check the model result, actual data of erosion and deposition
from the watershed for the study period were needed. Erosion and
deposition was measured experimentally on the study area. These
109
-------
Erosion (5to40Ibs)
Deposition (5 to 500 Ibs)
Deposition (>5001bs)
B
i i i
\
j i
1
l i
i ix-i fxl
i i
1 1 1
T I
Pxi
!x>>3 I I
I l I i T
llxj
I I I I
f£
I I
I I I
I I
I I I I
c-:-
I I
I I
I T f
I I I
?x! 1
mm i
vx!
P
t-X-X!
SO
yj
1 l K-:C-:H
i o ••*:::*'rsir
Fxi
I !
I l M I I I I
i i i i i i i
i m m
i i i l i
i i i
i
iTTT
I i i
l
\ffl
IT
\ \ \
J Li, I.
I
1
i I I
D
Figure 4.14. Areal distribution of erosion and deposition on each
subwatershed.
110
-------
measurements were accomplished by using metal rods (pind) distributed
over the area ABDC. These erosion pins were driven into the soil and
changes in their elevation were measured with a special micrometer
system accurate to 0.01 inches. Figure 4.15 presents the details of
the pin and its measuring technique. The distribution of the pins over
the area ABDC is shown in Figure 4.16.
To illustrate the model accuracy, sections MM' and GG' of Figure
4.16 were chosen for a more detailed study. Figures 4'. 17 and 4.18
show the changes in land elevation during Summer 1981. These results
indicated that model prediction tended to follow the experimental
variations at the measured sites. In the transects MM' and GG'
erosion predicted by the model is less than the values obtained from
the measuring pins, however, deposition predicted by the model is
more than that measured'using the pins. The underprediction of
erosion by the model could be related to the fact that most of the
parameters used in the model application were assumed to be constant
for the entire study period, and throughout the study area. The over-
prediction of deposition may be due to changes in rill transport
capacity not necessarily accounted for by the specific position of
erosion pins. There is also some evidence that erosion pins tend to
overpredict actual erosion and underpredict deposition (Sams, 1982).
To evaluate the model performance, a comparison between the pre-
dicted and recorded sediment load for the watershed was necessary.
Since the measured and predicted values coincided only at thirteen
points within the study area, a kriging procedure (Sampson, 1978) was
111
-------
Measuring points
use 6 inches
dial caliper
' 1
0.6 in
x 1
^^
.- —
-~ —
••••
• ••
,-— -
-*- •
/I I
^ i
^ I
/I C
I
K/
^
^
V
$
&
- — • —
^— -^
^40 in -JL ^L
0.2 in stove bolt
Bar stock - 1 in diameter
Bar stock drilled 3 in
Hole size approximately
0.6 in
Flat washer
rf Spoil surface
Erosion pin
Figure 4.15. Details of erosion pin.
112
-------
250
200
150
100
50
50
100
150
200
Figure 4.16. Erosion pins (•) distribution over the watershed ABDC.
113
-------
75
0)
0)
70
(0-
65
c
n
o
60L
Ground elevation before the storms
Predicted elevation after the storm
Measured elevation after the storm
M
.c
o
c
o
o
c
60 ~
co
40 j>
20 m
TJ
0 g
p
o
c
o>
D)
C
(0
.c
O
•
Figure 4.17. Erosion and deposition details of section MM.
-------
———— Ground elevation before the storms
Predicted elevation after the storms
Measured elevation after the storms
100
90
q
o
00
70
0)
W 60
•a
0 50
o
a ^o
- o
100
80
60
40
20
0
o
c
o
o
t-t
V.
1™^
CJ
0
•^«
d
0)
w
•a
a
3
0
O
c
bo
C
id
0
Figure 4.18. Erosion and deposition details of section GG.
-------
used to get best estimates of measured values at other points where
no recorded data was available and, according to the model, signif-
icant erosion or deposition was occurring. The resulting correlation
(r = 0.88) and regression lines between predicted and measured values
of the erosion and deposition (cumulative for all five storms) over
area ABDC are shown in Figure 4.19. In this figure the predicted
values are about 1.25 times the measured values. Thus, for the study
area the model was within plus or minus 25% of the measured erosion
and deposition.
The model actually goes one step further and predicts the sedi-
ment load at the stream or watershed outlet. However, because there
was no data available at the watershed outlet it was not possible to
check these values.
116
-------
0.15
0.10 -
w
03
r-i
'o
03
i_
3
rf
CD
0.05 -
0.00 -
-0.05
-0.10
-0.05
0.00
0.05
0.10
Predicted, inches
0.15
Figure 4.19.
Comparison of the predicted and measured erosion and
deposition for watershed ABDC. (Note: Curve a_ is for
measured/predicted values and curve b_ is for predicted/
measured values)
117
-------
-------
REFERENCES
Baver, L. D. 1965. Soil Physics. Third edition, 5th printing,
John Wiley and Sons, Inc., New York, London, Sydney, 489 pp.
Beasley, R. P. 1974. Erosion and sediment pollution control, The
Iowa State University Press, Ames, Iowa, 320 pp.
Bennet, J. P. 1974. Concepts of mathematical modeling of sediment
yield, Water Resources Res. 10(3):485-492, June,
Brant, G. H. , et_ al^. 1974. An economic analysis of erosion and
sediment control analysis for watersheds undergoing urbanization,
The Dow Chemical Co., Midland, Michigan.
Carnahan, B., H. A. Luther, and J. 0. Wilkes. 1969. Applied numerical
method, John Wiley and Sons, Inc.
Cooley, K. R. 1980. Erosivity values for individual design storms,
Journal of the Irrigation and Drainage Division 135-145.
David, W. P. and C. E. Beer. 1975. Simulation of soil erosion-
part I, Development of a mathematical erosion model, Transactions
ASAE 126-129, 133.
DuBoys, P. 1879. Etudes de regime du Rhone et L'action exercee'
par Les eaux sur un Lit a fond de graviers indefiniment
affouillable: Annales des points et Chaussees, Ser. 5, 18,
141-95, cited by Raudkivi, A. J., 1967, Loose boundary hydraulics,
Pergamon Press, New York, London, Sydney, 330 pp.
Ellison, W. D. 1945. Some effects of raindrops and surface-flow on
soil erosion and infiltration. Trans. Am. Geophy. Union 26(3):
415-429.
Einstein, H. A. 1950. The bed load function for sediment transpor-
tation in open channel flows, USDA, Technical Bulletin No. 1026.
Federal Aviation Agency. Airport Drainage, Department of Transporta-
tion, Advisory Circular, A/C 150-5320-5B, Washington, D.C., 1970.
Foster, R. L. and G. L. Martin. 1969. Effect of unit weight and slope
on erosion, J. Irrigation and Drainage Div., Proceeding of the
ASAE 95(IR4):551-561.
118
-------
Foster, G. R. and L. D. Meyer. 1972. A closed-form soil erosion
for upland areas, Sedimentation, edited by H. W. Shen, Water
Resources Publication, Fort Collins, Colorado.
Foster, G. R. and C. A. Onstad. 1973. Erosion equation derived from
modeling principles, Paper No. 73-277 presented at 1973 Annual
Meeting, Amer. Soc. Agricultural Engineers, Lexington, Kentucky.
Foster, G. R. and W. H. Wischmeier. 1974. Evaluating irregular slope
for soil loss prediction, Transactions ASAE 17 (2)-.305-309.
Green, W. H. and G. A. Ampt. 1911. Studies on soil physics: I. The
flow of air and water through soils, J. Agri. Sci. 4, 1-24.
Holtan, H. N. 1961. A concept for infiltration estimates in water-
shed engineering, U.S. Dept. Agr. ARS 41-51, 25 pp.
Horton, R. E. 1939. Analysis of runoff-plot experiments with varying
infiltration capacity, Trans. Am. Geophys, Union 20, 693-711.
Horton, R. E. 1940. An approach towards a physical interpretation of
infiltration capacity, Soil Sci. Soc. Am. Pro. 5, 399-417.
Kuh, H. C., Reddell, and E. A. Hiler. 1976. Two-dimensional model
of erosion from a watershed, ASAE Paper No. 76-2539.
Linsley, R. K. 1971. A critical review of currently available
hydrologic methods for analysis of urban storm water runoff.
Report to the Office of Water Resources Research, Hydrocomp
International.
Meyer, L. D., G. R. Foster, and S. Nikolor. 1975. Effect of flow
rate and canopy on rill erosion, Transactions ASAE 18(5):905-911.
Meyer, L. D. and W. H. Wischmeier. 1969. Mathematical simulation
of the process of soil erosion by water, Transactions ASAE
756-758, 762.
Meyer-Peter, E. and R. Muller. 1948. Formulae for bed-load trans-
port, Proc. 2nd Meeting IAHR, Stockholm, pp. 39-64.
Morris, H. M. and J. M. Wiggert. 1972. Applied Hydraulics in
Engineering, John Wiley and Sons, Inc., Second Edition.
Middleton, H. E. 1930. Properties of soils that influence soil
erosion. Dept. of Agr. Tech. Bull. 178.
Musgrave, G. W. 1947. The quantitative evaluation of factors in
water erosion: A first approximation, J. of Soil Water Conserv.,
2(3):133-138, July.
119
-------
Negeve, M. 1967. A sediment model on a digital computer, Department
of Civil Engineering, Stanford University, Technical Report
No. 76.
Onstad, C. A. and G. R. Foster. 1975. Erosion modeling on a water-
shed, Transactions ASAE 288-292.
Pedersen, T. A., A. S. Rogowski, and R. Pennock, Jr. 1978. Comparison
of some properties of minesoils and contiguous natural soils, U.S.
Environm. Protection Agency, Paper No. EPA-600/7-78-162.
Philip, J. R. 1957. The theory of infiltration: 4. Sorptivity and
algebraic infiltration equation, Soil Sci. Soc. Am. J.
Podmore, T. H. and G. E. Merva. 1971. Silt transport by thin film
flow, Transactions ASAE, pp. 1065-1067, 1072.
Robertson, A. F., et_ al. 1964. Runoff from impervious surface under
conditions of simulated rainfall, Transactions ASAE 9(3):343-346.
Ree, W. 0. and V. J. Palmer. 1949. Flow of water in channels pro-
tected by vegetative lining, U.S. Soil Conservation Service,
Technical Bulletin, 967.
Rogowski, A. S. 1979. Development of erosion pavement on strip mine
spoils, Amer. Soc. of Agr. Eng. 1979 Winter Meeting, Dec. 11-14,
New Orleans, Louisiana.
Rogwoski, A. S. 1980. Hydrologic parameters distribution on a mine
spoil, In Proceedings, Symposium on watershed management '80,
July 21-23, Boise, Idaho, pp. 764-780.
Rogowski, A. S. 1981. Factors affecting erosion on mine spoil, In
Proceedings of water forum '81, August 10-14, San Francisco,
California, pp. 656-663.
Rogowski, A. S., E. T. Engman, and E. L. Jacoby, Jr. 1974. Transient
response of a layered, sloping soil on natural rainfall in the
presence of a shallow water table: Experimental results, USDA,
Technical Note; ARS-NE-30.
Rogowski, A. S. and Tsuneo Tamura. 1970. Environmental mobility of
cesium-137, Radiation Botany 10:35-45.
Sampson, R. J. 1978. Surface II Graphics System, Kansas Geological
Survey, Lawrence, Kansas.
Sams, J. I. 1982. Erosion of Strip Mine Lands, Unpublished Paper
in Environmental Pollution Control, The Pennsylvania State
University, University Park, Pennsylvania.
120
-------
Shen, H. W. 1972. Sedimentation. Colorado State University, Fort
Collins, Colorado.
Shiao, L. Y. 1978. Water and sediment yield from small watersheds,
Ph.D. Dissertation, Colorado State University, Fort Collins,
Colorado.
Shields, A. 1936. Anwendung der Aehnlichkeitsmechanlk und der
turbulenzforschung auf die geschiebebewegung, Mitteilung der
Preussischen Versuchsanstalt fuer Wasserbau and Schiffbau,
Heft 26, Berlin.
Simons, D. B., R. M. Li, and M. A. Stevens. 1975. Development of
models for predicting water and sediment routing and yield from
storms on small watersheds, USDA, Forest Service, Rocky Mountain
Forest and Range Experiment Station.
Simons, D. B., R. M. Li, and T. J. Ward. 1977. Simple procedural
method for estimating on-site erosion, Prepared for USDA
Forest Service, Rocky Mountain Forest and Range Experiment
Station, Flagstaff, Arizona.
Smith, D. D. and A. W. Zingg. 1945. Investigation in the erosion
control and reclamation of eroded Shelby and related soils at
the Conservation Experiment Station, Bethany, Missouri, 1930-1942,
USDA Technical Bull. No. 833, 175 pp.
Smith, D. D. and W. H. Wischmeier. 1962. Rainfall erosion: Advances
in Agronomy, Academic Press, Inc., New York, Vol. 14, pp. 109-148.
Soil Conservation Service (SCS). 1970. Estimating peak discharge
for watershed evaluation storms and preliminary designs,
Technical Service Center, Technical Note Hydrology PO-2, USDA,
June, pp. 1-6.
Tricker, A. S. 1978. The infiltration cylinder: Some comments on its
use, J. Hydrology 36:383-391.
Venard, J. K. 1961. Elementary Fluid Mechanics, Fourth Edition,
John Wiley and Sons, Inc., New York.
Williams, J. R. 1974. Predicting Sediment Yield Frequency for
Rural Basins to Determine Man's Effect on Long-Term Sedimentation,
_In Effects of Man on the Interface of the Hydrological Cycle with
the Physical Environment, Inter. Assoc. of Hydro. Sci., Publica-
tion #113, pp. 105-108.
Williams, J. R. 1975. Sediment routing for agricultural watersheds,
Water Resour. Bull. 11(5) :965-974, October.
121
-------
Wischmeier, W. H. and D. D. Smith. 1958. Rainfall energy and its
relationship to soil loss, Transactions of AGU 29(2):285-291,
April.
Wischmeier, W. H. and D. D. Smith. 1962. Soil loss estimation as
a tool in soil and water management planning, Int. Assoc. Sci.
Hydrol. Publ. 59, pp. 148-159.
Wischmeier, W. H. and D. D. Smith. 1965. Predicting rainfall erosion
losses from cropland east of the Rocky Mountains, Guide for
selection of practices for soil and water conservation,
Agricultural Handbook 282, USDA, 47 pp.
Wischmeier, W. H. and D. D. Smith. 1978. Predicting rainfall erosion
losses—A guide to conservation planning, Agric. Handbook No. 537,
Science and Education Administration, U.S. Department of Agricul-
ture in cooperation with Purdue Agricultural Experiment Station,
pp. 1-58.
Yalin, Y. S. 1963. An expression for bed-load transportation, Pro-
ceedings of the ASCE, J. of the Hydraulic Division 89(HY3) :
221-250.
Young, R. A. and C. K. Mutchler. 1969. Soil movement on irregular
slopes, Water Resour. Res. 5(1) :184-187.
Young, R. A. and J. L. Wiersma. 1973. The role of rainfall impact
in soil detachment and transport, Water Resour. Res. 9(6):
1629-1636.
Zingg, A. W. 1940. Degree and length of land slope as it affects
soil loss in runoff, Agricultural Engineering 21(2):59-64.
122
-------
-------
APPENDIX A
EROSION SEDIMENT YIELD MODEL "KEM"
Note: KEM would be available from the Civil Engineering Department
at The Pennsylvania State University and from USDA-ARS,
University Park.
123
-------
c **************************************
C *THIS PROGRAM DETERMINES THE RILLS*
C * PATTERN,CONTRIBUTING ZONES,EROS ION *
C * AND DEPOSITION OVER THE WATERSHED *
C **************************************
C
c
c
C INPUT PARAMETERS AND THEIR DEFINITION
C *************************************************************
C
c
C E=ELEVATION OF EACH POINT (OF THE GRID SYSTEM) IN FT.
C S = SORPTIVITY (IN/SQ ROOT OF SEC.)
C A = A - VALUE FACTOR (INCH/SEC.) TO BE USED IN PHI LIP'S(1957) EQ
C L1=THE DISTANCE BETWEEN THE GRID POINTS(IN FT.)
C RR=AVERAGE RAINFALL INTENSITY (IN./SEC.)
C KK=NUMBER OF GRID POINTS IN VERTICAL DIRECTION (NO. OF ROWS)
C LL=NUMBER OF GRID POINTS IN HORIZONTAL DIRtfCTION(NO. OF COL.)
C T1=DESIRED TIME AFTER THE START OF RAINFALL'(SEC.)
C T=TOTAL DURATION OF RAINFALL(HOURS)
C P=TOTAL AMOUNT OF RAINFALL (INCHES)
C D & B - COEFFICIENTS DEPENDING ON THE TYPE OF STORM(SCS TYPES)
C (SCS STORM TYPES ARE: I, IA, II, IIA)
C C= CROP - COVER FACTOR IN UNIVERSAL SOIL LOSS EQUATION(USLE)
C PP=PRACTICE MANAGEMENT FACTOR IN USLE
C KE-SOIL ERODIBILITY FACTOR TO BE USED IN USLE(TONS/ACRE/EI )
C TIM-THE TIME AFTER THE START OF RUNOFF WHICH EROSION IS WANTED '
C US=SPECIFIC WEIGHT OF SOIL
C VIS=VISCOSITY OF RUNOFF(WATER)
C SG=SPECIFIC GRAVITY OF SOIL
C
-------
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
OUTPUT PARAMETERS AND THEIR DEFINITION
**************************************************************
MMM-THE MATRIX SHOWING RILL SOURCE,RILL PATTERN,CONTRIBUTING AREA
USLE-THE ERODED SOIL AVAILABLE ALONG THE R'lLL BEFORE ROUTINC(LBS)
RUNO-THE RUNOFF AVAILABLE ALONG THE RILL BEFORE ROUTING(CUB 1C IN.)
CFAA=THE C-COEFFICIENT FOR FED. AVIATION AGENCY'S EQUATION
NET=THE ERODED SOIL ALONG THE RILL AFTER ROUTING (LBS)
FNET-THE FINAL SCOUR AND DEPOSITION (LBS)
DIMENSION
****************************************************************
DIMENSION £(100,100),8(100,100),A(100,100),MMM(100,100)
DIMENSION CFAA(100,100),FINAL(100,100),RAIN(100,100)
DIMENSION EBP(100,100),NET(100, 100),OF(100,100) ,RUNO(100,100)
DIMENSION ASS(100,100),RKCP(100,100),USLE(100,100),KE(100,100)
DIMENSION BUD GET (100, 100) ,FNET(100 , 100 ) , C( 1 00 , 100 ) , PP ( 1 00 , 100 )
REAL L1.L2.R,RR,A1,A2,A3,AA,A5,T1,T,KE, LS . LBAR ,MIN , NM, LE , NET
REAL I , 11,12,13,14
-------
c
c
C READ THE INPUT DATA
C
C
C
READ,L1,T1,RR,KK,LL,P,T,D,B,TIM,US,VIS,SG
DO I K-l.KK
1 READ,(E(K,J),J-1,LL)
DO 2 K-l,KK
2 READ,(S(K,L),L-l,LL)
DO 3 M = l,KK
3 READ,(A(M,N),N-1,LL)
DO 44 K-l.KK
44 READ,(C(K,L),L-l,LL)
DO 55 K-l,KK
55 READ,(KE(K,L),L-i,LL)
DO 65 K-l ,KK
65 READ,(PP(K,L),L-l,LL)
C
C PRINT THE INPUT INFORMATION
C
PRINT 6000
PRINT,'THE INPUT INFORMATION CAN BE SUMMARIZED AS
PRINT,'THE RAINFALL DURATION (HOURS) -',T
PRINT,'THE RAINFALL AMOUNT (INCHES) -',P
PRINT,'THE DISTANCE BETWEEN NODES (FT)-',LI
PRINT,'THE SOIL PARTICLE DENSITY -',US
-------
to
c
c
c
c
c
45
46
C
C
C
C
c
c
c
c
33
C
C
C
DETERMINE THE TOTAL PRECIPITATION(R)
IF(T1 .LT.
IF(P .EQ.
R-P
GO TO 46
R-RR*T1
MM-KK/2
NN-LL/2
T)
0.)
GO TO 45
GO TO 45
FIND THE DISTANCE(L2) BETWEEN TWO ADJACENT GRID POINTS
AT 45 DEGREE ANGLE
L2-ABS(L1*2**0.5)
PRINT 6000
PRINT,'
PRINT,'
PRINT,'
PRINT,'
************************************
*THIS PART OF OUTPUT SHOWS THE RILL*1
* SOURCES AND THEIR PATTERNS *'
************************************
111
SET ALL THE ELEMENTS IN MATRIX "MMM" EQUAL TO ZERO
THIS MATRIX WILL SHOW THE RILL SOURCES AND THEIR PATHS
DO 33 K-l.KK
DO 33 L-l.LL
MMM(K,L)-0.
SET THE PROCEDURE FOR CHECKING ALL THE POINTS
Jl-0
IF(Jl.EQ.KK) GO TO 112
J1-J1+1
J-0
J-J + 1
-------
DO 4 K-J1.KK
DO 4 L-J.LL
C
C CHECK IF THE POINT HAS NOT ALREADY BEEN IN ANY RILL
C
IF(MMM(K,L) .EQ. 1) GO TO 4
C
C
C FIND THE INFILTRATED WATER AT ALL POINTS
C
C
I-S(K,L)*T1**0.50+A(K,L)*T1
C
C CORRECT THE ACTUAL INFILTRATION-USE TRICKER'S EQUATION(1978)
C USE SUBROUTINE "TRICK" FOR THIS PURPOSE
C
CALL TRICK(I,T1,I)
R2-R-I
C CONVERT INCHES TO FEET
C
R2-R2/12.
I-I/12.
C
C FIND THE POINTS WITH EXCESS RAINFALL
C
IF(R2 .LT. 0.) GO TO 998
C
C
C CALCULATE THE AVERAGE SLOPE FOR THOSE POINTS WITH RUNOFF
C SUBROUTINE SLOPES WOULD BE USED FOR THIS PURPOSE
C
C
CALL SLOPES(K,L,KK,LL,E,L1,L2,B1,B2,B3,B4,B5,AS)
-------
c
C FIND THE POINTS THAT FLOW STARTS TO CONCENTRATE(RILL SOURCE)
C
C
IF(L1 .LE. 30) GO TO 999
C
C
C FIND RUNOFF-INFILTRATION RATIO
C
C
RIR-R2/I
IF(RIR .GT. 1) GO TO 5
IF(KE(K,L) .GT. 0.50) GO TO 6
RIR1-1.25-(2.5*KE(K,L))
IF(RIR1 .GT. RIR) GO TO 5
GO TO 998
C
C FIND AND CHECK THE REYNOLD'S NUMBER(REY)
r
C
999 REY-32.2*AS*R2**3./(3*VIS**2.)
IF(REY .GT. 500.) GO TO 5
998 IF(L.EQ.LL) GO TO 114
GO TO 4
C
C CHECK IF THE SOIL IS ERODIBLE ENOUGH
C
5 IF(KE(K,L) .GT. 0.10) GO TO 6
GO TO 998
C
C CHECK IF THERE IS A POSITIVE SLOPE
C
6 IF(AS .GE. 0.)GO TO 7
GO TO 998
4 CONTINUE
GO TO 112
-------
c
c
C PRINT THE RILL SOURCE
C
C
7 PRINT,* '
PRINT,' THE START OF RILL IS AT K-',K,' AND L-',L
PRINT,'
C
C SPECIFY THE RILL SOURCES WITH 10 IN MATRIX "MMM"
C
MMM(K,L)-10.
Jl-K
J-L
PRINT,' THE RILL WOULD PASS THROUGH THE FOLLOWING POINTS'
100 CONTINUE
C
C CALCULATE THE SLOPES OF ALL THE DIRECTION THAT RILL COULD GO
C "SUBROUTINE SLOPES WOULD BE USED"
H" £
0 CALL SLOPES(K,L,KK,LL,E,L1,L2,A1,A2,A3,A4,A5,AV)
C
C
C FIND THE AVERAGE SLOPES (ASS) FOR EACH POINT
C
C
ASS(K,L)-AV
C
C FIND THE LARGEST SLOPE (RILL WOULD MOVE IN THAT DIRECTION)
C "SUBROUTINE BIG WOULD FIND THE LARGEST SLOPE
C
CALL BIG(A1,A2,A3,A4,A5,K,L,KK,LL,III)
IF(III .EQ. 1) GO TO 101
IF(III .EQ. 2) GO TO 102
IF(III .EQ. 3) GO TO 103
IF(III .EQ. 4) GO TO 104
IF(III .EQ. 5) GO TO 105
-------
c
C PRINT THE PATTERN OF RILL MOVEMENT
C
101 L-L+1
PRINT,'K-',K,*L-',L
GO TO 110
102 K-K+1
L-L+1
PRINT,'K-',K,'L-',L
GO TO 110
103 K-K+1
PRINT,'K-',K,'L-',L
GO TO 110
104 K-K+1
L-L-1
PRINT,*K-',K,'L-' ,L
GO TO 110
105 L-L-1
PRINT,'K-',K,'L-',L
GO TO 110
C
C SPECIFY THE RILL PATHS WITH 1 IN MATRIX MMMM"
C
110 MMM(K,L)-1.
IF(K .EQ. KK) ASS(K,L)-ABS((E(K,L)-E(K-1,L))/L1)
IF(K .EQ. KK) GO TO 113
GO TO 100
113 IF(J.EQ.LL) GO TO 114
GO TO 111
112 CONTINUE
C
C
C
C FIND THE EROSION CONTRIBUTING AREAS TO THE RILLS (PARTIAL AREAS)
C
-------
DO 200 K-l.KK
DO 200 L-l.LL
c
c
c
201
202
206
207
208
C
C
C
203
C
C
C
C
CHECK ALL THE INTERRILL POINTS
IF(MMM(K,L) .EQ. 10) GO TO 201
IF(MMM(K,L) .EQ. 1) GO TO 201
GO TO 200
IF(K .EQ. 1) GO TO 202
GO TO 206
IF(L .EQ. 1) GO TO 203
IF(L .EQ. LL) GO TO 204
GO TO 205
IF(K .EQ.KK)GO TO 207
GO TO 208
IF(L .EQ.DGO TO 209
IF(L .EQ. LL)GO TO 231
GO TO 232
IF(L .EQ. 1) GO TO 233
IF(L .EQ. LL) GO TO 223
GO TO 234
DETERMINE THE INFILTRATION ON
I1-S(K,L+1)*T1**0,50+A(K,L+1 )*T
ADJACENT TO RILLS
•
INTERRILL AREAS
1
MAKE THE CORRECTION FOR ACTUALL INFILTRATION
USE SUBROUTINE "TRICK"
CALL TRICK(I1,T1,II)
IF(MMM(K,L-H ) .EQ. 10) GO TO 230
IF(MMM(K,L+1).EQ. I) CO TO 230
-------
c
C CHECK IF THERE IS A RAINFALL EXCESS
C SET THESE CONTRIBUTING AREAS AS 8 IN MATRIX "MMM1
C
IF(R .GT. II) MMM(K,L+1 )-8 .
230 I2-S(K+1 ,L)*T1**0. 50+A(K+1 ,L)*T1
CALL TRICK(I2,Tl,12)
IF(MMM(K+1,L).EQ. 10) GO TO 200
IF(MMM(K+1,L).EQ. 1) GO TO 200
IF(R .GT. 12) MMM(K+1,L)-8.
GO TO 200
204 I3-S(K,L-1)*T1**0.50+A(K,L-1)*T1
CALL TRICK(I3,T1,I3)
IF(MMM(K,L-1).EQ. 10) GO TO 230
IF(MMM(K,L-1).EQ. 1) GO TO 230
IF(R .GT. 13) MMM(K,L-l)-8.
GO TO 230
G 205 I3-S(K,L-1)*T1**0.50+A(K,L-1)*T1
CALL TRICK(I3,T1,13)
IF(MMM(K,L-1).EQ. 10) GO TO 203
IF(MMM(K,L-1).EQ. 1) GO TO 203
IF(R .GT.I3) MMM(K,L-1)-8.
GO TO 203
209 I4-S(K-1,L)*T1**0.50+A(K-1,L)*T1
CALL TRICK(I4,T1,IA)
IF(MMM(K-1,L).EQ. 10) GO TO 210
IF(MMM(K-1,L).EQ. 1) GO TO 210
IF(R .GT. IA) MMM(K-1,L)-8.
A1-A(K,L+1)
210 I1-S(K,L+1 )*T1*0.50-I-A(K,L-H)*T1
CALL TRICK(I1,T1,II)
IF(MMM(K,L4-1 ).EQ. 10) GO TO 200
IP(MMM(K,L+1).EQ. 1) GO TO 200
IF(R .GT. II) MMM(K,L+1)-8.
GO TO 200
-------
H
*
H
*
CM i— I
O
oo
CM O
H
*
co ,*
CM CO
H
*
O • - O
HOoo^ H O •
HII.Hoo
O — » US O I
O
H
O eo
HI
o a\
CM O
CM
O«
H O 00
HI
O'-^
UOi-J
H
*
co 1-4
CM CO
CM
O»
HOoo
HI
H
«
co
CM co
CM CM
HOoo
HI
» O '"x O
O x"s O *"^
O ^^ O *"^
O **"s O ^^
O *
* H W W
*^ •* * • •*"* 00
f-t i—* rf-'N x-* ^-s g_t CM
«i_i^4.^<-i * 1-1
^~» VM" 4" 4- I— I ^N V_^
* -* cr o"z.
* H W Ui S
*— o-o-s
-It H U W S
*^ nCO
{^ CM *^^ x^ s-+k
* H
W X
• -.* —»v^4-4->-l .-»_,.., H4 —> -w . . >-l ^ v^, + + _,
(Ji-it^^Hi-ii-'iisirfE-iO.JH-i^^HCM.-ii-iiiHO.-iMS^isHcOi-Ji-ijrfi^E-i
>ce >-''>-' U I ofi >-- ^ O CM • as v_* N^ o + «>-'>-'OcM + as«^N^C3CM «os«-'N-'O
«HX:s»s^ ZZ O>^ Z Z H>^ ZZ O'w ZZ ON-- ZZ
crtJZZoScrtHJZZaSHWJZZa: crtJZXceHWHJZraiHWJJZZai
I »J N^ N^r V«^ I »J >»^ X.X X^ I *-} Vw' ^X ^IMI' | |-3 ^X >^^ ^^ | |^ >,^ **** \^f | hJ ^^ X^1 ^^
co
-------
GO TO 223
200 CONTINUE
C
C
C PRINT THE RILLS PATTERNS
C
C
PRINT 6000
6000 FORMAT('l" T40 '**********************************************)
PRINT 6002
6002 FORMATC *,T39,'* THE FOLLOWING CHART SHOWS THE START OF **)
PRINT 6003
6003 FORMATC ',T39,'* RILLS (AS *) AND THE RILLS PATTERNS (AS 1)*')
PRINT 6006
6006 FORMATC ',T39,'* AND THE EROSION CONTRIBUTING AREAS TO THE *')
PRINT 6007
6007 FORMATC ',T39,'* RILLS (AS 8). ZEROS ARE THE OTHER ZONES *')
w PRINT 6008
01 6008 FORMATC ' TAG '**********************************************)
PRINT,' '
PRINT,* '
DO 222 K-l.KK
WRITE(6,6001)(MMM(K,L),L-l,LL)
6001 FORMATC ',10011)
222 CONTINUE
C
C
C CALCULATE THE OVERLAND FLOW AT ALL POINTS OF THE WATERSHED
C
C
DO 298 K-l.KK
DO 298 L-l.LL
I-S(K,L)*T1**0.50+A(K,L)*T1
-------
c
c
c
c
298
C
C
C
C
C
C
299
300
C
C
C
C
CORRECT THE INFILTRATION
USE SUBROUTINE "TRICK"
CALL TRICK(I,Tl,I)
OF(K,L)-R-I
IF(OF(K,L) .LT. 0.0) OF(K,L)-0.
CONTINUE
DETERMINE THE RAINFALL ALONG THE RILL POINTS
DO 301 K-l.KK
DO 301 L-l.LL
IF(MMM(K,L) .EQ,
IF(MMM(K,L) .EQ,
10) MMM(K,L)-1.
1) GO TO 299
SET THE INITIAL MATRICES CONDITION
RKCP(K,L)-0.
USLE(K,L)-0.
RUNO(K,L)-0.
RAIN(K,L)-R
NET(K,L)-0.
GO TO 301
IF(P .NE. 0.) GO TO 300
P-RR*T*3600.
CONTINUE
USE COOLEY'S APPROACH FOR FINDING R TO BE USED IN USLE
" SUBROUTINE COOL WOULD BE USED TO FIND "R" VALUE IN USLE'
CALL COOL(P,T,D,B,RST)
SS-ASS(K,L)*100
-------
c
c
C USE SUBROUTINE MEY TO FIND "LS" VALUE IN USLE
C
C
CALL MEY(SS,L1,LS)
RUNO(K,L)-OF(K,L)
RAIN(K,L)-R
C
C FIND THE MATRIX RKCP TO BE USED IN USLE
C
RKCP(K,L)-RST*C(K,L)*KE(K,L)*PP(K,L)
C
C FIND THE ERODED SOIL AT EACH POINT(USING THE USLE)
C
USLE(K,L)-LS*RKCP(K,L)
301 CONTINUE
M C
OJ £
C DETERMINE THE ERODED SOIL AND RUNOFF VALUES CORRESPONDING TO
C THE PARTIAL AREAS AND DISTRIBUTE THEM ALONG THE RILLS PATTERN
C
C
DO 390 K-l.KK
DO 390 L-l.LL
IF(MMM(K,L) .EQ. 0.) ASS(K,L)-0.
IF(MMM(K,L) .NE. 8) GO TO 390
LBAR-ABS(L1**0.50)
C
C USE "USLE" FOR CALCULATING INTERRILL EROSION
C
RKCP(K,L)-RST*KE(K,L)*C(K,L)*PP(K,L)
IF(L .EQ. 1) GO TO 302
GO TO 310
-------
302 IF(K .EQ. KK) GO TO 303
GO TO 305
303 IF(MMM(K,L+l) .EQ. 1) GO TO 385
IF(MMM(K-1,L+1) .EQ. 1) GO TO 304
GO TO 382
304 IF(E(K-1,L) .GT. E(K-l.L-fl)) GO TO 383
GO TO 382
305 IF(MMM(K+1,L+1) .EQ. 1) GO TO 306
GO TO 308
306 IF(MMM(K+1,L) .EQ. 1) GO TO 307
GO TO 388
307 IF(E(K+1,L) .GT. E(K+1,L+l)) GO TO 388
GO TO 387
308 IF(MMM(K+1,L) ,EQ. 1) GO TO 387
IF(MMM(K,L+1) .EQ. 1) GO TO 385
IF(MMM(K-1,L+1) .EQ. 1) GO TO 309
£ GO TO 382
00 309 IF(E(K-1 ,L+1 ) .GT. E(K-1,L)) GO TO 382
GO TO 383
310 IF(L .EQ. LL) GO TO 311
GO TO 341
311 IF(K .EQ. KK) GO TO 312
GO TO 314
312 IF(MMM(K,L-1) .EQ. 1) CO TO 384
IF(MMM(K-1,L-1) .EQ. 1) GO TO 313
GO TO 382
313 IF(E(K-1,L) .GT. E(K-1,L-1)) GO TO 381
GO TO 382
314 IF(MMM(K+1,L-1) .EQ. 1) GO TO 315
GO TO 317
315 IF(MMM(K+1,L) ,EQ.i) GO TO 316
GO TO 386
-------
c
C FIND THE POINT THAT INTERRILL ERODED SOIL IS TRANSPORTED TO
C
316 IF(E(K+l,L) .GT. E(K-H.L-l)) GO TO 386
GO TO 387
317 IF(MMM(K+1,L) .EQ. 1) GO TO 387
IF(MMM(K,L-1) .EQ. 1) GO TO 384
IF(MMM(K-1,L-1) .EQ. 1) CO TO 318
GO TO 382
318 IF(E(K-1,L-1) .GT. E(K-1,L)) GO TO 382
GO TO 381
341 IF(K .EQ. KK) GO TO 340
GO TO 319
319 IF(MMM(K+1,L-1) .EQ. 1) GO TO 320
GO TO 328
320 IF(MMM(K+1,L) .EQ. 1) GO TO 321
GO TO 326
321 IF(MMM(K+1,L+1) .EQ. 1) GO TO 322
GO TO 325
322 IF(E(K+1,L-1) ..GT. E(K + 1,L+1)) GO TO 323
GO TO 324
323 IF(E(K+1,L+1) .GE. E(K+1,L)) GO TO 387
GO TO 388
324 IF(E(K + 1 ,L-1) .GE. E(K-fl.L)) GO TO 387
GO TO 386
325 IF(E(K+1,L-1) .GE. E(K+1,L)) GO TO 387
GO TO 386
326 IF(MMM(K+1,L+1) .EQ. 1) GO TO 327
GO TO 386
327 IF(E(K+l,L-l)-E(K+l,L+l)) 386,337,388
328 IF(MMM(K+1,L) .EQ. 1) GO TO 329
GO TO 331
-------
329 IF(MMM(K-H ,L+l) .EQ. 1) GO TO 330
GO TO 387
330 IF(E(K+1,L+1) .GE. E(K+1,L)) GO TO 387
GO TO 388
331 IF(MMM(K+1,L+1) .EQ. 1) GO TO 388
340 IF(MMM(K,L-i) .EQ. I) GO TO 332
GO TO 334
332 IF(MMM(K,L+1) .EQ. 1)GO TO 333
GO TO 384
333 IF(E(K,L-1)-E(K,L+1)) 384,338,385
334 IF(MMM(K,L+1 ) .EQ. 1) GO TO 385
IF(MMM(K-1,L-1) .EQ. 1) GO TO 335
GO TO 383
335 IF(MMM(K-1,L+1) .EQ. 1) GO TO 336
GO TO 381
336 IF(E(K-1,L-1)-E(K-1,L+1)) 381,339,383
337 IF(L .GT. NN) GO TO 386
GO TO 388
338 IF(L .GT. NN) GO TO 384
GO TO 385
339 IF(L .GT. NN) GO TO 381
GO TO 383
381 ASS(K,L)-(E(KSL)-E(K-1,L-1))/L2
SS-ABS(ASS(K,L)*100)
LS-LBAR*(0.0076+0.0053*88+0.00076*88**2)
C
C SET THE EROSION AND RUNOFF AT RILL POINTS BEFORE ROUTING
C
USLE(K-1,L-1)-USLE(K-1,L-1)+LS*RKCP(K,L)
RUNO(K-1,L-1)-RUNO(K-l,L-l)+OF(K,L)
RAIN(K-1,L-1)-RAIN(K-1,L-
GO TO 390
-------
382 ASS(K,L)-(E(K,L)-E(K-1,L))/L1
SS-ABS(ASS(K,L)*100)
LS-LBAR*(0.0076+0.0053*SS-M).00076*SS**2)
USLE(K-1,L)-USLE(K-1,L)+LS*RKCP(K,L)
RUNO(K-1,L)-RUNO(K-1,L)+OF(K,L)
RAIN(K-1,L)-RAIN(K-1, L) + R
GO TO 390
383 ASS(K,L)-(E(K,L)-E(K-1,L+1))/L2
SS-ABS(ASS(K,L)*100)
LS-LBAR*(0.0076+0.0053*SS+0.00076*SS**2)
USLE(K-1,L+1 )-USLE(K-l , L-f 1 )-f LS *RKCP(K , L)
RUNO(K-1 , L+l )-RUNO(K-l , L-f 1 )+OF(K,L)
RAIN(K-1 ,L+1 )-RAIN(K-l,L+l)+R
GO TO 390
38A ASS(K,L)-(E(K,L)-E(K,L-1))/L1
SS-ABS(ASS(K,L)*100)
LS-LBAR*(0.0076+0.0053*SS-f0.00076*SS**2)
USLE(K,L-1)-USLE(K,L-1)+LS*RKCP(K,L)
RUNO(K,L-1)-RUNO(K,L-1)+OF(K,L)
RAIN(K,L-1)-RAIN(K,L-l)+R
GO TO 390
385 ASS(K,L)-(E(K,L)-E(K,L+1))/L1
SS-ABS(ASS(K,L)*100)
LS-LBAR*(0.0076+0.0053*SS+0.00076*SS**2)
USLE(K,L-H)-USLE(K,L+1)+LS*RKCP(K,L)
RUNG(K,L+1)-RUNO(K.L+1)+OF(K,L)
RAIN(K,L+l)-RAIN(K,L-fl)+R
GO TO 390
386 ASS(K,L)-(E(K,L)-E(K+1,L-1))/L2
SS-ABS(ASS(K,L)*100)
LS-LBAR*(0.0076+0.0053*SS+0.00076*SS**2)
USLE(K-fl ,L-1 )-USLE(K-H,L-l)+LS*RKCP(K,L)
RUNO(K + 1 ,L-1 )-RUNO(K-fl ,L-1 )+OF(K,L)
RAIN(K-H ,L-1 )-RAIN(K-H ,L-
-------
CO TO 390
387 ASS(K,L)-(E(K,L)-E(K+1,L))/L1
SS-ABS(ASS(K,L)*100)
LS-LBAR*(0.0076+0.0053*SS+0.00076*SS**2)
USLE(K+1,L)-USLE(K+1,L)+LS*RKCP(K,L)
RUNO(K+1,L)-RUNO(K+1,L)+OF(K,L)
RAIN(K-H ,L)-RAIN(K+1, L)+R
GO TO 390
388 ASS(K,L)-(E(K,L)-E(K+1,L+l))/L2
SS-ABS(ASS(K,L)*100)
LS-LBAR*(0.0076+0.0053*88+0.00076*SS**2)
USLE(K+1, L+l)-USLE(K+l,L+l)+LS*RKCP(K,L)
RUNO(K+l,L+l)-RUNG(K+l,L+l)+OF(K,L)
RAIN(K+1,L+1)-RAIN(K+1,L+l)+R
GO TO 390
390 CONTINUE
C ADJUST THE UNITS OF VARIABLES
C
DO 392 K-l.KK
DO 392 L-l.LL
USLE(K,L)-((USLE(K,L)*L1**2)/43539.24)*2000
RUNO(K,L)-(RUNO(K,L)*L1**2)*144.
RAIN(K,L)-(RAIN(K,L)*L1**2)*144
C
C FIND THE "C-COEFFICIENT" FROM RARIONAL FORMULA APPROACH
C
CFAA(K,L)-RUNO(K,L)/RAIN(K,L)
C
C SET AVAILABLE EROSION BEFORE POINTS(EBP)TO ZERO
C
EBP(K,L)-0.
392 CONTINUE
-------
-p-
u>
C
C
C
6009
6010
6004
391
C
C
C
6011
6012
PRINT THE DETACHED SOIL AT EACH POINT BEFORE ROUTING
PRINT 6000
PRINT 6009
FORMATC ',T40,'* THE ERODED SOIL AVAILABLE ALONG THE RILLS*')
PRINT 6010
FORMATC ',T40,'* {{{{THE UNIT IS IN LBS }}} *')
PRINT 6008
PRINT,' '
PRINT,' '
DO 391 K-l.KK
PRINT,' FOR ROW NO.' ,K
WRITE(6,6004)(USLE(K,L)(L-lfLL)
FORMATC '.10F13.0)
PRINT,'
PRINT,' '
PRINT,' '
PRINT,'
PRINT THE AVAILABLE RUNOFF AT EACH POINT OF RILL BEFOR ROUTING
6005
PRINT 6000
PRINT 6011
FORMATC ',T40,'* THE
PRINT 6012
FORMATC ',T40,'*
PRINT 6008
PRINT,' '
PRINT,' '
PRINT,' '
DO 393 K-l.KK
PRINT,' FOR ROW NO.' , K
WRITE(6,6005)(RUNG(K,L),L-l,LL)
FORMATC '.10F10.2)
RUNOFF AVAILABLE ALONG THE RILLS *')
THE UNIT IS IN CUBIC INCHES *')
-------
PRINT,' '
PRINT,' '
PRINT,' '
393 PRINT,' '
C
C PRINT THE "C-COEFFICIENT" TO BE USED IN " FAA " EQUATION
C
PRINT 6000
PRINT 6013
6013 FORMATC ',T40,'* THE C-COEFFICIENT FOR FED, AV. AGENCY
PRINT 6008
DO 394 K-l.KK
PRINT,'ROW NO.',K
WRITE(6,6014)(CFAA(K,L),L-1,LL)
6014 FORMATC '.10F10.1)
PRINT,' '
PRINT,' '
PRINT,' '
394 CONTINUE
DO 395 K-l.KK
DO 395 L-l.LL
NET(K,L)-USLE(K,L)
395 CONTINUE
C
C
C START THE SEDIMENT AND RUNOFF ROUTING PROCEDURE
C -ROUTING WOULD BE DONE FROM DOWNSLOPE TO UPSLOPE-
C
C
OUTE-0.
DO 502 KKK-1,KK
DO 502 LLL-l.LL
K-(KK+1)-KKK
L-(LL+1)-LLL
-------
c
c
c
c
c
c
c
400
401
t_n
IF(MMM(K,L) .NE. 1) GO TO 502
RUN-RUNO(K.L)
CFAC-CFAA(K.L)
IF(K .EQ. KK) GO TO 502
SET TIME EQUAL TO ZERO FOR ROUTING
TFLOW-0.
TIME-0.
FIND THE SLOPES AND DIRECTION THAT ROUTING WOULD BE CONDUCTED
ALONG THE RILL PATHS
IF(L .EQ. 1) GO TO 401
GO TO 402
IF(A1 .LT.O.)A1-0.
A2-(E(K,L)-E(K+1,L+1))/L2
402
403
404
IF(A2
A3-(E(
IF(A3
A4-0.
A5-0.
GO TO
IF(L .
GO TO
Al-0.
A2-0.
GO TO
A1-(E(
IF(A1
.LT
K,L
.LT
406
EQ.
404
405
.LT
. 0.
)-E(
. 0.
LL)
)-E(
. 0.
)A2-0.
) A3-0
GO TO
K.L+1 )
) Al-0
403
A2-(E(K,L)-E(K+1,
IF(A2 .LT. 0.) A2
))/L2
-------
405 A3-(E(K,L)-E(K-H,L))/L1
IF(A3 .LT. 0.) A3-0.
A4-(E(K,L)-E(K-H,L-1))/L2
IF(A4 .LT. 0.) A4-0.
A5-(E(K,L)-E(K,L-1))/L1
1F(A5 .LT. 0.) A5-0.
406 CONTINUE
CALL BIG (A1,A2,A3,A4,A5,K,L,KK,LL,II)
IF(II .EQ. 1) GO TO 481
IF(II .EQ. 2) GO TO 482
IF(II .EQ. 3) GO TO 483
IF(II .EQ. 4) GO TO 484
IF(II .EQ. 5) GO TO 485
481 SL-A1
LE-L1
GO TO 486
482 SL-A2
LE-L2
GO TO 486
483 SL-A3
LE-L1
GO TO 486
484 SL-A4
LE-L2
CO TO 486
485 SL-A5
LE-L1
GO TO 486
486 IF(SL .LE. 0.0) SL-0.00001
C
C FIND THE FLOW TRAVEL TIME(TFLOW)FOR EACH SEGMENT OF RILL
C -USE SUBROUTINE TRTI FOR THIS PURPOSE-
C
CALL TRTI(CFAC,LE,SL,TFLOW)
-------
c
C DETERMINE THE FLOW RATE IN RILL SEGMENTS
C
Q-(RUN)/(TFLOW*103680.)
IF(Q .EQ. 0.)GO TO 888
C
C DETERMINE Tllti DEPTH OF RILL FLOW
C -USE SUBROURINE "NEW" FOR THIS PURPOSE-
C --SET INITIAL DEPTH EQUAL TO 0.001 FT. —
C
Yl-0.001
CALL NEW(Y1,Q,SL,DD)
GO TO 470
888 DD-0.
A70 DEPTH-DD
WIDTH-2*DEPTH
C DETERMINE THE TRANSPORT CAPACITY OF RILL FLOW
C -USE SUBROUTINE "YAL" FOR THIS PURPOSE-
C
CALL YAL(DEPTH,SL,SG,US,VIS,TF1,Y,YCR)
TF-TF1*TFLOW*60.*WIDTH
C
C DETERMINE THE DETACHMENT CAPACITY OF RILL FLOW
C -USE SUBROUTINE "DUB" FOR THIS PURPOSE-
C
CALL DUB(DEPTH,SL,US,Y,YCR,DF1)
DF-DF1*WIDTH*60.*TFLOW
C
C FIND THE TOTAL DETACHED SOIL AVAILABLE FOR TRANSPORT
C
DET-DF+EBP(K,L)+NET(K,L)
-------
C COMPARE THE TOTAL DETACHMENT WITH THE TOTAL TRANSPORT CAPACITY IN RILL
C
IF(TF-DET) 489,489,490
489 TTF-TF
GO TO 491
490 TTF-DET
GO TO 491
491 IF(II .EQ. 1) GO TO 492
IF(II .EQ. 2) GO TO 493
IF(II .EQ. 3) GO TO 494
IF(II .EQ. 4) GO TO 495
IF(II .EQ. 5) GO TO 496
C
C TRANSPOST THE NET RESULT TO THE NEXT SEGMENT OF RILL
C
492 EBP(K,L+1)-EBP(K,L+l)+TTF
MK-K
£ ML-L+l
00 GO TO 497
493 EBP(K+1, L+l)-EBP(K+l,L+l)+TTF
MK-K+1
ML-L-fl
GO TO 497
494 EBP(K+1,L)-EBP(K+1,L)+TTF
MK-K+1
ML-L
GO TO 497
495 EBP(K+1,L-1)-EBP(K+1, L-1)+TTF
MK-K+1
ML-L-1
GO TO 497
496 EBP(K,L-1)-EBP(K,L-l)+TTF
MK-K
ML-L-1
-------
GO TO 497
497 IF(TF-DF) 501,501,498
498 IF(EBP(K,L) .GT. (TF-DF)) GO TO 500
NET(K,L)-NET(K,L)-(TF-DF-EBP(K,L))
IF(NET(K,L) .LT. 0.) NET(K,L)-0
EBP(K,L)-0.
GO TO 501
500 EBP(K,L)-EBP(K,L)-(TF-DF)
GO TO 501
501 K-MK
L-ML
TIME-TFLOW+TIME
IF(TIME .GE. TIM) GO TO 502
IF(K .EQ. KK) GO TO 502
GO TO 400
502 CONTINUE
DO 504 K-l.KK
DO 504 L-l.LL
NET(K,L)-ABS(NET(K,L)+EBP(K,L))
504 CONTINUE
C
C PRINT THE ROUTED ERODED SOIL
C
PRINT 6000
PRINT 6015
6015 FORMATC ',T40,'* DTHE ERODED SOIL AFTER ROUTING *')
PRINT 6016
6016 FORMATC ',T40,'* THE UNIT IS IN LBS *')
PRINT 6008
DO 503 K-l.KK
PRINT,' ROW NUMBER - '.K
WRITE(6,6017)(NET(K,L),L-l,LL)
6017 FORMATC '.10F10.0)
-------
PRINT,' '
PRINT,* '
PRINT,' '
503 CONTINUE
C
C FIND AND PRINT THE FINAL SCOUR OR DEPOSITION
C
DO 505 K-l.KK
DO 505 L-l.LL
505 FNET(K,L)-NET(K,L)-USLE(K,L)
PRINT 6000
PRINT 6018
6018 FORMATC ',T40,'* THE FINAL SCOURE OR DEPOSITION *')
PRINT 6019
6019 FORMATC ',T40,'* NEGATIVE-SCOURE , POSITIVE-DEPOSITION *')
PRINT 6020
M 6020 FORMATC ',T40,'* THE UNIT IS IN LBS * *)
o PRINT 6008
DO 506 K-l.KK
PRINT,' ROW NUMBER - ',K
WRITE(6,6021)(FNET(K,L),L-l,LL)
6021 FORMATC '.10F9.0)
PRINT,' '
PRINT,' '
PRINT,' '
506 CONTINUE
STOP
END
C
C
C THE END OF THE MAIN PROGRAM
C ****************************************************************
-------
c
c
c
C SUBROUTINES USED IN THIS MODEL ARE;
C ****************************************************************
C
C
SUBROUTINE BIG(Al,A2,A3,AA,A5,K,L,KK,LL,III)
C
C THIS SUBROUTINE SELECTS THE BIGGEST SLOPE AMONG THE FIVE SLOPES
C —VARIABLES ARE,
C A1,A2,A3,A4,& A5-SLOPES
C K.L-THE ROW AND COLUMN POSITION OF THE GRID POINT
C KK.LL-THE TOTAL ROWS AND TOTAL COLUMNS RESPECTIVELY
C III-A PARAMETER WHICH INDICATES THE BIGGEST SLOPES
C
IF(A1-A2) 10,70,38
10 IF(A2-A3) 11,18,28
11 IF(A3-A4) 12,15,15
12 IF(A4-A5) 105,104,104
15 IF(A3-A5) 105,103,103
18 IF(A3-A4) 19,22,25
19 IF(A4-A5) 105,104,104
22 IF(A3-A5) 105,103,103
25 IF(A3-A5) 105,103,103
28 IF(A2-A4) 29,32,35
29 IF(A4-A5) 105,104,104
32 IF(A4-A5) 105,34,34
34 IF(K .LT. KK/2)THEN DO
IF(L .LT. LL/2)THEN DO
GO TO 104
ELSE DO
GO TO 102
END IF
ELSE DO
IF(L .LT. LL/2) THEN DO
GO TO 102
-------
ELSE DO
GO TO 104
END IF
END IF
35 IF(A2-A5) 105,102,102
38 IF(A1-A3) 39,49,59
39 IF(A3-A4) 40,43,46
40 IF(A4-A5) 105,104,104
43 IF(A3-A5) 105,103,103
46 IF(A3-A5) 105,103,103
49 IF(A3-A4) 50,53,56
50 IF(A4-A5) 105,104,104
53 IF(A3-A5) 105,103,103
56 IF(A3-A5) 105,103,103
59 IF(A1-A4) 60,63,66
60 IF(A4-A5) 105,104,104
63 IF(A1-A5) 105,104,104
66 IF(A1-A5) 105,68,101
68 IF(L .LT. LL/2)THEN DO
GO TO 101
ELSE DO
GO TO 105
END IF
70 IF(A2-A3) 71,91,81
71 IF(A3-A4) 72,75,78
72 IF(A4-A5) 105,104,104
75 IF(A3-A5) 105,103,103
78 IF(A4-A5) 105,103,103
81 IF(A2-A4) 82,85,88
82 IF(A4-A5) 105,104,104
85 IF(A2-A5) 105,34,34
88 IF(A2-A5) 105,102,102
91 IF(A3-A4) 92,95,98
92 IF(A4-A5) 105,104,104
95 IF(A3-A5) 105,103,103
-------
98 IF(A3-A5) 105,103,103
101 III-l
GO TO 110
102 III-2
GO TO 110
103 III-3
GO TO 110
lO'i 111-4
GO TO 110
105 III-5
GO TO 110
110 CONTINUE
RETURN
END
C
C
Ln
UJ
-------
c ****************************************************************
SUBROUTINE COOL(P1,T1,D1,B1,RST1)
C
C THIS SUBROUTINE DETERMINES THE "R" VALUE TO BE USE IN USLE
C COOLEY'S APPROACH IS USED IN THIS SUBROUTINE
C —VARIABLES ARE,
C PI-TOTAL AMOUNT OF STORM(INCHES)
C Tl-TOTAL STORM DURATION(HOURS)
C D1.B1-STORM DEPENDENT COEFFICIENTS
C RST1-"R" FACTOR IN USEL
C
POW-2. 119*T1**0.0086
RN-D1*P1*POW
RD-T1**B1
RST1-RN/RD
RETURN
END
C
C
C
C
-------
c ***************************************************************
SUBROUTINE MEY(SS1,DL1,SLS1 )
C
C THIS SUBROUTINE DETERMINES THE "LS" FACTOR TO BE USED IN USLE
C
C —VARIABLES ARE,
C SSI-SLOPE STEEPNESS
C DL1-LENGTH (FT.)
C SLS1--LS" FACTOR
C
SH-0.0076+0.0053*SS1+0.00076*SS1**2
SLS1-ABS(SH*DL1**0.50)
RETURN
END
-------
Ul
C
C
C
C
C
C
C
C
C
C
****************************************************************
SUBROUTINE TRTI(CI,LI,SI,TI)
THIS SUBROUTINE DETERMINES THE FLOW TRAVEL TIME IN RILL SEGMENT
-FEDERAL AVIATION AGENCY'S EQUATION IS USED-
--VARIABLES ARE,
CI-C-COEFFICIENT
LI-RILL LENGTH (FT.)
SI-RILL SLOPE
TI-TRAVEL TIME (MINUTES)
REAL LI
PSI-SI*100.
TI-(1.1-CI)*1.8*LI**0.50/(PSI**0.3334)
RETURN
END
C
C
C
C
-------
c ***************************************************************
SUBROUTINE NEW(Y1,Q,SL,D)
C
C THIS SUBROUTINE HELPS TO SOLVE THE NONLINEAR EQUATION
C RESULTING FROM MANNING'S EQUATION
C -NEWTON'S APPROACH IS USED IN THIS PART-
C
C —VARIABLES ARE;
C Yl-INITIAL ESTIMATE OF FOLW DEPTH(FT.)
C Q-RILL FLOW RATE (CFS)
C SL-SLOPE STEEPNESS
C D-THE RILL FLOW DEPTH (FT.)
C
Q1-0.745*SL**0.50
487 FA-(Q1*Y1**2.6667)-(Q*(0.01*Y1+0.018))
G FDA-2.6667*Q1*Y1**1.6667-0.01*Q
^ FAFDA-FA/FDA
Y2-Y1-FAFDA
IF(ABS(Y2-Y1) .LE. 0.02) GO TO 488
Y1-ABS(Y2)
GO TO 487
488 D-Y1
RETURN
END
C
C
C
C
-------
c ***************************************************************
SUBROUTINE TRICK(I,Tl,IF)
C THIS SUBROUTINE FINDS THE CORRECTED INFILTRATION
C FROM THE EXPERIMENTAL INFILTRATION
C --THICKER EQUATION IS USED HERE
C VARIABLES ARE —
C I-THE CUMULATIVE INFILTRATION(IN)
C Tl-TIME DURATION OF THE RAINFALL(SEC)
C IF-THE CORRECTED CUMULATIVE INFILTRATION(IN)
REAL I,Tl,IF,1C,IE
IE-(I*91440)/T1
IC-(3.74)*(IE**0.46)*((Tl/3600)**(-.64))
IF-(IC*T1)/91440.
RETURN
END
Ln
00
-------
Ul
VD
£ *** *************** AAA A ************ A AAAAA A ***********************
SUBROUTINE YAL(DEP,SL,SG,US,VIS,TCF,Y,YCR)
C
C THIS SUBROUTINE DETERMINES THE TRANSPORT CAPACITY OF RILL FLOW
C -YALIN'S EQUATION IS USED FOR THIS PURPOSE-
C
C --VARIABLES ARE;
C DEP=RILL FLOW DEPTH(FT.)
C SL=SLOPE STEEPNESS
C SG=SPECIFIC GRAVITY OF SOIL
C US=PARTICLE DENSITY
C VIS=RUNOFF VISCOSITY
C TCF=TRANSPORT CAPACITY OF RILL FLOW
C Y=SHEAR STRESS
C YCR=CRITICAL SHEAR STRESS
REAL NM
WID=2*DEP
NM=0.018+0.01*DEP
DIA=(NM/0.034)**6.
11R = DEP/(8**0. 50)
C
C GRAVITY ACCELERATION IS CONSIDERED TO BE 32.2 FT/SQ. SEC.
C RUNOFF DENSITYIS SET TO BE 1.98
C SPECIFIC WEIGHT OF RUNOFF IS CONSIDERED TO BE 1.0
C
V=(32.2*HR*SL)**0.50
Y=(V**2.)/((!.6500)*32.2*DIA)
XX-DIA*(((US*1.98-1.98)*32.2)/(l.98*VIS**2.))**0.3334
YY=(XX**0.92)/13.2
YYC-YY*((VIS**2.)*(((US*1.98-1.98)*32.2)**2.)*1.98)**0.3334
-------
YCR-(YYC)/(1.65*32.2*DIA)
1F( Y . LT. YCR)TIIEN DO
l)EL=0.
ELSE DO
DEL=(Y-YCR)/YCR
END IF
AA=(2.45)*((SG)**(-.40))*(YCR**0.5)
SIG=AA*DEL
1F(DEL .EQ. 0.) SIG=0.001
TCF=1.98*(SG)*DIA*V*32.2*0.635*DEL*(1-(ALOGIO(1+SIG))/SIG)
RETURN
END
C
C
C
-------
c ***************************************************************
SUBROUTINE DUB(DEP,SL,US,Y,YCR,DCF)
C
C THIS SUBROUTINE DETERMINES THE DETACHMENT CAPACITY OF RILL FLOW
C -DUBOY'S PRINCIPAL IS USED HERE-
C
C —VARIABLES ARE;
C DEP-RILL FLOW DEPTH (FT.)
C SL-SLOPE STEEPNESS
C US-PARTICLE DENSITY
C Y-FLOW SHEAR STRESS
C YCR-CRITICAL SHEAR STRESS
C DCF-DETACHMENT CAPACITY OF RILL FLOW
C
REAL NM
NM-0.018+0.01*DEP
DIA-(NM/0.034)**6.
HR-DEP/(8**0.50)
C
C FIND THE RILL FLOW VELOCITY FROM MANNING'S FORMULA
C
VEL-(1.49/NM)*(HR**0.6667)*(SL**0.5)
DCF-(10*VEL)*(Y**2.-Y*YCR)/((DIA)*(US-1.0)*62.5)
RETURN
END
C
C
C
C
-------
c ***************************************************************
SUBROUTINE SLOPES(K,L,KK,LL,E,L1,L2,A1,A2,A3,A4,A5,AS)
C
C THIS SUBROUTINE CALCULATES THE SLOPES IN THE FIVE DIRECTIONS
C SURROUNDING A POINTS AND THEIR AVERAGE
C
C —VARIABLES ARE;
C K.L-THE ROW AND COLUMN POSITION OF THE GRID POINT
C KK.LL-THE TOTOAL NUMBER OF ROWS AND COLUMNS
C E-ELEVATION
C Ll-THE HORIZONTAL DISTANCE BETWEEN GRID POINTS(FT.)
C L2-DISTANCE BETWEEN GRID POINTS(FT)AT 45 DEGREE ANGLE
C Al, A2, A3..AA, A5-SLOPES STEEPNESS
C AS-AVERAGE OF SLOPES
REAL L1.L2
DIMENSION E(50,50)
IF(K .EQ. KK) GO TO 12
IF(L .EQ. 1) THEN DO
A1-(E(K,L)-E(K,L+1))/L1
IF(A1 .LT. 0.') Al-0.
A2-(E(K,L)-E(K+1,L+1))/L2
IF(A2 .LT. 0.) A2-0.
A4-0.
A5-0.
GO TO 11
ELSE DO
IF(L .EQ. LL) THEN DO
Al-0.
A2-0.
10 A4-(E(K,L)-E(K-H9L-1))/L2
IF(AA .LT. 0.) A4-0.
A5-(E(K,L)-E(K,L-1))/L1
IF(A5 .LT. 0.) A5-0.
-------
GO TO 11
ELSE DO
A1-(E(K,L)-E(K,L+1
IF(A1 .LT. 0.) Al-0.
A2-(E(K,L)-E(K+1 ,L+1))/L2
IF(A2 .LT. 0.) A2-0.
GO TO 10
END IF
11 A3-(E(K,L)-E(K+1,L))/L1
IF(A3 .LT. 0.) A3-0.
END IF
C
C FIND THE AVERAGE
C
AS-(Al+A2+A3+A4+A5)/5.
GO TO 14
12 AS-ABS(E(K,L)-E(K-1,L)/L1)
14 CONTINUE
RETURN
END
C
C
-------
C START OF INPUT DATA
£ ****AAA**AA*AA*A*******A***A*AA**AAAAA*AAAAAA*AA*AA
//DATA.INPUT DD *
LI Tl RR KK LL P T D B TIM US VIS SG
E(l,l).E(l,2) E(l,3) E(1,LL)
!!!:!!.!!!:!!.:::::::;::::::::::::::::::::::.!!!:!'!'!
E(KK, 1) E(KK,2) E(KK,3) E(KK,LL)
S(l,l) S(l,2) S(1,LL) I
j ^ Sorptlvity Data
S(KK,1) S(KK,2) S(KK,LL) J
A(l,l) A(l,2) A(1,LL) -I
• * • • « • • j I* A-Values Data
A(KK,l) A(KK,2) A(KK,LL) J
C(l,l)C(l,2) C(l,LL) 1
j ^ Cover Factor Data
C(KK,1) C(KK,2) C(KK,LL) J
KE(l,l) KE(1,2) KE(1,LL) 1
. . . . . —^- Soil Erodlblllty Data
KE(KK,1) KE(KK,2) KE(KK.LL) J
PP(l,l) PP(1,2) PP(1,LL) I
a I ^ Practice Management Factor
PP(KK,1) PP(KK,2) PP(KK.LL) J Data
-------
APPENDIX B
SAMPLE PROBLEM OF INPUT AND OUTPUT FOR KEM
165
-------
I. Required Input Data for KEM
1. Elevation
2. K (Soil Erodibility)
3. C (Ground Cover Factor)
4. P (Practice Management Factor)
5. S (Soil Sorptivity)
6. A-value (A Soil Property Factor)
7. Subwatershed Size
8. Total Amount of Rain (Or Average Rainfall Intensity)
9. Storm Duration
10. Air Temperature
166
-------
II. Example Problem for a 10 x 10 Grid System (100 Sufawatersheds)
Li = 10 ft
Ti = 1800 sec
RR » 0.03 in/sec
KK - 10
LL - 10
P = 0 (zero if total amount of rain is not known)
T = 0.5 hrs
D = 15.03
B = 0.5780
TIM = 30 minutes
US = 2.65
VIS = 12 x 10~6 ft2/sec
SG = 2.65
The above information and other data for grid points in the
computer format follow.
Note: Free format has been used for all input data.
167
-------
OS
00
C
C START OF INPUT DATA
c *********************************
//DATA. INPUT DD *
10. 1800.0 0.03 10 10 0. 0.5 1
70 70 70 70 70 70 70 70 70
70 65 60 70 70 70 65 70 70
60 65 55 60 60 55 55 60 60
60 55 50 45 60 50 60 60 50
60
50
50
45
40
30
. 1
.8
. 1
.2
.8
.1
.8
.8
.8
.2
.01
.08
.01
.02
.08
.01
.08
.08
.08
.02
60
50
45
45
40
30
.2
.8
.8
.8
.2
. 2
.8
.8
.2
.8
.02
.08
.08
.08
.02
.02
.08
.08
.02
.08
60
40
40
40
40
30
.2
.8
.8
.8
.2
.8
.8
. 1
.8
.8
.02
.08
.08
.08
.02
.08
.08
.01
.08
.08
40
3
3
3
3
5
5
5
0
30
•
•
•
•
•
•
•
•
•
•
*
•
•
•
*
•
•
•
»
*
8
2
2
8
8
8
8
2
8
8
08
02
02
08
08
08
08
02
08
08
50
40
30
30
25
25
.8
.1
.8
.8
.8
.2
.8
.8
.8
.8
.08
.01
.08
.08
.08
.02
.08
.08
.08
.08
40
40
30
30
20
20
.8
.2
.8
.8
.8
.2
.2
.2
.8
.8
o08
.02
.08
.08
.08
.02
.02
.02
,08
.08
50
3
3
2
5
0
5
30
2
•
•
•
•
•
•
»
•
•
•
•
•
•
•
•
•
•
•
•
«
0
2
8
8
8
2
8
8
8
2
8
02
08
08
08
02
08
08
08
02
08
40
40
35
30
30
20
. 1
.8
.8
.2
.8
.8
.8
.8
.2
.8
.01
.08
.08
.02
.08
.08
.08
.08
.02
.08
50
40
35
35
30
20
.8
.8
.8
.8
.8
. 1
.8
.2
.8
.8
.08
.08
.08
.08
.08
.01
.08
.02
.08
.08
lt****4
.5.03
70 -1
70
70
60
50
40
35
35
30
20 -
.8 •
. 1
.2
.2
.8
.2
.2
.8
.8
.8 -
.08 •
.01
.02
.02
.08
.02
.02
.08
.08
.08 .
0.5780 30. 2.65 0.000012 2,65
Elevation Data
Sorptivity Data
A-Values Data
-------
•
•
•
•
*
*
•
•
•
•
*
•
•
•
•
*
1
1
1
1
1
1
1
1
1
1
1
3
3
3
3
3
3
3
3
3
/,
4
/,
/,
1
1
1
1
1
1
1
1
1
1
3
.3
.3
.3
.3
.3
.3
.3
.3
.3
/,
.4
/,
/,
1
1
1
1
1
1
1
1
1
1
•
•
•
•
*
•
•
•
•
•
1
1
1
1
1
1
1
1
1
1
1
3
3
3
3
3
3
3
3
3
4
4
/,
/,
1
1
1
I
1
1
1
1
1
1
1
.3
.3
.3
.3
.3
.3
.3
.3
.3
/,
.4
/,
/,
I
1
1
1
1
1
1
1
1
1
•
•
•
•
•
*
•
•
•
•
1
1
1
1
1
1
1
1
1
I
1
3
3
3
3
3
3
3
3
3
/,
4
4
1
1
1
1
1
I
1
1
1
1
.3 .3 .3 .3 .3 1
.3 .3 .3 .3 .3
.3 .3 .3 .3 .3
.3 .3 .3 .3 .3
.3 .3 .3 .3 .3
.3 .3 .3 .3 .3
.3 .3 .3 .3 .3
.3 .3 .3 .3 .3
.3 .3 .3 .3 .3
.3 .3 .3 .3 .3 J
.4.4.4.4.4 1
.4 .4 .4 .4 .4
.4 .4 .4 .4 .4
1 1 1
1 1
1 1
1 1
1 1
11 " ^
1 1
1 1
1 1
1 1 J
Cover Factor Data
Soil Erodibility Data
Practice Management Factor Data
-------
-------
III. Selected Output for Example Problems
170
-------
*********************************************
* THE FOILOWING CHART SHOWS THE START CF *
* RILLS (AS *) AND THE PILLS PATTERNS (AS 1) *
* 1ND THE EROSION CONTRIBUTING ARE*S TO THE *
* RILLS (AS 8) . ZERCS APE THE OTHEP ZONES *
*********************************************
****3***0g
*1 1*** •)***
1* 1** 1 1*1*
*H1*1**1*
** * 1*1 * 1**
1*11**1*1 1
81111 11811
08881811*1
00008 1681 1
OOOOH18811
-------
NO
rep sew NC.
370.
FCP ROV NO.
5554.
PtP PCH NC.
370.
tea RCW NO.
608.
PCR RON NO.
4780.
PCR 1CW NO.
170.
FCI ecu NC.
0.
rcR sew NO.
0.
FCP RCU KC.
0.
FCf> ICW NO.
0.
• THE ERODPD SOU AVAILABLE ALONG 1HE RILLS*
• ((((THE UNIT IS IH LBS (I) »
1
2574. 1150. 1041. 0. 170. 608.
2
4064. 608. 21009. 11731. 21004. 6748.
1
22599. 2574. 11221. 8058. 608. 370.
4
608. 2574. 608. 57179. 1858. 13839.
5
14681. 59781. 608. 11021. 170. 35010.
6
76h2. 1422. 170. 10109. 12891. 1422.
7
11087. 6515. 6515. 34. 370. 608.
8
0. 0. 0. '12239. 0. 18975.
9
0. 0. 0. 0. 16214. 0.
10
0. 0. 0. 0. 11432. 0.
170.
16016. 4811. 1077.
257*. 1858. 17408.
28466. 1043. 11221.
170. 17408. 4780.
6748. 257«.
370.
1422.
14.
608. 4780. 1422.
9054. 4780.
S4M40. 36201.
-------
ROW NUKBIR =
0.
0.
*•**
*
*
•*••****»#****»***********#
THE ERODED SOIL HFTER ROOTING »
THE UNIT IS IN tES *
3. 0. 0. 0.
0.
RCV HJPBER
0.
0.
0. 0. 0. 0.
0.
ROW NUMBER
nu.
0.
0. 10131. 127H. 0.
RCU NUMBER
0.
0. 24810.
37l»5.
POW MOIBfB
0.
5
10089.
0. 11618. 11012.
0. 51397.
RGW NIIMBFR
239.
6
HHUBh,
5<»219.
9663.
ROW DUMBER
0.
0. 150117. H206.
0. 3230. 6«26.
RCU NUMBER -
0.
0.
0. 0. 0. 16308.
0.
RCW NUMBER =
0.
0. 2016-iO.
FCW NUMBER =
0.
10
0. 33U32.
0. 6890M. 43069.
-------
*********************************************
* TH° PIHUL SCOUBE OR DEPOSITION *
• NECATIVE=SCOORr , POSITI»E=DEPOSrTION *
* THE UNIT IS IN IBS *
*********************************************
pew NniiPER = i
-370. -257H. -3153. -10«3. 0. -373. -608. -370. 0. 0.
ROH MOBBER = 2
-5551*. -<*06U. -608. -21009. -11733. -21009. -67UB. -16016. -U81U. -1077.
BCW NOHBEfi = 3
6U3. -21«06. -257I*. -11221. -8058. 9523. 90*. -257H. -1P58. -17H08.
BCW NUMBER = U
-6C8. -60ft. -257"*. -608. -323U9. -1858. -3009U. -28«66. -10«3. -11221.
8CW NOB3FE = 5
-U783. -11*681. -H969U. -606. -21U06. 10662. -35010. 51028. -17U06. -U780.
RCH NDHBER = 6
-131. -7662. U736S. 58880. -10109. -12893. -1H22. -67U8. "»090. -1Q22.
BOW NUHBER = 7
0. -11087. -6515. -6515. 150083. 7836. -608. 0. 2860. 6392.
BOH NnnBER = 9
0. 0. 0. 0. -12239. 0. -2667. -608. -H780. -1«22.
ROW KONBEB = 9
0. 0. 0. 0. 0. 165U1S. 0. 0. -9051. -«780.
RCH MIHBPB - 10
0. 0. 0. 0. 0. 0. 0. 0. 1U06U. 6468.
-------
TECHNICAL REPORT DATA
(Please read liiaruciions on ^ht reverse before completing}
1. REPORT NO
4. TITLE ANOSUBTITLE
Rill-interrill erosion ant
stripmine hydrology
7. AUTHOR(S)
M. R. Khanbilvardi , A. S.
2.
i deposition model of
Rogowski, and A. C. Miller
9. PERFORMING ORGANIZATION NAME AND ADDRESS
Northeast Watershed Research Center
USDA-ARS, 110 Research Building A
University Park, Pennsylvania 16802
12. SPONSORING AGENCY NAME AND ADC
U.S. Environmental Prot
Office of Research & De
Office of Energy, Miner
Washington, D.C. 20'
mess
.ection Agency
•velopment
als & Industry
160
3. RECIPIENT'S ACCESSIOfVNO.
5. REPORT DATE
6. PERFORMING ORGANIZATION CODE
8. PERFORMING ORGANIZATION REPORT N<
2
10. PROGRAM ELEMENT NO.
11. CONTRACT /GRANT NO.
EPA-IAG-D5-E763
13. TYPE OF REPORT AND PERIOD COVEREC
Interim 9/1/75-8/31/80
14. SPONSORING AGENCY CODE
EPA-ORD
15. SUPPLEMENTARY NOTES
This project is part of the EPA-planned and coordinated Federal Interagency
Energy/Environment R&D Program.
16. ABSTRACT
An erosion-sediment yield model (labeled KEM) was developed from the continuity
consideration for sediment transport and from equations describing rill and interrill
erosion. This computerized model is based on dividing the upland areas into a grid
containing rill and interrill zones and on the Universal Soil Loss Equation (USLE).
The USLE is used to estimate the sediment contribution from the interrill areas. Pre-
diction of soil loss from the interrill areas is based on the premise that both raindro
impact and overland flow energy can create soil erosion. The rill flow carries the
interrill erosion along with the rill scour. Rill transport capacity governs the
amount of removed soil from the site. If the flow transport capacity is less than the
available eroded soil, net erosion equals the transport capacity and the excess sedi-
ment is deposited in the flow paths. Otherwise, all eroded soil will move downslope
and out of the watershed.
The model was tested by simulating actual events on a small watershed in Central
Pennsylvania for summer storms during 1981. Applying the model to this stripmined and
reclaimed area created a set of information about the location and amount of watershed
erosion and deposition. The areal distribution of erosion and deposition was compared
with measured data. The model performed satisfactorily in predicting soil loss from
the site.
17. (Circle One or More)
KEY WORDS AND DOCUMENT ANALYSIS
a. DESCRIPTORS
feiib^ <^^^mnol°VY
U.nvironmeniai Engineering Combustion
Other;
Energy Conversion
Physical Chemistry
Matervals Handling
Inorganic Chemistry
Organic Chemistry
Chemical Engineering
13. DISTRIBUTION STATEMENT
b,lDENTIFlERS/OPEN ENDED TERMS '
I*.*. t™.« >~ „,.„. ^
s^r«^r.:r ir^0,r: 6"^H ^5W
*^i-i ^. r .^ i-*.
19. SECURITY CLASS (This Report)
20. SECURITY CLASS ,'Tlns pafej
c. cos AT i Field/Group
6F 8A 8F
8H 10A 10B
7B 7C 13B
21. NO. OF PAGES
22. PRICE
EPA Form 2220-1 J9-73)
------- |