Urban Erosion Potential Risk Mapping with GIS

With increased regulatory focus on eroded sediment and its bound pollutants, methods are needed to predict areas with high erosive potential (EP) in urbanized areas. Using EP to prioritize urban areas for maintenance, implementation of Stormwater Control Measures (SCMs), stream restoration or monitoring is crucial. This study utilizes commonly available geospatial layers in conjunction with a computational procedure for prioritizing the contribution of site specificand transport-erosion to compute relative EP risk throughout a target urban watershed. Factors that contribute to erosion were evaluated: local cell slope, soil erodibility, land cover, runoff volume, distance and slope to nearest stormwater conveyance point along a surface flow travel path. A case study of the developed methodology was performed on a 1.6 square kilometer urban watershed in Blacksburg, VA, to generate EP risk maps. Results of the study indicate areas of erosive potential within the target watershed and provide a methodology for creating erosion potential risk maps for use by MS4 planners, engineers and other individuals that manage erosion control programs.


Introduction
Changes in biodiversity, hydrology and biochemistry of steams, referred collectively as the 'urban stream syndrome,' represent the negative consequences imparted by urbanization (Walsh et al., 2005).Observable effects in urban areas include modified landscapes and stream hydrology by way of sedimentation and erosion.Channel morphology and water quality are strongly influenced by sediment originating from erosion of bare land or streambanks, or during construction in urban areas (Hogan, Jarnagin, Loperfido, & Van Ness, 2014;Krug & Goddard, 1986;Paul & Meyer, 2001;"Sediment," n. d.).Both coarse and fine sediments are eroded and have various negative consequences downstream of the point of erosion (POE).Heavier particles directly contribute to sediment deposition along the travel path, typically closer to the POE.Fine sediments, which remain in suspension until reaching receiving waterbodies, block sunlight from traversing the water column to reach bottom-dwelling plants and can smother certain aquatic organisms in waterways (Chesapeake Bay Program, 2017;Reshetiloff, 2004;"Sediment," n. d.).Aggregation of this suspended sediment in channels affects the morphology and dimensions of streams, and can landlock ports (Reshetiloff, 2004;"Sediment," n. d.;Walsh et al., 2005).Further environmental degradation results from binding or absorption of nutrients and chemicals (toxins, metals, etc.) to the surface of sediment, allowing their transport and/or deposition with the sediment particles (Hogarth et al., 2004;Reshetiloff, 2004;"Sediment," n. d.).
Sediment is transported primarily via surface stormwater runoff to receiving waterways, where it often becomes a regulated pollutant.Municipal governments are required to manage stormwater to comply with federal and state regulations regarding water quality and quantity.The Clean Water Act (CWA 1972) initially placed controls on point source pollutant discharges, through permits via the National Pollutant Discharge Elimination System (NPDES).The CWA specifically targeted discharges from publicly owned treatment works (POTWs) (33 U.S.C §1251 et seq., 1972).The Water Quality Act (WQA) of 1987, an amendment to the CWA, worked to ensure water quality standards are met through monitoring and assessment.The CWA required the establishment of Total Maximum Daily Loads (TMDLs) for the Nation's waters, which guides restoration efforts to improve water quality to acceptable levels.Nutrients, sediments and metals are some common TMDLs that are often found in stormwater.Nonpoint source pollutants, including stormwater, contribute to 40% of the impaired waterways (Carle et al., 2005).However, definition of stormwater runoff as a point source through Municipal Separate Storm Sewer Systems (MS4s) has led to their inclusion in NPDES regulation (Clarke, 2003).
Erosion control regulations did exist prior to the WQA (1987); Franklin D. Roosevelt created the Soil Erosion Service in 1933 (now National Resources Conservation Service (NRCS), formerly Soil Conservation Service, under the Department of Agriculture) in response to the Dust Bowl, to conserve soil (Virginia Department of Environmental Quality, n. d.).Both Maryland and Virginia established laws for soil and water conservation in the 1930s.Maryland authorized sediment control regulations in 1957 and Virginia created the Virginia Erosion and Sediment Control Law (VESCL) in 1973 (Maryland Department of the Environment, n.d.;Virginia Department of Environmental Quality, n. d.).In state regulatory codes, erosion and sediment control and stormwater management may be addressed in separate sections.Both sediment and stormwater are considered pollutants, where sediment can be discharged with stormwater.The creation of MS4s allows for concerted efforts in addressing community-scale stormwater management through the NPDES permit system or applicable state-level versions of the NPDES (such as the Virginia Pollution Discharge Elimination System (VPDES)).To meet regulatory requirements, statewide and/or local stormwater management programs (SWMP) have been implemented to meet thresholds of pollutant discharges in stormwater, including sediment.
Development of local TMDLs for impaired waters, as required by the CWA, set some of the pollutant discharge thresholds that MS4s and municipalities are required to meet by their SWMP.The sensitivity of the Chesapeake Bay to nutrients and sediments has made these pollutants the focus of Virginia's nonpoint source pollution prevention program (Chesapeake Bay Program, 2017;Reshetiloff, 2004).To meet sediment TMDLs, the sediment load is often controlled through use of sediment and erosion control practices during construction, and with stormwater best management practices (BMPs) in other cases (Hogan et al., 2014).Many of these BMPs are coordinated through United States Department of Agriculture (USDA) local Soil and Water Conservation Districts (SWCDs) with a primary focus on agricultural erosion reduction.While the control of sediment runoff from agricultural practices is paramount in meeting sediment reduction goals, many agricultural practices are exempt from regulatory erosion control guidelines, thus limiting their regulation by local jurisdictions.Due to this restriction, it is crucial that MS4s can understand, predict, and optimally control sediment migration from regulated urban areas in order to make progress in meeting their TMDL goals.Typically, the highest probability of mass soil migration is during land disturbing activities where surface vegetation and topsoil has been stripped to allow for mass grading.Entrapment and removal of sediment prior to reaching downstream receiving waters may also remove bound pollutants, such as nutrients and chemicals, as well (Liu et al., 2003;Reshetiloff, 2004).

