Gully Erosion Susceptibility Assessment in the Kondoran Watershed Using Machine Learning Algorithms and the Boruta Feature Selection

: Gully erosion susceptibility mapping is an essential land management tool to reduce soil erosion damages. This study investigates gully susceptibility based on multiple diagnostic analysis, support vector machine and random forest algorithms, and also a combination of these models, namely the ensemble model. Thus, a gully susceptibility map in the Kondoran watershed of Iran was generated by applying these models on the occurrence and non-occurrence points (as the target variable) and several predictors (slope, aspect, elevation, topographic wetness index, drainage density, plan curvature, distance to streams, lithology, soil texture and land use). The Boruta algorithm was used to select the most effective variables in modeling gully erosion susceptibility. The area under the receiver operating characteristic curve (AUC), the receiver operating characteristics, and true skill statistics (TSS) were used to assess the model performance. The results indicated that the ensemble model had the best performance (AUC = 0.982, TSS = 0.93) compared to the others. The most effective factors in gully erosion susceptibility mapping of the study region were topological, anthropogenic, and geological. The methodology and variables of this study can be used in other regions to control and mitigate the gully erosion phenomenon by applying biophilic and regenerative techniques at the locations of the most inﬂuential factors.


Introduction
Gully erosion is the most destructive type of water erosion, commonly occurring in arid and semi-arid regions. Gully erosion is considered when the cross-section of the erosion channel is greater than 929 cm 2 [1]. Gullies are usually created on plain and plateau areas depending on several factors such as land use, soil type, lithology, topography and vegetation cover [2,3]. Gullies are the most destructive type of water erosion, and have many effects such as land degradation, reduced soil fertility, and damage to infrastructures [4]. Generally, gullies were observed more in semi-arid and arid regions of Iran. A study of gully data in Iran indicated that the average length and depth of gullies in Iran were about 570 m and 2.8 m, respectively, which produced 21 m 3 /m sediment per unit gully length [5]. To this end, the assessment of the gully erosion status is necessary to mitigate its destructive effects in Iran and other countries. Application of physically-based The Kondoran watershed (KW) is located along 54 • 19 to 54 • 30 E longitude and 26 • 40 to 26 • 56 N latitude in Hormozgan province, south of Iran. Its area covers 258.86 km 2 of the region, varying by 3 m to 1406 m above sea level [28]. The annual precipitation average in the study region is about 83.62 mm and the climate is hot and humid [29]. In the study area, there are several land use types including bare rock, barren lands, poor rangeland, residential lands, salt land and shrublands. The study region with an estimated population of 1628 people includes four rural regions of Armak (971 people), Kondoran (485 people), Glango (34 people) and Lavaran (138 people) [30]. Runoff and gully erosion are the main issues of the region. The soil of the region is erosion-sensitive and mostly is loamy and silty.
Soil sensitivity, heavy rains, land-use change, vegetation degradation and human interventions have caused deep gullies with depths of 1-10 m. The KW gullies contain the widest area of soil erosion in the south of Iran, which significantly has been spread toward infrastructure and settlements and has caused significant damage to roads, electrical installations, buildings, fertile lands of rural regions, and other infrastructure, and which is increasing every year. Figure 1 displays the gullies' locations in the study region in Hormozgan province, Iran.

Data Set
The present study applied data of 140 gullies in the Kondoran watershed, collected through a field study using a GPS (Global Positioning System). Furthermore, the information of 100 gully absence places was gathered. Of these, 70% of data were randomly selected for training the learning models, and the rest (30%) were used to test and validate the fitted models [31].

Predictor Variables of Gully Erosion Susceptibility
Gullies are created under various conditions. Therefore, different variables influence gully formation and development. Selecting appropriate variables is one of the critical and effective issues in modeling gully erosion susceptibility. Various researchers have considered different variables as influence factors on gully erosion. Previous studies have considered 11 factors for gully erosion susceptibility modeling [32,33]. Figure 3 presents the methodological diagram of the present study. Using the Boruta algorithm, the most critical variables were selected for modeling gully erosion susceptibility. The evidential belief function (EBF) algorithm was applied to describe the correlation among predictor variables and gully erosion status. The gully susceptibility map in the Kondoran watershed of Hormozgan was generated by applying the multiple diagnostic analysis (MDA), support vector machine (SVM), and random forest (RF) algorithms. The models were applied to the occurrence and non-occurrence points. The AUC (area-underthe-curve) and ROC (receiver operating characteristics) have been used to assess model performance. The ensemble model was used to combine these models and produce the ensemble gully susceptibility map of the region.

