Gully Head-Cut Distribution Modeling Using Machine Learning Methods—A Case Study of N.W. Iran

: To more e ﬀ ectively prevent and manage the scourge of gully erosion in arid and semi-arid regions, we present a novel-ensemble intelligence approach—bagging-based alternating decision-tree classiﬁer (bagging-ADTree)—and use it to model a landscape’s susceptibility to gully erosion based on 18 gully-erosion conditioning factors. The model’s goodness-of-ﬁt and prediction performance are compared to three other machine learning algorithms (single alternating decision tree, rotational-forest-based alternating decision tree (RF-ADTree), and benchmark logistic regression). To achieve this, a gully-erosion inventory was created for the study area, the Chah Mousi watershed, Iran by combining archival records containing reports of gully erosion, remotely sensed data from Google Earth, and geolocated sites of gully head-cuts gathered in a ﬁeld survey. A total of 119 gully head-cuts were identiﬁed and mapped. To train the models’ analysis and prediction capabilities, 83 head-cuts (70% of the total) and the corresponding measures of the conditioning factors were input into each model. The results from the models were validated using the data pertaining to the remaining 36 gully locations (30%). Next, the frequency ratio is used to identify which conditioning-factor classes have the strongest correlation with gully erosion. Using random-forest modeling, the relative importance of each of the conditioning factors was determined. Based on the random-forest results, the top eight factors in this study area are distance-to-road, drainage density, distance-to-stream, LU / LC, annual precipitation, topographic wetness index, NDVI, and elevation. Finally, based on goodness-of-ﬁt and AUROC of the success rate curve (SRC) and prediction rate curve (PRC), the results indicate that the bagging-ADTree ensemble model had the best performance, with SRC (0.964) and PRC (0.978). RF-ADTree (SRC = 0.952 and PRC = 0.971), ADTree (SRC = 0.926 and PRC = 0.965), and LR (SRC = 0.867 and PRC = 0.870) were the subsequent best performers. The results also indicate that bagging and RF, as meta-classiﬁers, improved the performance of the ADTree model as a base classiﬁer. The bagging-ADTree model’s results indicate that 24.28% of the study area is classiﬁed as having high and very high susceptibility to gully erosion. The new ensemble model accurately identiﬁed the areas that are susceptible to gully erosion based on the past patterns of formation, but it also provides highly accurate predictions of future gully development. The novel ensemble method introduced in this research is recommended for use to evaluate the patterns of gullying in arid and semi-arid environments and can e ﬀ ectively identify the most salient conditioning factors that promote the development and expansion of gullies in erosion-susceptible environments.


Introduction
Gullies are common features in arid and semi-arid regions, and they are major causes of sediment erosion; they supply from 10 to 94% of the total sediment yield in some watersheds [1]. High erosion rates undercut agricultural sustainability and necessitate the search for (usually expensive) solutions in the context of costly governmental policies. However, studying and predicting gully erosion is difficult [2][3][4]. In terms of the ecosystem effects and environmental damages from gully erosion, studies have focused on the influential factors and on identification of susceptible areas using geographic information systems (GIS) and remote sensing (RS) [5][6][7][8]. This study develops a new model to detect and predict gully locations with high spatial accuracy to reduce gully erosion damages.
One of the difficulties in the regional GESM process is that the factors influencing gully erosion require data usually derived from various sources at different spatial scales, which may contain uncertainties and imprecisions. Traditional data-driven approaches cannot be used to determine the relationships between geo-environmental factors and gully erosion occurrence because of the limitations caused by imbedded statistical assumptions about variables' independence and data distributions in susceptibility analyses [20,21]. New modeling methods are needed that go beyond traditional data-driven approaches, and methods that can deal with the above issues and can enhance model performance.
Ensemble models have been used in GESM due to their novelty and their ability to comprehensively assess gully-erosion parameters for discrete classes of independent factors [51,52]. Although some studies have been conducted on the spatial prediction of gullies, a standard framework considering all influential factors for achieving a reasonable and reliable prediction has not been established. Some studies and techniques should be used in different hydro-geomorphological environments to devise a global framework for gully-erosion modeling. Additionally, some factors contribute to gullying that are either difficult to recognize (and measure), or they are difficult to convert to raster formats for modeling. Therefore, one of the future fields of gully modeling should focus on the detection and Water 2020, 12, 16 3 of 26 application of the unknown factors that influence gully formation. This may be achieved by combining gullying research with GIS and data-mining tools to create a tool or technique that can map future, unknown factors. This could help planners, decision makers, and environmental managers to prepare gully erosion maps of the highest quality with the best possible accuracy to better manage gullying in erosion-susceptible areas.
The main difference between this study and previously conducted studies is that this study explores a new ensemble-intelligence approach that employs bagging as a meta-classifier with an alternating decision tree (ADTree) as a base classifier to spatially predict gully erosion. The results produced by this new ensemble-intelligence approach are compared to the results generated with a single alternating decision tree, a rotational-forest-based alternating decision tree (RF-ADTree), and benchmark logistic regression (LR) to assess and improve the accuracy of GESM. These ML modes haven't been used for GESM, so we assess the performance of the new ensemble model using a variety of statistical metrics and the area under the curve (AUC).
The Chah Mousi watershed (northeastern Iran) is an arid region very prone to gully erosion. Gullies are widespread throughout the region and cause land degradation and economic damages every year. This study illustrates and compares individual and ensemble machine learning models to assess gullying susceptibility. We test the efficacy of these models and compare them to find the most suitable model for land use planning. The main objectives of this study are identifying and mapping the extant gullies in the Chah Mousi watershed by (a) creating an inventory; (b) mapping, modeling, and predicting the locations of gully head-cuts; (c) characterizing the roles of various geo-environmental features as factors that control the distribution of gullies; and (d) evaluating gully erosion susceptibility in the study area.

