**Alireza Arabameri 1,\*, Artemi Cerda <sup>2</sup> and John P. Tiefenbacher <sup>3</sup>**


Received: 24 April 2019; Accepted: 24 May 2019; Published: 29 May 2019

**Abstract:** Gully erosion is an environmental problem in arid and semi-arid areas. Gullies threaten the soil and water resources and cause off- and on-site problems. In this research, a new hybrid model combines the index-of-entropy (IoE) model with the weight-of-evidence (WoE) model. Remote sensing and GIS techniques are used to map gully-erosion susceptibility in the watershed of the Bastam district of Semnan Province in northern Iran. The performance of the hybrid model is assessed by comparing the results with from models that use only IoE or WoE. Three hundred and three gullies were mapped in the study area and were randomly classified into two groups for training (70% or 212 gullies) and validation (30% or 91 gullies). Eighteen topographical, hydrological, geological, and environmental conditioning factors were considered in the modeling process. Prediction-rate curves (PRCs) and success-rate curves (SRCs) were used for validation. Results from the IoE model indicate that drainage density, slope, and rainfall factors are the most important factors promoting gullying in the study area. Validation results indicate that the ensemble model performed better than either the IoE or WoE models. The hybrid model predicted that 38.02 percent of the study area has either high or very high susceptible to gullying. Given the high accuracy of the novel hybrid model, this scientific methodology may be very useful for land use management decisions and for land use planning in gully-prone regions. Our research contributes to achieve Land Degradation Neutrality as will help to design remediation programs to control non-sustainable soil erosion rates.

**Keywords:** soil; natural resources; modeling; hybrid model; Bastam watershed

## **1. Introduction**

Soil erosion in arid and semi-arid regions is one of the important factors that should be taken into consideration in land use planning [1–3]. Soil erosion is a major consequence of environmental and ecological change [4–6]. Water soil erosion has long been of interest to soil conservation researchers [1]. Topography, types of landforms, and resistance to erosion are the main determinants of erosion type and severity, and this is especially the case with linear erosion processes such as gullies [7].

Among the types of erosion, gully erosion most threatens numerous environmental resources and land use sustainability. This threat is not limited to soil degradation, changes in the landscape and/or landcover, limitations of agricultural activities, or the economic exploitation of natural resources. Gully erosion also promotes the initiation and expansion of badlands, promotes floods, lowering of water tables, desertification, and production and transportation of significant volumes of sediments along the watersheds to the coastal and lowlands [8–10].

A gully is a drainage channel with steeply sloping sidewalls and an active, steep eroding head-cut caused by fluctuating surface flows (i.e., during or after severe rainfall) [11]. Most scholars regard the destruction of natural ecosystems, non-sustainable land uses, degradation of vegetative cover, overgrazing, climate change, geological conditions, and human disturbance of natural systems to be the main causes of gullying [12–19].

Some studies have found that the range of sediment production attributable to gully erosion is 10% to 94% in different ecosystems and that spatial and temporal factors, soil types, land uses, topography, climate, and others can also influence sediment production [11]. A study of 22 watersheds in Spain determined that annual sediment production in areas without gullying is 0.74 Mg/ha, but in similar environments where gullying is occurring the annual production is 2.97 Mg/ha [20]. A single gully can generate sediments annually at a rate of up to 93,750 Mg/km [21].

Because the erosion is so fast at the head of a gully, there are no measures that can be taken to control or stop head-cut development; the measures that are taken are likely to only slow down their growth [22]. To prevent the rapid growth of gullies or to minimize the damage they cause, it is vital to understand the morphology of a gully, the manner in which it forms, and the causes for its development. Studying the amount of development of gully head-cut erosion in a specific timeframe can enable the prediction of the rate of expansion and growth, and this can also improve damage estimates. Identification of areas that are more susceptible to gully erosion is therefore possible. Accurate and adequate information about the type and extent of headward head-cut erosion can enable better land use decisions and more effective environmental management [23].

For more than a half-century, prediction and modeling of soil erosion have been a valuable component of soil conservation and engineering planning [24]. Mathematical simulations can be used to estimate erosion that is caused by an array of independent factors. Geomorphologic models are empirically based and must not include theoretical assumptions of relationships. They also must avoid scenarios and should not ignore multicollinearity among conditioning factors for the phenomenon [25]. Without a clear explanation of a geomorphological conceptual model, it is impossible to identify the appropriate sub-models that can numerically express the phenomenon [25]. The solution is to define the modeling goals and to identify the unknowns about the phenomenon; only then can a conceptual model for the desired phenomenon be developed, tested, and validated [26]. Modeling of soil erosion requires many empirical measurements of a landscape. Remote sensing techniques that produce aerial photographs, satellite data, or radar data, have eased the difficulty of this otherwise laborious problem [27].