Slope Angle
Slope angle is an essential factor that affects gully erosion. Gully erosion occurs when macropores become larger by infiltration of water into the soil [34]. The larger macropores result in a greater risk of soil collapse and gully erosion. Gentle slopes have a higher penetration rate than steeper slopes. Thus, gentle slopes are more prone to gully erosion. In this study, the slope angle map (degree) was obtained based on the digital elevation model (DEM) (Figure 4a).

Slope Aspect
The slope aspect is another factor that influences solar radiation reaching the ground, soil moisture, and vegetation cover [35]. Thus, it is considered an effective factor in mapping gully erosion susceptibility. The study region's slope aspect map was generated using DEM 30 m and classified into nine classes (Figure 4b).

Elevation
Elevation affects the microclimate and vegetation cover. This factor affects the runoff production process as one of the main stimuli of gully erosion. Therefore, this study considered elevation as one predictor for gully erosion susceptibility modeling [36]. In the present study, a DEM of 30 m was used, extracted from a topographic map on a scale of 1:50,000 ( Figure 4c).

Soil Texture
Soil texture and soil erodibility significantly affect gully erosions [37]. The study region's soil texture map was provided at a 1:50,000 scale using the database of the Agricultural Research Center of Hormozgan Province [38] In this study, the soil texture is classified into four classes of loamy skeletal, coarse loamy, fine silty and fine loamy ( Figure 4d).

Length-Slope Factor
The length-slope factor determines the influence of topography on erosion; therefore, it is considered one of the critical factors in predicting gully erosion status [39]. This factor is computed as below:  (1) in which the FAG is the flow accumulation grid and 0.01745 is a coefficient that converts the radian to degree measure. The LS index was obtained using SAGA-GIS 7.3.6 software based on 30 m DEM (Figure 4e).

Plan Curvature
Plan curvature is defined as the curvature of contour lines in the horizontal plane [40]. This factor affects converging and diverging surface flow, flow rate, and soil erosion. Thus, plan curvature can be applied as an important predictor of gully erosion. The plan curvature map was generated using DEM in an ArcGIS 10.7 environment. The plan curvature map of the study area was categorized into three classes, including concave, convex and flat categories (Figure 4f). The positive and negative values of plan curvature indicate the convexity and concavity of curvature slope, and the values close to zero determine that the surface is flat [1].

Topographic Wetness Index
The topographic wetness index (TWI) is obtained by combining the upstream catchment area and slope, and is commonly used to quantify the topography's influence on hydrological processes. It also determines the effect of topography on the level of saturation surface to produce runoff [41]. TWI calculates the probability of flow accumulation in soil due to the upstream catchment area and slope [42]. Thus, TWI was applied as a topographical variable for modeling gully erosion. The TW index is obtained using the following equation: where As is the specific area of the upstream watershed (m 2 /m) and b is the slope angle in degree scale. In this study, the TWI was generated in SAGA-GIS software based on the 30 m DEM (Figure 4g).

Land Use
Land use is another factor that affects gully erosion and land degradation [13]. The land use map of the study region was provided at a 1:50,000 scale from the Agricultural Research Center of the Hormozgan Province [43]. According to Figure 4h, barren land (without vegetation), poor rangeland, bare rock, salt land, residential lands and shrublands are the main land use types of the study area. The map of the predictor variables was generated in raster format and then classified using the existing function GIS software. Table 1 indicates information on predictor variables that affect gully erosion susceptibility. The lithology factor has an important role in the gully erosion phenomenon since it depends on lithological properties and weathering of materials at the ground surface [1]. In this study, the lithological map of the study region at a 1:50,000 scale was provided from the Geological Survey of Iran [38].

