Content

1. Digital landscape characterization and analysis

Land use management and landscape process modeling supported by Geographic information system (GIS) technology require representation of landscape and its features in a digital form as point, vector or raster data. Digital landscape characterization data are usually obtained from direct field measurements, remote sensing imagery or digitized maps. Recently, large volumes of digital spatial data have become available on Internet, for example through National Spatial Data Clearinghouse or Microsoft TerraServer. To effectively use these digital data for modeling and simulations needed for land use management we have introduced the concept of spatial phenomena representation by multivariate functions f(r), where r=(x,y,z,t), with effective discretization as multidimensional rasters (Mitas and Mitasova 1997), as illustrated by FIGURE 1. In the following sections we discuss the specific issues related to digital landscape characterization needed for erosion modeling.


1.1 Phenomena influencing erosion and their digital representation

Erosion, sediment transport and deposition involves complex interactions between rainfall, surface and subsurface hydrology, soil properties, land cover and topography. Modeling erosion and deposition in complex terrain within a GIS therefore requires a high resolution digital elevation model (DEM), as well as digital data describing spatial distributions of rainfall, soil and land cover properties, including the anthropogenic features.

Rainfall data are obtained from weather stations or simulated by weather simulators based on statistical analysis of long term rainfall data in a given area. For smaller areas, only uniform rainfall is used, for larger areas (regional level) the data from several stations are interpolated to create spatial representation of rainfall distribution in the raster format. For areas with negligible influence of terrain we use 2D spline interpolation (Mitasova and Mitas 1993), for areas with significant terrain impact on the distribution of rainfall, the use of more complex 3D spline interpolation is needed (Mitasova et al. 1995).

Elevation data with resolutions 30m and lower (e.g., 100m) are readily available for most areas from USGS. Horizontal and vertical resolution of these data is suitable for rough erosion risk assessment at a regional level. Resolution and accuracy of DEMs available from USGS is often insufficient for more detailed, installation/training area level modeling (Mitasova et al. 1996), because systematic errors create various artificial patterns in terrain surfaces, such as stripes or waves along contour lines leading to unrealistic predictions of erosion/deposition distribution. New type of elevation data produced by IFSARE provides significantly higher resolution as well as information on man made features, however, the data are distorted by speckle noise and represent surface including the vegetation. Therefore the use of these data is not straightforward and will require significant investment in new research and development before they could be used effectively. Recent advances in laser altimetry provide more accurate estimation of elevations thank the IFSARE and tools to remove the vegetation from the measured surface are already commercially available. Unfortunately, laser data are not yet widely available and cost is still relatively high, however their future use for high resolution erosion modeling is very promising. Currently, it seems that the most reliable and accurate source of elevation data are contours generated from orthophotomaps and interpolated to high resolution DEM (2-10m, depending on the contour interval) by the RST method, as described later. For small areas (fields, small watersheds) field measured elevation data are the most appropriate.

Land cover data. Vegetation and land use data can be often derived from satelite imagery at 30 - 20m resolution (Warren et al. 1989). Higher resolution data can be obtained from field surveys, orthophotomaps or digitized from thematic maps. We expect that the availability of high resolution land cover data will greatly increase with the rapid development of new sensors and remote sensing technology. Land cover is a highly dynamic phenomenon, so temporal series of land cover data is extremely useful for evaluation of erosion risk. Because temporal data are not currently readily available, temporal variability in land cover due to the vegetation growth is often simulated by plant growth models (Flanagan and Nearing 1995). Impact of land cover on erosion is complex and difficult to represent by a few parameters. Therefore transformation procedures or tables are needed to derive erosion model parameters such as C-factor, surface roughness factor (e.g., Mannings, Chezy coefficients), erodibility, transportability and critical shear stress adjustment factors. While extensive sets of these parameters are available for USLE, CREAMS and WEPP models, an analysis of the suitability of the published values of these parameters to conditions generated by military use needs to be performed to ensure the validity of the simulation results.

Soil data. Soil data are usually derived from thematic maps and soil data databases (e.g. STRATSGO, WEPP), however, high accuracy simulations may require field measurements of soil properties. Because the impact of soils properties is less dramatic than the impact of cover lower resolution data are often sufficient.

The WEPP manual (Flanagan and Nearing 1995) describes the procedures used for the development of soil and cover parameters, and gives empirical equations representing relations between erodibility and soil texture and organic matter content. It also presents numerous adjustment factors reflecting the impact of vegetation and various erosion prevention measures and agricultural practices. WEPP program includes a data base which automatically assigns the values of parameters for the given soil, vegetation characteristics and practices. Similar database can be created for military land use, and thus simplify the task of selecting the proper parameters.

To prepare the input representing soil and cover properties for erosion simulation it is necessary to provide the data in the form of a raster map with the resolution identical to the resolution of DEM. This is straightforward, if the data are in a raster or vector format, as usually simple GIS transformation programs for resampling and vector-to-raster transformation can be used. However, if the data are given in the form of measurements at irregularly distributed points, spatial interpolation should be used. While there exist a large number of methods for spatial interpolation, the spline and kriging methods offer the most flexibility and the highest accuracy.


1.2 Multivariate spline interpolation


To model erosion using a GIS an adequate discrete representation of terrain, soil and cover properties, and climatic phenomena must be created in a GIS database. These phenomena are usually measured at irregularly distributed points or digitized from isoline maps. However, models of landscape processes such as water and sediment flux are most often designed for grid data. Therefore, the following problem must be solved:

Given site data representing a continuous phenomenon (x,y,z)i, i=1,...N, where (x,y) are a geographically referenced coordinates and z is a measured value (elevation, precipitation, concentration, potential etc.) create a raster map representing the spatial distribution of this phenomenon, that means interpolate or approximate the values of z at regular grid points. This task is difficult to perform especially for geoscientific data because of the following reasons:

Various methods have been developed to solve this problem (Mitas and Mitasova 1999), however general solution giving perfect results for any data does not exist. Illustration of the results obtained by various methods available in GIS is given in the FIGURE 2.

