USPED model for erosion/deposition

developed with GRASS5.0

Erosion/deposition modeling with USPED using GIS

Helena Mitasova, Lubos Mitas, University of Illinois at Urbana-Champaign

Copyright © 1999 Helena Mitasova and Lubos Mitas
(note: this document contains unpublished material, please cite this document and papers Mitasova et al 1996, Mitas and Mitasova 1998, when using this material in publications)

1. Method

USPED (Unit Stream Power - based Erosion Deposition) is a simple model which predicts the spatial distribution of erosion and deposition rates for a steady state overland flow with uniform rainfall excess conditions for transport capacity limited case of erosion process. The model is based on the theory originally outlined by Moore and Burch 1986 with numerous improvements. For the transport capacity limted case, we assume that the sediment flow rate qs(r) is at the sediment transport capacity T(r), r=(x,y) which is approximated by (Julien and Simons 1985)

|qs(r)| = T(r) = K_t(r) |q(r)|^m sin b(r)ⁿ

where b(r) [deg]is slope, q(r) is water flow rate, Kt(r) is transportability coefficient dependent on soil and cover, m, n are constants depending on the type of flow and soil properties. For overland flow the constants are usually set to m=1.6, n=1.3 (Foster 1993). Steadu state water flow can be expressed as a function of upslope contributing area per unit contour width A(r)[m] |q(r)| = A(r) i

where i[m] is uniform rainfall intensity (note: approximation by upslope area neglects the change in flow velocity due to cover). For the uniform soil and cover properties represented by Kt=const., the net erosion/deposition rate is estimated as a divergence of the sediment flow (see Appendix in Mitas and Mitasova 1998):

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

where s(r) is the unit vector in the steepest slope direction, h(r) [m] is the water depth estimated from the upslope area A(r), kp(r) is the profile curvature (terrain curvature in the direction of the steepest slope), kt(r) is the tangential curvature (curvature in 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 e.g., by the RST (Mitasova and Mitas, 1993; Mitasova and Hofierka, 1993; Krcho 1991). According to the 2D formulation, the spatial distribution of erosion and 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 interplay between the magnitude of water flow change and both terrain curvatures reflected in the bivariate formulation determines whether erosion or deposition will occur.

No experimental work was performed to develop parameters needed for USPED, therefore we
use the USLE or RUSLE parameters to incorporate the impact of soil and cover 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

T = R K C P A^m (sin b)ⁿ note:explain normalization/units here
where R~i^m, KCP~Kt and LS=A^m sin b ⁿ (add units here), and m=1.6, n=1.3 for prevailing rill erosion while m=n=1 for prevailing sheet erosion. Then the net erosion/deposition is estimated as ED = div (T . s) = d(T*cos a)/dx + d(T*sin a)/dy

where a [deg] is aspect of the terrain surface.This equation is equivalent to the relationship with curvatures presented above, however the computational procedure is simpler.
Caution should be used when interpreting the results because the USLE parameters were developed for simple plane fields and detachment limited erosion therefore to obtain accurate quantitative predictions for complex terrain conditions they need to be re-calibrated ( Foster 1990, Mitasova et al 1997 reply).

2. GIS implementation

2.1 GRASS GIS

Given data:

raster: elevation, K, C, (P)

constant: R=120,resolution=10

Computation

1. r.flow elevation dsout=flowacc

2. r.slope.aspect elevation slope=slope aspect=aspect

3. r.mapcalc

      explain normalization/units here

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

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

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

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

4. r.slope.aspect qsx dx = qsx.dx

5. r.slope.aspect qsy dy = qsy.dy

6. r.mapcalc

        erdep = qsx.dx + qsy.dy , for prevailing rill erosion

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

7. Optional: create new color table, reclass to erosion/deposition risk classes,

run statistics, ...

2.2 ArcView-Spatial Analyst

Given data

grid: elevation, K, C, (P), Important! K and C should be multiplied at least by 100

      to avoid problems when displaying the results

constants: R=120,resolution=10

Computation

1. select elevation

   DERIVE SLOPE

   name it slope

   DERIVE ASPECT

   name it aspect

2. MAP CALCULATOR

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

       name it flowacc

3. MAP CALCULATOR add normalization in this step

   build an expression (prevailing rill erosion):

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

    Evaluate

   name it sflowtopo

  or for prevailing sheet erosion

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

    Evaluate

    name it sflowtopo

    build an expression

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

    Evaluate

    name it qsx

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

    Evaluate

    name qsy

4. select qsx

   DERIVE SLOPE

   name it qsx.slope

   DERIVE ASPECT

   name it qsx.aspect

5. select qsy

   DERIVE SLOPE

   name it qsy.slope

   DERIVE ASPECT

   name it qsx.aspect

6. MAP CALCULATOR

       build an expression

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

       name it qsx.dx

       build an expression

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

       name it qsy.dy

       build an expression

       [qsx.dx] + [qsy.dy] for prevailing rill erosion

       ([qsx.dx] + [qsy.dy]) * 10. for prevailing sheet erosion

       name it erdep

(keep in mind that the values are multiplied by the values used to multiply the C and K factors. Do not divide to get the correct values, because the pattern with values in interval <-1,1>, typical for well preserved areas as well as upland, relatively flat areas, will be lost.)

7. Optional: create a new color table, reclass to erosion/deposition risk classes,

   run statistics, ...

2.3 ArcGIS 8.1, ArcMap


Given data 

grid: elevation, K, C, (P)

constant: R=120, resolution=10


Computation 

1. Enable Spatial Analyst
   
   under View... Toolbars  select Spatial Analyst

2. Calculate Slope and Aspect

   from the Spatial Analyst toolbar, select Surface Analysis... Slope

   give the new theme name slope

   from the Spatial Analyst toolbar, select Surface Analysis... Aspect

   give the new theme name aspect

3. Raster Calculator 

   from the Spatial Analyst toolbar, select Raster Calculator

   build an expression:
   
   FlowAccumulation(FlowDirection([elevation]))

   Evaluate
   
   make the new theme permanent and change the name to flowacc

   build an expression in the Raster Calculator:

   (for prevailing rill erosion):
   Pow([flowacc] * resolution , 0.6) * Pow(Sin([slope] * 0.01745) , 1.3))

   (for prevailing sheet erosion):
   [flowacc] * resolution * Sin([slope] * 0.01745)

   Evaluate

   make the new theme permanent and change the name to sflowtopo

   build an expression in the Raster Calculator:
   # note : the original 280 was changed to 120 to reflect R=120 constant

   [sflowtopo] * [kfac] * [cfac] * 120 * Cos((([aspect] *  (-1)) + 450) * .01745)

   Evaluate

   make the new theme permanent and change the name to qsx

   build an expression in the Raster Calculator:

   [sflowtopo] * [kfac] * [cfac] * 120 * Sin((([aspect] *  (-1)) + 450) * .01745)

   Evaluate

   make the new theme permanent and change the name to qsy

