Comparative assessment using boosted regression trees, binary logistic regression, frequency ratio and numerical risk factor for gully erosion susceptibility modelling

Catena 183 (2019) 104223

Contents lists available at ScienceDirect

Catena journal homepage: www.elsevier.com/locate/catena

Comparative assessment using boosted regression trees, binary logistic regression, frequency ratio and numerical risk factor for gully erosion susceptibility modelling ⁎

T

⁎⁎

Alireza Arabameria, , Biswajeet Pradhanb,c, , Luigi Lombardod a

Department of Geomorphology, Tarbiat Modares University, Tehran 36581-17994, Iran Centre for Advanced Modelling and Geospatial Information Systems (CAMGIS), Faculty of Engineering and Information Technology, University of Technology Sydney, 2007, New South Wales, Australia c Department of Energy and Mineral Resources Engineering, Sejong University, Choongmu Gwan, 209 Neungdong-ro Gwangjin-gu, Seoul 05006, Republic of Korea d University of Twente, Faculty of Geo-Information Science and Earth Observation (ITC), the Netherlands b

A R T I C LE I N FO

A B S T R A C T

Keywords: Gully erosion susceptibility GIS Boosted regression trees Binary logistic regression Bayazeh watershed

The initiation and development of gullies as worldwide features in landscape have resulted in land degradation, soil erosion, desertification, flooding and groundwater level decrease, which in turn, cause severe destruction to infrastructure. Gully erosion susceptibility mapping is the first and most important step in managing these effects and achieving sustainable development. This paper attempts to generate a reliable map using four state-of-theart models to investigate the Bayazeh Watershed in Iran. These models consists of boosted regression trees (BRT), binary logistic regression (BLR), numerical risk factor (NRF) and frequency ratio (FR), which are based on a geographic information system (GIS). The gully erosion inventory map accounts for 362 gully locations, which were randomly divided into two groups (70% for training and 30% for validation). Sixteen topographical, geological, hydrological and environmental gully-related conditioning factors were selected for modelling. The threshold-independent area under receiver operating characteristic (AUROC) and seed cell area index (SCAI) approaches were used for validation. According to the results of BLR and BRT, the conditioning parameters namely, NDVI and lithology, played a key role in gully occurrence. Validation results showed that the BRT model with AUROC = 0.834 (83.4%) had higher prediction accuracy than other models, followed by FR 0.823 (82.3%), NRF 0.746 (74.6%) and BLR 0.659 (65.9%). SCAI results indicated that the BRT, FR and BLR models had acceptable classification accuracy. The findings, in terms of model and predictor choice, can be used by decisionmakers for hazard management and implementation of protective measures in gully erosion-prone areas.

1. Introduction Land degradation is a worldwide process induced by human activity on landscapes (Solomun et al., 2018; Hălbac-Cotoară-Zamfir et al., 2019; Arabameri and Pourghasemi, 2019). In the last two decades, soil erosion has received increasing attention within the scientific community (Romshoo et al., 2012; Issaka and Aqeel Ashraf, 2017; Arabameri et al., 2018c, 2018d; Keesstra et al., 2016, 2019; Kropacek et al., 2016). As the global population rapidly increases, the need for food sources will proportionally increase, and soil represents the natural support required in agriculture and farming (Magliulo, 2012). However, natural and anthropogenic activities contribute to an increment in estimated

erosion rates (Arabameri et al., 2017a, 2017b; Vaezi et al., 2017; Antoneli et al., 2018). Soil is not renewable (Rodrigo-Comino et al., 2018) within a human's lifespan (Ding et al., 2015; Kumar Samanta et al., 2016; Vanmaercke et al., 2016). Thus, an understanding of the erosional processes under current and future conditions (time) is as crucial as it is understanding the erosional processes occurring around the world (space) (Keesstra et al., 2018). To tackle these questions, the community usually implements physics- (Bout et al., 2018) or statistically- (Castro Camilo et al., 2017) based predictive models defined over time, space or both (Borrelli et al., 2017). The statistical component is often referred to as susceptibility modelling (Lombardo et al., 2018, 2019).

⁎

Corresponding author. Correspondence to: B. Pradhan, Centre for Advanced Modelling and Geospatial Information Systems (CAMGIS), Faculty of Engineering and Information Technology, University of Technology Sydney, 2007, New South Wales, Australia. E-mail addresses: [email protected] (A. Arabameri), [email protected] (B. Pradhan). ⁎⁎

https://doi.org/10.1016/j.catena.2019.104223 Received 19 March 2019; Received in revised form 11 June 2019; Accepted 19 August 2019 0341-8162/ © 2019 Elsevier B.V. All rights reserved.

Catena 183 (2019) 104223

A. Arabameri, et al.

32°51′56″–33°43′42″ N latitude. The minimum and maximum elevation of the study area is bounded between 644 and 2291 m.a.s.l, and the mean elevation is 937.49 m.a.s.l. The average annual temperature and rainfall during a 30-year period (1987–2017) are 17.23C° and 87.64 mm, respectively (IRIMO, 2012). The study area lies across arid and semi-arid regions. The mean slope is 4.42°, which shows a gentle topography. Poor range (37.38%), sand dune (13.58%) and salt lands (13.01%) are the most important land use/land cover (LU/LC) classes. Among the 103 villages in the study area, the most prominent are Bayazeh, Anarak, Asghar Aabad, Bahin, Hemat abad, Eraj, Eshgh abad, Ali abad and Garmeh. Aridisol (32.71) is the most important type of soil in the study area. The lithology of the study study area is listed in Table 1.

Soil erosion susceptibility modelling is the tool used for predicting potential areas subjected to soil erosion and is often a mandatory requirement in implementing actions and limiting land degradation processes at the site scale (Rahmati et al., 2017a; Pourghasemi et al., 2017). In terms of model structure, the community has already established a shared series of actions for assessing soil erosion susceptibility (such as the binary structure of the data and spatial partition). In recent years, quantitative approaches spanning from pure statistics to data mining have been continuously developing (Märker et al., 2011). This advancement thus calls into question the choice of algorithms used for prediction modelling. Hence, comparative studies have become popular in the literature over the years (Conoscenti et al., 2014). Gully erosion, as the main representative of erosional processes in spatial predictive models, aligns well with the current trend. Many algorithms have been separately tested (Märker et al., 2011; Conforti et al., 2011; Dewitte et al., 2015; Rahmati et al., 2016; Rahmati et al., 2017a; Pourghasemi et al., 2017; Zabihi et al., 2018; Arabameri et al., 2018a; Azareh et al., 2019; Arabameri and Pourghasemi, 2019; Arabameri and Pourghasemi, 2019) to produce susceptibility models. Each one is characterised by advantages and drawbacks (Rahmati et al., 2017a; Arabameri et al., 2018c). Therefore, even for gully erosion susceptibility, comparative studies are currently the sole means of assessing which method performs better in such a variable context. The scientific effort may appear redundant at times, but comparative studies can help the community converge on a common approach (or a small set of approaches) and ultimately lead to meta-analytic studies, given that comparative works are common in many other scientific branches (Shit et al., 2015) that already have a good agreement or basic reference. Our contribution well aligns with other comparative studies in natural hazard (e.g. Arabameri et al., 2019b) for it addresses this issue by comparing four state-of-the-art methods in the framework of gully erosion susceptibility. Our comparison included one heuristic, one empirical, one statistical and one machine learning approach (Mohammady et al., 2012; Zare et al., 2013; Dehnavi et al., 2015; Aghdam et al., 2016; Tien Bui et al., 2016; Hong et al., 2017; Mojaddadi et al., 2017; Chen et al., 2019). Specifically, we tested i) numerical risk factor (NRF), a simple heuristic method; ii) frequency ratio (FR, Rahmati et al., 2016; Meliho et al., 2018), an empirical method with a well-established literature; iii) binary logistic regression (BLR, Chaplot et al., 2005; Akgun and Türk, 2011; Dewitte et al., 2015), a multivariate statistical method and iv) boosted regression tree (BRT, Arabameri et al., 2018a), a data mining method. We selected these tools simply because we only searched for approaches with increasing complexity. Thus, we attempted to verify whether highly complex and recent methods would provide improved and accurate spatial information. In this regard, we considered enhancement in information—not only the final prediction, be it numerical or in map form, but also the ability to interpret causal relationships between gullies and the chosen set of causative factors. In our view, understanding such causal relationships should be part of the toolbox represented by any spatial model, although this ideal is not always the case. In the present study, we focused on a relatively small area located at the centre of Iran, which has well-exposed gullies developed in recent years. The structure of this manuscript includes the following. A description of the study area and processes is given in Section 2. An overview of the data and the models we adopted is presented in Section 3, and results, discussions and the conclusion are presented in Sections 4, 5 respectively.

2.2. Methodology This research consists of five stages (Fig. 2). First, two databases were provided. The first database included the determined location of the gullies and the preparation of a gully erosion inventory map (GEIM) from different sources. These gullies were subsequently divided randomly into two groups for training and validation purposes. The second database included 16 gully-related conditioning factors (GRCFs) that were identified using various resources, including the physiographic characteristics of the study area, literature review and a multi-collinearity test. The list of conditioning factors of GESM is listed in Table 2. The characteristics were later classified into four classes, namely, topographical, geological, hydrological and environmental. The second stage analysed the spatial relationship between the GRCFs and the location of the gullies using FR (Fanos and Pradhan, 2019; Pradhan and Lee, 2010) and NRF models. The third stage determined the importance of the GRCFs in gully occurrence in the study area by using BLR and BRT. The fourth stage involved the preparation of gully erosion susceptibility mapping (GESM) using FR, NRF, BLR and BR models (Youssef et al., 2015). The fifth stage validated the models based on area under receiver operating characteristic (AUROC, the curves) and seed cell area index (SCAI) metrics (Rahmati et al., 2019). 2.2.1. Gully erosion data In assessing the sensitivity of natural events, past records of events should be used as the basis for work and guidance in the modelling process (Van Westen et al., 2008). The more data is collected in relation to the location of the gullies, the higher the accuracy and quality; therefore, the prepared sensitivity map would have increased predictive accuracy (Conoscenti et al., 2018). The locations of the gullies (362 in total) are shown in Fig. 1. These areas were randomly spilt into two groups as follows: 70% for training (253 gully locations) and 30% for validation (109 gully location). The locations of the gullies were determined using different resources. We initially examined the gully inventory compiled by the Isfahan Agricultural and Natural Resources, Research and Education Centre (http://esfahan.areeo.ac.ir/). In the second step, we further validated the inventory using Google Earth images and field surveys. During the latter, we used Global Positioning Systems (GPS) measures to correct the incorrect gully positions. Fig. 3 shows the head cut-off of some of the surveyed gullies, where the main incision is several meters deep. Fig. 3 shows a sample of the mapped gullies in the study area. 2.2.2. Conditioning factors Geo-environmental features can be used as GRCFs to predict gully occurrence (Arabameri et al., 2018c). In this research, we selected 16 conditioning factors (Fig. 4) based on study area characteristics, analysis scale and a multicollinearity test to explain the spatial distribution of gullies. Part of the factors were derived by the Advanced Land Observing Satellite digital elevation model (DEM) with a 30 m resolution (http://www.eorc.jaxa.jp/ALOS/en/aw3d30). In addition, we used remotely sensed and thematic properties. Overall, we considered 1)

2. Materials and methods 2.1. Study area As shown in Fig. 1, Bayazeh Watershed has an area of 5179.34 km2 and is located at 54°43′24″–55°47′12″ E longitude to 2

Catena 183 (2019) 104223

A. Arabameri, et al.

Fig. 1. Study area. a) Location of study area in Iran, b) location of study area in Yazd, Southern Khorasen, and Isfahan Provinces, c) location of training and validation gullies in the study area. Table 1 Lithology of the study area. Group

