1. Introduction.............................................3 2. Methods..................................................4 2.1 Input data processing................................5 2.2 Erosion models.......................................8 2.2.1 USLE...........................................9 2.2.2 USPED.........................................10 2.2.3 Process based simulation of water erosion.....11 2.3 GIS implementation of models........................14 2.3.1 GRASS5.0......................................14 2.3.2 ArcView Spatial Analyst.......................16 3. Applications............................................18 3.1 Fort Hood subwatershed: comparison with Cs137.......18 3.2 Fort Polk erosion potential....................21 4. Conclusion and future directions........................22 5. References..............................................23 6. Appendix................................................24 6.1 Manuals and tutorials 6.2 Data 6.3 On-line documents

Research and experience have delivered a rich set of tools and techniques for evaluating and managing soil erosion from various sources within watersheds. While significant gains in the control of soil erosion have been made, there remain significant problems with sedimentation of water supply and navigation systems, transport of soil attached contaminants and the cost of some management techniques. Further, while much research exists for individual methods in individual locations little work has been devoted to integrated approaches which characterize interaction between erosion processes and management methods across scales and within a system of land management strategies

To better support the multilevel management we propose a methodology for watershed characterization and erosion modeling at multiple scales and levels of complexity. Because the spatial unit, modeled at the local level is a part of a larger watershed, the evaluation of the impact of numerous, locally implemented conservation practices on the entire watershed requires multiscale approach which links the high resolution, local level simulation with low resolution/regional simulation. Models and simulations of watershed behavior demonstrating how the lands within the watershed are connected, how they interact and influence each other help to better understand the impact of local actions within the watershed necessary for successful conservation and sustainable land management programs.

Recent advances in Geographic Information Systems (GIS) technology, linkage of numerous models with GIS, (Moore et al. 1993, Saghafian et al. 1995, Srinivasan and Arnold 1995) create a potential to develop an environment for coordination of conservation efforts at a hierarchical system of management levels by providing the tools to design and evaluate the impact of land use alternatives at both local and watershed/regional levels for a given time horizon. They facilitate evaluation of prevention practices not just depending on the type of the prevention measure, but also on their location within the watershed. Spatial analysis and simulations can also provide supporting information for allocation of resources to those areas and those types of practices which will provide the most effective protection. Distributed models facilitating multiscale approach have been recently identified as promising for studies of landscapes with significant human impact, such as agricultural and urbanized watersheds. However advances in models, algorithms and GIS tools are needed to fully support this approach.

Simulation modeling which can give decision makers interactive tools
for both understanding the physical system and judging how management actions
might affect that system was identified in a report "New Strategies for
America's Watersheds" (Committee on Watershed Management, National Research
Council 1999) as one area of special promise for watershed management.
This report analyzes the current status in watershed modeling for decision
making and concludes that the available models and methods are outdated
and "a major modeling effort is needed to develop and implement state-of-the-art
models for watershed evaluations"(pp.160-161). The development of distributed
models supported by SERDP and this project is one step in the direction
identified in the report as crucial for the land use management in America's
Watersheds.

*Local level* of management requires detailed spatial representation
(10-1m resolution) and models capable to simulate effects of spatially
variable land cover. This level of detail is necessary for implementation
of conservation practices at the most effective locations, as well as for
selecting the projects that save the most soil or benefit the most acres
per dollar cost. The empirical models such as USLE/RUSLE, which have been
used at this level for many years, are now being replaced by process-based
models, such as WEPP (Flanagan and Nearing 1995). Both RUSLE and WEPP are
based on routing of water and sediment over hillslope segments and pre-defined
channels, a concept which requires substantial manual input and is not
very practical for large areas typical for military installations.

New generation of models based on raster representation of multivariate functions has been developed recently with the support from SERDP and ARO (SIMWE: Mitas and Mitasova 1998, CASC2d: Julien, Ogden, Saghafian, Johnson 1995-99, MIT models: Brass, Willgoose, Moglen, SIBERIA, CHILD 1992-1998)). These models have a potential, after undergoing a substantial experimental testing and calibration, to be useful for both local and a regional scale modeling.

The watershed models based on homogeneous spatial units have been described in detail in literature (e.g., Arnold et al. 1993) and are widely available, including the possibility to run the model online. Therefore, we focus on distributed modeling of erosion and deposition patterns for spatially variable conditions. Within the spatially continuous approach, model inputs and outputs are represented by multivariate functions discretized as grids. Flows of water and sediment are described as bivariate vector fields rather than commonly used systems of 1D flows (e.g., Moore et al., 1993). To support modeling at different levels of complexity we have developed a set of tools which range from modifications of relatively simple empirical models to a more complex, process-based approach. GIS technology is used to support the processing, analysis and visualization of the data and simulation results (Mitas et al., 1997). The models are described in detail in our previous report (Mitasova et al., 1998) therefore in the following sections we focus on the issues related to practical applications of the models and their use with GIS.