Significant improvement in interpolation of digital elevation models and other scattered data has been achieved by new spline interpolation and approximation methods. Spline methods are based on the assumption that the approximation function should a) pass as closely as possible to the data points and b) should be at the same time as smooth as possible. These two requirements can be combined into a single condition of minimizing the deviation from the measured points and the smoothness seminorm of the function. The solution of this problem can be expressed as a sum of two components (Talmi and Gilat 1977) a 'trend' function and a radial basis function with an explicit form which depends on the choice of the smooth seminorm. For our choice of the smooth seminorm (Mitasova and Mitas 1993) the `trend' function is a constant and the general form of radial basis function for d-variate case (Mitasova et al. 1995) includes incomplete gamma function (Abramowitz and Stegun 1964). We call this function Regularized Spline with Tension (RST). For the the special case d=2 the radial basis function can be expressed through a logarithmic and exponential integral function (Abramowitz and Stegun, 1964). The explicit form for the function for 3D (volume) and 4D (volume in time) data has been presented by Mitasova and Mitas (1993) and Mitasova et al. (1995). The coefficients of the interpolation function are obtained by solving the system of N linear equations, where N is the number of given points.

Besides the high accuracy (Mitasova and Mitas 1993, McCauley 1995, Rohling 1998), the Regularized Spline with Tension has several useful properties. The generalized tension parameter controls the distance over which the given point influences the resulting surface. For the bivariate case tuning the tension can be interpreted as tuning the character of the resulting surface between a membrane and a thin plate. Smoothing parameter can be used to interpolate smooth surfaces from noisy data. This parameter can be spatially variable depending on the noise/accuracy for each data point. The proper choice of smoothing and tension parameters is important for successful interpolation or approximation. By tuning the tension and smoothing it is possible to minimize the overshoots (Mitasova and Mitas 1993), artificial pits or banding effect of the elevation values around the contours, observed for the less general forms of splines such as the thin plate spline or thin plate spline with fixed tension. GIS implementation of this function provides several tools for controlling the quality of resulting DEMs, such as computation of deviations, predictive errors based on ordinary crossvalidation procedure, and computation of curvatures useful for detecting artificial waves around contours caused by improperly chosen tension parameter. Optimal smoothing parameter can be found by ordinary or general cross-validation scheme (Wahba 1990, Mitasova et a. 1995).

Due to the solution of system of linear equations the computational demands for the presented method are proportional to N^3 and their application to large data sets becomes problematic. A solution to interpolation of large data sets using an explicit interpolation function was proposed by Franke (1982) and further developed and applied in geosciences by Mitasova and Mitas (1993). The approach is based on the division of given region into rectangular segments and computation of interpolation function for each segment using the data points from this segment and from its neighborhood. Because the method works with equally sized segments, the approach is not very efficient for strongly heterogeneous data like digitized contours or clustered drill hole data. A more effective approach to decomposition of the given region with heterogeneous spatial distribution of data into segments with approximately the same number of points for each segment was designed with quadtrees by Gerdes and Kosinovsky, (Mitasova et al. 1995). The interpolation function is then computed for each segment using the points from this segment and the points located in the window which adjusts its size to the density of points in the neighborhood of each segment. The segmentation is sensitive to tension, for very low tension parameter used in flat areas, the number of points needed for smooth connection of segments can be very large. In certain extreme situations where empty flat segment is next to a segment with rapidly changing terrain, additional smoothing may be necessary. To improve the performance of program for such situations the segmentation procedure was further enhanced by using automatically adjustable number of data points dependent on the size of each segment.

RST has been specially designed to meet the demands of topographic analysis. It has continuous derivatives of all orders allowing us to use these derivatives to compute slope, aspect and curvatures simultaneously with interpolation, as described in the next section.

1.3 Topographic analysis

Topographic parameters needed for erosion modeling include parameters describing the local geometry of elevation surface: slope, aspect, curvatures and parameters related to the flow over the terrain surface: upslope contributing area, slope length (FIGURE 3).

Local geometry parameters. Regularized spline with tension was specially designed to support direct computation of topographic parameters which are functions of first and second orders derivatives. Using the partial derivatives of the RST fucntion, slope angle is computed as arctan of the magnitude of the elevation surface gradient, and aspect angle is computed as arctan of the direction of minus gradient (Mitasova and Hofierka 1993). Computation of curvatures is more complicated because, in general, the surface has different curvatures in different directions and which one is important must be determined according to the type of processes under study. For applications in geosciences, the curvature in gradient direction (profile curvature) is important because it reflects the change in slope angle and thus controls the change of velocity of mass flowing down along the slope curve. The curvature in a direction perpendicular to the gradient reflects the change in aspect angle and influences the divergence/convergence of water flow. This curvature is usually measured in the horizontal plane as the curvature of contours and is called plan curvature (Zevenbergen and Thorne, 1987, Moore et al., 1991). However, for the study of flow divergence/convergence, it is more appropriate to introduce a curvature measured in the normal plane in the direction perpendicular to gradient (Krcho 1991). This curvature will be called here a tangential curvature because the direction perpendicular to gradient is, in fact, the direction of tangent to contour at a given point. Equations for these curvatures can be derived using the general equation for curvature of a plane section through a point on a surface as described by Mitasova and Hofierka 1993. The positive and negative values of profile and tangential curvature can be combined to define the basic geometric relief forms ( Krcho, 1991; Dikau, 1989). Each form has a different type of flow. Convex and concave forms in gradient direction have accelerated and decelerated flow, respectively, and convex and concave forms in tangential direction have converging and diverging flow, respectively.

Flow-related parameters. Derivation of flow-related parameters, such as upslope contributing area and hillslope length requires implementation of effective flow-tracing algorithm. Standard algorithms for flow tracing, which use only a limited number of flow directions (most often 8) from each grid cell, can lead to various unrealistic situations, such as prevailing flow in the direction parallel to x or y axes or diagonals (Fairfield and Leymarie 1991). Several new algorithms for flow tracing help to overcome deficiencies of standard algorithms by using the random-eight node approach (Fairfield and Leymarie 1991), multiple nearest neighbor nodes (Freeman 1991), and by using 360 directions of flow with the vector-grid algorithm (Mitasova and Hofierka 1993, Mitasova et al. 1996).

Upslope contributing area is the area from which the water flows into a given grid cell. Upslope contributing area per unit contour width Aj, for the given grid cell j is computed from the sum of grid cells from which the water flows into the cell j (Moore et al. 1991). This approximation is acceptable if the DEM is interpolated with the adequate resolution which depends on the curvature of terrain surface. Our experience, supported by several recent studies (Zhang and Montgomery 1994), is that the resolution 2 m-20 m is appropriate for models using upslope contributing areas in regions with complex terrain. For the computation of the number of cells draining into each grid cell, flow lines are constructed downhill from each grid cell in the minus gradient direction, until they reach a cell with a slope lower than the specified minimum, boundary line, singular point, or barrier (e.g. road). An improved algorithm for the construction of flow lines based on vector-grid approach (Mitasova and Hofierka 1993) is used. Flow lines constructed by this algorithm are represented in vector format. The points defining the flow line are computed as the points of intersection of a line constructed in the flow direction given by aspect angle and a grid cell edge. Downhill flow lines merge in valleys and they can be used also for the extraction of channels (Jenson and Domingue 1988).

Flow path length is used in the standard form of USLE/RUSLE and is appropriate for hillslopes with negligible water flow convergence/divergence. For the computation of flow path length, flow lines are generated uphill from each grid cell in the gradient direction until they reach a grid cell with the slope lower than the specified minimum, boundary line, singular point or barrier. Flow path length for each grid cell is then computed from the coordinates of points of their intersection with grid cells. Uphill flow lines are merging on ridge lines.

1.4 GIS implementation of interpolation and topographic analysis

GRASS GIS.
RST for bivariate interpolation was first released with GRASS4.1 as s.surf.tps. The original program was completely redesigned and enhanced for the GRASS5.0 release which includes support of floating point raster data, necessary for efficient erosion modeling. The new version of the RST interpolation program was renamed to s.surf.rst (Appendix) and it is currently being released with the GRASS5.0 produced by Baylor University, TX. The most important enhancements include: output of floating point rasters representing elevation, slope, aspect, profile and tangential curvatures as well as first and second order partial derivatives of the modeled surface, output of a site file with deviations in given points ; enhanced segmentation including an option to output the quadtree segments, spatially variable smoothing supporting approximation of surfaces from data with heterogeneous accuracy and easier to modify source code. If the raster DEM is available, local geometry parameters (slope, aspect, curvatures) can be computed in GRASS by r.slope.aspect (Appendix). The enhancements for the GRASS5.0 release include support for floating point, computation of curvatures and estimation of partial derivatives. Vector-grid flowtracing algorithm for computation of upslope contributing area, hillslope length and flowlines was implemented as a program r.flow (Appendix ) and is being released with GRASS5.0.

ARC/INFO, ArcView.
Spatial interpolation can be performed by TOPOGRID, which is a different version of thin plate spline with tension, developed by Hutchinson, 1989. The DEM derived by TOPOGRID should be used with caution for erosion/deposition modeling, because it tends to create waves along contours leading to artificial spatial pattern. Thin plate spline with tension and regularized spline developed by Mitas and Mitasova (1988) are available in ARCGRID and Arcview Spatial Analyst as SPLINE command. These functions do not include smoothing capabilities and topographic analysis which have to be performed separately after the surface is interpolated, using the functions SLOPE, ASPECT, CURVATURE. Flow related parameters are derived using the functions based on D8 algorithm (Jenson and Domingue 1988) FLOWACCUMULATION, FLOWLENGTH. The D8 algorithm works well for low resolution data, at high resolutions it tends to create artificial flow patterns biased towards grid cell diagonals.

1.5 Multidimensional dynamic visualization

Implementation of landscape process simulation tools within a GIS has stimulated integration of GIS and computer cartography with scientific visualization (Brown and Gerdes 1992 (SG3d) Brown and Astley 1995 (NVIZ), Brown et al. 1995, Mitasova, Brown, and Hofierka 1994). The integration provides methods and tools for creating dynamic cartographic models representing landscape phenomena and processes. These tools display multiple dynamic surfaces and isosurfaces, together with draped raster, vector and point data in an appropriate projection of 3D space with adequately selected parameters of the visualization environment. Visual analysis of input data is facilitated by interactive manipulations of visualization environment parameters such as viewing position, z-scale, cutting planes, light position and brightness etc., and by animating the sequences of images created by changing the viewing parameters or by displaying evolving series of data (Mitasova, Brown, and Hofierka 1994, Brown et al. 1995). Rendering shaded surfaces of input data prior to processing reveals the types of noise and other common landscape data anomalies discussed in section 1.1.2.

The new tools are used to create dynamic cartographic models of multi-variate landscape phenomena characterized by sets of discrete sampling points and by interpolated surfaces/hypersurfaces. The ability to visualize the data points in 3D space, fly through the point space, select certain values (e.g., time) as a third dimension and assign various graphical variables (size, color, shape, numerical label) to selected attributes (Brown 1996) facilitates visual analysis of spatial relationship between the location of point data and their measured values and helps in evaluation of interpolation techniques used to generate the surfaces.

Better visualization of point data for validation of source data. For visualization purposes, point data is often used to derive an artificial surface or a volume from which isosurfaces may be calculated. Such data extrapolation and interpolation is also necessary for modeling purposes to produce input data for models that work with grids or volumes. Better visualization of point data for comparison to derived surfaces at various stages of the modeling process allows the researcher to better understand and verify the model.

Building on recent improvements to the GIS library that work intelligently with spatial site data containing multiple attributes and more than two dimensions, visualization tools have been improved to allow multiple attributes/dimensions of point data to determine the shape, size, color and orientation of solid 3D markers as well as the 3D position. Additional tools have been developed to allow selection and display of a subset of sites from a larger file based on ranges of values for various attributes of each site (see Appendix). Given a value for an attribute of the data at a particular site, visualization is accomplished by mapping the value to one or more characteristics of the marker. In the case of color, three attributes (representing red, green, and blue components) may be used to map a single characteristic, color. For each data attribute, the visualization tools allow a transformation consisting of an addvalue and a multiplyvalue which converts the data value to visualization units, preventing duplication of data in the database solely for the purpose of visualization and allowing the researcher to more freely explore visualization combinations. If a single value is being used to map color, the existing format in the GIS for specifying color tables for raster data is used. Variable shape is accomplished by having an "alternative" size field which is used to scale the marker vertically. Plans to allow scaling the marker differently in each dimension according to attributes of the site data will further enhance the variability of shape. Using time stamps representing data measurement times, it is possible to visualize data collection sites and measured values in compressed time, drawing the site markers at intervals proportionally representing the time of measurement, and perhaps using derived data as a background surface. Such methods could highlight, for example, a measurement interval in which instruments were improperly calibrated, a period of sparse measurement, inadvertant influence due to spatial measurement patterns, or dramatic change that would need to be investigated further.

Animation of data is used to represent a change in time, change in a modeling parameter, or change in visible data. Producing an animation of the data involves redrawing the scene repeatedly with slightly changing data or viewpoint, saving each resulting image, and then combining them with existing multimedia tools into a movie. Our visualization tools incorporate a scripting mechanism which records the actions of the user while also allowing the use of loops. An integrated tool in the GUI provides methods for building the script, including functions such as Open Loop, Close Loop, Open File Sequence Loop, and Close File Sequence Loop, which help in the creation of animations by allowing a sequence of actions to be applied to a sequence of data files. Similarly, a File Sequence Tool allows iteration of data associated with the raster, vector, and site maps for producing the most common type of data visualization animation, where data has been precalculated to represent a timestep or change in modeling parameter and there will be no change in viewer position. This tool also makes it easy to add a somewhat constant dataset such as roads for better spatial reference. For more details on the use of animations for erosion modeling see Mitas et al 1996 and Mitas et al 1997.

Initially, tools were implemented on the Silicon Graphics platform using a proprietary graphics library, IRIS GL. While the interface language, Tcl, is publicly available for use on multiple platforms, the underlying graphics routines needed to be converted to OpenGL, a widely available graphic standard. Most of the tools have now been converted to use the OpenGL graphics API.


2 Distributed models of water and sediment transport

Erosion modeling within GIS focuses on describing spatial distributions or patterns of water and sediment flow as well as net erosion/deposition, along with estimating the statistical averages of magnitudes of soil loss. Capability to predict the pattern of erosion and deposition and to indentify the location of high risk areas is extremely important from the point of view of land use management because it provids the information necessary for effective implementation of erosion prevention measures. Recent advances in GIS technology, especially support for modeling with multivariate fields (Mitasova et al., 1995; Mitas et al., 1997), along with the exponential growth of computational power, stimulate the shift from empirical, lumped models to physically-based, distributed simulations, capable to predict the spatial distributions of landscape phenomena resulting from processes acting under various conditions. A typical example of this trend is the current development in hydrologic and erosion modeling (Maidment, 1996, Moore et al., 1993). Integration of distributed hydrologic models within GIS (Vieux et al., 1996; Saghafian, 1996) improved capabilities to estimate overland flow in complex landscapes. The development of analogous distributed simulation tools for erosion, sediment transport and deposition by overland flow has been slower, mostly because of the enormous complexity of the processes involved.

Several erosion prediction tools have been interfaced with GIS, including various modifications of an empirical model represented by the Universal Soil Loss Equation (USLE) (e.g., Warren et al. 1989, Flacke et al. 1990, Desmet and Govers 1996, Huang and Ferng 1990), and watershed models for non-point source pollution such as AGNPS or ANSWERS (de Roo et al. 1989, Rewerts and Engel 1991, Srinivasan and Engel 1991) which are based on the USLE concepts.

While important progress towards physically based erosion simulations was made by the introduction of a new generation erosion prediction tool within the Water Erosion Prediction Project (WEPP) (Foster et al., 1995), and by linking the existing simulation programs with GIS (Engel, 1995; Rewerts and Engel, 1991; Srinivasan and Engel, 1991), several problems persisted. Most models, including the distributed ones, are based on a 1 D sediment routing over planar hillslope segments (Hairsine and Rose, 1992; Govindaraju and Kavvas, 1991; Foster et al., 1995; Flacke et al., 1990) or through homogeneous (lumped) subwatersheds (Arnold et al., 1993). Due to these simplifications, most of the currently used models could only partially explain the impact of a complex spatial variability in terrain and land cover at a landscape scale, where even small variations can have a dramatic impact on location and rates of erosion and deposition. Unit stream power based methods and related approaches (Moore and Burch, 1986; Moore and Wilson, 1992; Mitasova et al., 1996a; Desmet and Govers, 1995; Hofierka and Suri, 1996) partially incorporated the influence of complex terrain. While this simple approach provides useful estimates of topographic effects on erosion/deposition, it does not simulate the impact of a wide range of soil and cover properties, typical for anthropogenic landscapes (Desmet and Govers, 1995). Several distributed landscape evolution models also describe erosion processes (e.g., Kirkby, 1986; Kramer and Marder, 1992; Howard, 1994), however, their application to land management is rather limited because these models are designed for time scales of thousands of years and focus on the evolution of topography and stream networks, rather than on predictions of erosion risk locations and erosion prevention.

Many existing methods require large number of empirical inputs, important effects of terrain shape are not taken fully into account, numerical methods used for the model implementations exhibit instabilities when used at high resolutions, and detailed spatially variable outputs are limited. Combination of these factors often make applications of traditional models either laborious or simply not practical, especially for large areas. Some of the key issues which restrict the use of the current models to increase our knowledge about spatio-temporal aspects of erosion processes and limit their application to effective erosion prevention for large, complex areas can be identified as follows:

It is clear that resolving these complex issues requires considerable time and effort of many research teams, therefore we have focused our effirt on a selected subset of the problems, especially those which we believe were relevant to erosion modeling at military installations.

2.1 Universal Soil Loss Equation (USLE)

The Universal Soil Loss Equation is an empirical equation designed for the computation of soil loss in agricultural fields. This equation was developed for the detachment capacity limited erosion in fields with negligible curvature and no deposition. The equation has the following form (Wischmeier and Smith 1978, Renard et al. 1991)

A = RKLSCP(1)

where R is the rainfall intensity factor, K is the soil factor, LS is the topographic (slope-length) factor, C is the cover factor and P is the prevention practices factor. Various modifications of this equation are often applied to the estimation of soil loss using GIS.

The topographic factor for USLE has been recently improved by incorporation of the influence of profile convexity/concavity using segmentation of irregular slopes and by improving the empirical equations for the computation of LS factor (Foster and Wischmeier 1974, Desmet and Govers 1996, Renard et al. 1991) as a part of the Revised Universal Soil Loss Equation (RUSLE). To incorporate the impact of flow convergence, the hillslope length factor was replaced by upslope contributing area (Moore and Burch 1996, Desmet and Govers 1996). The modified equation for computation of the LS factor in GIS in finite difference form for erosion in a grid cell representing hillslope segment was derived by Desmet and Govers (1996). A simpler, continuous form of equation for computation of the LS factor at a point on a hillslope, (Mitasova et al. 1996) is

LS(r) = (m+1) [A(r)/a]m [sin b(r) /b0]n (2)

where m and n are parameters, a is the length and b0 is the slope of the standard USLE plot. Impact of replacing the slope length by upslope area is illustrated in FIGURE 4 which shows that the upslope area better reflects the impact of concentrated flow on increased erosion. However, both USLE and RUSLE consider erosion only along the flow line without the full influence of flow convergence/divergence and both the standard and modified (e.g., Desmet and Govers 1996) equations can be properly applied only to areas experiencing net erosion. Depositional areas should be excluded from the study area. Therefore, direct application of USLE/RUSLE to complex terrain within GIS is rather restricted.

It has been 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 (Moore and Wilson 1992), for slopes with negligible tangential curvature. It is therefore possible to use the detachment capacity index as the LS factor for the computation of soil loss using the RUSLE with the influence of rainfall intensity, soil properties, and land cover, incorporated as RUSLE R, K, C factors respectively, with the assumption that transport capacity exceeds detachment capacity everywhere and erosion and sediment transport is detachment capacity limited.

2.2 Unit Stream Power based Erosion/Deposition model (USPED).

USPED 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. For this case, we assume that the critical shear stress is negligible and the sediment flow rate qs(r) is at the sediment transport capacity T(r), r=(x,y) which is approximated by a power function of slope b(r), water flow q(r) and transportability coefficient Kt(r) dependent on soil and cover (Julien and Simons 1985)

|qs(r)| = T(r) = Kt |q(r)|m sinb(r)n (3)

where m, n are constants depending on the type of flow and soil properties. For overland flow these parameters are usually set to m=1.5, n=1.3 (Foster 1993). Water flow can be expressed as a function of steady state water depth or as a function of upslope contributing area

|q(r)| = 1/n(r) h(r)5/3 sinb(r)1/2 = A(r) i(r) (4)

where n(r) is Manning's coefficient dependent on land cover and i(r) is rainfall intensity. Within the original USPED the water and sediment flow was modeled as 1D flow along a flow line generated over 3D terrain (Mitasova et al. 1996). The net erosion/deposition rate was computed as a change in the sediment flow rate along the flow line, approximated by a directional derivative of the sediment flow rate. For this univariate case, the net erosion/deposition rate D(r) is

D(r) = dT(r)/ds = Kt {[ grad h(r) . s(r) ] sin b(r) - h(r) kp(r)} (5)

where s(r) is the unit vector in the steepest slope direction, h(r) is the water depth estimated from the upslope area A(r), kp(r) is the profile curvature (terrain curvature in the direction of the minus gradient, i.e., the direction of the steepest slope). This 1D flow-based formulation includes the impact of water flow, slope and profile curvature, however, the impact of tangential curvature is incorporated only through the water flow term. The predicted pattern is in a good agreement with observations except for heads of valleys where it predicts only erosion while the soil maps and field experiments indicate that both deposition and erosion is observed. Also, for shoulders the 1D flow formulation does not predict commonly reported local maximum in erosion rates (Martz 1987, 1991, Sutherland 1991, Quine et al. 1994).

The USPED model was improved by using a 2D flow formulation, which was derived as a special case of the more general model SIMWE (Mitas and Mitasova 1998a, see next section ). Within this formulation we represent the water an sediment flow as a bivariate vector field q(r)=q(x,y), qs(r)=qs(x,y). Then, the net erosion/deposition rate is estimated as a divergence of the sediment flow. If we assume uniform rainfall, soil and cover conditions, and a transport capacity limiting case with sediment flow close to sediment transport capacity, the net erosion/deposition can be written as (see Appendix in Mitas and Mitasova 1998):

D(r) = div qs(r) = Kt { [grad h(r)] . s(r) sin b(r) - h(r) [kp(r) + kt(r)] } (6)

where kt(r) is the tangential curvature (curvature in the direction perpendicular to the gradient, i.e., the direction tangential to a contourline projected to the normal plane). Topographic parameters s(r), kp(r), kt(r) are computed from the first and second order derivatives of terrain surface approximated by the RST (Mitasova and Mitas, 1993; Mitasova and Hofierka, 1993; Krcho 1991). According to the 2D formulation, the spatial distribution of erosion/deposition is controlled by the change in the overland flow depth (first term) and by the local geometry of terrain (second term), including both profile and tangential curvatures. The bivariate formulation thus demonstrates that the local acceleration of flow in both the gradient and tangential directions (related to the profile and tangential curvatures) play equally important roles in spatial distribution of erosion/deposition. The impact of the tangential curvature kt(r) is therefore twofold. First, kt(r) influences the water depth h(r) through its control of water flow convergence/divergence, with tangential concavity leading to rapid increase in water depth and an increase in the potential for erosion. Second, kt(r) causes a local change in sediment flow velocity which for tangential concavity has an opposite effect (reduction in sediment transport), creating thus potential for deposition. The interplay between the magnitude of water flow change and both terrain curvatures reflected in the bivariate formulation therefore determines whether erosion or deposition will occur.

When the results of the 1D flow and the 2D flow models are compared with the observed pattern of colluvial deposits, it is clear, that the 1D model fails to predict deposition observed in areas where profile curvature is close to zero but there is a significant tangential concavity (FIGURE 5). It also underestimates erosion in areas with tangential convexity (shoulders) . The prediction by the 2D flow model in these areas is in significantly better agreement with the observed pattern of deposition. The 1D flow model underestimates the overall extent of deposition as only a 18% of the total area, while the 2D flow model, without taking into account the spatial variability in surface cover, predicts deposition at 26% of the total area. Colluvial deposits thicker than 10cm, indicating long term deposition, were observed in 40% of the sampling sites. The sedimentation area observed on a steep hillslope covered by grass, which contributes to the spatial extent of deposition cannot be predicted by using only the elevation data and the incorporation of spatially variable cover is necessary.

The 2D flow equation also provides a sound theoretical explanation for the results of field experiments reported by several authors, e.g., Busacca et al., (1993) and Sutherland, (1991); where `the highest erosion was observed on divergent shoulder elements and deposition on convergent footslope elements', as well as by Quine et al. (1994) observing `maximum soil loss from the slope convexities and maximum gain in both the slope concavities and the main thalwegs'.