Unit

Description

Age

A

E1c Ed.avs Egr Jub Jugr Jurb K2lm K2l2 Kdzsh K2m,l Ktl Kns Kbsh K1c K2d.asv Murm mb Pj pCmt2 Plc pCgn P Pel Pz pCdi Qft2 Qsf Qs,d Qft1 Qsl TRJs

Pale-red, polygenic conglomerate and sandstone Dacitic to Andesitic volcano sediment Granite Sandstone, siltstone, Pectinid limestone, marl, gypsum (Bidou series) Upper Jurassic granite including Shir Kuh Granite and Shah Kuh Granite Sandstone, siltstone, and fine grained conglomerate (Garedu red beds) Pale - red marl, gypsiferous marl and limestone Thick - bedded to massive limestone (maastrichtian) Marl, shale, sandstone and limestone (Darreh - Zanjir Fm) Marl, shale and detritic limestone Thin to medium bedded argillaceous limestone and thick bedded to massive, grey orbitolina bearing limestone (Taft Fm) Red sandstone and conglomeratic sandstone Dark grey slightly phyllitized shale with intercalations of sandstone and limestone (Biabanak Shale) Red conglomerate and sandstone Dacitic to andesitic subvolcanic rocks Ligth - red to brown marl and gypsiferous marl with sandstone intercalations Marble Massive to thick - bedded, dark - grey, partly reef type limestone and a thick yellow dolomite band in the upper part (JAMAL FOR) Low - grade, regional metamorphic rocks (Green Schist Facies) Polymictic conglomerate and sandstone Gneiss, granite gneiss and locally including migmatite Undifferentiated Permian rocks Medium to thick - bedded limestone Undifferentiated lower Paleozoic rocks Precambrian diorite Low level piedmont fan and valley terrace deposits Salt flat Unconsolidated windblown sand deposits including sand dunes High level piedmont fan and valley terrace deposits salt lake Dark grey shale and sandstone (SHEMSHAK FM)

Paleocene-Eocene Eocene Eocene Late Jurassic Late Jurassic Late Jurassic Late Cretaceous Late Cretaceous Cretaceous Late Cretaceous Early Cretaceous Early Cretaceous Cretaceous Cretaceous Late Cretaceous Miocene Triassic Permian Pre-Cambrian Pliocene Pre-Cambrian Permian Paleocene-Eocene Early Paleozoic Pre-Cambrian Quaternary Quaternary Quaternary Quaternary Quaternary Triassic-Jurassic

B

C

D E

F

G

3

Catena 183 (2019) 104223

A. Arabameri, et al.

Fig. 2. Flowchart of research.