In the last decades, several gully-erosion models have been developed to aid gully- erosion susceptibility mapping (GESM). These models can be classified into three groups: knowledge based models like the analytic hierarchy process (AHP) [16]; bivariate and multivariate statistical models like conditional probability (CP) [28], information value (IV) [29], frequency ratio (FR) [13], index of entropy (IoE) [30], evidential belief function (EBF) [14], weights-of-evidence (WOE) [31], certainty factor (CF) [32], and logistic regression (LR) [33]; and machine-learning models like maximum entropy (ME) [19], multivariate adaptive regression spline (MARS) [15], artificial neural network (ANN) [34], boosted regression tree (BRT) [35], linear discriminant analysis (ADA) [17], bagging best-first decision tree (Bag-BFTree) [36], random forest (RF) [37], flexible discriminant analysis (FDA) [38], support vector machine (SVM) [39], and classification and regression trees (CART) [15].

Each model has disadvantages and advantages. Therefore, improve performance of gully-erosion models; this study uses a hybrid approach to identify gully-erosion susceptible areas. A region that is heavily impacted by gully erosion is the Bastam district watershed in Semnan Province, of northern Iran. Gully erosion occurs in this area due to the highly erodible soils of these lands and misuse of soil and water resources. This has led to the destruction and transformation of agricultural lands into wastelands, the destruction of communication infrastructure, and damage to residential areas. In this study, two models of WoE and IoE are combined to map gully-erosion susceptibility in the Bastam district watershed. These two models have been previously used in studies of other fields of

destructive geomorphological and hydrological phenomena, such as landslides [40–42], floods [43], and groundwater depletion [44,45].

#### **2. Materials and Methods**

#### *2.1. Study Area*

The Bastam watershed covers a 1329 km<sup>2</sup> area (36◦27 02" to 36◦47 13" N, and 54◦24 23" to 55◦11 08" E (Figure 1). Elevation in the study area ranges from 1357 m a.s.l. to 3893 m a.s.l. The mean elevation is 278.34 m a.s.l.). The minimum and maximum slopes are 0◦ and 70.66◦; the mean is 13.55◦. The central and south parts of the watershed have relatively smooth topography with gentle slopes. The rest of the study area is mountainous; is part of the Alborz Mountains. The watershed experiences a mean annual rainfall amount of 262 mm. The mean annual temperature is 12.8 ◦C [46]. More than 25% of the land is poor rangeland. Another 22.84% is covered by a relatively sparse forest, while 0.01% is covered in dense forest and 0.01% is land on which dryland-farming is practiced. The main lithographic units include high elevation piedmont fans and valley terrace deposits, stream channels, braided channels, floodplain deposits, and low elevation piedmont fans, and valley terrace deposits [47]. Rock outcrops, entisols, entisols/inceptisols, and mollisols are the most common geological and pedological surfaces in the study area [48]. Morphometric analysis of gullies in the study area shows that approximately 9.3% of the study area is affected by gully erosion. This indicates that the study area is very susceptible to gully erosion. There are more gullies in the south-central portion of the study area where slopes are lower than in the northern parts of the study area where slopes are steep and rocky outcrops are abundant and few gullies form. Gullies range in length from several meters to several hundred meters and their depths can be several meters. Gully width also varies, ranging from several centimeters to several meters. Gullies in the northern parts of the study area have V-shaped cross-sections, whereas in the central and southern parts of the region, due to more erodible soils, concentrated runoff because of low slope, and the more resistant sediments on the floor of the watershed prevent erosion, so the valleys are U-shaped.

**Figure 1.** Location of the study area in Iran and Semnan province and locations of gullies used for training (black dots) and validation (green dots).

#### *2.2. Methodology*

This study involved six steps (Figure 2): data preparation, including the creation of a GEIM, compilation of thematic layers depicting the spatial distribution of gully-erosion conditioning factors (GECFs), and division of GEIM into training and validation groups; multi-collinearity analyses of the GECFs using indices of tolerance (TOL) and variance inflation factor (VIF); application of the WoE model to analyze the spatial relationship between GECFs and gully locations; determination of the

importance of the GECFs using the IoE model; preparation of the GESM using the WoE and IoE models separately; and integration of the hybrid WoE-IoE model; and validation using the area-under-the-curve receiver-operating characteristic (AUROC).

**Figure 2.** Flowchart of research in the present study for gully erosion susceptibility mapping.

## 2.2.1. Gully-erosion Inventory Map (GEIM)

A GEIM shows the spatial distribution of gullies and this is necessary to determine the spatial distribution of susceptibility to gullying [15]. Data acquired from the Agricultural and Natural Resources Research Center of Semnan Province were combined with extensive field surveys (Figure 3) and interpretation of satellite imagery to map gully erosion in the study area. Three hundred and three gullies were identified in the study area (Figure 1) and this set was randomly separated into two groups for modeling (70% or 212 gullies) and validation (30% or 91 gullies). In addition, an equivalent number (303) of locations without gullies were randomly selected using the random-point tool in ArcGIS to support calibration and validation [49]. The maximum and minimum lengths of the 303 gullies are 346 m and 0.74 m, respectively, whereas, the maximum depth is 7.2 m and the minimum depth is 0.65 m. The maximum width is 16.6 m and the minimum width is 0.74 m.