Study Area
The Chah Mousi watershed is in Semnan province, Iran, and is located between 35 • 15 05" and 35 • 37 12" N and 54 • 35 44" and 55 • 23 05" E ( Figure 1). It is a relatively small area of approximately 2176.02 km 2 . The greatest change in elevation is along a NE to SE axis. The average elevation in the northeastern quadrant is 2123 m.a.s.l. In the southeastern quadrant, it is 672 m.a.s.l. As the region is relatively small, the slope degree varies significantly from flat to 67.8 • , although the average is about 3 • . Due to the predominance of flat landscapes, standing and slow-moving water is more typical than runoff. The mean annual precipitation ranges from 48 to 206 mm, principally during the wet season from January to March [59]. Temperatures typically reach a peak of 41 • C during summer, especially in the south, and a low below 0 • C during winter in the northern parts of the region; though average temperatures during the rest of the year range from 13 to 23 • C [59]. Together, these numbers indicate the potential for meteorological stress on the land surface with high thermal and precipitation variations and local spikes that may cause freezing and thawing of soils and expansion and contraction within the regolith [13].
The land covers include agriculture, bare land, kavir (barren sandy and rocky desert), rangeland, rock outcrops, salt lakes, wetlands, and salt lands. The latter are particularly vulnerable to dissolution processes during the wet season as the salt crust is easily weathered, giving rise to pores that promote changing groundwater levels and erosion of soils [60]. The distribution of salt crusts is evident in the regional soil map (primarily in areas featuring aridisols and entisols and where the outcropping lithologies are also reported). The main lithological units in the study area are marl, gypsiferous marl and limestone, shale, sandstone, granite, conglomerate, and salt flat [60].

Gully Mapping
Archival records containing reports of gully erosion that have been compiled by the Semnan Agricultural and Natural Resources Research and Education Center were used as the first source of locational data. Upon this historical foundation, gully locations and dimensions were identified and measured using remotely sensed data viewed through Google Earth. Finally, a field survey was conducted in the study region to update and refine the inventory ( Figure 2). Sites of gully head-cuts were geolocated with a DGPS (Differential Global Positioning System) device. The survey yielded 119 gully head-cuts ( Figure 2) to be used for modeling. Of the overall dataset, 75 gullies (63.02%) were identified from archives, 19 gullies (15.96%) collected using Google Earth, and 25 gullies (21.008%) were collected in a field survey. All gullies were checked and mapped using DGPS with millimeter accuracy. The universal transverse Mercator (UTM) coordinate system was used. The models described above were applied to the locations of 83 head-cuts (70% of the total). The models were tested (or validated) with the remaining 36 gully locations (30% of the total). As the models selected in this study correspond to a family that predicts the presence or absence of a phenomenon, an equal number of locations (36 no gully locations as validation data and 83 no gully locations as calibration data) were selected and tested as well [52]. In turn, this procedure creates a balanced dataset for the subsequent analyses, although it should be noted that the geomorphological features still debates whether balanced or unbalanced datasets should be created prior to a susceptibility analysis [19,58]. Some of mapped gullies are shown in figure 3.

Gully Mapping
Archival records containing reports of gully erosion that have been compiled by the Semnan Agricultural and Natural Resources Research and Education Center were used as the first source of locational data. Upon this historical foundation, gully locations and dimensions were identified and measured using remotely sensed data viewed through Google Earth. Finally, a field survey was conducted in the study region to update and refine the inventory ( Figure 2). Sites of gully head-cuts were geolocated with a DGPS (Differential Global Positioning System) device. The survey yielded 119 gully head-cuts ( Figure 2) to be used for modeling. Of the overall dataset, 75 gullies (63.02%) were identified from archives, 19 gullies (15.96%) collected using Google Earth, and 25 gullies (21.008%) were collected in a field survey. All gullies were checked and mapped using DGPS with millimeter accuracy. The universal transverse Mercator (UTM) coordinate system was used. The models described above were applied to the locations of 83 head-cuts (70% of the total). The models were tested (or validated) with the remaining 36 gully locations (30% of the total). As the models selected in this study correspond to a family that predicts the presence or absence of a phenomenon, an equal number of locations (36 no gully locations as validation data and 83 no gully locations as calibration data) were selected and tested as well [52]. In turn, this procedure creates a balanced dataset for the subsequent analyses, although it should be noted that the geomorphological features still debates whether balanced or unbalanced datasets should be created prior to a susceptibility analysis [19,58]. Some of mapped gullies are shown in Figure 3.

