7

CHAPTER 2

LITERATURE REVIEW

A broad understanding of various topics in environmental science and modeling

technology was required to complete the studies presented in this thesis, and it is
important to thoroughly review each of these. The first topic that will be reviewed is
erosion science and modeling using the most popular and current erosion prediction
models: USLE, RUSLE, and WEPP. Because the major research of this thesis deals with
topography, a large portion of the review of erosion models will focus on slope and
length factors. Subsequently, the subject of watershed modeling together with hillslope
definition and grid based hydrological modeling will be reviewed and discussed. Finally,
the use of digital elevation models for flow routing and the topic of channel network
identification will be presented. This literature review is intended to describe the
potential topographic effects of slope and length on erosion modeling and give the reader
a comprehensive view into the steps involved in modeling and interfacing environmental
models using DEMs and GIS algorithms.

2.1. Erosion research and models

The three main models used today to predict erosion in small agricultural fields are the
USLE, RUSLE, and WEPP models. A description and comparison of all three will show
how slope and length affect the prediction of erosion. Due to possible interaction
between soils, management the way slope and length are defined has an effect on the
modeling with GIS and DEMs. It is fundamental to understand the relationships between
 8

slope and length as they have an effect on soil loss. This understanding is essential to
using GIS and DEM to discretize the watershed.

2.1.1. The USLE model

The USLE (Wishmeier and Smith, 1978) is an empirical model that was developed based
on 10,000 plot-years of basic runoff and soil loss data from 49 locations across the United
States. It is designed to compute long term average soil erosion rates from sheet and rill
erosion under specified conditions.

The soil loss equation is the following:

A=RKLSCP [2.1]

A= computed soil loss per unit area (usually in tons/acre/year)

n
R = rainfall and runoff factor (usually R = ∑ EI 30 , where EI30 is the kinetic
i =1

energy of a particular storm times a maximum 30 minute intensity)

K = soil erodibility factor (soil loss rate per erosion index unit for a specified soil
as measured on a unit plot of 72.6 feet length, 6 feet width, 9 percent uniform
slope, continuously in clean-tilled fallow)
L = slope length factor (ratio of soil loss from field slope length to a 72.6 feet
length plot under identical conditions.)
S = slope-steepness factor (ratio of soil loss from field slope gradient to a 9
percent slope under identical conditions.)
C = cover and management factor (ratio of soil loss from an area with specific
cover and management to an identical area in tilled continuous fallow)
P = support practice factor (ratio of soil loss with a support practice factor such as
contouring, terracing, etc, to a straight-row farming up and down slope)
 9

The L (slope length) and S (slope steepness) factors

The L and S factors account for the runoff concentration, velocity, and erosive potential
in the USLE. The effect of slope length on annual runoff per unit area is considered
negligible, however the soil loss per unit area increases as the slope length increases.
Runoff generally increases as slope gradient increases, but soil loss increases more than
runoff. The effects of slope steepness on soil loss are greater than the slope length effects
(Wischmeier and Smith, 1978).

The L and the S factor are defined in reference to the “unit plot” which is 6 feet wide by
72.6 feet long with a uniform slope of 9%. The relationships in the equations were
derived from natural rainfall plots on cropland with slopes ranging from 3 to 18 percent
and 30 to 300 feet lengths. For erosion prediction purposes, both the L and the S factors
are combined in a single LS factor. A variety of charts have been created to compute LS
values based on slope length and gradient (Wischmeier and Smith, 1978).

The slope length is defined as the distance from the point of origin of overland flow to
the point where either the slope decreases to the extent that deposition begins, or the
runoff water enters a well defined channel (Wischmeier and Smith, 1978). This is
because USLE is not applicable to areas where net deposition occurs. Slope length can
be expressed by the following equation:

L = (λ / 72.6 m ) [2.2]

Where 72.6 represents the unit plot length in feet, λ is the field slope length in feet, and m
is the exponent related to the slope. Suggested values for m are the following:
m = 0.5 if slope >=5%
m = 0.4 if slope is greater than 3 and less than 5%
m = 0.3 if slope is between 1 and 3 %
m = 0.2 for uniform slopes < 1%
 10

These m values were further improved in the revised universal soil loss equation
(RUSLE) when enough experimental data had been gathered to better quantify m for
slopes greater than 10 percent. Additional studies had shown that soil, crop, and
management variables could affect the L and S values for higher slopes because these
factors affect the water runoff and infiltration. Snowmelt and thawing could also affect
the calculation of L and S values as snowmelt on saturated soil surfaces could cause
increased runoff and soil losses.

The slope steepness factor for USLE can be computed with the following equation:

S = 65.41sin 2 θ + 4.56 sin θ + 0.065 [2.3]

Where θ is the angle of slope in degrees.

As mentioned earlier, the S factor was also believed to be influenced by interactions with
soil properties and surface conditions. Improved S factor equations were included in the
RUSLE model when interaction effects had been adequately quantified by research data.

2.1.2. The RUSLE model

Since the introduction of USLE, additional research and improvement have been
completed and are incorporated as the RUSLE (Renard et al., 1997). Like USLE,
RUSLE is a lumped model that should be used to provide average soil loss values over a
significant time span of usually one year. RUSLE still maintains the same format as
USLE, using the R, K, L, S, C, and P factors, but substantial improvements have been
made in the way those factors are determined.

Slope length for RUSLE is defined similarly as the USLE. The slope length is the
horizontal distance between the origin of overland flow to the point where either runoff
becomes concentrated or deposition begins. It is believed that surface runoff will usually
 11

concentrate after 100 meters, but longer slopes can occur. One important problem that
arises in the USLE and RUSLE is that the final selection of slope length and even
gradient is left up to the user. It has been shown that different users can select different
length and slopes for the same area of interest resulting in different erosion estimates
(mainly because determination of the start of a deposition region can be a judgement
call). A guideline showing the best methods of selecting lengths and gradients is
presented in the USDA Agricultural Handbook Number 703 to minimize this type of
discrepancy between users for the most common type of field situations.

New equations have been implemented to calculate the LS factors which reflect the ratio
of rill to interrill erosion and other concepts. For example, the equation to calculate the
slope length factor has been kept (equation 2.2), but the following new equation to
calculate the m exponent has been created:

m = β /(1 + β ) [2.4]

The m exponent is now related to the ratio β of rill erosion to interrill erosion. Rill
erosion is mostly affected by flow and interrill erosion by raindrop impact. McCool et al.
(1989) computed the value of β for soil moderately susceptible to both rill and interrill
erosion as:

(sin θ / 0.0896)
β= [2.5]
3.0(sin θ ) 0.8 + 0.56