4. select qsx in the Spatial Analysis Toolbar

   select Surface Analysis... Slope

   give the new theme name qsx_slope

   select Surface Analysis... Aspect

   give the new theme name qsx_aspect

5. select qsy in the Spatial Analysis Toolbar

   select Surface Analysis... Slope

   give the new theme name qsy_slope

   select Surface Analysis... Aspect

   give the new theme name qsy_aspect

6. Raster Calculator

   build an expression

   Cos((([qsx_aspect] * (-1)) + 450) * .01745) * Tan([qsx_slope] * .01745)

   give the new theme name qsx_dx

   build an expression

   Sin((([qsy_aspect] * (-1)) + 450) * .01745) * Tan([qsy_slope] * .01745)

   give the new theme name qsy_dy

   build an expression

   [qsx_dx] + [qsy_dy] for prevailing rill erosion
   ([qsx_dx] + [qsy_dy]) * 10. for prevailing sheet erosion

   give the new theme name erdep

2.4 Notes

We were not able to visualize the erosion/deposition map when using the true values of K and C, because ArcView has troubles handling floating point numbers in the interval <-1,1> . According to the Spatial Analyst manual floating point grid themes can only use equal interval or standard deviation for classification which are unsuitable for erosion/deposition maps. Displaying the maps, some map calculator operations as well as reclassification operations with erosion/deposition maps then yield incorrect results. Here are some examples of what we tried to do to make erosion/deposition maps from the computation with true values of K and C:

erdep1: changing legend (values <-1,1> are completely ignored)
erdep2: translating the values to positive territory
erdep3: translating to positive territory and multiplying by 100 (the i-table shows the value at its upper left corner - the color is obviously incorrect)
erdep.reclas: extracting only the values between (-1,1> yields a map with 0 and 1 with information about the spatial distribution of erosion/deposition between <-1, 1> lost.
erdep100000: thanks to Tomas Cebecauer we found the way how to avoid loss of small erosion/deposition values by using K*100 and C*1000 instead of their original values. It is important to remember, when performing quantitative analysis that the values in the erdep map are multiplied.
erdep4a: output the resulting grid into GRASS5.0 and viewing it with SG3d (or NVIZ), it can be smoothed by r.neighbors: erdep4b

It is important to note that the algorithms available in ARC are different from those in GRASS, therefore to obtain acceptable result more attention should be paid to the proper selection of resolution and to the quality of DEM.

3. References and links to related sites

Mitas, L., Mitasova, H., 1998, Distributed erosion modeling for effective erosion prevention. Water Resources Research Vol. 34, No. 3, pp. 505-516.

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

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

Mitasova 1996 Mitasova, H., L. Mitas, B.M. Brown, D.P. Gerdes, I. Kosinovsky, 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 et al. 1997,98,99 Applications presented at this WWW site:

Zhang www
Geomodel www

Back to erosion
Spatial modeling
GMSlab