Drainage Density
Drainage density is defined as the entire length of the channel per unit of watershed area (km/km 2 ) [44]. Regions with higher drainage density contain lower infiltration and higher runoff. Higher drainage density facilitates the evacuation of sediments from upstream areas in watersheds, causing expanding gully erosion. Thus, drainage density can be considered a predictor variable of gully erosion hazard [36]. The study area's drainage map was created using the line density function in an ArcGIS 10.7 environment (Figure 4j).

Distance to River
Distance to the river is the most important factor that affects the gully erosion process. The drainage network affects the stability and slope saturation degree, and consequently, the soil erosion process [45]. Common surface materials near the rivers contain higher moisture, which results in an acceleration of the processes of shedding, separation and soil transfer. Soil saturation depends on proximity to the stream and therefore, areas closer to streams are more susceptible to gully erosion [32]. The distance to the stream map of the study region was provided using the Euclidean distance method in ArcGIS 10.7 software. Finally, the map was classified into five categories of less than 276.5 m, 276.5 to 611.8 m, 611.8 to 1055.13 m, 1055.13 to 1806.48 m, and 1806.48 to 3118.41 m, based on the natural break approach (Figure 4k). The map of the predictor variables was generated in raster format and then classified using existing functions in GIS software. Table 1 indicates the information of predictor variables that affect gully erosion susceptibility.

Collinearity between Independent Variables
Existing collinearity between variables reduces the model's accuracy and performance [33]. The tolerance coefficient and variance inflation factor (VIF) are two important indices for investigating dependency between variables [46]. Variables often have collinearity problems when the tolerance coefficient or VIF between them are respectively less than 0.1 or greater than 5 [7].

Boruta Variable Selection Algorithm
Feature selection is an important step in classification and regression modeling. Feature selection algorithms are used to remove the unimportant predictors to improve the model accuracy [47,48]. The Boruta algorithm is a feature selection algorithm that is placed under the RF classification method [49]. It uses shadow features which are copies of original features. The shadow features are randomly assigned to objects; therefore, decision trees are generated based on the shadow features.
The importance of an attribute is measured by the loss of accuracy of the classification model caused by attributes that are randomly assigned to objects. Then, the mean and the standard deviation of the accuracy loss are calculated, and the Z-scores are obtained by dividing the mean and its standard deviation. In the Boruta algorithm, the Z-score is considered as the measure of importance. Therefore, the set of the importance of shadow attributes is used as a reference for detecting the importance of original attributes. Then, the importance of original features is compared with the highest importance of shadow features [27]. The Boruta algorithm has the following steps: (1) The information system is extended by generating the shadow attributes (at least five for each attribute). The attributes with importance significantly lower than MZS are removed (considered as unimportant), and those with importance significantly higher than MZS are considered as important. (6) All shadow attributes are removed and the procedure is repeated until the importance is assigned to all attributes.
The Boruta algorithm has advantages in comparison with other feature selection algorithms [27]. First, this algorithm follows an all-relevant variable selection approach in which it considers all features that are related to the output variable; whereas, most of the other variable selection algorithms follow a minimum optimum approach where they depend on a small subset of features which causes a minimum error on a selected classifier. Second, it considers multi-variable relationships and also can investigate interactions between variables. In this study, the Boruta package in R software was used to determine the importance of gully erosion susceptibility factors.

Investigating the Relationship between Gully Erosion and Conditioning Factors
This study applies the evidential belief function (EBF) algorithm to describe the correlation among predictor variables and gully erosion status. The evidential belief function (EBF) is a robust approach for obtaining reliable models through incorporating various factors to reduce uncertainty [50]. The statistical EBF model is computed based on Dempster-Shafer theory to combine the representations of several independent variables to achieve a combined measure of belief.
The degree of uncertainty (Unc), degree of disbelief (Dis), degree of belief (Bel), and degree of plausibility (Pls) are the main parameters of the EBF approach, and each describes particular information of a dataset [51]. The Bel and Pls are respectively defined as the lower and upper bound of probabilities, and the uncertainty degree (Unc) is described as the difference between Bel and Pls, which shows the ignorance or doubt that the evidence supports a hypothesis. Additionally, Dis is a degree of disbelief in evidence concerning the hypothesis, and the Dis value is obtained using (1-Pls) or (1-Unc-Bel) equations.

