Assessing Watershed Vulnerability in Bernalillo County, New Mexico Using GIS-Based Fuzzy Inference ()
1. Introduction
Urban sprawl is a leading cause of non-point source (NPS) pollution in the world [1] , with urban runoff being a major contributor to elevated levels of harmful pollutants in streams, lakes, rivers, and oceans. Common runoff pollutants include sediments, nutrients, organic materials, pathogens, hydrocarbons, pesticides, metals, chlorides, grease, trash, and toxic substances. These pollutants pose a serious threat to surface water quality and sustainable development. Additionally, the impact of runoff is much more realized in arid watersheds that are often characterized by scant but high-intensity episodic rainfall events with increased runoff volume [2] . A search for stormwater management practices to reduce stormwater runoff volume and improve overall runoff quality is a concern shared by many community planners and local governments for all land use classifications (urban and non-urban) within their jurisdictional boundary.
To address this concern, scientists, planners, and stakeholders are usually required to evaluate several competing watershed-specific factors in order to develop procedures and practices that can best minimize the issues associated with stormwater runoff. Common watershed factors include topography, precipitation, soil, vegetation, and drainage area. In most cases, selecting the best combination of methods and alternative conservation practices in the presence of competing factors poses a great challenge for decision makers [3] . A responsible approach to watershed protection is by adopting Best Management Practices (BMPs) [4] , also known as Stormwater Control Measures (SCMs).
Structural SCMs and non-structural SCMs involve techniques, methods, activities, and practices that have the primary objective of preventing or reducing the quantity of pollutants present in runoff water before reaching natural water systems. The increase in SCM awareness has been partly attributed to their capability to control runoff volume while preserving existing hydrological function [4] , and also by the number of stringent laws imposed by the federal government [5] . Many implementations of SCMs have been established in literature, but there is still no single SCM that can exclusively solve every watershed problem [3] . This implies that a good implementation of SCM strategies should be based on project and watershed-specific needs. Gautam et al. [6] indicated that due to rapid urbanization and limited water availability in the desert Southwest, the best implementation of SCMs for such areas should be based on water conservation and water reuse strategies.
One extensively used technique for improving SCM efficiency is by incorporating multi-criteria decision support (MCDS) systems. MCDS involves the use of scientific and mathematical tools to evaluate complex trade-offs between different environmental and socio-economic factors [2] . Commonly adopted MCDS techniques include analytic hierarchy process (AHP) and fuzzy inference system (FIS). The former is a measurement method to derive a scale of relative importance of alternatives or criteria from pair-wise comparisons [7] . Input can be actual measurements or subjective judgement. The latter uses fuzzy set theory [8] to map inputs to outputs. A more evolving concept is the integration of MCDS and GIS processes into SCM suitability studies. One potential benefit of GIS-based MCDS models is that they can be easily implemented in regions where precise spatial data are not readily available without having to rely on expensive and time-consuming field surveys. Information from MCDS models can prove useful for overall watershed management and resource planning.
This study combines specific soil erosion and infiltration criteria within a GIS-based FIS framework using Fuzzycell software to develop a watershed vulnerability map for Bernalillo County, New Mexico. The area is characterized by intermittent precipitation and limited water availability. Model output could be beneficial for exploring structural and non-structural SCM approaches to watershed management for targeted vulnerable areas within the county. The premise herein is that areas with high soil erosion potential and low infiltration potential are most vulnerable. An AHP based vulnerability map was developed concurrently [9] [10] .
2. SCMs in the Desert Southwest
Irrespective of the wide acceptance of SCM approaches in stormwater control projects, there has been little representation of such approaches in the desert Southwest [6] . This lack of representation in SCMs has often been attributed to the perception that these regions experience little or no rainfall and hence the expense and burden of SCMs are unjustifiable. Contrary to such belief, Gautam et al. [6] have argued that this perception is false and there is much need of SCM approaches for such areas. One important reason for this argument is that the limited but high-intensity episodic rainfall experienced in this region is a major contributing factor to the elevated levels of harmful pollutants found in natural watercourses. Another contention is that this region often suffers from water availability and sustainability problems, which is an indicator for a greater need of SCMs. Finally, Gautam et al. [6] suggest that, due to the rapid urbanization and severe water availability issues in the desert Southwest, the best strategy for implementing SCMs in such areas should be based on sustainable landuse, water conservation, and water reuse techniques. The authors note that one critical SCM for long-term sustainability is through low-impact development by reducing impervious areas and limiting the directly connected impervious surfaces. Water harvesting is another sustainable practice to address water quantity and quality in an integrated manner.
3. Economic Justification of SCMs
Use of economic argument to justify the expense and burden of SCM projects is a difficult task. One reason for this perception can be attributed to the fact that direct economic benefits of SCM projects are not readily quantifiable. Contrary to this belief, a study of the economic benefits of conservation practices in five states (Florida, Nebraska, New Mexico, Maine and Oregon) indicated a direct relationship between economic growth and watershed health [11] . Maintaining a clean and healthy watershed can greatly reduce healthcare cost and the risk of exposure to harmful chemicals and pollutants [11] . Additionally, maintaining a forested watershed has been shown to greatly reduce the capital investment and operational and maintenance cost of drinking water treatment facilities [12] . Considering watershed protection as an asset and a national capital could provide numerous health, social, and economic incentives over time.
4. MCDS Approach to SCM and Watershed Vulnerability Using FIS
Fuzzy Inference Systems (FIS)
Fuzzy inference systems (FIS) uses fuzzy set theory to classify data based on the degree of memberships. In classical set theory, an element can either exist as a member of a set (1) or a non-member (0). Contrary to classical set theory, fuzzy set theory classifies data based on the degree of memberships ranging from 0 to 1 (Figure 1). The ability to classify environmental factors to reflect spatial variability makes fuzzy logic a suitable system for watershed studies. Another important component of FIS is the ability to incorporate human thinking and linguistic variables into fuzzy membership functions. Incorporating human knowledge is an important step in recognizing the gray areas in memberships and minimizing the errors associated with modelling uncertainties in environmental data [13] .
Numerous attempts have been made to investigate the significance of fuzzy logic in watershed studies. Previous approaches, such as that by Nielsen and Hjelmfelt [14] , have focused on the placement of soil into hydrologic soil groups. The study concluded that fuzzy logic can be an effective technique in predicting hydrologic soil groups that correlates well with existing placements. Several authors have also conclusively shown that a fuzzy logic approach can be an in-expensive method of estimating soil erosion risk at considerable accuracies [15] [16] [17] . The need to incorporate water quality data, runoff data, and landuse information into fuzzy-based soil erosion risk models have been addressed by Black [18] . Other studies have highlighted the importance of using rule-based decision approach in developing groundwater susceptibility and contamination maps [19] [20] .
In recent years, there has been an increasing amount of literature on fuzzy based water conservation techniques. One example of such approaches is that by Shamsudin et al. [21] , involving the use of fuzzy logic to develop a rainwater
![]()
Figure 1. Fuzzy set versus conventional set.
harvesting system as an alternative source of potable water. A more innovative example has been presented by Ki and Ray [22] , using fuzzy logic to analyze infiltration trenches in order to determine the optimal locations for BMPs.
5. Factors Influencing Watershed Vulnerability
Watershed vulnerability, defined herein as high soil erosion and low infiltration potential, is a complex process to model. Several watershed factors such as slope gradient, length slope factor, soil texture, soil erodibility, soil permeability, landuse and landcover, vegetation cover, and drainage density influence the development of the watershed vulnerability model. The characteristics and effect of these factors on soil erosion and infiltration are outlined in this section.
Slope Gradient (S)
Slope gradient (S), usually expressed as percent slope or degree slope, is used to represent the degree of variability in elevation and is found to strongly influence the velocity and direction of stormwater runoff [1] . Generally, areas of watersheds with steeper slopes are characterized by high soil erosion susceptibility, whereas areas with mild or gentle slopes are suitable for infiltration and water harvesting practices.
Length Slope Factor (LS)
The length slope factor (LS) is an important topographic factor used in the Revised Universal Soil Loss Equation (RUSLE) computations to predict long term soil loss risks [23] . It is the ratio of soil loss from the field slope length to that from a 22.1 m (72.6 ft) length on the same soil type and gradient. Slope length is the distance from the origin of overland flow along its flow path to the location of either concentrated flow or deposition. High LS values represent watershed regions that are highly susceptible to soil erosion.
Soil Texture
Soil texture, classified with respect to percentages of sand and clay, is known to have a strong influence on the rate of infiltration and runoff potential of watersheds [24] . The National Resources Conservation Service (NRCS) have identified four distinct soil groups (A, B, C and D) based on their soil textures. Group A soils generally have higher infiltration potential, while Group D soils have lower infiltration potential.
Soil Erodibility (K)
The soil erodibility (K) factor is a quantitative soil property that plays a major role in determining the extent to which soil particles are eroded or transported by runoff water. Parameters such soil moisture, soil permeability, and organic matter content influence the degree of soil erodibility and stability. Soils with lower K values are highly stable compared to those with higher K values. The K factor is also an important parameter in the RUSLE soil erosion prediction calculations. Values for the K factor are available via the NRCS soil map STATSGO [25] .
Vegetation Cover
Normalized Difference Vegetation Index (NDVI) is a widely used standardization parameter to represent vegetation cover and land surface. NDVI values typically range from +1.0 to −1.0 with non-vegetated areas such as streams and rivers producing small or slightly negative values, and vegetated areas producing values ranging from 0.4 to 1.0 [26] . However, most processes for generating NDVI parameter values in GIS create a single band 8-bit raster which scales NDVI values of −1 to 1 to conform to a 0 to 255 range. The NDVI map used herein ranged from 115 to 160 based on New Mexico average maximum annual NDVI values from 1995-2009 [26] .
Soil Permeability (Ksat)
Soil permeability, or more correctly saturated hydraulic conductivity (Ksat), is an important physical soil property that refers to the ability of soils to conduct water. Ksat values give an indication of the degree of perviousness or imperviousness of soils. High Ksat values represent the portions of watersheds that are highly drained and are suitable for infiltration practices. Conversely, lower Ksat values correspond to watershed regions that experience high runoff and are prone to soil erosion.
Landuse and Landcover (LULC)
Landcover refers to the physical materials such as forest, crops, grass, and water covering the earth’s surface, while landuse refers to the use of land for human activities. LULC directly translates to the amount of water penetration into soil layers or being produced as runoff water. For example vegetative surfaces act as a protective layer which slows down runoff and allows adequate detention time for infiltration and groundwater recharge to occur [27] . A dense vegetation cover reduces soil erosion in an area even though the slope and rainfall may be relatively high; in contrast, fallow land results in increased soil erosion compared to other types of land covers [28] . Impervious surfaces such as concrete pavements, roads, roof tops, and frozen land, however, act as an impermeable layer that allows stormwater to flow directly to natural watercourses downstream.
Precipitation
Area precipitation characteristics (type, duration, intensity, and amount) are highly important in drainage projects due to their direct influence on the volume of generated runoff. Precipitation that falls on impervious surfaces eventually is carried to natural drainage. On porous surfaces, however, precipitated water may serve as a significant source of recharge for groundwater aquifer zones and rainwater harvesting sites. Precipitation characteristics play a significant role in soil erosion. Precipitation intensity and duration on soil erosion are important factors. In semi-arid watersheds, high intensity precipitation storms are believed to be main component for soil erosion. However, an event based erosion model applied to a catchment in southeast Spain indicated that large magnitude, low frequency events potentially contribute significantly to total soil erosion [29] . Regardless of the site characteristics, healthy watersheds are expected to be more resilient to the effects of precipitation changes [11] .
Drainage Density
Drainage density is a quantitative indicator that represents the ratio of the total length of streams to the total watershed area. Low drainage density values are favorable for infiltration while larger values indicate high runoff potential. Recent findings indicate that adequately drained watersheds have drainage densities less than 1 km/km2 while poorly drained watersheds have drainage densities greater than 5 km/km2 [27] . High drainage density values are an indication of extreme soil erosion vulnerability whereas lower values signify minor risk.
6. Modelling
The overall modelling approach used herein is based on using a MCDS FIS approach (Figure 1). The principal benefit of the method as it relates to developing watershed vulnerability is its ability to objectively and simultaneously consider an unlimited number of relative criteria. In this case, two sub-models are used to assess soil erosion potential and infiltration potential, respectively, each having multiple criteria.
Soil Erosion Sub-Model
Soil erosion is considered as one of the leading causes of environmental degradation in the world, contributing to long-term poor land productivity, surface and groundwater contamination, and loss of valuable natural resources. Nearly 99% of the total suspended solids (TSS) found in rivers and streams can be attributed to soil erosion in the United States [30] . Several models and methods have been suggested as a means of predicting soil erosion risks. Empirical methods like RUSLE continue to play a major role in soil erosion predictions and soil conversation programs. However, according to Black [18] , the best strategy for modelling environmental data and watershed attributes should be based on a decision support approach. Rainfall characteristics, vegetative cover, soil characteristics, slope gradient, and land-use are considered as main factors affecting soil erosion. Therefore, aspects of these factors, for example rainfall frequency and duration, and soil erodibility, may be used within a MCDS FIS approach to assess soil erosion vulnerability.
Infiltration Sub-Model
In order to fully address water availability concerns, it important to understand the impact of infiltration in water conservation programs. About 48% of irrigation water and 98% of domestic water in the United States are derived from groundwater [30] . Water planning and water conservation issues are predominantly at the forefront of regional priorities especially in arid regions which tend to experience limited water availability. Gautam et al. [6] proposed that the best strategy for addressing water availability concerns in the desert Southwest must be based on appropriate water harvesting and infiltration practices. Rokus [31] explained that well-intentioned infiltration policies must sustain groundwater recharge requirements and improve overall surface quality.
To examine infiltration potential within a MCDS FIS framework, attributes such as soil characteristics, vegetative cover, slope gradient, drainage density, and land-use are needed as input. For example, soil hydraulic conductivity significantly influences infiltration, as does vegetative cover with respect to surface entry of water. Slope gradient and drainage density affect surface runoff, which impacts overall infiltration. Thus, these interrelated factors may be used to assess infiltration potential using a MCDS FIS approach.
Soil Erosion Infiltration Potential (SIP) Model
Watershed protection programs should not only focus on soil erosion, but also infiltration as it is an integral component of the watershed surface and groundwater hydrology. Soil erosion potential is reduced with increased infiltration as overall surface runoff decreases. An innovative approach to watershed management is to develop a watershed vulnerability model that incorporates soil erosion and infiltration models. Mosadeghi et al. [32] explained that intersecting the areas of outcomes of multiple decision support models can be a useful technique in identifying spatial extents with greater confidence. In this study, a unique watershed model is presented on the basis of overlapping the areas of risks, or vulnerability, predicted by the soil erosion potential and the infiltration potential sub-models. An integrated watershed model could prove more beneficial especially in arid regions where water availability is a huge concern.
7. Site Description
This study was conducted using maps and digital data from Bernalillo County, New Mexico. The county has a total area of about 3023 km2 (1167 mi2) and is centrally located at longitude 106.669 and latitude 35.054. Having a population density of 220 persons/km2 (571 persons/mi2), Bernalillo County is the most populous county in the State of New Mexico, largely due to the City of Albuquerque. The region is covered by shrubs and scrubs (71%), deciduous, mixed, and evergreen forests (18%), associated urban (9%), and open water (1%). The two dominant soil orders are Entisols (47%) and Aridisols (32%). Bernalillo County is selected for this study because it falls within urban growth and high residential areas within New Mexico. Much of the region has an elevation of more than 1400 m. The 100 year 24 hour (100P24) storm event ranges from 96 to 164 mm (2.45 to 4.16 in) across the area. Limited water availability, unique soil characteristics, dominant shrubland vegetation, and intricate climate, coupled with a forecasted 35% population growth over the next 20 years, makes Bernalillo County a suitable case study area for the current research.
8. Materials and Methods
Tools and Data
A variety of software tools and databases are needed to facilitate implementation of an FIS approach to assess watershed vulnerability. The following sections summarize the tools and databases used for this research, as well as the extensive data preparation and processing required.
ArcGIS®: A Geographic Information System (GIS) provides an extensive toolset to manage multiple layers of spatially-distributed geographic data. Modelling approaches based on GIS have been shown to be successful in assessing regional scale soil erosion risks [15] and evaluating groundwater vulnerability in order to aid policy formulation [18] . ESRI’s ArcGIS® environment provides an infrastructure for making maps and working with large datasets of geographic and remotely sensed data. In this study, ArcMap® was used for a range of geo-processing tasks including referencing geographic data, creating erosion and infiltration models, analyzing data from decision support tools, and creating custom watershed vulnerability maps.
TauDEM: Terrain Analysis Using Digital Elevation Models (TauDEM) are a collection of Digital Elevation Model (DEM) processing tools [33] . TauDEM provides a variety of functions to extract hydrologic and topographic information from a DEM. Some common functions include developing hydrologically correct DEMs (pit removed), calculating flow paths and slopes, calculating slope contribution areas using the multiple flow direction method (D-Infinity), and delineating stream networks and sub-watersheds. TauDEM v 5.3 was used in this study to create the upslope contributing area map necessary for the generating a LS-Factor map.
Soil Data Viewer: The NRCS Soil Data Viewer is an ArcGIS® extension for analyzing soil data and building soil-based thematic maps. The program enables map makers to interpret soil attributes and properties without having to deal with the complexity associated with querying and processing soil databases. The Soil Data Viewer v 6.2 was used to create thematic layers such as soil texture, hydrologic soil groups, soil depth, and saturated hydraulic conductivity (Ksat) maps for the study area.
Thematic Layers: Several factors such as slope, soil erodibility, precipitation, LULC, soil type and texture, and vegetation were used in developing the watershed vulnerability model. The majority of the thematic layers were derived from base map layers such as LULC, DEM, soil data, and national hydrographic dataset (NHD) maps. The LULC map was provided by the National Land Cover Database [34] . The LULC layer was reclassified into seven (7) classes in order to fit the requirements of the model (water (11), urban (21), barren (31), forest (41), shrub (51), cultivated (81), and wetlands (91)). The DEM map of the study area was obtained from the National Map Viewer Platform [35] . Length slope (LS) and slope gradient (S) factors were extracted based on the DEM map. A drainage density factor (km/km2) was derived from the NHD map. The NRCS Soil Survey Geographic Database (SSURGO) soil map [36] , hydrographic datasets, and county boundary map were obtained from the New Mexico Resource Geographic Information System [37] . The soil texture and soil permeability factors were derived from the soil base map. Other input factors such as precipitation depth as 100P24, soil erodibility, NDVI, and rainfall erosivity maps were provided by Bulut et al. [26] , who developed a state-wide soil erosion risk factor map for New Mexico using the ArcGIS® environment.
Data Preparation and Processing: A grid size of 30 × 30 m pixel was chosen for all spatial layers to cover the study area extent. All acquired data with cell size different than 30 m pixel size were geometrically corrected using the cubic resampling technique in the Data Management toolbox of ArcGIS®. The Data management toolbox was then utilized to clip the input thematic layers to conform to the spatial extent of the study area. The Define Projection tool was used to project maps to the appropriate coordinate system in order for the maps to line up and display correctly in ArcGIS®. All layers were projected using the NAD 1983 Albers projection and the North American 1983 datum. In order to produce the soil texture and Ksat maps, the soil map together with the soil survey database (SSURGO) was loaded into the NRCS’s Soil Data Viewer program. Fuzzy membership functions, rule generation, and defuzzification of crisp inputs were conducted using the Fuzzy Rule Builder and FuzzyCell programs. In order for FuzzyCell to function properly, all input maps for the program were converted to the GRID raster format. All final map manipulations and overlay analysis were completed using the ARCGIS10 platform.
Derived Model Attributes
For the watershed vulnerability analysis herein, multiple objective criteria may be derived using the ArcGIS® environment.
Slope Gradient: A slope gradient map (degrees) was extracted from a digital elevation model (DEM) of the study area using the Slope tool in ArcGIS®. The resulting slope gradient map revealed that the eastern portions of the study area had high slope gradient values, while slope values varied from low to moderate towards the west.
D-Infinity Flow Accumulation: A flow accumulation (flowacc) map is needed to determine LS factor values for the RUSLE application. The flowacc map was created based on a three-step procedure within the TauDEM environment. The first step involved the development of a depressionless DEM raster using the Pit Remove tool. The purpose of creating a depressionless DEM raster was to ensure that areas enclosed by higher elevation values are filled in order to avoid erroneous flow-direction computations. Secondly, the D-Infinity Flow-Direction tool was then used to create a flow direction raster layer. The final flowacc raster was created using the depressionless DEM and the flow-direction maps by means of the D-Infinity Contributing Area tool in TauDEM.
Length Slope Factor (LS): The LS map for the study area was generated according to the overland flow method [5] . The LS model was developed in ArcGIS® using the Model Builder window. The model required a DEM in TIFF format and a D-Infinity flow accumulation map as model inputs. The output of the model was a 30 × 30 m 2D LS map of the study area.
Drainage Density: The Line Density function in Spatial Analyst toolbox of ArcGIS® was used to extract drainage density values from a hydrographic dataset. Drainage density values for the study area ranged from 0 to 25.54 km/km2.
Support Practice Factor (P): The support practice factor (P) value thresholds were determined based on the slope gradient (S) map expressed as percentages as proposed by Ahamed et al. [17] . The P factor is a quantitative representation of the soil loss ratio, and it is mostly used as a discretionary parameter in the RUSLE to determine long term average soil loss. P values for the study area ranged from 0.5 to 1 depending on their relationship with the percent slope gradient map. For example, for slopes between 5% and 8%, P is set at 0.5; for slopes greater than 25%, P equals 1 [17] .
Development of FIS Model
Fuzzy logic is a decision support tool that contributes to a better representation and modeling of GIS-based data [38] . Two common types of fuzzy inference systems (Mamdani and Sugeno) are identified in literature [39] . Regardless of the type inference, the basic of structure of any FIS involves developing membership classes, aggregating fuzzy rules, and the defuzzification of crisp input parameters.
FuzzyCell: FuzzyCell is among the most widely used fuzzy logic tools for processing geographic data. It allows users to incorporate human knowledge spatial analyses, classify data into membership functions, generate “if-then” rule aggregation, and provide inference methods for defuzzification. FuzzyCell was used as the main fuzzy inference system for classifying GIS raster data in fuzzy memberships and defuzzifying inputs. All generated fuzzy output maps are based on the Mamdani inference system and the Center of Area (COA) defuzzification method.
Fuzzy Rule Builder (FRB: Fuzzy Rule Builder (FRB) is a custom fuzzy inference tool developed in Microsoft Visual Studio 2010® express in order to handle fuzzy standardization and generating “if-then” fuzzy rule aggregation. The fundamental concept behind fuzzy inference processing involves the construction of a collection of “if-then” rules. The role of the “if then” rules are to establish a framework for converting crisp input data into defuzzified outputs. The total number of rules generated was determined as Fn, where n is the total number of input factors and F is the number of membership functions. A total of 243 fuzzy rules were generated for each model. Rules from the program were subsequently exported to FuzzyCell environment for fuzzy inference and defuzzification.
Fuzzy Membership Classifications for Watershed Vulnerability
Each input variable for fuzzy logic processing was coded into 3 membership classes (low, moderate, and high). The relationships between each input parameter and the output risk maps were used to establish membership classification threshold values. The trapezoidal shaped membership function (with a partial overlap to account for uncertainty) was used as the basis of establishing membership classes. Trapezoidal membership function curves were chosen due to its unique potential to optimize the performance of fuzzy inference processing [15] . The shape of membership function curves were established based on published threshold values and expert knowledge. All input variables were subsequently processed and defuzzified for decision making using the Mamdani fuzzy inference engine available in FuzzyCell. The quantitative description of the membership functions are given in Table 1 & Table 2. The output risk memberships were categorized based on the five suitability classes (Figure 2).
Fuzzy Standardization and Weighting
To be able to utilize input layers of variable spatial representation in fuzzy analysis, it was important that layers be associated with a common scale of measurement. This process of transforming input layers to a common scale (0 - 1) is termed standardization. All fuzzy input variables were transformed into a standardized input layer either using Equation (1) or Equation (2) for positive and negative correlations, respectively, with soil erosion or infiltration thematic layer inputs.
(1)
![]()
Table 1. Fuzzy membership classes for soil erosion models.
![]()
Table 2. Fuzzy membership classes for infiltration models.
![]()
Figure 2. Suitability classes for fuzzy logic models.
(2)
where wi is a standardized weight, xi is an initial value of criteria i, xi max is the maximum value, and xi min is the minimum value [40] .
9. Vulnerability Assessment
The development of the watershed vulnerability model presented in this study involved 3 major steps: 1) preparation of thematic layers, 2) extraction of fuzzy memberships using Fuzzy Cell, and 3) overlaying FIS results to produce the watershed output map. The watershed model depended on the determination of soil erosion and infiltration potential and identifying factors involved. Two sub-models were created including a soil erosion model and an infiltration model.
Vulnerability Map: The final watershed vulnerability map is a 30 × 30 m cell raster produced by overlaying soil erosion and infiltration potential maps obtained from the FIS process. The overlay process was conducted in ArcGIS® using the Weighted Overlay Tool. Equal weight was assigned to each map in the overlay analysis. A detailed illustration of the vulnerability assessment procedure is given in Figure 3. The various GIS Model Builder tools implemented herein are provided in Amankwatia [9] .
Model Comparison with RUSLE: Empirical soil erosion risk prediction models such as RUSLE continue to play a vital role in soil conservation and protection programs. RUSLE is a universally accepted empirical method developed by United States Agricultural Research Service [23] for assessing long-term soil loss. In this study, the resulting map of watershed vulnerability model was contrasted against annual soil erosion risk predicted using the RUSLE. The RUSLE map was developed in ArcGIS® based on annual soil loss as:
(3)
where E (tons/ac/yr = 2.242 tonnes/ha/yr) is the predicted average annual soil loss, R (hundreds of ft.tonf.in/ac/hr/yr) is the rainfall erosivity factor, K (ton/ac/unit R) is the soil erodibility factor, LS is the length slope factor, C is the cover management factor, and P is the support practice factor [5] .
The soil erosion map by merging five input raster layers, as defined in Equation (3). In order to compare the results to the watershed vulnerability model, the RUSLE output was subsequently divided into five vulnerability groups in tons/ac/yr, where 1 ton/ac/yr equals 2.242 tonnes/ha/yr: 1) not vulnerable (0 - 2); 2) slightly vulnerable (2 - 6); 3) moderately vulnerable (6 - 50); 4) highly vulnerable (50 - 150); and 5) extremely vulnerable (>150). These groups span across what is considered sustainable soil loss at 2 - 5 ton/ha/yr [41] to rates associated with extreme gully erosion or disturbed bare soil areas, such as construction sites, and are patterned after groupings of Rahman et al. [28] . A statistical comparison of the RUSLE vulnerability map with the FIS vulnerability map is now possible on a pixel by pixel basis using the Band Collection Statistics tool in ArcGIS® Spatial Analyst extension. The result is a coefficient of determination
![]()
Figure 3. Risk assessment process for the BMP-SIP model.
(r2) between the two maps.
10. Results and Discussion
Evaluation of Watershed Vulnerability
The assessment of the watershed vulnerability for this study was performed in ArcGIS® through a weighted overlay analysis of infiltration and soil erosion potential maps. The soil erosion maps were developed as a result of merging 5 input raster data: slope gradient, precipitation, K factor, NDVI, and LULC. Similarly, the infiltration map was developed by merging drainage density, soil permeability, soil texture, slope gradient, and LULC raster layers. The final output of each approach and subsequent overlay, which comprised an overall vulnerability map, is given in Figure 4.
As displayed in Figure 4, the study area was divided into five watershed vulnerability zones, namely, not vulnerable (N), slightly vulnerable (SV), moderately vulnerable (MV), highly vulnerable (HV), and extremely vulnerable (EV). Regions in red indicate extreme vulnerability, whilst regions in green represent no vulnerability. High to extreme vulnerability zones were mainly distributed in the eastern sections of the study area. The analysis indicated that about 76% of the study area is susceptible to moderate vulnerability (MV), and less than 12% of the area is subject to experience high or extreme vulnerability (HV/EV). The principal assumption in this study is that areas with low infiltration potential and high soil erosion potential are the most vulnerable and, therefore, require closer attention relative to other areas within the watershed. These results are summarized in Table 3.
Evaluation of Soil Erosion and Infiltration Risks
In this study, the distributions of soil erosion and infiltration risks were also considered. In order to compare the model outcomes, the five vulnerability classes (N, SV, MV, HV, and EV) used for the watershed vulnerability model were applied to both the soil erosion and infiltration risk models. The risk distribution of soil erosion and infiltration models followed a somewhat different pattern relative to the combined watershed vulnerability model, with soil erosion showing a higher percentage in the slightly vulnerable category and infiltration having a higher percentage in the moderately vulnerable as did the combined overlay.
![]()
Figure 4. Fuzzy logic watershed vulnerability risk map.
![]()
Table 3. Percentage of vulnerability distributions for the watershed model.
N = Not Vulnerable, SV = Slightly Vulnerable, MV = Moderately Vulnerable, HV = Highly Vulnerable, and EV = Extremely Vulnerable.
Figure 5 shows the FIS soil erosion vulnerability. As expected, the most problematic soil erosion prone areas were located near the eastern regions of the Bernalillo County. This area has a predominance of forested land cover; however, it is also characterized by high slope gradient coupled with a highly erodible soil and high 100P24 precipitation as depicted in these respective thematic layers developed for this study [9] .
Infiltration risk followed a similar pattern; however, the infiltration risk distribution also revealed high risk areas around the central and western portions of the study area. The vulnerability distributions for soil erosion and infiltration risk are summarized in Table 3.
Comparison with the RUSLE Model
The areas of watershed vulnerability for the study area were contrasted against soil erosion prediction from the RUSLE model. The spatial distribution of soil erosion vulnerability for Bernalillo County is given in Figure 6. For the two lower vulnerability groups (N and SV), the RUSLE predictions were considerably different with the results predicted by the soil erosion FIS vulnerability model, especially the not vulnerable class (Table 3). For both models, the percentage of the study area predicted to be susceptible to high and extreme vulnerability was small. For the slight vulnerability and moderate vulnerability groupings, the watershed vulnerability model predicted higher and lower percentage, respectively. RUSLE predictions were much higher for the not vulnerable category compared to the soil erosion vulnerability estimates. Soil erodibility was high within parts of the study area based on the magnitude of the derived K factor map. A not vulnerable vulnerability group from 0 to 1 ton/ac/yr versus 0 to 2 ton/ac/yr as modeled would better align the predictions of vulnerability between the FIS and RUSLE models for the two lower vulnerability classifications. For thin erodible soils, a tolerance level for sustainability of 1 ton/ac/yr soil loss is recommended [23] , and is taken herein as the lower bound of vulnerability.
![]()
Figure 5. Fuzzy logic soil erosion risk map.
![]()
Figure 6. RUSLE soil erosion (tons/ac/year) prediction.
![]()
Table 4. Risk distribution (%) for LULC classes.
N = Not Vulnerable, SV = Slightly Vulnerable, MV = Moderately Vulnerable, HV = Highly Vulnerable, and EV = Extremely Vulnerable.
Influence of Landuse and Landcover (LULC)
An attempt was made to quantify the influence of landcover and landuse activities on the watershed risk assessment using the 7 LULC classes given in Table 4. LULC vulnerability was evaluated by intersecting the watershed vulnerability FIS map with LULC classes using the Combine tool in ArcGIS®.
The results suggest that shrub, cultivated land, and forest zones are the most affected zones. About 65% and 9% of these LULC areas fell in the moderately vulnerable and highly vulnerable category, respectively, with slightly over 1% classed in the extremely vulnerable category. In general shrub lands were found to be the most vulnerable LULC category (highest total percentage in the top three vulnerability classes). This aligns with expectations based on area climate and literature. Ground cover, or vegetative litter, often may be sparse in semi-arid shrub lands. The role of this layer in modulating surface runoff and soil erosion remains poorly understood; however, in general, runoff volume and sediment yield are reduced in the presence of litter cover [42] . Langbein & Schumm [43] reported that mean annual sediment yield correlates with annual precipitation and is highest for semi-arid shrub land.
Use of Average Annual Maximum NDVI
NDVI input was a criterion for the soil erosion sub-model only. Hydrological partitioning at the watershed scale is a function of vegetative cover. Seasonal precipitation in New Mexico is distinctly bimodal with less intense spring rains and late summer monsoonal rains. As such a bimodal NDVI seasonal cycle is observed with a late spring-early summer initial peak, a mid-late summer pre-monsoon dip, and a second maximum peak in early autumn [44] . For the watershed vulnerability study herein use of an average pre-monsoonal NDVI would be more appropriate to simulate the effect of vegetative cover on infiltration potential and soil erosion potential in lieu of an average maximum NDVI. Such data was not available at the time of study.
11. Advantages of the MCDS Approach
Multi-Criteria Decision Support approaches such as FIS provide a relatively easy way of determining watershed vulnerability distributions on a regional scale. From the viewpoint of decision makers, these methods can provide several strategic benefits including the ability to incorporate project objectives and a variety of quantitative data. One possible explanation for its gain in popularity is the lack of statistically meaningful data and methodology to deal with issues associated with uncertainty. Overlapping memberships addresses some uncertainty in the FIS approach.
As demonstrated in this study, FIS can be an effective means of generating regional scale vulnerability map for watershed investigations. The FIS approach integrated both soil erosion potential and infiltration potential using specific criteria to characterize watershed vulnerability, whereas the RUSLE model only accounts for soil erosion. The FIS approach offers several benefits in developing regional scale watershed vulnerability model relative to conventional field surveys of watershed vulnerability, in terms of time and effort and degree of complexity. It is important to acknowledge, however, that detailed field investigations might be required in order to extend the applicability of the methodology discussed herein.
12. Conclusion and Recommendation for Future Work
Quantifying a regional scale risk index map can provide spatially relevant information for engineering management and decision making. In view of this, a 30 × 30 m watershed vulnerability risk map was generated for the Bernalillo County, New Mexico using an FIS methodology. The approach incorporated drainage density, slope, precipitation, vegetative cover, landuse and landcover, soil erodibility, soil texture, and soil permeability as the primary factors contributing to watershed vulnerability. The final map was produced by merging soil erosion and infiltration potential using the Weighted Overlay method in ArcGIS®.
Analysis of the vulnerability map identified that the most problematic areas were present in the eastern portions of the study area. Based on the analysis, less than 12% of the study area was identified to fall within the high vulnerability (HV) or extreme vulnerability (EV) category, whereas less than 1% of the total area was associated with no vulnerability (N). Around 88% of the area was predicted to experience slight (SV) to moderate (MV) vulnerability. About 9% of shrubs, cultivated lands, and forest zones were predicted to have high vulnerability (HV) with shrub land being most vulnerable class overall.
In terms of spatial extent and general patterns, the FIS soil erosion vulnerability output map visually compared with the RUSLE prediction. The AHP approach showed a higher correlation with RUSLE than the FIS method [9] . However, the AHP and FIS models revealed a strong correlation with each other (r2 = 0.99). Data analysis on a pixel by pixel basis using Band Collection Statistics tool in Spatial Analyst indicated that the correlation coefficient (r2) between FIS watershed vulnerability and RUSLE vulnerability predictions was 0.66. However, further analysis and field investigations may be necessary in order to validate the usefulness of the FIS model. For example, watershed areas assessed as being potentially vulnerable to soil erosion could be validated through ground-truthing of dominant erosion features such rill erosion, gully erosion, and cut-bank erosion. Gopal et al. [45] used two experts in remote sensing and aerial photography to assess accuracy in their fuzzy inference classification of urban landscapes based on a method by Zadeh [8] .
The current lack of high resolution spatial data and current site-specific landuse and land cover (LULC) information may limit the potential to field-test the watershed vulnerability model. It is, therefore, important that future extensions of the model focus on generating high resolution (e.g. a 10 × 10 m raster data) especially for the highly vulnerable regions and also obtaining detailed site-specific LULC information. NDVI data appropriate to the critical pre-monsoonal period should be developed to better address soil erosion potential. In addition, photographic documentation from various field sites is needed to understand how the landscape looks for different classes of soil erosion risk and infiltration potential. In order to improve on the functionality and applicability of the model, other factors such public awareness and economic indicators could be incorporated into the model.
Information from the FIS watershed vulnerability model can be useful for managing new structural or non-structural SCM projects or improving existing ones to mitigate impact to stormwater quality and quantity for identified vulnerable areas. A recent study showed that runoff water originating from nearby urbanized watersheds in Albuquerque, NM carried high levels of polychlorinated biphenyls (PCBs) into the Rio Grande [46] . These findings have inspired the Bernalillo County Public Works Division (BCPWD) to adapt a proactive approach toward identifying potential source areas of PCBs in their stormwater discharges. As noted, the most vulnerable area in Bernalillo County is the eastern portion with its high slope and higher precipitation. This large area may benefit from strategically placed structural SCMs, such as dry storage basins and filter strips, to manage and reduce runoff and sediment transport offsite. Having a MCDS FIS model that incorporates a multiplicity of technical, as well as non-technical, criteria can be beneficial in facilitating development of sustainable drainage programs. Use of this approach to identify potentially vulnerable watershed sites, pending validation by field investigations, is illustrated through its application to this study area. However, a sensitivity analysis of model criteria impact on watershed vulnerability should be conducted to evaluate the robustness of the levels of watershed vulnerability indicated.