The NDVI was computed from a LANDSAT-8 image archived by USGS (https://earthexplorer.usgs.gov/). ST, LIT and LU/LC derived from three thematic maps, namely, soil map of Isfahan (http://esfahan. areeo.ac.ir/) at 1:100,000 scale, geological map (http://www.gsi.ir/) at 1:100,000 scale and land use map (https://www.scwmri.ac.ir) at 1:100,000 scale. A set of 453 ground control points were used to validate the LU/LC map. The Kappa coefficient of the produced map was 0.945 which is acceptable. Eqs. (1)–(4) were used for the calculation of LS, TWI, SPI and DD (Moore et al., 1991; Horton, 1932; Rouse et al., 1974).

elevation; 2) slope and 3) aspect (ASP, Zevenbergen and Thorne, 1987; 4) plan curvature (PLC, Heedegen and Beran, 1982; 5) stream power index (SPI, Moore et al., 1991; 6) topographic wetness index (TWI, Beven and Kirkby, 1979; 7) LS (Moore et al., 1991; 8) convergence index (CI, Köthe and Lehmeier, 1993; 9) drainage density (DD, Horton, 1945); distance to 10) streams (DtS), 11) roads (DtR) and 12) faults (DtF); 13) normalised difference vegetation index (NDVI, Rouse et al., 1974); 14) soil type (ST); 15) lithology (LIT) and 16) LU/LC (Fig. 4). DtS, DtR and DtF were computed as Euclidean distances from each pixel centroid to the respective nearest polyline (Arabameri et al., 2018a). The vector-based information was then rasterised to enable it to coincide with the DEM resolution. Road and fault factor maps with scales of 1:50,000 and 1:100,000 were obtained from the National Geographic Organization of Iran (www.ngo-org.ir) and the Geological Society of Iran (http://www.gsi.ir/), respectively. A factor map of streams was extracted from DEM in ArcHydrov10.4 and ArcGISv10.4. 4

0.6 Sinβ ⎞1.3 A LS = ⎛ S ⎞ × ⎛ , 0.0896 ⎝ 22.13 ⎠ ⎝ ⎠

(1)

TWI = In (AS/tanβ),

(2)

SPI = AS × tanβ,

(3)

Catena 183 (2019) 104223

A. Arabameri, et al.

Table 2 List of factors used for GESM. Factor

Range

Source

Resolution

Classes

Method

1. (< 822 m), 2. (822 m – 1004 m), 3. (1004 m – 1209 m), 4. (1209 m – 1477 m), 5. (> 1477 m) 1. (< 5°), 2- (5–10°), 3. (10–15°), 4. (15–20°), 5. (20–30°), 6. (> 5°) 1. Flat (−1°), 2. North (337.5–360°, 0–22.5°), 3. Northeast (22.5–67.5°), 4. East (67.5–112.5°), 5. Southeast (112.5–157.5°), 6. South (157.5–202.5°), 7. Southwest (202.5–247.5°), 8. West (247.4–292.5°), and 9. Northwest (292.5–337.5°) 1. Concave (< 0), 2. Flat (0), 3. Convex (> 0) 1. (< 8.46), 2. (8.46–10.14), 3. (10.14–12.11), 4. (12.11–14.88), 5. (> 14.88) 1. (< 5.59), 2. (5.59–8.02), 3. (8.02–11.56), 4. (> 11.56) 1. (< 7.67 m), 2. (7.67–11.09 m), 3. (11.09–15.12 m), 4. (15.12–20.77 m), 5. (> 20.77 m) 1. (< −38.03), 2. (−38.03 to -12.15), 3. (−12.15–11.37), 4. (11.37–38.03), 5. (> 38.03) 1. (< 0.81), 2. (0.81–1.3), 3. (1.3–1.88), 4. (> 1.88)

Natural break

Min

Max

Elevation (m)

644

2291

PALSAR DEM

30 m

Slope (degree) Aspect

0.00° −1

72.70° 360

PALSAR DEM PALSAR DEM

30 m 30 m

Plan curvature SPI TWI LS

−9.29 6.27 1.10 3.64

11.08 24.89 23.15 55.03

PALSAR PALSAR PALSAR PALSAR

DEM DEM DEM DEM

30 m 30 m 30 m 30 m

Convergence

−100

100

PALSAR DEM

30 m

Drainage density (km/k2) Distance to river (m) Distance to road (m)

0.031

3.11

PALSAR DEM

30 m

0 0

2336.3 41,640

PALSAR DEM Topographic map

30 m 1: 50,000

Distance to fault (m)

0

44,668

Geological map

1: 100,000

NDVI Soil type

0.045 –

0.43 –

Landsat-8 image Soil map

30 m 1: 100,000

Lithology LU/LC

– –

– –

Geological map land use map

1: 100,000 1: 100,000

1. (< 100 m), 2. (100–200 m), 3. (200–300 m), 4. (300–400), 5. (> 400 m) 1. (< 500 m), 2. (500–1000 m), 3. (1000–1500 m), 4. (1500–2000 m), 5. (> 2000 m) 1. (< 500 m), 2. (500–1000 m), 3. (1000–1500 m), 4. (1500–2000 m), 5. (> 2000 m) 1. (< 0.036), 2. (0.036–0.064), 3. (> 0.064) 1. (Dune Lands), 2. (Playa), 3. (Rocky Lands-Aridisols), 4. (Salt Flats), 5. (Aridisols), 6. (Entisols/Aridisols) 1. (A), 2. (B), 3. (C), 4. (D), 5. (E), 6. (F), 7. (G) 1. (Agriculture), 2. (bareland), 3. (Saltland-poorrange), 4. (sanddune), 5. (wetland), 6. (midrange), 7. (poorrange), 8. (saltlake), 9. (saltland), 10. (sanddune), 11. (Urban), 12. (wetland), 13. (woodland)

Fig. 3. Examples of mapped gullies in the study area. 5

Natural break Equal interval

Natural Natural Natural Natural

break break break break

Natural break Natural break Natural break Natural break Natural break

Supervised classification Lithological units Supervised classification

Catena 183 (2019) 104223

A. Arabameri, et al.

Fig. 4. Gully erosion conditioning factors shown in map form. From top left to bottom right: Elevation, Slope, Aspect, Plan Curvature, Stream Power Index, Topographic Wetness Index, Slope Length and Slope Steepness Factor, Convergence Index, Drainage Density, Distance to Stream, Roads and Faults, NDVI, Soil, Lithology and Land Use maps. Each factor is represented by several classes, being obtained from reclassifying the original continuous information with the exception of the last three which were natively categorical. n

DD =

∑i = 1 Si a

,

2.2.3. Predictive models 2.2.3.1. Natural risk factor (NRF). NRF is a heuristic approach that relies on a subjective weighting procedure. For each factor, the proportion of (gully) occurrences in a reference class with respect to the average occurrences (gully) in all classes must be computed (Gupta and Joshi, 1990). The ratio between the two occurrence values is then used to assign certain weights to each pixel, where the reference class is observed. These weights can have values of 0 for classes with proportion < 1, 1 for classes with 1 < proportion < 2 and 2 for classes with proportion > 2. Once every pixel is assigned a set of weights, each one corresponding to a conditioning factor, a susceptibility map can be obtained by summing all values for each pixel. This method is clearly subjective for the weight values that are arbitrarily decided. Despite this evident flaw, NRF has recently been providing good results in predicting geomorphological processes over

(4)

where AS is the specific catchment area of the basin (m2/m), β is slope n steepness (degrees), ∑i = 1 Si is the total length of drainages in (km) and A is area of the drainage watershed (km2). The LR and BRT can handle linear and categorical covariates; by contrast, FR and NRF can use categorical properties only. Therefore, ST, LIT and LU/LC, which were categorical by nature, did not undergo any transformation. Meanwhile, all other continuous predictors were transformed into their ordinal counterpart. Each new ordinal factor was produced by reclassifying its continuous raster into several classes. Finally, all layers were unified with respect to PALSAR DEM spatial resolution (30 m) and the UTM Zone39N geographic coordinate system.

6

Catena 183 (2019) 104223

A. Arabameri, et al.

possible splits (for it would be an extremely demanding operation) and thus does not guarantee the finding of the absolute best solution. However, the recursive approach forces the tree to test a large number of sub-sub-splits towards a satisfactory performance.

Table 3 Multi-collinearity among conditioning factors. Factors

Multicollinearity TOL

Soil type LU/LC Slope aspect Convergence index Elevation Distance to stream Distance to fault Distance to road Drainage density LS SPI TWI Lithology NDVI Plan curvature Slope Rainfall Profile curvature Soil texture

0.329 0.566 0.907 0.764 0.387 0.740 0.248 0.362 0.527 0.211 0.196 0.148 0.587 0.801 0.676 0.294 0.021 0.018 0.027

VIF

2.2.4. Model building strategy and validation The original 362 gully head cuts represented our target variable. To calibrate our models, we randomly extracted 30% of the head cuts and merged them with an equivalent number of pixels denoting the absence conditions of gullies (Arabameri et al., 2019d, 2019e). These absences were also randomly extracted within the study area. Four predictive equations obtained for NRF, FR, BLR and BRT were tested in predicting the unknown subset of gully presences. This step allowed us to produce four susceptibility maps that depict areas prone to future gully erosion under the assumption that the four sets of causative relations obtained will not change significantly over time. Validation performance was evaluated based on receiver operating characteristics curves and their area under the curve (AUC) (Meliho et al., 2018; Azareh et al., 2019; Arabameri et al., 2019f, 2019h, 2019i, 2019j) and SCAI (Arabameri et al., 2018a, 2018b, 2018c). The values of AUROC varies between 0.5 and 1, and the closer the AUROC to 1, the better the performance of the method (Yesilnacar, 2005). SCAI is the ratio of the percentage area of each of the zoning classes to the percentage of phenomena occurring in each class (Yilmaz et al., 2012). In SCAI, the classification accuracy is acceptable when the values of the susceptibility classes increase from very low to very high (Süzen and Doyuran, 2004).

3.04 1.76 1.10 1.30 2.58 1.35 4.03 2.76 1.89 4.34 4.23 3.91 1.70 1.24 1.47 1.39 39.23 55.09 36.66

space (Mohammadi et al., 2014; Abedini and Tulabi, 2018). 2.2.3.2. Frequency ratio (FR). The FR is calculated as the ratio between (gully) occurrences and (gully) absences within a given conditioning factor class (Rahmati et al., 2016). Larger ratios imply larger contributions towards erosion-prone conditions at pixels where the factor class is assigned (Meliho et al., 2018). This simple case can easily be extended to situations where more than one factor is considered. In this case, the proneness or susceptibility to gully erosion at a given pixel is computed by summing all the ratios obtained for the other predisposing factors (Rahmati et al., 2016). The FR is an empirical approach. It is not classified as statistical because FR does not rely on an underlying statistical distribution (Lombardo and Mai, 2018).

3. Results 3.1. Multicollinearity test (MCT) Performing an MCT among GRCFs is important in preparing GESM (Arabameri et al., 2018a). A linear correlation between the GRCFs would result in a decreased prediction accuracy of the model (Arabameri et al., 2018b). In this study, two parameters, namely, tolerance (TOL) and variance inflation factor (VIF), were used for the linear test. Coefficient of TOL is ≤1 and VIF factor ≥ 10 indicate the existence of linearity between the GRCFs (Chen et al., 2019; Arabameri et al., 2018c). The results of the MCT among the GRCFs (Table 3) showed that among the 16 GRCFs, three conditioning factors, namely, rainfall, profile curvature and soil texture, with TOL values of 0.021, 0.018 and 0.027 and VIF values of 39.23, 55.09 and 36.66, respectively, had collinearity. They would decrease the accuracy of the model and were hence eliminated. The modelling was conducted using 16 factors, namely, elevation, slope, ASP, PLC, SPI, TWI, LS, CI, DD, DtS, DtR, DtF, NDVI, ST, LIT and LU/LC.

2.2.3.3. Binary logistic regression (BLR). BLR is a technique that is used in statistics and machine learning (Dewitte et al., 2015) or in numerous classification problems. This technique allows the modelling of dependent variables that are constrained to assume binary values. In contrast to classic linear regression, LR fits a logistic function to the data, where the function domain is bounded between 0 and 1. The function itself represents the probability of presence or absence given a set of covariates. The LR technique can be used in various cases. In this manuscript, we referred to BLR to represent the implementation of this technique in a generalised linear model (GLM) framework. Whereas FR does not rely on a statistical distribution, the BLR we used here assumed that the gullies occurred over space according to a Bernoulli distribution.

3.2. Application of FR and NRF models In this study, the spatial relationship between the GRCFs (Fig. 4) and the gully locations was analysed using FR and NRF models (Table 4). Based on the results (Table 4) of the elevation factor, areas with 822–1004 m elevation with FR = 2.32 and NRF = 2.8 had the greatest effect on gully occurrence in the study area but decreased in importance with increasing elevation. Hence, areas with the highest elevation (FR = 0 and NRF = 0) had the least effect on the gully occurrence. These results are consistent with those of Nazari Samani et al. (2009) and Arabameri et al. (2018b). Slope analysis using FR and NRF methods showed that classes of 5°–10° (FR = 1.99) and > 5° (FR = 1) in the FR method and classes of > 5 (NRF = 4.74) and 5–10 (NRF = 1.02) in the NRF method had the most correlation with the gully occurrence in the study area. In the two methods, an increase in slope decreased the correlation, given that surface runoff concentration and accumulation occur more in areas with low slopes than in places with high slopes. This finding is in line with that of Rahmati et al. (2016). Meanwhile, based on the ASP factor, areas with eastern

2.2.3.4. Boosted regression trees (BRT). BRT is a statistical method whereby a regression tree is used to explain discrete data (in which case it is denoted as a classification tree) (Breiman, 1998). BRT essentially aims at reconstructing the logistic function by splitting the predictors' domain and using the observations contained in each split (Therneau et al., 2014). The splitting procedure is performed recursively, and at each iteration, BRT tries to find the best split that minimises the error or loss function (Lombardo et al., 2015). This procedure terminates when the gain in reducing the error becomes negligible. At this point, each predictor is already divided into an optimal set of sub-classes, each one characterised by a specific rule (or, intuitively, ‘correlation’) with respect to the presence or absence of the process (gully) being modelled (Arabameri et al., 2018a). Interestingly, BRT is a ‘greedy’ model. In other words, BRT does not explore all the 7

Catena 183 (2019) 104223

A. Arabameri, et al.

Table 4 Spatial relation between gully erosion conditioning factors and gully locations using frequency ratio and LNRF. Factors

Elevation (m)

Slope (°)

Aspect

Curvature (100/m)

SPI (100/m)

TWI (100/m)

LS (100/m)

Convergence index (100/m)

Drainage density (km/km2)

Distance to stream (m)

Distance to road (m)

Distance to fault (m)

NDVI

Soil type

Classes

< 822 822–1004 1004–1209 1209–1477 > 1477 <5 5–10 10–15 15–20 20–30 > 30 F N NE E SE S SW W NW Concave Flat Convex < 8.46 8.46–10.14 10.14–12.11 12.11–14.88 > 14.88 < 5.59 5.59–8.02 8.02–11.56 > 11.56 < 7.67 7.67–11.09 11.09–15.12 15.12–20.77 > 20.77 < −38.03 −38.03 to −12.15 −12.15–11.37 11.37–38.03 > 38.03 < 0.81 0.81–1.3 1.3–1.88 > 1.88 < 100 100–200 200–300 300–400 > 400 < 500 500–1000 1000–1500 1500–2000 > 2000 < 500 500–1000 1000–1500 1500–2000 > 2000 < 0.036 0.036–0.064 > 0.064 Dune Lands Playa Rocky Lands-Aridisols Salt Flats Aridisols Entisols/Aridisols

Pixels in domain

Gully pixels

FR

No

%

No

%

2,553,993 1,385,671 939,836 683,904 192,919 4,553,985 492,780 229,940 153,667 195,627 130,324 856,168 1,301,109 660,265 475,926 545,910 450,903 433,990 464,222 567,829 1,839,028 2,033,078 1,884,215 1,744,947 1,853,109 1,351,701 616,428 190,138 1,524,577 2,662,377 1,205,062 364,307 1,825,323 1,901,468 1,229,891 610,621 189,020 623,362 1,296,345 1,914,404 1,297,319 620,399 1,279,836 2,329,890 1,690,059 456,538 1,324,206 1,036,330 913,801 625,674 1,856,312 205,295 190,368 182,048 176,457 5,002,155 270,009 261,075 231,726 212,200 4,781,313 1,496,061 3,148,317 1,111,934 743,296 297,794 1,891,604 862,219 1,883,409 78,001

44.37 24.07 16.33 11.88 3.35 79.11 8.56 3.99 2.67 3.40 2.26 14.87 22.60 11.47 8.27 9.48 7.83 7.54 8.06 9.86 31.95 35.32 32.73 30.31 32.19 23.48 10.71 3.30 26.49 46.25 20.93 6.33 31.71 33.03 21.37 10.61 3.28 10.84 22.54 33.28 22.55 10.79 22.23 40.48 29.36 7.93 23.00 18.00 15.87 10.87 32.25 3.57 3.31 3.16 3.07 86.90 4.69 4.54 4.03 3.69 83.06 25.99 54.69 19.32 12.91 5.17 32.86 14.98 32.72 1.36

47 141 44 20 0 199 43 8 2 0 0 26 43 35 51 26 19 12 13 27 81 90 81 64 101 50 26 11 74 123 35 20 60 100 53 31 8 20 56 111 49 16 13 166 67 6 81 30 50 18 73 44 36 36 22 114 18 14 8 4 208 33 169 50 0 8 116 19 109 0

18.65 55.95 17.46 7.94 0.00 78.97 17.06 3.17 0.79 0.00 0.00 10.32 17.06 13.89 20.24 10.32 7.54 4.76 5.16 10.71 32.14 35.71 32.14 25.40 40.08 19.84 10.32 4.37 29.37 48.81 13.89 7.94 23.81 39.68 21.03 12.30 3.17 7.94 22.22 44.05 19.44 6.35 5.16 65.87 26.59 2.38 32.14 11.90 19.84 7.14 28.97 17.46 14.29 14.29 8.73 45.24 7.14 5.56 3.17 1.59 82.54 13.10 67.06 19.84 0.00 3.17 46.03 7.54 43.25 0.00

0.42 2.32 1.07 0.67 0.00 1.00 1.99 0.79 0.30 0.00 0.00 0.69 0.75 1.21 2.45 1.09 0.96 0.63 0.64 1.09 1.01 1.1 0.98 0.84 1.24 0.84 0.96 1.32 1.11 1.06 0.66 1.25 0.75 1.20 0.98 1.16 0.97 0.73 0.99 1.32 0.86 0.59 0.23 1.63 0.91 0.30 1.40 0.66 1.25 0.66 0.90 4.90 4.32 4.52 2.85 0.52 1.52 1.22 0.79 0.43 0.99 0.50 1.23 1.03 0.00 0.61 1.40 0.50 1.32 0.00

NRF NRF

W

0.93 2.80 0.87 0.40 0 4.74 1.02 0.19 0.05 0.00 0.00 1.03 1.71 1.39 2.02 1.03 0.75 0.48 0.52 1.07 0.96 1.07 0.96 1.27 2.00 0.99 0.52 0.22 1.17 1.95 0.56 0.32 1.19 1.98 1.05 0.62 0.16 0.40 1.11 2.20 0.97 0.32 0.21 2.63 1.06 0.10 1.61 0.60 0.99 0.36 1.45 0.87 0.71 0.71 0.44 2.26 0.36 0.28 0.16 0.08 4.13 0.39 2.01 0.60 0.00 0.19 2.76 0.45 2.60 0.00

0 2 0 0 0 2 1 0 0 0 1 1 1 2 1 0 0 0 1 0 1 0 1 2 0 0 0 1 1 0 0 1 1 1 0 0 0 1 2 0 0 0 2 1 0 1 0 0 0 1 0 0 0 0 2 0 0 0 0 2 0 2 0 0 0 2 0 2

(continued on next page)

8

Catena 183 (2019) 104223

A. Arabameri, et al.

Table 4 (continued) Factors

Lithology

Land use/land cover

Classes

A B C D E F G Agriculture Bare land Salt land - poor range Sand dune Wetland Mod range Poor range Salt lake Salt land Sand dune Urban Wetland Woodland

Pixels in domain

Gully pixels

FR

No

%

No

%

58,984 110,211 1,220,702 214,436 556,507 3,563,985 31,498 9400 469,019 782,134 749,165 50,152 724,928 2,151,961 216,998 97,194 32 633 141,198 363,506

1.02 1.91 21.21 3.73 9.67 61.91 0.55 0.16 8.15 13.59 13.01 0.87 12.59 37.38 3.77 1.69 0.00 0.01 2.45 6.31

0 0 170 1 6 75 0 0 10 31 0 0 9 184 0 11 0 0 0 2

0.00 0.00 67.46 0.40 2.38 29.76 0.00 0.00 4.05 12.55 0.00 0.00 3.64 74.49 0.00 4.45 0.00 0.00 0.00 0.81

0.00 0.00 3.18 0.11 0.25 0.48 0.00 0.00 0.50 0.92 0.00 0.00 0.29 1.99 0.00 2.64 0.00 0.00 0.00 0.13

NRF NRF

W

0.00 0.00 4.72 0.03 0.17 2.08 0.00 0.00 0.53 1.63 0.00 0.00 0.47 9.68 0.00 0.58 0.00 0.00 0.00 0.11

0 0 2 0 0 2 0 0 0 1 0 0 0 0 0 0 0 0 0 1

most sensitive to gully erosion. This result is in line with that of Arabameri et al. (2018b). With regard to the LIT factor, class of C with FR (3.18) and NRF (4.72) had the most importance in gully occurrence. Based on the LU/LC factor, in the FR model, salt land (2.64) and in NRF model, poor-range (9.68) classes had the most important role in gully occurrence. After the relationship between GECFs and gully locations in the study area and the weight of each class of GRCFs were determined using FR and NRF models, the weights of each of the classes were added to the maps in the software ArcGIS10.5. Then, with use of the weighted sum tool, the GESM using FR (Eq. (6)) and NRF (Eq. (7)) was obtained.

orientation (FR = 2.45 and NRF = 2.02) had greater impact on gully occurrence in the study area compared with places with orientations in other geographic directions. According to the PLC factor, flat areas (FR = 1.1 and NRF = 1.07) had a strong relationship with gully occurrence, followed by concave and convex areas. This result is in line with those of Wilkinson and Humphreys (2006), Conforti et al. (2011), Conoscenti et al. (2014), Rahmati et al. (2017b) and Arabameri et al. (2018c). Based on the SPI, TWI, LS and CI factors, FR model classes of > 14.88 (FR = 1.24), > 1.25 (FR = 1.25), 7.67–11.09 m (FR = 1.2) and −12.15–11.37 (FR = 1.32) and NRF model classes of 8.46–10.14 (NRF = 2), 5.59–8.02 (NRF = 1.95), 7.67–11.09 m (NRF = 1.98) and −12.15–11.37 (NRF = 2.2) had the most important values in gully occurrence in the study area. In the case of DD, DtS, DtR and DtF in the FR model, classes of 0.81–1.3 km/km2 (FR = 1.63), < 100 m (FR = 1.40), < 500 m (FR = 4.9) and < 500 m (FR = 1.52) and classes of < 81 km/km2 (FR = 0.23), 300–400 m (FR = 0.66), > 2000 (FR = 0.52) and 1500–2000 (FR = 0.43) had the most and least importance in gullying, respectively. The results showed that the areas with the least distances with the streams, roads and faults were most sensitive to gully erosion, which indicated the considerable impact of linear parameters on gully occurrence. This result is in line with those in several studies (Croke and Mockler, 2001; Croke and Hairsine, 2006; Dube et al., 2014; Conoscenti et al., 2014; Rahmati et al., 2016; Arabameri et al., 2018b, Meliho et al., 2018). Vandekerckhove et al. (2000) stated that drainage network is a key factor in determining the location and trajectory of gully erosion. Meanwhile, roads generate impervious surface and runoff accumulation and thus disturb the natural process of drainage networks; therefore, this factor played a key role in gully occurrence (Croke and Hairsine, 2006; Rahmati et al., 2017b). By contrast, in the NRF model, classes of 0.81–1.3 (NRF = 2.63), < 100 m (NRF = 1.61), > 2000 m (NRF = 2.26) and < 2000 m (NRF = 4.13) had the most importance in gully occurrence. With regard to NDVI, previous research (Conforti et al., 2011; Dewitte et al., 2015; Pourghasemi et al., 2017; Zabihi et al., 2018) has indicated that vegetation substantially affects soil conservation, and areas with the lowest vegetation cover have the highest sensitivity to gully erosion. However, according to the results of this research, the areas with the lowest vegetation (< 0.036 class with FR = 0.50 and NRF = 0.39) had least sensitivity to gully erosion. The reason is that areas with low vegetation were covered with sandy hills. Meanwhile, according to the ST factor, rocky aridisol areas with FR = 1.40 and NRF = 2.76 were the

GESMFR = (NDVIFR) + (STFR) + (SPIFR) + (LITFR) + (DDFR) + (LU. LCFR) + (ASPFR) + (DtFFR) + (DtRFR) + (DtSFR) + (ElevationFR) + (CIFR) + (LSFR) + (SlopeFR) + (PLCFR) + (TWIFR)

(6)

GESMNRF = (NDVINRF ) + (STNRF) + (SPINRF) + (LITNRF) + (DDNRF ) + (LU. LCNRF) + (ASPNRF) + (DtFNRF) + (DtRNRF) + (DtSNRF) + (ElevationNRF) + (CINRF) + (LSNRF) + (SlopeNRF) + (PLCNRF) + (TWINRF)

(7)

Values resulting from GESMs using FR (9.90–31.01) and NRF (9.39–52.11) based on natural break method were divided into five classes of very low (9.90–14.20 in FR and 9.39–23.96 in NRF), low (14.20–16.69 in FR and 23.96–28.65 in NRF), moderate (16.69–19.25 in FR and 28.65–33.85 in NRF), high (19.25–22.23 in FR and 33.85–39.54 in NRF) and very high (22.23–31.04 in FR and 39.54–52.11 in NRF) (Fig. 5a and 5b). Results of susceptibility classes (Table 5) indicated that in FR (45.07%, 25.18%, 13.75%, 13.03% and 2.98%) and in NRF (11.18%, 28.49%, 23.95%, 19.98% and 16.39%) areas were located in very low, low, moderate, high and very high susceptibility classes. By contrast, in FR (9.75%, 12.57%, 7.5%, 24.43% and 45.75%) and in NRF (3.88%, 6.31%, 15.43%, 13.29% and 61.09%) of gullies were located in very low to very high susceptibility classes. 3.3. Application of BLR and BRT models The results of determining the importance of GRCFs in the gully occurrence using BLR (Fig. 6) and BRT (Fig. 7) methods showed that in 9

Catena 183 (2019) 104223

A. Arabameri, et al.

Fig. 5. Gully erosion susceptibility mapping: a) Frequency ratio (RF), b) Numerical Risk Factor (NRF), c) Binary Logistic regression (BLR), and d) Boosted regression tree (BRT).