Background
Soil erosion occurs through two main processes: 1) detachment and 2) transport of soil particles during rainfall events (Young & Wiersma, 1973), resulting in impact (splash) erosion and transport erosion, respectively.Detachment, or splash erosion, which is caused by the kinetic energy of rainfall, which results in the dislodgment of a soil particle (Quansah, 1981;Young & Wiersma, 1973).Rainfall intensity is often seen as a main parameter in determining the erosivity of the rainfall.Studies show that soil erosion is more sensitive to changes in rainfall intensity versus changes strictly based on total rainfall volume, due to a larger transfer of energy over a shorter period of time (Nearing, 2001;Nearing et al., 2005).Slope steepness is another crucial factor affecting soil erosion, especially for flow (transport) erosion.Laboratory studies have found that slope steepness was positively correlated with both soil detachment and transport (Quansah, 1981).Erosion by rainfall impact and surface water flow, in conjunction with various other physical parameters, result in varying potential rates of soil migration throughout a watershed.Dominant downstream migration through a watershed is via rills; erosion on surfaces between rills is considered interrill erosion (Foster et al., 1977).These two components of erosion are driven by different detachment and transport processes.Rill erosion is dominated by concentrated flow, and interrill erosion is dominated by rainfall (raindrop impact and flow from rainfall) (G.Zhang, Liu, Liu, He, & Nearing, 2003).
The broad range of variables influencing soil erosion has produced several models and methods to model and predict soil loss or erosion potential (DeRoo et al., 1994;Downer et al., 2010;M. K. Jain et al., 2004;Leh et al., 2011;Morgan et al., 1984;Renard & Foster, 1991;Williams, 1982;X. Zhang & Zhang, 2012;S. Zhang et al., 2017).These models integrate the processes of detachment and transport to estimate soil loss from a watershed.Models are empirically, conceptually or physically based.The USLE and its derivations (RUSLE, MUSLE, etc.) are empirical models based on observational research (Merritt et al., 2003).Other common erosion models are physical (LISEM, WEPP), based in physical equations such as conservation of mass and momentum or sediment transport equations.Models can also be grouped by scale (catchment, hillslope, plot) or event type (long term or single event).They vary based on input, quality of input and processes modeled (Merritt et al., 2003).The Limburg Soil Erosion Model (LISEM) simulates an erosion process that falls between sheet and rill erosion using GIS maps and rainfall data (DeRoo et al., 1994).The Watershed Erosion Prediction Project (WEPP) focuses on sheet and rill erosion on hillslopes using physics based equations (Merritt et al., 2003).The Gridded Surface Subsurface Hydrologic Analysis (GSSHA) developed by the U.S Army Corps of Engineers uses detachment by rainfall and surface runoff, and runoff transport capacity to model overland soil erosion (Downer & Ogden, 2004;Downer et al., 2010).Models are limited by availability and quality of data and their outputs reflect those limitations and varying operational assumptions.
The most commonly used equation to estimate annual soil loss is through the Universal Soil Loss Equation (USLE).Its simplicity, relatively low data requirements, and generally wide application makes the USLE and its variations common choice for evaluation of soil loss.Developed by the USDA, the USLE predicts long-term annual soil loss by combined rill and interrill erosion on agricultural hillslopes (Haan et al., 1994;Merritt et al., 2003;Wischmeier & Smith, 1978).The USLE was intended for agricultural watersheds, but applications have been extended to forests, construction sites and mixed watersheds (Fraser et al., 1995;Mattheus & Norton, 2013;Renard & Foster, 1991).Refinements of the equation due to improving research over the years has resulted in the Modified USLE, MUSLE (Williams, 1982), and the Revised USLE, RUSLE (Renard & Foster, 1991), which make adjustments to the determination of the USLE parameters.
For all the models and their simplification of soil erosion processes, certain factors consistently seem to drive erosion.Detachment and transport erosion can both be influenced by rainfall, water flow, slope and soil type (Merritt et al., 2003;Quansah, 1981).Other factors including land cover/land use, rainfall interception, rainfall intensity and water regimes have effects on soil erosion processes (Leh et al., 2011;Morgan et al., 1984;Römkens et al., 2002;S. Zhang et al., 2017).Research has also examined the relationship between erosion and slope gradient and shear stress (Fox & Bryan, 2000;Onstad & Foster, 1975;G. Zhang et al., 2003).Although a variety of factors influencing erosion have been investigated, no universal value of influence for individual factors on soil erosion was found in the literature.Literature indicated that factors may vary by other preexisting conditions.

