Chapter 1

MULTISCALE SOIL EROSION SIMULATIONS FOR LAND USE MANAGEMENT

Helena Mitasova

Geographic Modeling Systems Laboratory University of Illinois at Urbana-Champaign Urbana, Illinois, USA

helena@gis.uiuc.edu

Lubos Mitas

Department of Physics North Carolina State University Raleigh, North Carolina, USA

lmitas@unity.ncsu.edu

1. INTRODUCTION

Increasing pressures on the land and an improved understanding of human impacts on the environment are leading to profound changes in land management, with emphasis on integration of local actions with watershed-scale approaches. This trend has a significant impact on the development of supporting GIS and modeling tools. Complex, distributed, physics-based models are needed to improve understanding and prediction of landscape processes at any point in space and time. At the same time, land owners and managers working in the watersheds and fields need fast and easy to use models for which the input data are readily available.

Recent advances in Geographic Information Systems (GIS) technology and linkage of numerous models with GIS (e.g., Wilson and Lorang, 1999; Moore et al., 1993; Doe et al., 1996; Saghafian, 1996; Srinivasan and Arnold, 1994; Vieux et al., 1996; Johnston and Srivastava, 1999) have created a potential to develop an environment for coordination of conservation efforts at different management levels. These advances facilitate evaluation of prevention practices based not only on the type of the prevention measure, but also on their location within the watershed. Spatial analysis and simulation can also provide supporting information for allocation of resources to those areas and those types of practices which will provide the most effective protection.

To reflect the need for modeling at different levels of complexity, from fast, approximate estimates for risk assessment to more detailed simulations for predictions and land use design, a set of models with increasing complexity was developed (Mitas and Mitasova, 1998; Mitasova et al., 1999). The simple models RUSLE3D (Revised Universal Soil Loss Equation for Complex Terrain) and USPED (Unit Stream Power-based Erosion Deposition) are based on modifications of well established equations representing special cases of erosion regimes. The basic empirical parameters for these models are available, however their applicability to a wide range of conditions is limited. The new distributed, process-based model SIMWE provides capabilities to simulate more complex effects, however both the experimental and theoretical research are still very active and underlying equations, as well as the input parameters, are under continuing development.

This chapter presents theoretical basis for the distributed, processbased SIMWE model, describes the relation between the processes modeled by SIMWE and the RUSLE3D and USPED models and illustrates the use of the presented models for evaluation and design of different conservation measures.

2. METHODS

To model spatio-temporal distribution of sediment transport and erosion/deposition at any point and time, a complex system of interacting processes has to be simulated, including rainfall events, vegetation growth, surface, subsurface and ground water flow, soil detachment, transport and deposition. Excellent examples of continuous time simulation systems, which integrate a wide range of interacting processes important for land use management are WEPP (Flanagan et al., chapter 7), LISEM (Jetten and de Roo, chapter 8) or SWAT (Srinivasan and Arnold, 1994). Spatial components of these systems are usually based on 1D routing of water and sediment through homogeneous hydrologic units (e.g., hillslope segments, subwatersheds), limiting the range of spatial effects that can be simulated. To more accuratelly capture the impacts of spatially variable conditions a new generation of hydrologic and sediment transport models introduced 2D flow routing capabilities (SIBERIA: Willgoose and Gyasi-Agyei, 1995; CASC2d: Julien et al., 1995; Ogden and Heilig, chapter 14; CHILD: Tucker et al., chapter 12; SIMWE: Mitas and Mitasova, 1998). SIMWE uses spatially continuous approach to modeling of erosion and deposition with the modeled phenomena represented by multivariate functions and the flow of water and sediment described as multivariate vector fields rather than systems of 1D flows. The model was developed for subsystems describing overland water and sediment flow with focus on simulations of phenomena important for land use management. The concept can be extended to other processes, including the 3D subsurface flows.

2.1 Process-based Overland Water and Sediment Flow Model

The spatially continuous approach uses a bivariate form of continuity equations to describe water flow and sediment transport over a complex terrain with spatially variable rainfall excess, land cover and soil properties during a rainfall event.

Shallow overland flow. Shallow overland water flow is described by the bivariate form of the St Venant equations (e.g., Julien et al., 1995). The continuity relation is given by:

\frac{∂h(r, t)}{∂t} = i_{e}(r, t) − ∇ · q(r, t)  (1)

while the momentum conservation in the diffusive wave approximation has the form:

s_{f} (r, t) = s(r) − ∇h(r, t)  (2)

where r = (x, y) [m] is the position, t [s] is the time, h(r, t) [m] is the depth of overland flow, i_{e}(r, t) [m/s] is the rainfall excess = rainfall − infiltration, q(r, t) [m^{2}/s] is the unit flow discharge (water flow per unit width), s(r) = −∇z(r) is the negative elevation gradient, z(r) [m] is the elevation, and s_{f} (r, t) is the negative gradient of overland flow surface (friction slope).

For a shallow water overland flow, with the hydraulic radius approximated by the normal flow depth h(r, t) (Moore and Foster, 1990), the unit discharge is given by:

q(r, t) = v(r, t)h(r, t)  (3)

where v(r, t) [m/s] is the flow velocity. The system of equations (1-3) is closed using the Manning’s relation between h(r, t) and v(r, t)

v(r, t) = \frac{C}{n(r)}h(r, t)^{2/3}|s_{f} (r)|^{1/2}s_{f0}(r)  (4)