Preparation of data for modeling and land use management still requires a substantial effort, in spite of the great progress in the measurement techniques and distribution of data. While the spatial data are currently much easier to exchange and import/export between different systems, with higher resolutions the data sets are substantially larger and often more complex. The growth of the size of data sets is "keeping up or even exceeding" the growth of computational power and efficient GIS tools for processing, analyzing and visualizing spatial data for land use management are as important as ever. When processing the spatial data, we are focusing on extensive use of standard GIS tools available in both commercial and public license systems, such as GRASS5.0 and ArcView/Spatial Analyst, however enhancements or modifications of existing tools are often needed to fulfill the needs of distributed, process-based modeling.

**Spatial interpolation** and smoothing is necessary for transformation
of scattered data to regular grids, resampling of grid data and smoothing
of data with noise. The selection of an adequate method and its parameters
can have a profound impact on the predicted patterns, as we have shown
e.g. in Mitas and Mitasova 1999. The current project has again demonstrated
the crucial role of spatial interpolation for processing of DEMs from various
sources (ISARE, contours, field measurements) as well as for the interpolation
of Cesium data. We have used and enhanced the bivariate regularized spline
with tension and smoothing (RST, Mitas and Mitasova 1999) and its implementation
in GRASS5.0 to perform spatial interpolation and topographic analysis for
this project applications, in particular for smoothing of the IFSARE data
(see section 3.1), interpolating field measured elevations with fault,
interpolating very large combined point and contour data set and deriving
surface gradients for water and erosion models.

**Interpolation of data with a fault by RST.** Spline is by definition
a smooth function and has always been regarded as a tool for representation
of smooth surfaces. Therefore its application to surfaces with faults has
been often described as problematic. However, RST combined with GIS masking
tools can be very effectively used to create models of surfaces with complex
pre-defined faults by using multiple-pass approach similar to zonal kriging.
First, the points defining the fault are extracted from the elevation data
set based on their attributes and transformed from site through vector
to raster masks representing continuous (smooth) subsets of the data. The
masks are then used to create separate elevation data subsets for each
smooth subregion which is then interpolated by RST. The masked surfaces
are then merged using map algebra to create a single surface with faults
(Figure 4a)

**Surface analysis** provides basic topographic or surface geometry
parameters needed for modeling of fluxes in landscape, including water
and sediment. Estimation of slope and aspect (direction of flow) are standard
functions in GIS software and their computation using the RST methods is
described by Mitasova and Hofierka 1993. For several hydrologic and erosion
models which simulate water and sediment flow as a bivariate function,
partial derivatives *dz/dx, dz/dy*
representing the surface gradient grad *z*, are needed:

We have used two approaches to estimate partial derivatives. First is
based on a local polynomial approximation of a surface at a grid point
from the elevations in its 3x3 neighborhood (Shapiro and Westervelt 1992).
This approach can be applied only to raster data and has been implemented
as an option for the command r.slope.aspect in GRASS5.0. It can also be
implemented using map algebra functions as shown e.g. by Shapiro and Westervelt
1992. The second approach is based on the derivatives of the RST function
(Mitasova and Hofierka 1993) and it can be used for both scattered points
and raster data. This approach is implemented as an option for the s.surf.rst
command in GRASS 5.0. To compute the partial derivatives within a GIS where
the direct computation of derivatives is not available we use the known
relationship between the slope *b* and aspect *a* angles
and the partial derivatives of the function *z*:

where *b *is slope and *a* is aspect in degrees, computed
by standard GIS functions.

Flowtracing for computation of slope length and upslope area is performed by r.flow. The program has undergone some modifications in 1995 to reduce the stripes caused by 1D flow on convex hillslopes where dispersal flow occurs. This modification introduces a small level of noise into the results which can be smoothed out by r.neighbors. Better solution is possible by using multiple direction flow. The problems with dispersal flow have been fully resolved by simulation of flow as a bivariate vector field in SIMWE, however a simple flowtracing implemented as a GIS command is still a useful alternative to full process based simulation.

**Terrain modification** tools are useful for incorporation of data
about streams and rivers into terrain model, which usually captures only
the approximate water surface elevation. We have developed a methodology
for modification of a DEM by "carving-in" streams and rivers using their
vector representation and attributes such as width and depth. The method
was implemented as r.enforce in GRASS5.0, an example of terrain model with
"carved-in" river is in Figure 5.

**Visualization** is an integral part of input data processing as
well as of communication of the results. After porting the underlying
libraries to OPEN GL at GMSLab (Brown 1998) and as a result of collaboration
with a research group in Europe, the 3D surface visualization tool NVIZ2.2
is now available for PC under LINUX and for SUN workstations within the
new GRASS5.0 release. The capability to visualize multiple surfaces
has been invaluable in the development of erosion models.

New visualization tools were developed for tracking the computation progress during interpolation of very large DEMs (millions of cells from hundreds of thousands of points) using the RST method, which often takes several hours even on fast, state of the art computers.