New tables to calculate this parameter have been created which corrects β for conditions
when the soil is highly susceptible to rill erosion and when conditions favor less rill
erosion. The factors that can influence the susceptibility to rill erosion can be surface
cover, soils, snow effects, and others. These tables can be found in the USDA
Agricultural Handbook No. 703.
 12

New slope steepness equations have also been created which reflect current research
observations. These equations are the following:

For slopes less than or equal to 9%:

S = 10.8 sin θ + 0.03 [2.6]

For slopes greater than 9%:

S = 16.8 sin θ - 0.50 [2.7]

For slopes shorter than 15 ft:

S = 3.0 (sin θ)0.8 + 0.56 [2.8]

For conditions where the soil is thawing after being recently tilled (weakened state):

S = 10.8 sin θ + 0.03 s<9% [2.9]

S = (sin θ/0.0896)0.6 s≥9% [2.10]

These slope steepness equations were proposed by McCool et al. (1987). Additional
research of these factors has been conducted for steeper slopes by Liu et al. (1994). He
proposed the following equation:

S = 21.91 sin θ -0.96 [2.11]

One equation that unifies the equations for the slope steepness factor from RUSLE and
the equation for steep slopes from Liu et al. (1994) was proposed by Nearing (1997). It
has the following form:
 13

17
S = −1.5 + [2.12]
1 + exp( 2.3 − 6.1sin θ )

This logistic equation maintains a close fit with the RUSLE equations and provides
statistically representative results for steeper slopes.

As mentioned earlier, the selection of slope length in RUSLE does not take into account
the area of deposition at the end of a hillslope. However, equations for the calculation of
deposition due to management practices such as strip cropping are present in the RUSLE
model integrated into the calculation of the P factor. The WEPP model, on the other
hand, does simulate deposition directly related to the slope values along the hillslope
profile.

2.1.3. WEPP model

The WEPP hillslope model is a physically based continuous simulation model (Flanagan
and Nearing, 1995). It uses fundamental hydrologic and erosion mechanics as opposed to
the empirically based USLE and RUSLE models. It is based on a two-dimensional
hillslope profile approach and is able to predict deposition and erosion along the profile
as well as sediment delivery from the profile.

The WEPP hillslope model does not have a slope length and slope gradient factor as in
USLE or RUSLE. Instead, the slope gradient and length inputs of WEPP are deeply
integrated into the hydrology and erosion components of the model. This means that the
measurements of slope length and slope gradient are not limited to affecting L and S
factors as in USLE and RUSLE, but instead affect the calculations of runoff, friction,
transport capacity, and various other factors.

The slope length of a WEPP profile is defined as the distance from the top of the hillslope
to the end of the hillslope (usually ending in a channel or impoundment). Selection of the
length and slope profile in WEPP is thus easier than in USLE or RUSLE for the case of a
 14

single hillslope. Since WEPP is able to calculate both detachment and deposition along
the hillslope profile, it is therefore important that an accurate representation of the slope
profile (length and gradients) be used.

2.1.4. Slope gradient and slope length in the WEPP model

Slope gradient and slope length factors are deeply integrated in two major components of
the WEPP model: the erosion component and the hydrology component. A description of
both of these components in presented as follows.

Erosion Component:

The erosion component of WEPP is based on the continuity equation:

dG
= Dr + Di [2.13]
dx

where,
G= sediment load (kg·s-1· m-1)
x = distance downslope (m)
Dr = rill erosion rate ( + for detachment, - for deposition)
Di = interrill sediment delivery (kg·s-1·m-2)

The interrill component of the continuity equation delivers sediment from the interrill
parts of the hillslopes to the rills and is presented as follows:

R 
Di = K iadj I eσ ir SDRRR Fnozzle  s  [2.14]
w

where
Kiadj = interrill erodibility
Ie = effective rainfall intensity (m·s-1)
σir = interrill runoff rate (m·s-1)
 15

SDRRR = sediment delivery ratio (function of random roughness, side slope, and particle
size distribution)
Fnozzle = irrigation adjustment factor
Rs = rill spacing (m)
w = rill width (m)

Slope and length terms can be found in the runoff rate and the sediment delivery ratio. In
the case of the sediment delivery ratio, if the average profile slope is greater than the
interrill side slope, the profile slope steepness is used for the interrill delivery slope value.

Important parts of the erosion component are the rill equations governing detachment,
transport, and deposition down the hillslope. Rill erosion is divided into two parts,
detachment and deposition. Rill detachment is predicted when the flow shear stress
exceeds the critical shear stress of the soil and when the sediment load is below the
calculated sediment transport capacity.

G
Dr = Dc (1 − ) [2.15]
Tc

where

Dc = detachment capacity by rill flow (kg·s-1 m-2)

Tc = sediment transport capacity (kg·s-1·m-1)

Dc = K r (τ f − τ c ) [2.16]

Kr = rill erodibility parameter

τf = flow shear stress acting on soil (Pa)
τc = critical shear stress of soil (Pa)
 16

Deposition in a rill will be predicted when the sediment load is greater than the sediment
transport capacity.

βV f
Dr = (Tc − G ) [2.17]
q

where,

Vf = effective fall velocity for the sediment (m·s-1)

Q = flow discharge per unit width (m2·s-1)
β = rainfall turbulence coefficient

For shallow flows, the rainfall induced turbulence factor can influence the effect of
steepness on soil loss. On steep slopes (>5%) the effect seems to be negligible, but on
smaller inclines, the rainfall induced turbulence seems to have an observable effect on the
rate of deposition (particle size and density dependent) therefore affecting the overall soil
loss (Meyer et al., 1983; Cochrane, 1995; and ongoing studies).

Slope is also used to calculate the shear stress acting on the soil (flow shear stress) and
the transport capacity of the flow.

The shear stress factor for the rill flow is computed using an equation of the following
form:

fs
τ = γR sin(α )( ) [2.18]
ft

γ = specific weight of water (kg·m-2·s-2)

R = hydraulic radius (m)
α = average slope angle for uniform segment
 17

fs/ft = ratio of friction factor for the soil to total rill friction factor (accounts for shear
stress acting on surface cover like residue, etc)

In this case the cover will influence the effect that the slope gradient and length will have
on the overall soil loss. The sensitivity to soil losses of slope gradient and length would
be diminished under conditions where the shear stress acting on the soil is decreased by
residue cover or other management practices.

For the calculations of hydraulic radius, a rectangular rill is assumed and the flow depth
in the rill is calculated using the Darcy-Weisbach equation.

8 gRS
f = [2.19]
V2

where

f = calculated friction coefficient

g = gravity (m·s-2)
R = hydraulic radius (m)
S= average slope (%)
V = velocity of flow (m·s-1)

A uniform slope gradient is assumed for each length segment in the calculations of the
flow depth in the rill. The actual friction coefficient (f) used is computed using equations
developed by Gilley and Weltz in Flanagan and Nearing (1995) using experimental data.