where n(r) is the dimensionless Manning’s coefficient, C = 1 is the corresponding dimension constant [m^{1/3}/s] (Dingman, 1984), and s_{f0}(r) = s_{f} (r)/|s_{f} (r)| is the unit vector in the friction slope direction. To account for spatially variable cover necessary for land use management, we consider n(r) as explicitly location dependent.

In this chapter, the solution of continuity and momentum equations for a steady state ∂h(r, t)/∂t = 0, is considered to be an adequate estimate of overland flow for the land management applications (Flanagan and Nearing, 1995). In addition, the flow is considered to be close to the kinematic wave approximation for which s_{f} (r, t) ≈ s(r) and after using Eq. (3), the Eq. (1) is given by:

∇ · [h(r)v(r)] = i_{e}(r)  (5)

In order to incorporate the diffusive wave effects at least in an approximate way, a diffusion-like term ∝ ∇^{2}[h^{5/3}(r)] can be incorporated into (5) (Mitas and Mitasova, 1998):

−\frac{ε(r)}{2}∇^{2}[h^{5/3}(r)] + ∇ · [h(r)v(r)] = i_{e}(r)  (6)

where ε(r) is a spatially variable diffusion coefficient. 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 Dingman, (1984) and Lettenmeier and Wood, (1992). The diffusion term, which depends on h^{5/3}(r) instead of h(r), makes Eq. (6) linear in the function h^{5/3}(r) which enables it to be solved using the path sampling method.

Erosion and sediment transport by overland flow. The basic relationship describing sediment transport by overland flow is the sediment continuity equation, 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 bivariate form of the continuity of sediment mass equation is (e.g., Hong and Mostaghimi, 1995):

\frac{∂[ρ_{s}c(r, t)h(r, t)]}{∂t} + ∇ · q_{s}(r, t) = sources − sinks = D(r, t) (7)

where q_{s}(r, t) [kg/(ms)] is the sediment flow rate per unit width, c(r, t) [particle/m^{3}] is sediment concentration, ρ_{s} [kg/particle] is mass per sediment particle, ρ_{s}c(r, t) [kg/m^{3}] is sediment mass density, and D(r, t)

[kg/(m^{2}s)] is the net erosion or deposition rate. The sediment flow rate q_{s}(r, t) is a function of water flow and sediment concentration:

q_{s}(r, t) = ρ_{s}c(r, t)q(r, t)  (8)

Again, we assume a steady state form of the continuity equation:

\frac{∂[ρ_{s}c(r, t)h(r, t)]}{∂t} = 0 ⟶ ∇ · q_{s}(r) = D(r).  (9)

The sources and sinks term is derived from the assumption that the detachment and deposition rates are proportional to the difference between the sediment transport capacity and the actual sediment flow rate (Foster and Meyer, 1972):

D(r) = σ(r)[T (r) − |q_{s}(r)|]  (10)

where T (r) [kg/(ms)] is the sediment transport capacity, σ(r) [m^{−1}] is the first order reaction term dependent on soil and cover properties. The expression for σ(r) = D_{c}(r)/T (r) is obtained from the following relationship (Foster and Meyer, 1972):

D(r)/D_{c}(r) + |q_{s}(r)|/T (r) = 1  (11)