Gully Erosion Conditioning Factors
Several factors affect a location's susceptibility to gully erosion [17,19]. After completing a study of the gully-erosion literature, and considering local conditions and data availability, 18 variables were selected for inclusion in the modeling process. These include elements of topographical, geological, and hydrological conditions.
The following topographical factors were considered: elevation, slope gradient, aspect, plan curvature, convergence index (CI), slope length (LS), topographic wetness index (TWI), topographic position index (TPI), and terrain ruggedness index (TRI). Each was calculated using PALSAR DEM with 12.5 m spatial resolution applying the basic terrain analyses in SAGA GIS. A detailed explanation of the equations used to calculate LS, TWI, and SPI is available in Arabameri et al. [19].
The description of the lithology was acquired from a geological map at a scale of 1:100,000 (Geological Survey Department of Iran, [59]). The map was digitized and 6 geological classes were identified in the study area: A (including marl, gypsiferous marl, and limestone; dacitic to andesitic volcano sediment; well-bedded green tuff and tuffaceous shale; dacitic to andesitic volcanic; dacitic to andesitic volcano breccia; andesitic volcano breccia, sandstone, marl, and limestone; granite, pale-red polygenic conglomerate, and sandstone), B (including phyllite, slate, and meta-sandstone; Jurassic dacite to andesite lava flows), C (including Cretaceous rocks, in general), D (including light red to brown marl and gypsiferous marl with sandstone intercalations; red marl, gypsiferous marl, sandstone, and conglomerate), E (including fluvial conglomerate, piedmont conglomerate, and sandstone), and F (salt flat, high-level piedmont fan and valley terrace deposits, low-level piedmont fan and valley terrace deposits, and salt lake) ( Figure 4p).

Gully Erosion Conditioning Factors
Several factors affect a location's susceptibility to gully erosion [17,19]. After completing a study of the gully-erosion literature, and considering local conditions and data availability, 18 variables were selected for inclusion in the modeling process. These include elements of topographical, geological, and hydrological conditions.
The following topographical factors were considered: elevation, slope gradient, aspect, plan curvature, convergence index (CI), slope length (LS), topographic wetness index (TWI), topographic position index (TPI), and terrain ruggedness index (TRI). Each was calculated using PALSAR DEM with 12.5 m spatial resolution applying the basic terrain analyses in SAGA GIS. A detailed explanation of the equations used to calculate LS, TWI, and SPI is available in Arabameri et al. [19].
The description of the lithology was acquired from a geological map at a scale of 1:100,000 (Geological Survey Department of Iran, [59]). The map was digitized and 6 geological classes were identified in the study area: A (including marl, gypsiferous marl, and limestone; dacitic to andesitic volcano sediment; well-bedded green tuff and tuffaceous shale; dacitic to andesitic volcanic; dacitic to andesitic volcano breccia; andesitic volcano breccia, sandstone, marl, and limestone; granite, pale-red polygenic conglomerate, and sandstone), B (including phyllite, slate, and meta-sandstone; Jurassic dacite to andesite lava flows), C (including Cretaceous rocks, in general), D (including light red to brown marl and gypsiferous marl with sandstone intercalations; red marl, gypsiferous marl, sandstone, and conglomerate), E (including fluvial conglomerate, piedmont conglomerate, and sandstone), and F (salt flat, high-level piedmont fan and valley terrace deposits, low-level piedmont fan and valley terrace deposits, and salt lake) ( Figure 4p).
The hydrological gully erosion factors that were included in the modeling process are drainage density, distance-to-stream, mean annual rainfall, and stream-power index (SPI). Drainage density and distance-to-stream were calculated using the stream network information developed from the PALSAR DEM in ArcGIS 10.5. Raster maps of these factors were prepared using line-density and Euclidean-distance tools in ArcGIS 10.5. The SPI was calculated as follows: where As is the specific catchment area, and β is slope ( • ). Annual precipitation data were obtained for the period from 1984 to 2014 recorded at the Toroud, Razveh, Moalleman, and Hosseinan weather stations operated by the Iran Meteorological Organization (IRIMO, 2014). The rainfall data were interpolated using the kriging interpolation tool in ArcGIS 10.5. Gully erosion is also influenced by soils, land use, and vegetation [19]. Therefore, these factors are represented by soil types, land use/land cover (LU/LC), and normalized difference vegetation index (NDVI), and were used as conditioning factors. Soil type data were based on the information from the Soil Conservation Section of Agricultural and Natural Resources Research Centre of Semnan Province. LU/LC and NDVI data were obtained from Landsat 8 images (15 August 2017) with a 30 m resolution. The LU/LC map containing eight classes (agriculture, bare land, kavir, poor range, rock, salt lake, salt land, and wetland) was prepared using the supervised classification method and maximum likelihood in ENVI4.8 software. The map was verified using the kappa coefficient with 459 ground control points (GCP). The kappa value of the resulting map was 0.976. The NDVI was calculated using Landsat 8 bands 4 (red) and 5 (infrared) data in ArcGIS 10.5.

