GIS-based gully erosion susceptibility modeling, adapting bivariate statistical method and AHP approach in Gombe town and environs Northeast Nigeria

Gully erosion is a major environmental problem in Gombe town, a large area of land is becoming unsuitable for human settlement, hence the need for a gully erosion susceptibility map of the study area. To generate a gully inventory map, a detailed field exercise was carried out, during this investigation one hundred gullies were identified and studied extensively within the study area of about 550 km2. In addition to the mapped gullies, Google EarthPro with high-resolution imagery was used to locate the spatial extents of fifty (50) more gullies. Ten gully erosion predisposing factors were carefully selected considering the information obtained from literature, and multiple field survey of the study area, the factors include elevation, slope angle, curvature, aspect, topographic wetness index (TWI), soil texture, geology, drainage buffer, road buffer and landuse. In this study, a GIS-based Frequency Ratio (FR) and Analytical Hierarchy Process (AHP) models were employed to predict areas prone to gully erosion in Gombe town and environs. The result obtained from FR shows that drainage, soil texture, and slope have the highest correlation with gully occurrence, while the AHP model revealed that drainage buffer, soil texture, geology have a high correlation with the formation of a gully. Gully erosion susceptibility maps (GESM) were produced and reclassified into very high, high, moderate, and low zones. The overall accuracies of both models were tested utilizing area under the curve (AUC) values and gully density distribution.FR and AHP model have AUC values of 0.73 and 0.72 respectively, the outcome indicates that both models have high prediction accuracy. The gully erosion density distribution values revealed that gullies are concentrated in the very high susceptibility class and it decreases towards the low class, therefore the GESM produced using these models in this study area is reliable and can be used for land management and future planning.


Introduction
One of the factors that endanger water and soil is soil erosion (Magliulo 2012). The major cause of land degradation around the world is soil erosion by gully erosion (Nampak et al. 2018;Rizeei et al. 2016). Gully erosion is defined as a deep channel that is formed by concentrated water flow, which in the process removes surface soils and materials (Kirkby and Bracken 2009).
The changes in the quantity of moisture content resulting from dry and wet seasons are a major factor contributing to cracks and grooves in clay formations consequently forming rilled erosion and gullies (Torri et al. 2012). Runoff accumulates within these cracks at the first sudden rainfall and consequently, gully develops. Under natural conditions, vegetation helps to hold soil in place and protect it from the direct impact of rainfall thereby controlling run-off. The effect caused by excessive clearing, inappropriate landuse, and compaction of the soil caused by grazing is that the soil becomes exposed and unable to absorb excess water, this implies that there becomes an elevated level of Surface run-off that concentrates in drainage lines, making it possible for the occurrence of gullies in susceptible areas, (Ligonja and Shrestha 2015). The threat of gully erosion on the socio-economic development of Gombe state ( Fig. 1) includes the destruction of houses, loss of lives, displacement of people, land depreciation, destruction of roads and culverts (Mbaya 2017). To better identify areas prone to gully erosion and realize the mechanism involved it becomes expedient to invent a gully erosion susceptibility map (GESM) (Arabameri et al. 2018). The demarcation of an area into zones of varying susceptibility is made possible by the estimation of the input of each conditioning factor. Several models have been created and employed to analyze both quantitatively and qualitatively, the rate of gully erosion. Lately, researchers have executed statistical and machine learning methods in other to determine the statistical relationship that exists between gully erosion factors and the spatial distribution of gullies (Magliulo 2012). Several researchers have assessed control measures and the effect of gully erosion in Gombe town, but no work has been done to produce a reliable gully erosion susceptibility map (GESM) adapting Frequency Ratio (FR) and Analytical Hierarchy Process (AHP) models. The objective of this research includes: To develop a gully inventory map of the study area, secondly, to produce all relevant thematic maps for the conditioning factors and establish the correlation between predisposing factors and gully occurrence, to apply Frequency Ratio and Analytical Hierarchy Process (AHP) models to produce GESM and finally, to validate the produced susceptibility maps using AUC and gully density distribution.

