Morphometric Analysis for Soil Erosion Susceptibility Mapping Using Novel GIS-Based Ensemble Model

The morphometric characteristics of the Kalvārı̄ basin were analyzed to prioritize sub-basins based on their susceptibility to erosion by water using a remote sensing-based data and a GIS. The morphometric parameters (MPs)—linear, relief, and shape—of the drainage network were calculated using data from the Advanced Land-observing Satellite (ALOS) phased-array L-type synthetic-aperture radar (PALSAR) digital elevation model (DEM) with a spatial resolution of 12.5 m. Interferometric synthetic aperture radar (InSAR) was used to generate the DEM. These parameters revealed the network’s texture, morpho-tectonics, geometry, and relief characteristics. A complex proportional assessment of alternatives (COPRAS)-analytical hierarchy process (AHP) novel-ensemble multiple-criteria decision-making (MCDM) model was used to rank sub-basins and to identify the major MPs that significantly influence erosion landforms of the Kalvārı̄ drainage basin. The results show that in evolutionary terms this is a youthful landscape. Rejuvenation has influenced the erosional development of the basin, but lithology and relief, structure, and tectonics have determined the drainage patterns of the catchment. Results of the AHP model indicate that slope and drainage density influence erosion in the study area. The COPRAS-AHP ensemble model results reveal that sub-basin 1 is the most susceptible to soil erosion (SE) and that sub-basin 5 is least susceptible. The ensemble model was compared to the two individual models using the Spearman correlation coefficient test (SCCT) and the Kendall Tau correlation coefficient test (KTCCT). To evaluate the prediction accuracy of the ensemble model, its results were compared to results generated by the modified Pacific Southwest Inter-Agency Committee (MPSIAC) model in each sub-basin. Based on SCCT and KTCCT, the ensemble model was better at ranking sub-basins than the MPSIAC model, which indicated that sub-basins 1 and 4, with mean sediment yields of 943.7 and 456.3 m3km−2 year−1, respectively, have the highest and lowest SE susceptibility in the study area. The sensitivity analysis revealed that the most sensitive parameters of the MPSIAC model are slope (R2 = 0.96), followed by runoff (R2 = 0.95). The MPSIAC shows that the ensemble model has a high prediction accuracy. The method tested here has been shown to be an effective tool to improve sustainable soil management.

. Location of the study area in Iran.

Methodology
There are seven steps to this method (Figure 2), which are: (1) determine the boundaries of 11 sub-basins using topographic maps (1:50,000 scale and 20 m contours) and PALSAR DEM (12.5 m resolution); 2) extract basic MPs (basin area, basin length, numbers and lengths of streams of each order, basin perimeter, and bifurcation ratio) of the Kalvārī Basin from drainage networks, sub-basin polygons, and an elevation model generated from the PALSAR DEM; 3) extract linear, shape, and relief factors, including mean bifurcation ratio (Rbm), drainage density (Dd), stream frequency (Fu), texture ratio (T), length of overland flow (Lo), infiltration number (If), constant of channel maintenance (C), form factor (Rf), shape factor (Bs), elongation ratio (Re), compactness coefficient (Cc), circularity ratio (Rc), ruggedness number (Rn), basin relief (Bh), relief ratio (Rh), and slope (S) from the basic parameters and formula (Table 1); 4) apply the AHP model to determine the relative importance of the morphometric factors in SE; 5) calculate the relative weight of each alternative (for the sub-basins) and prioritize them using the COPRAS model; 6) prepare an SE susceptibility map using the AHP-COPRAS ensemble model; and 7) validate the results using the MPSIAC model and non-parametric correlation tests such as the Spearman correlation coefficient test (SCCT) and the Kendall Tau correlation coefficient test (KTCCT).
Morphometric (elevation and slope) and hydrologic (drainage networks) parameters were estimated using an ALOS PALSAR DEM (12.5 m resolution). Detailed descriptions of the methodology used to produce the ALOS PALSAR DEM using InSAR are discussed in [51,52]. The key step in DEM generation is the transformation of phase-change measurements to elevations [53]. InSAR was developed by Graham [54]. In terms of economics, efficiency, and resolution, InSAR is the best DEM production technique [55]. InSAR calculates phase changes from two radar image pairs taken at different times and reveals changes in the surface quantitatively and qualitatively [56]. There are 6 steps in InSAR generation of DEMs from PALSAR data [52] (Figure 3), which are: (1) register radar images-in this study, two images (Slave (2008/07/11) and Master (2010/08/16)) were entered chronologically; (2) produce interferogram-the images iSRTMn SARScape4.8 are converted into