Table 5 Gully erosion susceptibility classes and values of seed cell area index. Models

FR

LNRF

BLR

BRT

Value

Very low Low Moderate High Very high Very low Low Moderate High Very high very low Low Moderate High Very high Very low Low Moderate High Very high

Total area

Training gully

Validation gully

Count

%

Count

%

Count

%

2,593,488 1,449,064 791,083 749,992 171,262 643,471 1,639,488 1,378,449 1,149,933 943,487 4,663,920 49,863 441,798 498,424 102,317 1,283,630 1,616,464 1,451,204 1,009,102 394,432

45.07 25.18 13.75 13.03 2.98 11.18 28.49 23.95 19.98 16.39 81.02 0.87 7.68 8.66 1.78 22.31 28.09 25.22 17.53 6.85

26 31 17 63 115 8 11 50 30 153 79 3 25 52 93 4 34 19 46 149

10.32 12.30 6.75 25.00 45.63 3.17 4.37 19.84 11.90 60.71 31.35 1.19 9.92 20.63 36.90 1.59 13.49 7.54 18.25 59.13

10 14 9 26 50 5 9 12 16 67 31 1 13 28 36 1 14 10 19 65

9.17 12.84 8.26 23.85 45.87 4.59 8.26 11.01 14.68 61.47 28.44 0.92 11.93 25.69 33.03 0.92 12.84 9.17 17.43 59.63

10

Total gully

SCAI

9.75 12.57 7.50 24.43 45.75 3.88 6.31 15.43 13.29 61.09 29.89 1.05 10.92 23.16 34.97 1.25 13.17 8.36 17.84 59.38

4.624 2.003 1.832 0.534 0.065 2.881 4.514 1.553 1.503 0.268 2.710 0.822 0.703 0.374 0.051 17.810 3.133 2.017 0.983 0.115

Catena 183 (2019) 104223

Conditioning factor

A. Arabameri, et al.