Study Objective
Recent studies apply soil erosion models on mostly agricultural watersheds (Balousek et al., 2000;Chang et al., 2015;S.K. Jain et al., 2001;Mattheus & Norton, 2013;Nearing et al., 2005;Pandey et al., 2007;Prasannakumar et al., 2012;Šurda et al., 2007).Input parameters vary based on the model chosen, but common factors that influence the total soil loss are topography, rainfall, slope, soil type and land cover (DeRoo et al., 1994;M. K. Jain et al., 2004;Leh et al., 2011;Renard & Foster, 1991;S. Zhang et al., 2017).USLE is considered one of the most robust models with few input parameters and has been applied to watershed types that extend beyond its intended application.However, other models (LISEM, SedNet, WEPP, etc) do exist that leverage datasets related to erosion that are now commonly available.Long term rainfall data, digital elevation models (DEMs), sediment delivery ratios and knowledge of plant coverage are additional parameters required by other soil erosion models (SedNet, WEPP etc).Results for existing models of rural watersheds are often soil erosion loss, annually or by single storm event, or as erosion potential (Merritt et al., 2003).
Despite this, almost no literature exists describing a simplified process for utilizing commonly available geospatial datasets for predicting erosion potential in urban areas.With increased focus and established thresholds on eroded sediment and its bound pollutants, caused by the effects of urbanization, methods are needed to predict land surface areas with high relative erosion risk in urban areas.MS4s, municipal planners, engineers, and state erosion control programs and researchers will benefit from determining locations of high risk areas.The objective of this study is to use commonly available geospatial layers to evaluate land surface erosion potential in urban areas.A procedure is developed to create a theoretical raster based GIS model of urban erosion potential and test it on a 1.6 square kilometer watershed case study in the Town of Blacksburg.The ability to prioritize urban areas for planning overlays, additional erosion control measures during construction, BMP implementation, restoration or monitoring is crucial to ensure that MS4s continue to maximize their efforts in meeting their community TMDL goals.

Methods
The Central Stroubles watershed, located in the Town of Blacksburg, Virginia was selected as a case study application of the methodology described in this section.Central Stroubles is a 1.6 square kilometer, primarily urban watershed located within the corporate boundaries of the south-west Virginian town.The watershed is a mix of commercial, civic, and residential areas with multiple land cover types and slope conditions.Notable areas of interest within the watershed include the eastern side of Main Street, the old Blacksburg High School campus, Westview Cemetery and Wong Park.

Data Sources
The methodology is based on tools (Spatial Analyst and Hydrology, specifically) available in ArcGIS (ESRI, 2016) and a Python script writer.Urban areas are moving towards documenting their stormwater assets (infrastructure and BMPs) on geographic information systems.The methodology assumes the watershed has stormwater infrastructure data, location of catch basins and other inlets, available in shapefiles.For the case study application, data used for the analysis was obtained from the Town's geodatabase ("Blacksburg GIS Database," 2015) or from online data repositories managed by the USDA and the U.S. Geological Survey (USGS).Specific datasets include: • A polyline feature class representing 0.6-meter (2 foot) contours for the Town of Blacksburg.
• A polygon feature class representing the Town of Blacksburg's Detailed Land Cover Database (DLCD).
• A polygon feature class generated using the USDA Web Soil Survey (WSS) for Montgomery County and the Soil Survey Geographic Database (SSURGO) ("Web Soil Survey," 2017).
• A line feature class representing National Hydrography Database (NHD) Flowlines from USGS