**Figure 3.** Sample of mapped gullies in the study area.

#### 2.2.2. Gully Erosion Conditioning Factors (GECFs)

Gullying results from the interaction of several environmental factors that include hydrology, topography, land use/land cover, soil characteristics, human activities, and rainfall [50]. In this study, 18 GECFs were selected and mapped: elevation (m), slope (◦), plan curvature (100/m), aspect, convergence index (100/m), slope length (LS), profile curvature (100/m), topography position index (TPI), terrain ruggedness index (TRI), topography wetness index (TWI), stream power index (SPI), rainfall (mm), distance to road (m), distance to stream (m), drainage density (km/km2), land use/land cover (LU/LC), soil type, and lithology (Figure 4a–r). These factors were obtained from several sources. The geological map at the scale of 1:100,000 and made by the Geological Society of Iran (GSI) (http://www.gsi.ir/) was used to produce the lithology map. Ten lithological units in the study area were created based upon formation and susceptibility to gully erosion (Table 1). The topographic map at the scale of 1:50,000 were acquired from the National Geographic Organization of Iran (www.ngo-org.ir), and Google Earth images were used to digitize the road network on the topographic map. The ALOS DEM with 12.5 m resolution, downloaded from the Alaska Satellite Facility (ASF) Distributed Active Archive Center (DAAC) was used to extract topographic and hydrological data to map elevation, slope, aspect, plan curvature, profile curvature, slope length (LS), TWI, SPI, TPI, TRI, CI, distance to stream, and drainage density [29,51]. The reproduction of complex morphology and features depends both on accuracy and gridding techniques [52,53]. The quality of the reproduction influence the value of some of the topographical and hydrological gully-erosion conditioning factors, Therefore, in this research, ALOS DEM with a vertical accuracy of 0.3 m was used. Similar to the accuracy assessment procedures implemented by [54], vertical accuracies of the ALOS DEM were assessed by comparing the ALOS DEM elevations with those of the ground control points (GCPs). At each point, the DEM elevations were extracted using ArcGIS 10.5 software. Then, the differences in elevation were computed by subtracting the GCP elevation from its corresponding DEM elevations, and these differences are the measured errors in the ALOS DEM. For a particular DEM, positive errors represent locations where the DEM was above the GCP elevation, and negative errors occur at locations where the DEM was below the control point elevation. From these measured errors, the mean error and RMSE for each DEM were calculated, including standard deviations of the mean errors. The mean error (or bias) indicates if a DEM has an overall vertical offset (either positive or negative) from true ground level [54]. Finally, accuracy assessment results were analyzed. Details on how the ALOS DEM was produced using interferometric synthetic aperture radar (InSAR) are discussed in the papers of References [55,56]. The most important step in InSAR for DEM generation is the phase measurement, then the transformation of phase to height [55].

LU/LC was extracted from a land-use map prepared by the Iranian Soil Conservation and Watershed Management Research Institute (https://www.scwmri.ac.ir). The map of soil types was prepared by the Agricultural and Natural Resources Research Center of Semnan Province. To map annual rainfall, 30 years of annual data (1986 to 2016) were acquired for the meteorological stations at Mojen, Farahzad, Bastam, Abr, Karkhaneh, Semnan, Shahroud, Tarzeh, and Shah Kuh-e Bala. Kriging was used to prepare the final annual rainfall map in ArcGIS 10.5. ArcGIS V10.5, Arc Hydro V10.4, Google Earth Pro 7.8, and ENVI v4.8 toolboxes were used to prepare the conditioning factors for gully-erosion susceptibility assessment. Microsoft Excel v2016 and SPSS v24 were also used for exploratory, statistical, and validation analyses. Calculation of the GECFs is explained elsewhere [57–61].

**Figure 4.** Gully erosion conditioning factors. (**a**) elevation, (**b**) slope, (**c**) plan curvature, (**d**) slope aspect, (**e**) convergence index (CI), (**f**) slope length (LS), (**g**) profile curvature, (**h**) topography position index (TPI), (**i**) terrain ruggedness index (TRI), (**j**) topography wetness index (TWI), (**k**) stream power index (SPI), (**l**) rainfall, (**m**) distance to road, (**n**) distance to stream, (**o**) drainage density, (**p**) land use/land cover, (**q**) soil type, and (**r**) lithology.


#### **Table 1.** The lithology of the study area.

## *2.3. Multi-collinearity Analysis*

Collinearity indicates that an independent variable is the linear function of another independent variable. If collinearity in a regression equation is high, there is a high correlation between independent variables and that reduces the accuracy of the model [15]. In this study, the indicators tolerance (TOL) and variance inflation factor (VIF) were used to assess collinearity among the independent variables. These indices do not have standardized thresholds, however, the literature reports that the following intervals have been widely used by several researchers: VIF ≤ 5 or 10 and TOL ≤ 0.1 or 0.2, imply that no collinearity is present and the gully conditioning factors are independent [15].

#### *2.4. Models*

#### 2.4.1. Index of Entropy (IoE)

Entropy is a measurement of the instability, imbalance, disorder, and uncertainty of a system [62]. The value of entropy of a system has a one-to-one relationship with the degree of disorder. This relationship, called the Boltzmann principle, has been used to describe the thermodynamic status of a system [62]. Shannon improved upon the Boltzmann principle and established an entropy model for information theory. The IoE distinguishes the most important factors from less significant effective factors of a target; it identifies the variables that have the greatest impact on the occurrence of an event. While gully erosion is affected by several factors and gully-erosion susceptibility can be determined with a bivariate statistical model like the WoE, where all factors are weighted the same. Factors having a greater impact might be ignored. Therefore, the IoE can provide an important understanding of the conditioning factors and their impacts and the results can inform appropriate management [63]. The entropy of gully erosion refers to the extent that various factors influence the development of a gully. Several important factors provide additional entropy into the index system. As a result, the entropy value can be used to calculate the objective weights of the index system. To determine the relative importance of GECFs in the gully occurrence and GESM using IoE model, following Equations (1)–(6) have been used [63]:

$$P\_{\text{ij}} = \frac{\mathbf{b}}{\mathbf{a}} \tag{1}$$

$$\begin{pmatrix} \mathbf{P\_{ij}} \end{pmatrix} = \frac{\mathbf{P\_{ij}}}{\sum\_{\mathbf{j=1}}^{\mathbf{S}\mathbf{j}} \mathbf{P\_{ij}}} \tag{2}$$

$$\mathbf{H}\_{\mathbf{j}} = -\sum\_{\mathbf{i}=1}^{\text{S}\dagger} (\mathbf{P}\_{\overline{\mathbf{i}}}) \log\_2 \left( \mathbf{P}\_{\overline{\mathbf{i}}} \right)\_{\mathbf{i}} \qquad \mathbf{j} = 1, \; \dots, \mathbf{n} \tag{3}$$

$$\mathbf{H}\_{\mathbf{j}\text{ max}} = \log\_2 \mathbf{S}\_{\mathbf{j}} \tag{4}$$

$$\mathbf{I}\_{\mathbf{j}} = \frac{\mathbf{H}\_{\mathbf{j}\text{ max}} - \mathbf{H}\_{\mathbf{j}}}{\mathbf{H}\_{\mathbf{j}\text{ max}}}, \text{ I} = (0, 1), \qquad \mathbf{j} = 1, \dots, \mathbf{n} \tag{5}$$

$$\mathcal{W}\_{\vec{\text{ij}}} = \mathbf{I}\_{\vec{\text{j}}} \times \mathbf{P}\_{\vec{\text{ij}}} \tag{6}$$

where a and b are the domain and gully erosion percentages, respectively, Pij is the probability density, Hj and Hj max indicate the entropy values and maximum entropy, respectively, Ij is the information value, Sj is the number of classes and Wj indicates the resulting weight of each factor. Wij's values range from 0 to 1. After calculation of the final weight of each GECF and of their classes, these values were added to each related thematic layer and then each was weighted. Finally, they were summed to produce the final gully-erosion susceptibility map using Equation (7) and the weighted-sum tool in ArcGIS [43].

$$\text{GESM} = \sum\_{\mathbf{l}=1}^{n} \left( \mathbf{W}\_{\mathbf{l}} \times \mathbf{P}\_{\mathbf{l}} \right) \tag{7}$$

where Wj and Pj are the final weight and the probability density for the jth feature.

#### 2.4.2. Weight of Evidence (WoE) Model

WoE is a bivariate statistical method, based on the Bayesian probability framework, to statistically estimate the relative importance of conditioning factors [41]. By overlaying gully erosion locations on each gully-related conditioning factor, the spatial relationship between them can be determined and the explanatory significance of the effective variable for past gullying can be evaluated. The WoE model calculates the positive (W+) and negative (W−) weights for each gully conditioning factor (A) based on the presence or absence of gully sites (B). This model measurement the conditioning factors within the study area as follows [64]:

$$\mathbf{W}\_{\rm i}^{+} = \text{In}\left(\frac{\text{p}\langle\text{B}|\text{AL}\rangle}{\text{P}\left\{\text{B}\middle|\overline{\text{A}}\right\}}\right) \tag{8}$$

$$\mathbf{W}\_{\mathbf{i}}^{-} = \mathrm{Im} \left( \frac{\mathrm{P} \{ \overline{\mathbf{B}} | \mathbf{A} \}}{\mathrm{P} \left\{ \overline{\mathbf{B}} | \overline{\mathbf{A}} \right\}} \right) \tag{9}$$

P is the probability and In is the natural log function. B and B indicate the presence and absence of the gully conditioning factors, respectively. A is the presence of gully, and A is the absence of a gully. A positive weight (W+) designates the fact that the conditioning factor is present at the gully locations and its value is an indication of the positive correlation between the presence of the gully conditioning factor and the gullies. Similarly, a negative weight (W−) explains the absence of the

gully conditioning factor and reflects the level of negative correlation. In GESM, the weight contrast (C) measures and specifies the spatial association between the effective factors and gully erosion occurrences. C is negative for a negative spatial relationship and positive for a positive relationship. The standard deviation S(C) of W is calculated by Equation (10):

$$\mathbf{S(C)} = \sqrt{\mathbf{S}^2} \mathbf{(W^+)} + \mathbf{S^2} (\mathbf{W^-}) \tag{10}$$

where S<sup>2</sup> W+ is the variance of the W<sup>+</sup> and S2(W−) is the variance of W−. The variances of the weights can be determined as follows:

$$\mathrm{S}^{2}\mathrm{\{W}^{+}\} = \frac{1}{\mathrm{N}\{\mathrm{B}\cap\mathrm{L}\}} + \frac{1}{\{\mathrm{B}\cap\mathrm{L}\}}\tag{11}$$

$$\mathbf{S}^2(\mathbf{W}^-) = \frac{1}{\left\{\overline{\mathbf{B}} \cap \mathbf{L}\right\}} + \frac{1}{\left\{\overline{\mathbf{B}} \cap \overline{\mathbf{L}}\right\}}\tag{12}$$

The studentized contrast (G*final*) is a measure of confidence and is calculated, using the following equation:

$$\mathbf{G}\_{final} = \left(\frac{\mathbf{C}}{\mathbf{S(C)}}\right) \tag{13}$$

C indicates the overall association between a geo-environmental factor and gully occurrence. S(C) is the standard deviation of the contrast and W is the final weight.

After determining the relative weights of the GECFs using IoE model and the spatial relationships between GECFs and gullies in the study area using the WoE model, the two models are integrated to improve performance and decrease the disadvantages of each so that, relative weight of GECFs obtained by IoE multiple with weight of GECF classes obtained by WoE using Equation (14):

$$\begin{aligned} \text{GESM}\_{\text{WoE}-\text{I}\text{oE}} &= \left( \text{WoE}\_{\text{Elevation}} \times \text{Elevation}\_{\text{IoE}} + \text{WoE}\_{\text{Slope}} \times \text{Slope}\_{\text{IoE}} + \text{V}\_{\text{OE}} \times \text{s}\_{\text{T}} \right) \\ \text{WoE}\_{\text{Aspect}} &\times \text{Aspert}\_{\text{IoE}} + \text{WoE}\_{\text{SPI}} \times \text{SPI}\_{\text{IoE}} + \text{WoE}\_{\text{TVI}} \times \text{T} \text{WI}\_{\text{IoE}} + \text{WoE}\_{\text{TPI}} \times \text{TPI}\_{\text{IoE}} + \\ &\dots + \text{WoE}\_{\text{Lithology}} \times \text{Lithology}\_{\text{IoE}} \text{)}. \end{aligned} \tag{14}$$

#### *2.5. Validation Method*

Validation of the results is an important step in GESM [65]. The AUROC (area under the receiver operating characteristic) method is an efficient and accurate validation tool of GESMs [66]. This approach helps to visualize the quality of models by showing the incremental percentage of the model's prediction accuracy [67]. In the AUROC approach, the number of pixels correctly predicted by the model is plotted against the number of pixels incorrectly predicted. The AUROC is usually between 0.5 and 1; as the value approaches 1, the higher is the performance of the model [68]. AUROC values can be classified as performance ratings: 0.5–0.6 is poor, 0.6–0.7 is average, 0.7–0.8 is good, 0.8–0.9 is very good, and 0.9–1 is excellent [69].

#### **3. Results**

## *3.1. Multi-collinearity Analysis Among GECFs*

Among the 21 GECFs studied, three factors (NDVI with tolerance (or TOL) = 0.095 and variance inflation factor (or VIF) = 10.47), distance to fault (TOL = 0.022 and VIF = 45.04), and soil texture (TOL = 0.03 and VIF = 33.76) have collinearity and cannot be used for modeling (Table 2). The values of TOL and VIF for the other nineteen factors ranged from 0.271 to 0.975 and 1.02 to 4.17, respectively, which indicates non-linearity between them. Therefore, nineteen factors were included in the modeling.


**Table 2.** Multi-collinearity test among gully erosion conditioning factors.

#### *3.2. Determine the Spatial Relationship between GECFs and Occurred Gullies*

The analysis of the spatial relationship between GECFs and the gullies occurred in the study area was conducted first using the IoE and WoE models (Table 3 and Figure 5). Values less than 1 in the IoE model indicate low correlations and values larger than 1 indicate higher correlations. In general, high values of IoE indicate higher probability of gully occurrence and lower values indicate less likelihood that gullying will occur [16]. Negative values generated by the WoE model indicate negative spatial relationships and positive values reflect positive relationships. Zero indicates that the specific class of conditioning factors is not useful in the analysis [70]. Classes of the factors elevation, slope, plan curvature, aspect, and CI have strong relationships with gullying in the IoE model: classes of <1659 m = 2.495, <5◦ = 2.499, flat curvature = 1.565, east-facing = 1.941, and CI < −39.6 = 2.831. In the WoE model, these same classes also show strong relationships with sites of gullying: <1659 m = 12.43, <5◦ = 5.746, flat curvature = 4.156, east-facing = 5.914, and CI < −39.6 = 6.463. With regard to the factors in the IoE model, LS, profile curvature, TPI, TRI, and TWI have strongly predictive classes: 50.6 to 69.5 m = 1.490, −0.43 to 0.28 = 1.425, −3.3 to 2.9 = 1.607, <3.22 = 2.134, and > 10.85 = 3.009. These same factors and class are similarly strongly correlated with gullying in the WoE model: 50.6 to 69.5 m = 2.880, −0.43 to 0.28 = 7.096, −3.3 to 2.9 = 4.849, <3.22 = 11.011, and >10.85 = 6.226. In the case of the IoE model, the factors SPI, rainfall, distance to road, distance to stream, and drainage density have classes that strongly indicate locations of highest susceptibility to gully erosion: > 15.35 = 1.910, <248.5 mm = 2.461, <500 m=2.214, <100 m = 1.627, and >2.19 km/km2 = 3.385. In the WoE model these same classes indicated high susceptibility as well: > 15.35 = 2.636, <248.5 mm = 9.641, < 500 m = 7.139, <100 m = 5.500, and >2.19 km/km2 = 15.005). The IoE model aligns classes of LU/LC, soil type, and lithology with higher susceptibility: poor rangeland = 2.583, entisols/inceptisols = 2.933, and Group 9 (consisting of quaternary formations that include high elevation piedmont fan and valley terrace deposits, stream channel, braided channel, and floodplain deposits, low elevation piedmont fan and valley terrace deposits = 2.237. The same classes in the WoE model are highly predictive of gullying, as well: poor rangeland = 11.992, entisols/inceptisols = 11.522, and group 9 = 11.297.

**Figure 5.** The spatial relationship between conditioning factors and gully locations.


**Table 3.** The spatial relationship between gully erosion conditioning factors and gully locations.


**Table 3.** *Cont.*

#### *3.3. Determining the Relative Importance of GECFs in the IoE Model*

The order of importance of the GECFs in the IoE model (Figure 6) is: drainage density (0.564), slope (0.48), rainfall (0.448), TRI (0.427), soil type (0.391), elevation (0.383), TWI (0.373), TPI (0.309), LU/LC (0.248), profile curvature (0.246), lithology (0.24), distance to road (0.136), CI (0.132), distance to stream (0.13), slope aspect (0.117), plan curvature (0.107), SPI (0.105), and LS (0.024).

**Figure 6.** The relative importance of conditioning factors in the IoE model. A: landuse/land cover, B: Soil type, C: lithology, D: topography wetness index (TWI), E: stream power index (SPI), F: aspect, G: convergence index (CI), H: elevation, I: drainage density, J: distance to stream, K: distance to stream, L: slope length (LS), M: plan curvature, N: profile curvature, O: rainfall, P: slope, Q: terrain ruggedness index (TRI), R) topography position index (TPI).