Rotational Forest (RF)
RF modeling is a relatively new ensemble algorithm that increases the accuracy and diversity of base classifiers, and it was first proposed by Rodriguez et al. [39]. The success of RF modeling depends on the rotation matrix generated by transformations and base classifiers [61,62]. The basis of RF modeling is principal component analysis (PCA), which can extract features and create training datasets for learning base classifiers [63]. RF has been applied to classification problems, such as landslide-susceptibility research, land use mapping, and flash flood susceptibility research [64][65][66]. (1) Ri is produced by the following four steps: (i) Divide B into K subsets, and the number of gully conditioning factors of each subset is Q = n/K.
(ii) In case of the classifier Ai, let Bi,j be the jth, where j = 1, 2, 3,… and K is the subset of gully conditioning factors. Xi,j is the gully conditioning factor of Bi,j from X. Bi,j is randomly selected from the Xi,j with the 75% size by bootstrap algorithm. Then, Xi,j' would be transformed to achieve coefficient ai,1 (1) , ai,1 (2) , …, ai,l (Qi) , the size of ai,1′ is Q × 1. (iii) Arrange a sparse rotation matrix Ri with the obtained coefficients. (iv) The confidence of each class is calculated by the average combination method in the given test sample χ, Roads also affect gully erosion as they intercept and concentrate overland flow [17]. This factor is represented by the distance to road in gully and non-gully locations, which is determined by vectorizing topographic maps and Google Earth images, and then transforming the data to a raster map using line density tools in ArcGIS 10.5.

Rotational Forest (RF)
RF modeling is a relatively new ensemble algorithm that increases the accuracy and diversity of base classifiers, and it was first proposed by Rodriguez et al. [39]. The success of RF modeling depends on the rotation matrix generated by transformations and base classifiers [61,62]. The basis of RF modeling is principal component analysis (PCA), which can extract features and create training datasets for learning base classifiers [63]. RF has been applied to classification problems, such as landslide-susceptibility research, land use mapping, and flash flood susceptibility research [64][65][66]. Suppose x = (x 1 , x 2 , x 3 , . . . , x n ) is the vector of the landslide conditioning factor, y = (y 1 , y 2 ) is the vector of landslide or non-landslide class, X is the training dataset, A 1 , A 2 , A 3 , . . . , A L are the classifiers in the ensemble, and B is the landslide conditioning factor set. The steps of training classifier A i are as follows. The rotation matrix R i a generated by the matrix of R i is shown in Equation (2).
R i is produced by the following four steps: (i) Divide B into K subsets, and the number of gully conditioning factors of each subset is Q = n/K.
(ii) In case of the classifier A i , let B i,j be the jth, where j = 1, 2, 3, . . . and K is the subset of gully conditioning factors. X i,j is the gully conditioning factor of B i,j from X. B i,j is randomly selected from the X i,j with the 75% size by bootstrap algorithm. Then, X i,j ' would be transformed to achieve coefficient a i,1 (1) , a i,1 (2) , . . . , a i,l (Qi) , the size of a i,1 is Q × 1.
(iii) Arrange a sparse rotation matrix R i with the obtained coefficients.
Water 2020, 12, 16 9 of 26 (iv) The confidence of each class is calculated by the average combination method in the given test sample χ, where γ i,k (ηR i a ) is the probability produced by the classifier A i to the hypothesis that η belongs to the class k.

Alternating Decision Tree
The alternating decision tree (ADTree) model is an ensemble model that consists of a boosting algorithm and a decision tree [67]. It is a generalization of a decision tree in which each node is replaced by a splitter node and a prediction node [68,69]. The base rule mapping from an instance to real number involves a precondition c 1 , a base condition c 2 , and two real numbers a and b. If c 1 ∩ c 2 , the prediction is a, and the prediction is b when c 1 ∩ − c 2 ; − means negation. The values of a and b are determined by Equations (4) and (5), respectively.
where W(p) is the total weight of training instance. The best c 1 and c 2 values are obtained by minimizing the Z t (c 1 , c 2 ), which is defined as Equation (6).
Suppose that R is a set of base rules. Then, a new rule can be defined as R t+1 = R t + r t , r t (x), which shows two prediction values (a and b) at every layer of the tree. x is a set of instances. The classification of instances is the sign of the sum of all predicted values in R t+1 : The algorithm first finds the best constant prediction for the whole data set [70]. Cross validation is often used for selection [71].

Bagging
Bootstrap aggregation or bagging (BAG) was introduced by Breiman in 1996 [72]. The bootstrap technique randomly selects and replaces samples to generate multiple samples to form a training dataset. Every subset generated is used to build a decision tree, and they are later aggregated in the final model. The accuracy of classification is improved by reducing the variance of classification error [73,74]. In recent years, BAG has been widely applied in landslide susceptibility research and has performed well [75][76][77].

Logistic Regression
Logistic regression (LR) is one of the most popular multivariate statistical analysis methods [78][79][80]. It can make a multivariate regression correlation between a dependent variable and several independent variables [81,82]. The advantage of LR is that the variables can be continuous, discontinuous, or a combination of the two [83,84]. In this study, the main purpose of using an LR model is to determine the relationships between landslide occurrence and gully conditioning factors, calculated using Equation (8).
where P is the probability of gully occurrence and ranges from 0 to 1. Z is a linear sum of constants, and its range is (−∞, +∞). The calculation equation of Z can be defined as Equation (9).