Methodology
There are seven steps to this method (Figure 2), which are: (1) determine the boundaries of 11 sub-basins using topographic maps (1:50,000 scale and 20 m contours) and PALSAR DEM (12.5 m resolution); (2) extract basic MPs (basin area, basin length, numbers and lengths of streams of each order, basin perimeter, and bifurcation ratio) of the Kalvārī Basin from drainage networks, sub-basin polygons, and an elevation model generated from the PALSAR DEM; (3) extract linear, shape, and relief factors, including mean bifurcation ratio (Rbm), drainage density (Dd), stream frequency (Fu), texture ratio (T), length of overland flow (Lo), infiltration number (If), constant of channel maintenance (C), form factor (Rf), shape factor (Bs), elongation ratio (Re), compactness coefficient (Cc), circularity ratio (Rc), ruggedness number (Rn), basin relief (Bh), relief ratio (Rh), and slope (S) from the basic parameters and formula (Table 1); (4) apply the AHP model to determine the relative importance of the morphometric factors in SE; (5) calculate the relative weight of each alternative (for the sub-basins) and prioritize them using the COPRAS model; (6) prepare an SE susceptibility map using the AHP-COPRAS ensemble model; and (7) validate the results using the MPSIAC model and non-parametric correlation tests such as the Spearman correlation coefficient test (SCCT) and the Kendall Tau correlation coefficient test (KTCCT).    Morphometric (elevation and slope) and hydrologic (drainage networks) parameters were estimated using an ALOS PALSAR DEM (12.5 m resolution). Detailed descriptions of the methodology used to produce the ALOS PALSAR DEM using InSAR are discussed in [51,52]. The key step in DEM generation is the transformation of phase-change measurements to elevations [53]. InSAR was developed by Graham [54]. In terms of economics, efficiency, and resolution, InSAR is the best DEM production technique [55]. InSAR calculates phase changes from two radar image pairs taken at different times and reveals changes in the surface quantitatively and qualitatively [56]. There are 6 steps in InSAR generation of DEMs from PALSAR data [52] (Figure 3), which are: (1) register radar images-in this study, two images (Slave (2008/07/11) and Master (2010/08/16)) were entered chronologically; (2) produce interferogram-the images iSRTMn SARScape4.8 are converted into single-look complex (SLC) format with the following settings, TB (745 days), SB (114.8 m), and critical baseline (2623.5 m); (3) remove flat effect-the DEM derived from the Shuttle Radar Topography Mission (SRTM) and its Doppler Orbitography and Radiopositioning Integrated by Satellite (DORIS) files are used to remove the topographic effect; (4) filter noise-adaptive filters are used to eliminate the effect of noise from the interferogram, as they may reduce its quality (this step will also produce a coherent image, coherence indicates interferogram quality and should be less than 0.5; the Goldstein filter was used in this study and the coherence value of the image was 0.20); (5) unwrap phase-region-growing and minimum cost-flow algorithms are the most common methods for phase correction, a region-growing algorithm was used in this study; and (6) convert phase to height-the ALOS PALSAR DEM is created by phase-to-height conversion, the PALSAR DEM of the study area was produced with the mosaicking tools in ENVI v4.8.   An ALOS PALSAR DEM enhanced reproduction of the complex topography of the study area to improve approximations of morphometric and hydrological factors [57,58]. The vertical accuracy of the ALOS DEM was assessed by ground-truthing the model's elevation values in ArcGIS 10.5 at selected 230 ground control points (GCPs) (Figure 4), following procedures used in [59]. The root mean square error (RMSE) of DEM generated based on comparisons with GCPs was 1.2 m.   The boundaries of the sub-basins ( Figure 5) were extracted by determining the pour point, the location where water drained from the basin flows into the main river. For drainage network extraction, the Arc Hydro package was used. Compared to the manual approach, this package provides a rational, effective, and consistent algorithm [46]. Generating the drainage network using Arc Hydro has been explained by Ahmad Rather et al. [47]. To do this, DEM sinks were specified and filled to designate the flow direction and locations of accumulation. Stream networks in the sub-basin were defined according to the cumulative number of upstream cells draining into each cell. A threshold of greater than 500 was used to extract the drainage; this critical threshold is the minimum upstream area necessary to produce a stream. The areas and perimeters of the sub-basins were extracted by computing the geometry of the sub-basin polygons. Strahler's scheme was used for stream ordering [16].      [60]. The COPRAS method assumes direct and commensurate affiliations of the levels of magnitude and usefulness of alternatives in the presence of conflicting criteria [61]. The COPRAS procedure consists of the following steps [62]: Step 1: Prepare the primary matrix; Step 2: Normalize the primary matrix using Equation (1): where x ij is the normalized quantity of the j th criterion, x ij is the i th alternative performance of the j th criterion, and m denotes the alternative numbers; Step 3: Determine the normalized weighted decision-making matrix (Equation (2)): where x ij is the efficiency of the i th alternative, and w j is the criterion weight; Step 4: Compute the maximum and minimum indices for alternatives-in this step, alternatives are classified as maximising and minimising indices (Equations (3) and (4)): where y +ij and y −ij are the weighted normalized qualities for advantageous and non-advantageous adjectives, respectively. In fact, the highest value of parameters that have a direct relationship with SE, such as slope, and the lowest value of parameters that have an inverse relationship with SE, such as shape factors, is y +ij and vice versa; and Step 5: Calculate the relative weights of each alternative using (Equation (5)): where S − min is the minimum value of S − j . S + j and S − j are maximum and minimum indices, respectively.