Transport capacity is calculated by first applying the full Yalin (1963) equation to the end
of an overland flow element (OFE) or hillslope. A simplified equation, such as the one
presented below, is then applied to the interior points along the slope profile.
 18

Tc = k tτ 3f / 2 [2.20]

Tc = sediment transport capacity (kg·s-1·m-1)

kt = transport coefficient (m0.5·s2 kg –0.5)
τf = hydraulic shear stress acting on the soil (Pa)

The transport capacity coefficient, kt, is dependent on a variety of factors including slope
steepness (or slope of energy gradeline).

Hydrology component

The hydrology component of the WEPP model produces four inputs for the erosion
component:
Peak runoff rate, Pr (m·s-1)
Effective runoff duration, tr (s)
Effective rainfall intensity, Ie (m·s-1)
Effective rainfall duration, te (s)

Hydrology is also used to provide values for the amount of water that infiltrates into the
soil for the crop growth and decomposition components which will in turn affect the
infiltration, runoff routing, and erosion parameters. The water balance component also
makes use of the hydrology component and utilizes the surface slope to determine the
direction of flow of seepage.

Infiltration is calculated using the Green-Ampt Mein-Larson (GAML) equations (Mein

and Larson, 1973). Adjustments are made for soil saturated conditions and for
depressional storage. Depressional storage is related to random roughness and the slope
of flow and was derived from an analysis of over 1000 tilled plots by Onstad (1984):

Sd = 0.112rr + 3.1r2r – 1.2rrSo [2.21]

 19

Sd = depressional storage (m)

rr = random roughness (m)
So = slope of flow surface (m·m-1)

When the infiltration rate and depressional storage have been exceeded, rainfall excess
occurs. Once rainfall excess occurs, peak discharge (Pr) is calculated using a kinematic
wave model.

The kinematic wave equations are the following:

∂h ∂q
+ =v [2.22]
∂t ∂x

v = rainfall excess rate (m·s-1)

h = depth of flow (m)
q = discharge per unit width of plane (m3·m-1·s-1)
t = time (s)
x = distance (m)
and

q = αh m [2.23]

α = depth discharge coefficient (m0.5·s-1)

m= depth discharge exponent (value of 1.5)

In which the Chezy relationship is used:

α = C So0.5 [2.24]
 20

where,

C = Chezy coefficient (m0.5·s-1)

So = slope value (m·m-1)

Two methods are used to compute the peak discharge from the kinematic wave equations.
When WEPP is run in a single event mode a semi-analytical solution is always used, and
when WEPP is run in a continuous mode, the semi-analytical solution or an
approximation of the kinematic wave model may be used depending upon the runoff
event condictions.

The hydrographs created using these methods exhibit two properties, either partial
kinematic equilibrium will be reached or a flat topped hydrograph with infinite runoff
duration will be created. Because of this WEPP will terminate the routing of rainfall
excess when either the routed runoff volume equals 95% of the rainfall excess volume or
the discharge rate is 10% of the peak discharge rate (Stone et al., 1995).

If the rainfall excess duration is less than the time to kinematic equilibrium, then partial
kinematic equilibrium will be reached. The time to kinematic equilibrium can be
calculated in the following way:

1/ m
 L 
t e =  m −1  [2.25]
αv 

where

L = length of the overland flow surface (m)

v = rainfall excess rate (m·s-1)
m = 1.5
 21

This also shows that the slope and length play a role in the time computations.

The effective runoff duration (tr) is calculated as:

Vt
tr = [2.26]
Pr

where,

Vt = total runoff volume depth (m)

Pr = peak runoff rate (m·s-1)

The effective rainfall intensity is calculated over the time period when rainfall rate
exceeds infiltration rate as:

Ie =
∫ Idt [2.27]
te

where

Ie = effective rainfall intensity (m·s-1)

I = break point rainfall intensity (m·s-1)
te = time when rainfall rate exceeds infiltration rate (s)

As seen in the equations presented above, the slope and length inputs for WEPP are
deeply integrated into the erosion and hydrology components.
 22

2.1.5. WEPP slope gradients and slope length inputs

WEPP uses a slope file to input values for gradient and length. An overall profile length
is specified and slope values are assigned at user specified distances down the profile.
The profile input is processed to describe segments of the profile as simple slope shapes,
namely convex, concave, or uniform using the following equation for a slope at a point:

S*=ax*+b [2.28]

S* = normalized slope (m·m-1)

x* = distance downslope normalized to the slope length (m)
a, b = calculated values describing shape of slope
If a is positive then S is a convex slope
If a is negative then S is a concave slope
If a is zero then S is a uniform slope (b)

To facilitate the calculations of the erosion component in WEPP the values of erosion
have been normalized. Distance downslope (x) is normalized to the slope length (L)
x*=x/L. Slope values are normalized to a uniform average slope passing through the top
and end point on the profile. For calculation purposes, the model assumes a continuous
slope function and computes slope values at individual points.

2.1.6. WEPP slope length and gradient sensitivity

A variety of studies have been conducted to determine the sensitivity of the various input
parameters of WEPP with regard to the sediment delivery, runoff depth, peak runoff,
average soil detachment, average sediment deposition, and total sediment yield. These
include studies by Nearing et al. (1990), Chaves and Nearing (1991), Risse (1994), Zhang
et al. (1996), and others. A few of these have studied the effects of slope gradient, shape,
and length.
 23

In a study by Flanagan and Nearing (1991), the effects of slope steepness, length, and
shape from a single storm and in a continuous simulation were presented for WEPP
model version 91.2. It was found that for the single storm the sensitivity of detachment
and deposition were high for both steepness and length. Slope shape also influenced the
predicted detachment rates. For continuous simulation, runoff was slightly sensitive to
slope steepness and was attributed to effects of depressional storage on the hillslope.
Detachment and sediment delivery were sensitive to average slope steepness. For fall
plow corn conditions, detachment was sensitive to slope length whereas for no-till corn it
was not. This was attributed to limited rill detachment as a consequence of residue cover.

A comparison of WEPP and RUSLE for uniform slope plots with regard to soil loss was
conducted by Nearing et al. (1995). They found that for both models, cropping and
management practices have a strong influence on both slope gradient and length
relationships to soil loss. Since rills tend to facilitate soil loss, if the ratio of rill to
interrill erosion is diminished the soil loss is also diminished. The study showed that
there is a greater sensitivity to slope length for fallow and conventionally tilled soil,
where rilling is more severe, than for no-till and meadow conditions, where rilling is less
severe, for both WEPP and RUSLE.

2.1.7. The WEPP watershed model

