Use of Multi-Criteria Decision Analysis (MCDA) for Mapping Erosion Potential in Gulf of Mexico Watersheds

: The evaluation of soil erosion is often assessed using traditional soil-loss models such as the Revised Universal Soil-Loss Equation (RUSLE) and the Soil and Water Assessment Tool (SWAT). These models provide quantitative outputs for sediment yield and are often integrated with geographic information systems (GIS). The work described here is focused on transitioning towards a qualitative assessment of erosion potential using Multi-Criteria Decision Analysis (MCDA), for improved decision-support and watershed-management prioritization in a northern Gulf of Mexico coastal watershed. The foundation of this work conceptually deﬁned watershed erosion potential based on terrain slope, geomorphology, land cover, and soil erodibility (as deﬁned by the soil K-factor) with precipitation as a driver. These criteria were evaluated using a weighted linear combination (WLC) model to map generalized erosion potential. The sensitivity of individual criteria was accessed with the one-at-a-time (OAT) method, which simply removed one criterion and re-evaluated erosion potential. The soil erodibility and slope were found to have the most inﬂuence on erosion-potential modeling. Expert input was added through MCDA using the Analytical Hierarchy Process (AHP). The AHP allows for experts to rank criteria, providing a quantitative metric (weight) for the qualitative data. The individual AHP weights were altered in one-percent increments to help identify areas of alignment or commonality in erosion potential across the drainage basin. These areas were used to identify outliers and to develop an analysis mask for watershed management area prioritization. A comparison of the WLC, AHP, ensembled model (average of WLC and AHP models), and SWAT output data resulted in visual geographic alignment between the WLC and AHP erosion-potential output with the SWAT sediment-yield output. These observations yielded similar results between the qualitative and quantitative erosion-potential assessment approaches, with alignment in the upper and lower ranks of the mapped erosion potentials and sediment yields. The MCDA, using the AHP and ensembled modeling for mapping watershed potential, provided the advantage of more quickly mapping erosion potential in coastal watersheds for improved management of the environmental resources linked to erosion.