We have renewed our use of **p.vrml** which allows us to directly
serve the results as 3D interactive models on the Internet. While most
of the VRML viewers are not well designed for exploration of geospatial
data, their speed and convenience has improved sufficiently for practical
use.

The Universal Soil Loss Equation (USLE) is a well known empirical equation developed for simple, field-based estimates of average annual soil loss using factors derived from a very large experimental data set:

where *E [ton/(acre.year)] *is the average soil loss,* R [hundreds
of ft.tonsf.in/acre.hr.year]* (in SI: *[MJ.mm/ha.hr.year ]*,
*R[SI]=17.02R[EU]* or [N/h=kg m/s^{3}])
is the rainfall intensity factor, *K [tons per acre per unit R] = [tons.acre.hr/hundreds.acre.ft.tonsf.in]*
is the soil factor, *LS [dimensionless]* is the topographic (length-slope)
factor, *C [dimensionless] *is the cover factor and *P[dimensionless]*
is the prevention practices factor. Several approaches were developed to
adapt this field based model to GIS use. *R,K,C* and *P* factors
are usually derived from land use and soil maps (digitized or derived from
imagery), and slope and slope length for the LS-factor are derived from
a DEM. While GIS has relatively good tools for estimating slopes, deriving
a slope-length map is more problematic. Examples of L-factor
estimation and implementation include use of (i) average slope length for
the entire area, (ii) average slope length for each soil type, (iii) spatially
variable slope length estimated for each cell from a DEM using a flowtracing
program such as r.flow in GRASS (Mitasova and Hofierka 1993, Figure 9,
Appendix). The LS factor is then computed as

where *l[m] *is the slope length and *b[deg]* is the slope
angle.

To incorporate the impact of flow convergence, a replacement of the slope length by the upslope contributing area per unit contour width was suggested, e.g., by Moore and Wilson 1992, Mitasova et al. 1996, Desmet and Govers 1996. The modified LS factor at a point on a hillslope is

where *A[m]* is upslope contributing area per unit width, *b[deg]*
is the slope angle, *m* and *n* are parameters. Exponents
*m* and *n* should be calibrated for a specific prevailing type
of flow and soil conditions, if the data are available. Moore and Wilson
(1992) have shown that the values of *m=0.6,* *n=1.3* give results
consistent with RUSLE LS factor for the slope lengths <100m and the
slope angles <14 deg for slopes with negligible tangential curvature.
We use the equation (6) as an approximate *LS* factor for estimation
of soil loss using RUSLE, with the assumption that transport capacity exceeds
detachment capacity everywhere, erosion and sediment transport is detachment
capacity limited and therefore no deposition occurs. Impact of replacing
the slope length by upslope area is illustrated in Figure 9, which shows
that the upslope area better describes the increased erosion in the areas
of concentrated flow and allows to identify areas with potential for gully
erosion. More research is needed to quantitatively determine the erosion
rates in these areas as the type of flow changes dramatically and different
equations or at least parameters are needed to estimate the erosion rates
by concentrated flow.

Both the standard and modified equations can be properly applied only
to areas experiencing net erosion. Deposition areas should be excluded
from the study area because the model assumes that transport capacity
exceeds detachment capacity everywhere and erosion and sediment transport
is detachment capacity limited. Therefore, direct application of USLE/RUSLE
to complex terrain within GIS is rather restricted. The results can also
be interpreted as an extreme case with maximum spatial extent of erosion
possible.

**2.2.2 Unit Stream Power based Erosion/Deposition model**

USPED (Unit Stream Power - based Erosion Deposition) is a simple model
which predicts the spatial distribution of erosion and deposition rates
for a steady state overland flow with uniform rainfall excess conditions
for transport capacity limited case of erosion process. The model is based
on the theory originally outlined by Moore and Burch 1986 with numerous
improvements. For this case, we assume that the sediment flow rate **q***s[kg/ms]?*
is at the sediment transport capacity *T[kg/ms]*, (Julien and Simons
1985)

where **q***[m ^{2}/s]* is the water flow,

No experimental work was performed to develop parameters needed for USPED, therefore we use the USLE or RUSLE parameters to incorporate the impact of soil and cover to obtain at least a relative estimate of net erosion and deposition.We assume that we can estimate sediment flow at sediment transport capacity as

where* R ~ i ^{m}[N/h=kg m/s^{3}=c1*m/s, c1=kg/s^{2}],
KCP ~ K_{t}[kg/m^{3}]
* and