Modeling Gully Erosion Susceptibility
The support vector machine (SVM) is a supervised learning approach that is used for classification or regression modeling. This method is defined based on statistical learning theory and uses the structural risk minimization (SRM) method to obtain an optimized solution [52]. The SVM applies linear or non-linear (polynomial or radial) kernels to learn a classification model. This study applied the radial basis (RBF) kernel due to the high performance of this function. In this study, the SVM method was run using R software 3.5.3 and the SDM (Species Distribution Modelling) package [53].

Random Forest Model
Random forest (RF) is a classification approach that is obtained based on the improvement of bagging (bootstrap aggregation) trees [54]. The bagging tree is a decision tree that is built on bootstrap samples [55]. The main difference between random forest and bagging is the number of predictors in each bootstrapping step. On the other hand, for p predictors, the bagging trees use m = p, while the random forest splits m ≈ √ p predictors in each bootstrapping tree. Finally, the RF model is defined as the average of bootstrapping trees. The RF procedure is used to reduce the variance of the statistical learning method. Robustness to outliers, limiting overfitting and errors are further advantages of the RF approach [56].

Multiple Discriminant Analysis
Multiple discriminant analysis is a classification approach that is used to predict categorical responses. This model, also known as the Fisher discriminant analysis, is defined based on Bayes' theorem. The MDA attempts to estimate the conditional probability and the predictors are assumed to follow a multivariate normal distribution. The linear (LDA) and quadratic (QDA) discriminant analysis are special cases of MDA used for binary responses when respectively, a linear or non-linear boundary between classes is assumed. The normality assumption of the predictor variables is a limitation of using the MDA method. For the normal case, the conditional probability leads to the linear combination for which the coefficients should be estimated [57].

Ensemble Model
The ensemble model (EM) is a combination of several learning models and is commonly used to improve classification algorithms by decreasing variance and bias or improving prediction accuracy [58]. Incorporating the single models' predictions is done using weighted and unweighted averaging. In this study, the ensemble model was assembled using weighted averaging based on AUC statistics Equation (3): where EM is the ensemble model, and AUC i is the AUC value of the ith single model (Mi).
In the present study, the SVM, RF, MDA and ensemble methods have been run using R software 3.5.3 and the SDM (Species Distribution Modelling) package [53]. The probability maps were provided based on all models and were classified using the natural break method in ArcGIS 10.7. The classification map obtained from the EM is different from the maps obtained from the other three models. This is because the EM is the weighted average of SVM, RF and MDA models using the AUC weighting and the unclassified maps of every single model are merged to create an ensemble model. To merge the unclassified output map of single models, the weighted average is used.

Evaluation Model Performance
This study applied AUC (area under the receiver operating characteristic curve) and True Skill Statistics (TSS) to evaluate model performance. The AUC statistics describe the area under the receiving operator characteristics curve (ROC) and are used as a measure of classification accuracy. The greater AUC indicates a better classification result [58]. Each observation belongs to a positive (gully existence) or negative (non-existence) category in gully classification. The AUC determines the classification accuracy based on the ROC. The number of positive and negative pixels correctly classified is called true positive (TP) or true negative (TN). In contrast, it can be defined as the false positive (FP) and false-negative (FN) for wrongly classified pixels [59]. In the ROC curve, the X-axes represent the specificity and Y-axes indicate the sensitivity and are defined as below: The AUC varies between 0 and 1, such that a value close to 0 indicates that the learning model is not able to classify the observations, while an AUC close to 1 shows that the observations are classified very well [60]. In addition, an AUC less than 0.6, between 0.6 to 0.8, and greater than 0.8, describes low, moderate and high classification accuracy, respectively [33].
The true skill statistic (TSS), also known as Hanssen-Kuipers discriminant, is another common measure for evaluating classification accuracy [61]. TSS statistics explain the ratio of true predictions (true positive and true negative) [62]. TSS value changes between −1 to +1, where +1 determines high classification accuracy and values below zero show a performance no better than random [62]. The TSS is computed as TSS = Sensitivity + Specificity−1.