Climate and vegetation
The area is characterized by wet and dry seasons, having a mean annual rainfall and temperature of 850 mm and 32°C respectively. Rainfall within the study area occurs mostly between June and September. Precipitation is associated with a storm of high intensity, especially in July and August. The vegetation of the Gombe area can be described as Sudan savannah with open grassland and shrubs which dries up during the dry season, the vegetation comprises of scattered shrubs and trees.

Geology
The study area is located in the north-south trending Gongola basin of the upper Benue Trough, it is underlain by four geologic Formations, they include, Yolde Formation, Pindiga Formation, Gombe Sandstone, and Kerri-Kerri Formation. (Fig. 1) The Yolde Formation is Cenomanian in age, deposited in continental to the marine environment it comprises of sandstone and shale at the base while the top consists of sandstones, shales, and calcareous sandstone (Abubakar et al. 2008). The Cenomanian-Santonian marine Pindiga Formation, consisting of thick marine shale, with some limestone beds toward the base (Zaborski et al. 1997).while the Gombe sandstone is the last cretaceous deposit in the Gongola arm of the Upper Benue trough. It is underlain by the marine Pindiga Formation and overlain by the Kerri-Kerri Formation. The lacustrine to deltaic Gombe Formation consists of well-bedded fine to medium-grained friable, ferruginous sandstone, siltstone, and shale with ironstone (Orazulike 1988). Three major lithofacies characterize this formation, they include the basal transitional portion, bedded facies, and red sandstone facies. At its base, it comprises intercalation of silty shale occasionally, with plant remains (Zaborski et al. 1997). The continental Kerri-Kerri Formation ended sediment deposition in the Upper Benue Trough; it consists of sandstones, siltstones, and shales.

Geomorphology
The complex geologic crystalline bedrock formed the base on which the relief of Gombe was established. The elevation of the study area is within 330 m to 721 m (Fig. 2). The flat-topped to conical hills characterized the landscape; this landscape is the outcome of dissection and stream incision in the area after the deposition of sedimentary formation during the late cretaceous period. Aside from the Gombe hill and Liji hill (Fig. 3.) Gombe town is generally taken to be a flat area (Arabi et al. 2009).

Methodology
The methodology employed in this research is summarized in the flowchart in Fig. 4. The first approach to this study to generate a reliable gully inventory map was to embark on a detailed field survey. Ten gully erosion predisposing factors were carefully selected considering information obtained from the literature and field survey of the study area. Next was to produce thematic maps corresponding with the chosen conditioning factors. ASTER DEM with a resolution of 30 m × 30 m was used to extract Topographic related factors such as elevation, slope angle, curvature, aspect, topographic wetness index. The geological map was digitized from a previously existing geological map produced by Nigeria  Geological Survey Agency (NGSA) while soil map obtained from the Institute for Agricultural Research; Ahmadu Bello University Zaria Nigeria was obtained and digitized to produce the soil texture map of the area. Drainage buffer and road buffer were digitized from Google Earth imagery while Landsat 8™ 30 m resolution was employed to produce the landuse map in ArcGIS version 10.4. The modeling phase required that all the gully erosion training (75%) and testing (25%) data were selected randomly and overlain over all the produced thematic maps in the ArcGIS environment. After the rasterization of all produced thematic maps, in other to ensure uniformity of the areal extent and resolution, the thematic maps were masked with Digital Elevation Model (DEM) of a 30mx30m grid size. The prediction rate acquired from the Frequency Ratio approach together with the criteria weights from the AHP were subsequently integrated into the ArcGIS raster calculator, at this point each conditioning factor map was multiplied with their respective prediction rates and criteria weights, at the end of these process the susceptibility index map was produced. The Gully Erosion Susceptibility Index Maps generated were reclassified into four zones of different degrees of susceptibility to gully occurrence. To validate the accuracy of the prediction, the area under the curve (AUC) and gully density distribution technique was used. The area of gully erosions in a susceptibility zone was divided by the area of the zone to determine the gully erosion density distribution.