GIS implementation of the more complete USPED equation is rather simple, using the tools already developed for this project and implemented in GRASS GIS. Assuming that we have an adequate DEM the procedure involves computation of:

In Arcview Spatial Analyst or ArcGrid the following computations are performed:

It is important to note that the algorithms available in ARC for interpolation and topographic analysis are less sophisticated that 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.

Soil and cover parameters similar to those used in USLE or WEPP were not developed for the USPED model, because no systematic experimental work was performed. However, it was suggested (Moore and Wilson 1992) that under certain conditions, there is a relation between the USPED concept and USLE. It is therefore possible to combine the model with the USLE/RUSLE factors to predict the relative impact of land use change on erosion/deposition pattern assuming the conditions valid for USPED hold. This approach was used e.g., for the Camp Shelby study and for the Hohenfels installations (see section on applications for more details). Caution should be used when interpreting the results because the USLE parameters were developed for simple plane fields and to obtain accurate predictions for complex terrain conditions they need to be re-calibrated ( Foster 1990, Mitasova et al 1997 reply).

The GIS implementation (both for GRASS and Arcview) of the computation of erosion and deposition for uniform or spatially variable rainfall, soil and cover properties can be customized for routine use by creating a graphical user interface (GUI) for the entire computation procedure and by automated transformation of soil and cover data to relevant erosion parameters. The GUI implementation can be developed using Tcl/Tk for GRASS and Avenue for Arcview.

2.3 Process-based Simulation of Water Erosion (SIMWE)

The methodological framework for the simulation of erosion/deposition processes within the SIMWE model (SIMulation of Water Erosion) is based upon the description of water flow and sediment transport processes by first principles equations, a concept outlined previously, most often for a one dimensional case, for example by Foster and Meyer (1972) or Bennet (1974). Within our approach, inputs and outputs of the simulations are represented by multi-variate functions (scalar or vector fields) rather than homogeneous hillslope segments. The underlying continuity equations are solved by Green's function Monte Carlo method, to provide robustness necessary for spatially variable conditions and high resolutions. Detailed description of equations used in this model is given by Mitas and Mitasova 1998.

2.3.1 Formulation

