Study area
Geographically, the study area lies within the Southeastern part of Anambra state, Nigeria. The population density of this area is approximately one million (National Bureau of Statistics, 2020), which are scattered over many towns within the area, notably Ekwulubia, Agulu, Igboukwu, Nanka, Aguata, amongst others. Geomorphologically, it is characterized by undulating topography with the hilly areas standing at elevation of 400 m while the valley areas lies as low as 60 m above sea level (Fig. 1). The boundary between the high and the low areas is characterized by high angled slopes which have been steepened by geologic hazards, such as landslides and gullying (Egboka et al., 1990).
Tropical climate prevails within the study area, which is characterized by relatively high temperature throughout the year that ranges from 16 0C in the month of December to about 31 0C in the month of March. The total annual rainfall is about 1580 mm, with minimum rainfall of about 16 mm occurring in February and a maximum of about 350 mm in July (Eze, 2007). Tropical rain forest predominates in most part of the study area, apart from few high elevated areas that are comprised of shrubs and grasses.
Geologically, the study area lies in the Anambra basin which is part of the lower Benue Trough. The history and evolution of this Basin has been discussed extensively by several authors, as well as it’s sedimentological and stratigraphical composition (Nwajide, 1979; Oboh-Ikuenobe et al., 2005; Dim, 2021). Almost the entirety of the study area is been underlain by Ameki Formation which was deposited during the Eocene (Nwajide, 1979). This Formation comprises of successive layers of fine to coarse grained tidally influenced and fluvial sandstones at the basal part. The course grained sandstone is successively overlain by intercalations of thin layers of clay, shale, and limestone, with cross-bedded coarse grained sandstones and clays at the uppermost layer (Arua and Rao, 1987). Due to intensive tropical weathering, most part of the outcropping rock units of this Formation have been severely weathered into lateritic soils. Two major soil type, deep porous red soil derived from sandy deposits, and reddish brown soils derived from sandstone and shale deposits covered more than 95 percent of the study area. The remaining five percent of the area which is seen at the edges and at some low lying terrains of the study area is been underlain by dark brownish soils derived from shale, likely from the outcropping adjacent Imo Formation (Fig. 2).
Desk study
Prior to field sampling, the digital elevation model (DEM) data of the study area was acquired from the United States Geological Survey online archives, whereas the population density data was acquired from NASA Socioeconomic Data and Applications Center (SEDAC, 2020). With the help of surfer 11, the DEM data was analysed and some geomorphologic parameters such as elevation and slope were extracted. The extracted information was used in producing 2D and 3D contour maps which helped to identify potential gully erosion channels. The potential gully channels that were indentified from the DEM data were verified by visual confirmation using Google earth pro in order to ensure that the channels were neither man-made drainage systems nor stream channels. This helped in narrowing down the specific sites visited during field investigation and sampling.
Field investigation
At the gully sites, the depth and width of some of the gullies were measured using a long ranging pole and a measuring tape. The width and depth of some of the gullies that are relatively very wide and deep, respectively were estimated. In areas where the gullies could not be accessed, the width was measured with the help of the ruler function in Google earth pro. Because it is an arduous task to manually trace and measure the length of the gullies, their lengths were estimated using the ruler function in both surfer 11 and Google earth pro. Altogether, about 120 gullies were captured from both field investigation and through image analysis. Information from 60 gully erosion event scars were used to produce the gully hazard map, which is a 2D representation of the studied gully erosion sites on map. The remainder of the gully location points was used for model validation.
Soil Sampling
A total of sixty (60) disturbed soil samples were collected at different locations within the study area. Samples were collected considering the elevation (topographic) contour and the frequency of gully erosion. More samples (49) were collected in areas with high frequency of gully erosion, compared to 11 samples which were collected in areas with little or no presence of gully erosion. However, it was ensured that samples were collected at each contour interval which was set at 50 m. Topographic contour interval was considered because it was assumed that contour variation may represent change in physical properties of the underlying soil. Soil sampling was done by digging horizontally with a hand auger through the gully walls to a depth of about 1m. In areas where there was no gully, the digging was done vertically. Samples collection at the depth of 0.5-1m was to ensure that relatively fresh samples devoid of plant roots were collected. The soil samples were bagged in a cellophane bag and were properly labelled before transporting to the laboratory for analysis.
Soil geotechnical properties
Air dried soil samples that passed American Society for Testing and Materials (ASTM) sieve No. 40 (< 425µm) was used to determine the Atterberg limit of the soils following ASTM D2487 (2011) standard procedure, and the plasticity index was calculated from the results thereafter.
The soil samples which has been air dried in the laboratory was subjected to shear strength test using the direct simple shear (DSS) test method following ASTM D3080-04 (2004) standard procedures. This test was carried out under consolidated – drained (CD) condition. For each sample, the test was repeated for four different times with increase in the normal stress at every round; thus the shear stress required to cause failure was determined for each normal stress. The failure envelope was obtained by plotting the points corresponding to the shear strength (shear stress at failure) at different normal stresses and joining them afterwards with a straight line. The intercept of this straight line on the vertical axis gives the cohesion (c) of the soil sample, whereas the inclination of line to the horizontal axis gives the angle of internal friction (ϕ).
Values of the plasticity index, cohesion, and friction derived from the geotechnical tests were used in creating a geospatial map which distributed these values in space within the study area. This was done by plotting the values of the geotechnical parameter against the geographical reference point where they are sampled. The distribution of the geotechnical data across the study area followed kriging’s interpolation method as described in (Séguret and Huchon, 1990; Wackernagel, 2003). This resulted in three maps each of cohesion, friction, and plasticity of the soil samples.
Geomorphologic and geoenvironmental properties
The slope information was extracted from the DEM of the study area using the extract tool function in ArcGIS 10.0. The extracted information was used to produce a slope map of the study area. Similarly, the population density (POD) data was extracted from the population density raster map of the study area with the aid of ArcGIS 10.0 clip tool function.
Geospatial Analysis
The geotechnical, geomorphologic, and geo-environmental data was analysed following this stepwise procedure in order to generate the erosion susceptibility map.
1. Data rasterization and projection
All the polygon maps which include soil geotechnical properties maps (cohesion, friction and plasticity), gully event map, slope map, and population density map were converted to raster. During this process, the maps were set at a cell size of 30 by 30 m and were projected geographically using the geographic coordinate settings of the study area.
2. Map classification
The raster map was classified into several classes depending on the distribution of the values in the erodibility factors that was considered. The geotechnical parameters were classified into 9 classes each, while slope and population density were classified into five classes each (Tables 1–5). The classification of the geotechnical parameters, slope, and population density maps were done following the Jenks natural break method.
3. Reclassification (Class weightage assignment)
Weight contribution of the individual erodibility factor class (EFC) was assigned following the frequency ratio (FR) method. This method has been explained by (Lee and Sambath, 2006; Sharma and Mahajan, 2018). It relates the frequency of gully erosion occurrence to the erodibility factor under consideration. It can be mathematically represented as the area covered by gully erosion event (GEA) divided by the area covered by a particular erodibility factor class (EFCA) as shown in Eq. 1,
FR = \(\frac{GEA}{EFCA}\) (Eq. 1)
GEA was determined from the gully event raster map (Fig. 3a), while the EFCA was determined from the geotechnical and geoenvironmental parameter raster maps (Figs. 3b-f) using the Tabulate Area Function (TAF) in ArcGIS version 10.0.
The relative frequency (RF) of the individual erodibility factor class, which is the percentage contribution of the FR was calculated from Eq. 2,
\(RF= \frac{FR}{\sum FR} \times 100\) (Eq. 2)
The value of RF was used to assign weightage values to the Individual factor classes (Tables 2–5).
4. Ranking (Factor weightage assignment)
The erodibility factors which were considered for this analysis were ranked based on their percentage contribution (PFC) to the frequency of gully erosion. This was calculated from Eq. 3;
\(PFC= \frac{FW}{\sum FW} \times 100\) (Eq. 3)
Where FW is factor weightage, which is the summation of the frequency ratios of each class in an erodibility factor under consideration. ΣFW is the sum of all the factor weightage for the erodibility factors that was considered. Ranking was done using the weighted overlay function of the ArcGIS version 10.0.
Two scenario of the erosion susceptibility map was created. The first scenario considered only the soil geotechnical parameters (cohesion, friction, plasticity), and slope; while the second scenario considered all the factors in the first scenario in addition to population density. The resultant susceptibility map was classified into five classes very low, low, moderate, high, and severe.
Correlation Analysis
Pearson’s correlation was used to evaluate the relationships between the erodibility factors and the frequency of gully occurrence. This was achieved by plotting the mean of the erodibility factor class against their corresponding frequency ratios (Tables 1–5), using the linear model function in Genstat – a statistics analytical tool.
Model validation
Sixty (60) gully event points which were not used in building the susceptibility maps was used for the validation of the model. This was done by creating a post map with the gully event location points, which was thereafter superimposed on the generated susceptibility map. The number of gully event points that coincides with each susceptibility class was delineated using the identify object function in ArcGIS, and the percentages were thus calculated. The percentage of gully event for each susceptibility class was used in validating the gully erosion susceptibility model.