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 30m x 30m 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 8TM 30m 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 (Ouyanget al2017) 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 (Oh et al. 2017),
FR = (A∕B) ∕ (C∕D)
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(Poudyalet 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 (Bathrelloset 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 (Saaty1980). 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
CR is expressed as follows:
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. 2000) and more gullies were identified from high-resolution imagery from Google Earth. Fig.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.
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; Gutie´rrez et al. 2009a,). 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; Confortiet al. 2011)
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, O. H, Igwe, O. (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.5km, 0.5–1km, 1–2km, 2–3km, 3-5km.
Soil’s physical properties are a major contributing factor to runoff, soil infiltration, gully occurrence, and soil resistance to erosion (Xia, D., 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 (Chaplotet 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,
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 (Rahmati et al. 2015b; Ghorbani Nejad et al. 2017)
Through supervised classification, the landuse map (Fig.7d) was derived from Landsat 8TM 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 erosion-prone (Gutie´rrez et al. 2009).
The type of vegetation and the nature of precipitation are influenced by altitude. DEM was used to obtain the elevation of the area in the ArcGIS environment, five classes were obtained (Fig.7e): 330 - 414m, 414 - 477m, 477 - 544m, 544 - 611m, 611 - 721m
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).Go´mez-Gutie´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
The lithological properties of the exposed geological formation affect the occurrence of gullies (Golestaniet 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
Paved roads effectively concentrate surface runoff. Road construction and expansion projects tend to encourage slope instability. Five classes of the road buffer map corresponding to buffer distances of 0-1km, 1-2km, 2-3km, 3-6km, 6-9km (Fig.7h).
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)
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.