Overland water flow A 2D shallow water flow is described by the bivariate form of Saint Venant equations (e.g., Julien et al., 1995). The continuity of water flow relation is coupled with the momentum conservation equation and for a shallow water overland flow, the hydraulic radius is approximated by the normal flow depth The system of equations is closed using the Manning's relation. In this project, we assume that the solution of continuity and momentum equations for a steady state, provides an adequate estimate of overland flow for the land management applications (Flanagan and Nearing, 1995). In addition, we assume that the flow is close to the kinematic wave approximation, but we include a diffusion-like term to include the impact of diffusive wave effects. Such an incorporation of diffusion in the water flow simulation is not new and a similar term has been obtained in derivations of diffusion-advection equations for overland flow, e.g., by Lettenmeier and Wood, (1992). In our reformulation, we simplify the diffusion coefficient to a constant and we use a modified diffusion term. The diffusion constant which we have used is rather small (approximately one order of magnitude smaller than the reciprocal Manning's coefficient) and therefore the resulting flow is close to the kinematic regime. However, the diffusion term improves the kinematic solution, by overcoming small shallow pits common in digital elevation models (DEM) and by smoothing out the flow over slope discontinuities or abrupt changes in Manning's coefficient (e.g., due to a road, or other anthropogenic changes in elevations or cover).

Sediment flow. The basic relationship describing the sediment transport by overland flow is the continuity of sediment mass, which relates the change in sediment storage over time, and the change in sediment flow rate along the hillslope to effective sources and sinks (e.g., Haan et al., 1994; Govindaraju and Kavvas, 1991; Foster and Meyer, 1972; Bennet, 1974). The sediment flow rate is a function of water flow and sediment concentration. For shallow, gradually varied flow the storage term can be neglected leading to a steady state form of the continuity equation. The sources/sinks term is derived from the assumption that the erosion and deposition rates D(r) are proportional to the difference between the sediment transport capacity T(r) and the actual sediment flow rate |qs(r)|(Foster and Meyer, 1972)

D(r) = c(r) [T(r) - |qs(r)|](7)

with the first order reaction term c(r) dependent on soil and cover properties. The c(r) is obtained from the relationship (Foster and Meyer, 1972):

D(r)/Dm(r) + |qs(r)|/T(r) = 1 (8)

which states that the ratio of erosion rate D(r) to the detachment capacity Dm(r) plus the ratio of the sediment flow to the sediment transport capacity is a conserved quantity (unity). The sediment transport capacity and detachment capacity represent maximum potential sediment flow rate and maximum potential detachment rate, respectively, and are functions of a shear stress t(r) (Foster and Meyer, 1972)

t(r) = w h (r) sin b(r) (9)

where w is the hydrostatic pressure of water with the unit height. Then the simplified equation for transport capacity is

T(r) = Kt(r) t(r)p(10)


where Kt(r) is the effective soil transportability coefficient and p is exponent dependent on the type of flow. Detachment capacity estimated as a function of shear stress can be expressed as

Dm(r) = Kd(r) [t(r)-ts(r)] q (11)

where Kd(r) is the effective erodibility (detachment capacity coefficient), q is exponent and ts(r) is the critical shear stress. The parameters and adjustment factors for the estimation of detachment and transport capacity are functions of soil and cover properties, and their values for a wide range of soils, cover, agricultural and erosion prevention practices are being developed within the WEPP model (Flanagan and Nearing, 1995).

To solve the bivariate continuity equation we introduce a small diffusion term which represents local dispersion processes of the suspended flow (Bennet, 1974), for example, the impact of small, local slope changes below the DEM resolution, we rewrite the continuity equation with this term On the left hand side of the equation, the first term describes local diffusion, the second term is a drift driven by the water flow while the third term represents a velocity dependent 'potential' . The size of the diffusion constant used was again about one order of magnitude smaller than the reciprocal Manning's constant so that the impact of the diffusion term was relatively small.

Green's function stochastic method of solution. The equations can be solved, e.g., by projection methods (Rouhi and Wright, 1995). Another equivalent alternative is to interpret equations as a representation of stochastic processes with diffusion and drift components (Fokker-Planck equations) and carry out the actual simulation of the underlying process utilizing stochastic methods (Gardiner, 1985). This is very similar to Monte Carlo methods in computational fluid dynamics or to quantum Monte Carlo approaches for solving the Schrodinger equation (Schmidt and Ceperley, 1992, Hammond et al., 1994; Mitas, 1996). The solution by stochastic approach is illustrated by an animation in Mitas et al., (1997).

The Monte Carlo technique has several unique advantages which are becoming even more important due to new developments in computer technology. Perhaps one of the most significant Monte Carlo properties is robustness which enables us to solve the equations for complex cases, such as discontinuities in the coefficients of differential operators (in our case, abrupt slope or cover changes, etc). Also, rough solutions can be estimated rather quickly, which allows us to carry out preliminary quantitative studies or to rapidly extract qualitative trends by parameter scans. In addition, the stochastic methods are tailored to the new generation of computers as they provide scalability from a single workstation to large parallel machines due to the independence of sampling points. Therefore, the methods are useful both for everyday exploratory work using a desktop computer and for large, cutting-edge applications using high performance computing.

2.3.2 Properties

First we analyze the role of the SIMWE model parameters for complex terrain with spatially uniform rainfall excess, soil and cover properties.

First-order reaction coefficient c(r) is related to the soil detachability and transportability and controls the spatial extent of deposition. Depending on this coefficient, there are two limiting cases of erosion and sediment transport (Foster and Meyer, 1972; Hairsine and Rose, 1992): i) detachment limited and ii) sediment transport capacity limited. The first case is represented by c(r) close to 0, which results in the net erosion equal to the detachment capacity. Therefore, for conditions when c(r)<<1 (e.g., soils and cover with with substantially higher transportability than detachability), there is a prevailing detachment limited erosion, with almost all the detached sediment transported to the stream and with deposition restricted to small concave areas and channels. This case for a 1D hillslope is modeled by the Universal Soil Loss Equation (USLE). For the second limiting case represented by c(r) > 1, sediment flow is close to sediment transport capacity and the net erosion/deposition is modeled as a divergence of the sediment transport flow and the model predicts large extent of areas with deposition. Such a behavior is close to the observed distribution of colluvial deposits, suggesting the prevailing influence of the transport capacity limited case on a long term pattern of deposition. This case is modeled by the Unit Stream Power based erosion/deposition model (USPED) (Mitasova et al. 1996b). The results of simulations for the c(r) values increasing from 0.01 (fine soils) to 10.0 (sandy soils), assuming a rough cover (grass) with n=0.1 (Mitas and Mitasova 1998), demonstrate that with the changing c(r) the spatial distribution of erosion and deposition over the landscape changes dramatically. The erosion and sediment transport shows a transition between the two limiting cases (detachment and transport capacity limited). With the increasing c(r), the deposition areas expand while the sediment flow rate in the streams decreases until the pattern typical for a transport capacity limiting case is reached. Because the c(r) is dependent on the ratio between detachment capacity and sediment transport capacity both erodibility and transportability coefficients influence the spatial extent of deposition, however, as we demontstrate in the next section, each of these coefficients has different impact on the sediment loads in streams.

Erodibility or detachment capacity coefficient Kd(r) is a measure of soil susceptibility to detachment by water flow (Flanagan and Nearing, 1995). It is often defined as the increase in soil detachment per unit increase in shear stress of clear water flow. Change in erodibility factor while keeping the other parameters constant changes the ratio between transport capacity and detachment capacity . This leads to the change in the character of erosion process from detachment capacity limited for erodibility significantly lower than transportability to transport capacity limited case when erodibility is greater than transportability. Therefore, as shown in the FIGURE 6, change in erodibility will change the spatial pattern of erosion and deposition, however the impact on the magnitude of the sediment load in stream is small. This example is a good illustration of the fact that measuring sediment load in the watershed outlet does not provide enough information for understanding the erosion processes on watershed's hillslopes and quite different processes can result in the same sediment concentration levels. This makes the use of instream sediment concentrations as a measure reflecting the impact of erosion protection measurers needed for land use management problematic.

Transportability or transport capacity coefficient Kt(r) is a measure of soil capability to be transported by water flow. It depends on soil properties but can be also influenced by vegetation. This coefficient is not directly measured and provided in the WEPP, rather it is estimated indirectly, which is making the proper determination of this parameter problematic. However, the parameter can be derived at least for some types of soils using the published values of first order reaction coefficient or using the procedure suggested by Finkner et al. (1989). Our simulations show that this parameter has a profound impact on the erosion process as it influences both spatial distribution and magnitude of sediment flow and erosion/deposition rates. Recently, the importance of transport capacity for overland flow erosion processes has been fully recognized and more experimental and theoretical work is being performed (Guy et al. 1991, Govers 1991, Nearing et al. 1997). The FIGURE 6 illustrates how the change in transportability changes the erosion regime from detachment capacity limited to transport capacity limited while also reducing the magnitude of erosion rates.

Surface roughness is represented by the Manning's n(r) and has significant impact on the location and magnitude of deposition. For the same value of c = 1.0, the extent of deposition for smooth surfaces (e.g., bare soil with n=0.01) is smaller than for rough surfaces (e.g., grass with n=0.1), assuming also an increase in detachability and transportability for bare soil as compared to grass.

Rainfall excess i(r) is estimated as rainfall intensity - infiltration rate where infiltration rate can be estimated from a separate module or using map calculator based on the soil data. Rainfall excess influences the magnitude of erosion/deposition rates, with increasing rainfall excess both the erosion and deposition rates increase, however the spatial pattern of erosion and deposition does not change.

Critical shear stress ts(r) represents the soils resistance to the shearing forces of water flow. It depends on soil and cover properties and the values are available e.g., from the WEPP manual (Flanagan and Nearing, 1995). If shear stress at the given location is lower than the critical shear stress, no soil is detached. This parameter therefore has an impact on erosion/deposition pattern. Its high value reduces the spatial extent of erosion, on the other hand, it can also increase the magnitude of erosion rates on steeper hillslopes and in areas with concentrated flow due to the fact that clean water has higher potential to transport the sediment.

It is important to note that the parameters do not act independently, they are interrelated and it is their interaction which controls the pattern and magnitude of erosion. For example, the growth of vegetation reduces both Kd, Kt and 1/n and the resulting erosion/deposition pattern depends on the interaction between the rates of this change. If both Kd and Kt change at the same rate, the spatial distribution of erosion/deposition stays the same and only the magnitude of rates changes. If vegetation growth reduces Kd faster than Kt, the erosion/deposition pattern will change from transport capacity limited towards detachment limited. The interrelation between the parameters is an open research question and there is a lack of systematic experimental and theoretical/modeling work in this area

The presented analysis shows that for the uniform soil and cover there is a basic pattern of erosion and deposition in complex terrain which does not change significantly even if rainfall, soil and cover changes. High erosion risk areas are located on upper convex parts of hillslopes and in hollows and centers of valleys with concentrated flow. Deposition occurres in concave valleys and concave lower parts of hillslopes. With changing soil and cover conditions the line of starting deposition can move up- or down-slope, potentially influencing the sediment loads in streams. The analysis demonstrates that besides the obvious fact that the worst situation occurs for large events with high rainfall intensity when the soil surface is rather smooth (e.g. without a protective vegetation cover), sediment loads can be very high even for smaller events if soil is highly transportable. Change in the soil and cover from a uniform to spatially variable distribution has a great impact on the basic pattern of erosion/deposition (as we show later). Depending on the location this spatial variability can positively or negatively influence the overall soil loss as well as the sediment loads in streams.