which states that the ratio of the erosion rate to the detachment capacity D_{c}(r) [kg/(m^{2}s)] plus the ratio of the sediment flow to the sediment transport capacity is a conserved quantity (unity). To keep the model simple, σ(r is applied for both net erosion and deposition. Equation (11), proposed by Foster and Meyer (1972), is based on the observed relationship between soil detachment and transport described e.g., by Meyer and Wischmeier (1969). This concept is used in several erosion models including WEPP (Haan et al., 1994; Flanagan and Nearing, 1995). The qualitative arguments, experimental observations and values for σ(r) are discussed by Foster and Meyer, (1972) and Foster (1982). While it is possible to use other frameworks for estimation of σ(r), we have chosen the Foster and Meyer concept because of its simplicity and widespread use.

The sediment transport capacity T (r) and detachment capacity D_{c}(r) represent the maximum potential sediment flow rate and the maximum potential detachment rate, respectively, and there are numerous simplified equations representing these rates under different conditions. In the WEPP model they are expressed as functions of a shear stress (Foster and Meyer, 1972):

T (r) = K_{t}(r)[τ(r)]^{p}  (12)

D_{c}(r) = K_{d}(r)[τ(r) − τ_{cr}(r)]^{q}  (13)

where τ(r) = ρ_{w} gh(r) sin β(r) [P a] is the shear stress, β [deg] is the slope angle, p and q are exponents, K_{t}(r) [s] is the effective transport capacity coefficient, K_{d}(r) [s/m] is the effective erodibility (detachment capacity coefficient), ρ_{w} g is the hydrostatic pressure of water with the unit height, g = 9.81 [m/s^{2}] is the gravitational acceleration, ρ_{w} = 10^{3} [kg/m^{3}] is the mass density of water, and τ_{cr}(r) [P a] is the critical shear stress. The parameters and adjustment factors for the estimation of D_{c}(r), T (r) are functions of soil and cover properties, and their values for a wide range of soils, cover, agricultural and erosion prevention practices were developed within the WEPP model (Flanagan and Nearing, 1995).

The exponents p, q substantially influence model behavior. The WEPP model uses q = 1 and p = 1.5 which means that with increasing water flow, transport capacity increases faster than detachment. The comparison with the spatial extent of colluvial deposits in our previous study (Mitas and Mitasova, 1998) indicates that over long term period (decades) the pattern of erosion/deposition is closer to the results obtained with a lower value of p. This is in agreement with several studies which indicate that for the landscape scale modeling the Eqs. (12), (13) are not general enough for different types of flow and sediment transport processes present in complex landscapes (Willgoose et al., 1989; Kirkby, 1987; Willgoose and Gyasi-Agyei, 1995; Mitas and Mitasova 1998).

Similar to the water flow equation, the steady state sediment flow equation (Eq. 9) can be rewritten to include a diffusion term. First, a function representing the mass of water-carried sediment per unit area ϱ(r) [kg/m^{2}] is defined as:

ϱ(r) = ρ_{s}c(r)h(r)  (14)

and after introducing a small diffusion term ∝ ∇^{2}ϱ(r), the continuity equation is rewritten as:

−\frac{ω}{2} ∇^{2}ϱ(r) + ∇ · [ϱ(r)v(r)] + ϱ(r)σ(r)|v(r)| = σ(r)T (r) (15)

where ω[m^{2}/s] is the diffusion constant. On the left hand side of Eq. (15) 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’ acting on ϱ(r). The size of the diffusion constant is about one order of magnitude smaller than the reciprocal Manning’s constant so that the impact of the diffusion term is relatively small. It represents local dispersion processes of the suspended flow (Bennet, 1974), caused by microtopography which is not captured by the DEM. The diffusion term can be modified to reflect impact of various processes.

The water and sediment flow described by Eqs. (6) and (15) can be solved by the path sampling stochastic approach described in the following section.

2.2 Path Sampling Solution Method

Equations (6) and (15) have a similar form in which a linear differential operator O acts on a nonnegative unknown function γ(r) (either h(r) or ϱ(r)), while on the right hand side, there is a given source term S(r) (either i_{e}(r) or σ(r)T (r)):

Oγ(r) = S(r)  (16)

Denoting O^{−1} the inverse operator to O, the solution can be symbolically written as:

γ(r) = O^{−1}S(r)  (17)

or explicitly, using the Green’s function:

γ(r) = ∫_{0}^{∞} ∫ G(r, r′, p)S(r′)dr′dp  (18)

The Green’s functions are often used for expressing the solutions of linear differential equations in physical or mathematical applications. For the theory see, e.g., Karlin and Taylor, (1981); Glimm and Jaffe, (1972); Stakgold, (1979); Carslaw and Jaeger, (1947). G(r, r′, p) is given by the following equation and an initial condition:

\frac{∂G(r, r′, p)}{∂p} = −OG(r, r′, p) ; G(r, r′, 0) = δ(r − r′)  (19)

where r, r′ are locations, p is time and δ is the Dirac function. In addition, we assume that the spatial region is a delineated watershed with zero boundary condition which is fulfilled by G(r, r′, p). The corresponding equations can be solved, e.g., by projection methods (Rouhi and Wright, 1995). Another equivalent alternative is to interpret Eqs. (6),(15) and (16) as describing stochastic processes with diffusion and drift components (Fokker-Planck equations) and carry out the actual simulation of the underlying process using a path sampling method (Gardiner, 1985).

The method is based on duality between the particle and field representation of spatially distributed phenomena.Within this concept, density of particles in space defines a field and vice versa, field is represented by particles with corresponding spatial distribution of their densities. Using this duality, processes can be modeled as evolution of fields or

Figure 1.1. Path sampling solution of the continuity equation for water depth h(r) using duality between particle and field representation: a) water depth at 1 minute,

b) water depth after 24 minutes. The grid is 416x430 cells at 10m resolution.See the CDROM for animation.

evolution of spatially distributed particles (Fig. 1), with the solution obtained as follows.

First a selected number of particles, also called walkers or sampling points, is distributed according to the source S(r′). These walkers are then propagated according to the function G(r, r′, p), generating a number of sampling paths. Averaging of these path samples provides an estimation of the actual solution γ(r) with statistical accuracy proportional to 1/\sqrt{M} where M is the number of walkers (Fig. 2). The solution is not restricted to the steady state and the state of the modeled quantity at any given time p can be obtained by averaging the path samples at a given time p.

The path sampling technique has several unique advantages which are becoming even more important due to new developments in computer

Figure 1.2. Path sampling solution of the continuity equation for sediment flow and net erosion/deposition: a) results for 7000 walkers, b) results for 50 million walkers. The grid is 280x250(70,000 cells) at 2m resolution. See the CDROM for animation.

technology. Perhaps one of its most significant properties is robustness which makes it possible 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). In addition, the independence of sampling points makes the stochastic methods perfectly suited to the new generation of computers as they provide scalability from a single workstation to large parallel machines and computers distributed over diferent types of networks.

3. SIMPLIFIED SPECIAL CASES AND MODEL EXTENSIONS

Land use management poses specific challenges for hydrologic and erosion modeling because of the necessity to capture the spatial aspects of the modeled processes. In the following sections some of the modeling issues important from the point of view of landuse design are addressed, in particular the simplified erosion and deposition models, simulation of water flow in flat areas and depressions, and modeling with spatially variable resolution.

3.1 Simple erosion and deposition models

To satisfy the need for models which are easy to implement, simple to compute, and for which the data a readily available, it is useful to derive the models for special cases of sediment transport regimes. There are two limiting cases of erosion and sediment transport (Foster and Meyer, 1972; Hairsine and Rose, 1992; Tucker et al., chapter 12): (i) detachment limited, and (ii) sediment transport capacity limited.

The detachment limited case. This case is represented by σ → 0, which, after substituting into Eqs. (9)-(10) results in the net erosion equal to the detachment capacity D(r) = D_{c}(r), which can be directly computed from Eq. (13) using an estimate of water depth h(r) computed e.g. by the SIMWE water flow module. For this case, transport capacity exceeds detachment capacity everywhere, erosion and sediment transport is detachment capacity limited and therefore no deposition occurs.