#### *3.4. Gully Erosion Susceptibility Mapping (GESM)*

After the computation of the spatial relationship between gully erosion and conditioning factors, and determination of the relative importance of GECFs in the IoE model, GESM calculations based on the WoE model Equation (15), the IoE model Equation (16), and the hybrid WoE- IoE model Equation (17) were constructed.

GESMWoE = (WWoEElevation) + (WWoESlope ) + (WWoEPlan curvature) +(WWoESlope aspect) + (WWoECI) + (WWoELS) +(WWoEProfile curvature) + (WWoETPI) + (WWoETRI) +(WWoETWI) + (WWoESPI) + (WWoERainfall) +(WWoEDistance to road) + (WWoEDistance to stream) +(WWoEDrinage density) + (WWoELU/LC) + (WWoESoil type) +(WWoELithology) (15) GESMIoE = (WIoEElevation × 0.383) + (WIoESlope × 0.48) +(WIoEPlan curvature × 0.107) + (WIoESlope aspect × 0.117) +(WIoECI × 0.132) + (WIoELS × 0.024) +(WIoEProfile curvature × 0.246) + (WIoETPI × 0.309) +(WIoETRI × 0.427) + (WIoETWI × 0.373) + (WIoESPI × 0.105) +(WIoERainfall × 0.448) + (WIoEDistance to road × 0.136) +(WIoEDistance to stream × 0.13) + (WIoEDrinage density × 0.565) +(WIoELU/LC × 0.248) + (WIoESoil type × 0.391) +(WIoELithology × 0.24) (16) GESMWoE−IoE = (WWoEElevation × 0.383) + (WWoESlope × 0.48) +(WWoEPlan curvature × 0.107) + (WWoESlope aspect × 0.117) +(WWoECI × 0.132) + (WWoELS × 0.024) +(WWoEProfile curvature × 0.246) + (WWoETPI × 0.309) +(WWoETRI × 0.427) + (WWoETWI × 0.373) +(WWoESPI × 0.105) + (WWoERainfall × 0.448) +(WWoEDistance to road × 0.136) + (WWoEDistance to stream × 0.13) +(WWoEDrinage density × 0.565) + (WWoELU/LC × 0.248) + (WWoESoil type × 0.391) + (WWoELithology × 0.24). (17)