TWI Plan curvature Slope LS Aspect Convergence index Distance to stream Elevation Distance to road Distance to fault LU/LC Drainage density Lithology SPI Soil type NDVI

0.00 1.07 1.85 2.47 2.92

2.92 2.92 2.92 2.92 2.93 4.68

5.10 6.95 7.19 12.08 41.10

0

10

20

%

30

40

50

Conditioning factor

Fig. 6. Relative importance of conditioning factors using Binary logistic regression.

Soil type Plan curvature LS SPI Aspect Convergence index TWI Slope Drainage density Distance to stream LU/LC Distance to fault Elevation Distance to road NDVI Lithology

0.20 1.16 1.17 1.91 2.72 3.93 4.03 4.13 5.24 5.34 5.44

Fig. 8. Validation of results using area under curve (AUC).

resulting from the BRT was is classified based on natural break method and divided into five classes (Fig. 5d) from very low to very high. Based on the BRT results (Table 5), 22.31%, 28.09%, 25.22%, 17.53% and 6.85% of the study area were in very-low- to very-high-susceptibility classes.

7.75 9.77 10.98 13.80

3.4. Validation of results

22.45

0

5

10

%

15

20

25

The validation results of the four different models using the AUROC method are shown in Fig. 8 and using the SCAI indicator are shown in Table 5. Based on AUROC results, BRT and FR models with AUROC = 0.834 and 0.823 (83.4% and 82.3%) had higher prediction accuracy than NRF and BLR with AUROC = 0.746 and 0.659 (74.6% and 65.9%). Results of SCAI indicated that BRT, FR and BLR models have acceptable classification accuracy. The values of susceptibility classes decreased from very low to very high in these models but a similar trend cannot be seen in the NRF model.

Fig. 7. Relative importance of conditioning factors using boosted regression tree.

BLR, factors of NDVI (41.10%), ST (12.08%) and SPI (7.19%) and in BRT, factors of LIT (22.45%), NDVI (13.8%) and DtR (10.98%) had the largest impact on the gully occurrence in the study area. By contrast, in BLR, factors of TWI (0), PLC (1.07%) and slope (1.85%) and in BRT, factors of ST (0.201%), PLC (1.15%) and LS (1.16%) had the least impact on gully occurrence in the study area. BLR factors of LIT, DD, LU/LC, DtF, DtR, elevation, DtS, CI, ASP, and LS and BRT model factors of elevation, DtF, LU/LC, DtS, DD, slope, TWI, CI, ASP, and SPI ranked 4th to 13th, respectively. Eq. (8) was used for the calculation of GESM using the BLR model.

4. Discussion Soil erosion by water every year imposes extensive damage to humankind around the world (Arabameri et al., 2018b). Water erosion is commonly divided into four groups comprising sheet, splash or raindrop, rill and gully (Gong et al., 2011). Gully erosion is the most destructive form and can be exacerbated by factors that are dependent on human activity and change in its ecosystem equilibrium (Dymond et al., 2016; Goodwin et al., 2017). This kind of erosion can cause much damage to land, soil, ecosystems and human infrastructure (Pourghasemi et al., 2017). This phenomenon, by shortening the relationship between upstream and downstream areas, exacerbates surface runoff and decrease groundwater nutrition (Poesen et al., 2003). The results of the field data of the morphometric characteristics of the gullies in the study area showed that the longitudinal extension of the gullies is high and sometimes reaches several hundred meters. This result is more due to the backwardness of the forehead of the gully, which indicates the intense activity of gullying in the study area. The lateral and deep development of the gullies in the study area is due to lateral erosion, which is caused by the undercutting slope bank of the gullies and collapse caused by gravity. The gullies are mainly V and U-

GESMBLR = −16.20 + (8.82×NDVI) + (1.97 × ST) + (0.92 × SPI) + (0.86 + LIT) + (0.46 × DD) + (0.38 × LU. LC) + (0.001 × ASP) + (0.00007 × DtF) + (−0.0001 × DtR) + ( −0.0005 × DtS) + (−0.0021×Elevation) + ( −0.0023 × CI) + ( −0.1 × LS) + (−0.23×Slope) + (− 0.4 × PLC) + (− 0.63 × TWI)

(8)

Values of the resulted GESM by using BLR model ranged from 0 to 0.005. These values are based on natural break method and divided into five classes (Fig. 5c) of very low (0–0.00004), low (0.000044–0.000048), moderate (0.000048–0.000092), high (0.000092–0.00057) and very high (0.00057–0.0057). Based on the results, 81.02%, 0.87%, 7.68%, 8.66% and 1.78% of the study area was located in very-low- to very-high-susceptibility classes. The map 11

Catena 183 (2019) 104223

A. Arabameri, et al.

methods showed that the ranking and significance of the factors are different in each method, which indicates uncertainty in the data and the structural difference of the used models. Validation results showed that BRT and FR performed better than the NRF and BLR models. BRT and FR methods have been used in various environmental management fields of floods, landslides, underground waters and gullies and have shown high predictive and zoning capabilities (Arabameri et al., 2019c; Lee and Pradhan, 2007; Ozdemir, 2011; Naghibi et al., 2016; Chen et al., 2019). Naghibi et al. (2016) tested 3 machine learning models of BRT, CART and RF for groundwater potential mapping in Iran using 13 factors of slope, ASP, elevation, TWI, LS, PLC, profile curvature, DtS, DtF, LIT, LU/LC, DD and fault density. The BRT model had the best prediction accuracy. Ozdemir (2011) tested FR, weights of evidence (WoE) and BLR methods for groundwater spring potential mapping using 18 conditioning factors in the Sultan Mountains (Konya, Turkey) and stated FR and WoE models were relatively good estimators, whereas BLR was a relatively poor estimator. Lee and Pradhan (2007) tested FR and BLR for landslide hazard mapping at Selangor, Malaysia using nine geo-environmental factors, including slope, ASP, PLC, DtS, LIT, distance to lineaments, LU/ LC, NDVI and precipitation distribution and stated that FR is better in prediction than BLR. However, in some cases (Meten et al., 2016; Rasyid et al., 2016), BLR performed better than FR. Overall, one cannot conclusively state that BRT and FR should always be preferred to the other two methods we tested. However, we consistently used the same dataset among the four algorithms; the only differences are linked to the structure of the model itself. Therefore, we can assume that the greater classification capacity of BRT is responsible for the high performance. The FR case serves instead as a reminder that simple models can sometimes suffice and a greater complexity does not always lead to greater performance. As for BLR, the unsatisfying results we achieved are quite surprising. However, it is important to note that in order to respect the computational requirements of FR and NRF and to keep the same dataset structure across models; we reclassified the conditioning factors into a number of classes. Moreover, the number of gullies we modelled was not extremely large. Therefore, the BRT is a classifier, which creates simple binary rules in a tree structure, thus it can be adaptive even with a small dataset. The FR and NRF are two simple models, which essentially only rely on gully frequencies within factor classes. Therefore, even small datasets can still provide informative frequencies and associated relations to gully erosion. Conversely, BLR is a multivariate statistical approach and the regression procedure can be affected by limited information. This may be especially the case in the present contribution where we distributed/ nonlinearized the gully information into a number of factor classes rather than keeping it together in a single vector to be used as a linear covariate.

shaped cross-sections with retrograde and backstop with steep and unstable walls. Gullies around the Bayazeh Village has deepened buy have little width due to high density of loess deposits. The high density of ditches is mainly around the Bayazeh village. Other types of gullies are created on river terraces and on slopes leading to main networks, which are of frontal gullies. Conversely, retrogressive erosion could destroy agricultural land, rural roads and irrigation canals. A retrograding erosion contributes to the expansion and lengthening of a gully in the direction of the head-cut where the process started rather than the direction towards which the gully commonly evolves. Extreme rainfall, which could fall in arid areas, acts as a triggering factor. Piping plays an important role in the development of gullies. This phenomenon is due to the underground and undercutting of disintegrating and dissolvable materials and is sometimes referred to as tunnel erosion. These tunnels are mainly caused by permeable and dissolvable sediments and possibly due to the concentration of surface runoff. Underground flow causes this phenomenon due to the finegrained formations (silt and clay) and soluble (salt, gypsum and lime) elements. In the study area, in marl and silt formations at low slopes, tunnelling erosion is an important factor in the formation and development of gullies. This phenomenon from the upper part of the gullies is created in the form of holes and then leads to the absorption of surface water. Subsequently, the tunnel begins to expand and eventually, after the fall of the roof, become a new gully and cause the expansion of earlier gullies. The factors causing the formation of piping in the study area are clay and evaporative formations, such as clay, gypsum and salt, which are highly dissolved. The dry climatic condition of the study area results in cracks in the ground. During rainfall season, the water enters these gaps, penetrates the earth and moves towards the direction of a slope by dissolving and washing the colloid soils. Piping is subsequently created. Field observations showed that the gullies were mainly found in areas where the formation of the region is marl and silt and consisted of different parts of the gypsum, limite and salt with high solubility. Piping, tension crack development, dispersion and bank collapse and rill erosion are the most important mechanisms in the initiation and development of gullies in the study area. This result is in line with those of Gómez-Gutiérrez et al. (2015), Ollobarren et al. (2016) and Barnes et al. (2016). Finally, 362 gullies were used as inventory in the study area. Given that various environmental factors are effective in gully erosion, understanding and recognising the relationship among these factors and the occurrence of gullies, and providing a GESM constitute an important strategy for water and soil resource management (Conoscenti et al., 2014; Shit et al., 2015). Gully erosion is a result of the interaction of precipitation, topography, land use, soils and lithology (Conoscenti et al., 2014; Luffman et al., 2015). In this research, we tested four models (FR, NRF, BLE and BRT) for GESM. Each of these models has some disadvantages and advantages. FR and NRF models are unable to determine the importance of factors (Abedini and Tulabi, 2018). Simplicity of implementation and interpretation are the most important advantages of FR and NRF (Hong et al., 2017). By contrast, BLR and BRT require a more rigorous implementation and computational effort but they provide useful estimates on the spatial relationships between a given phenomenon and the factors affecting its occurrence (Dewitte et al., 2015; Arabameri et al., 2018a). The ability to use continuous and categorical variables as independent variables are also, another advantage of BRT and BLR compared to the other two simple model we tested, where continuous properties need to be reclassified prior to the analyses (Akgun and Türk, 2011). Based on the results of BLR, NDVI, ST and SPI and in BRT model, LIT, NDVI and DtR factors have key roles in gully occurrence. This result is in line with those of various studies (Croke and Mockler, 2001; Vanacker et al., 2003; Croke and Hairsine, 2006; McCloskey et al., 2016; Shellberg et al., 2016; Rahmati et al., 2017b; Arabameri et al., 2018a, 2018c). The results of analysing the importance of factors using BLR and BRT