The most common erosion model which represents this case is USLE and its revised version RUSLE (Lane et al., chapter 10). If we assume that τ_{cr}(r) = 0 and that the spatial distribution of steady state water flow is adequately represented by a function of upslope contributing area per unit width (e.g., Moore et al., 1993) the detachment limited erosion can be approximated by RUSLE3D, which for a point on a hillslope has the form (Moore et al., 1993; Mitasova et al., 1999):

D(r) = R(r)K(r)C(r)P (r)(m + 1)[A(r)/22.13]^{m}[sinβ(r)/0.09]^{n} (20)

where D(r)[ton/(acre.year) = 0.2242 kg/(m^{2}.year)] is the average annual soil detachment (soil loss) rate, A(r)[m^{2}/m] is the upslope contributing area per unit width, 22.13m is the length and 0.09 is the slope of the standard USLE plot, R(r)[(hundredsofft.tonf.in)/(acre.hr.year) = = 17.02 MJ.mm/(ha.hr.year)] is the rainfall energy factor, K(r) [ton.acre.hr/(hundredsofacreft.tonf.in) = 0.1317 t.ha.hr/(ha.MJ.mm)] is soil erodibility, C(r)[dimensionless] is the cover factor and P (r) [dimensionless] is the prevention measures factor (Haan et al., 1994; Lane et al., chapter 10). Single storm and monthly R(r) is also available, making Eq. (20) suitable for estimation of D(r) for single storms and for modeling of monthly soil loss distribution over a year (Haan et al., 1994). Exponents m, n depend on the prevailing type of erosion (sheet, rill) and the typical values are m = 0.4 − 0.6 and n = 1 − 1.3. Replacement of slope length used in the original formulation of USLE/RUSLE by the upslope area provides a better spatial description of increased erosion due to the concentrated flow without the need to a priori define these locations as inputs for the model.

Transport capacity limited erosion. For this case σ → ∞, which, after substituting into Eq. (10), leads to |q_{s}(r)| ≈ T (r) and net erosion/deposition can be computed directly as a divergence of the sediment transport capacity:

D(r) = ∇ · q_{s}(r) = ∇ · [T (r)s_{0}(r)]  (21)

where s_{0}(r) = s(r)/|s(r)| is the unit vector in the steepest slope direction. Mitas and Mitasova (1998) have shown that the results obtained for this case were close to the observed distribution of colluvial deposits in their study area, suggesting the prevailing influence of the transport capacity limited case on a long term pattern of deposition. Equation (21) was also used to demonstrate the impact of both tangential and profile curvatures and the importance of 2D flow routing for predicting net erosion/deposition pattern (Mitas and Mitasova, 1998).

The Unit Stream Power Based Erosion/Deposition model (USPED) estimates the transport capacity limited case of erosion/deposition using the idea originally proposed by Moore and Burch (1986). It combines the USLE/RUSLE parameters and upslope contributing area per unit width A(r) to estimate the sediment flow at sediment transport capacity:

q_{s}(r) = T (r) ≈ R(r)K(r)C(r)P (r)[A(r)]^{m}[sinβ(r)]^{n}  (22)

The net erosion/deposition D(r) is then estimated using Eq. (21) as:

D(r) = ∇ · [T (r)s_{0}(r)] = \frac{d[T (r)cosα(r)]}{dx} + \frac{d[T (r)sinα(r))}{dy}  (23)

where α(r)[deg] is the aspect of the terrain surface (direction of flow). The exponents m, n control the relative influence of water and slope terms and reflect the impact of different types of flow. The exponents leading to values of erosion rates on hillslopes consistent with the values from RUSLE are m = 1.6, n = 1.3 and they seem to reflect the pattern for prevailing rill erosion with erosion sharply increasing with the amount of water. The observed extent of colluvial deposits in our previous study (Mitas and Mitasova, 1998) indicated that a lower exponent m = 1 better reflects the pattern of compounded, long term impact of both rill and sheet erosion and averaging over a long term sequence of large and small events.

Models representing limiting cases of erosion are simple to compute in GIS by combining the flow-tracing and topographic analysis functions with map algebra. They can be applied to a single storm, monthly and annual estimates of soil detachment and net erosion/deposition. Caution should be used when interpreting the results from both RUSLE3D and USPED because the USLE/RUSLE parameters were developed for simple plane fields and detachment limited erosion. Therefore, to obtain accurate quantitative predictions for complex terrain and land cover conditions they need to be re-calibrated, especially in areas of concentrated flow (Foster, 1990; Mitasova et al., 1997 reply). While the capabilities of both RUSLE3D and USPED to accurately predict the rates of erosion and deposition at any point in the complex landscape are limited (fact which is true about almost any erosion model), they are useful tools for land management. Both models use readily available parameters and can provide valuable spatial information about: (i) the location of areas with high erosion risk from both shallow overland and concentrated flow, (ii) location of areas with deposition, and (iii) relative estimates of erosion and deposition rates for different land use alternatives and conservation strategies. Locations identified as high risk from both RUSLE3D and USPED should be primary targets for field erosion inventory (to validate the risk) and implementation of prevention/mitigation measures (if the high risk is confirmed in the field). Computation of net erosion and deposition is also useful for evaluation of the landscape’s capacity to deposit the eroded material before it can reach the streams.

3.2 Water Depth in Flat Areas and Depressions

