Proposing a Novel Predictive Technique for Gully Erosion Susceptibility Mapping in Arid and Semi-arid Regions (Iran)

: Gully erosion is considered to be one of the main causes of land degradation in arid and semi-arid territories around the world. In this research, gully erosion susceptibility mapping was carried out in Semnan province (Iran) as a case study in which we tested the e ﬃ ciency of the index of entropy (IoE), the Vlse Kriterijumska Optimizacija I Kompromisno Resenje (VIKOR) method, and their combination. Remote sensing and geographic information system (GIS) were used to reduce the time and costs needed for rapid assessment of gully erosion. Firstly, a gully erosion inventory map (GEIM) with 206 gully locations was obtained from various sources and randomly divided into two groups: A training dataset (70% of the data) and a validation dataset (30% of the data). Fifteen gully-related conditioning factors (GRCFs) including elevation, slope, aspect, plan curvature, stream power index, topographical wetness index, rainfall, soil type, drainage density, distance to river, distance to road, distance to fault, lithology, land use / land cover, and soil type, were used for modeling. The advanced land observing satellite (ALOS) digital elevation model with a spatial resolution of 30 m was used for the extraction of the above-mentioned topographic factors. The tolerance (TOL) and variance inﬂation factor (VIF) were also included for checking the multicollinearity among the GRCFs. Based on IoE, we concluded that soil type, lithology, and elevation were the most signiﬁcant in terms of gully formation. Validation results using the area under the receiver operating characteristic curve (AUROC) showed that IoE (0.941) reached a higher prediction accuracy than VIKOR (0.857) and VIKOR-IoE (0.868). Based on our results, the combination of statistical (IoE) models along with remote sensing and GIS can convert the multi-criteria decision-making (MCDM) models into e ﬃ cient and powerful tools for gully erosion prediction. We strongly suggest that decision-makers and managers should use these kinds of results to develop more consistent solutions to achieve sustainable development on degraded lands such as in the Semnan province. gully locations and 15 environmental factors, namely elevation, slope, aspect, plan curvature, stream power index (SPI), topography wetness index (TWI), rainfall, soil type, drainage density, distance to river, distance to road, distance to fault, lithology, land use / land cover (LU / LC), and normalized di ﬀ erence vegetation index (NDVI). Our results showed that the IoE and expert knowledge models considered soil type, lithology, elevation, and land use to have the biggest impact on gully occurrence in the study area. The validation of our results using the AUROC curve showed that IoE has a higher prediction accuracy than the other models. The scientiﬁc methodology presented in this study can be used in areas that do not have accurate data on the spatial distribution of past and current gullies. The results of this research could be used by environmental planners to manage and reduce the risks associated with the development of gully erosion in the study area as well as in the conservation of natural resources such as water and soil.

In arid and semi-arid climates, such as the Semnan province, Iran, gullies cause very extensive damage and are considered a major environmental problem [26,46]. Therefore, the identification of the causes of gully erosion development and its zoning for comprehensive management planning by executive agencies remain essential. This will be definitive to the design of restoration strategies that can be developed following nature-based solutions to achieve sustainability [47]. Our literature review to date shows that despite the high sensitivity of this study area to gully erosion, no comprehensive studies have been carried out to identify areas that are particularly susceptible to gully erosion at all; although the scientific community are working very hard and achieving interesting and promising results. In this study, a new hybrid approach has been used to identify areas susceptible to gully erosion, which includes the combination of MCDM (multi-criteria decision making) and bivariate statistical methods. The main objective of this research is the development of a gully erosion susceptibility mapping in a representative area with known soil erosion processes; in this case, the Semnan province in Iran. To this end, the efficiency of the index of entropy (IoE) as a statistical model, Vlse Kriterijumska Optimizacija I Kompromisno Resenje (VIKOR) as multi-criteria decision-making model, and the combination of the two methods were assessed in terms of their suitability for prediction of areas prone to gully erosion.