Values of gully erosion susceptibility for the WoE model, the IoE model, and the hybrid WoE-IoE model varied for the minimums (−77.558, 0.265, and −23.3362) and the maximums (143.404, 13.810, and 51.6535) of each, respectively. These values were classified into five classes using the natural breaks method: very low, low, moderate, high, and very high susceptibility classes (Figure 7a–c). The results for each category using the WoE model were 41.45%, 14.59%, 10.91%, 18.5%, and 14.52%, respectively (Figure 8). The results for each category using the IoE model were 40.03%, 15.43%, 11.10%, 17.18%, and 16.22%, respectively (Figure 8). And using WoE-IoE integrated model, the classes of susceptibility were 27.95%, 23.78%, 10.22%, 16.60%, and 21.42%, respectively.

**Figure 7.** Gully erosion susceptibility map. (**a**) weight-of-evidence (WoE) model, (**b**) index-of-entropy (IoE) model, and (**c**) WoE-IoE integrated model.

**Figure 8.** Percentage of each susceptibility classes in three models of the weight of evidence, index of entropy, and WoE-IoE.

#### *3.5. Validation of Results*

The AUROC graphs were created for the training dataset (success-rate curve) and for the validation dataset (prediction-rate curve) (Figure 9a,b). These demonstrate that the WoE model had better performance in terms of the success-rate curve (SRC = 0.899) than did the IoE model (SRC = 0.896). The results also indicate that integration of these models improve their performance: the WoE-IoE integrated model had a better prediction-rate curve (PRC = 0.903) and a better success-rate curve (SRC = 0.918) than either the WoE (PRC = 0.877 and SRC = 0.877) or IoE (PRC = 0.899 and SRC = 0.896).