While flat areas and depressions are not high erosion risk locations, they play an important role within watersheds by holding water and reducing water flow to neighboring steeper slopes or into the streams. Modeling water flow over terrain surfaces with ∇z > 0 can be successfully performed using kinematic wave approximation. However, flat areas and bottoms of depressions where ∇z → 0 pose a problem because the water flow direction becomes undefined. Incorporation of the spatially variable diffusion term ε(r) in the SIMWE model provides capabilities to approximately simulate water depth in these locations. By defining the diffusion term ε(r) as a function of water depth and the velocity of flow as a function of an approximate water flow momentum, the water fills the depressions or spreads in the flat area and flows out in the prevailing flow direction (Fig. 3b).

For the situations when the flat areas or depressions are drained by natural or man-made swales and channels, the water flow can be simulated by combining the gradient field derived from the DEM with the

MULTISCALE SOIL EROSION SIMULATIONS 13

Figure 1.3. Simulation of water flow through existing shallow depressions: a) 2D kinematic wave flow, b) approximate 2D diffusive wave flow with variable ε(r) and flow momentum, c) approximate 2D diffusive wave flow with a shallow, 1 degree slope channel running through the depressions. Simulation was performed by the SIMWE model using a 10m resolution DEM, 100x170 cells. See the CDROM for animation.

gradient of the drainage (Fig. 3c). The input vector field s(r) representing the flow direction is then defined as:

s(r) = (1 − δ_{ij})s_{e}(r) + δ_{ij}s_{d}(r)  (24)

where s_{e}(r) is the vector in the steepest slope direction derived from a DEM and s_{d}(r) is the vector representing the flow through the surface drainage which can be estimated, for example, from drainage line data using GIS tools and/or from field measurements.

3.3 Multiscale Water and Sediment Flow Simulation

Because the spatial unit modeled at the local, implementation level is part of a larger watershed, the evaluation of the impact of numerous, locally implemented conservation practices across the entire watershed requires multiscale approach. This approach links the high resolution, local level simulation with low resolution/regional simulation and is being implemented for the SIMWE model. It supports simulations with spatially variable resolution which can include the following cases: (i) study area is represented by several data sets with different resolutions and levels of accuracy, and the best available data are used for each subregion; (ii) study area is large, with spatially variable complexity and it is sufficient to run the more homogeneous areas at lower resolution while running the more complex areas or areas experiencing land use change at high resolution. Both spatially variable accuracy and resolution can be implemented by reformulating the solution through the Green’s function given by Eq. (18). The integral Eq. (18) can be multiplied by a reweighting function W (r):

W (r)γ(r) = ∫_{0}^{∞} ∫ W (r)G(r, r′, p)S(r′)dr′dp =
= ∫_{0}^{∞} ∫ G^{∗}(r, r′, p)S(r′)dr′dp  (25)

which is equal to the appropriate increase in accuracy (W (r) > 1) in the regions of interest while it is unity elsewhere. The function W (r) can change (abruptly or smoothly) between regions with unequal resolutions and in fact, can be optimally adapted to the quality of input data (terrain, soils, etc) so that the more accurate solution is calculated only in the regions with correspondingly accurate inputs.

The reweighted Green’s function G^{∗}(r, r′, p), in effect, introduces a higher density of sampling points in the region with large W (r). The statistical noise will be spatially variable as ≈ 1/[W (r)\sqrt{M}], where M

Figure 1.4. Simulation of steady state water flow at spatially variable resolution 10m and 2m: a) particle representation, b) corresponding water depth field. The low resolution area is 416x430 and the high resolution area is 660x375 rows and columns. See the CDROM for animation.

is the average number of samples resulting in the accuracy increase for the areas with W (r) > 1.

This approach provides an alternative to the finite element methods (see e.g. Tucker et al., chapter 12) because it uses multiple standard grids with the given resolutions instead of finite element meshes, which often lack adequate GIS support. The implementation is based on the multipass simulation. First, the entire area is simulated at lower resolution, and the walkers entering the high resolution area(s) are saved. The saved walkers are resampled according to the Eq. (25) by spliting each walker into a number of ”smaller” walkers which are randomly distributed in the neighborhood of the original walker. The model is then run at high resolution only for the given subarea, with the resampled walkers used as inputs (Fig. 4). If several different land use alternatives are considered for the given subarea, this approach can be used to perform simulations for each alternative only within the high resolution subarea. The approach also provides useful spatial information about the locations where water flows into the given subarea and where it flows out (Fig. 4).

4. LANDSCAPE SCALE EROSION PREVENTION PLANNING AND DESIGN

Interactions between different land covers and topography significantly influence the spatial distribution of surface water depth, sediment flow and net erosion/deposition. The capabilities to simulate these interactions at both watershed and field scale can support the design of sustainable, cost effective conservation strategies and erosion prevention measures. Mitasova and Mitas (1998) have demonstrated the use of SIMWE and GIS for finding an effective spatial distribution of protective grass cover for a small agricultural watershed. The following sections provide examples of a wide range of applications where simulations of water flow, sediment transport and erosion/deposition are used to support land use management at different scales and levels of complexity.

4.1 Watershed scale erosion risk assesment and evaluation of conservation strategies with simple distributed models