where **s **is a unit vector in the flow direction, *a* [deg]
is aspect of the terrain surface.This equation is equivalent to the relationship
with curvatures presented in the report by Mitasova et al. 1998, however,
the computational procedure is simpler. The exponents *m*
and *n* control the relative impact of water and slope terms and reflect
different erosion patterns for different types of flow. The exponents theoretically
consistent with the RUSLE are *m=1.6, n=1.3 *and they seem to
reflect pattern* *for prevailing rill erosion (net erosion is computed
as a derivative of T, which leads to an exponent 0.6 used in USLE) with
erosion sharply increasing with the amount of water (Figure 10a). The observed
extent of colluvial deposits in our previous study indicated lower exponents
reflecting compound, long term impact of rill and sheet erosion and both
large and small events. The exponents *m=n=1* predict greater spatial
extent of deposition suggesting prevailing sheet erosion (Figure 10b).
Comparison with SIMWE simulations shows that *m<1* leads to pattern
with small erosion/deposition rates but large spatial extent of deposition
a situation which may occur during a small event.

Caution should be used when interpreting the results because the USLE parameters were developed for simple plane fields and detachment limited erosion therefore to obtain accurate quantitative predictions for complex terrain conditions they need to be re-calibrated ( Foster 1990, Mitasova et al 1997 reply).

**2.2.3 Process-based simulation of water erosion**

The new generation erosion model SIMulation of Water Erosion model (SIMWE)
is based upon the description of water flow and sediment transport by first
principles equations (e.g., Foster and Meyer 1972, Bennet 1974) with the
core theoretical principles based on the WEPP (Flanagan and Nearing 1995).
The model is described by Mitas and Mitasova (1998) and in our previous
reports, therefore we briefly present only its principles and recent enhancements.
The steady state version SIMWE1.0 estimates the net erosion/deposition
*D( r)* at a point

(10)

under the assumption that the erosion and deposition are proportional to
the difference between the sediment transport capacity

(11)

where

(12)

where The model is being continually enhanced by including more complete and more general description of the simulated processes. The following enhancements were useful for the applications in this project:

- simulation of water depth in flat areas and depressions without a channel
- simulation of waterflow using combination of DEM and stream network to improve the accuracy of water flow in streams and rivers
- simulation of water flow and erosion/deposition for short events (before steady state)
- simulations at spatially variable resolutions (multi-scale)

We are also moving beyond the steady state modeling by introduction of real time and simulation of events with a given time period. This allows us to predict erosion/deposition patterns for small events when water flow has not reached the steady state, as illustrated by the following example of a field in Exeter, England, when center of the valley has depositional area for smaller events and gully erosion for large events.

The necessity for spatially heterogeneous modeling of fluxes can arise
due to various reasons:

(a) study area is represented by data with spatially variable resolution
and accuracy (study area has high resolution data however there is
mass coming in from an area for which we have only low resolution data);

(b) study area is large, with spatially variable complexity for which
simulation with spatially variable precision or detail is the most effective
(study area is large, but a certain part of it is homogeneous so it is
unnecessary to run it all at high resolution, only area of interest/highs
spatial variability will be simulated at high resolution while the rest
is low resolution.)

The system of continuity equations used in SIMWE describes the water
and sediment flow at a spatial scale equal or larger than an average distance
between rills (i.e., grid cell size more than 1m) and therefore the presented
approach allows us to perform landscape scale simulations at variable spatial
resolutions from one to hundreds of meters, depending on the complexity
and importance of studied subregions. The solution of continuity equations
through the Green's function can be reformulated for accommodation of spatially
variable accuracy and resolution by introducing a weighting function.
This weighting function can change (abruptly or smoothly) between
regions with unequal resolutions and in fact, can be optimally adapted
to the quality of input data (terrain, soils, etc.) so that the accurate
solution is calculated only in the regions with correspondingly accurate
inputs. The reweighted Green's function in effect, introduces
higher density of sampling points in the region with large weight.
The statistical noise will be spatially variable as a function
of the average number of samples resulting in the accuracy increase
for the areas with weight function > 1.

We have tested implementation in 2 systems :

- research oriented, GPL licensed, GRASS5.0,
- commercial system ArcView

In the following examples we assume that the given data are raster maps representing elevation, K, C, R and optionally also the P factors, and resolution as a constant. The erosion is always negative (including USLE), deposition is positive.

**r.flow** -u elevation lgu=length
**r.slope.aspect **elevation slope=slope aspect=aspect
**r.mapcalc**

slg1=if(K_mask,length)

slg2=if(K_mask,0,length)

length_av=if(K_mask,140,69)

usle_lg_av=-R*K*C*1.6*exp(length_av/22.13,0.6)*(65.4*sin(slope)*sin(slope)+4.56*sin(slope)+0.0654)

**r.flow** -u elevation lgu=length
**r.slope.aspect **elevation slope=slope aspect=aspect
**r.mapcalc**

usle_lg=-R*K*C*1.6*exp(length/22.13,0.6)*(65.4*sin(slope)*sin(slope)+4.56*sin(slope)+0.0654)

**r.flow **elevation dsout=flowacc
**r.slope.aspect **elevation slope=slope aspect=aspect
**r.mapcalc**