Frequency ratio (FR)
Frequency ratio (FR) is the ratio of the area where gullies occurred in the study area to the whole study area. For a given factor, it is the ratio of the possibility of gully occurrence to its nonoccurrence (Lee and Talib 2005). The frequency ratio model is a bivariate statistical method (Ouyang et al. 2017) employed for calculating the possible relationship between gully erosion and conditioning factors. It is necessary to have several independent variables that are required to determine the dependent variables (Abedini et al. 2017). According to , A represents the number of pixels within gully erosion for each class of geo-environmental factors, and B represents the total number of gullies in the study area, while C represents the total number of pixels within each class of the predisposing factors, and finally, D represents the total number of pixels in the study area.
The Frequency Ratio correlation value varies between < 1 to > 1. When FR value is observed to be less than 1, it indicates a low correlation between gully location and each class of the predisposing factor, on the other hand, a higher correlation exists if the FR value is observed to be higher than 1 (Poudyal et al. 2010).

Analytic hierarchy process (AHP)
AHP was introduced by Saaty (1980) and is considered as a multi-criteria and multi-objective method, in this method, weights, which are decisions made are assigned based on the knowledge and experience of an expert. These weights are assigned in a form of pairwise relative comparison (Bathrellos et al. 2017;Papadakis and Karimalis 2017).
In the Analytical Hierarchy Process (AHP) the numerical value for each conditioning factor must be between 1 and 9 in the comparison matrix as represented in (Table 1). In the matrix of this study, the conditioning factors for gully occurrence were organized hierarchically. A numerical value was allocated to each predisposing factor, based on their importance as compared with one another (Saaty 1980). The Eigenvalue (λmax) and the Consistency Ratio (CR) were gotten after calculations were executed using the average of the hierarchically arranged factors according to Saaty (1980), in other to achieve a consistent comparison matrix the Eigenvalue (λmax) and the total number of factors (n) have to be the same.
Mathematically Consistency Index is expressed as.
CI is the consistency index, λmax is the Eigenvalue, and n represents the total number of factors being compared. Consistency Ratio (CR) is used to check for the consistency of the comparison matrix. Table 1 The scale of preference between two factors in AHP (Saaty 1980) Preference factor Degree of preference Definition 1 Equally two factors have equal importance to the desired objective 3

Moderately
One factor slightly has more importance than the other to the desired objective 5 Strongly A factor strongly has higher importance than the other towards achieving the desired objective 7 Very strongly Here a factor has very strong importance more than the other towards achieving its objective 9 Extremely One factor is extremely more important than the other toward achieving the desired objective 2,4,6,8 Intermediate when factor importance is between 1,3,4,7 and 9 CR is expressed as follows: CR ¼ CI RI CR represents the Consistency Ratio, CI represents the Consistency Index, and RI is the Random consistency Index of the pairwise comparison matrix. Table 2 shows the Random consistency index. The rule for the consistency index states that a Consistency Ratio (CR) that is equal to or less than 0.1 signifies that the matrix is acceptable; while a ratio of more than 0.1 means that the matrix is not consistent (Saaty 1980). The susceptibility index map was finally generated for the AHP model in the ArcGIS environment integrating the Arc Map raster calculator, the gully susceptibility index map generated was subsequently divided into low, moderate, high, and very high susceptibility zones according to the natural break classification method (Tian et al. 2017).

Gully erosion inventory
To generate a comprehensive and dependable gully inventory map, as shown in Fig. 5, a detailed field survey was done, (Guzzetti et al. 2002) and more gullies were identified from high-resolution imagery from Google Earth. Figure 6 shows some photographs of some gullies in the study area. Next was the conversion of the gully areas represented as polygon to points and they were also overlain on a hill shade view of the area. Among all the 150 mapped gullies, 75% (113 gullies) of the mapped gullies in the study area were used for training. The training data was used in the modeling phases for the spatial prediction of gullies susceptible zones in the study area. While 25% (37 gullies) of the mapped gullies were used as testing data (Chung and Fabbri 2003;Dube et al. 2014) this was utilized in the validation of the gully erosion susceptibility models generated.