AHP
Several methods can characterize the weights of criteria, but in this study the AHP was used. Weights were calculated with a pair-wise variable comparison matrix developed from experts' opinions. For this purpose, an AHP questionnaire was designed and was administered to 18 geomorphology and 15 hydrology experts. Initially, due to the incompatibility of some of the paired-comparison matrices from the experts' votes, the questionnaire was redistributed to confirm the matrices' compatibilities and the validity of the questionnaire. The expert judgments reflect a blend of rational thinking and experience [63]. Based on the AHP method, Saaty's linguistic scales (Table 2) of pair-wise comparisons were converted to quantitative values [64]. Then, the weights of criteria were determined using Equations (6) and (7) [63]: where W j is the weight of criteria by AHP, n ij is normalized of pair-wise comparison matrix and a ij is matrix element in row i and column j. Table 2. Saaty's linguistic scales in the analytical hierarchy process (AHP) [64].

Preference Factor
Degree of Preference 1 Equally 3 Moderately 5 Strongly 7 Very strongly 9 Extremely 2, 4, 6, and 8 Intermediate between 2 adjacent judgments The consistency ratio (CR) is the mechanism by which the validity of the expert response is measured in the pair-wise comparison matrix [63]. A CR < 0.1 is acceptable. Equations (8) to (12) were used to calculate CR [4]: where CR is consistency ratio, CI is consistency index, RI is a random index (extracted from Table 3), n is the number of criteria, λ max is the largest special matrix value, λ is consistency vector, WSV is weighted sum vector, A is pair-wise comparison matrix, and W is weight of criteria vector.

Non-Parametric Correlation Tests to Comparing Models Ranking
The non-parametric Spearman correlation coefficient test (SCCT) and Kendall Tau correlation coefficient test (KTCCT) were used to compare the ranks of the observed values of two independent variables within the models instead of comparing their values [65,66] to determine whether the variables are statistically dependent. A reciprocal comparison is made between each random pair of variables. The number of comparisons is equal to where n is the number of alternatives.
Remote Sens. 2020, 12, 874 10 of 24 KTCC (Equation (13)) is used to compare two variables with dissimilar ranks. When two variables have similar rankings of observed values, Equation (14) is used: where C and D are the numbers of concordant pairs and the number of discordant pairs, respectively. T and U are the numbers of pairs having the same ranks within the two data sets. The non-parametric SCCT test compares the ranked values of two variables. Equation (15) is used when the observations of two variables are never similarly ranked and Equation (16) is used when one set of observations for two variables have the same rank: where d i is the difference between the ranks of models for each alternative, and x and y are the mean of x and y model, respectively.