Multi-Collinearity Test
Multi-collinearity is a serious problem in regression and classification, and refers to a situation in which two or more predictor variables are linearly correlated. In this study, we applied variance inflation factor (VIF) and tolerance (TOL) criteria to explore collinearity between explanatory variables. According to these criteria, a TOL less than 0.1 and a VIF above 10 indicate that there exists a collinearity problem among predictor variables. Table 2 provided VIF and TOL values between gully erosion predictors. As shown, there exists no collinearity between predictor variables.

Gully Erosion Susceptibility Maps (GESM)
The maps of gully erosion susceptibility were generated using effective parameters based on the SVM, RF and MAD models. Next, a new map was provided using the ensemble model, which was defined based on the RF, MAD and SVM models. The maps were classified into low, moderate, high, and very high susceptibility classes using the natural break classification method ( Figure 5). The results revealed that most parts of the study region are allocated in the very high susceptibility category. Further analysis of Table 3 indicated that the MDA model allocated the lowest part (429 km 2 , 1.67%) and the highest part (169.93 km 2 , 66.03%) of the study regions respectively assigned to the moderate and very high susceptibility, while according to the SVM model, the lowest part (31.06 km 2 , 12.07%) and the highest part (121.52 km 2 , 47.22%) of the study regions belong to the low and very high susceptibility classes. In this regard, based on the RF and ensemble models, the low susceptibility areas respectively covered 73.37 km 2 (28.51%) and 74.34 km 2 (28.88%) of the study regions, while very high susceptibility areas respectively covered 84.07 km 2 (32.67%) and 97.69 km 2 (37.96%) of the study region. According to the results, in all models, the study area's central and eastern parts belong to the very high susceptibility class and the northern parts of the region are in the low susceptibility class.

Evaluation of GESMs Performance
The performance of the fitted models for the gully erosion susceptibility maps was evaluated using area-under-the-curve (AUC) and true skill statistics (TSS). The results are reported in Table 4 and Figure 6. All models present high values in both AUC and TSS. The ensemble model had the highest values in both criteria. Thus, the accuracy among the fitted models in predicting gully erosion susceptibility is high.   6. Accuracy assessment of gully erosion susceptibility based on testing data (30%).

Computing Variable Importance Using Boruta Algorithm
Results and outputs of the Boruta algorithm (Table 5) indicate that distance to the stream (33.5), land use (17.41), elevation (12.18), and lithology (7.34) are the most important gully erosion factors, followed by soil type, drainage density, LS factor, and slope angle. In addition, the plan curvature, TWI, and slope aspects have the least rank among the predictor variables and should be ignored from the modeling process.

Evaluating the Relationship among Conditioning Factors and Gully Erosion Using the Evidential Belief Function (EBF) Model
The relationship between gully erosion and predictor variables was investigated and results are shown in Table 6. According to the results, the highest weight of distance to streams variable (0.51) has been assigned to class 1, [0-276.5], and is regarded as the most effective factor on gully erosion occurrence in the study area. For this factor, higher distances (>1055.3 m) have zero weight in Bel measure. In this regard, the class [1.8-3.69] km/km 2 of drainage density variable had the highest Bel (0.45) weight and is considered as the area with a high probability of gully occurrence. The highest correlation (Bel = 0.44) of gully erosion with TWI variable belongs to the [12-22.3] class, which denotes a positive correlation between these two variables.  In the case of elevation factor, the complete correlation (Bel = 1) between gully erosion and  m category denoted that the gully erosion generally occurs in lowlands and there is a low probability for the gully erosion occurrence in higher elevation lands. Additionally, the results of the relationship between gully erosion and slope length factor indicated that the highest Bel value (0.88) was observed at the [0-2.5] class, signifying that commonly, the gully erosion in the region occurs with low slope length (<5.91). In the case of slope degree, the highest Bel value was observed between gully erosion and the [<10 • ] classes. Further, the slope aspect results determined that most of the gullies have been formed (Bel = 0.25) in the north-east directions in the study area.
In the case of the lithology variable, the highest correlation (Bel = 0.71) was obtained between gully erosion and Qfp, followed respectively by Qal and Qaf classes. In other lithological units, the Bel value was zero. The plan curvature is another predictor variable that defines the curvature of a contour line in the horizontal plane and determines a region's position [45,63]. As the result indicated, the [−0.1 −0.1] curvature class had the highest impact (Bel = 1) on gully erosion occurrence, compared to the other curvature classes.
Concerning the land-use variable, most gully erosions have occurred in the barren lands (0.71), followed by shrublands (0.29). Finally, the soil type analysis indicated that the highest correlation (Bel = 0.94) was observed between gully erosion and the fine silty soil texture category, followed by the coarse loamy class (Bel = 0.06).