Study Area
The Semnan watershed with an area of 8700 km 2 is located in the western part of Semnan province (52 • 37 19" to 54 • 06 25" E longitude and 34 • 52 40" to 35 • 51 04" N latitude) ( Figure 1). The minimum elevation of the study area is 743 m.s.a.l. the maximum elevation is 3266 m.s.a.l. and the mean elevation lies at 1184 m.s.a.l. The mean annual rainfall is 120 mm/year, and the average annual temperature is 21.8 • C, which is considered a dry and semi-arid climate [48]. The inclination of the study area varies from a minimum of 0 • to a maximum of 73 • , with mean values of 51.5 • . The topography of the Semnan watershed can be considered mostly flat in the central and southern parts, and the northern area has a mountainous landscape. The most important settlements in the study area are Semnan, Sorkheh, Mahdishahr, and Shahmirzad. The soils of the study are mainly in the group of plateau soils and can be classified as Inceptisols, Entisols, and Aridisols [49]. The lithology units in the study are mostly marls, gypsum, conglomerate, sandstone, shale, tuff, dolomite, limestone, fluvial conglomerate, volcanic rocks, clays, and saline areas [50]. The main land uses are cultivated areas, bare lands, kavir, mid-range, poor-range, urban, and waterbodies.
Morphometric analysis of gullies in the study area shows that about 8% of the study area could be affected by severe soil erosion, which would indicate a high sensitivity of the study area to the development of gully erosion processes. The abundance of gullies in the center and south of the study area with a low slope is higher, while, in contrast, the northern parts of the study area features more steep slopes and the rocky outcrops. A first overview shows that gully lengths range from several to hundreds of meters with depths reaching some meters. The width of the gullies also can vary from several centimeters to meters. We observed that the cross-section of the detected gullies in the northern parts is mainly V-shape, but in the central and southern parts, possibly due to the more erodible and lower slopes, they are U-shaped. In terms of the position of the gully in the landscape, most of the gullies are on the floor valley type. Considering the evolution of their morphologies [15], the gullies are mostly of a continuous type (drainage network), and the rest are discontinuous ones. In terms of the head-cut plan [10], the gullies in the study area could be classified as a digitized type. They have a maximum width of 8 m, a maximum depth of 12 m, and a maximum length of 360 m. The expansion of these corridors and the fall of their roofs are causing an increase in their sizes.
Remote Sens. 2019, 11, 2577 4 of 25 In terms of the head-cut plan [10], the gullies in the study area could be classified as a digitized type. They have a maximum width of 8 m, a maximum depth of 12 m, and a maximum length of 360 m. The expansion of these corridors and the fall of their roofs are causing an increase in their sizes.

Methodology
As shown in Figure 2, this research consists of five different steps. Firstly, the data were prepared, which consisted of 3 steps: (i) Determining the location of the gullies in the study area, preparing a gully erosion inventory map (GEIM) and dividing the data into two groups, namely training data (70%) and validation data (30%); (ii) selecting 500 sample points (alternatives) in the study area for applying the VIKOR model using the tool "Create random point" in ArcGIS10.5 (ESRI, USA); and, (iii) including the gully-related conditioning factors (GRCFs), which was divided into the four groups of topographical (slope, aspect, elevation, plan curvature, and topography wetness index (TWI)), geological (lithology and distance to fault), hydrological (stream power index (SPI), rainfall, drainage density, and distance to stream), and environmental (distance to road, soil type, normalized difference vegetation index (NDVI), and land use/land cover (LU/LC)) factors. Secondly, an assessment of the multicollinearity among the conditioning factors was performed. In the third step, the index of entropy (IoE) and the VIKOR (VIseKriterijumska Optimizacija I Kompromisno Resenje) models were applied. The fourth step aimed to introduce a novel hybrid model through the integration of the IoE and VIKOR models (IoE-VIKOR) for gully erosion susceptibility mapping (GESM). Finally, in the fifth step, a validation of the models was performed using the seed cell area index (SCAI) and the area under the curve (AUC).