5. Conclusion GESM is the first and most important step in reducing the damage and managing gully erosion leading to waste of water and soil resources. Reducing mapping uncertainty towards an accurate identification of the areas prone to gully erosion can be done by selecting the most performing and robust model. In this study, a comparative evaluation of the FR, NRF, BLR and BRT methods was used to select the best model for GESM. Remote sensing data and GIS were applied as well during pre- and post- processing phases. This paper shows that BRT and FR outperform the NRF and BLR according to both AUROC and SCAI metrics. Being machine learning techniques constantly developed and improved, we assume the performance difference between such methods and statistically-based ones will stay the same if not increase. As a result, we envision similar uses in other regions where gully erosion represent a threat to local communities. Master planners and engineers can use GEMS information to support decision-making administrative and technical duties for land conservation and development 12

Catena 183 (2019) 104223

A. Arabameri, et al.

projects. In fact, our GESM is quite affordable for it does not require substantial computing facilities nor prolonged computational times. And, it converged to satisfactory and reliable predictive results by using a relatively small dataset which in turn, makes our approach even more appealing in data-scarse environments. However, we recall here that in order to ensure the comparison among the four models, we performed some pre-processing steps, which may have specifically promoted some algorithms. For instance, the logistic regression may have suffered from a small number of instances within factor classes, possibly resulting in larger uncertainty. Conversely, a linear use of each factor could have strengthened the regressors' estimates. However, building NRF and FR models would not have been possible under the linear use constraint. Taking aside the predictive component of this research, a better understanding can be drawn on local environmental factors which contribute to gully-erosion-prone conditions. In fact, despite the algorithmic architecture of the four models resulted in a different predictor subset deemed as relevant to promote gullying, the key parameters stands out by being commonly shared among predictive models. This is the case for the vegetation density expressed via the NDVI and the outcropping lithotypes. This is a striking result for there is no agreement among models on the role of terrain properties which are usually among the most relevant contributors. Keeping the local perspective, the GESM for the Bayazeh watershed largely showed areas prone to gullies and erosion. Hence, residential areas and infrastructure, such as roads and power lines, in areas with high susceptibility to erosion should adopt management measures, such as vegetation restoration. In addition, mechanical structures should be put in place to protect the soil from water erosion. These practical measures could reduce and control soil and gully erosion in the study area. Future applications may involve testing the same procedures on a larger dataset, possibly moving from the catchment to the regional or even national scales. This could also shed light on how transferable FR, NRF, BLR and BRT can be in a different geographic setting and associated data structures. This could be a step further in promoting a set of methods to be commonly shared among the community as a reference, irrespective of the scale of application. In turn, this could help to leave behind methods that become outdated as new algorithms are proposed and to converge to more standardized or even meta-analytical procedures like in other branches of science.

Total Environ. 613–614, 1385–1400. Arabameri, A., Pourghasemi, H.R., Yamani, M., 2017b. Applying different scenarios for landslide spatial modeling using computational intelligence methods. Environ. Earth Sci. 76, 832. Arabameri, A., Pradhan, B., Pourghasemi, H.R., Rezaei, K., Kerle, N., 2018a. Spatial modelling of gully erosion using GIS and R programing: a comparison among three data mining algorithms. Appl. Sci. 8 (8), 1369. https://doi.org/10.3390/ app8081369. Arabameri, A., Rezaei, K., Pourghasemi, H.R., Lee, S., Yamani, M., 2018b. GIS-based gully erosion susceptibility mapping: a comparison among three data-driven models and AHP knowledge-based technique. Environ. Earth Sci. 77, 628. https://doi.org/10. 1007/s12665-018-7808-5. Arabameri, A., Pradhan, B., Rezaei, K., Yamani, M., Pourghasemi, H.R., Lombardo, L., 2018c. Spatial modelling of gully erosion using evidential belief function, logistic regression and a new ensemble EBF–LR algorithm. Land Degrad. Dev. 29 (11), 4035–4049. https://doi.org/10.1002/ldr.3151. Arabameri, A., Pradhan, B., Pourghasemi, H.R., Rezaei, K., 2018d. Identification of erosion-prone areas using different multi-criteria decision-making techniques and GIS. Geomat. Nat. Haz. Risk 9, 1129–1155. https://doi.org/10.1080/19475705.2018. 1513084. Arabameri, A., Pradhan, B., Rezaei, K., 2019b. Gully erosion zonation mapping using integrated geographically weighted regression with certainty factor and random forest models in GIS. J. Environ. Manag. 232, 928–942. https://doi.org/10.1016/j. jenvman.2018.11.110. Arabameri, A., Rezaei, K., Cerda, A., Lombardo, L., Rodrigo-Comino, J., 2019c. GIS-based groundwater potential mapping in Shahroud plain, Iran. A comparison among statistical (bivariate and multivariate), data mining and MCDM approaches. Sci. Total Environ. 658, 160–177. Arabameri, A., Cerda, A., Tiefenbacher, J.P., 2019d. Spatial pattern analysis and prediction of gully erosion using novel hybrid model of entropy-weight of evidence. Water 11 (6), 1129. Arabameri, A., Pradhan, B., Rezaei, K., Conoscenti, C., 2019e. Gully erosion susceptibility mapping using GIS-based multi-criteria decision analysis techniques. CATENA 180, 282–297. Arabameri, A., Pradhan, B., Rezaei, K., Lee, C.-W., 2019f. Assessment of landslide susceptibility using statistical-and artificial intelligence-based FR–RF integrated model and multiresolution DEMs. Remote Sens. 11 (9), 999. Arabameri, A., Pradhan, B., Rezaei, K., Saro, L., Sohrabi, M., 2019h. An ensemble model for landslide susceptibility mapping in a forested area. Geocarto Int. https://doi.org/ 10.1080/10106049.2019.1585484. Arabameri, A., Pradhan, B., Rezaei, K., Sohrabi, M., Kalantari, Z., 2019i. GIS-based landslide susceptibility mapping using numerical risk factor bivariate model and its ensemble with linear multivariate regression and boosted regression tree algorithms. J. Mt. Sci. 16 (3), 595–618. Arabameri, A., Rezaei, K., Cerda, A., Conoscenti, C., Kalantari, Z., 2019j. A comparison of statistical methods and multi-criteria decision making to map flood hazard susceptibility in northern Iran. Sci. Total Environ. 660, 443–458. Azareh, A., Rahmati, O., Rafiei-Sardooi, E., Sankey, J.B., Lee, S., Shahabi, H., Bin Ahmad, B., 2019. Modelling gully-erosion susceptibility in a semi-arid region, Iran: investigation of applicability of certainty factor and maximum entropy models. Sci. Total Environ. 655, 684–696. https://doi.org/10.1016/j.scitotenv.2018.11.235. Barnes, N., Luffman, I., Nandi, A., 2016. Gully erosion and freeze-thaw processes in clayrich soils, northeast Tennessee, USA. GeoRes J. 9–12, 67–76. https://doi.org/10. 1016/j.grj.2016.09.001. Beven, K.J., Kirkby, M.J., 1979. A physically based, variable contributing area model of basin hydrology/un modèle à base physique de zone d'appel variable de l'hydrologie du bassin versant. J. Hydrol. Sci. Bull. (1), 43–69. https://doi.org/10.1080/ 02626667909491834. Borrelli, P., Robinson, D.A., Fleischer, L.R., Lugato, E., Ballabio, C., Alewell, C., Meusburger, K., Modugno, S., Schütt, B., Ferro, V., Bagarello, V., Van Oost, L., Montanarella, L., Panagos, P., 2017. An assessment of the global impact of 21st century land use change on soil erosion. Nat. Commun. 8, 1–15. https://doi.org/10. 1038/s41467-017-02142-7. Bout, B., Lombardo, L., van Westen, C.J., Jetten, V.G., 2018. Integration of two-phase solid fluid equations in a catchment model for flashfloods, debris flows and shallow slope failures. Environ. Model. Softw. 105, 1–16. Breiman, L., 1998. Arcing classifiers. Ann. Stat. 26, 801–849. https://doi.org/10.1214/ aos/1024691079. Castro Camilo, D., Lombardo, L., Mai, P.M., Dou, J., Huser, R., 2017. Handling high predictor dimensionality in slope-unit-based landslide susceptibility models through LASSO-penalized generalized linear model. Environ. Model. Softw. 97, 145–156. Chaplot, V., Coadou le Brozec, E., Silvera, N., Valentin, C., 2005. Spatial and temporal assessment of linear erosion in catchments under sloping lands of northern Laos. Catena 63, 167–184. https://doi.org/10.1016/j.catena.2005.06.003. Chen, W., Pradhan, B., Li, S., Shahabi, H., Rizeei, H.M., Hou, E., Wang, S., 2019. Novelhybrid integration approach of bagging-based fisher’s linear discriminant function forgroundwater potential analysis. Nat. Resour. Res. 1–20. https://doi.org/ 10.1007/s11053-019-09465-w. Conforti, M., Aucelli, P.P., Robustelli, G., Scarciglia, F., 2011. Geomorphology and GIS analysis for mapping gully erosion susceptibility in the Turbolo stream catchment (Northern Calabria, Italy). Nat. Hazards 56 (3), 881–898. https://doi.org/10.1007/ s11069-010-9598-2. Conoscenti, C., Angileri, S., Cappadonia, C., Rotigliano, E., Agnesi, V., Märker, M., 2014. Gully erosion susceptibility assessment by means of GIS-based logistic regression: a case of Sicily (Italy). Geomorphology 204 (1), 399–411. https://doi.org/10.1016/j. geomorph.2013.08.021.

Funding This research is supported by the Centre for Advanced Modelling and Geospatial Information Systems (CAMGIS), UTS under grant numbers 323930, and 321740.2232335 and 321740.2232357, 321740.2232424; and 321740.2232452. References Abedini, M., Tulabi, S., 2018. Assessing LNRF, FR, and AHP models in landslide susceptibility mapping index: a comparative study of Nojian watershed in Lorestan province, Iran. Environ. Earth Sci. 77, 405. https://doi.org/10.1007/s12665-0187524-1. Aghdam, I.N., Varzandeh, M.H.M., Pradhan, B., 2016. Landslide susceptibility mappingusing an ensemble statistical index (Wi) and adaptive neuro-fuzzy inference system (ANFIS) model at Alborz Mountains (Iran). Environ. Earth Sci. 75 (7), 553. https:// doi.org/10.1007/s12665-015-5233-6. Akgun, A., Türk, N., 2011. Mapping erosion susceptibility by a multivariate statistical method: a case study from the Ayvalik region, NW Turkey. Comput. Geosci. 37 (9), 1515–1524. https://doi.org/10.1016/j.cageo.2010.09.006. Antoneli, V., Rebinski, E.A., Bednarz, J.A., Rodrigo-Comino, J., Keesstra, S.D., Cerdà, A., Pulido Fernández, M., 2018. Soil Erosion induced by the introduction of new pasture species in a faxinal farm of Southern Brazil. Geosciences 8 (5), 166. Arabameri, A., Pourghasemi, H.R., 2019. Spatial modeling of gully erosion using linear and quadratic discriminant analyses in GIS and R. In: Pourghasemi, H.R., Gokceoglu, C. (Eds.), Spatial Modeling in GIS and R for Earth and Environmental Sciences, First edition. Elsevier publication. https://doi.org/10.1016/B978-0-12-815226-3.000132. 796 p. Arabameri, A., Pourghasemi, H.R., Cerda, A., 2017a. Erodibility prioritization of subwatersheds using morphometric parameters analysis and its mapping: a comparison among TOPSIS, VIKOR, SAW, and CF multi-criteria decision making models. Sci.