The WEPP watershed model is an extension of the hillslope model that is applicable to
small agricultural watersheds. The model was originally designed for use on field-sized
watersheds with a size limit of about 260 ha (Foster and Lane, 1987). The WEPP
watershed model is not applicable to areas containing classical gullies or stream channels
which may have (1) headcut erosion, (2) sloughing of gully sidewalls, (3) seepage effects
on erosion in channels, (4) perennial stream channels, or (5) partial area hydrology
(Ascough et al., 1997).

The three main components of the WEPP watershed model are hillslopes, channels, and
impoundments. Each watershed can therefore be divided into these components (Figure
 24

2.1). Hillslopes can drain to the left, right, or top of individual channels. Each hillslope
can only contribute sediment and runoff to a single channel. All runoff and sediment
yield from adjacent hillslopes is distributed evenly along the channel length. Runoff and
sediment yield from top hillslopes drain directly to the top of the channel. Hillslopes are
represented as rectangular areas with a width and length representing the area of the
hillslope as can be seen in Figure 2.2. The topography of each hillslope is represented by
a single slope profile with length equal to the length of the representative rectangular
hillslope. Hillslopes can further be divided into overland flow elements (OFE's) as long
as they are divided perpendicular to the direction of runoff on the hillslope (Figure 2.1).
An OFE represents a unique area within a hillslope consisting of a single management
practice, crop, and soil type. In many cases, the whole hillslope will have a uniform soil
type, crop, and management practice, which will eliminate the need to divide the hillslope
into OFE's. The WEPP watershed model can also simulate various types of
impoundments. These include farm ponds, culverts, filter fences, and check dams.
Impoundments can be located between channels or at the bottom of a hillslope.
 25

Figure 2.1. WEPP watershed description.

Figure 2.2. WEPP watershed formats.

 26

In a typical small agricultural watershed with no impoundments, the main components

are the hillslopes and channels. The runoff and sediment yield from the hillslopes is then
routed through the channels using an adaptation from the CREAMS model channel
component (Knisel, 1980). The channel erosion routines are similar to the hillslope
model except that flow shear stress is calculated from regression equations developed by
Foster et al. (1980) and only entrainment, transport, and deposition of concentrated flow
are simulated (Ascough et al., 1997).

2.1.8. Erosion models summary

Erosion on hillslope profiles can be divided into two major components, rill erosion and
interrill erosion. Interrill erosion is assumed to be constant down a profile on uniform
slopes, whereas rill erosion may increase down the slope. The sediment load is a
summation of the rill and its constant interrill contribution.

In RUSLE both interrill areas and rills are lumped together, but the effect of rills and
interrill areas are adjusted by changing the equations for slope factor and the ratio of rill
and interrill erosion in the length factor. The slope gradient effect is greater for the long
slopes. Rill erosion is more prominent than interrill erosion in longer slopes than in
shorter slopes.

WEPP partitions runoff between rill and interrill areas. It calculates parameters such as
shear stresses based on rill flow and rill hydraulics rather than sheet flow. WEPP models
the soil loss by adding the sediment loss from the interrill area to the rill erosion, which is
dependent on detachment capacity by rill flow and sediment transport capacity.

WEPP uses an empirical equation to estimate sediment yield from interrill areas. The
individual processes of the detachment and transport over the interrill area are not
present, but the empirical relations have been based on numerous studies and are
considered simple, stable, and to give reasonable results (Nearing and Parker, 1994).
Interrill contribution to rills is considered independent of distance downslope and
 27

sediment delivery to rills is adjusted for ground cover, dead roots, live roots, canopy
cover, and plant and soil management practices affecting infiltration. Interrill processes
may dominate in situations where gradients are low and slope lengths short. Other
examples where interrill processes may dominate are rangelands and no-till areas (areas
in which the formation of rills is limited). For these conditions, interrill cover and
erodibility factors have an important influence on soil loss.

In most other cases rills are the major carriers of sediment. As slope length increases, the
occurrence of rills increases. As slope steepness increases rill erosion also increases. For
these situations, hydraulic factors and factors affecting friction in the rills become very
important.

In WEPP, rills transport the sediment downslope. Any factor that will influence the
formation of rills, quantity of rills, or transport and detachment in a rill will have a
significant effect on the sediment loss. One major influence can be surface cover, which
may have a strong effect on both slope gradient and length relationships to soil loss.
Similarly, management practices and soil factors, such as texture, that affect the hydraulic
properties in the rills and erodibility in the rill will have an influence on the effects of
gradient and slope length on soil loss.

The effect of topographical input such as gradient and slope length on soil loss varies
under different management practices or soil types. This can be illustrated with the help
of Figure 2.3. Assuming all slopes (A,B,C, and D) have the same gradient of 5 percent,
the soil loss ratio of short to long slope length will not be the same under different
management practices or soil types. For this case, the soil loss ratio of A/B was 0.90 and
for C/D it was 0.50 because of the different tillage practices. However, since the model
is very complex, and takes into account a variety of climatic, soil, crop, and management
factors, the exact outcome will differ from case to case.
 28

A C

B D

no-till system conventional

till system

Where,
A, C = 20 meter plot at 5% slope
B, D = 50 meter plot at 5% slope *using Monona soil series and 100 year Cligen generated climate

Figure 2.3. Slope length effect affected by management system.

The simulation and prediction of soil erosion is very complex. The empirically-based
USLE and RUSLE have proven to be great advancements in the science of soil erosion
providing predictions of average soil loss over long periods of time. The physically-
based WEPP model has taken erosion prediction further by providing estimates of
erosion and deposition along a slope profile over a continuous simulation or for an event
storm. Future modeling should include more features and increase reliability of
predictions, but it is important to note that even if we could identify all the functions and
parameters related to soil erosion, simulations will always be subject to uncertainties
because of natural variability in soils, topography, cover, and other factors.

When using the WEPP watershed model, the factors that will influence choosing the
appropriate slope profile are the selection of the hillslopes within the watershed and the
water flow routing within each hillslope. In a 3-dimensional landscape, a hillslope can be
of irregular shape either converging, diverging, or uniformly draining into a channel or
impoundment. For these cases, different slope profiles (with different lengths and
gradients) exist. As in USLE and RUSLE, selection of the appropriate length and
gradient (or slope profile for WEPP) for this situation is a problem since it is subject to
much discrepancy between users. Therefore, creating an algorithm that minimizes the
errors in representing a slope profile for irregularly shaped hillslopes is important.
 29

2.2. Defining hillslopes

In erosion modeling, the main topographical modeling unit is the hillslope. A hillslope is
usually defined as a transect of topography starting at the top of a ridge or a point where
overland flow begins and ending at a point where deposition starts or at the edge of a
channel. Where the hillslope ends will depend on whether the model being used can
simulate deposition or not. The width of the hillslope depends on the uniformity of the
hillslope. Combinations of hillslopes can then been used to define watersheds.

Some work on defining hillslopes has been done for the purpose of using the USLE and
RUSLE models. However, hillslopes for these models are defined as starting from the
point of origin of overland flow to the point where either the slope decreases to the extent
that deposition begins or the runoff water enters a well defined channel (Wischmeier and
Smith, 1978; Renard et al., 1997).

In hydrological modeling, irregular shaped hillslopes are often approximated by a

rectangular plane draining into a channel segment (Wooding, 1965). Such a
representation, called a Wooding plane, uses the width (W), length (L), and slope (S) to
reproduce the hydraulic runoff characteristics of the original hillslope.

The watershed model KINEROS makes use of this concept to simulate runoff (Smith et
al., 1995). The flood hydrograph package HEC-1 also makes use of the concept of
overland flow being characterized by a rectangular plane draining into a channel link to
simulate flooding (U.S. Army Corps of Engineers, 1985).

DEM’s and GIS are ideal for representing hillslopes using the Wooding representation.
Garbretch and Martz (1997), provide algorithms to calculate the length and the slope of
the representative hillslopes. The length (L) can be represented by the following
expression:
 30

∑l i * ai
L= i =1
n
[2.29]
∑a
i =1
i

l = flow path length (m)

a = upstream drainage area (m2)
i = flow path counter
n = total number of flow paths in hillslope

The assumptions of this are:

(1) that a flow path is the route traveled by the water from the upstream source (top of the
hillslope) to the end of the hillslope where it drains into the channel and (2) that flow
paths of large discharge contribute proportionally more to the runoff hydrograph
characteristics of the hillslope than flow paths with small discharge.

The width of the hillslope can easily be obtained by dividing the total area of the irregular
hillslope by the calculated length.

Slope can be calculated in a variety of ways. One approach, also presented by Garbretch
and Martz (1997), is the following:
n

∑ s *k i i
S= i =1
n
[2.30]
∑k
i =1
i

Where, s is the average flowpath slope, k a weighing factor, and n is the total number of
flowpaths in the hillslopes.
The weighing can be done in the following ways:
1. upstream drainage area ki=ai
2. flow path length ki=li
3. upstream drainage area times flow path length ki=ai*li
 31

Another option is using one of the many algorithms available to determine a slope value
for each raster cell in an irregular hillslope. A mean value of all of these could also give
a representative slope and a good idea of distribution of slope within the irregular
hillslope. One problem with this approach is that the slope values for each cell are
dependent on the chosen algorithm and can differ significantly (Srinivasan and Engel,
1991).

Four widely used algorithms to predict slope for individual raster cells from a digital
elevation model were tested by Srinivasan and Engel (1991). Two 10 by 10 cell data sets
from the DEM were selected and slopes were calculated using the neighborhood method,
the quadratic surface method, the best fit plane method, and the maximum slope method.
The neighborhood method uses the eight neighboring cells to predict slope for the center
cell, however, it does not take into account the elevation of the center cell. This can
cause problems if the elevation data has small sinkholes (pits) or small ridges.
Smoothing of the elevation data has to be conducted before using this method to
eliminate these pits and ridges. The quadratic surface method (Zevenbergen and Thorne,
1987) uses a partial quadratic equation to fit the surface that passes through nine
elevation data points. The slope is calculated by the derivative of this polynomial which
takes into account four neighbors, but not the center cell elevation. Pits and ridges also
pose a problem to this method and have to be eliminated prior to the calculation of slope
values. The best fit plane method used in ANSWERS (Beasley and Huggins, 1982) uses
a linear equation to fit a plane through the four corner points of a cell. Slope is calculated
from the linear fit of the four diagonal neighbors. Slope values for individual cells using
the maximum slope method are calculated by selecting the maximum slope between the
eight neighbors using the distance between the midpoints of the cells and the center cell
as the reference slope distance. This method is normally used to calculate slope using
contour maps, but it tends to over-predict slope values when using square grid cells.

Results from the application of these four methods were compared to an observed slope
estimation derived from USGS topographic quad-sheets and site observations. Results
 32

showed that the neighborhood method most closely approximated the observed slope
values. The maximum slope method always over-predicts the slope. Compared to the
neighborhood method, the maximum slope method over-predicts by 286% for flat areas
and 80% for steep areas. This suggests that more importance be given to the methods of
estimating slope and other topographic parameters, particularly for runoff and erosion
models where topography plays a major role in the simulation of these processes
(Srinivasan and Engel, 1991).

2.2.1. Irregular topography in erosion modeling

The next step in defining hillslopes is being able to create them using non-uniform
topography. Foster and Wischmeier (1974) proposed an equation to calculate the LS
factor for irregular slopes. In this equation, irregular slopes are divided up into segments
of uniform slope. The LS factors from these segments can then be combined to create an
effective LS value for a complex profile.

∑S
j =1
j (λ j
m +1
− λ j −1
m +1
)
LS e = [2.31]
(λ j − λ j −1 ) ⋅ (22.13) m

where
LSe = effective LS factor
Sj = slope factor for j-th segment
λj = distance from top of slope to lower end of segment (m)
λj-1 = slope length above segment (m)
m = USLE coefficient related to rill/interrill erosion

This procedure was incorporated by Desmet and Govers (1996a) in a computer algorithm
to calculate USLE and RUSLE LS factors over a two-dimensional landscape using GIS,
but physically-separable hillslopes were not defined. Instead, individual grid cells, where
deposition was not expected to occur, had LS factors assigned to them. This method
 33

provided continuity of hydrological properties throughout the watershed and effects of

flow convergence were accounted for to some degree. This study found that the
automation of the calculation of the LS factors using GIS eliminates the subjectivity that
can arise in the selection of length and slope parameters by different users, but that the
precision of the final simulation is determined by the accuracy of the DEM. A similar
approach was used by Srinivasan and Engel (1991) to calculate the LS factors for each
grid cell in a GIS simulation.

Moore and Burch (1986a) developed a physically-based analytical framework to predict

the effects of topography on erosion and deposition on two dimensional and quasi three-
dimensional non-planar hillslopes. They started by calculating sediment transport
capacity by using Yang’s (1972, 1973) and Yang and Song’s (1979) equation for an
alluvial channel:

C = γ ( P − PC ) β [2.32]

C = sediment transport capacity of an alluvial channel (g·m-3)

P = unit stream power (time rate of potential energy dissipation per unit weight of water)
PC = critical unit stream power at incipient sediment motion
γ and β = functions of the median sediment size, the kinematic viscosity of water, and the
terminal fall velocity of sediment particles in water

Although other sediment transport equations could also be used, this equation was
selected because it is simple, conceptually attractive, and has already been found to be
applicable to sediment transport in shallow sheet and rill flow (Moore and Burch,
1986a,b and Alonso et al. 1981). When applying this equation to overland flow, a
simplification of the equation is done by neglecting critical unit stream power:

C = γ (P) β [2.33]
 34

This is an acceptable assumption for many situations. Unit stream power can then be
written in the following form (Moore and Burch, 1986a,b):

q x0.4 ∂y 1.3
P= ( ) [2.34]
n 0.6 ∂x

Where qx is discharge per unit width at time t and at point x on a slope and n is
Manning’s surface roughness coefficient.

The sediment flux per unit width is:

Yb = ρq x C [2.35]

Where,
Yb = sediment flux per unit width (kg·m-1·s-1)
ρ = density of water (g·m3)

The erosion/deposition rate (kg·m-2·s-1) is:

Yr = ∂Yb / ∂x [2.36]

If the effects of cover and soil detachment are ignored and only the topography is
accounted for, the following equations apply for a hillslope profile:

1.3 t
 ∂y 
Y∝x  
1.4
 ∂x  ∫i
0
1.4
dt [2.37]

Where Y is the cumulative sediment yield per unit width (kg·m-1) and i is the rainfall
excess.
 35

1.3
 ∂y 
Yb ∝ x   i 1.4
1.4
[2.38]
 ∂x 

 ∂y ∂ 2 y  1.4
0.3
 ∂y 
Yr ∝ x  
0.4
 14
. + 13
. x 2 i [2.39]
 ∂x   ∂x ∂x 

Where Yb is the sediment flux per unit width (kg·m-1·s-1) and Yr is the erosion and
deposition rate (kg·m-2·s-1).

Experiments were conducted on 3 sinusoidal slopes and 1 uniform slope using these
equations. Results showed that hillslopes having convex profiles near the end of their
outlets had a greater sediment flux.

The same equations were expanded to deal with a three dimensional case. For example,
the cumulative sediment yield per unit width (kg·m-1) was computed as:

1.4 1.3 t
 x (1 + xZ / 2)   ∂y 
Y ∝   
 (1 + xZ )   ∂x  ∫i
0
1.4
dt [2.40]

Where Z = θ/Wo and Wo is the top width of the hillslope and θ is the convergence or
divergence angle in radians. Theta is positive for diverging and negative for converging
surfaces. A similar approach is used to estimate the sediment flux per unit width and the
erosion/deposition rate.

The application of these equations to a non-convergent slope, a 20 degree divergent

slope, and a 20 degree convergent slope showed that convergence can have a major
impact on erosion whereas divergence has a lesser impact. Convergence concentrates the
flow, increasing flow velocity and therefore increasing the sediment transport capacity of
the runoff (Moore and Burch, 1986b).
 36

2.3. Watershed modeling with DEMs

There are three ways by which data can be manipulated using digital elevation models.

1. TINs (triangular irregular networks)

2. Vector or contour lines (vector based GIS)
3. Grids (raster based grid cells)

Triangular irregular networks are used primarily to represent terrain surfaces or

landscapes. TINs have also been used by some to model hydrological processes in small
watersheds. Vieux (1991), for example, used the TIN facets to provide land surface slope
in a finite element solution of overland flow. Other studies have also used TINs to apply
finite element methods to solve the kinematic wave equations for runoff. However, no
studies have been conducted to determine well defined hillslopes using TINs.

Vector based DEMs are based on the familiar contour topographic maps. Information on
elevation is stored as contour lines each having a series of x and y coordinates.
Analyzing this type of representation of elevation is more computationally demanding
that using a standard grid based system, but one of the advantages of this system is that
better accuracy of topography can be maintained at larger scales than in grid based
systems for very complex and dissected terrain. Vector based DEMs can be used to
model hydrology and erosion within watersheds. One model that uses the vector
approach to modeling is TAPES-C, Topographic Analysis Programs for the
Environmental Science – Contour (Moore and Grayson, 1991). This program
automatically partitions a catchment into interconnected, irregularly shaped elements
using a “stream tube” approach. The “stream tube” concept consists of assuming that
contour lines are equipotential lines and pairs of orthogonals to the equipotential lines
form the “stream tubes.” The distance between the pairs of orthogonals to the
equipotential lines can be arbitrarily defined, but they start at a specific section of a
channel and end at the edge of the catchment. The result is a flowpath with a certain
thickness. These flowpaths could be grouped together to create hillslopes. Modeling
continuous hydrology or erosion using vector based hillslopes is no trivial task and
 37

development of vector based applications for this purpose has been limited as compared
to the more popular raster based format.

For the present study we have chosen to use the grid approach. The following are the
main advantages of a raster based system (Burrough and McDonnell, 1998):

1. Simple data structures (each grid cell is represented by a single data value).
2. Location-specific manipulation of attribute data is easy (each grid cell has a X and Y
coordinate that represents its geographical location).
3. Many kinds of spatial analysis and filtering may be used.
4. Mathematical modeling is easy because all spatial entities have a simple, regular
square shape.
5. The technology is inexpensive.
6. Many forms of data are available and commonly used.

Some of the disadvantages of the raster based system have been overcome by advances in
computing storage capabilities and data processing speeds of today's computers. For
example, one of the disadvantages of the raster based systems in the past was large data
volumes, which now is not a big problem because storage capacity is no longer a major
concern.

Hillslopes can be create using raster based DEMs by first extracting a drainage channel
network from a defined watershed. For each channel link, the hillslopes can be defined
as the areas that drain into the left, right, and top of the channel (if a first order channel).
Hillslopes defined in this manner maintain their overland flow characteristics and include
areas where deposition might occur. The size, shape, and number of these hillslopes are
dependent on the density of the channel network extracted by flow routing algorithms.
Flow routing algorithms are therefore explained in the next section.
 38

2.3.1. Flow routing algorithms

There have been many proposed algorithms for flow routing and to identify channel
networks within a watershed by using raster based DEMs. These can be divided into
single flow and multiple flow algorithms. Single flow algorithms transfer all the flow
from the source cell to a single neighboring cell downslope. Multiple flow algorithms
divide the flow from the source cell between many neighboring cells. Single flow
algorithms can only simulate parallel and convergent flow, where as multiple flow
methods can also simulate divergent flows.

Single flow:

Flow routing to identify the channel network is usually obtained by using an algorithm
that determines the steepest descent direction between cells. Examples of these
algorithms are present in the work of O’Callaghan and Mark (1984), Jenson and
Domingue (1988), Martz and Garbretch (1992), and Zenenbergen and Thorne (1987).
Some of these algorithms have been adopted by popular geographic information systems
such as GRASS, Arc ViewTM (Spatial AnalystTM), and IdrisiTM.

The steepest slope criterion can be illustrated in the following way:

z i, j − z i +a, j +b
S i , j = max( ) [2.41]
df a ,b

Where,
Si,j = local slope at cell(i,j)
d = the cell size
a,b = coefficients to define the neighboring cells (-1,0,1)
fa,b = a coefficient related to the cardinal or diagonal direction of neighboring cells
(sqrt(2) if a and b ≠0, 1 otherwise)

The flow direction is thus assumed to be from the cell(i,j) to the cell(i+a,j+b) for which
local slope is the maximum.
 39

Other less popular algorithms rely on aspect direction for flow routing (eg. Gardner et al.
1990; Drayton et al. 1992; and Lea 1992).

Multiple flow:

Two of the most talked about multiple flow algorithms in the literature are the ones used
in the ANSWERS model (Beasley et al. 1980) and an algorithm developed by Quinn et
al. (1991). The ANSWERS algorithms divide a grid cell into two parts using a line
passing through the cell center in the direction of the aspect of the cell. Each of the parts
then contributes its flow to the corresponding downslope neighbor. The algorithm
proposed by Quinn et al. (1991) uses the following equation to divide the flow passed to
the downslope cells:

tan
∆Ai
i
A n

∑ (tan
j
j Lj
 40

algorithms are believed to be superior in zones of convergent flow and along well defined
valleys (Freeman, 1991; Quinn et al. 1991). However, the accuracy of the DEM and its
grid size may be more influential than the flow detection method. (Pilotti et al. 1996).
The calculation of topographic attributes and all applications in which a routing algorithm
is present, are greatly affected by the spatial scale of the elevation data. Increasing the
element size of the DEM implies a loss of topographic information, especially in regions
of pronounced curvature. This will lead to a modification of flow distribution. (Desmet
and Govers, 1996b).

2.3.2. Removal of depressions and flow routing through flat areas

One important requirement to run the flow algorithms is the removal of all small
depressions (pits or sink holes). Small depressions prevent a continuous flow routing to
be determined. The most common algorithm to accomplish this is presented by Jenson
and Domingue (1988) where the depressions are first identified and then eliminated by
raising the depression cells to the height of the lowest pour point. Garbretch and Martz
(1997) have made significant improvements to this method by creating an algorithm that
is able to distinguish between two types of depressions, sink depressions and
impoundment depressions. Sink depressions are caused by a group of cells with a lower
elevation than the surrounding cells, whereas, impoundment depressions are caused by a
narrow strip of cells across a drainage path. Sink depressions are treated in the same
manner as described before, but impoundment depressions are treated by simulating a
breach in the impoundment.

Flow routing through flat areas can also be a problem. Flat areas represented by clusters
of cells of equal elevation can be real topographic features, the result of insufficient detail
in the data, or a result of filling depressions. Whichever the case, flow routing through
these areas has to be done by assignment of flow directions based on the surrounding
topography. Early algorithms tended to create unrealistic parallel flow lines through flat
areas. Tribe (1992) suggested an algorithm that creates converging flow patterns on flat
areas. An inlet and an outlet cell are identified based on the surrounding topography and
 41

a flow line is interpolated between the two cells. Other flow paths converge to this flow
line eliminating parallel flow patterns on flat areas. Algorithms developed for TOPAZ
(Garbretch and Martz, 1997) to deal with flat areas also utilize a similar approach of
converging flow paths in flat areas. The treatment of these flat areas is of particular
importance when applying models that require continuous flow down the hillslope such
as the erosion models presented earlier.

2.3.3. Channel initialization

One of the primary assumptions in watershed erosion modeling is that there is a

difference in processes between hillslope erosion and channel erosion and that by
identifying the channels one can make that distinction. In larger watersheds, channel
erosion and sediment transport in channels rather than the hillslope processes may be
responsible for the bulk of the erosion from the watershed. In smaller watersheds,
hillslope processes (rill and interill) are the main source of erosion. The problem lies in
the fact that it is not exactly clear where the transition between hillslope to channel
processes occur, or how to identify precisely a channel network or drainage density.

Tarbonton et al. (1991) has suggested a way to select an appropriate drainage density
based on traditional geographic scaling laws. It is therefore important to first review the
traditional laws regarding channel networks as characterized by Horton and Strahler
(Horton 1932,1945 and Strahler 1952,1957).

Horton’s laws:
Law of stream numbers Nw-1/Nw = RB (bifurcation ratio)
Law of stream lengths Lw/Lw-1 = RL (length ratio)
Law of stream areas Aw/Aw-1 = RA (area ratio)
Law of stream slopes Sw-1/Sw = RS (slope ratio)

Nw is number of streams of order w

Lw is mean length of streams of order w
A is mean area contributing to streams of order w
w is the mean slope of streams of order w

R , RL A, and R represent characteristics of a particular network.

Ordering (w) of channels is important. The most widely used ordering system is the

1. All exterior links have order 1.

When two upstream links of the same order join, the downstream link has its order

3. When two upstream links of different order join, the downstream link takes the higher

d = L /A

D = drainage density
T = total length of channels
A = contributing area

Constant drop law:

rop law refers to the change in slope with channel order. It was noted that
the average drop (H which on average is equal to Sw·Lw) along streams of order w was
approximately constant, in other words, independent of order (Broscoe, 1959). This
implies that typically the slope ratio is close to the length ratio:

Sw-1/Sw ≈ Lw/Lw-1 or that RS ≈ RL [2.44]

 43

Power law scaling of slope with area:

The power law scaling of slope with area states that slope is empirically related to
contributing area (Tarbonton et al., 1991):

S ∝ A-θ [2.45]

where
ln RS
θ= [2.46]
ln R A

In extracting channel networks from digital elevation models it is important that the
networks extracted be close to what traditional workers using maps or fieldwork would
regard as channel networks. (Tarboton et al., 1991) It was also suggested that the channel
networks should have as fine a resolution as possible to the point where the key elevation
scaling laws break. The procedure they propose is to look for the smallest scale, in terms
of the threshold area, where the constant drop law and the power law of scaling of slope
with area break. It is argued that the breaking point represents a physical transition from
channel erosion and transport mechanisms at large scale to hillslope erosion and transport
mechanisms at small scale (Tarboton et al., 1991).

Fractal geometry, introduced by Mendelbrot (1983), can also be used to study channel
networks and define channel initiation. Fractals are used to define the length to area
relationships between channel networks and basin areas. Veltri et al. (1996) and Puente
and Castillo (1996) conducted such studies on the fractal nature of natural channels.
These studies showed that the laws pertaining to channel networks could be applied at all
scales. Moussa and Bocquillon (1996) later defined “the scale of observation” as S/So (S
is the threshold area for channel initiation and So is the catchment area). Moussa and
Bocquillon (1996) studied the effects of changing the threshold area and calculated three
 44

constants that are independent of the observation scale, the bifurcation ratio, the channel
network shape index, and the elongation index. It was suggested that these could be used
to define the scale of observation and therefore the threshold area for the channel
initialization.

Some studies suggest a more physically-based approach to solving the problem of the
initialization of a channel. Once depressions in a DEM are removed and the routing
algorithms are applied, the channel network density can be determined by specifying a
threshold value for channel initiation (critical source area or a threshold area). The
threshold area to initialize a channel is usually determined arbitrarily even though it is
known that by choosing different thresholds, the density of channel network created will
vary significantly.

However, a few studies have been conducted that allow for the calculation of a critical
threshold area. For example, exceedance of a critical threshold area at cell (i, j) can be
formulated in the following way (Pilotti et al., 1996):

Ai,j ≥ (As)i,j = Ci,jSi,j-α [2.47]

where
A = contributing area at cell (i,j)
As = critical threshold area
α = an exponent whose value depends on flow conditions
C = coefficient that depends on both climatic and soil properties
S = slope at cell (i,j)

Most studies and algorithms use C = As and α=0. This means that a constant threshold
area is chosen for the whole watershed, without any regard to possible spatial variability
of the threshold area within the watershed. The threshold value of the upstream area
depends, among other things, on slope characteristics, soil properties, surface cover, and
 45

climatic conditions (Martz and Garbrecht, 1992). Algorithms in TOPAZ, for example,
allow the user to introduce this variability into the computations. This is done by adding
a second map layer with weighted values for the threshold areas to the original elevation
map layer.

Some scientists believe that channel initiation is controlled by four processes: incision by
Hortonian overland flow, incision by saturated overland flow, seepage erosion, and
shallow landsliding (Dunne, 1980; Montgomery and Dietrich, 1989). It is also believed
that threshold area (As) to define the initialization point of a channel can be expressed as
a function of the ground slope and of some soil and climatic properties (Montgomery and
Dietrich, 1989). Montgomery and Dietrich (1994) propose the following equations to
calculate the threshold area needed to start a channel under laminar and turbulent flow
conditions.
Laminar flow:
τ c5 / 3 1
as = [2.48]
pn γ 5/3
S 7/6

Turbulent flow:
2τ c3 g 1
as = [2.49]
pkvγ 3 S 2

where
S = slope (m·m-1)
as = support area per unit contour length (m2·m-1)
p = uniform rainfall rate in excess of infiltration capacity (m·s-1)
k = dimensionless coefficient dependent on geometry of the erodible bed
v = kinematic viscosity of the fluid (m2·s-1)
γ = sediment-laden water specific gravity (N·m-3)
g = gravitational acceleration constant (m·s-2)
n = Manning roughness coefficient (m-1/3·s-1)
 46

τc = critical shear stress (N·m-2)

Prediction of channel initiation is not a simple task. All of the methods presented here
have had some degree of success, but in many cases comparisons of predicted and
observed channel initiation show large errors. Some of these errors can be attributed to
low accuracy DEMs, while others to natural variability in local slope, soil parameters,
climate, and cropping factors. Channel heads may also advance and retreat depending on
discharge cycles (Montgomery and Dietrich, 1994). This implies that the actual location
of initiation of a channel can vary over time too.

The identification of the channels within a watershed is a crucial step in watershed

modeling. Integration of channel networks with GIS is the next logical step as many of
these procedures already rely on DEMs to extract channels. The logic and levels for the
integration of watershed modeling with GIS is described next.

2.3.4. Watershed models using GIS

There are at least two purposes for developing an environmental model. The first one is
to assist in the understanding of a physical system by having a framework for hypothesis
testing (Grayson et al., 1992). The second one is to provide a predictive tool (Beven,
1989). The integration of environmental models with GIS fits both of these purposes.
GIS can be used as a tool to help our understanding of spatially distributed watershed
modeling. It can be used to test different hypothesis on issues such as data resolution,
effects of topography, and scale. Integration of GIS and environmental models also help
facilitate the implementation of the watershed models as predictive tools.

There are three levels of integration between a GIS and an environmental model (Savabi
et al., 1995; Tim and Jolly, 1994). The first level of integration is an "ad hoc integration"
where the model and the GIS are developed separately and the GIS is only used as a
means to create input data for the model. This type of integration was used by Savabi et
al. (1995) to produce soil files and some of the slope files for the WEPP model. Similar
 47

systems have been used to create input files for the USLE and RUSLE models (Huang
and Ferng, 1990; Flacke et al., 1990). The second level of integration is a "partial-
integration." A partial integration consists of an interface that accesses GIS maps and
databases to parameterize the model. The development of an interface to exchange data
between the GIS and the model is a crucial part of the partial integration system. Most of
the current GIS and environmental model integrations fall into this category. A few
examples of models that have been integrated with GIS are SWAT (Manguerra and
Engel, 1998), ANSWERS (Rewerts and Engel, 1991), AGNPS (Line et al., 1997) and
others. The third level of integration is a "complete integration". In this level, the
environmental model is completely integrated into the GIS. Due to the structure and
complexity of the current WEPP watershed model, a full integration would not be
possible at this time without completely modifying the original code. However, a partial
integration is possible without modifying the WEPP source code. Since it is not the
intention to modify the source of WEPP in this thesis, a partial integration system will be
developed to test methods of applying WEPP with GIS and DEMs.

2.4. Literature review summary

A general description of the USLE, RUSLE, and WEPP erosion models as well as a
comprehensive look into the details and challenges of watershed modeling with GIS and
DEMs were presented in this chapter. The review of the three erosion models described
our current understanding of erosion processes and the need for further knowledge. Of
particular interest for the research conducted in this thesis was the description of the
effects of slope and slope length on the soil loss from hillslopes for each of the models.
The understanding of these effects is fundamental for the development of algorithms to
integrate erosion models and GIS. Furthermore, the advancement in the integration of
these models with GIS technology using DEMs is an obvious step in increasing our
knowledge of watershed modeling.
 48

When using DEMs to model watersheds, a few general steps have to be followed. First,
the DEM has to be treated to remove all inconsistencies such as sink holes or depressions.
Once this is done, the desired watershed can be delineated. Subsequently a channel
network has to be extracted by using flow routing algorithms. Special care has to be
taken so that the density of the simulated channel network matches the real channels in
the watershed. After the channels are identified, hillslopes draining to the channels can
be extracted. Finally, gradients and slope lengths for each of the hillslopes have to be
defined. This process can be used to interface the WEPP watershed model to a GIS.
However, apart from the variety of interfacing problems that have to be solved, there are
also many theoretical questions that have to be answered. The main question is whether
channels, hillslopes, gradients, and slope lengths can accurately be extracted from DEMs
for the WEPP input and if the simulations conducted with them are comparable to a
manual application of WEPP. There is also the question of DEM resolutions and how
accurate the data has to be to run the model and obtain acceptable results. These
problems and questions will be addressed and answered in the following chapters.