Methodology
As shown in Figure 2, this research consists of five different steps. Firstly, the data were prepared, which consisted of 3 steps: (i) Determining the location of the gullies in the study area, preparing a gully erosion inventory map (GEIM) and dividing the data into two groups, namely training data (70%) and validation data (30%); (ii) selecting 500 sample points (alternatives) in the study area for applying the VIKOR model using the tool "Create random point" in ArcGIS10.5 (ESRI, USA); and, (iii) including the gully-related conditioning factors (GRCFs), which was divided into the four groups of topographical (slope, aspect, elevation, plan curvature, and topography wetness index (TWI)), geological (lithology and distance to fault), hydrological (stream power index (SPI), rainfall, drainage density, and distance to stream), and environmental (distance to road, soil type, normalized difference vegetation index (NDVI), and land use/land cover (LU/LC)) factors. Secondly, an assessment of the multicollinearity among the conditioning factors was performed. In the third step, the index of entropy (IoE) and the VIKOR (VIseKriterijumska Optimizacija I Kompromisno Resenje) models were applied. The fourth step aimed to introduce a novel hybrid model through the integration of the IoE and VIKOR models (IoE-VIKOR) for gully erosion susceptibility mapping (GESM). Finally, in the fifth step, a validation of the models was performed using the seed cell area index (SCAI) and the area under the curve (AUC). Remote Sens. 2019, 11, 2577 5 of 25 Figure 2. Methodology of research conducted for the study area.

Gully Erosion Inventory Map
For gully erosion susceptibility mapping (GESM), a map of the spatial distribution of the considered gullies is a prerequisite [51][52][53]. The gully erosion inventory map was prepared from the data provided by the Agricultural and Natural Resources Research Center of Semnan Province. In order to complete inventory map, the interpretation of satellite images and extensive field surveys were conducted in the study area. A total of 206 gully plots were identified in the study area ( Figure  1). In the next step, the head-cuts of the gullies were converted to points and were used to build the gully erosion susceptibility model [7,9,29,31,40,51,53]. Finally, the gullies were randomly divided into two groups, with 70% of points used for modeling (145 gullies) and 30% used for validation (61 gullies) [41,52]. One of the identified gullies in the study area is shown in Figure 3.

Gully Erosion Inventory Map
For gully erosion susceptibility mapping (GESM), a map of the spatial distribution of the considered gullies is a prerequisite [51][52][53]. The gully erosion inventory map was prepared from the data provided by the Agricultural and Natural Resources Research Center of Semnan Province. In order to complete inventory map, the interpretation of satellite images and extensive field surveys were conducted in the study area. A total of 206 gully plots were identified in the study area ( Figure 1). In the next step, the head-cuts of the gullies were converted to points and were used to build the gully erosion susceptibility model [7,9,29,31,40,51,53]. Finally, the gullies were randomly divided into two groups, with 70% of points used for modeling (145 gullies) and 30% used for validation (61 gullies) [41,52]. One of the identified gullies in the study area is shown in Figure 3.