• High resolution Aerial Imagery of the Town of Blacksburg
The USDA WSS enables download of soil data for an area of interest.The information downloaded for Montgomery County, VA includes polygons of soil type, along with their corresponding map unit symbols and keys.Additional information about each map unit is found in the geodatabase included with the exported WSS download.

Parameter Characterization
A number of important parameters, as determined by literature review, are used to characterize soil erosion.These include: • Slope along surface flow travel path to nearest stormwater conveyance point These parameters are used to determine the erosion potential (EP) within a raster cell.The watershed is analyzed using a 10 x 10-meter raster grid.Each of the parameters of interest are ranked from 1 -10 for each 10-meter cell.For the target watershed, ranges corresponding to the various ranks are shown in Table 1.b Mass per area per erosivity unit (G.R. Foster et al., 1981) This analysis focuses on land surface erosion risk, so cells indicating stream locations are excluded from the raster image (and the analysis).These cells indicate areas of channelized flow and can skew the risk maps by presumed stream bank erosion.Streams within the watershed were estimated using a flow network.The flow network was identified using a non-weighted flow accumulation.For purposes of this study, cells exhibiting a DEM flow accumulation greater than 300 are removed from the analysis.The set of cells with accumulation values exceeding 300 generally correspond to stream locations based on visual examination of NHD flowlines and high-resolution aerial photos and contain the subset of points indicating channel cross-sections as catalogued in the Town of Blacksburg stormwater infrastructure geodatabase.The 0.6 meter (2 foot) contours for the Town of Blacksburg are converted in ArcGIS (ESRI 2016) to a digital elevation model (DEM) for the watershed extents using the default settings.Depressions within the DEM are filled using default settings of ArcGIS Hydrology Fill tool; filled version of the DEM is used throughout the analysis.
Rainfall intensity was determined to be an important factor in soil erosion by the literature review but is not used in the study.It is assumed that the rainfall intensity is uniform across the area of analysis.As the aim is to create a relative erosion risk map, the effect of rainfall intensity will cancel out as the value is the same for each cell.

Erosion Runoff Volume
Total runoff volume represents the accumulated water volume in a cell from all contributing upstream cells, including water from the cell in question.The NRCS runoff equation is used to calculate the volume of water generated by each contributing cell based on a total storm precipitation depth and the ArcMap Flow Accumulation function is used to aggregate upstream contributing runoff.The volume of runoff is influenced by curve numbers (CN) based on land cover and hydrologic soil group (HSG) of a cell.A runoff volume raster image was generated for the 2-year 24-hour storm (7.01 cm), based on NOAA Atlas 14 data for Blacksburg (3 SE 44-0766) ("NOAA Atlas 14," n.d.).The 2-year storm is the typical return period used to evaluate erosion in manmade channels in Virginia and many other states (New Jersey, Connecticut, etc.).Note that while the NOAA Atlas 14 dataset is sampled at a (~0.93 km latitude) resolution which would yield similar scale, but varying rainfall depths across the watershed, the storm total was kept constant using the NOAA precipitation depth value at the centroid of the watershed.This was done to preserve the assumption of homogeneity of rainfall across the watershed, thus preserving the ability to compute relative risk across watershed cells.
Curve numbers are determined by land cover and hydrologic soil group.The SSURGO polygons, with the HSG value field, are intersected with the DLCD layer to create a raster which identifies both land cover and HSG for a polygon.Using Urban Hydrology for Small Watershed Technical Release 55 (SCS, 1986), the following curve numbers (Table 2) were assigned to land cover types by hydrologic soil group; unless specified fair condition was assumed.The NRCS runoff equation is used to compute a frequency storm runoff volume for each cell (Equation 1).
where V is frequency storm runoff volume for a 10-m raster cell in cubic meters and Q m-NRCS is runoff depth from TR-55 in meters.
Once a frequency storm runoff volume raster has been established for the watershed, all of the contributing upstream cells need to be accounted for.Using the ArcToolbox Hydrology tool, a flow direction raster is created using the filled DEM.This flow direction raster is used to create two flow accumulations rasters.The first is an unweighted flow accumulation used for flow length and flow network calculations.The second is weighted by the frequency storm runoff volume raster.This second raster, representing cubic meters of accumulated runoff volume, is ranked from 1 -10 using Natural Jenks Breaks based on relative values (Table 1).