**Figure 9.** Area under curve values for weight-of-evidence (WoE), index-of-entropy (IoE), and WoE-IoE models: (**a**) validation dataset (prediction rate curve), (**b**) training dataset (success rate curve).

#### **4. Discussion**

In terms of the position of the gully in the landscape [71], most of the gullies are on the floor of the valley and in terms of the evolution of the gullies [72], most are continuous. The continuous gully is part of the drainage network, and a discontinuous gully is separated from the drainage network. In terms of gully head-cut plan [73,74], the gullies in the study area are the digitized type, usually

found adjacent to rivers and often located at the intersection of their branches. These gullies form in areas of 0% to 5% slope. The most important factors that led to the development of gullies in the study area are: the undercutting of the boundaries of flood flows and gravitational formation of mass fractures; formation of a groove in gullies due to rill erosion and the resulting surface runoff; the influences of tunnel erosion (piping) and underground corridors which have a maximum diameter of 8 meters and a maximum length of 12 m—the expansion of these corridors and the collapse of their roofs cause gully development; and head-cut retrogradation and upward development of gullies.

Determination of the susceptible areas affected by different kinds of soil erosion such as gully erosion is useful for land use planning, mitigating soil erosion, and conservation practices [29]. So far, various methods have been applied for gully erosion susceptibility assessment globally [28–30]. Given the shortcomings and limitations of each of the models, scientists have proposed and developed integrated methods to overcome their disadvantages and increase their efficiency [14]. In this study, two types of statistical methods, namely, IoE and WoE, and their ensemble were applied to produce GESM. The IoE–WoE integrated method eliminates the disadvantages of IoE and WoE individual models and improves their advantages. The main advantages of the WoE method are that it calculates the weighted value of the factors based on a statistical formula and thus avoids the subjective choice of weighting factors. In addition, input maps with missing data (incomplete coverage) can be accommodated in the model and under-sampled data do not significantly impact the results. The main shortcomings of the model are the failure to calculate the relative importance of the parameters, and that the weighted values calculated for different areas are not comparable in terms of the degree of hazard [31]. The main advantage of IoE model is that this model does not require that assumptions be made about the proper distribution of explanatory variables; therefore, several properties can be used and tested. This method also examines the statistical relationships between independent and dependent variables and provides metrics for the significance of variables [30].