A large number of watershed associations are being organized EPA, 2000 with the goal of improving the management of America’s watersheds. The basis for this work are watershed management plans which identify the problems and set priorities in funding and implementation of conservation measures. GIS combined with simple erosion models and free spatial data available through National Spatial Data Clearinghouse (NSDI, 2000) provide a cost effective way to assess the current state of watersheds, as well as evaluate the impact and prioritize various conservation strategies. The simplified models RUSLE3D and USPED were applied to the Court Creek watershed which serves as a pilot area for the Illinois Department of Natural Resources (IDNR) program aimed at demonstration of community based watershed management with strong scientific support. First, the current erosion risk areas were identified and their pattern was analyzed using RUSLE3D and GIS tools (Fig. 5a). The analysis has shown a favorable land use pattern with protective forested buffers along the bigger streams and on steep slopes. However, the analysis has also shown that headwater areas and areas with lower slopes but convergent water flow are not sufficiently protected. These sensitive areas are relatively small and scattered, and the results from RUSLE3D indicate that only 16% area (10,000 acres) produces 87% of total detached soil available for transport.

The impact of several conservation strategies was then evaluated, with the following two alternatives presented here: (i) buffers along the streams extended from 90m to 300m; (ii) conservation areas based on the erosion risk map.

The comparison of the strategies in terms of loss of agricultural land and reduction of erosion is shown in Figure 5a. Elimination of high erosion potential would require reduction of agricultural land by 16% and extension of the criteria currently used for the conservation program by including the headwater areas and areas with convergent water flow. The results obtained from the USPED model (Fig. 5b) indicate that a substantial portion of the eroded soil moves only for a short distance and there is enough concave areas to deposit the sediment before it can enter the streams. These results support some recent observations and hypotheses (Roseboom and Mollahan, 1999; Trimble, 1999) that in the Midwestern watersheds most of the sediment observed in the streams originates within the streams and from erosion by concentrated flow rather than from hillslope erosion by shallow overland flow.

This application demonstrates that the simple models used with widely available data can be useful for preliminary assessment of erosion and sedimentation risk, identification of ”hot spots” in the watersheds and approximate evaluation of different conservation strategies. Design and implementation of best management practices with reduced maintenance cost and increased effectiveness requires high resolution data and process based modeling. Application of SIMWE to different land use design

Figure 1.5. Court Creek watershed erosion potential for different land use alternatives estimated by RUSLE3D: a) current land use and erosion potential,

b) buffers along the streams c) grass in all areas where erosion potential exceeds 10 t/acre year. The maps are shown for a subarea 5x6km at 10m resolution, numbers are for the entire watershed (see maps for entire watershed at http://www2.gis.uiuc.edu:2280/modviz/courtcreek/cc.html)

tasks illustrates the type of problems that the simulations can help to understand and solve.

4.2 Wetlands and Drainage

Topographic potential for wetlands. Preservation and restoration of wetlands is among the most important and popular best management practices. Their success depends on many factors, including a sufficient supply of water. The SIMWE hydrologic submodel was used to identify the locations within the Court Creek Pilot Watershed which have topographic conditions favorable for wetlands. Several simulations were performed for various rainfall intensities and uniform land cover and soil conditions assuming that the flow velocity is controlled only by the terrain gradient - the existing drainage and channels were not considered. Comparison of the resulting simulated water depth with existing wetlands has shown that the areas with existing wetlands have accumulated depth from one event overland flow of at least 0.3m (except for wetlands in the upland areas which were not captured by the 30m DEM used in this study). Using this threshold, a map for topographic potential for wetlands was computed using map algebra (Fig. 6). While the simulation was very simplified, the map can serve as a useful starting point for identification of land owners with suitable land for wetlands and for evaluation of the proposals for wetland locations.

Drainage location design. Simulation of spatial distribution of water depth provides valuable information also for an ”opposite” task - identification of locations which require drainage to prevent negative impact of standing water on yields. Using a high accuracy DEM ( 6m resolution, 0.05m vertical accuracy) interpolated from rapid kinematic survey data by the RST method (Mitas and Mitasova, 1999), the water depth distribution was simulated for a typical rainfall for Midwestern agricultural fields (9mm/hr) under saturated conditions. The resulting water depth map was used to evaluate suitability of the locations of current drainage under saturated conditions. The resulting water depth map was used to evaluate suitability of the locations of current drainage and to plan the location of new drainage network in the negatively affected field (Fig. 7). While the model was very useful for evaluating and planning a suitable spatial pattern of the drainage network, detailed soil data and more complex dynamic simulations are needed to design the size, depth and other parameters of the drainage.

4.3 Concentrated flow erosion and grassways

The suitability of the SIMWE model for spatial design of vegetation based best management practices was evaluated by application to small experimental watersheds and fields with planned or installed erosion prevention measures.

Concentrated flow erosion. Development of high erosion in areas of concentrated flow was studied by performing simulations of water flow

Figure 1.6. Map of topographic potential for wetlands representing locations which can hold water depth 0.3m and more without being drained by a channel. Existing wetlands are displayed as polygons. The area is 4x6km, simulated at 10m resolution (resampled from 30m resolution DEM

.

and net erosion deposition for an experimental field with uniform land cover (350x270m, modeled at 2m resolution; Zhang, 1999). For a short rainfall event ending before the flow has reached steady stte, the maximum erosion rate was on the upper convex part of the hillslope and there was only deposition in the center of the valley (Fig. 8a). As the duration of the rainfall increased, water depth in the center of the valley has grown rapidly until it reached a threshold when linear features with very high erosion rates developed within the depositional area, indicating potential for gully formation (Fig. 8b). This effect is modeled by both USPED (Mitasova et al., 1996, 1999) and SIMWE (Mitas and Mitasova, 1998),

Figure 1.7. Simulated spatial distribution of water depth for agricultural fields (2.5x4.5km) draped as color over the DEM (6m resolution, 30-times vertical exaggeration), with the existing drainage drawn as lines.

however, a smooth, high resolution DEM without artifacts is needed to realistically capture this commonly observed phenomenon (see Fig. 2c in Mitas and Mitasova, 1999). This example also demonstrates that for a dynamic event modeling incorporation of re-entrainment process is important (Hairsine and Rose, 1992).

Grassways. The common practice for prevention of erosion by concentrated flow are grassways. Their design is guided by the topographic conditions and roughness within the grassway, represented by Mannings coefficient (SCS, 1988). To investigate the actual impact of a

Figure 1.8. Water depth and net erosion/deposition pattern for 18mm/hr rainfall excess for a) short event, with only deposition in the valley center, b) long event leading to steady state flow, with both high erosion and deposition in the valley center, indicating a potential for gully formation. The 350x270m field is modeled at 2m resolution.