usle_area=-R*K*C*1.6*exp(flowacc*resolution./22.13,0.6)*exp(sin(slope)/0.0896),1.6)

**USPED (can be written as a script)**

**r.flow **elevation dsout=flowacc
**r.slope.aspect **elevation slope=slope aspect=aspect
**r.mapcalc**

sflowtopo=exp(flowacc*resolution,1.6)*exp(sin(slope),1.3) (prevailing rills)

or

sflowtopo=flowacc*resolution*sin(slope) (prevailing sheet)

qsx = R * K * C * sflowtopo * cos(aspect)

qsy = R * K * C * sflowtopo * sin(aspect)

erdep = qsx.dx + qsy.dy (prevailing rills)

erdep = 10. * (qsx.dx + qsy.dy) (prevailing sheet)

**Note: **Similarly as in the case of modified USLE, the model may
produce very high erosion in areas of concentrated flow, especially for
the prevailing rill erosion. These areas represent only a very small percentage
of modeled area and should be reclassed to values typical for gully erosion.

**2.3.2 ArcView/Spatial Analyst**

There are several aml scripts and other implementations of USLE in the ESRI products or as field based, on-line calculators available over the Internet. We present an implementation for raster data using the standard commands in ArcView with Spatial Analyst extension.

**USLE ***based on upslope area*

select elevation

under Analysis select DERIVE SLOPE give the
new theme name slope

MAP CALCULATOR

build an expression

([elevation].FlowDirection(FALSE)).FlowAccumulation(NIL)

Evaluate, give the new theme name flowacc

build an expression :

(([flowacc] * resolution/22.1).Pow(0.6))* ((([slope]*0.01745).Sin)/0.09).Pow(1.3))*1.6

Evaluate, give the new theme name lsfac

build an expression

-1.*R*[K]*[C]*[P]*[lsfac]

Evaluate, give the new theme name usle_area

**USPED**

select elevation

DERIVE SLOPE, name it slope

DERIVE ASPECT, name it aspect

MAP CALCULATOR

build en expression

([elevation].FlowDirection(FALSE)).FlowAccumulation(NIL)

Evaluate, name it flowacc

MAP CALCULATOR

build an expression (*prevailing rill erosion)*:

(([flowacc] * resolution).Pow(1.6))* ((([slope]*0.01745).Sin).Pow(1.3))

Evaluate, name it sflowtopo

or for *prevailing sheet erosion*

([flowacc] * resolution)* (([slope]*0.01745).Sin)

Evaluate, name it sflowtopo

*Important : the values of K and C should be multiplied by 100 (see
note below)*

build an expression

(((([aspect]*(-1))+450)*0.01745).Cos)*[sflowtopo]*[K]*[C]*R

Evaluate, name it qsx

(((([aspect]*(-1))+450)*0.01745).Sin)*[sflowtopo]*[K]*[C]*R

Evaluate, name qsy

select qsx

DERIVE SLOPE, name it qsx.slope

DERIVE ASPECT, name it qsx.aspect

select qsy

DERIVE SLOPE, name it qsy.slope

DERIVE ASPECT, name it qsx.aspect

MAP CALCULATOR

build an expression

(((([qsx.aspect]*(-1))+450)*0.01745).Cos)*([qsx.slope]*0.01745).Tan

Evaluate, name it qsx.dx

build an expression

(((([qsy.aspect]*(-1))+450)*0.01745).Sin)*([qsy.slope]*0.01745).Tan

Evaluate, name it qsy.dy

build an expression

[qsx.dx] + [qsy.dy]

Evaluate, name it erdep

It is necessary to keep in mind that the resulting values for erosion
and deposition are multiplied by the values used to multiply the C and
K factors. **Do not divide **to get the correct values, because the
pattern with values in interval <-1,1>, typical for well preserved areas
as well as for the relatively flat areas, will be lost. We were not able
to visualize the erosion/deposition map when using the true values of K
and C, because ArcView has troubles handling floating point numbers in
the interval <-1,1> . According to the Spatial Analyst manual, floating
point grid themes can only use equal interval or standard deviation for
classification which are unsuitable for erosion/deposition maps. Several
examples of problems encountered when computing erosion/deposition maps
with true values of K and C are presented in the online Appendix. It is
also important to note that the algorithms available in ARC are different
from those in GRASS, therefore to obtain acceptable result more attention
should be paid to the proper selection of resolution and to the quality
of DEM.

Because of profound changes in ESRI software which involves merging of ArcINFO and ArcView, release of ArcModel and our troubles with displaying the maps in Arcview3.1,3.2 we propose testing the model again under the new release of ArcInfo8 for which the GMSLab will be a beta tester.

**Process-based simulation **implemented as** SIMWE1.0 **is designed
as a research tool linked to GRASS and ArcInfo/ArcView through import/export
of raster data. Most of preprocessing, visualization and analysis is done
within GIS only the simulation of the water flow and sediment flow processes
is run by the stand alone program.