Validation of Results Using MPSIAC Model
To validate the model, the sediment-delivery ratio (SDR) was calculated for each sub-basin using the modified Pacific Southwest Inter-Agency Committee (MPSIAC) method. The PSIAC method was created in 1968 by the US Water Management Committee [67]. MPSIAC was used in Walnut Gulch basin in southeastern Arizona. In 1982, Johnson and Gembhart quantified the descriptive concepts of the first model and presented each of the factors mathematically. In this method, the effects of nine important effective parameters-surface geology (X1), soil (X2), climate (X3), runoff (X4), topography (X5), land cover (X6), land use (X7), upland SE (X8) and channel SE (X9)-are evaluated. Depending on the relative importance of each parameter, values are attributed and the sum of the values of each is used to estimate SE severity and sediment yield [68]. Details of the MPSIAC model have been reported elsewhere [69][70][71].

Estimation of Sediment Yield and Total Sediment Production
Each effective parameter is divided into classes. Based on the estimated impact of each class, augmentation values are assigned by consulting model tables. The SE severity and annual sediment yield are estimated by summation of the values. This is signified by R: where R is ranking value (m 3 km −2 year −1 ) and Xi is each factor in the model. To manage the accuracy of interpolations and extrapolations of the nine parameter values in the MPSIAC model, this equation was developed to estimate sediment yield [70]: where Qs is rate of sediment yield (m 3 km −2 year −1 ) and R is ranking value (m 3 km −2 year −1 ). To calculate the total sediment production of a study area, the rate of sediment yield is multiplied by surface area: where S is the total sediment production based on sediment yield (m 3 year −1 ), Qs is the sediment rate (m 3 km −2 year −1 ), and A is surface area.

Analysis of MPs
The characteristics of physical processes in a basin drainage system significantly impact its infiltration capacity and runoff dimensions [14]. Basin morphometry describes the relationships between the geomorphic, hydrologic, and geologic surface processes and a landscape [16]. Quantitative analysis of the Kalvārī Basin and its 11 sub-basins was carried out to evaluate the basin's morphometric characteristics and the characteristics of each sub-basin drainage network. This analysis enables prioritization of the variables in terms of conservation and management efforts. In this regard, 23 MPs that represent basic (Table 4), linear, shape, and relief characteristics of the basin were examined.

Basic Parameters
The Kalvārī Basin is a 5th-order ( Figure 6) basin with an area of 70.71 km 2 . The basin contains 2204 stream segments. The total length is 327.93 km. The number and length were determined by the drainage threshold defined during the extraction of the stream network from the DEM. These characteristics are indicative of Horton's First Law [15], which states that the number of streams of different ranks in the basin tends to have an inverse geometric ratio. This inverse geometric relationship is shown in the form of straight lines when the logarithm of the number of streams is plotted on a regular graph (Figure 7a). Changes in stream rankings are strongly dependent on the morphological and structural features of the basin. The number and length of streams vary directly with the size of the sub-basins. The lengths of the stream (Lu) were calculated according to Horton's law. The length indicates the temporal development of a stream interacting with tectonic disturbances. The higher the stream ranking, the greater the length of the stream. The characteristics of the basin length of the Kalvārī Basin conform to Horton's Second Law [15], which states that the average flow length of each stream tends to have a straight geometric ratio relative to basin length. This linear geometrical relationship is shown when the logarithms of the basin length values are plotted on a regular graph (Figure 7b). The lengths of streams in the basin vary from 4.62 km (sub-basin 6) to a maximum of 91.62 km (sub-basin 7). Mean stream length ranges from 0.1216 (sub-basin 6) to 0.1877 (sub-basin 8). The minimum and maximum heights of the basin are 1100 m and 1999 m. A comparison of the areas and lengths of sub-basin waterways also reveals a direct relationship (Figure 7c).

Linear Parameters
Stream density reflects landscape dissection and a basin's runoff potential. Slope angle and relative relief are the primary morphological factors that control drainage density. Strahler [16] concluded that Dd is low when basin relief is high. Lower Dd in a basin indicates highly permeable subsurface material, good vegetation, and low roughness, and the opposite conditions produce high