Local Cell Slope
The slope of the site (a raster cell) influences erosion potential and is modeled as the local slope parameter.The local slope (in degrees) is calculated from the filled DEM using the 3D Spatial Analyst tool.Local slope is ranked into 10 categories using the soil detachment rate presented by Zhang (2003) for shallow flow (Table 1).The equation provided by Zhang includes a flow rate variable along with the slope variable to estimate soil detachment.
Assuming a near unit flow rate (1 m/s), resulted in a range of soil detachment rate between 1 and 10 for slopes angles 1 to 45 degrees.Angles greater than 45 degrees are assigned a value of 10; values up to 45 degrees are good assumed angles of repose for soil (United States Department of Agriculture, 1994Agriculture, , 2007;;Upadhyaya, 2009).

Land Cover
The DLCD for the Town of Blacksburg classifies land covers for the town under twelve different types.These types were ranked on a scale from 1 -10 based on their erosive potential, where impervious surfaces are least erodible (1) and dirt is the most erodible (10) (Table 3).In the ranking of land cover types, rank 2 and 9 were left empty to allow for other land cover types.If a MS4s are using another land cover classification, such as the National Land Cover Database (NLCD), other land cover types may be used.Rank 2 was left empty to represent values that may exhibit low erosivity but not complete imperviousness (i.e.water based on NLCD) while rank 9 is highly erosive with some protection, such as agricultural land, based on NLCD cover types.

Soil Erodibility
Polygons imported from the WSS outline areas of specific soil types with few soil parameter attributes included with the feature class.Although this information is not included in the layer by default, both soil erodibility, K f , and the hydrologic soil group, necessary for the erosion potential calculation, can be found in the SSURGO database.The RUSLE2 Related Attributes Table, found in the SSURGO database, contains K f and HSG fields for all the soil types and is joined with the polygon layer using the map unit symbol/key (mukey).Udorthents and urban land soil types were assumed to have a soil erodibility of 0.27 (Renard, Foster, Weesies, McCool, & Yoder, 1997) as SSURGO places a "null" in this field; these values were manually updated.Polygons which represent impoundments were updated to have a soil erodibility of 0, as no surface soil erosion is expected within these polygons.For HSG, the worst-case scenario in terms of soil erosion assumes a soil group of "D" for urban land/udorthents and in any case where two HSGs are given (e.g.HSG B/D or C/D).A land cover type of "water" is also assumed to belong to HSG "D".
The soil polygons were converted to a raster image using the Polygon to Raster tool within ArcGIS using the soil erodibility, K f , values.The K f values are ranked on a scale of 1 -10 based on the range of soil erodibility values consistent with the USLE nomograph (Wischmeier & Smith, 1978).According to Wischmeier and Smith (1978), soil erodibility ranges from 0 to 0.092 ton hectacre hour/hectacre megajoule millimeter (0 and 0.7 ton acre hour/hundreds of acre-foot tonf inch), so the ranks 1 -10 as established by this procedure are defined as equal intervals spanning this erodibility range (Table 1).

Distance and Slope to Nearest Stormwater Conveyance Point
Once surface flow reaches a receiving channel, other factors such as bank stabilization, main channel flow, and other factors related to larger volumes of typically rapidly moving water begin to drive the erosion process.These factors, including stream bank erosion, are separate contributors to sediment loads that are not considered in this erosion potential calculation since the goal here is to find the highest EP risk in upslope areas.Due to this limitation, the distance to a receiving channel or point of entry to a conveyance network (storm sewer inlet) is considered important as it correlates with the risk that a particular cell will contribute to downslope erosion.
In urban areas, stormwater infrastructure typically conveys runoff directly to receiving channels rendering distances to these inlets important.For this study, nearest stormwater conveyance points include stormwater infrastructure (inlets, catch basins etc.) and points along a stream network.Due to the complexity of calculating a distance-to-outlet for each 10-meter raster cell, the distance and slope to nearest stormwater conveyance point was calculated using a Python script, Distance to Inlet (Python Software Foundation, n. d.).The Python script (described in detail in subsequent sections) evaluates the distance using a 'Near' Analysis on the points representing the nearest stormwater conveyance ("Near Analysis," 2017).Other factors such as elevation of the cell, rim elevation and flow lengths are considered when determining the most representative distance.This ensures that surface flow is not attributed to a nearest inlet that is upgrade from the elevation of the cell.
Developing Stormwater Conveyance Point Layer.The Town of Blacksburg's stormwater geodatabase contains a stormwater node point layer that includes catch basins, manholes, junctions, inlets and other related stormwater infrastructure nodes.The 'Node ID' and Rim Elevation are required inputs for the Python analysis.Node types were limited to nodes representing catch basins, headwalls/endwalls, cross sections and pond outlet structures as they represent inlets receiving stormwater runoff.The flow network was used to supplement the known infrastructure locations with stream networks based on upper and lower threshold network values.For the Central Stroubles case study watershed, the upper threshold was determined by where the stream entered an underground system in the aerial photograph.Those cells that fall within that range of network values were assigned elevations using the DEM and then converted to points and used as additional outlet points for the analysis.Infrastructure nodes and stream network points are combined into one stormwater conveyance point layer with two fields, Node ID and rim elevations.
Master Point Layer.The master point layer contains the fields: OBJECTID, Shape, pointid and grid_code, where pointid is the ID of the raster cell point and grid_code is the elevation of that cell point.Prior to running the Distance to Inlet script, the following fields are added to the master point layer, with pre-population of some fields (Table 4).