The approach used for SIMWE which is based on modeling with multivariate fields and the Monte Carlo solution of continuity equations opened new possibilities for simulation of a broad range of phenomena related to flow of water and sediment and their interaction with terrain and land cover. Therefore SIMWE is being continuously enhanced and further developed by adding more complete description of modeled processes and modifications reflecting the recent results of experimental research. The proposed multiscale approach can be implemented within any standard raster GIS without the need for developing new data structures. For each resolution a separate set of inputs is created with overlapping region, simulation algorithm combines the data from different files by going through several passes of computation. Algorithm produces a useful byproduct - identification of borderline segments with inflow of water and sediment and segments with outflow - important for identification of areas which require protection from e.g. inflow of harmful chemicals or areas where outflow of contaminated material should be controlled (Figure 13).

To perform the routine modeling efficiently, a customized interface
can be created for particular applications.

**Computation of Digital Elevation Model**. We have used 2 different
sources of elevation data : a) IFSARE from TEC and b) point data from direct
field measurement from the local NRCS. The IFSARE data from TEC were provided
as integer values creating artificial steps along 1m contourlines and making
the use of data for modeling of water and sediment flow problematic (Figure
15a). Upon our request, we were also provided with the original floating
point IFSARE data where the artificial steps were not present and data
were more suitable for modeling. However, both data sets contained speckle
noise typical for radar data which has to be removed to make the data useful
for applications. The noise was relatively homogeneous and appeared as
"bumps" about 1-3m high/wide on top of the measured terrain surface with
vegetation (Figure 15a,b). In spite of the fact that the RST method has
not been designed specifically to deal with the speckle noise, we have
found that the RST method performed very well as a smoothing tool allowing
us to create usable representation of terrain for the studied watershed
(Figure 15c).

High quality field data were provided by NRCS with coded fault and border lines which allowed us to test the RST method for surface interpolation with fault (Figure 4). The gully was mapped with a substantially higher density of points therefore we will be able to use this data set to test it for multiscale representation of terrain and modeling of flow on the hillslopes and inside the gully. Because of their higher accuracy and detail we present here the results of modeling using the NRCS field data with slightly smoothed gully. The model results for the IFSARE data are in the online Appendix.

**Erosion modeling**. We have focused on comparison of different
approaches based on USLE parameters and we have also run some simulations
with SIMWE. The USLE and USPED parameters were determined based on the
soil data and field survey and the following values were used: *R
= 160, K* for the following soils *Slidel and Topsey K=0.32,* *Brackett
K=0.17, * C for the upper part of the watershed was set to *C=0.2
*(20% cover) and for the central and lower part *C=0.04* (60% cover)

The Cs data were sampled at 100 points and erosion/deposition rates were derived using the models presented by Walling et al. 199?. ( see Warren et al., 1999). To visually compare the pattern of erosion and deposition derived from models and from the Cs data, we have interpolated the Cs-derived rates using RST (Figure 16)

We have computed the maps of erosion and erosion/deposition using all
the USLE methods described in section 2.2.1 and 2.3.1 with the results
presented in Figure 16 a, b, c. By definition, the results do not show
any deposition thus overestimating the spatial extent of erosion. The methods
based on slope length predict the lowest erosion rates in the area where
the erosion rates are the highest and gully is present. This demonstrates
that when using slope-length approach the areas of concentrated flow should
be identified (manually, from imagery or using flowtracing GIS tools) and
modeled separately and the estimates added to the total and average soil
loss. Such an approach would be extremely time consuming for large areas
typical for military installations. The method based on the upslope area
predicts the highest erosion rates in the center of the valley where the
gully is located, however, it does not predict the deposition so the total
erosion is overestimated.

Comparison with Cs data shows the missing deposition in the USLE results
(Figure 9), however the values were the closest to the averaged slope length
model - these results predicted relatively homogeneous distribution of
erosion rates similar to the Cs distribution and also has very little relation
with the amount of flowing water.

We have tested the USPED model for various values of *m,n* parameters
(Figure 10). The model predicts both erosion and deposition, approximately
in the areas indicated by the Cs measurements. However, the pattern from
the USPED is more complex than the pattern derived from the Cs data (partially
due to the fact that Cs was sampled by 100 points while the terrain was
sampled by over thousand points) and the erosion rates predicted by USPED
in the gully are substantially higher than the erosion rates on the ridges,
while the Cs estimates in these areas are about the same. We have found
that the lowest differences between the Cs and USPED were for the *m=0.6*
while the spatial extent of deposition was closest for *m=1 *(Table
1). This indicates impact of relatively local processes (in case of water
erosion, short events with little runoff) which is in agreement with the
published Cs research results.

erosion [%area] | deposition [%area] | |

Cs137 | 72 | 28 |

p=0.6 | 56 | 41 |

p=1.0 | 74 | 25 |

p=1.6 | 81 | 18 |