While terrain plays an important role in spatial distribution of erosion/deposition, and in general we find deposition in concave areas and erosion in convex areas, this pattern can change significantly due to the spatial variability of land cover and soil properties. Borders between different land covers (e.g., bare soil and dense grass) cause abrupt changes in flow velocities, as well as in transport and detachment capacities, creating effects important for erosion prevention. We illustrate these effects for the original, conventional land use design in a subset of the study area, which combines arable land and meadows (FIGURE 7). We performed a simulation with the parameters set for dense grass in the meadow area (n=0.1, Kt = Kd = 0.0003) and bare soil in the arable area (n=0.01, Kt = Kd = 0.03), and compared the results with the observed spatial patterns of erosion/deposition.

The results of the simulation show several important phenomena. The highest rates of both the net erosion and the net deposition are predicted in valleys with high concentrated sediment flow (FIGURE 7), a phenomenon typical for thalweg erosion. Field measurements confirm that this area has the thickest layers of colluvial deposits with large linear erosion features observed after an unusually strong storm (FIGURE 7), (Auerswald et al., 1996). About a magnitude lower, but still relatively high erosion rates were predicted on bare upper parts of hillslopes which had the highest loss of radio-tracers and the lowest yields and where initiation of dense rilling occurred (FIGURE 8). Another area of increased erosion is predicted for narrow stripes below the grass areas, where water accelerates after depositing the sediment and leaving the grass. This leads to an increase in the difference between the sediment transport capacity and the actual sediment flow, creating the conditions for higher net erosion.

Deposition is predicted at the lower, concave parts of hillslopes which is in agreement with the spatial distributions of colluvial deposits and radio-tracers (FIGURES 7, 8). The sediment flow rate sharply decreases and eroded material deposits at the upper edges of meadows, explaining the cut-off rills with unusual shapes, which follow the borderline between the grass and bare soil areas (FIGURE 7). The influence of grass cover on flow velocity, transport and detachment capacities, also explains the observed location of deposited material on the convex, upper part of the hillslope with a meadow, however, to properly simulate the impact on the depth of colluvial deposits, long term data on rainfall and especially land cover would be needed.

2.4 Role of visualization in simulations

Scientific visualization is used as both a process of research and discovery and a method of communicating measured or modeled geographic phenomena. As a process of discovery, the visualization process is cyclical in nature, with visualizations feeding a refinement of the model. As a method of communication, visualization is used to demonstrate complicated processes through the use of images and animations. Combinations of various graphical representations of raster, vector, and point data displayed simultaneously allow researchers to study spatial relations in 3D space. At the same time, visual analysis of data requires the capability to distort this spatial relation by changing vertical scale, separating surfaces, performing simple transformations on point or vector data for scenario development, etc.

Visual analysis of multi-variate data represented by multi-dimensional rasters is also supported by capabilities to create fly-by's and to animate series of surfaces and isosurfaces. Vertical relationships between multiple surfaces, typical for representation of soil data, are analyzed by moving cutting planes (FIGURE 9). Surface/hypersurface representation is not restricted to globally continuous fields. A method for effective representation and visualization of surfaces with discontinuities using splines and GIS tools was developed by Cox et al. 1994. Phenomena represented by classes, such as vegetation or land use, can be represented as surfaces with faults as well. Proper interpolation to high resolution 2D,3D and 4D rasters, adequate scaling in 3rd and 4th dimension, selection of proper color schemes and visualization parameters play a key role in creating dynamic cartographic models from land characterization data. The efficiency and suitability of visualization tools for exploring multi-variate land characterization data is ensured by a high level of interactivity and by a combination of advanced visualization capabilities with the traditional spatial query and analysis functions of GIS (Brown et al. 1995, Brown and Astley 1995).

Multiple surfaces have proven useful to visualize boundaries of layers when examining soil core data and soil horizon crossections. Achieving intuitive visual representations of these horizons presents a technical challenge in terms of dimensional scale, as demonstrated by FIGURE 9. The vertical dimension is often quite small relative to the other two (eg, soil depth of a few meters vs. region dimensions of several kilometers), so is often exaggerated when a single surface is displayed. This exaggeration is usually adequate to add relief to an otherwise featureless surface, but in order to separate close stratified layers, the required exaggeration grossly distorts the modeled layers. If vertical translation of a surface is used to separate surfaces enough that they may be viewed separately, intersections between horizons and relative distances are misrepresented. This is unacceptable since these may actually be the features we are interested in viewing. To study differences between two similar surfaces, we use a scaled difference approach where only the spatial distance between surfaces is exaggerated, maintaining correct surface intersections.

A post-processed terrain surface rendered in the same 3-space as the original data surface help to reveal the successes and failures of the model; localized overshoots and undershoots, excess smoothing of significant detail or revelation of detail hidden by noise is easily compared in a qualitative manner. (Fig?)

Animation is an important tool for exploring large and complex data sets, especially when they involve both spatial and temporal dimensions. Recent progress in graphics technology and emerging standards for animation file formats have made desktop animations easier to produce and share with colleagues. Animation is useful for representation of change in time, change in a modeling parameter (Ehlschlaeger 1994), change in viewer position (fly-bys) or change in visible data (fence cuts, slices). FIGURE 10 shows several frames from an animation where a fence cut is moved through data to better view underlying surface structure. We implemented animation capabilities in both 2D and 3D tools.

Querying a 2D data set displayed as a raster image can be thought of as a scale operation and a translation operation. When a user clicks on a pixel, the relative position in the image is scaled by the resolution of a pixel and the north and east offsets are added to obtain the geographical position. When displaying surfaces in 3D space with perspective, however, clicking on a pixel on the image of the surface really represents a ray through 3D space. The point being queried is the intersection of this ray with the closest visible (unmasked) part of one of the surfaces in the display. One method for 3D querying provided by some graphics libraries is a "selection" method, where objects are "redrawn" without actually drawing to the screen, and any objects drawn at the query point are returned as the "selected" objects. This method is slow and at best the returned object is a polygon rather than a specific point on the polygon. Therefore, we require the user to specify the type of data they are querying (surface, point, or vector) and then use our knowledge of the geometry of that data to perform a geometric query in 3D. FIGURE 11 shows that using cutting planes can allow the user to query a specific location on a surface that is covered by another surface. This specialized point-on-surface algorithm can be outlined as follows: 1. Transform point on view plane to a view ray. 2. Intersect view ray with convex polyhedron defined by the intersection of the parallelepiped view region with any user defined cutting planes. 3. If ray enters this polyhedron, trace ray to find any intersections with visible (unmasked) parts of any surfaces. 4. Choose closest intersection to viewer (or return an ordered list). Such point-on-surface functionality is useful for 3D data querying, setting center of view and setting center of rotation for vector transformations.

GRASS, as an open system with a full range of GIS capabilities has provided a sound basis for testbed development of visualization tools. Each GRASS data type (raster, vector, and site) plus our own 3D grid format may be used for visualization in a single 3D space. In our implementation, there are four object types and various ways to represent each.

Surfaces. A surface requires at least one raster data set to represent topography and may use other raster data sets to represent attributes of color, transparency, masking, and shininess. These data sets may have been derived from vector (e.g., countour) or scattered point data using tools within the GIS. Users are allowed to use a constant value in the place of any raster data set to produce, for example, a flat surface for reference purposes or a constant color surface.

Vector sets. 3D vector sets are not currently supported, so in order to display 2D vector sets in 3 dimensions, they must be associated with one or more surfaces. The 2D vector sets are then draped over any of these surfaces.

Site Objects. Point data is represented by 3D solid glyphs. Attributes from the database may be used to define the color, size, and shape of the glyphs. 2D site data must be associated with one or more surfaces, and 3D site data may be associated with a surface (e.g., sampling sites measured as depth below terrain).

Volumes. 3D grid data may be represented by isosurfaces or cross sections at user-defined intervals. Color of the isosurfaces may be determined by threshold value or by draping color representing a second 3D grid data set over the isosurfaces. Implementation and initial testing of a 3D grid data file format for managing volumetric spatial data was completed. The storage format and programmer's interface routines we developed allows random access to compressed floating point double precision 3D data with caching. It is fully integrated within the GIS, using the established database hierarchy for header and data files. See the html documentation for the specification used. The resulting library was used to write several utility applications such as r3.in.ascii, r3.out.ascii, r3.in.grid3, r3.mask, r3.null, r3.info, g3.region . In addition, a program r3.mkdspf , reads the 3D grid data and creates a "display" file containing geometry for drawing isosurfaces to represent the data for visualization purposes. Future work will incorporate the 3D grid file and display file formats for use within visualization tools. Appendix 2 contains the specification and design documents for our 3D grid API, and the function prototypes.
http://www.cecer.army.mil/grass/viz/htdoc/g3d/specification.html
http://www.cecer.army.mil/grass/viz/htdoc/g3d/protos.html

For interactive viewer positioning, scaling, zooming, etc., we use custom GUI widgets and a "fast display mode" where only wire mesh representations of the data are drawn. When rendering a scene, the user may select various preset resolutions for better control over rendering time. For positioning we also chose to use a paradigm of moving viewer rather than moving object because we think it is more intuitive when modeling a reality of generally immobile geography. To focus on a particular object, the user simply clicks on the object to set a new center of view. Scripting is used to create animations from series of data (e.g. time series or a changing modeling parameter), automatically loading the data sets and rendering frames for the animation. A keyframe technique is used to generate animations when there is no change in data, e.g., to create fly-bys or show a series of isosurfaces in volumetric data. (Also see nviz Tutorial in Appendix )


3. Applications and results

The methods and algorithms developed for this project were applied and tested for several different sites. High resolution data from Scheyern experimetal farm (Technical University Munchen) were used for testing, evaluating and analysis of the developed algorithms and tools. Applicability of the proposed methdology to military land management was tested at various military installations using different types of input data.

3.1 Erosion prevention design for experimental farm

The properties and functioning of the developed methods for interpolation, topographic analysis and erosion modeling were illustrated in previous chapters using the data from Scheyern experimetal farm (Technical University Munchen). We have also used these data to investigate the possibility to use the SIMWE model for analyzing and designing the placement of selected erosion protection measures based on land cover, as illustrated by the following example.

We have applied SIMWE to the experimental farm for the traditional land use management scenario (FIGURE 7A). This land use has resulted in severe erosion when the large storm event occurred during the time when the agricultural fields were bare. The field data showed that the areas with the highest water and sediment flow potential were unprotected and significant rilling and gully formation occurred. Simulations performed for these conditions confirm high potential for erosion on convex parts of the bare hillslopes and high sediment flows in the areas of concentrated water flow. We have then used the SIMWE model to investigate the possibilities to redesign the land use so that the net soil loss as well as sediment loads are minimized.

First, we used the model to identify locations with the highest erosion risk, assuming a uniform land use. Then, the protective grass cover was distributed to the high risk areas while preserving the extent of grass cover at the original 30% of the area (FIGURE 7C). We performed a simulation with the new land use to evaluate its effectiveness. The results demonstrate that the new design has a potential to dramatically reduce soil loss and sediment loads in the ephemeral streams (FIGURE 7C). The crest in sediment flow in the valley disappears and is replaced by light deposition within the grassways, while the maximum and total rates of erosion are significantly reduced (Mitas and Mitasova 1998). We have found, that the effectiveness of this design depends on differences in roughness, combination of very smooth bare soil and very dense grassway resulted in predictions of higher erosion along the borders of the grassway - more study is needed to decide whether this is an artifact of the method or a realistic prediction.