Conditioning factors
Quite a several factors contribute to gully erosion as discussed in the literature (Dube et al. 2014), the occurrence  of gully erosion and its behavior depends on the following factors: topography characteristics, climate, lithology, soil characteristics, and land use (Poesen et al. 2003;Gutiérrez et al. 2009). The erodibility of the geologic formation in an area and the erosivity of surface runoff determine how susceptible the environment is to the formation of the gully (Conoscenti et al. 2008;Conforti et al. 2011).

Drainage buffer
The movement of eroded sediments is usually facilitated by drainage, therefore it could be said that gullies are associated with drainage (Conoscenti et al. 2014). The erosive actions in the drainage channels consequently reduce the shear strength of the slope material thereby increasing their chances of failure by sliding (Ozioko and Igwe 2020). The distance from the drainage network was calculated in ArcGIS 10.4.1. Five buffer zones were created within the study area to determine the effect of drainages on gullies occurrence (Fig. 7a), they include 0-0.5 km, 0.5-1 km, 1-2 km, 2-3 km, 3-5 km.

Soil texture
Soil's physical properties are a major contributing factor to runoff, soil infiltration, gully occurrence, and soil resistance to erosion (Xia et al. 2015). The occurrence of subsurface flow and piping is a function of the soil texture, the formation of gullies follows the collapse of the top of these pipes Chaplot et al. 2005). It is expedient that the soil's texture is put into consideration in other to assess gully erosion susceptibility (Conoscenti et al. 2013). Four soil texture types were identified in the study areas (Fig. 7b), they include sandy loam, stony, sandy clay, sandy,

Slope degree
ASTER DEM was used to obtain the slope information of the area the slope of the area was divided into five cases: 2-4,4-9, 9-16, 16-43 (Fig. 7c). Surface flow accumulation is common in flat areas and consequently initiation the formation of the gully Ghorbani Nejad et al. 2017).

Land use
Through supervised classification, the landuse map ( Fig. 7d) was derived from Landsat 8™ Main landuse type in the study area include urban area, bare surface, shrubs and grasses, and farmlands. The stability of slope and gully occurrence is influenced by the landuse type (Zakerinejad and Maerker 2015). Infrastructural development such as roads, housing, and industries, prevents infiltration and increases the amount of surface runoff; this runoff concentrates and consequently produces gullies. Generally, bare surfaces and low vegetated areas are gully erosionprone (Gutiérrez et al. 2009).

Elevation
The type

TWI
The topographic wetness index (TWI) which is a secondary topographic factor within the runoff model is usually employed to determine the extent to which topography control hydrological processes. TWI estimates the possibility of water accumulating in soil due to slope and upstream catchment area (Rahmati et al. 2016). Gómez- Gutiérrez et al. (2015), noted that TWI has to be considered in evaluating gully erosion susceptibility. TWI was derived from DEM. The TWI map was divided into four classes (Fig. 7f), they include: < 7, 7.0-8.5, 8.5-10.6, > 10.6.

Lithology
The lithological properties of the exposed geological formation affect the occurrence of gullies (Golestani et al. 2014). The geology thematic map (Fig. 7g) represents the distribution of the gully erosion within the study area for the underlying lithology. The lithology thematic map was digitized from a geologic map produced by Nigeria Geological Survey Agency (NGSA) and was classified into four, according to the lithologic unit of the area; they include Kerri-Kerri Formation, Gombe Sandstone, Yolde Formation, and Pindiga Formation.