NEAR_NodeID
NodeID from the stormwater conveyance point layer that corresponds to the NEAR_FID as determined by a search cursor in the script.

NEAR_RimElev
The rim elevation corresponding to the NodeID in the stormwater conveyance point layer as determined by a search cursor in the script.

FlowLength
Flow length as determined by the flow length tool.Populated prior to running the script.
dElevToOutlet Change in elevation between the cell point and its outlet, as determined by the flow length weighted by decimal slope.Populated prior to running the script.

Avg_Slope
Average slope between the point and the nearest stormwater conveyance point, calculated as part of the script.
Flow Length Determination.The watershed boundary is infrastructure-corrected, taking into account storm inlets and pipes that convey stormwater outside of elevation-based drainage areas.Therefore, the watershed boundary does not necessarily follow the DEM which was created using surface contours.Due to these boundary deviations, for some raster cells, the nearest available stormwater conveyance may be outside of the boundaries of the watershed.Some of the points near the watershed boundary will therefore flow out of the watershed.To account for those situations, flow length for each cell was determined using the flow direction hydrology tool.The shortest distance, when comparing the NEAR distance and flow length, is used to represent the distance to nearest stormwater conveyance point.Because the average surface flow slope along this travel path (if flow length is the shortest distance) is a parameter of interest, the flow length weighted by decimal slope was also created.This weighted flow length raster represents the change in elevation between the current cell and its outlet, as determined by the flow length.If flow length is shorter, the near rim elevation is calculated as the difference between the cell elevation and this change in elevation.These flow lengths and change in elevation values are added to fields in the master point layer.
Python Script.The script uses the Near Analysis to determine the nearest distance between a master point layer, representing each 10-meter cell, and the stormwater conveyance points described above.The master point layer is updated with the Near analysis results, while a copy point layer is used in conjunction with the search cursor.
Iterating through each point in the watershed, the script steps through several layers to determine the nearest stormwater conveyance.The distance and average slope to the nearest node are returned from the script to be used as an erosion potential risk factor.The distance to the nearest node can be estimated using one of two methods, described above.The script compares the two distances, as flow length distance is pre-populated in the master file and the NEAR_Dist is determined during the current iteration (Figure 1).Final values for each parameter resulting from computations with the Python script were used in generation of raster grids for both distance to inlet and average slope.The distance to the nearest inlet was reclassified to rank 1 -10 based on Natural Jenks Breaks applied to the range of calculated values.Average slope along the travel path to the nearest inlet is ranked 1 -10 using the same scale as discussed previously for classification of local slope.For ranking thresholds, see Table 1.