Analyses of the spatial relationships between GECFs and gully locations using IoE, and WoE show that areas of low elevation and slope and flat topography, where surface runoff concentrates and where erosion-sensitive evaporation deposits of gypsum or salt have formed, have a high susceptibility to gully erosion. This has been demonstrated in other local studies [69,75]. Areas nearer to streams and roads, with sparse vegetation and higher drainage density than other areas, have more potential for gully occurrence. These findings are in line with References [22,76]. Limestone, sandstone, marl, shale, and red conglomerate geological units have a high susceptibility to gully erosion. These geomorphological features are generally known to promote gully erosion; this is confirmed within the study area and has been reported in other research [77,78].

The IoE model has shown that drainage density, slope degree, and rainfall are key conditional factors for gully erosion in the study area. This is in line with References [16,63,79]. Yesilnacar et al. [79] proved that gently sloping areas are susceptible to surface flow accumulation and gully erosion. Pourghasemi et al. [63] assessed the capability of WoE and FR models for spatial prediction of gully erosion susceptibility in the Chavar region of Ilam Province, Iran. Sixty-three gullies and ten GECFs were used in that study. Their results indicate that the distance to roads, drainage density, and LU/LC are key conditions affecting gully occurrence. Arabameri et al. [16] compared three data-driven models and an AHP knowledge-based technique for gully-erosion susceptibility mapping in the Toroud watershed in Semnan Province, Iran using 80 gully locations and 13 GECFs (including elevation, slope degree, slope aspect, plan curvature, distance from river, drainage density, distance from road, lithology, LU/LC, TWI, SPI, NDVI, and LS). Their results show that lithology, slope, and NDVI are the primary determinants of gully occurrence in the Toroud watershed.