grassway, the water and sediment flow as well as net erosion/deposition pattern was simulated for a field within the Scheyern experimental farm (Auerswald et al., 1996; Mitas and Mitasova, 1998) for the bare field and after the installation of grassway with different values of roughness in the field. For the bare field, there is a potential for gully formation (Fig. 9a). After the installation of grassway the center of the valley becomes a depositional area. However, if the roughness in the field is several times smaller than in the grassway, high erosion develops around the grassway, potentially replacing one big gully with two smaller ones.

This is a known problem called double channeling which can substantially increase the cost of grassway maintenance (Fig. 9b). Increasing the roughness in the field reduces the risk of double channeling and the transition from erosion in the field to deposition in the grassway is relatively smooth (Fig. 9c).

Figure 1.9. Impact of grassway and differences in roughness on sediment flow and net erosion/deposition pattern: a) bare field without grassway with gully potential in the center, b) grassway (n=0.1) and the bare field (n=0.01) with sediment flow and erosion around the grassway, c) grassway (n=0.1) and the field with increased roughness (n=0.05) with smooth transition from erosion to deposition.

5. CONCLUSIONS

This chapter is focused on methodology and applications of simulation methods for prediction and solution of land management problems related to overland flow erosion. The presented approach aims at keeping the models, and particularly the number of required input parameters, as simple as possible while capturing the effects important for sustainable land use design. The applications demonstrate the need for a set of modeling tools with different levels of complexity to support land use management from strategic planning to design and implementation. To satisfy this need three interrelated models were prepsented.

The first model SIMWE is based on generalization of hillslope erosion model used in WEPP (Flanagan and Nearing, 1995). It models erosion regimes from detachment to transport capacity limiting cases, includes approximate diffusive wave effect and prefferential flow. It supports multiscale modeling which which can be further extended to incorporate multiscale-multiprocess simulations. It is possible to derive simple, explicit equations for the computation of limiting erosion regimes - for these cases SIMWE differs from RUSLE3D/USPED mainly by the use the WEPP parameters.

RUSLE3D/USPED differ from SIMWE in that they model only the limited cases of erosion and sediment transport, however they use readily available parameters and therefore are easy to implement and use. All of the presented models can be used for single storms as well as for long term averages.

The difference between the models can be illustrated by the prediction of erosion/deposition pattern on the upper border of a grass strip. While RUSLE3D will predict only reduction of erosion on the upper border of a field strip (no deposition), USPED will predict all sediment deposited on the edge, SIMWE allows to predict some deposition above the strip (due to approx. backwater effect) and deposition in grass strip spread over a longer distance depending on the size of the soil particles.

The applications of spatially continuous simulations revealed gaps in the theory of erosion processes in complex landscapes, especially in relation to the mathematical description of transport capacity suitable for practical application while at the same time general enough to capture the variability of flow at the landscape scale. Spatially distributed field experiments based on new technologies for field data collection, monitoring and remote sensing closely coupled with modeling are needed to improve our understanding of complex interactions involved in erosion processes and bring the quantitative accuracy of predictions (which is currently at about 50-150%) to acceptable and useful levels.

The report ”New Strategies for America’s Watersheds” (Committee on Watershed Management, National Research Council 1999) identifies simulation modeling as one area of special promise for watershed management. At the same time, this report analyzes the current status in watershed modeling for decision making and concludes that the available models and methods are outdated and ”a major modeling effort is needed to develop and implement state-of-the-art models for watershed evaluations”(pp.160-161). The presented approach along with other models presented in several chapters in this book, are a step towards the development and implementation of such tools.

Acknowledgments

We would like to acknowledge the long term support for this research from Geographic Modeling Systems Laboratory director Douglas M. Johnston as well as GIS assistance by William M. Brown. The funding was provided by the USArmy CERL, Strategic Environmental Research and Development Program - SERDP, Illinois Council on Food and Agricultural Research - CFAR and Illinois Department of Natural Resources. We greatly appreciate the sharing of data by K. Auerswald, S. Warren, K. Drackett, Zhang Yusheng, and D. Timlin. Our special thanks goes to two reviewers for their thorough review of the paper, as well as stimulating questions and comments which had a a very positive impact on the revised version of this chapter.

References

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 Universitaet Muenchen, Freising-Weihenstephan, Germany.

Bennet, J. P. (1974). Concepts of Mathematical Modeling of Sediment

Yield, Water Resources Research, 10, 485-496. Carslaw, H. S., and J.C. Jaeger. (1947). Conduction of Heat in Solids,

Oxford University, London. Dingman, S. L. (1984). Fluvial hydrology, Freeman, New York. Doe, W.W., B. Saghafian, and P.Y. Julien. (1996). Land Use Impact

on Watershed Response: The Integration of Two-dimensional Hydrological Modeling and Geographical Information Systems. Hydrological Processes, 10, 1503-1511.

