NCSU GIS/MEA582:
Geospatial Modeling and Analysis

Erosion modeling

Resources: ESRI Academy

Launch ArcGIS Pro and login with ArcGIS Online credentials if required. Select New> Map> Create New Project. Create a new project in your preferred workspace, in the instructions below we are using the default C:\Users\myname\Documents\ArcGIS\Projects\ folder.

Use your judgement to decide which maps to include in your report.

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 flow accumulation)

For units and explanation of constants see Erosion modeling notes

First compute slope, and flow accumulation for the elevation raster map e_lid792 (NOT elid_792): --Leave other values defaul if not specified--

In the Map tab, in the Layer group, click 'Add Data'
Browse to C:\Users\myname\Documents\ArcGIS\582data\ncrast.gdb
select 'e_lid792'  (NOT elid_792)
"Ok"

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

In the Geoprocessing pane, search for and select the 'Slope' tool
input'e_lid792
output=slope_1m
Keep 'DEGREES' for output measurement

Compute the flow accumulation map Flowacc:

In the Geoprocessing pane, search for and select the 'Fill' tool
input=e_lid792 
output=Fill_e_lid

In the Geoprocessing pane, search for and select the 'Flow Direction' tool
input=Fill_e_lid 
ouput=FlowDir_fil

In the Geoprocessing pane, search for and select the 'Flow Accumulation' tool
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.

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:

In the Geoprocessing pane, search for and select the 'Raster Calculator' tool
expression=
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 'Run'

Soil loss

E = RKLSCP

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

Add ncrast.gdb/soils_kf
Add ncrast.gdb/cfactb
In the Geoprocessing pane, search for and select the 'Raster Calculator' tool
expression=
270.0 * "cfactb" * "lsfac" * "soils_kf"
Set the 'output raster' to '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:

In the Geoprocessing pane, search for and select the '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 Geoprocessing pane, search for and select the 'Raster Calculator' tool
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 Geoprocessing pane, search for and select the 'Raster Calculator' tool
expression=
(Power("Flowacc" * 1. , 1.3) * Power((Sin("slope_1m" * math.pi/180.0)), 1.2))
output raster = sflowtopo2

Compute divergence of sediment flow

Be sure the 'cfactb' layer from earlier is in your map or you will recieve an error-
In the Geoprocessing pane, search for and select the 'Raster Calculator' tool
expression=
"sflowtopo" * "soils_kf" * "cfactb" * 270. * Cos((- "aspect_1m" + 450.) * math.pi / 180.0)
Output raster = qsx

In the Geoprocessing pane, search for and select the 'Raster Calculator' tool
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 the Geoprocessing pane, search for and select the 'Slope' tool
input=qsx
output=Slope_qsx
In the Geoprocessing pane, search for and select the 'Aspect' tool
input=qsx
output=Aspect_qsx
In the Geoprocessing pane, search for and select the 'Slope' tool
input=qsy
output=Slope_qsy
In the Geoprocessing pane, search for and select the 'Aspect' tool
input=qsy
output=Aspect_qsy

In the Geoprocessing pane, search for and select the 'Raster Calculator' tool
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:

In the Geoprocessing pane, search for and select the 'Raster Calculator' tool
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:

RMC the erdep layer->Symbology
Select 'Classify' from 'Primary symbology'
Set Classes to 11
Choose a broad, divergent color ramp from the Color Scheme drop down field
that ranges from brown grading through white to green (brown-green continuous)
Edit the eleven break values (located in the classes tab) 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)