It is interesting to note, that the land use design obtained by this rather simple computational procedure, using only the elevation data, had several common features with the sustainable land use design proposed and implemented in 1993 at the farm, based on extensive experimental work (Auerswald et al., 1996). A simplified version of that landuse design, along with prediction of water and sediment flows and net erosion/deposition is presented in FIGURE 7B. It uses significantly higher proportion of permanent grass cover and fallow and increases roughness in the hop field. The results show that this design keeps lots of moisture and, at a higher cost, reduces the erosion even further.

3.2 Topographic potential for erosion at Yakima

In the early stages of the project, we have tested the possibilities to use a standard USGS 30m DEM product to predict topographic potential for erosion/deposition (Mitasova et al. 1996). We have found that the horizontal and vertical resolution, 30 m and 1 m respectively were insufficient for the analysis of water flow and erosion/deposition modeling. To get reasonable results from these data, we reinterpolated the elevations to the new DEM with 10m horizontal and 0.01m vertical resolution using regularized spline with tension, with simultaneous smoothing and computation of slope and aspect. Detailed analysis of the original and smoothed DEM, using 3 dimensional shaded views, curvatures and histograms, revealed that two levels of systematic errors were significantly reduced - dense strips in the direction parallel to y-axis (see Mitasova et al. 1995) and less frequent strips in the north-western direction. The strong waves along the 20m contour values, present in the original DEM, were also reduced but profile curvature and histogram analysis proved that it was not possible to remove the waves completely. Resampling and smoothing has improved the results for the upslope contributing areas (FIGURE 12) and the detachment capacity index (FIGURE 12). However, the erosion/deposition model still predicted waves of deposition along 20m contours (FIGURE 13a). To confirm that these waves are artificial, we have extracted 10m contours from the original grid DEM and interpolated a new DEM from these contours using our spline function with properly set smoothing and tension. The new DEM has r.m.s deviation 2.3m which is well within the the vertical accuracy of data given as 7m and the bounding effect around each of the given contours is not present. The erosion/deposition model predicts the highest erosion in strongly convergent areas indicating creation of channels, and on the steep upper part of hillslopes. Several depositional areas are predicted in concave parts of hillslopes, mostly with dispersal flow. Comparison of the results from the original, smoothed and recomputed DEMs (FIGURE 13) demonstrates that the accuracy of USGS DEM was not sufficient for erosion/deposition modeling and more accurate data are needed to obtain reliable results. Therefore, for the entire Yakima range, we have computed only topographic potential for detachment limited erosion (FIGURE 14) which does not require the computation of deposition and is thus less sensitive to artifacts in DEM.

3.3 Analysis of topographic potential for erosion/deposition at Ft. McCoy

Topographic potential for soil detachment and for net erosion and deposition at Fort McCoy was estimated by the modified USLE and USPED models (Mitasova et al. 1996). Even with rather crude elevation data, it was possible to identify some areas with high topographic risk for soil detachment and erosion/deposition, as illustrated by the following sequence of maps showing the input 30m DEM (FIGURE 15) and results of computations (FIGURE 16). The plateaus in low elevation areas due to the insufficient vertical resolution of 1m create an artificial pattern of steeper slopes along 1m contours and upslope contributing areas computed from the original 30m DEM by r.flow are underestimated. Pits and plateaus in the DEM cause problems in flowtracing with almost no flow generated in lower elevation areas and disrupted flow in valleys. Modified LS factor predicts well some areas of high potential for soil detachment, especially in hilly parts of the region. However, it does not predict the location of areas with concentrated flow with high potential for gully formation . Also, in lower elevation areas the artificial pattern of increased erosion is predicted along 1m contours. The net erosion/deposition map computed by the USPED model predicts high net erosion in hilly regions and on slopes along the streams. It also shows that a significant portion of the material eroded from hillslopes is deposited before it could reach the main streams. However, this map lacks the prediction of high erosion due to concentrated flow in valleys and shows artificial waves of erosion and deposition along the 1m contours in flatter areas.

To reduce the negative impact of lower resolution of the available elevation data the given DEM was reinterpolated to 10m horizontal and 0.01m vertical resolution using the RST method. Simultaneously with interpolation, maps of slope, aspect, profile and tangential curvatures were generated. Steady state water flow computed by r.flow from the smoothed 10m DEM shows potential for channel formation in valleys and predicts water flow also in low elevation region (FIGURE 17). Flow in the two main streams is not adequately described because of lack of data, but it can be incorporated if stream data are available. Modified LS factors were computed from the 10m DEM using different exponents for the water flow term (p=0.6 and p=1.5). Value of this exponent is still a subject of research and discussion in erosion research community. The first result for p=0.6 puts more weight on the influence of slope and theoretically represents the detachment limited erosion (detachment capacity of overland flow). The second result with p=1.5 puts more weight on the influence of water flow and theoretically represents the sediment transport capacity of overland flow (sediment flow for transport capacity limited case). Because the exponent depends on the conditions of water flow in a particular area it should be calibrated to reflect the type of flow typical for the modeled area and time during the year. Topographic potential for net erosion/deposition was computed from the 10m DEM using the USPED model. Similarly as for the result of 30m DEM, model shows high erosion in hilly area and along main streams and deposition in concave areas. It also indicates location of areas with concentrated flow which can reach streams.

3.4 Impact of land use change for a proposed Camp Shelby extension

We have used the modified USLE and the USPED model to evaluate spatial and temporal distribution of erosion and deposition potential at a planned military installation extension undergoing significant changes in land use during the year. The high resolution (5m) DEM with slope and aspect was created from digitized contours using the regularized spline with tension. The upslope contributing area for each grid cell was computed by the vector-grid based flow tracing algorithm. To estimate the spatio-temporal change in detachment and transport capacity limited erosion due to the change in rainfall and cover we have computed the modified USLE and USPED models using the raster map calculator r.mapcalc (Shapiro and Westervelt 1992). The series of raster maps representing the spatial and temporal distribution of detachment capacity index D(r) with the influence of rainfall, soils and landuse during the year was then computed as

Di(r,ti) = R(ti) . C(r, ti) . K(r) . LS(r),         i=1,...,12 (12)

where R is the rainfall factor uniform over the area, but significantly changing during the year, C is the cover factor, which varies in space and time, K is the soil factor and LS is the modified topographic factor (FIGURE 18). The RUSLE factors were used here as weights, to model the relative increase/decrease in detachment due to the changes in rainfall and cover. The net erosion and deposition index was computed using the USPED model for the current undisturbed landuse and for several land use alternatives with a significant portion of the area exposed to intensive landuse during training (FIGURE 19).

The study has shown that both spatial and temporal distribution of erosion are important for land use management to minimize the soil loss. Spatial analysis allowed us to evaluate various land use alternatives and assess the impact of planned extension of training areas. Clearing the current forest and training without any erosion proptection measures will lead to a 30-times increase in average soil loss per acre and about the same dramatic increase in net erosion which will be entering the streams. Combination of protective stream buffers and exclusion of steeper slopes (>10%) from training can reduce the potential for erosion by about 40%, however it will still be more than 10 times higher than before clearing (FIGURE 19).

Temporal analysis showed that the highest intensity of land use occurres during the time of the highest rainfall intensity, compounding thus the negative effect of rainfall and land cover disturbance (FIGURE 18). Rescheduling the timing of intensive use so that it does not coincide with the highest rainfall intensity is one of the possibilities to reduce the soil erosion risks for this area.

3.5 Multiscale analysis of topographic potential for erosion/deposition at Ft. Irwin

To illustrate the issues associated with simulations for large areas we use an example of a mountainous region in California. The standard 30m DEM available for the entire study area (3000 square km) represents 4 million grid cells, a challenging data set for the current process-based simulation tools and workstations, at the resolution hardly sufficient for rough identification of high erosion risk areas. The DEM at 5m resolution, needed e.g. for capturing at least approximately the effects of roads or grassways, would require to use data sets with 121 million grid cells and run simulations which would be prohibitively expensive, if not impossible with the current computational resources. It is clear, that for such a large area, modeling at different scales and resolutions is needed, depending on the importance and complexity of area's watersheds. It is important to note that our aim in this application is to illustrate the possibilities to use standard elevation data for erosion simulations in a large area. We are not taking into account the climate, soil and cover properties which certainly have a profound impact on water flow in this area. In spite of the fact that this is only a hypothetical example, this test allows us to assess the capabilities of our tools and prepare the environment for more realistic applications (e.g. Hohenfels installation).

To demonstrate the impact of resolution, noise and systematic errors in standard elevation data we have computed topographic parameters from the 30m DEM available for Ft. Irwin, with more detailed analysis performed for the subareas A, B, C, D (FIGURE 20). Topographic parameters are useful for evaluating the quality of a DEM and for identifying possible noise and systematic errors, as illustrated by the following example of terrain with draped tangential curvature. Tangential curvature for 30m DEM shows acceptable structure in mountainous area while significant noise and systematic errors (stripes) are present in lowland (FIGURE 21a). After smoothing and resampling to 10m resolution using the RST interpolation method the noise is reduced and the major topographic features become more visible (FIGURE 21b). The analysis clearly demonstrates that the need for precision and accuracy in elevation data is spatially variable. Areas with flat terrain are much more sensitive to noise and systematic errors in elevation data than mountains. Moreover, the artificial structure in flat areas continuously transforms into the real terrain structure in mountains creating a risk that, if this DEM is going to be used for simulations, the artificial structure can be mistaken for the real topographic feature. Smoothing and resampling the original 30m data to 10m resolution has also improved flowtracing by reducing the artificial pits in the original DEM leading to creation of continuous stream network for a test subarea C.

Topographic potential for detachment limited erosion and transport capacity limited erosion/deposition was estimated by the modified USLE and the USPED models. To illustrate the impact of smoothing and resampling on erosion modeling we have computed detachment capacity of water flow (modified USLE LS-factor, FIGURE 22) and topographic potential for net erosion/deposition (FIGURE 23) at three resolutions: a) downscaled to 90m for the entire region, b) original 30m for a 500 square km subregion, c) smoothed and resampled to 10m for a 150 square km subregion. To demonstrate the differences in detail which can be achieved at various resolutions we have visualized the results as color maps draped over the 10m resolution DEM for a small subregion (36 square km) (FIGURES 24, 25). The results show that the 90m resolution is inadequate for capturing the erosion by concentrated flow and does not allow to predict a realistic net erosion/deposition pattern. The 30m resolution starts to reveal the structure of erosion pattern, including the concentrated flow erosion in the valleys, however it is not sufficient for predicting net erosion/deposition pattern. The 10m resolution does not improve the accuracy of the original elevation model, however, the smoothing reduces the noise in the elevation data and higher resolution allows better description of terrain geometry leading eventually to more realistic results both for the detachment limited erosion and net erosion and deposition (FIGURES 24, 25).

