Terrain Modeling
and Soil Erosion Simulation
Final Report
Prepared for:
U.S. Army Engineering Research and Development Center (ERDC),
under DACA8899D0002 Task Order No. 05,
Project Title: Terrain Modeling and Soil Erosion Simulation,
Matthew Hohmann, POC/COR ERDC/USACERL,
under contract to the University of Illinois Board of Trustees,
Douglas M. Johnston, Principal Investigator.
1.1 Simulation and modeling
1.2 Types of models
2.1. Exponents in RUSLE3d
2.2 Exponents in USPED
5.1 Comparison with NRCS erosion inventory
5.2 Planning prevention measures with GIS
One of the key challenges in environmental research is to model interacting physical processes with sufficient accuracy and efficiency. Rapid development of computer technology offers new opportunities to tackle extremely complex environmental problems. Computational approaches belong to "young" methodologies which were developed only over the past few decades and as such, they have their own rules, challenges, successes and limitations. The role of algorithms, data structures, computationally efficient methods, advanced visualization and exploration of parallelism are crucial for new advances in environmental research and require close collaboration between traditional research disciplines and computational science.
Originally, GIS applications were focused on static spatial data processing, analysis and computer cartography. However, development of new mapping technologies and computer capabilities together with acute environmental problems have pushed the GIS applications into more sophisticated levels. Advanced geoscientific applications involve modeling landscape processes at an unprecedented level of detail. Nevertheless, the processbased modeling of geospatial phenomena involves substantially more uncertainty than modeling in physics or chemistry because of the above mentioned complexity. Practical solutions have to rely, then, on the best possible combination of physical models, empirical evidence, intuition and available measured data. In physics, accuracy is usually understood in a much stricter sense, because many fundamental laws are known over a broad range of scales in energy, distance or time. For example, the Schrodinger equation describes matter at the electronic level virtually exactly, that means, within spectroscopic accuracy of 6 to 12 digits. This is seldom the case in complex geoscientific applications where 50% differences between measurements and model predictions can be in many instances considered satisfactory.
Computational approaches to investigations of physical systems are based on simulation and modeling. By simulation we mean a computer representation of reality in which the simulated system is governed by a set of known physical laws expressed in mathematical language. In this case the model is already in place and its range of validity and accuracy is supposed to be wellknown and verified. The task of simulation is to solve it for a particular realistic situation. The fact that the fundamental laws are known does not mean that simulations are straightforward or easy. The corresponding equations are often difficult to solve and barriers in computational feasibility and efficiency often limit accuracy, resolution or size of the modeled system.
On the other hand, by modeling we mean a process by which the scientist is trying to build a simplified version of reality for a phenomenon for which the fundamental laws are either unknown, impractical to use or simply do not exist. This typically involves systems which are very complex, with many constituents and a variety of interactions between them and with limited amount of available experimental information. Typical examples are many ecological models, and systems which involve anthropogenic activities. The modeling process often involves trial and error and in some cases its predictive power, accuracy and relation to reality may be a research problem on its own.
Physical systems are often described by a combination of deterministic (physically based) and empirical (observation based) models. Empirical models are based on observations and statistical analysis of observed data and their applicability is limited to the conditions for which they were developed. They can provide a rough picture of the phenomenon under study, but they cannot explain how the system works. Because of their simplicity, they are widely used for practical applications and as components of more complex models for the subprocesses for which the physical model is unknown or too complicated.
Current research and development in physical systems modeling is focused on distributed, processbased models, often dynamic in three dimensional (3D) space. This trend has been stimulated by the availability of geospatial data and supporting geographic information system (GIS) tools. GIS has greatly reduced the time for preparation of inputs, however, this task can still be rather tedious and time consuming. Artifacts from interpolation remain a perennial problem and the support for high precision floating point, temporal and multidimensional data is still inadequate. Therefore coupling of GIS and models is done at various levels and incorporation of GIS functionality within the modeling systems is now quite common. The models of physical processes are often core modules for integrated modeling systems, due to their impact on ecosystems and society. Models of physical systems are also important components of decision support systems.
A typical example of modeling of a complex landscape process which requires combination of empirical and physicsbased methods is soil erosion, and sediment transport and deposition modeling. This combination is necessary because there are still gaps in understanding of the relevant processes. Also, some of the processes are too complex or require difficult to measure parameters to be incorporated using the known physics principles.
A number of erosion models explicitly distinguish and simulate various types of erosion processes, such as sheet erosion, rills, gullies etc. However such an approach requires a large number of empirical input parameters which cannot be obtained at sufficient accuracy for landscape scale applications without substantial cost. Therefore in simpler models, the different impact of these various processes is averaged and incorporated through empirical parameters, which is the case for RUSLE3d and USPED models. These two models and their various applications have been covered in previous reports and publications (Mitasova et al. 19972000) and in the Landscape Erosion Modeling Tutorial. Here we provide some explanations regarding the types of erosion modeled and selection of parameters (Engel. 1999, Desmet and Govers.1996.).
Both models
replace slope length by upslope area as a measure of water flow. This
makes the models applicable to complex topography. It also means that
the model captures impact of a wider range of types of flow than the
original USLE. It includes the combined, averaged impact of sheet and
rill flow on hillslopes as well as concentrated flow erosion and
potential for gully formation that has not been covered by
traditional USLE. While the sheet and rill flow rates are comparable
to USLE for the range of slopes and hillslope length typical for USLE
applications erosion rate estimates for
concentrated flow erosion and very long hillslopes are not sufficiently verified because of lack of
empirical data. However,
preliminary comparisons show that the models correctly predict over one to two
magnitudes of increase in erosion rates due to concentrated flow and
these areas should be classified as areas of high erosion risk
(rather than excluded from the results). When making erosion risk assessment using RUSLE3d and USPED it is therefore not necessary to
add impact of gullies from field observations because they are
already incorporated.
The basic equation for RUSLE estimates erosion rate by multiplying the following empirical factors:
E = R K LS C P
where E is the average soil loss, R is the rainfall intensity, K is the soil, C is the cover, P is the prevention practices and LS is the topographic factor:
LS(r) = (m+1) [ A(r) / a_{0 }]^{m } [ sin b(r) / b_{0 }]^{n}
where A is upslope contributing area per unit contour width, b is the slope, a_{0 } = 22.1m, b_{0} = 0.09. Typical values for m=0.40.6 and n=11.3. Exponents for water and slope terms in the sediment transport and detachment equations reflect the interaction between different types of flow and soil detachment and transport. We explain their application and impact on the result in the following sections.
For sheet
flow, detachment and sediment transport increases relatively
slowly with the amount of water. Geometric properties of topography
(slope, curvatures) play a more important role in the evolution of the
pattern of soil detachment and net erosion/deposition than the
pattern of water flow. This type of flow is typical for areas with
good vegetation cover but also for a severely compacted, smooth soil
cover where compaction prevents soil detachment and formation of
rills. This type of flow is reflected by the lower value of exponent
m for the water term represented by the upslope area (Figure
1.).
Figure. 1 RUSLE3D LS factor for prevailing sheet flow reflected by the value of m=0.1. In this case the impact of steeper slope on the erosion pattern is stronger than the impact of concentrated flow.
SUM = 334314 Number of cells: 735318 Minimum: 0 Maximum: 3.1872103214 Range: 3.18721 Arithmetic Mean: 0.454652 Variance: 0.137468 Standard deviation: 0.370766  MAP: 1.1*exp(house3.dsd*3./22.13,0.1)*exp(sin(sl3.rst50s2)/0.0896,1.3)    Category Information  %    #description  cover acres  01stable . . . . . . . . . . . . . . . . . . . . . . . .  89.911476.09046 15low. . . . . . . . . . . . . . . . . . . . . . . . . .  9.62 157.96604  *no data. . . . . . . . . . . . . . . . . . . . . . . .  0.47 7.64897 
If both types
of flow are present in the given area, which is usually the case due
to spatial variability in land cover and soil properties, the value
m=0.4 provides reasonable, averaged results (Figure 2). It has
also a theoretical foundation (Moore and Burch 1986) and is being
used, for example, by Engel (1999) This exponent balances the impact
of turbulent and sheet overland flow.
Figure 2. RUSLE3D LS factor for averaged sheet, rill, and gully erosion including both impact of concentrated flow and slope pattern using m=0.4.
SUM = 666252 Number of cells: 735318 Minimum: 0 Maximum: 24.6981582642 Range: 24.6982 Arithmetic Mean: 0.906074 Variance: 0.779073 Standard deviation: 0.882651  MAP: 1.4*exp(house3.dsd*3./22.13,0.4)*exp(sin(sl3.rst50s2)/0.0896,1.3)    Category Information  %    #description  cover acres   01stable . . . . . . . . . . . . . . . . . . . . . . .  68.451123.70468  15low. . . . . . . . . . . . . . . . . . . . . . . . .  30.77 505.07843  510moderate . . . . . . . . . . . . . . . . . . . . . .  0.28 4.51560 1050high . . . . . . . . . . . . . . . . . . . . . . . .  0.05 0.75779  *no data. . . . . . . . . . . . . . . . . . . . . . .  0.47 7.64897 
For
prevailing rill and gully erosion that appears on disturbed
land, with soils vulnerable to rilling and gully formation, creating
conditions for highly turbulent water flow, the impact of the water term
is much higher (Figure 3), reflected by the higher exponent m.
Figure 3. RUSLE3D LS factor for prevailing rill and gully erosion reflected by the value of m=0.6. In this case the impact of flow on the erosion pattern is stronger than the impact of slope and overall erosion rates are much higher.
SUM = 1112308 Number of cells: 735318 Minimum: 0 Maximum: 120.2604675293 Range: 120.26 Arithmetic Mean: 1.51269 Variance: 5.92545 Standard deviation: 2.43422  MAP: 1.6*exp(house3.dsd*3./22.13,0.6)*exp(sin(sl3.rst50s2)/0.0896,1.3)    Category Information  %    #description  cover acres   01stable . . . . . . . . . . . . . . . . . . . . . .  51.96 853.10417  15low. . . . . . . . . . . . . . . . . . . . . . . .  43.84 719.65615  510moderate . . . . . . . . . . . . . . . . . . . . .  3.12 51.25830  1050high . . . . . . . . . . . . . . . . . . . . . . .  0.57 9.30009 505000severe . . . . . . . . . . . . . . . . . . . . . .  0.04 0.73779  *no data. . . . . . . . . . . . . . . . . . . . . .  0.47 7.64897 
If we assume
that a dense vegetation cover prevents creation of rills and keeps
water flow in dispersed, sheet flow while the loose, bare soil has an
opposite impact, increasing the flow turbulence and formation of
rills, we can use a spatially variable exponent based on the cover. In
the following example (Figure 4), we have used the exponents m=0.2 for forest,
m=0.4 for grass, m=0.5 for dormant sparse grass and m=0.6 for
disturbed areas. The result increases the negative impact of
disturbed areas and reduces the impact of vegetated areas.
Figure 4. RUSLE3D LS for variable m exponent
SUM = 792530.019225 Number of cells: 735318 Minimum: 0 Maximum: 86.2473754883 Range: 86.2474 Arithmetic Mean: 1.07781 Variance: 1.7033 Standard deviation: 1.30511  MAP: (1.0+mvar)*exp(house3.dsd*3./22.13,mvar)*exp(sin(sl3.rst50s2)/0.089... (   Category Information  %    #description  cover acres   01stable . . . . . . . . . . . . . . . . . . . . . .  62.891032.51487  15low. . . . . . . . . . . . . . . . . . . . . . . .  35.44 581.81698  510moderate . . . . . . . . . . . . . . . . . . . . .  1.06 17.44907  1050high . . . . . . . . . . . . . . . . . . . . . . .  0.14 2.22891 505000severe . . . . . . . . . . . . . . . . . . . . . .  0.00 0.04667  *no data. . . . . . . . . . . . . . . . . . . . . .  0.47 7.64897 
Note, that RUSLE provides a variable
m exponent expressed as a function of
slope angle, which reflects the fact that planar hillslopes with low
slope have prevailing dispersed flow, while flow on steeper slopes
can be more turbulent. However, this exponent gives reasonable
results only for shorter slopes in RUSLE3d. The RUSLE equation for variable m
was developed for slope length and traditional field applications, so it
is not recommended for use with upslope area without evaluating the
results using field measurements. The slopebased formula for m can
result in values of 0.8 or greater on steeper slopes. For slopes hundreds of meters long,
or for concentrated flow, this exponent predicts extremely high
erosion rates when using upslope area in RUSLE3d.
Figure 5. Full RUSLE3D estimate with a) constant m=0.4; b) variable m. The result with variable m predicts higher risk for gullies.
When different LS factors are used in the full RUSLE3D equation, the quantitative results and pattern (Figure 5) are as follows:
m=0.4 SUM = 6377646 Number of cells: 722067 Minimum: 734.1213378906 Maximum: 0 Range: 734.121 Arithmetic Mean: 8.83249 Variance: 328.88 Standard deviation: 18.135 m=variable SUM = 9074719 Number of cells: 722067 Minimum: 3863.8823242188 Maximum: 0 Range: 3863.88 Arithmetic Mean: 12.5677 Variance: 1269.59 Standard deviation: 35.6313
where b(r) [deg] is slope, q(r)is water flow
rate, Kt(r) is transportability coefficient, which is
dependent on soil and cover; m, n are constants that vary according to
type of flow and soil properties. For overland flow the constants are usually
set to m=1.6, n=1.3 (Foster 1990). Steady state water flow can be
expressed as a function of upslope contributing area per unit contour
width A(r)[m]
where i[m] is uniform rainfall intensity (note: approximation by upslope area neglects the change in flow velocity due to cover). No experimental work was performed to develop parameters needed for USPED, therefore we use the USLE or RUSLE parameters to incorporate the approximate impact of soil and cover and obtain at least a relative estimate of net erosion and deposition. We assume that we can estimate sediment flow at sediment transport capacity as:
ED = div (T . s) = d(T*cos a)/dx + d(T*sin a)/dy
where a [deg] is aspect of the elevation surface (or direction of flow, steepest slope direction, minus gradient direction).
Because USPED computes divergence of sediment flow, the impact of the exponents is more complex. The water flow term exponent here controls the ratio between the extent of erosion and deposition as it is illustrated by the following examples. It reflects the fact that turbulent flow can carry sediment farther and the impact of concentrated erosion will be wider than if flow is dispersed by vegetation. The following figures (Figures 6, 7, 8,9,10) demonstrate the influence of the exponent m on the erosion/deposition pattern and gully erosion risk estimates.
Figure 6. Spatial pattern
of topographic potential for erosion and deposition by USPED with
m=1, prevailing dispersed sheet flow, deposition extends high onto
hillslopes.
Figure 7. Spatial pattern
of topographic potential for erosion and deposition by USPED with
m=1.4, both sheet flow and rill flow influence erosion and deposition.
Deposition starts lower on the hillslope, gullies start farther towards the headwater
area and they are wider.
Figure 8. Spatial pattern of topographic potential for erosion and deposition by USPED with m=1.6, with prevailing rill and concentrated flow. Extent of deposition is further reduced and the potential for gullies starts very high in the headwater area, so gullies are even longer and wider.
Figure 9. Spatial
pattern of topographic potential for erosion and deposition by USPED
with variable m.
Figure 10. Full USPED with variable m
The RUSLE3d and USPED models as described above were implemented in ArcGIS 8.1 as well as in ArcView3.x and GRASS GIS. For ArcView, still the most widely used GIS on installations, we developed sample Avenue scripts to make it easier to run the models. We examined implementation in ArcGIS 8 using the ArcModel interface, but that interface was not released with ArcGIS 8.1 and is reportedly undergoing major revision with no release date set. Detailed procedures for running both models in ArcGIS are provided in the online tutorial developed as part of this work (see link below).
Using Soil Erosion Modeling for Improved Conservation Planning: A GISbased Tutorial (http://www.gis.uiuc.edu/erosion) is a richly illustrated hypertext tutorial, available online and on CDROM , produced as part of this contract. The tutorial covers erosion theory, data collection and evaluation, methodology, running the models using GIS, and interpreting results.
In the previous reports (Mitasova et al. 2000) we have applied the models to a predominately natural watershed at Fort Hood (Owl creek). Here we focus work on highly disturbed areas in the House Creek watershed. While the hydrology of this watershed has been investigated using various models, high resolution analysis of erosion and deposition by overland flow provides new insights about the distribution of sediment sources and their relation to land cover and topography. We obtained data from an NRCS erosion inventory for this area and entered them into our database. While the inventory data cannot be directly used for model calibration and validation because of the approach that was used, it provides some indication of the consistency of both the estimates of input parameters, and erosion rates estimated in the field at a few locations and from the GIS data.
SITES FILENAME: erosionf  Header Information:  name erosionhouse description erosion inventory data for House creek at Ft. Hood labels XY#ID%K%L%S%LS%C%P%ET%EA%EG%AG%ER Number of DIMENSIONS: 2    MIN     MAX   dim 1 604010.000000 613355.000000 Easting dim 2 3445785.000000 3452260.000000 Northing Type of CATEGORY information: CELL_TYPE    MIN     MAX   158 214 ID number Number of DOUBLE attributes: 11    MIN     MAX   dbl 1 0.2 0.39 Kfactor dbl 2 40 300 Slope Length [m] dbl 3 1 14 Slope angle dbl 4 0.12 1.45 LS dbl 5 0.003 0.45 Cfactor dbl 6 1 1 Pfactor dbl 7 0 28 dbl 8 0 13 Soil loss (USLE)(t/a/y) dbl 9 0 89 gully erosion (t/y) dbl 10 0 0.5 area gully affected (acre) dbl 11 0 152 road erosion (t/y) TOTAL SITES COUNTED: 50 
The inventory data is also available as a GRASS sites database file (http://www.gis.uiuc.edu/erosion/finalreport/erosionf.txt).
Notes about the soil erosion inventory data :
X  Y  ID  K  dsd/L  S  LS1.1  LS1.4  LS1.6  LSVAR  C  ET  EA  EG  AG  ER 
605462.

3451852



0.17

23

1.4dg

0.27

0.27

0.32

0.44  0.5  8/10    
605466  3451855  169  0.32  200  1.5  0.20        0.1  4  1.9  0  0  0 
605725

3452206

  0.17  14  1.9dg  0.30  0.39  0.45  0.40  0.1  2/2  +  
605721  3452202  160  0.32  200  3.5  0.46  0.1  9  4.3  24  0.2  32  
605797

3452086

0.17  37  2.0dg  0.38  0.85  1.42  1.1  0.1  4/5  +    
605799  3452086  162  0.32  150  4  0.47  0.2  17  8.2  39  0.5  3?  
606091  3452263  0.17  17  1.5  0.23  0.32  0.40  0.35  0.1  2      
606102  3452260  158  0.2  300  1.5  0.23  0.1  3  1.4  7  0.1  9  
606121 
3451954

0.17  20  2.3  0.43  0.67  0.88  0.7  0.5  14/17      
606132  3451953  164  0.2  200  4  0.53  0.4  28  13  86  0.4  0  
606220 
3451780

0.28  27  1.2  0.17  0.26  0.34  0.32  0.01  0.2  +  
606232  3451784  171  0.32  200  3  0.35  0.2  13  6.3  1  0.1  0 
From the point of view of conservation measures, contour filter strips, spreaders and, to some extent, grassed waterways, are effective in changing the turbulent flow typical for rills and gullies to dispersed sheet flow. Erosion prevention in areas with sheet flow should focus on upper convex parts of the slope (that are prone to increased net erosion by creating accelerated flow). Dense vegetation should be preserved in these areas. Areas with more turbulent flow require flow dispersal measures and prevention of concentrated flow development which can cause severe and irreparable soil loss. Once concentrated flow fully develops it is more difficult to control erosion and more expensive control measures, such as sedimentation ponds, need to be added.
The first step in identifying erosionprone areas is to select a numerical threshold from the results that identifies "hot spots", or areas to investigate for remediation (Figures 11, 12). Expert knowledge of the study area and field observations should be compared to values obtained by the model to help identify an appropriate threshold. Once a threshold value (or a series of priority values) is chosen, the hot spot areas are selected in the GIS. Appropriate vegetation may be added to the model in selected areas within the hotspots by changing the Cfactor for those areas. The model run is then repeated and the results compared to those from the initial condition. Repeating this process and visualizing the results for multiple scenarios, using the GIS to analyze area, types, and cost of remediation, should help in the development of a satisfactory remediation plan (Figures 13, 14).
Figure 11. Hot spots from RUSLE3d
Figure 12. Hot spots from USPED
Figure 13. Current (dark green) and needed new (light green) protective cover based on RUSLE3d
%area acres sparse cover 76.351253.40376 current forest 7.73 126.97239 proposed new forest/dense cover 13.66 224.23337
Figure 14. Current (dark green) and needed new (light green) protective cover based on USPED.
% area  acres sparse cover 79.581306.47097 current forest 6.94 113.91671 new dense vegetation 8.78 144.49688
The work summarized in this report provides a current application of stateoftheart modeling of erosion and sediment transport deposition processes. Modeling has been improved by incorporating more accurate representations of physical processes such as soil detachment, and by the availability of higher resolution (and higher accuracy) data sources. Because of the lack of appropriate field data, however, validation and calibration of model parameters and results remains a research need.
The analysis of the RUSLE3D and USPED exponents explains how the flow term exponent controls the magnitudes and spatial pattern of detachment and net erosion/deposition. Because there are no extensive experimental data available to assign the values of this exponent based on local conditions several approaches can be taken for its determination. The most accurate results can be obtained if the exponent is calibrated using spatially distributed measurements. If no field observations are available, it is possible to assign a value to m based on the land cover. Cfactor reflects the impact of cover on soil detachability, while m may incorporate the impact of land cover on the type of flow. The same cover may have a different impact on C and on m. Other models such as WEPP and SIMWE use the impact of land cover in several parameters as well. Other indicators of type of flow, such as surface roughness, propensity for seepage (which increases rilling and gullies, suggesting higher m), and soil properties can also be used. For many applications, especially where totals or averages are used, or the results are classified into a small number of classes, it is sufficient to use m=1.4. Because the range of values for slope is much smaller than the range of values for upslope area, the impact of its exponent on the results is smaller, but still significant, so if the experimental data are available it is useful to include a range of values into calibration.
The tutorial developed as part of this effort provides a more complete discussion of work accomplished to date, as well as examples and instructions for application.
Desmet, P. J. J., and G. Govers (1996) A GIS procedure for automatically calculating the USLE LS factor on topographically complex landscape units, J. Soil and Water Cons., 51(5), 427433.
Engel, B. (1999) Estimating Soil Erosion Using RUSLE (Revised Universal Soil Loss Equation) Using ArcView. http://danpatch.ecn.purdue.edu/~engelb/abe526/gisrusle/gisrusle.html.
Foster, G. R., 1990, Processbased 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. 429445.
Julien, P.Y. and Simons, D.B. (1985) Sediment transport capacity of overland flow. Transactions of the ASAE, 28, 755762.
Mitasova, H., Mitas, L., Brown, W. M., Johnston, D., 2000, Terrain modeling
and Soil Erosion Simulation: applications for evaluation and design of
conservation strategies Report for USA CERL. University of Illinois,
UrbanaChampaign, IL.
Mitasova, H., Mitas, L., Brown, W. M., Johnston, D., 1999, Terrain modeling and
Soil Erosion Simulations for Fort Hood and Fort Polk test areas. Report for USA
CERL. University of Illinois, UrbanaChampaign, IL.
Mitasova, H., Mitas, L., Brown, W. M., Johnston, D., 1998, Multidimensional Soil
Erosion/deposition Modeling and visualization using GIS. Final report for USA
CERL. University of Illinois, UrbanaChampaign, IL.
Mitasova, H., Mitas, L., Brown, W. M., Johnston, D., 1997, Multidimensional Soil
Erosion/deposition Modeling. PART IV and V: Process based erosion simulation for
spatially complex conditions and its applications. Report for USA CERL.
University of Illinois, UrbanaChampaign, IL,
Moore, I.D., and Burch, G.J. (1986) Physical basis of the lengthslope factor
in the Universal Soil Loss Equation. Soil Sciences Society America Journal, 50,
12941298.
Wilson, JP, and Lorang, MS, 1999, Spatial Models of Soil Erosion and GIS:
Spatial Models and GIS:New Potential and New Models (M Wegener, and AS
Fotheringham, eds.), Taylor and Francis, London: 83108.