Discussion
Gully erosion constitutes a serious problem for land degradation in a wide range of environments. The Kondoran watershed is one of the most susceptible regions to gullying in the Hormozgan province due to environmental factors and human interventions. Indeed, this study provides an understandable platform to inform decision-makers about the current situation of this study area. To this end, we evaluated gully erosion susceptibility in the study area based on multiple diagnostic analysis (MDA), support vector machine, and random forest (RF) algorithms, and also a combination of these models, namely the ensemble model.
Each of the models contains specific weaknesses in modeling gully susceptibility mapping, which can be resolved using the ensemble modeling method. The ensemble method is a procedure that combines a specific number of learning models in such a way that its predicted response is the weighted average of the predicted response obtained from the incorporated models [26]. The evaluation of the performance of models showed that the ensemble model had the highest accuracy in predicting gully erosion susceptibility. These results are consistent with the findings of Pourghasemi et al. [64] in using the ensemble model for predicting gully erosion susceptibility.
The outputs of fitted models were spatially different because it depends on the different structures of the models, and also the importance of input variables to each model. For instance, the MDA model is highly correlated with altitude, the SVM model has the highest correlation with LS factor, while the RF and the ensemble models showed the highest correlation with the distance to the stream predictors. The Boruta feature selection algorithm was applied to determine the most important gully erosion factors. Results and outputs of the Boruta algorithm indicated that distance to the stream was the most influential factor on gully erosion susceptibility. This result corroborates the findings of Conoscenti et al. [45], and Rahmati et al. [61].
The relationship between gully erosion and predictor variables was investigated using the EBF model. The results showed that the highest weight of distance to streams variable had been assigned to the regions with the nearest distance to streams. Therefore, most gullies in study regions have been formed close to the streams because the soils near streams are water-saturated, and as a result, the disintegration and collapse of soil particles intensifies; subsequently, gully erosion susceptibility increases [64]. The regions with the most drainage density had a high probability of gully occurrence due to lower infiltration and more runoff events [31]. The highest correlation of gully erosion with the TWI variable belongs to the regions with the higher topographic wetness. Generally, the higher topographic wetness results in a higher probability of occurrence of gully erosion, corroborating with Pourghasemi et al. [64] and Hembram et al. [65].
The study of elevation factor indicated that gully erosion generally occurs in lowlands and there is a low probability for the gully erosion occurrence in higher elevation lands. The probability of washing and disintegrating the soil particles increases in lowland regions because the runoff accumulation is higher than in elevated regions. As a result, the gully erosion hazard increases in lowlands, in line with the findings of Dickson et al. [66] and Amiri et al. [7]. According to the results, the regions with lower slope angles and slope length are more prone to gully erosion than the regions with higher slope angles and slope length, and this is similar to the findings of Azareh et al. [32]. When the slope increases, the water velocity increases; subsequently, the soil erosion rate increases. However, gully erosion generally forms on low slopes and flat areas. In other words, to form a gully erosion, water flow should accumulate at a point and create a hole; a phenomena which commonly happens to low slopes.
Investigating the lithology variable indicated that the highest gully erosion has occurred in Quaternary floodplain deposits (Qfp). These sediments contain fine-grained floodplain deposits, silt, and fine-grained sand and clay. Therefore, these regions are the most susceptible to gully erosion. In the case of the land-use variable, most gully erosions have occurred in the barren lands that, due to the lack of vegetation, are more susceptible to gully erosion. This is consistent with the finding of Fernández-Raga et al. [62]. Soil type analysis showed that most gully erosions have occurred in the fine silty category, followed by the coarse loamy class. These soils are more susceptible to erosion because their particles are easily separated and transferred after moisture absorption, and thereby cause the gully erosion phenomenon [67,68].
Raising the level of public awareness about this phenomenon, not interfering and misappropriating in nature, preventing the destruction of natural resources, preventing land-use change and runoff management in watersheds by implementing watershed management projects, especially rainwater harvesting systems and vegetation restoration, can be considered as some of the most important preventive measures in controlling gully erosion and land degradation in the Hormozgan province [69].
The study of gully erosion susceptibility maps indicated that, in all models, the central and eastern parts of the studied region belong to a very high sensitivity class, while the northern part of the region is located in the low sensitivity class. Most of the gullies are located in lowlands, near the streams, in silty soils, in barren lands and shrublands, and in proximity to the rural settlements. In the past few years, the gully erosion was expanded toward the Kondoran and Armak villages lands in the middle part of the study region. Therefore, experts and planners should prioritize the management of these lands to reduce damages caused by gully erosion. Distance to streams and elevation cannot change; meanwhile, restoration and development of the degraded land might change the land use and decrease the susceptibility to gully erosion [70,71]. Additionally, sustainable soil management could be achieved by afforestation of the barren lands and regenerative agriculture techniques in areas most susceptible to gully erosion [72,73].