Frequency Ratio
The ratio between the frequency of occurrences and non-occurrences at a location within a given causative factor class is called the FR [19]. Larger ratios suggest that those factor classes are more important determinants of the occurrence (in this case, gully-erosion proneness or susceptibility. As there are numerous pertinent factors at play in each location (or area defined by a pixel in our digital map), the potential for gully erosion can be computed as the sum of all ratios for the predisposing factor classes [19]. FR is empirical. It is, in fact, not a statistical method; it is not based on statistical distributions.

Random Forest (RAF)
RAF uses multiple trees to classify locations based on a single conditioning factor [85]. The RAF algorithm continuously replaces the factors affecting each pixel space, thereby creating numerous decision trees. A combination of all decision trees in a study area provide the information to support decision making [85]. An RAF contains 3 user-defined parameters: (1) the number of variables used to construct each decision tree, which indicates the power of each independent tree; (2) the number of trees included in the RF; and (3) the minimum number of nodes within the trees. The prediction power of RAFs increase as the strength of independent trees increases and as the correlation between them decreases. Sixty-six percent of the data (the testing data) are used to grow a tree, and the result is called a bootstrap. A randomly introduced predictor variable splits a node in the tree's construction during the growing process. The remaining third of the data is used to evaluate (or validate) the fitted tree. The average of all predicted values produced during several iterations of the algorithm creates the final modeled prediction. In this model, two factors-the mean decrease accuracy and the mean decrease Gini index-are used to prioritize the effective factors. Comparing the mean decrease accuracy to the mean decrease Gini index determines the relative importance of the effective factors, especially the relationships between environmental factors. RAF analyses were carried out in R 3.3.1 using the "Randomforest" package [85].

Multicollinearity Assessment
In GESM, testing for collinearity among the effective factors in gullying is very important, because the collinearity reduces the accuracy of the GESM [86][87][88][89]. The variance inflation factor (VIF) and Tolerance (TOL) are very commonly used indicators for checking multicollinearity among parameters [90,91]. TOL values less than 0.1 or 0.2 and VIF values greater than 5 or 10 indicate collinearity between the parameters [17,19,86,89,92]. In the present study, the multicollinearity test of gully erosion conditioning factors (GECFs) was done using Equations (10) and (11) in SPSS software: where R 2 J is the regression coefficient for determining independent variable j.

Methodology
As described above, an inventory of gullies was created, and the gully-erosion conditioning data were compiled in a GIS to provide input for the modeling process ( Figure 5). The gully sites were divided into two datasets: 70% were used for training, and 30% were used for validation of the models. An assessment of multicollinearity among the conditioning factors was performed. The relative weights of the GECFs were determined using an RAF model, and an analysis of the spatial relationships between GECFs and gullies was conducted with FR. GESMs were created using each of the four models: ADTree, RF-ADTree, Bagging-ADTree, and LR. Finally, the models were evaluated and validated using the receiver operating characteristic (ROC) curves and by calculating the area under the ROC curve (AUC) for each model [93][94][95]. The AUC values are between 0 and 1, which can be interpreted following these categories: 0.6-0.7 have poor, 0.6-0.7 medium, 0.7-0.8 good, 0.8-0.9 very good, and 0.9-1 excellent accuracy [9,17,19]. The four models used were objectively compared to determine the most effective approach.

Multicollinearity Assessment
A multicollinearity analysis of the GECFs was performed ( Table 1). The analysis revealed that TOL and VIF values for all factors are >0.1 and <5, respectively, indicating that the variables are not significantly correlated and that they can be used in further analyses.

Spatial Relationship between Gully Locations and Conditioning Factors by Applying FR Model
Analyses of the spatial relationships between gully locations and GECFs (Table 2) showed that classes of conditioning factors with FR values greater than 1 are susceptible to gully erosion [17].

Multicollinearity Assessment
A multicollinearity analysis of the GECFs was performed ( Table 1). The analysis revealed that TOL and VIF values for all factors are >0.1 and <5, respectively, indicating that the variables are not significantly correlated and that they can be used in further analyses.

Spatial Relationship between Gully Locations and Conditioning Factors by Applying FR Model
Analyses of the spatial relationships between gully locations and GECFs (Table 2) showed that classes of conditioning factors with FR values greater than 1 are susceptible to gully erosion [17].  According to LS factor, areas with the lowest slope length have the highest susceptibility to gully occurrence, so that class of <15.2 m, with FR = 1.244, showed the strongest relationship to gullying in the study area.
Generally, TPI values > 0 indicate ridges, 0-flat areas (or constant slopes), and <0-valleys. This is confirmed with the statistical relationships between gully locations and TPI values in the study area. Most of the study area is flat and classes of TPI < 1.96 are those with the gully locations. This is in accordance with TRI values that show terrain heterogeneity. Higher TRI values show increased local relief heterogeneity. In contrast, lower TRI values indicate more level surfaces (e.g., planar surfaces or various depositional landforms). The results showed that gullies occur in areas belonging to classes of TRI values < 7.84, and the most susceptible are areas with TRI < 1.47. Despite the occurrence of gullies, the terrain is quite homogenous; most of the study area is flat. TWI reveals the areas with drainage depressions where water is likely to accumulate. Thus, the areas with high values of TWI should be more susceptible to gully formation, which is in accordance with the results that showed that higher TWI values (>11.8) have a higher occurrence of gullies in the study area. SPI values indicate potential flow-erosion at a point in the topographic surface. Most of the gullies occur in areas where SPI values are <14.9 (FR = 4.66).
Distance-to-stream and drainage density are important factors conditioning gully occurrence [17]. Gullies occur mainly in the areas close to streams (<100 m) [13]. In addition, most of the gullies occur in areas receiving 68 to 85 mm of precipitation annually [16] (Table 2).
In lithological units, class of B (phyllite, slate and meta-sandstone, and Jurassic dacite to andesite lava flows) showed the strongest correlation with gully occurrence in the study area.
According to NDVI, class of 0.043 to 0.132 had the highest FR (1.34) and therefore the strongest relationship to gully formation. Moreover, most of the gullies occur in areas of kavir and poor rangeland, which had FR values of 1.961 and 0.672, respectively. Gully erosion occurs mainly in areas with entisols/aridisols (Table 3).
Roads may intercept overland flow and promote gully formation. Most of the gullies occur near roads (<1000 m) [16]; the strongest relationship is <500 m (FR = 6.43).

Validation of Results
The results were validated using AUC values both in SRC (success rate curve) and PRC (prediction rate curve) (Table 4, Figure 8a,b). Generally, the models tested achieved excellent accuracy. The success rate curves, a degree-of-fit measure (i.e., comparison of the susceptibility maps with training dataset), indicated that bagging-ADTree (0.964) was most accurate, and LR least accurate (0.867). The AUC values computed for prediction rate curves, indicating the predictive power of the susceptibility maps, confirmed that Bagging-ADTree was most accurate (0.978) and LR least (0.870).

Validation of Results
The results were validated using AUC values both in SRC (success rate curve) and PRC (prediction rate curve) (Table 4, Figure 8a,b). Generally, the models tested achieved excellent accuracy. The success rate curves, a degree-of-fit measure (i.e., comparison of the susceptibility maps with training dataset), indicated that bagging-ADTree (0.964) was most accurate, and LR least accurate (0.867). The AUC values computed for prediction rate curves, indicating the predictive power of the susceptibility maps, confirmed that Bagging-ADTree was most accurate (0.978) and LR least (0.870). (prediction rate curve) (Table 4, Figure 8a,b). Generally, the models tested achieved excellent accuracy. The success rate curves, a degree-of-fit measure (i.e., comparison of the susceptibility maps with training dataset), indicated that bagging-ADTree (0.964) was most accurate, and LR least accurate (0.867). The AUC values computed for prediction rate curves, indicating the predictive power of the susceptibility maps, confirmed that Bagging-ADTree was most accurate (0.978) and LR least (0.870).

Discussion
Different sources were used to prepare the input dataset. Because many factors used in GESM were extracted from a digital elevation model (DEM), the quality of the DEM significantly influences the accuracy of the results [96,97]. The Advanced Land Observing Satellite (ALOS) DEM with 12.5 m spatial resolution was used as it has been shown to provide better accuracy than both the Shuttle Radar Topography Mission (SRTM) and Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) and DEMs [98].
In this study, we developed and explored a new ensemble intelligence approach using bagging and RF as a meta-classifier and with ADTree as a base classifier, to spatially predict gully head-cut erosion in the Chah Mousi watershed. We produced GESMs based on a modeling procedure including training and validation datasets, and 18 conditioning factors (elevation, slope angle, aspect, plan curvature, CI, LS, SPI, TPI, TRI, TWI, distance to stream, drainage density, rainfall, distance to road, NDVI, lithology, land use/land cover, and soil type). These factors were checked for collinearity with statistical metrics, including TOL and VIF. The results reveal that all GECFs influenced gully erosion occurrence.
Based on FR analysis, the relationship between the factors and gully locations were assessed. Conditioning-factor classes with FR values >1 indicated areas with greater gully-erosion susceptibility [82]. Elevation plays an important role in vegetation and precipitation type and, therefore, controls the spatial distribution and gully erosion processes [99]. Elevations in the study region below 1000 m a.s.l. are more susceptible to gully erosion. Thus, the higher occurrence of gully head cut erosion in the lowland areas agrees with Dickson et al. [100]. However, Arabamiri et al. [19] determined that elevations below 829 m were most prone to gullying. In terms of slope angle and curvature, the FR analysis showed that slopes of less than 5 • (including flat areas) were most likely to be sites of gully occurrence. Because lower slope angles have greater soil depth, intensive rainfall impaction and greater runoff from upslope will decrease soil strength resulting in the development and extension of the gully channel [9]. Curvature causes accumulation of runoff and enhances the velocity and volume of flow, so this variable positively correlates to locations of gully erosion. The slope aspect that controls several climate conditions, such as the intensity of precipitation, moisture, evapotranspiration, and vegetation cover [101], indirectly influences gully erosion. Among the slope-aspect classes, east-and southeast-facing slopes are the most highly correlated to gully erosion. These two slope aspect classes get more solar radiation in the northern hemisphere and, as a result, they experience more evaporation, higher soil porosity (total pore space), lower soil strength, and lower vegetation density. This is in accordance with Zabihi et al. [9], who reported that southward slope aspects are more susceptible to gully erosion. CI values below −39.6 100/m were most predictive of gully formation: the lower the CI value, the greater is the potential for gully erosion. Arabameri et al. [17] concluded, based on the WoE method, that CI values between 0 and 10 signify locations that are more susceptible to gully occurrence in their study area. LS less than 15 m indicate a more likely formation of gullies and reflects that gullies are more likely formed in flat areas with lower slope angles. This confirms the findings of Gayen et al. [102], but conflicts with the results of Zabihi et al. [9], who shows a direct relationship between LS and gully erosion locations. Zabihi et al. also implied that the higher the LS, the higher the probability of gully erosion occurrence due to increasing runoff velocity and a decreasing detachment and transport threshold of soil particles [103,104].
The most susceptible classes for the other GECFs were SPI between >14.9, TPI less than −7, TRI less than 1.4, and TWI more than 11.8. These results are confirmed by the findings of Arabameri et al. [17] who reported that, for example, the greater the TWI factor, the greater is the potential for gully occurrence. High values of TWI increase the filtration rate and provide the conditions for piping and roof collapse, resulting in the development of gully tunnels and, eventually, the appearance of gullies on the surface [105].
Moreover, the nearer sites are to a river, the higher the susceptibility to gully erosion. In this study, locations at distances less than 100 m from a stream were more likely to see gully formation. Some researchers have confirmed these results [9,13,16,42]. The sheer force of flow can overcome and decrease the strength of soil along the sides of gully forms and lead to the development of gullies of greater dimensions.
Areas with drainage densities exceeding 1.75 km/km 2 were most correlated to gully erosion. The role of this factor can be made clearer when other factors are considered. For example, a location with a lower slope angle and higher drainage density has a higher TWI, and if the soil at that location was loose and erodible, gully erosion is easier to achieve. In the study area, the lower classes of annual precipitation amounts (between 68.3 and 85.7 mm) were most susceptible to gully incidence. This suggests that though rainfall has a positive role in gully formation, it is not the most important factor. In other words, lower rainfall values are positively related to gullying.
Distances from roads are important to gully erosion and, like distances from rivers, the nearer the site, the higher the potential for gully erosion. Distances of less than 500 m from a road were positively correlated to gully locations, which underscores the importance of the roles of development and disturbance of ground surfaces in promoting landscape degradation.
Results of the NDVI factor show that vegetation plays a very important role in protecting soil against erosion, so that, with increasing vegetation, the sensitivity of an area to gully erosion decreases. Vegetation cover greatly reduces the erosion of runoff through the increase in infiltration and protection of soil through roots [106]. The findings agree with those of Arabameri et al. [13], Arabameri et al. [19], and Chaplot et al. [107] stating that low values of NDVI have a positive association with gully erosion and that it is easier for a gully to develop in areas with lower NDVI values. Generally, barren lands and sparsely vegetated areas are more susceptible to erosion than forests, where vegetation cover strongly reduces the erosive action of surface runoff.
Because gully erosion depends on the lithological properties of materials at Earth's surface, lithology is a vital factor in gullying [104]. As for lithology, Quaternary lithotypes have a high susceptibility to gully erosion. The result is in agreement with findings reported by Arabameri et al. [13], who found that Quaternary lithotypes have a strong effect on gully occurrences. In terms of land use, which plays a key role in geomorphological and hydrological processes by controlling overland flow runoff generation and sediment dynamics [108], the areas of kavir are most susceptible to gully erosion. In these regions, the complete lack of vegetation leaves the soil exposed, and it is easily eroded by precipitation. These results are in line with [13]. The entisol/aridisol soils are the most susceptible soils to gully erosion occurring in the study area, which is in accordance with [19].
In terms of the FR values, the most important GECFs in the study area were the distance to nearest road and drainage density. This is confirmed by the RAF algorithm analysis, which was used to rank the importance of the GECFs for the spatial prediction of gullies in the study area. This result is consistent with [17,109,110]. Roads are impervious surfaces, and they disrupt natural drainage systems due to improper culverts, concentration of surface runoffs, and by altering the hydrological functions of hillslopes, which significantly contribute to overland flow and allow rapid run-off, easily eroding bare soil and causing gullying [111,112]. An example of the effect of roads on gullying is shown in Figure 9. Distance to a road is the most important factor. It is followed in importance by drainage density, distance to stream, land use, rainfall, NDVI, elevation, SPI, TPI, CI, lithology, soil type, plan curvature, TRI, aspect, and LS. Though other factors affect gully erosion, the above are the most meaningful in the study area.
gully erosion based on the past patterns of formation, but it also provides excellent predictions of future development. The RF and bagging as a meta-classifier can decrease over-fitting and noise problems in the training dataset. Some researchers have confirmed the prediction power of RF in applications to some environmental problems [42,[115][116][117].
For example, Tien Bui et al. [21] predicted gully locations in a semi-arid watershed of Iran using ADTtree and its ensembles using RF meta-classifier. They concluded that the RF model could well enhance the prediction power of ADTree as a base classifier. However, the RF-ADTree ensemble model outperformed some benchmark models, including SVM based on the polynomial and RBF kernels, LR, naïve Bayes, and ADTtree. Additionally, Shirzadi et al. [42] used four meta-classifiers, namely, multiboost, bagging, RF, and random subspace (RS), for the spatial prediction of shallow landslides in Bijar City, Kurdistan province, Iran. They used ADTree as a base classifier for the modeling process. The four ensemble models were combined with the ADTree under two scenarios of different sample sizes and raster resolutions. They reported that the RS model was more capable for sample sizes of 60%/40% and 70%/30% with a raster resolution of 10 m. According to the results, the new proposed ensemble model can spatially predict gully erosion occurrences with reasonably good accuracy.

Conclusions
Soil erosion is an important environmental challenge to ecosystem's condition and function. Land degradation and decreasing land productivity are a result of on-site and off-site erosion in a gully-prone area. However, detection, prediction, and management of gully-prone areas using protective measures and mitigation techniques are important efforts. Some quantitative and qualitative methods and techniques have been developed and explored for modeling and preparing A novel ensemble intelligence approach, bagging-ADTree, and other ML algorithms-ADTree, RF-ADTree and LR-were used to create gully erosion susceptibility maps. The goodness-of-fit and the performance of the models were checked by AUROC of success and prediction rate curves. The results illustrate that bagging ADTree and RF-ADTree outperformed ADTree and LR. These results are in line with [42,113,114]. The new model accurately identified the areas that are susceptible to gully erosion based on the past patterns of formation, but it also provides excellent predictions of future development. The RF and bagging as a meta-classifier can decrease over-fitting and noise problems in the training dataset. Some researchers have confirmed the prediction power of RF in applications to some environmental problems [42,[115][116][117].
For example, Tien Bui et al. [21] predicted gully locations in a semi-arid watershed of Iran using ADTtree and its ensembles using RF meta-classifier. They concluded that the RF model could well enhance the prediction power of ADTree as a base classifier. However, the RF-ADTree ensemble model outperformed some benchmark models, including SVM based on the polynomial and RBF kernels, LR, naïve Bayes, and ADTtree. Additionally, Shirzadi et al. [42] used four meta-classifiers, namely, multiboost, bagging, RF, and random subspace (RS), for the spatial prediction of shallow landslides in Bijar City, Kurdistan province, Iran. They used ADTree as a base classifier for the modeling process. The four ensemble models were combined with the ADTree under two scenarios of different sample sizes and raster resolutions. They reported that the RS model was more capable for sample sizes of 60%/40% and 70%/30% with a raster resolution of 10 m. According to the results, the new proposed ensemble model can spatially predict gully erosion occurrences with reasonably good accuracy.

Conclusions
Soil erosion is an important environmental challenge to ecosystem's condition and function. Land degradation and decreasing land productivity are a result of on-site and off-site erosion in a gully-prone area. However, detection, prediction, and management of gully-prone areas using protective measures and mitigation techniques are important efforts. Some quantitative and qualitative methods and techniques have been developed and explored for modeling and preparing the susceptibility assessments. However, due to differences in their probability distribution functions, their performances are also different. For example, some of them do not fit the data that are available. All models present advantages and disadvantages, so one of the most important aspects of the modeling strategy is selecting the appropriate model. Machine-learning models are more often used because of their ability to overcome over-fitting and noise challenges during the modeling process and because they have higher goodness-of-fits and perform better compared to other conventional models. Moreover, among the machine-learning classifiers, ensemble models are more powerful than single classifiers. They randomly divide a training dataset into subsets and perform a single classifier, which provides an output with the lowest error and the highest performance rather than the single classifier. This process overcomes the weakness of the single classifier and achieves a more powerful classifier. In response to the advantage of ensemble classifiers, a novel ensemble intelligence approach, namely bagging-ADTree, was performed and gully erosion maps were obtained. Some other machine-learning algorithms (including ADTree, Bagging-ADTree, and LR) were used for comparison and validation of the results of the new model. The random forest model is used to determine the relative importance of conditioning factors. The results indicate that distance-to-road and drainage density are very important to gully occurrence in the study area. The validation indicated that although the models achieved high goodness-of-fit scores and were powerfully predictive, the ensemble model was better than others at spatially predicting gully erosion and produced a more accurate gully-susceptibility map of the study area. Based on these results, we can recommend the new model, bagging-ADTree, for gully modeling in other zones of potential gully erosion susceptibility, but offer one caution: there may be other conditioning factors responsible for gully erosion in other areas. Finally, the results from a case study of the Chah Mousi watershed show that selecting suitable predisposing factors and combining machine-learning ensemble models with GISs can be used to efficiently predict an area's susceptibility to gully formation with high accuracy. Therefore, the gully-erosion susceptibility map generated by the method can aid decision makers, planners, and engineers in their quests to identify and develop the most effective protective measures to sustainably prevent and mitigate gully-erosion damage.