Stream Density (Dd)
Stream density reflects landscape dissection and a basin's runoff potential. Slope angle and relative relief are the primary morphological factors that control drainage density. Strahler [16] concluded that Dd is low when basin relief is high. Lower Dd in a basin indicates highly permeable subsurface material, good vegetation, and low roughness, and the opposite conditions produce high Dd [72]. In the Kalvārī Basin, the lowest stream density was observed in sub-basin 4 (3.94 km). This sub-basin, therefore, has the highest infiltration among all sub-basins. If the condensation parameter alone is considered, it has the highest SE resistance and is the least susceptible to SE, and sub-basin 1 has the highest value. Sub-basins 10 (5.12), 3 (5.08), 9 (4.9), 5 (4.78), 11 (4.55), 7 (4.43), 8 (4.13), 2 (4.06), and 4 (3.94) are in the following categories of SE susceptibility (Table 5).

Stream Frequency (Fu)
Fu is the ratio of the number of streams in a basin to that basin's area [15,73]. The Fu is inversely related to infiltration and is directly related to basin roughness [73]. High Fu indicates that the basin has a rocky surface and low permeability that contributes to further SE, and vice versa. Values of Fu in the study area vary from 22.15 streams/km2 for basin 8 to 41.06 for basin 3. Thus, sub-basin 3 has the lowest absorption capacity and is the most susceptible to SE, whereas basin 8 is least susceptible to SE. Sub-basins 9 (34), 1 (33.09), 5 (33.01), 7 (32.34), 10 (32.05), 2 (29.57), 6 (27.08) 11 (25.41) and 4 (24.4) rank second to tenth in susceptibility.

Mean Bifurcation Ratio (Rbm)
Rbm indicates that the infiltration of the basin are inversely correlated. A high Rbm value is the peak of the initial hydrograph when flooding results in high soil degradation. Rbm values are too high for all sub-basins in this study area, indicating that they are structurally complex and have low infiltration rates. Rbm values range from 3.66 for sub-basin 3 to 9.96 for sub-basin 5. Sub-basin 5 is the basin most susceptible to SE and sub-basins 6, 4, 7, 2, 1, 10, 8, 11 and 9 with values (8.96, 6.48, 6.01, 5.87, 5.81, 5.51, 4.56, 4.17, and 3.96) follow in rank.

Constant of Channel Maintenance (C)
This indicator reflects infiltration and the control of flow to the basin outlet [19]. The relationship between this parameter and SE is analogous to the relationship of drainage density to stream frequency. The values of C range from a minimum (0.13) for sub-basin 1 to a maximum (0.26) for sub-basin 6. Sub-basin 6 is the most erodible sub-basins, and sub-basin 1 is not susceptible to SE. Sub-basins 4 (0.21), 2 (0.

Elongation Ratio (Re)
The values of Re range from 0.6 to 1.0 based on climate and geological conditions [20]. Values of about 1 are typical of areas with very low roughness [15], mild topography, and little frictional resistance to flow, while values of 0.6 to 0.8 are associated with high roughness and steep terrain [75]. In the study area, sub-basin 9 is most the elongated sub-basin (0.78) and is therefore least susceptible to SE. Sub-basin 7 has the smallest elongation ratio (0.66) and is the most susceptible to SE. Sub-basins 5 (0.67), 10 (0.69), 2 (0.7), 11 (0.73), 4 (0.75), 1 (0.77), 8 (0.77), 6 (0.8) and 3 (0.79) fill out the top ten in terms of susceptibility.

Circularity Ratio (Rc)
Rc relates to several basin characteristics: stream length and frequency, geological structure, climate, roughness, and slope. High Rc indicates a circular basin with moderate to high roughness, high infiltration, less elongated, lower roughness, and low infiltration. Sub-basin 9, with the lowest Rc value (0.34), is least susceptible to SE due to high infiltration. Sub-basin 6, having the highest Rc value (0.67), is most susceptible to SE. In terms of Rc and susceptibility to SE, sub-basins 5 (0.

Compactness Coefficient (Cc)
The Cc is directly linked to infiltration capacity. Therefore, the relationship between Cc and SE is the same as Rf and Bs. Sub-basin 6 has the lowest Cc value (1.14) and has surfaces with low permeability. Therefore sub-basin 6 is the most SE-susceptible basin. Sub-basin 9 has the highest Cc

Basin Relief (Bh)
Bh indicates height difference [76]. This parameter has a significant role in hydrological characteristics [6]. The Bh indicates the overall slope of a basin and therefore the intensity of the SE forces operating on the slopes. The relationships of this parameter to SE is the same as Dd, and

Ruggedness Number (Rn)
Rn is used to calculate the flood potential of streams [16]. This parameter reflects the geometrical characteristics of the basin. The Rn is directly related to erodibility: increasing Rn increases erosivity. Rn ranged from a minimum of 0.079 for sub-basin 4 to a maximum of 5.41 for sub-basin 1. Sub-basin 1 is therefore most susceptible, followed by sub-basins 2, 10,8,7,3,11,9,4, and 6.

Prioritization of Sub-Basins Using Novel AHP-COPRAS Ensemble Model
To determine the relative importance of the contribution that each MP makes to determining soil-SE potential, a set of academic experts were asked to express their informed opinions of expected importance of each variable. From their responses, a pairwise comparison matrix was created to determine the weight of each parameter using AHP ( Table 6). The consistency ratio from this matrix was 0.05, indicating that the opinions were consistent. Based on the AHP results (Table 6 and Figure 9). Accordingly, sub-basin 4 is in the very low-SE class. Sub-basins 5, 6, and 9 are in the low-SE group. Sub-basins 7, 10, and 11 have moderate SE potential. Sub-basin 2 is highly likely to be prone to SE, and sub-basins 1, 3, and 8 are very highly likely to experience SE.      The COPRAS algorithm indicated ( Table 7) that sub-basins 1 (0.319), 8 (0.257), and 3 (0.208) are the most susceptible to SE. Sub-basins 2, 10, 11, 7, and 9 rank next in terms of SE susceptibility. And sub-basins 4 (0.118), 6 (0.137), and 5 (0.147) are the least susceptible to SE.

Validation of Results
Results of mean sediment yield from the MPSIAC model (Table 8) Table 9). Comparing the results of sub-basin prioritization using the hybrid method and the MPSIAC method showed that the ensemble method prioritizes the sub-basins rather accurately (100%) (Figure 10). Non-parametric tests indicate that, compared to the AHP and COPRAS models, the AHP-COPRAS ensemble had the best correlation (Table 10).   Figure 9. Prioritization of sub-basins for conservation programs.

Validation of Results
Results of mean sediment yield from the MPSIAC model (Table 8) Table 9). Comparing the results of sub-basin prioritization using the hybrid method and the MPSIAC method showed that the ensemble method prioritizes the sub-basins rather accurately (100%) (Figure 10). Non-parametric tests indicate that, compared to the AHP and COPRAS models, the AHP-COPRAS ensemble had the best correlation (Table 10).