Conclusions
In recent years, gully erosion has caused losses of a large volumes of soil in the Kondoran watershed in Hormozgan province, Iran. Therefore, accurate prediction of gully erosion susceptibility is a fundamental issue to protect the soil and reduce the erosion rate. This study applied four classification models, including random forest (RF), support vector machine (SVM), multivariate discriminant analysis (MDA), and the ensemble model (as the combination of the three previous models) to provide gully erosion maps of the study region. The accuracy of gully erosion susceptibility maps was evaluated using the area under ROC curve (AUC) and TSS metrics. Results and outputs indicated that the ensemble method (AUC = 98.2%, TSS = 0.93) has better performance than the RF (AUC = 97.1%, TSS = 0.91), SVM (AUC = 93.2%, TSS = 0.84) and MDA (AUC = 91.4%, TSS = 0.82) models.
The Boruta feature selection algorithm was used to determine the most important factors in gully erosion susceptibility mapping. It was found that gully occurrence in the study area was mainly influenced by topological factors (elevation, and distance from streams), anthropogenic factors (land use and proximity to rural settlements), and geological factors (lithology and soil texture). Further, the results of the Boruta algorithm denoted that the plan curvature, TWI, and slope aspect could be ignored from the modeling process. The results indicated that the central, eastern, and southern parts of the study area are more susceptible to gully erosion. Models and predictors applied in this study can be used in similar areas. Regenerative agriculture techniques and afforestation of these lands should be prioritized to reduce land degradation in the present and future.
Appendix A   Table A1. The abbreviations included in the text are reported alphabetically.

AGNPS
Agricultural Non-point Source AGNPS is a computer-simulation model that simulates the behavior of runoff, sediment, and nutrient transport from watersheds that have agriculture as their prime use. The model operates on a cell basis and is a distributed parameter, event-based model.

AUC of ROC
Area-under-the-curve of Receiver Operating Characteristics The AUC statistics describe the area under the ROC curve and are used as a measure of classification accuracy. The greater AUC indicates a better classification result.

Bel Belief
Belief is one of the functions that is used in EBF model. It is the pessimistic measures of the spatial relationship of response variable, i.e., Bel indicates the lower probabilities of evidence that supports a hypothesis.

DEM Digital Elevation Model
A digital elevation model (DEM) is a representation of the bare ground (bare earth) topographic surface of the Earth excluding trees, mountains, buildings, and any other surface objects.

Dis Disbelief
Disbelief is one of the functions that is used in EBF model. It is a degree of disbelief in evidence for the hypothesis and the Dis value is obtained from 1-Pls or 1-Unc-Bel.