Water and sediment flow simulation by SIMWE. To test the applicability of the SIMWE model to areas significantly larger than the test Germany data we have computed the spatial distribution of steady state water flow, sediment flow and net erosion/deposition for a 36 square km subregion D using the reinterpolated and smoothed 10m DEM (FIGURE 26). The water flow simulation captured the complex pattern of overland water flow, including the dispersal flow features typical for this area, which are difficult to predict by more traditional approaches. Sediment flow rates were estimated by solution of continuity of mass equation predicting high sediment flow rates in centers of valleys with concentrated flow and dispersal of sediment flow with reduction of sediment loads in areas with alluvial cones (FIGURE 26a). Net erosion/deposition rates were estimated for prevailing transport capacity limiting regime due to the prevailing sandy soils (FIGURE 26b). The simulation shows the formation of split gullies and their disapearance as the terrain flattens and the alluvial cones are formed.

This application demonstrates that processing of original elevation data by adequate tools, such as in our case the RST method and application of a robust, process-based model can lead to a more realistic prediction of a complex erosion/deposition pattern. This application also proved that the SIMWE model can be applied to large areas, in spite of its relatively complex formulation and computational demands.

The resampled and smoothed DEM and the results of simulation by the SIMWE model were provided for a development and demonstration of virtual reality tools for hydrologic simulations for CAVE at NCSA (Johnston and Reez, 1997) as illustrated by FIGURE 27. http://www.gis.uiuc.edu/hpgis/visualization.htm.

3.6 Simulation of erosion/deposition at Hohenfels

Cost effective land maintenance and rehabilitation at military installations requires identification of critical areas which should be targeted for erosion prevention. High resolution, process-based simulations provide valuable tool for assessing the current situation, evaluating the impact of high erosion and cost effective planning of implementation of rehabilitation and mitigation measures. Using the Hohenfels installation as an example, we demonstrate the capabilities of different models developed for this project to provide maps supporting assessment of erosion risk and identification of priorities in implementation of erosion prevention measures.

Topographic analysis As we have demonstrated by the previous examples, computation of a DEM and topographic analysis are crucial components of data preparation for erosion modeling. For the Hohenfels installation, we have evaluated various sources of elevation data and we have found that the most suitable DEM can be created by interpolating from 5m contours using the RST method. This approach has allowed us to perform simultaneous topographic analysis including the computation of slope, aspect and curvatures. For a watershed in eastern part of the installation we have computed a 2m resolution DEM from the 5m contours and for a small subarea of this watershed we have also computed a 2m resolution DEM from 1m contours for multiscale simulations (FIGURE 28).

Application of the modified USLE and the USPED model We have performed modeling of detachment limited erosion and transport capacity limited erosion and deposition for the current land use, represented by a C-factor based on 1997 land cover data (FIGURE 29). The analysis was computed for the entire installation at 10m resolution and for an eastern part of installation at 2m resolution, using the modified USLE and the USPED models. We have assumed uniform rainfall and we have used the K-factor based on 1997 soil data.

Results show that compared to the USLE estimate from previous study (Warren and Kowalski 1989) which predicted that 14.6 % of total area has potential for severe erosion areas, only 6.6% area is at high erosion risk (FIGURE 30), however, the rates of erosion and sediment transport in these areas can be much higher than indicated in the previous study (Warren and Kowalski 1989) For comparison, we have also used modified USLE with m=1, however this exponent leads to prediction of extremely high erosion rates and of large spatial extent of erosion risk at 22% area. Field observation and/or remote sensing data would be useful for validation of these results and for assessing the validity of our hypothesis that the exponent m=1 leads to overprediction of potential erosion. Net erosion/deposition modeled by USPED (with transport capacity computed similarly as the detachment capacity in the previous case, but with m=1.6) predicts that only 6.6% of total area has potential for severe net erosion and a significant portion of the eroded material is deposited before it enters the streams.

To increase the detail of the predicted erosion/deposition pattern we have used a 2m DEM computed for an eastern part of the installation from the same 5m contours and we have performed the same analysis of erosion and deposition patterns as in the previous case, but at a higher resolution. The results are consistent with the assessment for the entire installation, with the maps providing more detailed information on location of high risk areas (FIGURE 31).

Application of the SIMWE model. Because the SIMWE model is computationally more demanding than USPED we have performed the simulations at 10m resolution for the eastern watershed, which required to run the simulations on a 400x400 grid. For comparison, the 2m resolution modeling by USLE/USPED for this area involved processing of 2000x2000 grid. We have performed the simulations for the current land use for 2 different conditions (a) fully saturated soil with all rainfall contributing to surface runoff (b) spatially variable rainfall excess with high portion of rainfall infiltrated in the undisturbed areas covered mostly by forest and with reduced infiltration in the disturbed areas. The condition (b) involves modeling with spatially variable rainfall excess which we haven't tested in our previous simulations. The input parameters (elevation gradient, rainfall excess, erodibility, transportability, roughness-Manning's n, critical shear stress) were estimated using the elevation, C-factor and K-factor data, as well as the information about the land cover and soils from the report by Warren and Kowalski, 1989.

SIMWE results are consistent with the USPED model, however there are some important differences. SIMWE predicts greater risk for erosion in hollows and centers of valleys and on top/upper parts of disturbed ridges. The area with severe erosion risk is about..... Spatially variable infiltration increases the difference between sediment loads from disturbed and undisturbed areas, with undisturbed areas being much more stable and producing very little sediment compared to disturbed areas where overland flow is much higher and faster due to the lack of infiltration on compacted soil. The new study shows several issues relevant to land management: (a) the importance of protection in areas of concentrated flow which were not identified in the previous study, (b) significant amount of eroded soil from hillslopes has a potential for being deposited before entering the streams.

The comparison of the results of erosion modeling from 10 years ago and using the current technology developed for this project demonstrated the advancements in distributed simulations technology. In spite of the fact that the spatial accuracy of the erosion risk assessment is limited by the resolution of cover data, (which is 20m, while the DEM is 10m and 2m resolution), the high resolution elevation data and new models allow us to predict erosion risk areas with greater detail and identify also the high risk areas with concentrated water flow which were missing in the results from the previous study. However, it is important to note that the models are not calibrated so the quantitative results should be used with caution.


4. Conclusion

Integration of topographic analysis and erosion modeling with a GIS provides an environment for effective evaluation of various approaches to erosion/deposition risk assessment for landscape scale applications. The analysis of the models developed for this project: Modified USLE, USPED and SIMWE demonstrates that the SIMWE model formulation and implementation provides the most general approach, capturing erosion, sediment transport and deposition for a wide range of conditions. However, it is also the most computationally intensive approach and requires new type of input data which are being provided for the WEPP model. Both modified USLE and USPED models can be derived as special, extreme cases of the SIMWE model, demonstrating thus that they represent erosion process only under special conditions. However, because of their simplicity both modified USLE and USPED are useful for erosion risk assessment, as they are quite consistent with the SIMWE in predicting the high erosion risk areas (although the erosion maxima predicted from USLE are slightly shifted downslope when compared to USPED.)

The comparison of the modified USLE and the USPED model demonstrated that the USPED approach is more appropriate for landscape scale erosion modeling in complex terrain, especially when the location of both areas with erosion risk and with deposition potential is important. We have enhanced the unit stream power based approach by improving the computation of DEM using the smoothing spline with tension, and by refining the flow tracing using the vector-grid algorithm. The originally proposed computation of the change in sediment transport capacity by finite difference technique was improved by using the divergence of the bivariate function representing the sediment transport capacity of flow. Several practical applications demonstrated that the use of the standard, widely available 30m DEM for erosion/deposition modeling is problematic and can be used only after significant reprocessing. We believe that more accurate data, such as the digitized contours from at least 1:24000 topographic maps field measured elevation data, or contours derived from digital orthophoto maps used along with reliable interpolation are necessary.

While several important phenomena were not studied in this work and are the subject of our current research, both the USPED and SIMWE models and the comparisons of their results with the observed data have demonstrated the important role of terrain in the development of erosion/deposition patterns in the landscape. We have also shown that spatially variable cover can significantly change the general pattern of erosion/deposition and that the capability to simulate the cover's impact can become a powerful tool for a computer aided design of cost effective erosion protection measures.

Our simulations have pointed out to some unsolved issues related to the interpretation and determination of the soil and cover parameters. Our effort to predict the observed extent of deposition indicates that the equation for estimating the sediment transport capacity may not be applicable to a wide range of situations and a more general equation suitable for overland flow over complex terrain is still needed, as suggested by several papers (e.g., Govers, 1991; Guy et al., 1991). It is also important to note that the parameters n, Kt and Kd incorporate the influence of numerous physical properties of soil and cover which are not mutually independent, however, their functional relationship is not known. The ongoing experimental research (e.g., Flanagan and Nearing, 1995), as well as the presented simulations, can provide a better insight into the physical meaning of these parameters and improve the quantitative accuracy of the predictions.

The proposed erosion simulations tools are supported by numerous GIS functions, such as interpolation and visualization, and by new data structures which had been significantly enhanced within this project. These functions have been developed and anhanced for the public domain GRASS GIS which is being relased as GRASS5.0 by Baylor University, Texas. We have also tested some commercial software tools for support of erosion modeling and while the necessary commands for USPED and modified USLE are available and the methodology developed for this project can be used, often the underlying algorithms are not sophisticated enough to produce the results comparable with the maps and images presented in this report.

5. References

Abramowitz, M., and Stegun, I.A., 1964, Handbook of Mathematical Functions: Dover, New York, p.297-300,228-231.

Arnold, J. G., P. M. Allen, and G. Bernhardt, 1993, A comprehensive surface-groundwater flow model. Journal of Hydrology, 142, 47-69.

Auerswald, K., A. Eicher, J. Filser, A. Kammerer, M. Kainz, R. Rackwitz, J. Schulein, H. Wommer, S. Weigland, and K. Weinfurtner, 1996, Development and implementation of soil conservation strategies for sustainable land use - the Scheyern project of the FAM, in Development and Implementation of Soil Conservation Strategies for Sustainable Land Use, edited by H. Stanjek, Int. Congress of ESSC, Tour Guide, II, pp. 25-68, Technische Universitat Munchen, Freising-Weihenstephan, Germany.

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

Bjorneberg, D. L., Aase, J.K., Trout, T.J., 1997, WEPP model erosion evaluation under furrow irrigation. Paper no. 972115, 1997 ASAE Ann. Int. Meeting, Minneapolis, MI.

Brown, W.M., and Gerdes, D.P., 1992, SG3d - supporting information: U.S. Army Corps of Engineers, Construction Engineering Research Laboratories, Champaign, Illinois. http://www.cecer.army.mil/grass/viz/sg3d.html

Brown, W.M., 1996, SG3d - supporting information. Enhacements for GRASS4.2: Geo- graphic Modeling and Systems Laboratory, University of Illinois at Urbana-Champaign.

Brown, W.M., and Astley, M., 1995, NVIZ tutorial:Geographic Modeling and Systems Laboratory, University of Illinois at Urbana-Champaign.

Brown, W. M., M. Astley, T. Baker, and H. Mitasova, 1995, GRASS as an integrated GIS and visualization environment for spatio-temporal modeling, in Proceedings of the Auto Carto 12, edited by D.J. Peuquet, pp. 89-99, ACSM/ ASPRS, Charlotte, NC.

Busacca, A. J., C. A. Cook, and D. J. Mulla, 1993, Comparing landscape scale estimation of soil erosion in the Palouse using Cs-137 and RUSLE, Journal of Soil and Water Conservation, 48, 361-67.

