Gully Erosion Susceptibility Mapping using Machine Learning Techniques in Ihioma Community, Imo State, Nigeria

639 Nigerian Research Journal of Engineering and Environmental Sciences 10(2) 2025 pp. 639-654 p ISSN: 2635-3342; e ISSN: 2635-3350 Original Research Article Gully Erosion Susceptibility Mapping using Machine Learning Techniques in Ihioma Community, Imo State, Nigeria *Abiodun, O.E., Ohikhueme A.I. and Alabi, A.O. Department of Surveying and Geoinformatics, University of Lagos, Akoka, Yaba, Lagos, Nigeria. *[email protected]; [email protected] http://doi.org/10.5281/zenodo.18062130 ARTICLE INFORMATION ABSTRACT Article history: Received 24 Nov. 2025 Revised 08 Dec. 2025 Accepted 09 Dec. 2025 Available online 30 Dec. 2025 Gully erosion poses significant environmental and socioeconomic challenges in many regions worldwide, including the Ihioma community in Imo State, Nigeria. This study employs Machine Learning techniques to predict and map gully erosion susceptibility in the Ihioma community, contributing to the understanding of erosion processes and informing targeted mitigation strategies. Environmental factors such as rainfall, elevation, slope, topographic indices, vegetation cover, and land use/land cover are analyzed to identify primary drivers of erosion vulnerability. The Extreme Gradient Boosting and Random Forest algorithms were compared for their effectiveness in predicting erosion susceptibility, with Random Forest demonstrating superior performance with accuracy of 79.8% and recall of 90.4%. Feature Importance analysis ranked elevation, Normalized Difference Vegetation Index (NDVI), and rainfall in that order as critical variables influencing erosion susceptibility. Using each algorithm, susceptibility maps were created that indicated locations that are likely to experience gully erosion. From the results presented, about 35% of the study area have high susceptibility to gully erosion. The susceptibility map produced could serve as basis for proper planning in a way to reduce the negative impact of gully erosion in the study area. © 2025 RJEES. All rights reserved. Keywords: Erosion Susceptibility Extreme gradient boosting Random forest Mapping 1. INTRODUCTION Soil stands as a valuable and irreplaceable natural asset, providing essential ecosystem services and fulfilling numerous ecological roles crucial for sustaining life on our planet. It contributes to the production of food, fodder, and timber, aids in water purification and storage, facilitates nutrient cycling, filters out harmful substances, and supports biodiversity preservation, among other functions (Lorenz, 2019). These roles and significance notwithstanding, almost one-third of the world's soils face degradation, with an annual loss of 25–40 billion tons due to erosion (FAO, 2019; UNCCD, 2019). This degradation has profound implications for the productivity, resilience, and sustainability of agricultural 640 O.E. Abiodun et al. / Nigerian Research Journal of Engineering and Environmental Sciences 10(2) 2025 pp. 639-654 and ecological systems (Omid Rahmati et al., 2017). While soil erosion is a global issue, Africa bears the brunt of its impact as a result of unplanned mitigation strategies. Gully erosion happens when water flows across exposed ground, eroding the soil near the drainage lines. In an ideal situation where the land surface has been left undisturbed, vegetation would reduce runoff and hold the soil together, shielding it from excessive runoff and direct downpour. However, it is a common occurrence now to leave the soil exposed as a result of various human activities thereby making the coil incapable of absorbing extra rainfall. The development of gully erosion in sensitive areas is thus made possible by the subsequent increase in surface runoff that concentrates in drainage lines. Gully erosion is one of the major kinds of erosion that brings about land degradation problems (Chandrakala, 2020). Gully erosion causes both a decrease in the amount and quality of arable soil in various regions of the world. Gully erosion also causes several environmental problems such as desertification, inundation, and sedimentation in lakes (Torri et al., 2012). It can also lead to the reduction of soil fertility and agricultural productivity, thereby impacting negatively on both global and national economy (Kirkby & Bracken, 2009). The production of gully erosion maps involves the use of the association between gully erosion occurrence and their geo-environmental conditioning factors (Lana, 2023). Many studies have used the Universal Soil Loss Equation (USLE) and the Revised Universal Soil Loss Equation (RUSLE) to compute the intensity of soil erosion in various regions of the world (Panditharathne, 2019; Thakuriah, 2023). However, these models may not necessarily capture the complexities of soil erosion because the model simplifies the process of erosion without accounting for all factors that unarguably contribute to soil loss (Borrelli et al., 2021; Ezeh et al., 2024). Hence, this makes it difficult to accurately quantify soil loss and subsequently determine locations susceptible to soil erosion. Despite the fact that Orlu is located in a highly erosion susceptible area (South East Nigeria), studies on soil erosion in Nigeria which are mostly USLE and RUSLE based have been carried out elsewhere in Nigeria [Katsina (Adediji et al., 2010); Uyo (Adeola et al., 2016); Edda-Afikpo (Amah et al., 2020); Osun (Fasinmirin & Olorunfemi, 2013) among others] but none in Orlu. It is therefore important to carry out studies on erosion susceptibility in Orlu using an approach that will present a more accurate result. Although there are recent developments in remote sensing, machine learning, and geographic information systems (GIS), there is a gap in their integration for gully erosion susceptibility mapping especially in Nigeria (Igbokwe, 2008). Leveraging these advanced technologies can enhance the accuracy, efficiency, and scalability of susceptibility mapping models. This study therefore focuses on the use of selected Machine Learning algorithms (Random Forest and Extreme Gradient Boosting (XGBoost)) in predicting and mapping locations that are vulnerable to gully erosion within Ihioma community, located at Orlu Local Government Area, Imo State. 2. MATERIALS AND METHODS 2.1. Study Area Ihioma is a community in Orlu Local Government Area, Imo state in Southeastern Nigeria and it is in the humid tropical rainforest belt of Nigeria. In the rainy season, the average daily temperature is 20°C, while in the dry season, it is usually 33°C. The average yearly temperature ranges from 26.5 to 27.3°C, while the relative humidity is between 65 and 75%. The area's vegetation is made up of shrubs and trees from Nigeria's rainforest belt. However, because of human activities like farming and building infrastructure, the majority of the vegetation has been cleared. The study area is shown in Figure 1. 641 O.E. Abiodun et al. / Nigerian Research Journal of Engineering and Environmental Sciences 10(2) 2025 pp. 639-654 Figure 1: Map of the study area 2.2. Material Collection and Preparation of Samples Based on an extensive literature review and expert review, the factors of soil erosion used in this study and their sources are as contained in Table 1. The data processing phase started with the determination of the erosion triggering factors. 10 geographic factors, 7 of which are continuous and 3 of which are categorical (Table 1), are taken into consideration in this study as potential gully erosion influences. And most of these factors were extracted using geospatial tools and techniques in ArcMap. The machine learning model will analyse the data to determine how each feature impacts suitability. Table 1: Factors and types of erosion triggering factors Data type Source Digitized map of the study area Author Rainfall NASA (gpm.nasa.gov) Elevation www.cgiar-csi.org/data/srtm30m-digital-elevationdatabase v4-1and https://lta.cr.usgs.gov/SRTM Normalized Difference Vegetative Index (NDVI) NASA (gpm.nasa.gov) Land Cover Esri Sentinel-2 Land, Cover Explorer Soil Data FAO Geological Data Imo State Ministry of Industry and Solid Minerals Slope Shuttle Radar Topography Mission (SRTM) Stream power index (SPI) Topographical Wetness index (TWI) Shuttle Radar Topography Mission (SRTM) Topographic Position Index (TPI) Shuttle Radar Topography Mission (SRTM) All variables were further preprocessed in ArcGIS (ArcMap) as a type of data preparation for the machine learning model development. The study area boundary was overlaid on all the conditioning factors for the clipping of the features required for the model. 1. Rainfall: According to Roy and Saha (2019), rainfall is a key factor in determining the likelihood of gully development in a certain area, indicating the climate conditions that are common in the research 642 O.E. Abiodun et al. / Nigerian Research Journal of Engineering and Environmental Sciences 10(2) 2025 pp. 639-654 area. After acquiring the study area's annual mean rainfall, "Multidimensional Tools" and "Make NetCDF Raster Layer" in ArcGIS 10.8 were used to convert the rainfall data to a raster layer. Using "Conversion Tools," "From Raster," and "Raster to point," the raster layer got converted to points. Afterwards, the acquired points were inputted into "Spatial Analyst Tools" from where Kriging interpolation was performed to get the map of rainfall. 2. DEM: elevation has a major influence in determining the vulnerability to gully erosion at different time, as observed in the occurrence and progression of gully erosion (Zabihi et al., 2018). This is attributed to its influence on factors such as vegetation, precipitation patterns, and the dynamics of gully erosion itself (Golestani et al., 2014). Satellite Imageries from SRTM were downloaded. These images were imported to ArcGIS 10.8 interface and clipped. Next, the data was reclassified and a DEM map was created. 3. Slope: Surface runoff on slopes and surface drainage patterns are important causes of erosion. According to Conforti et al. (2011) and Luca et al. (2011), they are widely recognized as crucial indicators influencing the processes of gully erosion. ArcMap 10.8 software was used to build the slope. Select the Spatial Analyst Tool >> Surface >> Slope option from the ArcTool Box. The Digital Elevation Model raster in a projected coordinate system was inputted and used to calculate the slope, reclassified into classes, and created a Slope map. 4. Land cover: After the LANDSAT imagery for the study area was downloaded. The Land Use Land Cover (LULC) was classified through supervised classification in ArcGIS 10.8. In order to create signature files, the supervised classification process involves choosing and digitizing polygons, then putting them in a "Area of Interest" layer. To make this LULC specific, numerous polygons for a particular LULC type were created. Although the supervised classification method takes a lot of time, Enderle and Weih (2005) found that it produced results with higher general accuracy than unsupervised classification. 5. Normalized differenced vegetation index: The Normalized Difference Vegetation Index (NDVI) serves as a reliable indicator of photosynthetic activity, as highlighted in prior research (Pourghasemi et al., 2014). The NDVI was estimated from the Landsat 8 imagery downloaded and was calculated in the raster calculation in ArcGIS (Arc toolbox) using Equation 1. NDVI = 𝑁𝐼𝑅−𝑅𝑒𝑑 𝑁𝐼𝑅+𝑅𝑒𝑑 (1) where NIR and Red values indicate the infrared and red sections of the electromagnetic spectrum, respectively. 6. Topographical Wetness Index (TWI): As one of the factors leading to gully erosion, the topographic wetness index measures the amount of water in the studied area (Moore & Wilson, 1992). According to Moore et al. (1991), Equation 2 can be used to describe TWI: TWI = ln 𝑎 tan 𝑏 (2) Where 𝑎 is the upslope area and tan 𝑏 is the slope gradient in radians. The following steps was carried out in the ArcGIS interface: DEM >> Fill DEM >> Flow Direction >> Flow Accumulation >> Slope in degree >> Radians of slope = (slope in degree * 1.570796) \ 90 >> Tan slope = con (slope>0, tan(slope), 0.001) >>flow Accumulation scaled = (flow Accumulation+1) * cell size >> TWI = Ln (Flow Accumulation scaled/Tan slope). 7. Topographic position index (TPI): The topographic position index is a commonly used method for automating zone ordination and assessing topographic slope placement. This function generates a singleband raster that quantifies various characteristics based on elevation measurements. The following steps was carried out in the ArcGIS interface: DEM >> Focal statistics >> Minus >> TPI. 8. Stream power index (SPI): The idea that flow changes based on the particular watershed is the basis for measuring the erosive power of water flow. It can be described by Equation (3): 643 O.E. Abiodun et al. / Nigerian Research Journal of Engineering and Environmental Sciences 10(2) 2025 pp. 639-654 SPI = ln (𝑎 * tan𝑏 * 100) (3) The following steps was carried out in the ArcGIS interface: DEM >> Fill DEM >> Flow Direction >> Flow Accumulation >> Slope in degree >> SPI 9. Soil and geology: From the processing carried out, it was found that all locations within the study area had the same soil type and Geology which might bring negative value, and this would not necessarily have any importance in the model; hence, this was discarded. 2.3. Machine Learning Models Within the research area, two models were assessed for their capacity to forecast gully erosion. Moreover, the area's gully erosion map was created using these models: 1. Random forest (RF) 2. Extreme gradient boosting (XGBoost) Initially, the model was trained with the algorithms' default parameters, and the test set was used to evaluate its performance. 2.3.1. Random forest algorithm Combining several decision trees, Random Forest is an approach to ensemble learning that enhances accuracy in prediction. According to Cutler et al. (2007), Random Forest produces reliable estimates of erosion susceptibility by combining the predictions of several trees, which lowers the possibility of overfitting. Within a Random Forest, every decision tree acquires knowledge autonomously from a randomized subset of both features and training data. Mathematically, the algorithm's complexity is primarily determined by the number of decision trees (n), the number of features selected at each node (m), and the size of the training dataset. High-dimensional feature spaces in large datasets can be handled using Random Forest, which is generally computationally efficient. This diversity among individual trees helps capture different aspects of the complex relationships between environmental factors and erosion susceptibility (Liaw & Wiener, 2002). A measure of feature relevance is provided by Random Forest, which shows how each predictor variable affects the prediction performance of the model. This feature importance analysis can help identify the most influential environmental factors driving gully erosion susceptibility in the study area (Cutler et al., 2007). Because Random Forest is resistant to noisy and overfitting data, it can handle big and complicated datasets that are frequently encountered in mapping erosion susceptibility (Cutler et al., 2007). 2.3.2. XGBoost classifier The XGBoost is a gradient boosting method well-known for its effectiveness and efficiency in tasks involving predictive modelling. In order to produce a strong learner with higher accuracy, it iteratively constructs an ensemble of weak learners (decision trees) and merges their predictions. It is particularly effective at capturing nonlinear relationships between environmental factors and erosion susceptibility (Chen & Guestrin, 2016a). XGBoost, also known as Extreme Gradient Boosting, is rooted in optimization methodologies, particularly gradient-based techniques and decision trees. The training process of XGBoost entails minimizing an objective function, which integrates a loss component and regularization factors. The training of XGBoost is an iterative procedure. A new tree is added at each iteration in order to correct the errors generated by the previous trees. The ultimate prediction of the model is computed by aggregating the predictions of all trees, with weights determined by the learning rate. It can model complex interactions among predictor variables, including topography, soil characteristics, land use/land cover, and hydrological processes, which influence gully erosion dynamics. XGBoost incorporates regularization techniques such as L1 and L2 regularization (also known as Lasso and Ridge regularization) to prevent overfitting (Chen & Guestrin, 2016b). This helps ensure that the model generalizes well to unseen data and avoids memorizing noise in the training dataset. It provides control parameters for tree pruning and constraints on tree complexity, such as maximum tree depth and minimum child weight. These parameters allow researchers to fine-tune the structure of individual decision trees, balancing model complexity with predictive accuracy (Chen & Guestrin, 2016a). In gully erosion susceptibility mapping, the occurrence of erosion events may be relatively rare compared to non-erosion areas, leading to class imbalance. XGBoost offers options to address 644 O.E. Abiodun et al. / Nigerian Research Journal of Engineering and Environmental Sciences 10(2) 2025 pp. 639-654 class imbalance, such as adjusting class weights or using evaluation metrics like area under the ROC curve (AUC) that are robust to imbalanced data (Chen & Guestrin, 2016b). XGBoost provides insights into feature importance that enable researchers to pinpoint the environmental factors that have the most impact on gully erosion susceptibility. This feature importance analysis helps prioritize mitigation efforts and inform land management strategies by highlighting key variables contributing to erosion occurrence (Chen & Guestrin, 2016a). XGBoost is a flexible and effective Machine Learning method for gully erosion susceptibility mapping. Its ability to handle complex relationships, prevent over-fitting, provide insights into feature importance, and scale efficiently makes it well-suited for accurately predicting erosion occurrence and informing land management decisions in the study area. 2.4. Model Validation 2.4.1. Confusion matrix A confusion matrix is a certain table arrangement that allows one to see the performance of an algorithm, usually one for supervised learning. The four different prediction types are distinguished by the confusion matrix. They can be explained in terms of a binary classification of gully erosion. • True positive (TP): This represents a case where the positive class was accurately predicted, denoting that an area or occurrence was designated as "vulnerable". • True negative (TN): This indicates that a place or event that is categorized as "not vulnerable" is in fact in the "not vulnerable" class, as the prediction successfully identified the negative class in this case. • False positive (FP): This indicates that a place or event that is labeled as "vulnerable" may not actually be in the "vulnerable" class, as the prediction may have identified the positive class incorrectly in this case. • False negative (FN): This indicates that a place or event that is categorized as "not vulnerable" is actually in the "vulnerable" class when the prediction incorrectly recognizes the negative class. 2.4.2. Accuracy The percentage of right predictions is known as accuracy. The confusion matrix counts the instances where class 0 items are labeled as class 1 and vice versa. Equation 4 provides the formula to calculate binary classification accuracy: Accuracy = 𝑁𝑢𝑚𝑏𝑒𝑟 𝑜𝑓 𝐶𝑜𝑟𝑟𝑒𝑐𝑡 𝑃𝑟𝑒𝑑𝑖𝑐𝑡𝑖𝑜𝑛𝑠 𝑇𝑜𝑡𝑎𝑙 𝑛𝑢𝑚𝑏𝑒𝑟 𝑜𝑓 𝑝𝑟𝑒𝑑𝑖𝑐𝑡𝑖𝑜𝑛𝑠 = 𝑇𝑃+𝑇𝑁 𝑇𝑃+𝐹𝑃+𝑇𝑁+𝐹𝑁 (4) Precision Precision refers to the accuracy of the positive predictions, indicating the fraction of positive predictions that were actually positive. This can be determined by Equation 5. Precision = 𝑇𝑃 (𝑇𝑃+𝐹𝑃) (5) Recall The percentage of positive cases that the classifier correctly detects is known as recall and this can be determined using Equation 6 Recall = 𝑇𝑃 (𝑇𝑃+𝐹𝑁) (6) F1-Score The F1 score is designed to simultaneously work effectively with unbalanced data. Equation 7 can be used to determine F1 Score. F1 Score = 2 x 1 (1 𝑃𝑟𝑒𝑐𝑖𝑠𝑖𝑜𝑛+1 𝑅𝑒𝑐𝑎𝑙𝑙) (7) 645 O.E. Abiodun et al. / Nigerian Research Journal of Engineering and Environmental Sciences 10(2) 2025 pp. 639-654 Area under curve The Area Under Curve (AUC) was computed using the ROC curve as a guide. Equations 8 and 9 were used to calculate the false positive rate (1-specificity) and true positive rate (recall, sensitivity): ROC = 𝑆𝑒𝑛𝑠𝑖𝑡𝑖𝑣𝑖𝑡𝑦 1−𝑠𝑝𝑒𝑐𝑖𝑓𝑖𝑐𝑖𝑡𝑦 (8) 1-specificity = 𝐹𝑃 (𝐹𝑃+𝑇𝑁) (9) 3. RESULTS AND DISCUSSION The remotely sensed data was used to extract the input variables. Each feature was obtained independently through several spatial analyses. 3.1. Rainfall The rainfall map was created, and it is shown in Figure 2. The analysis revealed distinct rainfall patterns categorized into five classes based on intensity ranges. Areas experiencing moderate to high rainfall levels (351.001 - 406 mm) demonstrated heightened susceptibility to gully erosion, likely due to increased runoff and soil saturation. Regions with lower rainfall intensity (336 - 351 mm) exhibited comparatively lower erosion rates. 3.2. Elevation According to Figure 3, the Digital Elevation Model (DEM) analysis categorized terrain characteristics into five classes based on elevation ranges. Areas with higher elevation levels (148.577 - 193 meters) exhibited increased susceptibility to gully erosion, often associated with rugged terrain features and steep slopes. Regions with moderate elevation (112.946 - 148.576 meters) showed intermediate erosion susceptibility depending on slope gradient and soil composition. Lower elevation areas (75 - 112.945 meters) generally displayed slower erosion rates due to flatter terrain. Figure 2: Rainfall Map of Ihioma Community Figure 3: Elevation Map of Ihioma Community 3.3. Slope The analysis of terrain slope as shown in Figure 4 categorized slope values into five classes based on slope gradients. Areas with higher slope values (13.445 - 24.333) indicated steep terrain, potentially experiencing higher erosion rates. Regions with moderate slope values (7.349 - 13.455) exhibited varied slope gradients, 646 O.E. Abiodun et al. / Nigerian Research Journal of Engineering and Environmental Sciences 10(2) 2025 pp. 639-654 influencing erosion susceptibility. Lower slope values (4.486 - 7.348) corresponded to areas with intermediate slope gradients. Areas with slope values ranging from (2.387 - 4.485) represented gentle slopes. The lowest slope values (0 - 2.386) indicated flat terrain, potentially facing lower erosion risk. 3.4. Topographic Wetness Index (TWI) The analysis of the topographic wetness index categorized TWI values into five classes based on wetness gradients (Figure 5). Areas with higher TWI values (17.072 - 20.332) indicated higher levels of topographic wetness, potentially leading to increased erosion rates. Regions with moderate TWI values (13.810 - 17.071) exhibited varied wetness gradients, influencing erosion susceptibility. Lower TWI values (12.143 - 13.809) corresponded to areas with intermediate wetness levels. Areas with TWI values ranging from (10.259 - 12.142) represented relatively lower wetness levels. The lowest TWI values (1.85 - 10.258) indicated drier terrain, potentially facing lower erosion risk. Figure 4: Slope Map of Ihioma Community Figure 5: Topographic Wetness Index map of Ihioma Community 3.5. Topographic Position Index (TPI) The topographic position index analysis categorized topographic position into five classes based on TPI values (Figure 6). Areas with higher TPI values (0.967 - 5.171) indicated elevated topographic positions, such as ridges, potentially experiencing lower erosion rates. Regions with moderate TPI values (0.051 - 0.966) exhibited varied topographic positions, influencing erosion susceptibility. Lower TPI values (-0.811 - 0.05) corresponded to areas with relatively flat terrain. Areas with TPI values ranging from (-3.399 - - 0.812) represented slight depressions or low-lying areas. The lowest TPI values (-8.574 - -3.4) indicated valley bottoms or depressions, potentially facing increased erosion risk. 3.6. Stream Power Index (SPI) The stream power index analysis categorized SPI values into five classes based on stream power gradients. Areas with higher SPI values (4,855,627,990.313 - 7,444,975,616) indicated higher stream power, potentially leading to increased erosion rates and gully formation (Figure 7). Regions with moderately high SPI values (2,380,149,771.912 - 4,855,627,990.312) exhibited significant stream power, influencing erosion susceptibility. Lower SPI values (13,533,428.479 - 2,380,149,771.911) corresponded to areas with intermediate stream power levels. Areas with SPI values ranging from (-2,249,008,360.566 - 13,533,428.478) represented relatively normal stream power conditions. The lowest SPI values (- 4,615,624,704 - -2,249,008,360.567) indicated areas with minimal stream power, potentially facing lower erosion risk. 647 O.E. Abiodun et al. / Nigerian Research Journal of Engineering and Environmental Sciences 10(2) 2025 pp. 639-654 Figure 6: Topographic power index map of Ihioma Community Figure 7: Stream Power Index map of Ihioma Community Figure 8: Normalized Difference Vegetation Index map of Ihioma community Figure 9: Land Cover map of Ihioma Community 3.7. Normalized Difference Vegetative Index (NDVI) The normalized difference vegetation index analysis categorized vegetation cover into five classes based on NDVI values (Figure 8). Areas with higher NDVI values (0.399 - 0.525) indicated denser vegetation cover and potentially lower erosion susceptibility. Regions with moderate NDVI values (0.253 - 0.398) exhibited intermediate vegetation cover and erosion rates. Lower NDVI values (0.028 - 0.252) corresponded to areas with sparse vegetation cover, potentially leading to increased erosion risk. Integrating NDVI data with other factors such as rainfall and topography enhances erosion prediction models, facilitating targeted erosion control strategies. 654 O.E. Abiodun et al. / Nigerian Research Journal of Engineering and Environmental Sciences 10(2) 2025 pp. 639-654 Roy, J., and Saha, S. (2019). Landslide susceptibility mapping using knowledge driven statistical models in Darjeeling District, West Bengal, India. Geoenvironmental Disasters, 6(1), . DOI: 10.1186/s40677-019-0126-8 Thakuriah, G. (2023). GIS-based revised universal soil loss equation for estimating annual soil erosion: a case of lower Kulsi basin, India. SN Applied. Science. 5(3) https://doi.org/10.1007/s42452-023-05303-0 Torri, D., Poesen, J., Borselli, L., Bryan, R., and Rossi, M. (2012). Spatial variation of bed roughness in eroding rills and gullies. CATENA , 90(3), pp. 76-86. DOI:10.1016/j.catena.2011.10.004 United Nations Convention to Combat Desertification (UNCCD), (2019). The Global Land Outlook, West Africa Thematic Report, Bonn, Germany Zabihi, M., Mirchooli, F., Motevalli, A., Darvishan, A. K., Pourghasemi, H. R., Zakeri, M. A., and Sadighi, F. (2018). Spatial modelling of gully erosion in Mazandaran Province, northern Iran. CATENA, 161, pp. 1-13. https://doi.org/10.1016/j.catena.2017.10.010.