EBF Evidential Belief Function
The evidential belief function (EBF) algorithm describes the correlation among predictor variables and response variable. The statistical EBF model is computed based on Dempster-Shafer's theory to combine the representations of several independent variables to achieve a combined measure of belief.

GIS
Geographic Information System A geographic information system (GIS) is a conceptualized framework that provides the ability to capture and analyze spatial and geographic data.

GLM Generalized Linear Model
GLM is the extension of the classic linear regression model. Contrasted with the normal linear model, the response variables of GLM are not confined to normal distribution, and these response variables can also obey binomial or Poisson distributions. In addition, the link function is introduced into GLM to establish the relationship between the expectation of the response variable and the linear combination of explanatory variables.

GPS
Global Positioning System GPS (Global Positioning System) is a radio wave receiver used to provide coordinates that give the exact position of an element in a certain space.

LiDAR Light Detection and Ranging
LiDAR is a method for determining ranges (variable distance) by targeting an object with a laser and measuring the time for the reflected light to return to the receiver.

LISEM Limburg Soil Erosion Model
The Limburg soil erosion model (LISEM) is a physically-based hydrological and soil erosion model which can be used for planning and conservation purposes.

MaxEnt Maximum Entropy
MaxEnt is a data mining method to predict the occurrence of one event based on maximum entropy that approximates the probability distribution of presence data based on environmental limitations.

MDA Multiple Diagnostic Analysis
Multiple discriminant analysis is a classification approach that is used to predict categorical responses. This model, also known as the Fisher discriminant analysis, is defined based on Bayes' theorem. The MDA attempts to estimate the conditional probability and the predictors are assumed to follow a multivariate normal distribution.

Pls Plausibility
Plausibility is one of the functions that is used in EBF model. It is the optimistic measures of the spatial relationship of response variable, i.e., Pls indicates the upper probabilities of evidence that supports a hypothesis.

RF Random Forest
Random forest (RF) is a classification approach that is obtained based on the improvement of bagging (bootstrap aggregation) trees.

RUSLE Revised Universal Soil Loss Equation
RUSLE is an easily and widely used model that estimates rates of soil erosion caused by rainfall and associated overland flow.

SAR Synthetic Aperture Radar
Synthetic aperture radar (SAR) is a form of radar that is used to create two-dimensional images or three-dimensional reconstructions of objects, such as landscapes. SAR uses the motion of the radar antenna over a target region to provide finer spatial resolution than conventional stationary beam-scanning radars.

SVM Support Vector Machine
The support vector machine (SVM) is a supervised learning approach that is used for classification or regression modeling. This method is defined based on statistical learning theory and uses the structural risk minimization (SRM) method to obtain an optimized solution.

TOL Tolerance
Tolerance is relevant and frequently used quantities that may be consulted to examine individual predictors for potentially strong contributions to (near) multicollinearity. This index reflects estimates of the degree of interrelationship of an independent variable with other explanatory variables in a regression model. The TOL less than 0.1 indicates that there exists a collinearity problem among predictor variables.

TSS True Skill Statistic
The true skill statistic (TSS) is known as Hanssen-Kuipers discriminant, and is commonly measure for evaluating classification accuracy. The true skill statistics is defined based on the components of the standard confusion matrix representing matches and mismatches between observations and predictions.

Topographic Wetness Index
The topographic wetness index (TWI) is a physically-based index of the effect of local topography on runoff flow direction and accumulation. The index is a function of both the slope and the upstream contributing area. The computation of TWI is performed using both geographic information systems (GIS) and Python, a programing software used to enhance computing capabilities. The indices help identify rainfall runoff patterns, areas of potential increased soil moisture, and ponding areas.

Unc Uncertainty
Unc value is one of the functions that is used in EBF model. It is the difference between the Pls and Bel function, which shows the ignorance or doubt that the evidence supports a hypothesis.

VIF Variance Inflation Factor
Variance inflation factor measures how much the behavior (variance) of an independent variable is influenced, or inflated, by its interaction/correlation with the other independent variables. Variance inflation factors allow a quick measure of how much a variable is contributing to the standard error in the regression. VIF is the reciprocal of Tolerance. The VIF above 10 indicates that there exists a collinearity problem among predictor variables.