Erosion Potential
Erosion potential (EP) of the watershed is evaluated using the following equation: where EP is erosion potential, V R is accumulated runoff volume, S L is local slope of the cell, LC is land cover, K f is soil erodibility, D is distance to inlet and S A is average slope to inlet.This equation is theoretical in nature based on the components of USLE that appeared to be valid for urban areas.
Erosion potential is divided into a site component (first term of equation 2) and a downstream transport component (second term).Lack of literature regarding the influence of these components on erosion potential prevent the assignment of relative weights to each.The individual parameters are equally weighted on a scale of 10.Each parameter has a different ranking scale; a single ranking could not be applied equally across all six parameters.Certain parameters (local slope, average slope, land cover, and soil erodibility) have ranking divisions independent of conditions in the watershed.Runoff volume and distance to nearest inlet are calculated on a watershed-by-watershed basis and their ranking divisions reflects that.The runoff volume and distance to nearest inlet ranking divisions are found using Natural Jenks breaks.The site component represents erosive processes occurring within the cell, influenced by the accumulated runoff volume, land cover, slope and soil erodibility.The distance to inlet and average slope to inlet are indicative of the potential of sediment to be transported from the cell and into a receiving channel (right side of Eq. 2).All factors are ranked into categories 1 (least erosive) to 10 (most erosive), as described in Table 1.Thus, the total EP is calculated as a unitless relative risk score.The resulting risk map generated from Equation 2is not intended to quantify the amount of sediment eroded but instead a relative risk of erosion, as compared to other locations within a target area.
Per the equation, transport erosion appears to have much less weight than the site component.Unranked EP scores can range from 2 to 10,100 with the transport erosion component contributing a maximum 100.Transport represents the potential for eroded sediment from a cell to move into a receiving channel.Once in a channel, sediment often becomes a pollutant controlled by federal, state and/or local regulations.While erosion can occur during transport, the site is assumed to contribute more to the overall EP.The site component focuses on land surface erosion and the four interacting factors that contribute to this surface erosion.Most of the erosion will occur here before being transported downstream; therefore, when equally weighted, the site component should have higher values than the transport component.Cases where the transport component contributes a higher score to the EP value are discussed after application of the methodology to the case study watershed.

Erosion Potential Risk Map
Erosion potential is calculated for each cell in the watershed based on Equation 2. Relative erosive potential risk scores range from 2 (low) to 618 (high) for the watershed.Results are classified into 10 bins by Natural Jenks breaks based on the range of scores from the calculated EP (Figure 2).Breaks determined by Natural Jenks are used throughout this analysis because of its ability to minimize deviation within a class while maximizing deviations between classes.By this, the values within a bin are more similar than the values outside a bin.The Natural Jenks breaks aim to accurately represent the data's spatial attributes.Erosion potential risk maps are the desired output of the model; a classification method that creates bins that are similar within and different between, and also enhances the data's spatial attributes is an ideal choice.
All of the maps in Figure 2 are classified using the same scale.Classification is performed using Natural Jenks rather than equal intervals to show more visual distinction of EP components.When using equal intervals of 60, the bins generated result in 99% of the downstream transport values falling in the lowest erosive potential category and all the values within the first two bins (Table 5).With classification according to the Natural Jenks breaks, the transport values fall within the first four bins.The site (cell) specific component has 75% of the values in the first bin using equal intervals; the same approximate percentage falls within the first three bins using Natural Jenks.Even using Natural Jenks, the range of values for downstream transport potential fall on the lower end of the ranking as the highest possible value of 100 falls in the fourth lowest bin.For the total EP and site component, values encompass the entire range.Looking at the EP risk map for the Central Stroubles watershed, the majority of the watershed falls in the lower risk category.The median EP score is 38 with a median site component value of 30 and transport component of 8.
A quick visual inspection of points across the watershed confirm that there are few visible signs of erosion.Ten sites with EP scores of 15 to 155 were chosen for a field check.The likelihood of finding major signs of erosion was considered small since active maintenance by landowners and/or the local government can often reduce or eliminate erosion completely (Figure 3).Of the sites chosen, only one showed visible signs of erosion.The remaining 9 locations were stabilized grass (lawn) or brush/bush surfaces, well maintained by the property owner.The potential for erosion, though, exists whether the location does or doesn't show visible signs of erosion.The EP risk map does not quantify sediment erosion or determine if erosion is occurring, simply the relative risk of erosion when compared to other cells within the analysis area.

Influence of EP Components
The potential for erosion is divided into two components: site (cell) erosion and downstream transport potential.
Although the components are equally weighted, site erosion carries more importance than downstream transport, per equation 2, based on its number of influencing parameters.Site erosion has a maximum theoretical value of 10, 000 and transport erosion has a maximum theoretical value of 100.Visually, the site component features a distribution more similar to the total EP risk map, indicating its influence on the total value.Certain areas within the watershed have the transport component contributing a majority of the score to the overall EP risk (Figure 4).
In approximately 14% of the cells within the watershed, the transport component contributes 50% or more of the score towards total EP.75% of the cells have up to 33% of the EP as transport potential, and 25% of cells contribute less than 10% from transport potential towards the EP.The transport component often has less influence than the site component on total EP values, but in certain areas it is the driver of the EP value.These areas that are highly influenced by transport generally fall in the lower EP risk categories, based on Figure 2, above.Impervious surfaces are prominent in the lower southwestern section of the watershed (downtown Blacksburg) and the northeastern corner (old Blacksburg High School campus).Both sections have several stormwater inlets within a small area.While the transport component has a median contribution of 18% towards the total EP value, it has enough influence on heavily impervious areas of the watershed that it warrants inclusion in Equation 2 despite its maximum theoretical contributing score is 100 times smaller than site erosion.