13

Catena 183 (2019) 104223

A. Arabameri, et al.

s40808-016-0170-2. Lee, S., Pradhan, B., 2007. Landslide hazard mapping at Selangor, Malaysia using frequency ratio and logistic regression model. Landslides 4, 33–41. https://doi.org/10. 1007/s10346-006-0047-y. Lombardo, L., Mai, P.M., 2018. Presenting logistic regression-based landslide susceptibility results. Eng. Geol. 244, 14–24. https://doi.org/10.1016/j.enggeo.2018.07. 019. Lombardo, L., Cama, M., Conoscenti, C., Märker, M., Rotigliano, E., 2015. Binary logistic regression versus stochastic gradient boosted decision trees in assessing landslide susceptibility for multiple-occurring landslide events: application to the 2009 storm event in Messina (Sicily, southern Italy). Nat. Hazards 79 (3), 1621–1648. https:// doi.org/10.1007/s11069-015-1915-3. Lombardo, L., Opitz, T., Huser, R., 2018. Point process-based modeling of multiple debris flow landslides using INLA: an application to the 2009 Messina disaster. Stoch. Env. Res. Risk A. 32 (7), 2179–2198. Lombardo, L., Opitz, T., Huser, R., 2019. Numerical recipes for landslide spatial prediction using R-INLA: a step-by-step tutorial. In: Spatial Modeling in GIS and R for Earth and Environmental Sciences. Elsevier, pp. 55–83. Luffman, I.E., Nandi, A., Spiegel, T., 2015. Gully morphology, hillslope erosion, and precipitation characteristics in the Appalachian Valley and Ridge province, southeastern USA. Catena 133, 221–232. https://doi.org/10.1016/j.catena.2015.05.015. Magliulo, P., 2012. Assessing the susceptibility to water-induced soil erosion using a geomorphological, bivariate statistics-based approach. Environ. Earth Sci. 67, 1801–1820. https://doi.org/10.1007/s12665-012-1634-y. Märker, M., Pelacani, S., Schröder, B., 2011. A functional entity approach to predict soil erosion processes in a small Plio-Pleistocene Mediterranean catchment in Northern Chianti, Italy. Geomorphology 125, 530–540. https://doi.org/10.1016/j.geomorph. 2010.10.022. McCloskey, G.L., Wasson, R.J., Boggs, G.S., Douglas, M., 2016. Timing and causes of gully erosion in the riparian zone of the semi-arid tropical Victoria River, Australia: management implications. Geomorphology 266, 96–104. https://doi.org/10.1016/j. geomorph.2016.05.009. Meliho, M., Khattabi, A., Mhammdi, N., 2018. A GIS-based approach for gully erosion susceptibility modelling using bivariate statistics methods in the Ourika watershed, Morocco. Environ. Earth Sci. 77, 655. https://doi.org/10.1007/s12665-018-7844-1. Meten, M., Bhandary, N.P., Yatabe, R., 2016. GIS-based frequency ratio and logistic regression modelling for landslide susceptibility mapping of Debre Sina area in central Ethiopia. J. Mt. Sci. 12 (6), 1355–1372. https://doi.org/10.1007/s11629-0153464-3. Mohammadi, A.T., Irani, A., Sorur, J., 2014. The preparation of landslide map by landslide numerical risk factor (LNRF) model and geographic information system (GIS). Egypt. J. Remote Sens. Space Sci. 17 (2), 159–170. https://doi.org/10.1016/j.ejrs. 2014.08.001. Mohammady, M., Pourghasemi, H.R., Pradhan, B., 2012. Landslide susceptibility mapping at Golestan Province, Iran: A comparison between frequency ratio, Dempster–Shafer, and weights-of-evidence models. J. Asian Earth Sci. 61 (15), 221–236. https://doi.org/10.1016/j.jseaes.2012.10.005. Mojaddadi, H., Pradhan, B., Nampak, H., Ahmad, N., Ghazali, A.H., 2017. Ensemble machine-learning-based geospatial approach for flood risk assessment using multisensor remote-sensing data and GIS. Geomatics, Nat. Hazards and Risk 8 (2), 1080–1102. https://doi.org/10.1080/19475705.2017.1294113. Moore, I.D., Grayson, R.B., Ladson, A.R., 1991. Digital terrain modeling: a review of hydrological, geomorphological, and biological applications. Hydrol. Process. 5, 3–30. https://doi.org/10.1002/hyp.3360050103. Naghibi, S.A., Pourghasemi, H.R., Dixon, B., 2016. GIS-based groundwater potential mapping using boosted regression tree, classification and regression tree, and random forest machine learning models in Iran. Environ. Monit. Assess. 188, 44. https://doi. org/10.1007/s10661-015-5049-6. Nazari Samani, A., Ahmadi, H., Jafari, M., Boggs, G., Ghoddousi, J., Malekian, A., 2009. Geomorphic threshold conditions for gully erosion in southwestern Iran (BoushehrSamal watershed). J. Asian Earth Sci. 35, 180–189. https://doi.org/10. 1016/j.jseaes.2009.02.004. Ollobarren, P., Capra, A., Gelsomino, A., La Spada, C., 2016. Effects of ephemeral gully erosion on soil degradation in a cultivated area in Sicily (Italy). Catena 145, 334–345. https://doi.org/10.1016/j.catena.2016.06.031. Ozdemir, A., 2011. GIS-based groundwater spring potential mapping in the Sultan Mountains (Konya, Turkey) using frequency ratio, weights of evidence and logistic regression methods and their comparison. J. Hydrol. 411 (3–4), 290–308. https:// doi.org/10.1016/j.jhydrol.2011.10.010. Poesen, J., Nachtergaele, J., Verstraeten, G., Valentin, C., 2003. Gully erosion and environmental change: importance and research needs. Catena 50 (2–4), 91–133. https://doi.org/10.1016/S0341-8162(02)00143-1. Pourghasemi, H.R., Yousefi, S., Kornejady, A., Cerdà, A., 2017. Performance assessment of individual and ensemble data-mining techniques for gully erosion modelling. Sci. Total Environ. 609, 764–775. https://doi.org/10.1016/j.scitotenv.2017.07.198. Pradhan, B., Lee, S., 2010. Regional landslide susceptibility analysis using backpropagation neural network model at Cameron Highland, Malaysia. Landslides 7 (1), 13–30. https://doi.org/10.1007/s10346-009-0183-2. Rahmati, O., Haghizadeh, A., Pourghasemi, H.R., Noormohamadi, F., 2016. Gully erosion susceptibility mapping: the role of GIS based bivariate statistical models and their comparison. Nat. Hazards 82, 1231–1258. https://doi.org/10.1007/s11069-0162239-7. Rahmati, O., Tahmasebipour, N., Haghizadeh, A., Pourghasemi, H.R., Feizizadeh, B., 2017a. Evaluation of different machine learning models for predicting and mapping the susceptibility of gully erosion. Geomorphology 298, 118–137. https://doi.org/10. 1016/j.geomorph.2017.09.006.

Conoscenti, C., Agnesi, V., Cama, M., Alamaru Caraballo-Arias, N., Rotigliano, E., 2018. Assessment of gully erosion susceptibility using multivariate adaptive regression splines and accounting for terrain connectivity. Land Degrad. Dev. 29 (3), 724–736. https://doi.org/10.1002/ldr.2772. Croke, J.C., Hairsine, P.B., 2006. Sediment delivery in managed forests: a review. Environ. Rev. 14, 59–87. https://doi.org/10.1139/a05-016. Croke, J., Mockler, S., 2001. Gully initiation and road-to-stream linkage in a forested catchment, southeastern Australia. Earth Surf. Process. Landf. 26, 205–217. https:// doi.org/10.1002/1096-9837. Dehnavi, A., Aghdam, I.N., Pradhan, B., Varzandeh, M.H.M., 2015. A new hybrid model using step-wise weight assessment ratio analysis (SWARA) technique and adaptive neuro-fuzzy inference system (ANFIS) for regional landslide hazard assessment in Iran. Catena 135, 122–148. https://doi.org/10.1016/j.catena.2015.07.020. Dewitte, O., Daoudi, M., Bosco, C., Eeckhaut, M., 2015. Predicting the susceptibility to gully initiation in data-poor regions. Geomorphology 228, 101–115. https://doi.org/ 10.1016/j.geomorph.2014.08.010. Ding, L., Chen, K.L., Cheng, S.G., Wang, X., 2015. Water ecological carrying capacity of urban lakes in the context of rapid urbanization: a case study of East Lake in Wuhan. Phys. Chem. Earth 89–90, 104–113. https://doi.org/10.1016/j.pce.2015.08.004. Dube, F., Nhapi, I., Murwira, A., Gumindoga, W., Goldin, J., Mashauri, D.A., 2014. Potential of weight of evidence modelling for gully erosion hazard assessment in Mbire District – Zimbabwe. Phys. Chem. Earth 67, 145–152. https://doi.org/10. 1016/j.pce.2014.02.002. Dymond, J.R., Herzig, A., Basher, L., Betts, H.D., Marden, M., Phillips, C.J., Ausseil, A.G.E., Palmer, D.J., Clark, M., Roygard, J., 2016. Development of a New Zealand SedNet model for assessment of catchment-wide soil-conservation works. Geomorphology 257, 85–93. https://doi.org/10.1016/j.geomorph.2015.12.022. Fanos, A.M., Pradhan, B., 2019. Spatial Ensemble Model for Rockfall Source Identification From High Resolution LiDAR Data and GIS. IEEE Access. 7, 74570–74585. https:// doi.org/10.1109/ACCESS.2019.2919977. Gómez-Gutiérrez, Á., Conoscenti, C., Angileri, S.E., Rotigliano, E., Schnabel, S., 2015. Using topographical attributes to evaluate gully erosion proneness (susceptibility) in two Mediterranean basins: advantages and limitations. Nat. Hazards 79, 291–314. https://doi.org/10.1007/s11069-015-1703-0. Gong, J.G., Jia, Y.W., Zhou, Z.H., Wang, Y., Wang, W.L., Peng, H., 2011. An experimental study on dynamic processes of ephemeral gully erosion in loess landscapes. Geomorphology 125, 203–213. https://doi.org/10.1016/j.geomorph.2010.09.016. Goodwin, N.R., Armston, J.D., Muir, J., Stiller, I., 2017. Monitoring gully change: a comparison of airborne and terrestrial laser scanning using a case study from Aratula, Queensland. Geomorphology 282, 195–208. https://doi.org/10.1016/j.geomorph. 2017.01.001. Gupta, R.P., Joshi, B.C., 1990. Landslide hazard zonation using the GIS approach—a case study from the Ramganga catchment, Himalayas. Eng. Geol. 28, 119–131. https://doi. org/10.1016/0013-7952(90)90037-2. Hălbac-Cotoară-Zamfir, R., Keesstra, S., Kalantari, Z., 2019. The impact of political, socioeconomic and cultural factors on implementing environment friendly techniques for sustainable land management and climate change mitigation in Romania. Sci. Total Environ. 654, 418–429. https://doi.org/10.1016/j.scitotenv.2018.11.160. Heedegen, R.G., Beran, M., 1982. Quantifying source areas through land surface curvature and shape. J. Hydrol. 57 (s 3–4), 359–373. https://doi.org/10.1016/00221694(82)90155-X. Hong, H., Liu, J., Zhu, A.X., Shahabi, H., Pham, B.T., Chen, W., Pradhan, B., Tien Bui, D., 2017. A novel hybrid integration model using support vector machines and random subspace for weather-triggered landslide susceptibility assessment in the Wuning area (China). Environ. Earth Sci. 76 (652). https://doi.org/10.1007/s12665-0176981-2. Horton, R.E., 1945. Erosional development of streams and their drainage basins: hydrophysical approach to quantitative morphology. Bull. Geol. Soc. Am. 56, 275–370. https://doi.org/10.1130/0016-7606. I.R. of Iran Meteorological Organization (IRIMO), 2012. Available online: http:// www.mazandaranmet.ir. Geology survey of Iran (GSI). 1997. http://www.gsi.ir/ Main/Lang_en/index.html. Issaka, S., Aqeel Ashraf, M., 2017. Impact of soil erosion and degradation on water quality: a review. Geol. Ecol. Landscapes 1, 1–12. https://doi.org/10.1080/ 24749508.2017.1301053. Keesstra, S.D., Bouma, J., Wallinga, J., Tittonell, P., Smith, P., Cerdà, A., Bardgett, R.D., 2016. The significance of soils and soil science towards realization of the United Nations Sustainable Development Goals. Soil 2 (2), 111–128. https://doi.org/10. 5194/soil-2-111-2016. Keesstra, S., Mol, G., de Leeuw, J., Okx, J., de Cleen, M., Visser, S., 2018. Soil-related sustainable development goals: four concepts to make land degradation neutrality and restoration work. Land 7 (4), 133. https://doi.org/10.3390/land7040133. Keesstra, S.D., Rodrigo-Comino, J., Novara, A., Giménez-Morera, A., Pulido, M., Di Prima, S., Cerdà, A., 2019. Straw mulch as a sustainable solution to decrease runoff and erosion in glyphosate-treated clementine plantations in Eastern Spain. An assessment using rainfall simulation experiments. Catena 174, 95–103. https://doi.org/10.1016/ j.catena.2018.11.007. Köthe, R., Lehmeier, F., 1993. SAGA – Ein programmsystem zur automatischen reliefanalyse. Z. Angew. Geogr. 1993 (4), 11–21. Kropacek, J., Schillaci, C., Salvini, R., Marker, M., 2016. Assessment of gully erosion in the Upper Awash, Central Ethiopian highlands based on a comparison of archived aerial photographs and very high resolution satellite images. Geogr. Fis. Din. Quat. 39 (2), 161–170. Kumar Samanta, R., Sankar Bhunia, G., Kumar shit, P., 2016. Spatial modelling of soil erosion susceptibility mapping in lower basin of Subarnarekha river (India) based on geospatial techniques. Model. Earth Sys. Environ. 2, 99. https://doi.org/10.1007/