We have also computed an average from model results with *m=0.6, 1.0,
1.6* approximating an impact of combination of smaller and larger events
and prevailing sheet and rill erosion. The resulting estimate was
relatively close to the Cs distribution and at the same time the pattern
captured the gully erosion in the valley center. This indicates that it
may be important to use continuous time simulations which capture the changes
in rainfall and vegetation during the year, to predict the erosion/deposition
pattern accurately.

The results show that calibration of water erosion models based on Cs is problematic and our experience from this project as well as several studies demonstrate that only a portion of the Cs movement can be explained by water erosion and that relatively local processes have a substantial impact in Cs distribution. For example, in agricultural areas, tillage erosion which which is not present in our area, is used to explain the Cs distribution (Quinne et al. 1998). Even more intriguing is the fact that Cs is bound to clay particles which can be carried by water over long distances and deposited into the lakes and rivers, so at least theoretically, the pattern of clay particle movement has a long range character as opposed to local character of Cs distribution.

SIMWE simulations were ran with the IFSARE data and will be run for the more accurate field data and reported in the future. The SIMWE model will allow us to explore the effects which are not captured in the USLE or USPED such as adding detachment by raindrops which is dependent from water flow (although the term seems to be relatively low) or diffusive term reflecting combined impact of various processes (weathering, wind, raindrop splash) which may be important for long term effects.

A complete set of maps computed from both the IFSARE and the NRCS field
data is available in online Appendix
1.4.1 and Appendix
1.4.2.

- 20m resolution 1000x2000 grid which can be handled by current upper end workstations (Figure 18.)
- 5m resolution (to preserve smaller detail) with 5000x7000 grid which requires high end workstation with large RAM and disk space.

The 5m resolution DEM is aimed at studies for smaller subareas for which the sub-DEM can be "clipped" from the large file which we have provided. The 5m resolution DEM also served as an excellent test for the capabilities of s.surf.rst to process such a large data set in a reasonable time without the necessity to split the area. We have computed the topographic parameters slope and aspect for both models simultaneously with interpolation. We are now processing the LCTA data and comparing the field measured slopes with the values derived from the DEM to ensure consistency in input data when the erosion estimates from models are compared to the LCTA estimates. Preliminary comparison shows significant differences in slope at several LCTA sites and further work is necessary to ensure that the data have been imported and compared correctly.

We have used the 20m DEM to perform analysis of topographic risk for detachment limited erosion and net erosion/deposition (Figure 19). These results can be combined with the land use data to estimate the actual erosion risk under the current or simulated/planned land use.

We have selected a subarea to simulate water and sediment flow as well as erosion/deposition under homogeneous conditions in a selected subarea using SIMWE (Figure 11). Water flow simulation shows topographic potential for wetlands in several locations along the streams and smaller "spots" in the upland locations.

Modeling effort for Fort Polk will continue with the new land use data when they become available and with studies of impact of various land use management alternatives and BMP implementations.

The current research has shown that there are gaps in the understanding of sediment transport processes in complex landscapes and further research which integrates experimental and modeling approaches is needed. For example, the simulations have shown that for all models more general equations for estimation of detachment and transport capacity is needed, so that numerous equations, each applicable only to very specific conditions, are not necessary. This is an area of active research with the development of a new, more general equation for slope factor for RUSLE by. Nearing et al. 1999 a good example. There is also a strong need for calibration of the new models and equations for computation of their components both quantitatively and spatially to increase the accuracy of predictions (currently at 50-150% error) to acceptable levels.

The implementation of models within a GIS demonstrates the capabilities
to compute maps of erosion risk for large areas efficiently, however improvements
in GIS capabilities are needed, especially in handling and displaying the
floating point data in ArcView. With the completely redesigned ArcInfo
and ArcView being released next year for which the GMSLab is being the
beta tester we recommend to focus on the improved implementation of erosion
models in the new system rather than in ArcView3.1 and Avenue.

Bennet, J. P., 1974, Concepts of Mathematical Modeling of Sediment Yield, Water Resources Research, 10, 485-496.

Desmet, P. J. J, and Govers, G., 1995, GIS-based simulation of erosion and deposition patterns in an agricultural landscape: a comparison of model results with soil map information. Catena 25, 389-401

Desmet, P. J. J., and G. Govers, 1996, A GIS procedure for automatically calculating the USLE LS factor on topographically complex landscape units, J. Soil and Water Cons., 51(5), 427-433.

Flanagan, D. C., and M. A. Nearing (eds.), 1995, USDA-Water Erosion Prediction Project, NSERL, report no. 10, pp. 1.1- A.1, National Soil Erosion Lab., USDA ARS, Laffayette, IN.

Foster, G. R., and L. D. Meyer, 1972, A closed-form erosion equation for upland areas, in Sedimentation: Symposium to Honor Prof. H.A.Einstein, edited by H. W. Shen, pp. 12.1-12.19, Colorado State University, Ft. Collins, CO