Adaptations to Methodology
In its current form, Equation 2is purely theoretical in nature with equal weights applied to all parameters, which yields an unbalanced influence between the local and transport components as previously described.Adaptations may need to be made to the methodology when applying the method to a larger or smaller scale than the target watershed in this study.Because of the relative rankings used for the 2-yr runoff volume and distance to inlet, neighboring watersheds cannot be directly compared unless they use a common ranking for each of those parameters.For MS4s who want to compare erosive potential across the community, absolute rankings must be established from aggregated data for runoff volume and distance to inlet, prior to creation of any risk maps.
To better represent erosion potential, further information should be gathered to determine if any of the six parameters or EP components need to be individually weighted.In this study, all parameters are equally weighted, as are the site and downstream transport components of EP.The current form of Equation 2 has the transport component delivering a maximum value of 100 towards the total EP, while the local component delivers a maximum value of 10,000.It is possible that this weight should be adjusted to represent the influence of transport on erosion potential through the contribution of sediment to receiving channels.In addition, within the site component, certain parameters, such as local slope or land cover, may have greater influence on surface erosion, and could be adjusted accordingly.No research in the literature was found to propose or justify any non-uniform weighting of the parameters used in this study; therefore, all parameter weightings remain equal.Because erosive risk cannot be quantified by field observation or measurement, absolute risk is unknown, which prevents use of a sensitivity analysis to determine potential non-equal weights of the parameters used in Equation 2.

Conclusions
This study proposes geospatial model for estimating erosive potential in urban areas using commonly available geospatial layers and Python scripting.A literature review brought forward many influential factors in determining land surface erosion.Six major components were selected as the basis for the proposed EP calculation.EP is broken up into site specific erosive potential and downstream transport potential components.This benefit of this model is that components making up site erosion and transport erosion potentials are focused on urban conditions.Distance to stormwater infrastructure is a component unique to urban areas that is accounted for in this EP model.Land cover, while not unique to urban areas, has a significant impact on the erosion potential when considering impervious surfaces.By focusing on urban areas, parameters were modeled using data commonly available within these regions.Data can be sourced from national websites, including USDA's WSS and USGS's The National Map viewer.Urban areas, particularly those under NPDES and MS4s, have mapped stormwater infrastructure that can easily be modeled in GIS if not already.
Application of the methodology to the Central Stroubles watershed in Blacksburg, Virginia indicates low erosive potential across majority of the watershed.Of the two components, site specific erosion potential seems to have the most influence on total EP.Within the Central Stroubles watershed, transport potential had high contribution to total EP when impervious cover was high and distance to inlet shorter.Visual inspection of ten locations across the watershed confirm little to no visible signs of erosive damage, possibly due to maintenance by landowners which can curtail most signs of impending surface erosion; therefore, visual inspection by itself is an insufficient method for assessing or confirming EP.
Erosion potential risk maps provide a tool for MS4s and engineers, giving them the ability to focus attention and resources on areas of an urban watershed that may be contributing to sediment concentrations in local channels.
Particularly, it highlights areas of extreme erosive risk that may need to be managed with larger and more efficient erosion control measures during land disturbance.While traditional zoning restrictions have typically been limited to slope overlays to aid in controlling erosion, the methodology used within this study provides many additional contributing parameters beyond slope that can be used by future planners in the preparation of overlay mapping focused on preserving and improving the quality of the communities receiving waters.

Figure 1 .
Figure 1.A sample Python code segment used for determining the nearest node distance based on one of two methods

Figure 3 .
Figure 3. Areas visually inspected for visible signs of erosion within the watershed

Figure 4 .
Figure 4.The percentage contribution of the transport component toward the combined EP

Table 1 .
Summary table of established ranking thresholds 1 (least erosive) to 10 (most erosive) for the six parameters used to determine final erosion potential values

Table 2 .
DLCD land cover types paired with a TR-55 cover type and their corresponding curve numbers

Table 3 .
Land covers organized by their erosive potential on a scale of 1 to 10.All impervious surfaces are placed into the same erosive category

Table 4 .
Field names and description for the master point layer of the Python script calculating distance and slope to nearest inlet

Table 5 .
Frequency distributions of two different classification methods, equal interval and Natural Jenks, for final erosion potential of Central Stroubles