14

Catena 183 (2019) 104223

A. Arabameri, et al.

Vaezi, A.R., Abbasi, M., Bussi, G., Keesstra, S., 2017. Modeling sediment yield in semi-arid pasture micro-catchments, NW Iran. Land Degrad. Dev. 28 (4), 1274–1286. https:// doi.org/10.1002/ldr.2526. Van Westen, C.J., Castellanos, E., Kuriakose, S.L., 2008. Spatial data for landslide susceptibility, hazard, and vulnerability assessment: an overview. Eng. Geol. 102, 112–131. https://doi.org/10.1016/j.enggeo.2008.03.010. Vanacker, V., Govers, G., Poesen, J., Deckers, J., Dercon, G., Loaiza, G., 2003. The impact of environmental change on the intensity and spatial pattern of water erosion in a semi-arid mountainous Andean environment. Catena 51, 329–347. https://doi.org/ 10.1016/S0341-8162(02)00172-8. Vandekerckhove, L., Poesen, J., Oostwoud Wijdenes, D., Nachtergaele, J., Kosmas, C., Roxo, M.J., Figueiredo, T.D.E., 2000. Thresholds for gully initiation and sedimentation in Mediterranean Europe. Earth Surf. Process. Landf. 25, 1201–1220. https:// doi.org/10.1002/1096-9837. Vanmaercke, M., Poesen, J., Van Mele, B., Demuzere, M., Bruynseels, A., Golosov, V., Bezerra, J.F.R., Bolysov, S., Dvinskih, A., Frankl, A., 2016. How fast do gully headcuts retreat? Earth Sci. Rev. 154, 336–355. https://doi.org/10.1016/j.earscirev.2016.01. 009. Wilkinson, M.T., Humphreys, G.S., 2006. Slope aspect, slope length and slope inclination controls of shallow soils vegetated by sclerophyllous heath—links to long-term landscape evolution. Geomorphology 76 (3–4), 347–362. https://doi.org/10.1016/j. geomorph.2005.11.011. Yesilnacar, E.K., 2005. The Application of Computational Intelligence to Landslide Susceptibility Mapping in Turkey. Ph.D Thesis. Department of Geomatics the University of Melbourne (423p). Yilmaz, C., Topal, T., Suzen, M.L., 2012. GIS-based landslide susceptibility mapping using bivariate statistical analysis in Devrek (Zonguldak-Turkey). Environ. Earth Sci. 65, 2161–2178. https://doi.org/10.1007/s12665-011-1196-4. Youssef, A.M., Al-Kathery, M., Pradhan, B., 2015. Landslide susceptibility mapping at AlHasher area, Jizan (Saudi Arabia) using GIS-based frequency ratio and index of entropy models. Geosci. J. 19 (1), 113–134. https://doi.org/10.1007/s12303-0140032-8. Zabihi, M., Mirchooli, F., Motevalli, A., Darvishan, A.K., Pourghasemi, H.R., Zakeri, M.A., Sadighi, F., 2018. Spatial modelling of gully erosion in Mazandaran Province, northern Iran. Catena 161, 1–13. https://doi.org/10.1016/j.catena.2017.10.010. Zare, M., Pourghasemi, H.R., Vafakhah, M., Pradhan, B., 2013. Landslide susceptibilitymapping at VazWatershed (Iran) using an artificial neural networkmodel: a comparison between multilayer perceptron (MLP) and radial basic function (RBF) algorithms. Arab. J. Geosci. 6 (8), 2873–2888. https://doi.org/10.1007/s12517-0120610-x. Zevenbergen, L.W., Thorne, C.R., 1987. Quantitative analysis of land surface topography. Earth Surf. Process. Landf. 12 (1), 47–56. https://doi.org/10.1002/esp.3290120107.

Rahmati, O., Tahmasebipour, N., Haghizadeh, A., Pourghasemi, H.R., Feizizadeh, B., 2017b. Evaluating the influence of geo-environmental factors on gully erosion in a semi-arid region of Iran: an integrated framework. Sci. Total Environ. 579, 913–927. https://doi.org/10.1016/j.scitotenv.2016.10.176. Rahmati, O., Kornejady, A., Samadi, M., Deo, R.C., Conoscenti, C., Lombardo, L., Dayal, K., Taghizadeh-Mehrjardi, R., Pourghasemi, H.R., Kumar, S., Bui, D.T., 2019. PMT: new analytical framework for automated evaluation of geo-environmental modelling approaches. Sci. Total Environ. 664, 296–311. Rasyid, A.R., Bhandary, N.P., Yatabe, R., 2016. Performance of frequency ratio and logistic regression model in creating GIS based landslides susceptibility map at Lompobattang Mountain, Indonesia. Geoenviron. Dis. 3, 19. https://doi.org/10. 1186/s40677-016-0053-x. Rodrigo-Comino, J., Keesstra, S., Cerdà, A., 2018. Soil erosion as an environmental concern in vineyards. The case study of Celler del Roure, Eastern Spain, by means of rainfall simulation experiments. Beverages 4 (2), 31. https://doi.org/10.3390/ beverages4020031. Romshoo, S.A., Bhat, S.A., Rashid, I., 2012. Geo informatics for assessing the morphometric control on hydrological response at watershed scale in the Upper Indus Basin. J. Earth Syst. Sci. 121 (3), 659–686. https://doi.org/10.1007/s12040-012-0192-8. Rouse, J.W., Haas, R.H., Schell, J.A., Deering, D.W., 1974. Monitoring vegetation systems in the Great Plains with ERTS. In: Freden, S.C., Mercanti, E.P., Becker, M. (Eds.), Third Earth Resources Technology Satellite–1 Syposium. Volume I: Technical Presentations, NASA SP-351. NASA, Washington, D.C., pp. 309–317. Shellberg, J.G., Spencer, J., Brooks, A.P., Pietsch, T.J., 2016. Degradation of the Mitchell River fluvialmegafan by alluvial gully erosion increased by post-European land use change, Queensland, Australia. Geomorphology 266, 105–120. https://doi.org/10. 1016/j.geomorph.2016.04.021. Shit, P.K., Paira, R., Bhunia, G., Maiti, R., 2015. Modeling of potential gully erosion hazard using geo-spatial technology at Garbheta block, West Bengal in India. Model. Earth Sys. Environ. 1 (1–2), 1–16. https://doi.org/10.1007/s40808-015-0001-x. Solomun, M.K., Barger, N., Cerda, A., Keesstra, S., Marković, M., 2018. Assessing land condition as a first step to achieving land degradation neutrality: a case study of the Republic of Srpska. Environ. Sci. Pol. 90, 19–27. https://doi.org/10.1016/j.envsci. 2018.09.014. Süzen, M.L., Doyuran, V., 2004. A comparison of the GIS based landslide susceptibility assessment methods: multivariate versus bivariate. Environ. Geol. 45 (5), 665–679. https://doi.org/10.1007/s00254-003-0917-8. Therneau, T.M., Atkinson, B., Ripley, B., 2014. RPART: Recursive Partitioning and Regression Trees. R Package Version 4. pp. 1–8. Tien Bui, D., Ho, T.C., Pradhan, B., Pham, B.-T., Nhu, V.-H., Revhaug, I., 2016. GIS-based modeling of rainfall-induced landslides using data mining-based functional trees classifier with AdaBoost, Bagging, and MultiBoost ensemble frameworks. Environ. Earth Sci. 75, 1101. https://doi.org/10.1007/s12665-016-5919-4.

15