Surface Erosion AML

Hoodsport Watershed Analysis

/* Surface Erosion AML
/* Created by Susan Seaberg 3/1/99

*******AML for Road Erosion*******

/* Isolate roads within project area
rd = linegrid(trans,road.class.cd,#,#,32.823)
/*gridclip rd rds box 1351927.517 811058.007 1367387.150 824055.915
roads = con(watershed_dem > 0,rds)
/* Identify bedrock for soil erosion rates
soil_par = polygrid(soil,prnt.mat.spf.cd,#,#,32.823)
prnt_mat = con(watershed_dem > 0,soil_par)
/*gridclip prnt_mat prtmat box 1351927.517 811058.007 1367387.150 824055.915
/* parent material identified by the above equations are:  0 = no data; 11 = basic igneous/basalt
/* 50 = glacial drift; 51 = glacial till; 52 = glacial outwash; 55 = volcanic ash over glacial 
/* drift
/* Assume that all roads are older than 2 years
rd_prtmat_cd = con(prtmat == 0,1,prtmat == 11,10,prtmat == 50 or prtmat == 51 or~
 prtmat == 52 or prtmat == 55,30)
/*gridclip rd_prtmat_cd rd_pmat box 1351927.517 811058.007 1367387.150 824055.915
/* Basic Erosion rates
/* Tread and cutslope/ditch = 40% * rd_pmat_cd, fillslope = 20% * rd_pmat_cd
tread = con(roads >= 0,rd_pmat * 0.4)
/* Add correction factor by assume some vegetation has grown on cut and fill slopes since 
/* these are older roads,vegetation > 80% = .18(factor)
cut_ditch = con(roads >= 0,rd_pmat * 0.4 * 0.18)
fill = con(roads >= 0,rd_pmat * 0.2 * 0.18)
/* Assume that roads = 4, it is >6" gravel with a correction factor of .2;
/* roads  = 0(trail) it is native soil/rock and factor = 1
tread1 = con(roads == 0,tread,roads == 4,tread * 0.2)
/* Traffic and precip factor;  Given the 160" per year average in the area, 
/* The #4 road and trail will have 4 months of no traffic and the other 8 months
/* with light traffic giving them a factor of .25*.1 + .75*1 ; (.025+.75 = .775)
tread2 = con(roads == 0 or roads == 4,tread1 * .775) 
/* Determine areas of cutslope and fillslope
/* Trails, rw0 is 5' with no cut or fill slope; #4 roads rw4 = 12 then
/* calculate the cut and fill slopes with last years aml 
&sv rw0 = 5
&sv rw4 = 12
slope = slope(watershed_dem,percentrise)
/* Cutslope = cslope ; Fillslope = fslope
cslope = con(slope < 50,100,133)
fslope = con(slope < 45,67,80)
cw = con(roads == 4 and slope > 80,0,roads == 4 and slope < 30,%rw4% / 2 * slope~
 / (fslope - slope) / 2,roads == 4 and slope < 50,%rw4% * slope / (cslope - slope) * 2 / 3,~
 roads == 4 and slope >= 50,%rw4% * slope / (cslope - slope),roads == 0,0)
fw = con(roads == 4 and slope < 30,%rw4% * slope / (fslope - slope) / 2,~
 roads == 4 and slope < 50,%rw4% * slope / (fslope - slope) / 3,roads == 0,0)
/* Estimate sediment production from each road cell in Tons per year
road1 = con(roads >= 0,1)
tsed_tread = con(roads == 0,road1 * 32.283 * 5 * tread2 / 43560,roads == 4,road1 *~
 32.823 * 12 * tread2 / 43560)
cw_null = con(isnull(cw),0,cw)
tsed_cut = cw_null * cut_ditch * 32.823 / 43560
fw_null = con(isnull(fw),0,fw)
tsed_fill = fw_null * fill * 32.283 / 43560

*******AML for Road Delivery*******

/* Identify road cells that will deliver to streams by ditch
road_stream = con(road1 == 1 and streams_shed == 1,1)
edistroad = eucdistance(road_stream, #, #, 200)
roaddist = con(edistroad == 0,1,edistroad)
roadcell_del = roaddist * road1
/* Find Roads within 200' around stream
edist200 = eucdistance(streams_shed, #, #, 200)
null_edist200 = con(isnull(edist200),0,1)
road_win_200 = null_edist200 * (con(isnull(road1),0,1))
/* Separate out the roads crossing streams and the roads near streams
/* roads crossing is roadcell_del, but need to change so that it is 1,0
road_cross = con(isnull(roadcell_del),0,1)
/* Turn roads into streams
rd_and_str = con(road_cross == 1 or streams_shed == 1,1,0)
null_rdstr = con(isnull(rd_and_str),1,setnull(rd_and_str == 1,1))
flowdir = flowdirection(watershed_dem,channel_grade)
r_s_flowdir = flowdir * null_rdstr
rands_length = flowlength(r_s_flowdir,#,downstream)
/* Estimate the fraction of sediment delivered to the stream (not filtered out on its way)
/* as an exponential function of the distance from the stream.  Assume 50% is delivered in per cell.
delivery = pow(2,- rands_length)
null_del = con(isnull(delivery),0,delivery)
/* Roads crossing for cut and tread will dump 100% into streams
c_and_t200 = con(road_cross == 1, tsed_tread + tsed_cut,0)
f_200 = con(road_win_200 == 1,tsed_fill * .1,0)
c_t_f_gt200 = con(road_cross == 0 or road_win_200 == 0,null_del * (tsed_tread + tsed_cut +~
 tsed_fill))
/*

*******AML for USLE (Hillslope Erosion)*******

/* USLE equation in tons/acre, A = S*L*R*CP*K, 
/* Create length grids
up = flowlength(flowdir,#,upstream)
down = flowlength(flowdir,#,downstream)
/* Rainfall erosivity was given in 160 in/year (I think) needs to be in in/hr 
/* From Susan Bolton's Hydro class, USLE handout and NOAA atlases, R can be found by
/* R = 27.38 * (P (2yr 6 hr)^2.17
/* We have 2yr 24 hr , changed to 6 hr and put into the above equation yields
/* R = 27.38 * (.475 / 24 * 6)^2.17 = 0.268775
&sv r = 0.268775
/* cover is 95 - 100% and canopy cover is 80% (c = .003). Assume that there is no erosion control 
/* being practiced (p = 1)
&sv cp = 0.003
/* conversion to cells from acres
&sv a = 0.024732537
usle = (((slope / 100) * ((up + down) * 32.823)) / 10) * %r% * %cp% * (float(soils.soil_k)) * %a%
/*
/* USLE equation in tons/acre, A = S*L*R*CP*K, Assuming that where we aren't
/* harvesting anything else but the 55 or older trees, then CP = 0
/* Create age grid from riu
age = polygrid(riu,fiu.pri.age,#,#,32.823)
/* Find all trees older than 55
tre50 = con(age >= 50,1,0)
gridclip tre50 trees_50 box 1351927.517 811058.007 1367387.150 824055.915
cut_area = con(isnull(trees_50),(rand() * 1),trees_50)
usle2 = con(cut_area == 1 and slope < 30,(((slope / 100) * ((up + down) * 32.823) / 10) *~
 %r% * .33 * (float(soils.soil_k)) * %a%), con(cut_area == 1 and slope >= 30,(((slope / 100) * ((up + down) *~
 32.823) / 10) * %r% * 1 * float(soils.soil_k)) * %a%),usle))
/*
/*

*******AML for Hillslope and Road Delivery*******

usle_sed = usle * null_del
hill_road_sed = c_and_t200 + f_200 + con(isnull(usle_sed),0,usle_sed)~
 + con(isnull(c_t_f_gt200),0,c_t_f_gt200)
usle_sed2 = usle2 * null_del
hr_sed2 = c_and_t200 + f_200 + con(isnull(usle_sed2),0,usle_sed2)~
 + con(isnull(c_t_f_gt200),0,c_t_f_gt200)
/*

*******AML for Soil Creep*******

/* Estimate background sediment rate
/* First average out the slopes using the focalmean function
mean_slope = focalmean(slope,circle,1)
/* Determine soil fines fraction
rk_mn = polygrid(soil,rock.frag.min,#,#,32.823)
rk_mx = polygrid(soil,rock.frag.max,#,#,32.823)
rock_min = con(watershed_dem > 0,rk_mn,0)
rock_max = con(watershed_dem > 0,rk_mx,0)
finefrac = (100 - ((float(rock_max) + float(rock_min)) / 2)) / 100
gridclip finefrac fines box 1351927.517 811058.007 1367387.150 824055.915
/* Determine the soil depth (depth is in inches on attribute table)
soil_dmin = polygrid(soil,soil.depth.min,#,#,32.823)
soil_dmax = polygrid(soil,soil.depth.max,#,#,32.823)
soil_min = con(watershed_dem > 0,soil_dmin,0)
soil_max = con(watershed_dem > 0,soil_dmax,0)
soil_d = ((float(soil_min) + float(soil_max)) / 2) / 12
gridclip soil_d depth box 1351927.517 811058.007 1367387.150 824055.915
/*
strm_length = streams_shed * up
/* Find the background rate in US tons
creep = con(mean_slope < 30,strm_length * 2 * depth * .001 * 3937 / 1200 * fines,~
 strm_length * 2 * depth * .002 * 3937 / 1200 * fines)
/*

*******Accumulated Delivery*******

/* Accumulate it down the basin
hr1_accum = flowaccumulation(flowdir,hill_road_sed)
hr2_accum = flowaccumulation(flowdir,hr_sed2)