Discussion
The Kalvārī Basin was selected for study because of the extent to which it is impacted by SE. Erosion processes in the basin reflect its morphology-visible in its linear relief pattern, its shape, and its spatial extent. After determining its drainage network and MPs using RS, GIS, and a DEM, the AHP-COPRAS ensemble model was used to develop a map of SE susceptibility and to rank the sub-basins by the intensity of SE. The results of the map were compared to the results from two individual models.
Digitally estimating MPs provided an easier, more accurate, and more quantitative method to evaluate morphometric characteristics and to analyze the variations within the region. PALSAR DEM has been used previously for morphometric studies, especially in mountainous regions, as it provides more accurate elevation measurements and better morphometric and geomorphic details (compared to ASTER and SRTM DEMs) [77].
First-order streams comprise nearly 67% of the streams in the study area. That first-order streams are the most numerous of all orders indicates that there is a structural weakness in the basin in the form of lineaments [78]. Stream frequency in all sub-basins is moderate, which indicates that there is moderate run-off intensity and high drainage density. The terrain is highly dissected. Sub-basins are finely textured, indicating highly developed drainage and high rates of SE. Stream frequency, drainage density, and texture are very high in sub-basins 1 and 3, which indicates that these sub-basins experience intensive SE. The values of these three characteristics are lowest in sub-basin 4, a unit that is less prone to SE. Relief ratios generally depict a high-energy basin with substantial SE and high sediment loads. This parameter is very high in sub-basins 1 and 8, where SE rates are high, and slopes are moderate. By contrast, the relief ratios are very low in sub-basins 4 and 5, low-energy drainages as exhibited by their gentle topographies.
Despite the low mean slope for the region overall, both the maximum and minimum slope angles of the Kalvārī Basin are in sub-basin 4, where the gentle topography has attracted greater human activity. This has created a plateau-like terrain-a flat platform with steeply sloping edges on its perimeter. This unit is a low-energy landscape with less sediment in streams. Sub-basins 1, 11, and 8 have higher mean slopes, and provide favorable topography for higher SE potential. The ruggedness numbers of all the sub-basins reflect the mountainous nature of the region's terrain, with significant amounts of SE and dissection. The lowest ruggedness value is found in sub-basin 4 and the highest is in sub-basin 8; these compliment the patterns of slope and relief ratios discussed above. Shape factors reflect the geomorphology of a landscape and provide evidence of run-off and infiltration processes [76]. Sub-basins 6, 3, 9, and 1 have the highest form factors of the sub-basins, reflecting their more ovate and less elongated shapes (these forms are consistent with the elongation measures and circulatory ratios for those units). The higher form-factor values of these sub-basins indicate more drainage development and more structural control.
Based on the results of the AHP model, we know that slope, drainage density, infiltration number, and texture ratio significantly impact SE in the study area, results which are consistent with others [12,[79][80][81][82]. Arabameri et al. [82] used four MCDM-based models to rank the sub-basins' SE susceptibilities by analyzing the MPs. The mean bifurcation ratio (Rbm), slope (S), and infiltration number (If) have key roles in SE rates. Slope is a morphometric factor associated with hydrology that indicates runoff volumes and runoff concentration time [76]. Soil infiltration-capacity and the initial resistance of a surface to SE depend on drainage density [73]. Infiltration number is very important for expressing a basin's infiltration characteristics and depends directly on the basin's runoff capacity [76]. Validation of the results by using non-parametric tests shows that the ensemble model performed better than the AHP and COPRAS individual models. This is consistent with the findings of Arabameri et al. [83] which indicated that ensemble MCDM-based models perform better than individual models.
Validation of the ensemble model with the MPSIAC model shows that the ensemble model very accurately prioritized sub-basins and could be used to prepare SE susceptibility maps. Sub-basins 1, 3, and 8 might be experiencing heavy SE and sedimentation due to either higher erodibility of hills Remote Sens. 2020, 12, 874 20 of 24 in these areas, greater slopes, or land-cover changes (from compact to less compact pastures). The sub-basins with slight and moderate SE and sedimentation may be less affected because of ultrabasic and crystallized limestone formations, rock outcrops, or cultivated lands. The results provide ways to identify very erosive areas and open new horizons in watershed management and sediment control by providing greater evidence to support conservation-project prioritization.
Compared to other commonly used approaches to produce SE susceptibility maps, this ensemble approach can be achieved with much simpler input data, specifically morphometric data, that can be easily extracted from DEMs. Our approach does not need other SE parameters that might necessitate a soil-inventory map, a source that is very time-consuming and expensive create. One limitation of this method is that it ignores human activities that interfere with hydrologic processes (e.g., reducing a river's flow capacity), actions that often amplify SE. Moreover, additional information derived from land use analyses that can be used to detail soil-erosion history would be beneficial, as they would clarify the connection between spatial patterns of activities and SE zones. But suitable results can be achieved without the best information.

Conclusions
The common approaches to assess the problem of SE, to develop management plans to address it, require quantitative SE rate data at regional and global scales. In this research, the morphometric characteristics of the Kalvārī Basin were determined for use in an MCDM-based ensemble model. A PALSAR DEM was input into a GIS, allowing for the development of quantitative and qualitative morphometric analyses and the extraction of the spatial patterns of MPs. An AHP model was used to evaluate the importance of each of the 16 parameters and revealed that slope, drainage density, and infiltration number were the most important predictors of SE potential. A COPRAS-AHP ensemble model indicated that sub-basins 1, 3, and 8 are highly susceptible to SE. Validation of the modeled results using MPSIAC and non-parametric tests show that the ensemble model and the ranking of MPs achieved strong prediction accuracy of SE susceptibility and this enables a faster and less expensive prioritization of sub-basins for management actions. The method presented here takes advantage of computer-assisted extraction and computation of morphometric characteristics to predict the spatial patterns of the intensity of SE potential. A significant advantage is that only high-quality empirical topographical data are employed. In this way, one can produce assessments for catchments that either lack stream gages and extensive, large-scale records of past SE patterns, or are remote and difficult to access. This, in fact, describes most of the watercourses in Iran. Given the predicted high SE-susceptibility of sub-basins 1, 3, and 8, it is recommended that protective measures be taken to arrest SE, reduce sedimentation in reservoirs, stabilize slopes against mass wasting, and reduce flood risk.