Slope aspect
The slope aspect has major control over the vegetation type because it influences the duration of sunlight. The slope aspect affects transpiration and moisture content also, evaporation, and indirectly affects the erosion process (Jaafari et al. 2014). Nine classes corresponding to flat, north, northeast, east, southeast, south, southwest, west and northwest are represented in the aspect map (Fig. 7i).

Curvature
Extensive analysis of profile curvature reveals useful geomorphological information (Davoodi Moghaddam et al. 2015). Three classes of the profile curvature (Fig. 7j) include, convex, flat, and concave.

Results and discussion
Prediction rate (frequency ratio) The prediction rate result shows that drainage buffer, soil texture, slope degree, and landuse, have the major contribution to gully occurrence with a prediction rate of 4.21, 3.25, 2.49, 2.27, respectively as shown in (

Interpretation of frequency ratio values for each factor class
The Frequency Ratio correlation value varies between < 1 to > 1, when FR value is observed to be less than 1, it indicates a low correlation between gully location and each class of the predisposing factor, on the other hand, a higher correlation exists if the FR value is observed to be higher than 1 (Poudyal et al. 2010). Analysis of result obtained from the drainage buffer shows that the 0-500 m class has the highest correlation with FR of (1.87) (Table 4), others include 500 m-1 km (0.553), 1 km-2 km (0.038), therefore it can be deduced that the probability of gully occurrence decreases as the distance from drainage increases.
For the soil texture, the FR value of 2.29, was observed for the sandy clay class while other classes have FR value < 1 (Table 5). There is a dominance of sand in all classes, which could be attributed to the high erodibility since sand has no cohesion. This finding is in agreement with Mbaya et al. (2012) who stated that the dominance of the sand portion in Nigeria savanna soils has accelerated the occurrence of gully erosion. The high sand content of the soil also indicates that infiltration and percolation are high, thereby bringing about a piping phenomenon and ceiling collapse caused by subsurface flow, this is also in sync with the outcome of work done by Mbaya et al. (2012). The presence of clay in the soil may have contributed to the formation of a gully since clay exhibits shrink and swell characteristics, this could result in the formation of cracks during the alternating warm and dry seasons. These cracks at the first sudden rainfall concentrate the runoff and therefore gully erosion is formed. This is similar to findings by McCloskey et al. (2016).
Assessment of FR for the relationship between slope degree and gully occurrence reveals that slope degree class 0-2 0 is the only class with FR value > 1, (1.25), indicating a high correlation. While slope degree greater than 2 all have FR value < 1 (Table 6). Golestani et al. (2014) proved that areas with gentle slope are susceptible to surface flow accumulation and gully erosion.
The landuse variable shows that urban class and farmland class have the FR value > 1(2.39, 1.31) respectively, while other classes have values less than 1, ( Table 7). The high FR value in the urban class could be attributed to the infrastructural development such as roads, housing, industries, which have sealed the ground surface thereby preventing infiltration and increasing the amount of surface runoff, this runoff concentrates consequently forming gullies. This land-use type also deprives the soil of vegetation cover which is supposed to hold the soils together and prevent the direct impact of a raindrop on the soil. The high FR value owned by the farmland class could be attributed to agricultural activities such as plowing, compaction by farm machines and overgrazing, etc. which might have caused the soil to lose its structure and cohesiveness and becomes easily eroded.
Analysis of the frequency ratio for TWI shows a strong correlation for factor class> 10.6 and 8.5-10.6, both having FR value > 1 (Table 9). It can be deduced that classes with high TWI value are pointers to gully occurrence. This result is in agreement with the finding of Rahmati et al. (2016) and Dube et al. (2014).
Assessment of geology shows that the Yolde with units of sandstone and shale and Kerri-Kerri Formation class consisting of sandstone, shale, and clay, have a major contribution with FR value of 1.68 and 1.26 respectively (Table 10). Pindiga formation and Gombe sandstone have FR value < 1.
Road buffer reveals that Factor class < 1 km and 1-2 km have a strong correlation with gully occurrence, with FR value of 1.119 and 1.079. Table 11, while other classes have FR value less than 1. This suggests that the probability of gully occurrence decreases with distance from the road.
The result from aspect shows that the following classes have a strong relation with gully occurrence: flat face, east, southeast, and north with the following FR values 3.70, 1.68, 1.20, and 1.24, Table 12. Other classes have FR value less than 1.
Analysis of the profile curvature reveals that the concave class has a major contribution to gully occurrence with FR of 1.45 Table 13. The flat class and convex class have values <1indicating a very low correlation with FR value of 0.89 and 0.87 respectively.

Gully susceptibility maps
The respective weights from the Frequency ratio and AHP were multiplied by each predisposing factor map. This method was employed to produce the area's gully erosion susceptibility index map. The GESM shows a range of susceptibility to the study area. The produced susceptibility index maps for the frequency ratio and AHP models were reclassified into four zones according to the quantile classification scheme ( Fig. 8a and b). That is low, medium, high, and very high susceptibility zones (Youssef et al. 2015).

Gully density distribution
Prediction rate curves coupled with gully density distribution within the susceptible areas were employed to check for the validity of the susceptibility map. Gully erosion density distribution for frequency ratio and AHP reveals that gullies are concentrated in the very high susceptibility class and it decreases as it gets to the low susceptibility class, density values for FR include: 0.729, 0.223, 0.097, 0.005 (Table 16), and AHP include: 0.74,10.23, 0.10, 0.006   (Table 17). Gupta et al. (2008) stated in his work that the distribution of density decreases from the high susceptibility class to low class, which is in agreement with this study.
The area under the curve (AUC) To determine the accuracy of the Frequency ratio and AHP model applied, the (AUC) area under the curve was employed (Mohammady et al. 2012); AUC analysis is an efficient technique for evaluating the correctness of a test (Razandi et al. 2015). This curve indicates the accuracy and reliability of a predicting system. Rahmati et al. (2016) classified the AUC values as: 0.5-0.6, poor; 0.6-0.7, average; 0.7-0.8, good; 0.8-0.9, very good; and 0.9-1, excellent. For this study, the FR and AHP model gave AUC values of 0.73 and 0.72 respectively (Fig. 9). The outcome indicates that both models have high prediction accuracy with the Frequency ratio model slightly higher than the AHP model. Therefore the GESM produced using these models in this study area is reliable and can be used for land management and future planning.

Conclusion
It is important to identify areas susceptible to gully erosion and to generate susceptibility maps. A comprehensive and dependable gully inventory map was successfully developed through detailed field surveys and from high-resolution Google Earth imagery. The spatial relationship between gully occurrences and its ten carefully selected causative factors which include: drainage buffer, soil texture, slope degree, land use, elevation, TWI, Lithology, Road buffer, and aspect were assessed using FR and AHP model, the FR model revealed that drainage buffer, soil texture, slope degree, and landuse, have the major contribution to gully occurrence and according to AHP model drainage buffer, soil texture, geology, and landcover are the most important contributors. These two methods were employed to generate two gully erosion susceptibility maps (GESM) for the study area. The GESM shows a range of susceptibility in the study area. The gully erosion susceptibility index maps generated from the FR and AHP models were reclassified into four zones according to the quantile classification scheme.  The performances of the models were subsequently validated using the AUC plot and gully density distribution. AUC values of 0.73 and 0.72 were obtained for FR and AHP models respectively. The obtained values indicated that both models have high prediction accuracy, though that of the FR is slightly higher than AHP. In addition to the AUC, the gully erosion density distribution for both methods revealed that gullies are concentrated in the very high susceptibility class and it decreases as it gets to the low susceptibility class, this outcome is as expected. The high accuracy obtained from AUC plots coupled with the alignment of the gully density distribution with the varying susceptibility classes proves that the GESM produced is reliable and accurate. The GESM produced shall be of enormous importance for town planners and decisionmakers when there is a need for future infrastructural developments.