Introduction
Sediment is the largest volumetric nonpoint-source pollutant to surface waters [1][2][3] and one of the most important water-quality problems in the United States [4][5][6]. Upland watershed erosion is a serious issue for estuaries and the coastal region of the southeastern United States. During precipitation events, overland and streambank erosion increases in the watershed, often resulting in degradation to downstream resources in the associated estuary. Erosive rates are amplified in areas experiencing active land use changes with agriculture, and increasing urbanization and industrialization [7,8]. The influence of the growing human population and unrestricted development in coastal watersheds is proving to be very detrimental to the overall integrity of the fragile, yet highly productive estuarine vegetation index), land use, soil texture and slope aspect were studied for soil-erosion risk assessment [25]; additionally, drainage density, slope, land use/cover, and runoff measurement were used for identifying potential zones for rainwater harvesting [30], etc. Therefore, based on the availability of data in the specified zone, the number of factors may be varied as an input to the models.
Multi-criteria decision analysis (MCDA) methods have become very popular for spatial planning and management issues and are a significant tool for decision makers, especially for multicriteria assessment [31]. The applications of MCDA are wide. It is applied to identify priority areas for soil-erosion risk measurement [25], to calculate landscape deformation index [31], to identify potential zones for rainwater harvesting [30], in the field of transport for determining suitable management [32], to generate the ranking of green bonds in corporate office management [33], etc.
Hence, expanding GIS utilization for MCDA has improved decision-support models for land-based suitability evaluations. These expanding efforts have increased the need for ways to evaluate the performance of the models and tools utilized, as well as the sensitivity of the variables or layers used [34,35]. There are numerous procedures that are used with GIS for MCDA; examples include Boolean overlay, weighted linear combination (WLC), ordered weighted averaging (OWA), and the analytical hierarchy process (AHP) [36]. The WLC is one of the most commonly used decision-support tools in the GIS environment [37,38]. Additionally, GIS coupled with the AHP [39] is proving to be an important tool for MCDA [40,41]. GIS utilizing AHP is an established and credited approach to MCDA for land-resource-management decisions [42,43], and is an important part of sustainable land-planning approaches [44,45]. The AHP has been found to be a robust method for determining criteria weights based on expert input [46] and works well with MCDA in the GIS environment. Additionally, the combination of GIS and AHP is useful for MCDA in the management of natural resources related to soil-erosion mapping [47]. While models and tools of this type do not allow for the quantification of sediment yields or soil-loss rates due to erosion, they do offer resource managers and decision makers the necessary information to better manage and prioritize watersheds and the related resources.
Therefore, qualitative modeling approaches are often driven by MCDA, typically with expert input. This makes them very useful in the decision-making process, specifically for tasks such as vulnerability assessments and other methods [48]. The qualitative nature of MCDA often requires nontraditional methods of uncertainty assessment. Sensitivity analysis is one of the common methods used to reduce uncertainty in the variable weights, which can assist with identifying stability in model performance with changing criteria weights [34]. Sensitivity analysis with GIS-based MCDA can also provide insights into the spatial aspects of the changing criteria weights. Feick and Hall suggested that efforts to analyze criteria weight sensitivity can help to geographically visualize the sensitivity of the results [49]. Another approach to reducing uncertainty in the modeling approach is the combination of different models using the average, which is termed as ensembled model [50,51]. The ensembled model is a useful combinational approach and can enhance the strength of prediction mapping while reducing the weakness of each source map [51]. Studies showed better/improved prediction in clay content mapping for soil-quality assessment and decision making in land use [50], and for combining digital soil-property maps derived from disaggregated legacy soil-class maps and scorpan-kriging (using soil pan data) [51] from the ensembled model, respectively.
The primary aim of this research is to develop a qualitative assessment to map erosion potential using WLC, AHP, and ensembled modeling approaches for watersheds associated with the northern Gulf of Mexico. This assessment will provide a process to resource managers for the identification and prioritization of watershed management areas. The mapped erosion potential found from these models will be compared to the sediment yield from the SWAT model (https://swat.tamu.edu/). The SWAT model is widely used across the globe in assessing soil-erosion prevention control, nonpoint-source pollution control, and regional management in watershed (https://swatplus.gitbook.io/docs/). The SWAT was used in the development of the Weeks Bay Watershed Management plan. The model delineated the Weeks Bay watershed into 237 sub-watersheds (197 for the Fish River and 40 for the Magnolia River); these are used to produce the computational hydrologic response units (HRU's) in SWAT. Sediment-yield results from this model were based on 2011 land use/cover, and it was reported that over half of the sediment yield was produced from about one-third of Weeks Bay watershed [52]. Therefore, the comparisons will be limited to basic observations between the qualitative and quantitative output of the data. The comparisons will show a general visual alignment in sub-basins of increasing development and headland drainage areas in the study area. The comparison against the output of the SWAT model will help the resource managers to look at scenarios or management priorities without the understanding and execution of more complex soil-loss models. The ensembled model will aid in the visualization of the priorities of the management areas. The models will serve as a base for the multi-criteria decision analysis (MCDA) of erosion potential by decision makers and resource managers.

Study Area
The Weeks Bay watershed is located on the eastern shore of Mobile Bay in Alabama. It is an ideal basin for the assessment of erosion potential, as it relatively small and secluded from surrounding watersheds. The watershed is limited to inputs from two major rivers (Fish and Magnolia) that both directly drain to Weeks Bay ( Figure 1). The Weeks Bay watershed is a diverse natural and anthropogenically influenced landscape with natural, forested, agricultural, and developed areas that are reflective of the region's natural resources and demographics [52]. The area is within the humid subtropical climate region, characterized by warm summers and relatively mild winters. Average annual precipitation averages about 165 cm due to winter storms (cold fronts), summer thunderstorms (including those from the sea breeze), and tropical systems. The abundant water resources in the area make for a range of very productive land uses from timber production, cash-grain crops, and forage production [53]. The Fish River provides nearly 75% of the total discharge to the bay itself and is made up of three sub-watersheds (Upper, Middle, and Lower Fish River). The Magnolia River provides the remaining discharge and consist of a single sub-watershed.

Data Collection and Processing
The process used to map soil-erosion potential was based on several geospatial variables. National data sources were used for these variables to ensure transition between different watersheds and scalability. Below, we describe the most relevant variables to be included in the model.

Slope
The slope for the study area was calculated using the USGS (United States Geological Survey) 30-m National Elevation Dataset (https://www.usgs.gov/3d-elevation-program; accessed on 9 October 2019). The slope calculation for each raster cell was based on the amount of descent between it and the surrounding eight cell neighborhoods using Horn's algorithm [54]. The maximum value of descent was, thus, recorded as the cell's slope, and could be calculated in percent or degrees. The slope raster was then normalized to 0-1.

Soil Erodibility
The soil erodibility (K-factor) was from the 30-m gridded USDA (United State Department of Agriculture) Soil Survey Geographic Database (https://www.nrcs.usda.gov/wps/ portal/nrcs/main/soils/survey/; accessed on 10 March 2020). K-factor accounts for both the susceptibility of a cell to soil erosion based on soil texture and rate of runoff. Values less than 0.2 are considered as low erodibility, 0.

Data Collection and Processing
The process used to map soil-erosion potential was based on several geospatial variables. National data sources were used for these variables to ensure transition between different watersheds and scalability. Below, we describe the most relevant variables to be included in the model.

Slope
The slope for the study area was calculated using the USGS (United States Geological Survey) 30-m National Elevation Dataset (https://www.usgs.gov/3d-elevation-program; accessed on 9 October 2019). The slope calculation for each raster cell was based on the amount of descent between it and the surrounding eight cell neighborhoods using Horn's algorithm [54]. The maximum value of descent was, thus, recorded as the cell's slope, and could be calculated in percent or degrees. The slope raster was then normalized to 0-1.

Stream Density
The stream density was a 30-m raster collected from the USGS National Hydrography Dataset (https://www.usgs.gov/national-hydrography/national-hydrography-dataset; accessed on 9 October 2019). Higher instances of stream density are associated with increased erosion rates, specifically as they relate to the dissection of the landscape and landdrainage system interactions [55]. Similar data layers are used in soil-erosion analysis [56] and are also used in numerous landscape evolution models that simulate erosion and deposition [57]. The density function used for calculation utilized a neighborhood area with a specified search radius; all stream segments intersecting the area were counted and a continuous surface with the specified cell size was returned. The default search radius used in commercial GIS software is based on the minimal spatial dimension of the dataset [58]. The derived density surface was normalized by the maximum value within the Weeks Bay Watershed. The resulting data layer was a continuous index of stream density with unitless values ranging from 0.055-1.

Soil Brightness
The soil brightness was calculated dynamically from the tasseled cap transformation of the 30-m Land cover raster from USGS Global Land Survey Dataset (https://www.usgs. gov/landsat-missions/global-land-survey-gls; accessed on 9 October 2019). Hence, the soil brightness band of the Tasseled Cap transformation provided an index of measure for soil reflectance/exposure and not just the lack of vegetation [59]. The soil brightness data from the GLS dataset were extracted, subset to the Weeks Bay Watershed, and normalized by the maximum value. The resulting data layer was a continuous index of soil brightness with unitless values ranging from 0.018-1.

Precipitation
The precipitation data were collected from the 4-km 30-year normal grid of PRISM dataset (https://prism.oregonstate.edu/normals/; accessed on 15 March 2020) for rainfall variation across the basin. These data were extracted from the database for the region of interest, normalized by the maximum value, and resampled to 30 m (from 4 km). The resulting data layer was a continuous index of annual precipitation climatology with unitless values ranging from 0.96-1.

Model Description
The concept of watershed erosion is summarized as the total erosion for the combination of physical erodibility, land sensitivity, and precipitation erosivity factors [20,23,41,56]. For this study area, the physical erodibility was measured by slope and stream density; land sensitivity included measures of soil K-factor and soil brightness (exposure); and rainfall erosivity was measured as average precipitation of the watershed (rainfall variation). A schematic diagram of the modeling approach for watershed erosion potential is indicated in Figure 2. The algorithm used for this WLC mapping was a standard weighted linear combination (WLC) for the summation of the five raster data layers, i.e., slope, stream density, soil brightness, soil erodibility, and precipitation [37]. The weights for the AHP model were calculated using Saaty's method of a continuous rating scale for pairwise comparison [46]. This procedure set each data layer with a possible data range of zero to one for a common scale of assessment. Zero would be a minimal impact on erosion potential, with values of one having the greatest impact. The scale for the AHP model is shown in Table 1. Each data layer was then compared individually with the other data layers as they relate to erosion potential. Weights for each data layer were assigned based on results from the pairwise comparison matrix (Table 1). With all layers standardized and weighted, the WLC was used to apply the weights from the expert input for the assessment of erosion potential. This allowed each data layer to be multiplied by the expert-defined weight and then summed for a continuous surface of overall erosion potential. The algorithm for the WLC and AHP model is shown in the Equations (1) and (2), respectively.
where EP is the watershed erosion potential, P is the average Precipitation, K is the soil erodibility, S is the slope length and steepness, SB is the soil brightness, and SD is the stream density of this region. a is the standard weight for the WLC model, and a 1 , a 2 , a 3 , a 4 , and a 5 are weights expected from running the AHP model.
where is the watershed erosion potential, is the average Precipitation, is the soil erodibility, is the slope length and steepness, is the soil brightness, and is the stream density of this region. is the standard weight for the WLC model, and , , , , and are weights expected from running the AHP model.    The erosion potential (EP EN ) from the ensembled model was calculated and is mentioned in Equation (3).
where EP EN is the ensembled erosion-potential model, n is the number of models, and EP i indicates the erosion potential from each of the i th models.

Sensitivity Analysis
A basic sensitivity assessment was performed for the WLC model variables. The assessment was a simplistic one-at-a-time (OAT) procedure wherein a single variable or layer was removed and the erosion-potential analysis was processed again. The procedure followed the method used by Chen et al. [34] and Romano et al. [36].

Ranking of the Management Area
The study area was divided into 18 management areas based on the smaller streams in the watershed. The management areas were ranked or prioritized based on the average EP value for each zone. The average EP values found in the WLC, AHP, and SWAT models were compared to visualize the similarities in the zones' priorities for management, to support improved decision making for management and prioritization of resources.

WLC Model
Potential watershed erosion cells were mapped with a standard weighted linear combination (WLC) to build a foundation for the qualitative assessment for watershed erosion. Criteria sensitivity was evaluated using the one-at-a-time (OAT) method to better understand the influence of the landscape layers. The WLC model was executed using the five raster data layers, weighted equally, for the initial assessment of generalized erosion potential. The output of the WLC model was a continuous surface of erosion potential based on physical erodibility, land sensitivity, and 30-year precipitation for the Weeks Bay watershed. Mean erosion potential was calculated for the entire watershed (0.527, S.D. = 0.057) and the four sub-basins within the Weeks Bay watershed. Field observations during site visits showed that the upland, head water areas of the watershed, and the areas dominated by cultivated agricultural areas are expected to have higher erosion potential. Lower erosion-potential values are expected in the densely vegetated riparian and marsh areas.
The erosion-potential data were classified based on the standard deviation spread from the mean erosion potential across the drainage basin, to better define areas based on the upper and lower ranks of erosion potential. This resulted in seven classes that were used to define the ranks for the cells, with class 1 having the lowest potential and class 7 having the highest potential. At the watershed level, 1.5% of the cells are in the upper-most erosion-potential ranks (classes 6 and 7), approximately 14% are in the moderate erosionpotential rank (class 5), and the lower erosion-potential ranks (classes 1-3) are similar in distribution to the upper ranks (classes 5-7).
Erosion potential across the four sub-basins level are varied in comparison. The two downstream basins (Magnolia and Lower Fish) have decreased proportions of cells around the mean erosion potential, and the two upstream basins (Upper and Middle Fish) have increased cell proportions around the mean erosion potential. The Magnolia River subbasin has the largest count of cells in the upper erosion-potential ranks (classes 5-7), at 23%. The Lower Fish River sub-basin has the next highest count, with almost 17% of the sub-basin in the upper erosion-potential ranks, and the Upper and Middle Fish sub-basins have 12% or less in the upper ranks. Table 2 provides a complete description of cell counts (with upper and lower ranks) at the basin and sub-basin levels.

WLC Model Criteria Sensitivity Assessment
The sensitivity of individual criteria was accessed using the one-at-a-time (OAT) method by removing one criterion from the WLC model and recalculating erosion potential. This resulted in five additional WLC model outputs of erosion potential produced by equally weighting four of the five criteria; the detailed results are mentioned in Table 3. The five additional OAT WLC models were compared to the initial WLC model to understand the influence of each variable in the WLC. The OAT WLC models all had moderate-tostrong correlations with the initial WLC model. The model without the precipitation variable had the strongest correlation (R = 1.00), followed by the model run without slope input (R = 0.94). The runs without stream density and soil brightness were moderately correlated, R = 0.88. The run with the weakest correlation was the one without K-factor, R = 0.79. The correlation results showed that the WLC model was most sensitive to the K-factor variable, moderately sensitive to the variables of stream density and soil brightness, and least sensitive to the precipitation and slope variables.

AHP Model
This Analytical Hierarchy Process (AHP) model was executed on the five raster data layers utilized in the WLC model. The AHP experts input defined weights for each of the five layers with the pairwise comparison. The weight assignments for each of layers were based on scores from the experts' qualitative ranking of factors. The most weight was given to terrain slope at 33.8% with lesser weights given to geomorphology, land cover, and soil erodibility (~15.0% for each). The 30-year precipitation was effectively left unchanged at 20.5%. The AHP mean erosion potential for the entire watershed (0.472, S.D. = 0.051) decreased as compared to the WLC erosion potential. The range of the AHP erosion potential increased slightly from that of the WLC and was strongly correlated, with a Pearson's R value of 0.923.
The AHP erosion-potential data were classified into seven classes based on a standard deviation spread (Figure 3b) of just the WLC erosion-potential data (Figure 3a). At the watershed level, 2.7% of the cells are in the upper-most erosion-potential ranks (classes 6 and 7) and approximately 12% are in the moderate erosion-potential rank (class 5). The lower erosion-potential ranks (classes 1 and 2) are similar to the upper ranks with 2.27% of the cells. The AHP model produces slightly more cells in the lower ranks than the upper ranks of erosion potential ( Table 4). The spatial differences between the AHP model and the WLC model erosion-potential class changes are seen in Figure 3c, highlighting the areas where the AHP increased or decreased erosion potential from the WLC model.     The AHP weights were altered in one-percent increments to identify outliers with AHP model runs. An example of weight alteration in the AHP model for 41 iterations for a single layer is displayed in Appendix A, Table A1. The percentage of outliers at each map (or grid) cell was calculated to look at variations spatially. A threshold of 25% was used to define areas with minimal alignment between AHP runs. About 37.5% of the mapped watershed cells were producing inconsistent or varying erosion-potential results. About 62.5% of the mapped watershed cells were in alignment with 75% of the model runs, with no outliers. This area of the watershed is where the AHP model runs were in alignment for the upper and lower ranks of erosion potential. Erosion cell counts in this area for the moderate and upper ranks (classes 5, 6, and 7) decreased proportionally. The proportional decrease shows that the higher instances of outliers are not clumped within the lower or upper ranks of erosion potential. This provides the definition of a focus area (an analysis mask) for increased erosion potential in the Weeks Bay watershed (Figure 4). This analysis mask serves as spatial filter for mapped erosion-potential cells for improved agreement in the modeled outputs.

WLC and AHP Model Comparison
Comparisons between the WLC and AHP models were expanded to include an evaluation of the outputs with a SWAT model for the Weeks Bay Watershed. This comparison went beyond the mapped potential erosion cells to evaluate the outputs for watershed management scenarios. The comparison took the classified erosion cells for the WLC and AHP models and mapped them with the sediment-yield data for 237 sub-watersheds from the SWAT model output. The WLC and AHP data were similar, with an apparent alignment in the southwestern and northeastern areas of the watershed, with the higher ranks of mapped erosion potential aligned with the SWAT sediment yield ( Figure 5). Additionally, the AHP data highlighted hillslope areas in the headland areas of the watershed due to the experts' emphasis on topography.

Management-Priority Area Ranking
The final comparison was an aggregate of the mapped erosion-potential cells from the WLC, AHP, and ensemble model runs and the SWAT model sediment-yield data to the management areas of the Weeks Bay watershed (Figure 6). The average erosion-potential values and the SWAT sediment-yield data were summarized for each of the management areas for ranked prioritization. Observational trends in the summarized data exhibited spatial alignment in the WLC and AHP upper erosion-potential ranks with higher SWAT sediment-yield data for the watershed management areas. Both the qualitative AHP model and the numerical SWAT model agreed in mapping the management area

Management-Priority Area Ranking
The final comparison was an aggregate of the mapped erosion-potential cells from the WLC, AHP, and ensemble model runs and the SWAT model sediment-yield data to the management areas of the Weeks Bay watershed (Figure 6). The average erosion-potential values and the SWAT sediment-yield data were summarized for each of the management areas for ranked prioritization. Observational trends in the summarized data exhibited spatial alignment in the WLC and AHP upper erosion-potential ranks with higher SWAT sediment-yield data for the watershed management areas. Both the qualitative AHP model and the numerical SWAT model agreed in mapping the management area ranked with the highest erosion. Figure 6 shows the management area rankings for the WLC, AHP, and SWAT models with the ranking of each area based on calculated erosion or sediment yield for the respective model output. The WLC, AHP, and SWAT ranks for the higher ranks of erosion aligned approximately 70% of the time.
Water 2022, 14, x FOR PEER REVIEW 18 of 28 ranked with the highest erosion. Figure 6 shows the management area rankings for the WLC, AHP, and SWAT models with the ranking of each area based on calculated erosion or sediment yield for the respective model output. The WLC, AHP, and SWAT ranks for the higher ranks of erosion aligned approximately 70% of the time.
(a) Similarly, the mapped erosion potential from the ensemble model was averaged and summarized for the management areas, and similar rankings were found to that of the SWAT model. Three of the four management areas-Pensacola Branch, Waterhole Branch, and Turkey Branch-aligned with areas ranked by the SWAT model for high sediment yield. In the lower ranks, alignment was found with the Weeks Branch and Weeks Bay management areas (Table 5). Overall, the ensembled model produced similar results Similarly, the mapped erosion potential from the ensemble model was averaged and summarized for the management areas, and similar rankings were found to that of the SWAT model. Three of the four management areas-Pensacola Branch, Waterhole Branch, and Turkey Branch-aligned with areas ranked by the SWAT model for high sediment yield. In the lower ranks, alignment was found with the Weeks Branch and Weeks Bay management areas (Table 5). Overall, the ensembled model produced similar results to that of the WLC and AHP models, with agreement in mapped erosion potential in areas of higher ranks.

Discussion
The erosion potential for a watershed along the northern Gulf of Mexico was mapped qualitatively using layers representative of physical erodibility, land sensitivity, and 30-year precipitation for the Weeks Bay watershed. The criteria used to define these layers were based on regional availability and input from watershed managers and stakeholders on erosional trends in this area. The approach used a WLC model and was set up with criteria similar to numerical models such as RUSLE [19]. Other approaches to soil-erosion mapping include empirical, conceptual, physically based, and hybrid models [60]. Some hydrodynamic numerical models, i.e., 1-D (one-dimensional) and 2-D hydrodynamic models showed better efficiency in urban-flood-risk mapping [61]. Compared to the conceptual models, i.e., the Hydrologic Simulation Program (HSPF) and SWAT [60], this study proposed a WLC model which is also applicable to larger areas with less complexity. Additionally, the AHP model in the study area considered the factors' importance in terms of being responsible for regional watershed erosion. However, this study did not show any physically based model, which could be a future demand for better estimation of erosion risk at a large scale. In addition, factors such as surface hydrology, slope aspect, and storm events could be used in 2-D physically based simulation models such as GSSHA (Grided Surface/Subsurface Hydrologic Analysis), DWSM (Dynamic Watershed Simulation Models), etc.
In this study, the northern headland and the southern agricultural regions of the watershed had the highest erosion potential, as expected [1,28]. Overall erosion potential in the Weeks Bay watershed tends to be lower in densely vegetated riparian and marsh areas. Many of these areas, especially in the southern area near the bay, are part of the Weeks Bay National Estuarine Research Reserve. Areas in the watershed with higher erosion potential are more associated with transitional-type lands that appear to be more agricultural or dynamic in terms of land practices. The southeast region of the watershed reflects this, as it is an area dominated by agricultural practices such as cultivated crops and turf-grass farms.
Land sensitivity accessed by soil erodibility (as defined by the soil K-factor and soil exposure or brightness) was the criterion that had the most influence on mapped erosion potential with the WLC model. The K-factor is used in USLE and RUSLE applications that represents soil texture and composition [19,20,24], and was the most influential of the layers. Soil brightness is indicative of disruptive land uses and increases in erosion potential from this variable are apparent in agriculture-dominated areas of the watershed. The physical erodibility (topographic) criteria for slope and stream density moderately influenced the WLC-modeled erosion potential. These areas of increased potential are indicative of higher concentrations of stream reaches, with more surface interaction with runoff waters and lower soil infiltration rates [55], and proximity to active stream and river channels.
The assessment was enhanced with expert input through the AHP model, allowing experts to rank the criteria. The ranking of criteria quantified the weights based on their relative importance in the study area. The physical erodibility criterion of slope was identified as the most important by the experts. In this, the AHP model differs from the traditional soil-erosion models, with the experts' minimal emphasis on land cover or land sensitivity [19,21]. This difference is also concerning due to the increased erosion rates due to development in the watersheds in this region of the Gulf of Mexico [52]. The AHP model proved to be beneficial in mapping areas of increased erosion potential, as defined by the upper ranks in the classified data. These shifts of the mapped erosion cells from the WLC model align with similar approaches using MCDA techniques [42,55]. The AHP model is not a typical numerical model and can be better-suited to qualitative geospatial assessment for mapping erosion potential with MCDA.
The variations in the AHP weights are normally used to identify shifts (increasing and decreasing) in the mapped erosion-potential cells. Data outliers were used for the identification of areas of model alignment for the high ranks of erosion potential. This approach allowed for an analysis (management) mask to be generated for these areas that were consistently high, irrespective of the variation in criteria weights. The areas identified were generally associated with higher slopes, which were associated with stream and channel networks. The outlier analysis was successful in helping to identify management areas; however, a suitable model approach would allow for a quantification shift between ranks of erosion potential [34,39,45].
The WLC and AHP erosion-potential models were run as an ensemble to improve the reliability of the output. The outputs (from models) were compared with the SWAT sediment yield; the comparisons were used to help identify if there were any visual alignments or trends between the qualitative and quantitative outputs. The comparisons were very similar, with the primary difference being the focus of higher erosion-potential values in the WLC and AHP models. The alignment of the data was most apparent in areas of transition with expanding development and agricultural land practices. The AHP model, however, produced some areas that were more focused on topographic features because of the experts' input placing more emphasis on physical erodibility (slope and other terrain measures). This was most apparent in the headland area of the watershed where the landscape has more dissection. The comparison of prioritization and ranking of the mapped erosion-potential cells for Weeks Bay Watershed management areas displayed alignment for select areas of higher erosion potential.
These management areas include the Pensacola Branch, Waterhole Branch, and Turkey Branch basins of the Weeks Bay Watershed. Each of these management areas have experienced increased erosion due to expanding land development associated with the surrounding communities. The management areas with the lowest ranks were also in agreement. This, coupled with the alignment in the upper ranks, indicates a generalized agreement between the qualitative and quantitative assessments for the prioritization of management areas. The approach has limitations in management areas that would be considered to be of intermediate concern with regard to erosion. This is, in part, due to various reasons, including the added emphasis by the experts on terrain characteristics and temporal generalizations in layers used in the WLC/AHP models and the SWAT model.
However, these areas of upper and lower erosion ranks are indicative of the agreement between the qualitative and quantitative modeling approaches, supporting the use of MCDA for management decisions and improved applications for watershed managers.

Conclusions
The watersheds that drain areas along the northern Gulf of Mexico have dynamic landscapes that experience erosion and contribute volumes of sediment to the associated estuaries. These watersheds and their estuaries are important for both their services and the resources they provide. To maintain the function and value of these services and resources, these areas require proper management in terms of soil-loss and erosion. Those involved in the management of these watersheds and the associated resources need to have information that is both accurate and timely. This study took an approach that focused on a qualitative assessment with GIS and MCDA for mapping erosion potential, to facilitate the prioritization of watershed management areas for improved management decisions. This aim of this approach was to provide watershed managers with a process to quickly map erosion potential and prioritize areas for management. This will allow for better and more efficient allocation of resources to be utilized in watershed management for areas such as the Weeks Bay watershed.
Field measurements and numerical models are critical to accurately measure and estimate sediment yield and erosion, as they provide quantitative information to facilitate management plans and decisions. Qualitative assessments can produce similar results and help in the management process. Qualitative assessments do not provide numerical sediment-yield information, but often provide more rapid assessments with a more simplistic approach and execution. The design of these assessments needs to be similar to numerical modeling approaches, using similar criteria and inputs. Their design can also allow for expert input for situation-specific applications due to their understanding of the processes and issues unique to the watershed of interest. A WLC model and an AHP model were used to map erosion potential based on terrain slope, geomorphology, land cover, soil erodibility, and long-term precipitation trends.
The WLC and AHP models developed for this study mapped erosion potential to cells as defined by the input layers for the Weeks Bay watershed. The mapped erosion potential aligned with the erosion trends described in the Weeks Bay Watershed Management Plan, with increased erosion occurring in areas associated with agricultural practices and expanding areas of urban and suburban development. The MCDA with the AHP model mapped areas that are most susceptible to erosion, as evidenced by the shifts of cells in the upper ranks. This, coupled with the analysis mask generated by the criteria weighting variations, identified areas of alignment or commonality in erosion potential in the Weeks Bay watershed. These areas, in both the WLC and AHP models, were in agreement with the SWAT sediment-yield data from the Weeks Bay Watershed Management Plan. The qualitative approach was effective in prioritizing management areas in the Weeks Bay watershed and offers a simplified approach to mapping erosion potential. This simplified approach provides watershed managers with the means to define and prioritize management actions. This does not replace the need for numerical modeling for quantitative soil-erosion metrics; it does, however, provide management alternatives when needed. There are numerous pathways for future research for this work. Refinement and adaptation of the qualitative approach will continue to improve reliability as compared with numerical model outputs such as sediment yield. This will involve more interactions with modelers and watershed managers (and stakeholders). Scaling the approach to a regional level in future efforts could be beneficial in prioritizing modeling needs and guidance in watershed management plan needs. This, coupled with more development of geospatial tools, would help to transition these types of approaches to a more operational environment, allowing for enhanced planning and management-scenario development. Funding: This paper was supported through the Geospatial Education and Outreach Project funded by the National Oceanic and Atmospheric Administration Regional Geospatial Modeling Grant, award # NA19NOS4730297.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript; or in the decision to publish the results.