Cox, S., Kohn, B.P., and Gleadow, A.J.W., 1994, Towards a fission track tectonic image of Australia: model based interpolation in the Snowey Mountains using GIS: http://www.ned.dem.csiro.au/AGCRC/papers/snowys/12agc.cox.html

DeRoo, A.P.J., Hazelhoff, L., and Burrough, P.A., 1989, Soil erosion modelling using 'ANSWERS' and Geographical information systems. Earth Surface Processes and Landforms, 14, 517-532.

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

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

Dikau, R, 1990, The application of a digital relief model to landform analysis in geomorphology: in Raper, J. (Ed.), Three dimensional applications in Geographic Information Systems: Taylor \& Francis, London, 1989, p.51-77.

Ehlschlaeger, C., Goodchild, M. 1994, Uncertainty in Spatial Data: Defining, Visualizing, and Managing Data Errors. Proceedings, GIS/LIS'94, pp. 246-53, Phoenix AZ, October 39 1994.

Fairfield, J. and Leymarie, P., 1991, Drainage networks from grid digital elevation models. Water Resources Research, 27, 709-717.

Engel, B., Hydrology Models in GRASS, WWW, 1995.


http://soils.ecn.purdue.edu:80/~aggrass/models/hydrology.html

Fikner 1989

Flacke, W., Auerswald, K., and Neufang, L., 1990, Combining a modified USLE with a digital terrain model for computing high resolution maps of soil loss resulting from rain wash. Catena, 17, 383-397.

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.


http://soils.ecn.purdue.edu/~wepp/wepp.html

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., and Wischmeier, W.H.,1974, Evaluating irregular slopes for soil loss prediction. Transactions of the ASAE, 12, 305-309.

Foster, G.R., Flanagan, D.C., Nearing, M.A., Lane L.J., Risse L.M., Finkner, S.C., 1995, Hillslope erosion component. {\sl WEPP: USDA-Water Erosion Prediction Project} (Flanagan, Nearing eds.), NSERL report No. 10 National Soil Erosion Lab., USDA ARS, Laffayette, IN.

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.

Freeman, T.G., 1991, Calculating catchment area with divergent flow based on regular grid, Computers and Geosciences}, 17, 413-422.

Franke, R., 1982, Smooth interpolation of scattered data by local thin plate splines: Comput.Math.Applic., v.8, p.273-281.

GRASS4.1 Reference Manual, 1993, U.S.Army Corps of Engineers, Construction Engineering Research Laboratories, Champaign, Illinois, 422-425.

Gardiner, C. W., 1985, Handbook of Stochastic Methods for Physics, Chemistry, and the Natural Sciences, Springer, Berlin.

Govers, G., 1991, Rill erosion on arable land in central Belgium: rates, controls and predictability, Catena, 18, 133-155.

Govindaraju, R. S., and M. L. Kavvas, 1991, Modeling the erosion process over steep slopes: approximate analytical solutions, Journal of Hydrology, 127, 279-305.

Guy, B. T., Rudra, R. P. Rudra, and W. T. Dickinson, 1991, Process-oriented research on soil erosion and overland flow, Overland Flow Hydraulics and Mechanics, edited by A. J. Parsons, and A. D. Abrahams, Chapman and Hall, New York, pp. 225 -242.

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

Hairsine, P. B., and C. W. Rose, 1992, Modeling water erosion due to overland flow using physical principles 1. Sheet flow, Water Resources Research, 28(1), 237-243.

Hammond, B. L., W. A. Lester, Jr., and P. J. Reynolds, 1994, Monte Carlo Methods in Ab Initio Quantum Chemistry, World Scientific, Singapore.

Hibbard, W. L., and Santek, D.A., 1989, Visualizing large data sets in the earth sciences: IEEE Computer Graphics and Applications, v. 22, no. 8, p. 53-57.

Hibbard W.L., Paul B.E., Battaiola A.L., Santek D.A., Voidrot-Martinez M-F., and Dyer C.R., 1994, Interactive Visualization of Earth and Space Science Computations Computer v. 27, no. 7, p. 65-72.
href="http://www.ssec.wisc.edu/ billh/vis.html"

Hofierka, J., Suri, M., 1996, Soil Water Erosion Modelling Using GIS and Aerial Photographs. Proceedings of the JEC-GI 96 Conference Barcelona, Spain: 376-381.

Hong, S., and S. Mostaghimi, 1995, Evaluation of selected management practices for nonpoint source pollution control using a two-dimensional simulation model, ASAE, paper no. 952700. Summer meeting of the ASAE, Chicago, IL.

Howard, A. D., 1994, A detachment-limited model of drainage basin evolution, Water Resources Research, 30(7), 2261-2285.

Huang, S.L., and Ferng, J.J., 1990, Applied Land Classification for Surface Water Quality Management:II. Land Process Classification. Journal of Environmental Management, 31, 127-141.

Hutchinson, M.F., 1989, A new procedure for gridding elevation and stream line data with automatic removal of spurious pits. Journal of Hydrology, 106, 211-232.

Jenson, S.K., and Domingue, J.O., 1988, Extracting topographic structure from digital elevation data for geographic information system analysis. Photogrammetric Engineering and Remote Sensing, 54, 1593-1600.

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.

Kirkby, M.J., 1986, A two-dimensional simulation for slope and stream evolution, in Hillslope Processes, edited by A. D. Abrahams, pp. 203-222,Allen and Unwin, Winchester Mass.

Kramer, S., and M. Marder, 1992, Evolution of River Networks, Physical Review Letters, 68, 205-208.

Krcho, J., 1991, Georelief as a subsystem of landscape and the influence of morphometric parameters of georelief on spatial differentiation of landscape-ecological processes: Ecology /CSFR/, v.10, p.115-157.

Lettenmaier, D. P., and E. F. Wood, 1992, Hydrologic forecasting, in Handbook of Hydrology, edited by D. R. Maidment, pp. 26.1-26.30, McGraw-Hill, Inc., New York.

MacEachren, A.M., and Ganter, J.H., 1990, A pattern indentification approach to cartographic visualization: Cartographica, v. 27, no. 2, p. 64-81.

MacEachren, A.M., 1994, Time as cartographic variable: in Hearnshaw, H.M., and Unwin, D.J., ed., Visualization In geographical Information Systems: John Wiley and Sons, Chichester, p. 115-130.

Maidment, D. R., 1996, Environmental modeling within GIS, 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. 315-324.

Martz, L. W., and E. de Jong, 1987, Using Cesium-137 to assess the variability %of net soil erosion and its association with topography in a Canadian prairie landscape, Catena}, 14, 439-451.

Martz, L.W., and E. de Jong, 1991, Using Cesium-137 and landform classification %to develop a net soil erosion budget for a small Canadian prairie watershed, Catena, 18, 289-308.

McCauley J D, Engel B A 1995 Approximation of noisy bivariate traverse data for precision mapping.

Mitas, L., H. Mitasova, 1988, General variational approach to the interpolation problem. Computers and Mathematics with Applications 16, p.983-992.

Mitas, L., 1996, Electronic structure by Quantum Monte Carlo: atoms, molecules and solids, Computer Physics Communications, 97, 107-117.

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

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 (in press).

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 (includes CDROM and WWW: www.elsevier.nl/locate/cgvis).

Mitasova, H., W.M. Brown, J. Hofierka, 1994, Multidimensional dynamic cartography. Kartograficke listy, 2, p. 37-50.

Mitasova, H., and Mitas, L., 1993, Interpolation by Regularized Spline with Tension: I. Theory and implementation. Mathematical Geology, 25, 641-655.

Mitasova, H., and Hofierka, J., 1993, Interpolation by Regularized Spline with Tension: II. Application to terrain modeling and surface geometry analysis. Mathematical Geology , 25, 657-669.

Mitasova, H., Mitas, L., Brown, W.M., Gerdes, D.P., Kosinovsky, I., 1995, Modeling spatially and temporally distributed phenomena: New methods and tools for GRASS GIS. International Journal of GIS, 9 (4),special issue on integration of Environmental modeling and GIS, p. 443-446.

Mitasova, H., J. Hofierka, M. Zlocha, and R. L. Iverson, 1996a, 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)

Moore, I.D., and Burch, G.J., 1986a, Physical basis of the length-slope factor in the Universal Soil Loss Equation. Soil Sciences Society America Journal, 50, 1294-1298.

Moore I.D., Burch G.J., 1986b, Sediment transport capacity of sheet and rill flow: Application of unit stream power theory, Water Resources Research, v.22, p.1350-1360.

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

Moore, I.D., Grayson, R.B., and Ladson, A.R., 1991, Digital terrain modelling: a review of hydrological, geomorphological and biological applications. Hydrological Processes, 5, 3-30.

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

Moore, I.D., and Wilson, J.P., 1992, Length-slope factors for the Revised Universal Soil Loss Equation: Simplified method of estimation. Journal of Soil and Water Conservation, 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.

Renard G.K., Foster G.R., Weesies G.A., Porter J.P., 1991, RUSLE - Revised universal soil loss equation. Journal of Soil and Water Conservation, v. 46, p.30-33.

Rewerts C.C. and Engel B.A., 1991, ANSWERS on GRASS: Integrating a watershed simulation with a GIS. ASAE Paper No.91-2621. American Society of Agricultural Engineers, St.Joseph, Missouri, 1-8.

Rohling 1998

Rouhi A., and J. Wright, 1995, Spectral implementation of a new operator splitting method for solving partial differential equations, Computers in Physics, 9(5), 554-563

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.

Shapiro, M., and Westervelt, J., 1992 R.MAPCALC, an algebra for GIS and Image processing. US Army Corps of Engineers, Construction Engineering Research Laboratories, Champaign, Illinois, 422-425

Srinivasan, R. and Engel B.A., 1991, A knowledge based approach to extract input data from GIS, ASAE Paper No. 91-7045, American Society of Agricultural Engineers, St.Joseph, Missouri, 1-8.

Stephan, E.M., 1995, Exploratory data visualization of environmental data for reliable integration: Proceedings Auto-Carto 12, Charlotte, NC, p. 100-109.

Sutherland, R. A., 1991, Caesium-137 and sediment budgeting within a partially closed drainage basin, Zeitschrift Fur Geomorphologie N.F. Bd., 35, Heft 1, 47-63.

Vieux, B. E., N. S. Farajalla, and N. Gaur, Integrated GIS and distributed storm water runoff modeling, 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. 199-205, 1996.

Talmi, A., and Gilat, G., 1977, Method for Smooth Approximation of Data. Journal of Computational Physics, 23, 93-123.

Wahba, G., 1990, Spline models for observational data. CNMS-NSF Regional Conference series in applied mathematics, 59, SIAM, Philadelphia, Pennsylvania.

Warren, S.D., Diersing, V.E., Thompson, P.J., and Goran, W.D., 1989, An erosion-based land classification system for military installations. Environmental Management, 13, 251-257.

Wischmeier, W.H., and Smith, D.D., 1978, Predicting rainfall erosion losses, a guide to conservation planning. Agriculture Handbook No. 537, US Department of Agriculture, Washington D.C.

Zevenbergen, L.W. and Thorne, C.R., 1987, Quantitative analysis of surface topography. Earth Surface Processes and Landforms, 12, 47-56

Zhang, W. and Montgomery, D.R., 1994, Digital elevation model grid size, landscape representation and hydrologic simulations. Water Resources Research, 30, 1019-1028.

6. Appendix