NCSU GIS/MEA582:
Geospatial Modeling and Analysis

Erosion modeling

Resources: ESRI virtual campus Task: Compute soil detachment and net erosion/deposition using simple GIS-based models

Create the assignment working directory (folder) .\Erosion
Start ArcMap:

Start->Programs->ArcGIS->ArcMap

Check out a Spatial Analyst Extension license:

Under Tools->Extensions make sure there is a check next to 'Spatial Analyst'
Select View->Toolbars->Spatial Analyst to activate the extension

 
Geoprocessing>Environments>Workspace

Compute soil detachment rate using USLE3D

E = RKLSCP

where

• E is the annual soil detachment rate (soil loss)
• R is the rainfall
• K is the soil
• C is the land cover
• P is the preventive measures factor
• LS = 1.4 * (U/22.1)**0.4 * (sin (slope)/0.09)**1.2
  LS is the topography factor, U is contributing area per unit width (function of flowaccumulation)

For units and explanation of constants see Erosion modeling notes

First compute slope, and flowaccumulation for the elevation raster map e_lid792 (NOT elid_792):

In ArcMap
Add the raster feature layer: ncrast.mdb/e_lid792  (NOT elid_792)
Open the ArcToolbox

Compute slope_1m map using the Slope function, we will use it later:

Select 'Spatial Analyst Tools->Surface->Slope'
input: e_lid792 output: slope_1m
keep DEGREES for units

Compute the flow accumulation map Flowacc:

Compute Fill input=e_lid792 output=Fill_e_lid
Compute FlowDirection input=Fill_e_lid ouput=FlowDir_fil
Compute FlowAccumulation input=FlowDir_fil output=Flowacc

LS factor

Note that the lsfac map represents the topographic factor LS used in the USLE3D model.
Flowacc map multiplied by resolution is the contributing area per unit width U.

In the ArcToolbox, Select 'Spatial Analyst->Map Algebra->Raster Calculator

Type in the expression for lsfac layer (make sure you have spaces around operators).
Slope should be in radians, flowaccumulation is in number of cells so we multiply by resolution:

1.4 * (Power("Flowacc" * 1. / 22.1 , 0.4) * Power(Sin("slope_1m" * math.pi/180.0) / 0.09, 1.2))
Set the output raster to lsfac
Click 'OK'

Soil loss

E = RKLSCP

Multiplying all the factors in Raster Calculator we will compute the raster map of erosion rate

In the ArcToolbox, Select 'Spatial Analyst->Map Algebra->Raster Calculator
expression=
270.0 * "cfactb" * "lsfac" * "soils_kf"
output=s_loss

USPED model

We have already computed slope_1m and flowaccumulation, but we will also need direction of flow.
Compute aspect_1m map (direction of flow, direction of gradient vector) using the Aspect function:

'Spatial Analyst Tools->Surface->Aspect'
input=e_lid792, output=aspect_1m

Build an expression for sflowtopo (topographic factor LS for sediment flow T) using the Raster Calculator.
For the exponents use m=n=1, resolution is 1.

In the ArcToolbox, Select 'Spatial Analyst->Map Algebra->Raster Calculator
expression=
"Flowacc" * 1. * Sin("slope_1m" * math.pi/180.0)
output raster = sflowtopo

OR use m=1.3 and n=1.2 for study areas with extensive rills.
In this case, channels/streams will have large erosion rates due to high values of flowaccumulation.

In the ArcToolbox, Select 'Spatial Analyst->Map Algebra->Raster Calculator
expression=
(Power("Flowacc" * 1. , 1.3) * Power((Sin("slope_1m" * math.pi/180.0)), 1.2))
output raster = sflowtopo2

Compute divergence of sediment flow

In the ArcToolbox, Select 'Spatial Analyst->Map Algebra->Raster Calculator
expression=
"sflowtopo" * "soils_kf" * "cfactb" * 270. * Cos((- "aspect_1m" + 450.) * math.pi / 180.0) 
Output raster = qsx

In the ArcToolbox, Select 'Spatial Analyst->Map Algebra->Raster Calculator
expression=
"sflowtopo" * "soils_kf" * "cfactb" * 270. * Sin((- "aspect_1m" + 450.) * math.pi / 180.0) 
Output raster = qsy

Compute components of change in sediment flow in x and y direction as partial derivatives of sediment flow field, derived from slope and aspect - see eqs 1,2,3 from here.

In ArcToolbox
'Spatial Analyst Tools->Surface->Slope'
input=qsx, output=Slope_qsx

In the same way compute Aspect_qsx, Slope_qsy and Aspect_qsy

In ArcToolbox
'Spatial Analyst Tools->Surface->Aspect'
input=qsx, output=Aspect_qsx
'Spatial Analyst Tools->Surface->Slope'
input=qsx, output=Slope_qsy
In ArcToolbox
'Spatial Analyst Tools->Surface->Aspect'
input=qsx, output=Aspect_qsy

'Spatial Analyst->Map Algebra->Raster Calculator:
expression=
Cos((- "Aspect_qsx" + 450) * math.pi / 180.0) * Tan("Slope_qsx" * math.pi/180.0)
output=qsx_dx
expression=
Sin((- "Aspect_qsy" + 450) * math.pi / 180.0) * Tan("Slope_qsy" * math.pi/180.0)
output=qsy_dy

Finally. compute net erosion deposition:

'Spatial Analyst->Map Algebra->Raster Calculator:
expression=
"qsx_dx" + "qsy_dy"
output=erdep

Assign an appropriate color classification scheme to the erdep raster and generate a map for the report.
Note: the differences in the erdep are subtle and the default color table will almost assuredly render
the erosion deposition map in a way that little detail is revealed.
You are free to experiment with the the tools in Arc in order to define classification scheme which will reveal the erosion/deposition pattern.
Here's an option that should yield pretty good contrast between erosion and deposition:

Open Layer Properties for the erdep layer
Under the Symbology tab‚
Select Classified from options (Unique Values, Classified, Stretched, Discrete) in the left column
Set Classes to 11
Choose a broad, divergent color ramp from the Color Ramp drop down field
that ranges from brown grading through white to green
Click on Classify
Edit the eleven break values (located in column on far right in Classification dialog to:
-250000.00  (data minimum)
-50.00
-5.00
-1.00
-0.10
0.10
1.00
5.00
50.00
330000.00    (data maximum)
Click OK
Click Apply