EPA. (2000). Surf your watershed; http://www.epa.org.surf2/ Flanagan, D. C., and M. A. Nearing (eds.). (1995). USDA-Water Erosion

Prediction Project, NSERL, report no. 10, pp. 1.1- A.1, National Soil Erosion Lab., USDA ARS, Laffayette, IN.

Foster, G. R. (1982). Modeling the erosion processes, in: Hydrologic modeling of small watersheds, edited by C. T. Haan, H. D. Johnson, and D. L. Brakensiek, ASAE Monograph No. 5, ASAE, St. Joseph, MI, pp. 197-380.

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

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

Gardiner, C. W. (1985). Handbook of Stochastic Methods for Physics,

Chemistry, and the Natural Sciences, Springer, Berlin. Glimm J., and A. Jaffe. (1972). Quantum Physics. A Functional Integral

Point of View, Springer, Berlin. Govindaraju, R. S., and M. L. Kavvas. (1991). Modeling the erosion

process over steep slopes: approximate analytical solutions, Journal of Hydrology, 127, 279-305.

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.

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

Johnston, D.M. and A. Srivastava. (1999). Decision Support Systems for Design and Planning: The Development of HydroPEDDS (Hydrologic Performance Evaluation and Design Decision Support) System for Urban Watershed Planning, in: 6th International Conference on Computers in Urban Planning and Urban Management (CUPUMS’99), Venice, Italy (CDROM).

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.

Karlin, S., and H. M. Taylor. (1981). A Second Course in Stochastic

Processes, Academic Press, New York. Kirkby, M. J. (1987). Modelling some influences of soil erosion, land-

slides and valley gradient on drainage density and hollow development. Catena Supplement, 10, 1-14.

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.

Meyer, L.D., and W. H. Wischmeier, W.H. (1969) Mathematical simulation of the process of soil erosion by water. Transactions of the ASAE, 12, 754-758.

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

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

Mitasova, H., and L. Mitas. (1999). Modeling soil detachment by RUSLE

3d using GIS. http://www2.gis.uiuc.edu:2280/modviz/erosion/usle.html Mitasova, H., and L. Mitas. (1999). Erosion/deposition modeling with

USPED using GIS. http://www2.gis.uiuc.edu:2280/modviz/erosion/usped.html Mitasova, H., Mitas, L., Brown, W. M., and Johnston, D. (1999). Ter-

rain modeling and Soil Erosion Simulations for Fort Hood and Fort Polk test areas. Report for USA CERL. University of Illinois, Urbana-Champaign, IL; http://www2.gis.uiuc.edu:2280/modviz/reports/cerl99/rep99.html

Mitasova, H., J. Hofierka, M. Zlocha, and L.R. Iverson. (1997). Modeling topographic potential for erosion and deposition using GIS. Reply to a comment. Int. Journal of Geographical Information Science, 11(6), 611-618.

Mitasova, H., J. Hofierka, M. Zlocha, and L.R. Iverson. (1996). Modeling topographic potential for erosion and deposition using GIS. Int. Journal of Geographical Information Science, 10(5), 629-641.

Moore I.D., and Burch G.J. (1986). Modeling erosion and deposition:

Topographic effects. Transactions ASAE, 29, 1624-1640. Moore, I. D., and G. R. Foster. (1990). Hydraulics and overland flow,

in Process Studies in Hillslope Hydrology, edited by M. G. Anderson and T. P. Burt, John Wiley, 215-54.

Moore, I. D., A. K. Turner, J. P. Wilson, S. K. Jensen, and L. E. Band. (1993). GIS and land surface-subsurface process modeling, in: Geographic Information Systems and Environmental Modeling, edited by M. F Goodchild, L. T. Steyaert, and B. O. Parks, Oxford University Press, New York, 196-230.

National Research Council. (1999). New Strategies for America’s Water-

sheds (Washington DC: National Academy Press) NSDI. (2000). National spatial data infrastructure. http://www.nsdi.org/ Roseboom, D. and Mollahan, R. (1999). Lake Pittsfield National Mon-

itoring Project. Report for Illinosi State Water Survey and Illinois Environmental Protection Agency, Peoria, Il.

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. (1996). 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.

SCS. 1988. Manual. SCS DOPLN Srinivasan, R., and J. G. Arnold. (1994). Integration of a basin scale

water quality model with GIS, Water Resources Bulletin, 30(3), 453-462.

Stakgold, I. (1979), Green’s Functions and Boundary Value Problems,

John Wiley, New York. Trimble, S.W. (1999), Decreased rates of alluvial sediment storage in the

Coon Creek basin, Wisconsin. Science 285(8), 1244-1246. Vieux, B. E., N. S. Farajalla, and N. Gaur. (1996). 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.

Willgoose, G. R., and Y. Gyasi-Agyei. (1995). New technology in hydrology and erosion assessment for mine rehabilitations, Proceedings of the APCOM XXV Conference, Brisbane, pp. 555-562.

Willgoose, G. R., R. L. Bras, and I. Rodriguez-Iturbe. (1989). A physically based channel network and catchment evolution model, technical report no. 322, Ralph Parsons Lab., MIT, Cambridge, Mass., USA.

Wilson, J.P. and M.S. Lorang. (1999). Spatial Models of Soil Erosion and GIS. In Wegener M. and A.S. Fotheringham (Eds.), Spatial Models and GIS: New Potential and New Models (London: Taylor and Francis), 83-108.

Zhang Yusheng. (1999). GIS, Erosion and Deposition Modelling, and Caesium Technique. http://www.ex.ac.uk/ yszhang/welcome.htm