Validation results show that the IoE model performed better than the WoE model. This is consistent with the work of others [16,80,81]. Arabameri et al. [16] States that IoE with an SRC = 0.939 and a PRC = 0.925 better predicts areas prone to gully erosion than does WoE with an SRC = 0.926 and a PRC = 0.921. Wang et al. [81] used IoE and FR models for groundwater qanat potential mapping in the Moghan watershed, Iran and their results depict the excellence of IoE model in qanat occurrence potential

estimation. The integration of WoE and IoE has improved upon their separate performances, which is also consistent with References [14,39,82–88]. Arabameri et al. [14] used EBF, LR, and a new ensemble EBF–LR algorithm to spatially model gully erosion at Semnan Province, Iran. Their results show that their ensemble method performed considerably better (AUC = 0.909) than did the individual LR (0.802) and EBF (0.821) methods. Pourghasemi et al. [39] assessed the individual and ensemble data-mining techniques for gully erosion modeling and stated that the ensemble models had a higher goodness-of-fit and predictive power than individual models. Arabameri et al. [84] compared the performance of individual and ensemble models for assessment of landslide susceptibility in Sangtarashan watershed, Mazandran Province, Iran and state that FR–RF integrated model (AUC = 0.917) achieved higher predictive accuracy than the individual FR (AUC = 0.865) and RF (AUC = 0.840) models. Du et al. [85] integrated IV and LR individual models for landslide susceptibility mapping in the Bailongjiang watershed, Gansu Province, China and state that the proposed integrated method was reliable to produce an accurate landslide susceptibility map. [88] Used ensemble RF and EBF models for landslide susceptibility assessment in Western Mazandaran Province, Iran and state that introduced ensemble model can be a powerful tool for landslide assessment at regional scales.

Our research will contribute to achieve a better knowledge of the landscape and to develop sustainable policies to achieve the Land Degradation Neutrality challenge and the United Nation Goals for Land Degradation [89,90]. This is a relevant issue to achieve sustainable land management where gullies must be restored when developed as a consequence of human mismanagement, and for this is necessary to use nature-based solutions [91]. To achieve success in gully erosion control the strategies must find a way to reduce the connectivity of the flows [92].

#### **5. Conclusions**

Gully erosion is a common geomorphological problem in arid and semi-arid regions, therefore, it is essential to develop methods for predicting gullying with highly accurate and effective models. Knowing the locations that are prone to or susceptible to gullying can enable ways to avoid casualties and financial losses caused by gully development, and can even promote sustainable development. In recent years, many quantitative and qualitative methods have been introduced for GESM. No approach is believed to be a best-approach, but the general consensus is that each method has its advantages and disadvantages. In this study, two models (the IoE and the WoE) and their integrated offspring (WoE-IoE) were tested for GESM to assess what their advantages and shortcomings were. The most important conclusions to be drawn from this study are:


**Author Contributions:** Conceptualization, A.A Data curation, A.A.; Formal analysis, A.A. Investigation, A.A and A.C; Methodology, A.A., and A.C; Resources, A.A and A.C.; Software, A.A.; Supervision, A.A and A.C; Validation, A.A; Writing–original draft, A.A.; Writing–review and editing, A.A., A.C., and J.P.T.

**Funding:** This research received no external funding.

**Acknowledgments:** We are grateful to the Editor, Miley Fan, and two anonymous referees for their constructive comments which were valuable to improve our manuscript.

**Conflicts of Interest:** The authors declare no conflict of interest.
