### 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 extensionSet your workspace and scratch workspace to C:\path\.Erosion

Geoprocessing>Environments>Workspace

### Compute soil detachment rate using USLE3D

Compute a map of soil detachment rate using USLE3D equation which assumes that there is no deposition.This represents detachment limited sediment transport by overland flow.

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):

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

Build an expression in the Spatial Analyst Raster Calculator for lsfac map using the Flowacc and slope_1m layers you just generated.
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

Now build the expression for soil detachment rate E (soil loss) s_loss. From the formulaE = RKLSCPWe computed LS factor R factor is considered to be constant and set to 270.0 K factor values can be added as a layer ncrast.mdb/soils_kf C factior values can be added as a layer ncrast.mdb/cfactb P factor (prevention measures) is here taken to = 1.

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_lossAdjust color table for the s_loss map - histogram equalized color ramp may be needed Find average soil loss rate in the area from the s_loss map (in Properties:Source).

### USPED model

Compute net erosion and deposition as a divergence of sediment flow vector field.Note that custom color tables that cover entire range of values including negative are important in this section.

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

First compute components of sediment flow in x and y direction: Calculate qsxIn 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 = qsxCalculate qsy

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=qsy, output=Slope_qsy In ArcToolbox 'Spatial Analyst Tools->Surface->Aspect' input=qsy, output=Aspect_qsyCompute qsx_dx and qsy_dy

'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 ApplyYou should now see the erosion and deposition areas rendered with pretty good contrast.