Gully-Related Conditioning Factors (GRCFs)
In general, various GRCFs affect the occurrence of gully erosion in different regions based on different environmental conditions [18]. In the current study, our 15 GRCFs were: Elevation, slope, aspect, plan curvature, stream power index (SPI), topography wetness index (TWI), rainfall, soil type, drainage density, distance to river, distance to road, distance to fault, lithology, land use/land cover (LU/LC), and normalized difference vegetation index (NDVI) (Figure 4a-o). These factors were determined based on a literature review [41,51,52], physiographic features of the study area, the expert knowledge of natural resources experts of Semnan and Isfahan provinces, and a multicollinearity test.
Topographic parameters, namely elevation, slope degree, slope aspect, plan curvature, and TWI, and the hydrological parameters of SPI, drainage density, and distance to stream, were extracted from ALOS (the advanced land observing satellite), PALSAR (the phased array type L-band synthetic aperture radar), and a DEM (digital elevation model) with a spatial resolution of 12.5 m.
The SPI is a factor considered for the measurement of the erosive power of water runoff [30]. The SPI was obtained as per Equation (1) [54]: where As is the specific catchment area in meters and is the slope gradient in degrees. The TWI is a parameter that plays a crucial role in gully occurrence, and was calculated as per Equation (2) [54]: where S is the cumulative upslope area draining through a point and a is the slope gradient in degrees. In this research, the stream network affected gully occurrence. The distance to streams and drainage density are calculated after the extraction of the stream network from the DEM in Arc Hydro using the Euclidean distance and line density in spatial analysis tools in ArcGIS [37].
The road and fault factor maps, with scales of 1:50,000 and 1:100,000, respectively, were obtained from the National Geographic Organization of Iran (www.ngo-org.ir) and the Geological Society of Iran (GSI) (http://www.gsi.ir/), respectively.
The NDVI was computed from the LANDSAT-8 image archived by USGS (https://earthexplorer.usgs.gov/) with 30 m spatial resolution. The NDVI map of the study area was produced using Equation (3):

Gully-Related Conditioning Factors (GRCFs)
In general, various GRCFs affect the occurrence of gully erosion in different regions based on different environmental conditions [18]. In the current study, our 15 GRCFs were: Elevation, slope, aspect, plan curvature, stream power index (SPI), topography wetness index (TWI), rainfall, soil type, drainage density, distance to river, distance to road, distance to fault, lithology, land use/land cover (LU/LC), and normalized difference vegetation index (NDVI) (Figure 4a-o). These factors were determined based on a literature review [41,51,52], physiographic features of the study area, the expert knowledge of natural resources experts of Semnan and Isfahan provinces, and a multicollinearity test.
Topographic parameters, namely elevation, slope degree, slope aspect, plan curvature, and TWI, and the hydrological parameters of SPI, drainage density, and distance to stream, were extracted from ALOS (the advanced land observing satellite), PALSAR (the phased array type L-band synthetic aperture radar), and a DEM (digital elevation model) with a spatial resolution of 12.5 m.
The SPI is a factor considered for the measurement of the erosive power of water runoff [30]. The SPI was obtained as per Equation (1) [54]: where As is the specific catchment area in meters and σ is the slope gradient in degrees. The TWI is a parameter that plays a crucial role in gully occurrence, and was calculated as per Equation (2) [54]: where S is the cumulative upslope area draining through a point and a is the slope gradient in degrees. In this research, the stream network affected gully occurrence. The distance to streams and drainage density are calculated after the extraction of the stream network from the DEM in Arc Hydro using the Euclidean distance and line density in spatial analysis tools in ArcGIS [37]. The road and fault factor maps, with scales of 1:50,000 and 1:100,000, respectively, were obtained from the National Geographic Organization of Iran (www.ngo-org.ir) and the Geological Society of Iran (GSI) (http://www.gsi.ir/), respectively.
The NDVI was computed from the LANDSAT-8 image archived by USGS (https://earthexplorer. usgs.gov/) with 30 m spatial resolution. The NDVI map of the study area was produced using Equation (3): where NIR is the infrared portion of the electromagnetic spectrum (band 5) and R is the red portion of the electromagnetic spectrum (band 4).   To prepare the annual rainfall data, the rainfall statistics of Gandab, Ghahvehkhaneh-Bahar, Mahdi-Shahr, Dar-Jazin, Istegah-Pajohesh-Semnan, Kheir-Abad-Semnan, Attary, Shahmirzad, Abkhorieh-Semnan, Semnan, Chah-Kashani, Lahoord, and Sorkheh climate stations for a 30-year period (1986 to 2016) were used, and a statistical interpolation method was applied to prepare the final rainfall map. Univariate kriging was applied because it includes the simple kriging, ordinary kriging, universal kriging, and block kriging. In this study, ordinary kriging was selected to minimize the variance of error. Ordinary kriging assumes that the mean is locally constant based on the estimated point, which was adjusted to the requirements of our study area and variables.
The types of land use/land cover (LU/LC) and soil types were prepared following the information provided by the soil map of the Isfahan Agricultural and Natural Resources, Research and Education Center (http://esfahan.areeo.ac.ir/) at 1:100,000 scale, the geological map of the Geological Society of Iran (GSI) (http://www.gsi.ir/) at 1:100,000 scale, and the land use map from the Iranian Soil Conservation and Watershed Management Research Institute (https://www.scwmri.ac.ir) at 1:100,000 scale.
Finally, all layers were unified in the 12.5-meter pixel size (the same spatial resolution as the DEM) and the UTM Zone39N geographic coordinate system. An overview of the factors used for gully erosion susceptibility mapping is shown in Table 1.

Multicollinearity Test (MT)
The collinearity is a model that indicates an independent variable is the linear function of other independent variables. If the collinearity in a regression equation obtains an elevated value, it indicates a strong correlation among the independent variables, which reduces the accuracy of the model [55]. In this research, the two indicators of TOL and VIF were used to study the collinearity among independent variables. If the value of TOL is ≤ 1 and that of VIF is ≥ 10, then the existence of collinearity among the independent factors is confirmed.

Index of Entropy (IoE)
Entropy is one of the management approaches used to deal with irregularities and uncertainties in a system [56]. In order to prioritize the effective factors in the occurrence of the gully and to prepare its susceptibility map using the method mentioned above, the following equations have been used [57]: H j max = log 2 S j S j − number of classes (6) W J = I j P ij (8) where (Pij) is the probability density, Hj and Hj max represent the entropy values and maximum entropy, respectively, Ij is the information value, Sj is the number of categories, and Wj represents the resultant weight value for each factor. The range of wj is between 0 and 1. After calculating the final weight of each factor and their classes, these values were added in each related map and then weighted. Finally, they were added up, and the final gully erosion susceptibility map was prepared using Equation (9) and the weighted sum tool in ArcGIS [58].
where Wj and Pj are the final weight and the probability density for the jth feature.

VIKOR
MCDM is a common decision methodology in environmental science because it can increase the accuracy of decisions by rendering the process in a highly rational and efficient manner [59][60][61][62][63][64][65][66]. The VIKOR model is based on the adaptive planning of multi-criteria decision-making problems. The emphasis is on ranking and choosing from a set of alternatives and identifying an agreeable solution to the problem with conflicting criteria [67]. In a situation where the decision-maker is unable to identify and express the advantages of an issue at the time of its initiation and design, this method can be considered as an effective and powerful tool for decision-making. This approach focuses on the selection of a set of alternatives and sets up compromise responses for an issue with opposite criteria. Thus, it can help decision-makers to reach a final decision. The main difference between this model and other models, such as AHP, is that in these models there are no paired comparisons between criteria and alternatives, and each alternative is independently measured and evaluated by a criterion [68]. This model consists of five steps, namely: 1. Creating the decision matrix according to the number of criteria and alternatives 2. Standardizing the decision matrix using Equation (10): 3. Determining the weight of each criterion using expert knowledge 4. Determining the best and worst of the values available for each criterion (Equations (12) and (13) 5. Calculating the utility and incompatibility indexes using Equations (14) and (15), respectively: 6. Calculating the VIKOR index (Q i value) using Equation (16): where V is a constant (0.5), S * is min SS − i , is max S i , R * is min R i , and R * is max R i . Values of Q i can vary between 0 and 1. The alternatives with a score of 1 are the best, and those with a score of 0 are the worst [68].

Model Validation
The verification of the results is a key step in GESM and gully erosion susceptibility mapping [30]. The area under the receiver operating characteristic (AUROC) curve is an efficient and accurate method for verifying GESMs [32,51]. In the current research, the AUROC curve was used to evaluate the efficiency of each model in preparing the gully erosion sensitivity map [52]. This curve is an effective method for assessing the quality of models and shows the percentage of the model's prediction [69]. In the AUROC curve approach, the number of pixels that are correctly predicted by the model is plotted against the number of pixels that are incorrectly predicted. The AUROC is usually between 0.5 and 1, with values closer to 1 indicating a better performance of the method [70]. AUROC values can be classified as follows: 0.5-0.6, poor; >0.6-0.7, average; >0.7-0.8, good; >0.8-0.9, very good; and >0.9-1, excellent [71].

Applying the Index of Entropy (IoE) Model
The results of our analysis regarding the importance of the GRCFs using the IoE method revealed that the factors of soil type (1.8), LU/LC (1.3), and SPI (0.53) played the most important roles in the development of gullies in the study area. In contrast, the factors of distance to fault (0.001), aspect (0.02), and distance to road (0.02) were shown to be of least importance.
The spatial relationship between the GRCFs and the existing gullies in the study area are shown in Table 3. The results of the elevation factor analysis using the IoE method revealed that areas with low elevation are more sensitive to gully formation than the highlands. As such, the classes of <890 m (IoE = 1.24) and 890-935 (IoE = 1.42) have the highest score. These results are in line with Arabameri et al. [17]. In terms of the slope, areas with a slope between 5 • and 10 • were more prone to gully erosion (IoE = 1.56) than other areas, and, in contrast, slopes steeper than 20 • did not affect gully erosion. This is possibly because the accumulation of surface runoff and human activities were higher on the gently inclined areas [72]. With regard to the aspect factor, east aspects with IoE = 1.57 were shown to have the most significant impact on gully occurrence. This was mainly due to lower vegetation cover density in the east-facing areas compared with other aspects. The findings are in agreement with those of Arabameri et al. [43]. In terms of plan curvature, concave (IoE = 1.65) and flat (IoE = 0.98) areas had the highest impact on gully occurrence compared to the convex (IoE = 0.92) areas. The findings are in agreement with those of Arabameri et al. [43] and Rahmati et al. [72]. In the case of SPI and TWI, the classes of <9.2 and >10.44 with IoE = 1.26 and 3.78, respectively, were shown to have the most significant effect on gully occurrence. Results indicated that areas with low values of SPI and high value of TWI have a higher susceptibility to gully erosion. In areas with high values of TWI, the flow velocity, which is related to energy level of surface runoff, exceeds the soil shear stress and gully development [72].
The results related to rainfall indicated that the class of 125-215 mm with IoE = 1.18 and >0.350 with IoE = 0.17 had the highest and the lowest effect on gully occurrence, respectively. In the study area, vegetation cover is sparse in areas with lower rainfall, but mountainous areas with higher rainfall have rocky outcrops and are not prone to gully erosion. These results are in line with Arabameri et al. [73]. Based on soil types, Aridisols (IoE = 4.10) showed a strong dependency on the development of gullies in the study area. These soils are characterized by the presence of gypsum and evaporative soluble salts, which are highly sensitive to erosion [74,75].
In the case of the drainage density factor, the class of 0.14-0.18 km/km 2 with IoE = 1.111 showed a strong correlation with gully occurrence in this study area. A high drainage density causes a larger surface runoff ratio and areas with a high drainage density have a high susceptibility to gully erosion [43]. It is important to highlight the importance of the results related to the distance to rivers, faults, and roads. Several authors confirmed that the roads are one of the main causes of rill and gully development because of the destabilization of the hillslope and the elimination of the vegetation cover [26,76,77]. Our results indicated that the classes of <650 m, <5000 m, and >23,000 m with IoE = 1.64, 1.39, and 1.46, respectively, have the most impact on gully occurrence. These results are in agreement with Conoscenti et al. and Maliho et al. [28,33]. In the case of the lithology type, volcanic rocks (IoE = 1.73) and clay accumulations (1.6) showed a higher susceptibility to gully occurrence than other lithology units. During the dry and wet seasons, areas with high clay contents are subject to cracks in the ground and piping processes can develop [78,79]. Also, agricultural areas are prone to an increase in connectivity processes, which can increase runoff peaks and sediment transport [80]. According to the LU/LC factor, Kavir areas (IoE = 1.63) that cover the evaporation deposits such as gypsum and salt have a strong correlation with the location of gullies in the study area. The findings are in agreement with those of Arabameri et al. [43].
The results of the NDVI factor showed that vegetation plays a significant role in soil conservation against gully erosion, and that with increasing vegetation cover the sensitivity of the areas to this erosion decreases and that areas with the lowest vegetation cover have the highest susceptibility to gully erosion. The findings are in agreement with those of Arabameri et al. [43].

Applying the VIKOR Model
Before applying the VIKOR model, we calculated the weight of the GRCFs using pairwise comparisons and expert knowledge. According to the results of the expert knowledge method (Figure 7), soil types (0.368), lithology (0.213), and elevation (0.097) had the biggest impact on gully erosion in this study area. In a second group of factors, we can also consider with lower influence rainfall (0.088), LU/LC (0.074), slope (0.48), distance to road (0.033), plan curvature (0.023), distance to river (0.018), distance to fault (0.008), drainage density (0.006), aspect (0.004), TWI (0.001), and SPI (0.0003).
The values obtained using the VIKOR model in each sample point are shown in Table 4. The results obtained from the GESM using the VIKOR model ranged from 0.055 to 0.984. These values were divided into five classes: 0.055-0.248 (very low), 0.248-0.404 (low), 0.404-0.575 (moderate), 0.0575-0.752 (high), and 0.752-0.984 (very high) (Figure 5b). Based on the results ( Figure 6), a total of 24.5% of the study area is located in the very low susceptibility class. Furthermore, 20.2% is situated in the low susceptibility and 17.9% in moderate susceptibility class, while 19.9% was in the high susceptibility and 17.5% in the very high susceptibility class.

Applying the VIKOR Model
Before applying the VIKOR model, we calculated the weight of the GRCFs using pairwise comparisons and expert knowledge. According to the results of the expert knowledge method (Figure 7), soil types (0.368), lithology (0.213), and elevation (0.097) had the biggest impact on gully erosion in this study area. In a second group of factors, we can also consider with lower influence rainfall (0.088), LU/LC (0.074), slope (0.48), distance to road (0.033), plan curvature (0.023), distance to river (0.018), distance to fault (0.008), drainage density (0.006), aspect (0.004), TWI (0.001), and SPI (0.0003).    Table 4. The results obtained from the GESM using the VIKOR model ranged from 0.055 to 0.984. These values were divided into five classes: 0.055-0.248 (very low), 0.248-0.404 (low), 0.404-0.575 (moderate), 0.0575-0.752 (high), and 0.752-0.984 (very high) (Figure 5b). Based on the results (Figure 6), a total of 24.5% of the study area is located in the very low susceptibility class. Furthermore, 20.2% is situated in the low susceptibility and 17.9% in moderate susceptibility class, while 19.9% was in the high susceptibility and 17.5% in the very high susceptibility class.

Integration of IoE and VIKOR Models and Validation
The values obtained from the GESM ranged from 0.068 and 0.975. These values were classified into five classes to range from very low to very high based on the natural break interval method (Figure 5c). According to our results (Figure 6), 19.5%, 21%, 19.9%, 19.9%, and 19.5% of the study area had a very low, low, moderate, high, and very high susceptibility, respectively. Finally, the results of GESMs were validated using the IoE, VIKOR, and IoE-VIKOR integrated model with the AUROC curve ( Figure 8). It was shown that the IoE model (AUROC = 0.941) had a higher prediction accuracy than the other models, and that the combination of the IoE with VIKOR models (AUROC = 0.868) increased the accuracy of the VIKOR MCDM model (AUROC = 0.857). The IoE model obtained an AUROC = 0.868 prediction accuracy, while the VIKOR model obtained an AUROC = 0.857.

Integration of IoE and VIKOR Models and Validation
The values obtained from the GESM ranged from 0.068 and 0.975. These values were classified into five classes to range from very low to very high based on the natural break interval method (Figure 5c). According to our results (Figure 6), 19.5%, 21%, 19.9%, 19.9%, and 19.5% of the study area had a very low, low, moderate, high, and very high susceptibility, respectively. Finally, the results of GESMs were validated using the IoE, VIKOR, and IoE-VIKOR integrated model with the AUROC curve ( Figure 8). It was shown that the IoE model (AUROC = 0.941) had a higher prediction accuracy than the other models, and that the combination of the IoE with VIKOR models (AUROC = 0.868) increased the accuracy of the VIKOR MCDM model (AUROC = 0.857). The IoE model obtained an AUROC = 0.868 prediction accuracy, while the VIKOR model obtained an AUROC = 0.857.

Discussion
The evaluation of susceptibility to gully erosion is crucial for planning erosion control measures and mitigating the effects of this process. We demonstrated that the spatial distribution of gully erosion can be predicted by exploiting GIS analysis tools and establishing statistical relationships between gully occurrence and variability of physical variables related to topography, bedrock, soils, land cover, and rainfall. We used different modeling techniques, including statistical and data mining/machine learning, to achieve these results. In this research, a new scientific approach was developed to provide a gully erosion sensitivity map that does not require the location of gullies. This is very useful for areas that lack accurate data [81,82]. For this purpose, we confirmed that VIKOR MCDM models could be used and improved its efficiency and accuracy if the IoE bivariate model is integrated with it without comparing gullies data with each other.
The variables that we found to play a key role in our models, such as soil type, LU/LC, SPI, lithology, or elevation, are in line with other traditional and recent investigations conducted on gully erosion [27,37,40]. Zakerinejad and Maerker [27] used the maximum entropy model to investigate the sensitivity of the Maziyegan area of Fars province in Iran to gully erosion and demonstrated that many factors influence the occurrence and development of this type of erosion, including those also highlighted in our research, namely soil type, lithology, and LU/LC. Dube et al. [37] assessed gully erosion susceptibility mapping using the weight of evidence model in New Zealand and demonstrated that soil type, distance from streams, and LU/LC were the main factors affecting gully erosion. Rahmati et al. [40] compared the performance of seven state-of-the-art machine learning models, including SVM with four kernel types, BP-ANN, RF, and BRT to model the occurrence of gully erosion in the Kashkan-Poldokhtar Watershed, Iran, and stated that the drainage density, elevation, and distance to roads played a crucial role in gully occurrence. Reference [43] proposed a new model by combining the geographically weighted regression (GWR) technique with the certainty factor (CF) and random forest (RF) models to produce gully erosion zonation maps in the Mahabia watershed (Iran), and indicated that the distance to streams, distance to roads, and LU/LC factors were of primary importance in terms of gully occurrence. Our results also agree with other types of investigation carried out with field experiments where parent material, land use changes, and slope positions are key factors of soil erosion processes at smaller scales, such as hillslopes and pedons ones [41].
The validation of our results showed that IoE had a higher prediction accuracy than the VIKOR model and that the combination of IoE with VIKOE increased the prediction accuracy of the VIKOR model. These results were also confirmed by Zabihi et al. [32] and Arabameri et al. [41]. One of the main advantages of the IoE method is that this method can determine the importance of effective factors and analyze the spatial relationship between effective factors and gullies occurring in the study area [52]. The IoE model demonstrated how the most important factor could be estimated from the factors affecting the phenomenon (gully erosion) occurrence [83]. In other words, it can identify the variables that have the highest impact on the event of a phenomenon, which is important because, depending on the physiographic conditions of the area, there are usually several factors affecting gully erosion. Also, in determining the risk susceptibility using bivariate and probable statistical models, all effective factors showed the same weight, so if one factor has more impact, the effect will be ignored [84]. Therefore, the aforementioned theory as a management approach can have a significant role in the identification of effective factors and their impact on phenomena occurrence [17]. The results also confirmed that the combination of the IoE method with the VIKOR method increased the prediction accuracy of the VIKOR model. Thus, we stated that it is vital to check the effectiveness of the combination of different models rather than using single models adapted to each study area.

Conclusions
For gully erosion mapping, three different models (IoE, VIKOR, and IoE-VIKOR) were successfully applied to assess 145 gully locations and 15 environmental factors, namely elevation, slope, aspect, plan curvature, stream power index (SPI), topography wetness index (TWI), rainfall, soil type, drainage density, distance to river, distance to road, distance to fault, lithology, land use/land cover (LU/LC), and normalized difference vegetation index (NDVI). Our results showed that the IoE and expert knowledge models considered soil type, lithology, elevation, and land use to have the biggest impact on gully occurrence in the study area. The validation of our results using the AUROC curve showed that IoE has a higher prediction accuracy than the other models. The scientific methodology presented in this study can be used in areas that do not have accurate data on the spatial distribution of past and current gullies. The results of this research could be used by environmental planners to manage and reduce the risks associated with the development of gully erosion in the study area as well as in the conservation of natural resources such as water and soil.