Foster, G. R., 1990, Process-based modelling of soil erosion by water on agricultural land, in Soil Erosion on Agricultural Land}, edited by J. Boardman, I. D. L. Foster and J. A. Dearing, John Wiley & Sons Ltd, pp. 429-445

Haan, C. T., B. J. Barfield, and J. C. Hayes, 1994, Design Hydrology and Sedimentology for Small Catchments, pp. 242-243, Academic Press.

Julien, P.Y., and Simons, D.B., 1985, Sediment transport capacity of overland flow. Transactions of the ASAE, 28, 755-762.

Julien, P. Y., B. Saghafian, and F. L. Ogden, 1995, Raster-based hydrologic modeling of spatially varied surface runoff, Water Resources Bulletin, 31(3), 523-536.

Mitas, L., Mitasova, H., 1999, Spatial Interpolation. In: P.Longley, M.F.Goodchild, D.J. Maguire, D.W.Rhind (Eds.), Geographical Information Systems: Principles, Techniques, Management and Applications, Wiley, 481-492.

Mitas, L., and Mitasova, H., 1998a, Distributed soil erosion simulation for effective erosion prevention. Water Resources Research, 34(3), 505-516.

Mitas, L., Mitasova, H., 1998b, Multi-scale Green's function Monte Carlo approach to erosion modelling and its application to land-use optimization In: Modelling Soil Erosion, Sediment Transport and Closely Related Hydrological Processes, (Eds: W.Summer, E. Klaghofer and W.Zhang), IAHS Publication no. 249, pp. 81-90.

Mitas L., Brown W. M., Mitasova H., 1997, Role of dynamic cartography in simulations of landscape processes based on multi-variate fields. Computers and Geosciences, Vol.23, No. 4, pp. 437-446

Mitasova, H., Mitas, L., Brown, W. M., Johnston, D., 1998, Multidimensional Soil Erosion/deposition Modeling and visualization using GIS. Final report for USA CERL. University of Illinois, Urbana-Champaign, IL.

Mitasova, H., J. Hofierka, M. Zlocha, L.R. Iverson, 1996, Modeling topographic potential for erosion and deposition using GIS. Int. Journal of Geographical Information Science, 10(5), 629-641. (reply to a comment to this paper appears in 1997 in Int. Journal of Geographical Information Science, Vol. 11, No. 6)

Mitasova, H., J. Hofierka, 1993, Interpolation by regularized spline with tension : II. Application to terrain modeling and surface geometry analysis. Mathematical Geology 25, p. 657-669.

Moore I.D., and Burch G.J., 1986b, Modeling erosion and deposition: Topographic effects. Transactions ASAE, 29, 1624-1640.

Moore, I.D., Turner, A.K., Wilson, J.P., Jensen, S.K., Band, L.E., 1993, GIS and land surface-subsurface process modeling. In GIS and Environmental Modeling, edited by Goodchild, M.F., Parks, B., Steyaert, L.T., Oxford University Press, New York, 196-203.

Moore, I.D., and Wilson, J.P., 1992, Length-slope factors for the Revised Universal Soil Loss Equation: Simplified method of estimation. J. Soil and Water Cons., 47, 423-428.

Nearing M.A., Norton L.D., Bulgakov D.A., Larionov G.A., West L.T., Dontsova K.M., 1997, Hydraulics and Erosion in eroding rills. Water Resources Research, 33, 865-876.

Quine, T. A., P. J. J. Desmet, G. Govers, K. Vandaele, and D. E. Walling,1994,
A comparison of the roles of tillage and water erosion in landform development
and sediment export on agricultural land near Leuven, Belgium,

in Variability in stream erosion and sediment transport, edited by
L. J. Olive, R. J. Loughran, and J. A. Kesby, IAHS Publication, no. 224,
Proc. of an Int. Symposium Canberra, Australia, pp. 77-86.

Saghafian, B., Implementation of a Distributed Hydrologic Model within GRASS, in GIS and Environmental Modeling: Progress and Research Issues, edited by M. F Goodchild, L. T. Steyaert, and B. O. Parks, GIS World, Inc., pp. 205-208, 1996.

Srinivasan, R., and Arnold, J.G., 1994, Integration of a basin scale water quality model with GIS, Water Resources Bulletin, 30(3), 453-462.

Westervelt, J, Shapiro, M., 1991, r.mapcalc, Tutorial, USA CERL.

A1.2 USPED for GRASS and ArcView

A1.3 GIS tools: enhanced s.surf.rst, r.flow, r.slope.aspect, r.enforce

A1.4 WWW documents

elevation: IFSARE: original, smoothed, field data smoothed

topographic analysis: slope, aspect, upslope area

erosion: K, LS, USLE, USPED, SIMWE

**Ft. Polk**

elevation: 20m and 5m from points+contours

topo analysis: slope, aspect, upslope area,

erosion: LS, topo erdep, SIMWE test