Mass Wasting AML

Hoodsport Watershed Analysis

/* Mass Wasting AML
/* Created by Susan Seaberg 3/1/99

*******AML for Infinite Hillslope Equation (FS)*******

/* created soil grid, soilgd, using ArcView
/* isolated the soils within the watershed:  soils = soilgd * watershed_1
/*
/* Make a grid of soil depth; soils coverage contains min and max depths in inches then
/* change to feet
soil_inch_min = polygrid(soil,soil.depth.min,#,#,32.823)
gridclip soil_inch_min min_s_d box 1351927.517 811058.007 1367387.150 824055.915
min_soil = setnull(soils == 9,min_s_d / 12)
soil_inch_max = polygrid(soil,soil.depth.max,#,#,32.823)
gridclip soil_inch_max max_s_d box 1351927.517 811058.007 1367387.150 824055.915
max_soil = setnull(soils == 9,max_s_d / 12)
/*
/* Create a grid from riu of mbf/acre.  Located in
/* fiu.tot.tim.sv6g or Board ft volume to a 6" top
/* per acre and then multiply by 2(from lab sheet?)
/* and then multiply by density.  I will set this at 
/* .45 for douglas firs
vegetation = polygrid(riu,fiu.tot.tim.sv6g,#,#,32.823)
gridclip vegetation mbf_acre box 1351927.517 811058.007 1367387.150 824055.915
/*
&sv tree_wt = 31100
veg_weight = con(isnull(mbf_acre),((0 + (rand() * %tree_wt%)) * 2 * .45),2 * mbf_acre * .45)
/* 
/* Make the slopes
degree = slope(watershed_dem,degree)
/* 
/* Make a grid of minimum and maximum tan(phi)
minphi_rad = float(soils.phi_min) * pi / 180
maxphi_rad = float(soils.phi_max) * pi / 180
/*
/* Make a grid of cos(theta)and sin(theta)
cos_angle = cos(degree * pi / 180)
sin_angle = sin(degree * pi / 180)
/*

/* FS minimums

/* Use Cr, root reinforcement = 50; using saturation as 0
fs_min1 = (soils.c_min + 50 + (tan(minphi_rad) * (sqr(cos_angle))) *~
 (veg_weight + (min_soil * soils.u_wt_sat_min) - (62.4 * 0)))~
 / ((veg_weight + (min_soil * soils.u_wt_sat_min)) * (sin_angle) * (cos_angle))
fs_11 = con(fs_min1 <= 1,1,0)
/*
/* Use Cr,root reinforcement = 50; use saturation as total soil thickness
fs_min2 = (soils.c_min + 50 + (tan(minphi_rad) * (sqr(cos_angle))) *~
 (veg_weight + (min_soil * soils.u_wt_sat_min) - (62.4 * min_soil)))~
 / ((veg_weight + (min_soil * soils.u_wt_sat_min)) * (sin_angle) * (cos_angle))
fs_12 = con(fs_min2 <= 1,1,0)
/*
/* Use Cr, root reinforcement = 250; using saturation as 0
fs_min5 = (soils.c_min + 250 + (tan(minphi_rad) * (sqr(cos_angle))) *~
 (veg_weight + (min_soil * soils.u_wt_sat_min) - (62.4 * 0)))~
 / ((veg_weight + (min_soil * soils.u_wt_sat_min)) * (sin_angle) * (cos_angle))
fs_15 = con(fs_min5 <= 1,1,0)
/*
/* Use Cr,root reinforcement = 250; use saturation as total soil thickness
fs_min6 = (soils.c_min + 250 + (tan(minphi_rad) * (sqr(cos_angle))) *~
 (veg_weight + (min_soil * soils.u_wt_sat_min) - (62.4 * min_soil)))~
 / ((veg_weight + (min_soil * soils.u_wt_sat_min)) * (sin_angle) * (cos_angle))
fs_16 = con(fs_min6 <= 1,1,0)
/*

/* FS maximums

/* Use Cr, root reinforcement = 250; use saturation as 0
fs_max1 = (soils.c_max + 250 + (tan(maxphi_rad) * (sqr(cos_angle))) *~
 (veg_weight + (max_soil * soils.u_wt_sat_max) - (62.4 * 0)))~
 / ((veg_weight + (max_soil * soils.u_wt_sat_max)) * (sin_angle) * (cos_angle))
fs_21 = con(fs_max1 <= 1,1,0)
/*
/* Use Cr, root reinforcement = 250; use saturation as total soil thickness
fs_max2 = (soils.c_max + 250 + (tan(maxphi_rad) * (sqr(cos_angle))) *~
 (veg_weight + (max_soil * soils.u_wt_sat_max) - (62.4 * max_soil)))~
 / ((veg_weight + (max_soil * soils.u_wt_sat_max)) * (sin_angle) * (cos_angle))
fs_22 = con(fs_max2 <= 1,1,0)
/*
/* Use Cr, root reinforcement = 50; use saturation as 0
fs_max5 = (soils.c_max + 250 + (tan(maxphi_rad) * (sqr(cos_angle))) *~
 (veg_weight + (max_soil * soils.u_wt_sat_max) - (62.4 * 0)))~
 / ((veg_weight + (max_soil * soils.u_wt_sat_max)) * (sin_angle) * (cos_angle))
fs_25 = con(fs_max5 <= 1,1,0)
/*
/* Use Cr, root reinforcement = 50; use saturation as total soil thickness
fs_max6 = (soils.c_max + 250 + (tan(maxphi_rad) * (sqr(cos_angle))) *~
 (veg_weight + (max_soil * soils.u_wt_sat_max) - (62.4 * max_soil)))~
 / ((veg_weight + (max_soil * soils.u_wt_sat_max)) * (sin_angle) * (cos_angle))
fs_26 = con(fs_max6 <= 1,1,0)
/*
/* FS max and min using average soil cohesion
/* Use Cr, root reinforcement = 50; using saturation as 0
fs_min3 = (((soils.c_min + soils.c_max) / 2 ) + 50 + (tan(minphi_rad) * (sqr(cos_angle))) *~
 (veg_weight + (min_soil * soils.u_wt_sat_min) - (62.4 * 0)))~
 / ((veg_weight + (min_soil * soils.u_wt_sat_min)) * (sin_angle) * (cos_angle))
fs_13 = con(fs_min3 <= 1,1,0)
/*
/* Use Cr,root reinforcement = 50; use saturation as total soil thickness
fs_min4 = (((soils.c_min + soils.c_max) / 2 ) + 50 + (tan(minphi_rad) * (sqr(cos_angle))) *~
 (veg_weight + (min_soil * soils.u_wt_sat_min) - (62.4 * min_soil)))~
 / ((veg_weight + (min_soil * soils.u_wt_sat_min)) * (sin_angle) * (cos_angle))
fs_14 = con(fs_min4 <= 1,1,0)
/*
/* Use Cr, root reinforcement = 250; use saturation as 0
fs_max3 = (((soils.c_min + soils.c_max) / 2 ) + 250 + (tan(maxphi_rad) * (sqr(cos_angle))) *~
 (veg_weight + (max_soil * soils.u_wt_sat_max) - (62.4 * 0)))~
 / ((veg_weight + (max_soil * soils.u_wt_sat_max)) * (sin_angle) * (cos_angle))
fs_23 = con(fs_max3 <= 1,1,0)
/*
/* Use Cr, root reinforcement = 250; use saturation as total soil thickness
fs_max4 = (((soils.c_min + soils.c_max) / 2 ) + 250 + (tan(maxphi_rad) * (sqr(cos_angle))) *~
 (veg_weight + (max_soil * soils.u_wt_sat_max) - (62.4 * max_soil)))~
 / ((veg_weight + (max_soil * soils.u_wt_sat_max)) * (sin_angle) * (cos_angle))
fs_24 = con(fs_max4 <= 1,1,0)
/*

/* FS grid using Random Variables

setwindow watershed_dem
setcell 32.283
rand_cr = 50 + (rand() * (250 - 50))
rand_c = soils.c_min + (rand() * (soils.c_max - soils.c_min))
rand_veg = rand() * veg_weight
rand_depth = min_soil + (rand() * (max_soil - min_soil))
rand_satwt = soils.u_wt_sat_min + (rand() * (soils.u_wt_sat_max - soils.u_wt_sat_min))
rand_z = 0 + (rand() * (max_soil))
rand_phi = minphi_rad + (rand() * (maxphi_rad - minphi_rad))
/*
fs_rand1 =  (rand_c + rand_cr + (tan(rand_phi) * (sqr(cos_angle))) *~
 (rand_veg + (rand_depth * rand_satwt) - (62.4 * rand_z)))~
 / ((rand_veg + (rand_depth * rand_satwt)) * (sin_angle) * (cos_angle))
fs_r1 = con(fs_rand1 <= 1,1,0)
/*
setwindow watershed_dem
setcell 32.823
rand_cr1 = 50 + (rand() * (250 - 50))
rand_c1 = soils.c_min + (rand() * (soils.c_max - soils.c_min))
rand_veg1 = rand() * veg_weight
rand_depth1 = min_soil + (rand() * (max_soil - min_soil))
rand_satwt1 = soils.u_wt_sat_min + (rand() * (soils.u_wt_sat_max - soils.u_wt_sat_min))
rand_z1 = 0 + (rand() * (max_soil))
rand_phi1 = minphi_rad + (rand() * (maxphi_rad - minphi_rad))
/*
fs_rand2 =  (rand_c1 + rand_cr1 + (tan(rand_phi1) * (sqr(cos_angle))) *~
 (rand_veg1 + (rand_depth1 * rand_satwt1) - (62.4 * rand_z1)))~
 / ((rand_veg1 + (rand_depth1 * rand_satwt1)) * (sin_angle) * (cos_angle))
fs_r2 = con(fs_rand2 <= 1,1,0)
/*
setwindow watershed_dem
setcell 32.823
rand_cr2 = 50 + (rand() * (250 - 50))
rand_c2 = soils.c_min + (rand() * (soils.c_max - soils.c_min))
rand_veg2 = rand() * veg_weight
rand_depth2 = min_soil + (rand() * (max_soil - min_soil))
rand_satwt2 = soils.u_wt_sat_min + (rand() * (soils.u_wt_sat_max - soils.u_wt_sat_min))
rand_z2 = 0 + (rand() * (max_soil))
rand_phi2 = minphi_rad + (rand() * (maxphi_rad - minphi_rad))
/*
fs_rand3 =  (rand_c2 + rand_cr2 + (tan(rand_phi2) * (sqr(cos_angle))) *~
 (rand_veg2 + (rand_depth2 * rand_satwt2) - (62.4 * rand_z2)))~
 / ((rand_veg2 + (rand_depth2 * rand_satwt2)) * (sin_angle) * (cos_angle))
fs_r3 = con(fs_rand3 <= 1,1,0)
/*
setwindow watershed_dem
setcell 32.823
rand_cr3 = 50 + (rand() * (250 - 50))
rand_c3 = soils.c_min + (rand() * (soils.c_max - soils.c_min))
rand_veg3 = rand() * veg_weight
rand_depth3 = min_soil + (rand() * (max_soil - min_soil))
rand_satwt3 = soils.u_wt_sat_min + (rand() * (soils.u_wt_sat_max - soils.u_wt_sat_min))
rand_z3 = 0 + (rand() * (max_soil))
rand_phi3 = minphi_rad + (rand() * (maxphi_rad - minphi_rad))
/*
fs_rand4 =  (rand_c3 + rand_cr3 + (tan(rand_phi3) * (sqr(cos_angle))) *~
 (rand_veg3 + (rand_depth3 * rand_satwt3) - (62.4 * rand_z3)))~
 / ((rand_veg3 + (rand_depth3 * rand_satwt3)) * (sin_angle) * (cos_angle))
fs_r4 = con(fs_rand4 <= 1,1,0)
/*
setwindow watershed_dem
setcell 32.823
rand_cr4 = 50 + (rand() * (250 - 50))
rand_c4 = soils.c_min + (rand() * (soils.c_max - soils.c_min))
rand_veg4 = rand() * veg_weight
rand_depth4 = min_soil + (rand() * (max_soil - min_soil))
rand_satwt4 = soils.u_wt_sat_min + (rand() * (soils.u_wt_sat_max - soils.u_wt_sat_min))
rand_z4 = 0 + (rand() * (max_soil))
rand_phi4 = minphi_rad + (rand() * (maxphi_rad - minphi_rad))
/*
fs_rand5 =  (rand_c4 + rand_cr4 + (tan(rand_phi4) * (sqr(cos_angle))) *~
 (rand_veg4 + (rand_depth4 * rand_satwt4) - (62.4 * rand_z4)))~
 / ((rand_veg4 + (rand_depth4 * rand_satwt4)) * (sin_angle) * (cos_angle))
fs_r5 = con(fs_rand5 <= 1,1,0)
/*
setwindow watershed_dem
setcell 32.823
rand_cr5 = 50 + (rand() * (250 - 50))
rand_c5 = soils.c_min + (rand() * (soils.c_max - soils.c_min))
rand_veg5 = rand() * veg_weight
rand_depth5 = min_soil + (rand() * (max_soil - min_soil))
rand_satwt5 = soils.u_wt_sat_min + (rand() * (soils.u_wt_sat_max - soils.u_wt_sat_min))
rand_z5 = 0 + (rand() * (max_soil))
rand_phi5 = minphi_rad + (rand() * (maxphi_rad - minphi_rad))
/*
fs_rand6 =  (rand_c5 + rand_cr5 + (tan(rand_phi5) * (sqr(cos_angle))) *~
 (rand_veg5 + (rand_depth5 * rand_satwt5) - (62.4 * rand_z5)))~
 / ((rand_veg5 + (rand_depth5 * rand_satwt5)) * (sin_angle) * (cos_angle))
fs_r6 = con(fs_rand6 <= 1,1,0)
/*
/*

/* Average of the FS grids

/* FS INTEGER AVERAGE
fs_ave = (fs_11 + fs_12 + fs_13 + fs_14 + fs_15 + fs_16 + fs_21 +~
 fs_22 + fs_23 + fs_24 + fs_25 + fs_26 + fs_r1 + fs_r2 + fs_r3 + fs_r4 + fs_r5 +~
 fs_r6) / 18
/*
/* FS FLOATING POINT AVERAGE
fs_ave2 = (float(fs_11) + float(fs_12) + float(fs_13) + float(fs_14) + float(fs_15)~
 + float(fs_16) + float(fs_21) + float(fs_22) + float(fs_23) + float(fs_24) + float(fs_25)~
 + float(fs_26) + float(fs_r1) + float(fs_r2) + float(fs_r3) + float(fs_r4) + float(fs_r5)~
 + float(fs_r6)) / 18
/* 
/*

/* Determine Probability of Failure based on Critical Soil Thickness

/* Soil creep rate (q) = 16 in^2 / year;  change to ft^2
/* equation:  q * total local curvature  * stand age / hc
/* Hc = (Cs + Cr)/((Ysat * sin(theta) * cos(theta))-((Ysat - Yw)*tan(phi)*cos^2(theta)))
/*
/* We are not getting the critical thickness for failure.....later try and play with the 
/* numbers in the hc grid
hc_min = ((float(soils.c_min) + 50) / ((float(soils.u_wt_sat_min) * sin_angle *~
 cos_angle) - ((float(soils.u_wt_sat_min) - 62.4) * (tan(minphi_rad) * (sqr(cos_angle))))))
hc_max = ((float(soils.c_max) + 50) / ((float(soils.u_wt_sat_max) * sin_angle *~
 cos_angle) - ((float(soils.u_wt_sat_max) - 62.4) * (tan(maxphi_rad) * (sqr(cos_angle))))))
/* need a curvature grid
curve = curvature(watershed_dem)
/* need grid of stand ages
ages = polygrid(riu,fiu.pri.age,#,#,32.823)
gridclip ages tree_age box 1351927.517 811058.007 1367387.150 824055.915
trees = con(isnull(tree_age),(0 + (rand() * 77)),tree_age)
pfail_min = ((16 / (12 * 12)) * curve * trees) / hc_min
pfail_1 = con(pfail_min <= 0,0,pfail_min <= 1,pfail_min)
pfail_max = ((16 / (12 * 12)) * curve * trees) / hc_max
pfail_2 = con(pfail_max <= 0,0,pfail_max < 1,pfail_max)
/*

/* Delivery to streams

/* 
/* Slides will stop if it reaches a slope that is 1/2 as steep as its internal friction angle
/* Make a grid of average friction angles (max + min)/2
half_iaf = con(soils.phi_min > 0,(((soils.phi_max + soils.phi_min) / 2) / 2))
/* Identify the slopes that equal these areas
if (degree <= half_iaf)
stop_slopes = 1
else
stop_slopes = 0
endif
/*

/* Find streams vs stop slopes

flowdir = flowdirection(watershed_dem,channel_grade)
del_streams = con(streams_shed == 1 and stop_slopes == 0,1)
nonstop_area = watershed(flowdir,del_streams)
delivery_area = con(isnull(nonstop_area),0,nonstop_area)
/* Modify hazard and probability maps
mw_delivery = fs_ave2 * delivery_area
mw_accum = flowaccumulation(flowdir,mw_delivery)
/*
/*

******* Shaw Johnson*******

pslope = slope(watershed_dem,percentrise)
sj_slope = con(slope <= 15,1,slope <= 25,2,slope <= 47,3,slope <= 70,4,slope > 70,5)
sj_curve = con(curve < 0,1,curve == 0,2,curve > 0,3)
sj_green = con(sj_slope == 1 or sj_slope == 2 and sj_curve == 1 or sj_slope == 2 and sj_curve == 2~
 or sj_slope = 3 and sj_curve == 1 or sj_slope == 3 and sj_curve == 2~
 or sj_slope == 4 and sj_curve == 1,1)
sj_yellow = con(sj_slope == 2 and sj_curve == 3 or sj_slope == 4 and sj_curve == 2~
 or sj_slope == 5 and sj_curve == 1,2)
sj_red = con(sj_slope == 3 and sj_curve == 3 or sj_slope == 4 and sj_curve == 3~
 or sj_slope == 5 and sj_curve == 2 or sj_slope == 5 and sj_curve == 3,3)