The nexus between land use land cover dynamics and soil erosion hotspot area of Girana Watershed, Awash Basin, Ethiopia.

Maintaining hilly agriculture and food security remains challenging due to the ongoing degradation of the land caused by soil erosion on Ethiopia's highlands. . Soil erosion is one of the major problems affecting land and water resources. With the increase of land-use change, erosion and soil degradation increase significantly, leading to a loss of fertile soil every year. This study was therefore designed to identify erosion hotspot areas and their spatial and temporal alteration with land use land cover (LU/LC) change in the Girana watershed to give an option to local government decisions makers towards watershed management strategies. An attempt was made to combine a set of factors such as topographic wetness index (TWI), soil type, land use (1989 and 2019), slope, rainfall, and gully locations using geographic information system (GIS) based multi-criteria decision analysis (MCDA) to achieve the stated objective. Criterion maps of each factor have been processed and the factors were weighted using analytical hierarchy process (AHP) based pair-wise comparison methods, and weights have been combined using weighted overlay analysis to obtain the final erosion hotspots areas of the two-time references (1989 and 2019). The result found that 0.01%, 8.01%, 84.06%, and 7.92% of the total area fall under highly sensitive, moderately sensitive, marginally sensitive, and currently not sensitive erosion risk zone respectively for the year 1989 and 0.06%, 17.42%, 80.88% and 1.63% of the total area fall under highly sensitive, moderately sensitive, marginally sensitive, and not sensitive erosion risk zone respectively for the year 2019. Parts of the area that are highly sensitive, and moderately sensitive to Soil erosion classes increased markedly for the last thirty years in the Girana watershed, as a result of the conversion of thousands of forest areas to cultivated land and residential area. Therefore an urgent soil conservation intervention in hotspot areas is compulsory in the Girana watershed.

Introduction

Soil erosion is one of the most critical global environmental challenges, with human activities and land-use change having a significant impact on the quality of soil, land, and water resources on which humans rely for existence (Hrissanthou et al., 2010; Pimentel, 1993, 2006). It is one of the results of landscape modification and/or change as a result of poor land management (Hurni et al., 2005).

It reduces soil productivity (on-site effect) and has negative effects on water quality, such as increased sedimentation, chemical contamination from agricultural fertilizers and pesticides, and the likelihood of flooding (off-site) (El-Swaify, 1994; Saha et al., 2019; Weifeng and Bingfang, 2008). Soil erosion rates are ten to fourteen times higher than soil formation rates worldwide (Pimentel, 2006): nearly 10 million hectares of agricultural land are lost to soil erosion each year (Amiri et al., 2019; Cerdà et al., 2017) and more than 60% of eroded soil ends up in water bodies (Pimentel, 2006).

Rainfall, land use land cover (LULC), soil properties, landform, conservation, and management methods all affect a watershed's soil erosion and sedimentation vulnerability (Chang and Bayes, 2013; Lambin and Geist, 1990; Sharma et al., 2011; Sun et al., 2014). The resistance of a terrain unit to soil erosion is determined by the sort of surface cover it has. Poor surface cover promotes soil erosion, land degradation, and the elimination of habitat and biodiversity, as well as a rapid reaction to rainfall and excessive runoff (Kiage, 2013; Ziadat and Taimeh, 2013). Cultivation of slope slopes in highland locations with significant rainfall causes soil erosion and increases soil loss by extending gullies and allowing new gullies to form in the watershed (Zeleke and Hurni, 2001).

LULC is the most important predictor of soil erosion susceptibility; it delays or accelerates soil erosion rates depending on the situation (Chen et al., 2001). Uneven patterns of LULC are the strong determining factor of hydrological responses of watersheds (Siriwardena et al., 2006). As a result, LULC management is critical in improving soil characteristics and vegetation quality, both of which help to prevent soil degradation (Kosmas et al., 2000b; Wei et al., 2007). Severe soil erosion degrades soil quality by removing topsoil and, causes flooding, disrupts agricultural activities, jeopardizes community and regional security, and can damage dams and deplete reservoirs and it may eventually clash with their original purposes (Lal, 2001; Saha et al., 2019). In Ethiopia, severe soil erosion and the resulting land degradation are seriously affecting agricultural productivity, farmers' livelihood, and the growth of the country's economy (Ananda and Herath, 2003; Bewket, 2003; Constable, 1985; Hengsdijk et al., 2005; Sonneveld, 2002; Tesema, 1997) and the problem is serious in highland areas. According to the Ethiopian highland reclamation research, about half of Ethiopia's highland region was severely eroded in the mid-1980s, with 14 million hectares seriously damaged and over 2 million hectares beyond reclamation. The Ethiopian highlands range in altitude from 1500 to 4260 masl, with an average slope of more than 25% (Zemadim et al., 2011) and relatively high rainfall (between 900 and 1500 mm) that is lost as runoff and resulting soil erosion rather than being retained as surface water and groundwater. In the area, a soil and water conservation program started in 1970. However, because to the ineffective use of soil and water conservation practices, it was only partially successful (Bewket, 2003; Taye, 2006; Tesema, 1997). In this regard, the Girana watershed is a part of these highland areas in which all the upper mentioned problems are found boldly. Watershed Management must detect the spatial distribution and magnitude of soil erosion to develop effective remedies (Worku, 2015; Shiferaw and Singh, 2011).

The analytical hierarchy process (AHP) based MCDA modeling can be used to document the quantitative correlations between soil erosion conditioning variables and their spatial distributions in a GIS setting (Arabameri et al., 2018; Svoray et al., 2012). In this circumstance, GIS and Remote Sensing (RS) become valuable tools to set AHP scale-based MCDA to detect and map erosion hotspot locations (Baigorria and Romero, 2007; Saha et al., 1992; Jain and Das, 2009).

Multi-temporal satellite images can be used to analyze erosional features and provide vital information about seasonal land-use trends (Alkharabsheh et al., 2013; Kabisch et al., 2019). Multi-temporal satellite images are effective for extracting relevant information related with seasonal land-use dynamics for LULC when considering seasonal variability. The contribution of different variables for soil erosion can be weighted and prioritized in soil erosion modeling using AHP scale-based MCDA to further propose mitigation measures.

In the Girana watershed, soil erosion and abundance of gullies are still there as a result of both anthropogenic activities and natural processes in the area. Inadequate understanding and information on soil erosion risk areas may lead policymakers and resource managers to make poor environmental decisions that are poorly justified. Furthermore, for proper watershed management, updated information on the stimulation of LULC change on erosion hotspot locations is required. Accordingly, this study was aimed (1) to identify LULC change dynamics, (2) modeling and mapping soil erosion hotspot areas of 1989 and 2019 of Girana watershed using GIS-based MCDA, and (3) to evaluate LULC change effect on the extent and magnitude of soil erosion hotspot area.

Materials and methods

Description of the study area

The Girana watershed is located in the Amhara region of the Northern Ethiopian highlands, as part of the Awash Basin, between the coordinates of 11.4870 N and 11.9060 N latitude and 39.2710 E and 40.0870 E longitude, and covers an area of 248.72 km2 (Figure 1). The Mersa River (also known as Mersa Wuha) is a tributary of the Awash River basin that flows through the Amhara and Afar regions of Ethiopia. The study area is 15.8 km long and 19 km wide, with significant altitudinal variation in the watershed, ranging from 1422 m near to the outlet to 3585 m at the watershed's Northwest ridge. The slope gradient over the entire watershed varies from zero to 73°, with an average slope of 50.2°. The average annual rainfall reported at five weather stations over the last 30 years (from 1989 to 2019) ranges from 613 mm to 1324.9 mm, with an average annual rainfall of 978.2 mm and a bimodal distribution. June to September is the main rainy season, and March to April is the short rainy season. In general, the lower altitude parts of the watershed receive less rainfall than the higher altitude parts of the watershed, showing that mountains have an impact on rainfall patterns. The yearly average temperature obtained from the five meteorological stations is 21.2 °C, with maximum and minimum temperatures of 29.6 °C and 13 °C, respectively.

Figure 1

Location map of Girana watershed.

Materials & data

The data required for the AHP scale-based MCDA method are spatial. There were four different categories of data collected (land use, DEM, soil type, and precipitation). Table 1 presents the data sources for each category.

Table 1

Data types and sources.

Data Types	Sources
Land use	USGS Websitea
DEM (30 m Resolution)	(Shuttle Radar Topography Mission (SRTM) websiteb
Soil Type	Ministry of Water Resource, Irrigation & Energy
Rainfall	National Metrological Agency, Kombolcha branch

http://earthexplorer.usgs.gov.

http://earthexplorer.usgs.gov/.

Land use land cover change analysis

Spatial data on the surface cover types enables in assessing the resistance of terrain units to erosion as a result of surface protection (Kiage, 2013; Ziadat and Taimeh, 2013). Initially, six broad LULC classes were identified in the area; bare land, residential area, grassland, cultivated land, forest area, and shrubland. Further investigations were carried out by taking 228 sample points in the six LULC types. Cloud-free Landsat satellite images with a 30-meter spatial resolution were obtained from the USGS (United States Geological Survey) data repository (http://earthexplorer.usgs.gov) for the dry seasons (January to February) of 1989 (Landsat 5) and 2019 (Landsat 8) (Table 2).

Table 2

The source of land use land cover materials & details.

No.	Attributes to description	Image type
1		Landsat TM 1989	Landsat 2019
2	Path/row	168/52	168/52
3	Acquisition date	1/29/1989	1/16/2019
4	Type	GEOTIFF	GEOTIFF
5	Location	Ethiopia	Ethiopia
6	Satellite	Landsat 5	Landsat 8
7	Sensor	TM	OLI_TIRS
8	Band number	1,2,3,4,5,6,& 7	1,2,3,4,5,6,& 7
9	Pixel size	30 m × 30 m	30 m × 30 m
10	Map projection	UTM	UTM
11	Datum	WGS84	WGS84
12	UTM Zone	37	37

The ERDAS Imagine 2014 software was used to do image pre-processing (layer stacking, contrast stretching, and subsetting) as well as image classification. A linear method for calculating contrast stretching was used (equation 1).Where; DN′ is the digital number assigned to a pixel in the new (output image), DN is the original digital number of the pixel in the input image, MIN is the minimum value of the input image digital number, and MAX is the maximum value of the input image digital number.

For Landsat image classification, a supervised classification approach was used: primarily, ground control points were collected to comprehend the attributes of different LCTs, to assist a visual interpretation of the images, to identify reference areas (areas of interest), and for accuracy assessment. Focus group discussions (FGD) and informal interviews with elders in a local community provided further information regarding the area's long-term experience with the LULC changing trend. In ERDAS Imagine 2014 software, ground control points were first transformed to a vector file, then to an area of interest (AOI). The spectral signature of each LULC type has been extracted using the AOI.

The reliability of the LULC classification and validation was assessed using a confusion matrix (error matrix). The matrix matched ground truth data from reference sites with classified image outputs from sample areas in the ERDAS Imagine software's accuracy assessment tab. Thus, the error matrix calculated the overall accuracy, producer and user accuracies, and kappa statistic (Table 3) provides information about the quality of land cover classification (Congalton and Green, 2019).

Table 3

Error matrix of field data versus landsat 8 OLI 2019 for accuracy assessment.

Automated classification result	Bare land	Residential area	Grass land	Forest area	Shrub land	Cultivated land	Row total	User's accuracy
Bare land	17	1	0	0	3	0	21	81.0
Residential area	0	14	1	0	0	0	15	93.3
Grassland	0	0	25	5	0	0	30	83.3
Forest area	0	0	0	86	1	0	87	98.9
Shrub land	0	0	4	2	34	1	41	82.9
Cultivated land	3	0	3	0	3	25	34	73.5
Column total	20	15	33	93	41	26
Producer's Accuracy	85	93.3	75.8	92.5	82.9	96.2
Total GCPts collected	228
Overall classification accuracy	88.2%
Kappa coefficient	84.5%

The kappa coefficient was indicated to show how strong the agreement using the formula given by Congalton and Green (2019). It is a discrete multivariate technique of use in accuracy assessment.Where, TS, is total referenced samples; TCS, total correctly referenced samples; CT, column total; and RT, row total.

The accuracy assessment was conducted for the 2019 image but not for the old image (1989 year) because the accuracy assessment of image classification should be based on ground and known reference pixels of that time. The overall accuracy and the kappa coefficient value of the LU classification was 88.2% and 84.5% respectively which shows a high level of agreement between the classified image and the referenced data as well as good accuracy (Congalton and Green, 2019; Sim and Wright, 2005; and Anderson, 1976).

According to the classification results (Table 3), the forest area has a user accuracy of 98.9%, while cultivated land is the least accurately classified class among the other LCTs. The lower accuracy of the latter could be related to the precision of GPS data collected in the field, and/or area covered in vegetative and growing crops was not recognized as cultivated land. Green crops found between grasslands had a spectral signature that was identical to the closest LULC class, i.e. grasslands were recognized as grasslands.

The rate of LULC change in hectare/year and percentage share of each class in the studied periods were performed using both ArcGIS 10.4 and ERDAS Imagine 2014. According to Hassen and Assen (2018), the rates of change of LULC classes (proportion) for each LULC overtimes were calculated using Eq. (6).Where, ΔA (%), is the % change in the area of land use type between the initial year (At1) and final year (At2). On the other hand, the rate of change of LULC type was calculated by Shiferaw and Singh (2011) the following formula:Where: RΔ, rate of change; Z, recent area of land use land cover type in ha; X, the previous area of land use land cover type in ha; and W, the time interval between Z and X in years.

Soil erosion influence factors

Even if the degree of influence by each decision factor varies based on an anthropogenic and morphological characteristic of the watershed, several factors are responsible for the spread of soil erosion in Ethiopian highlands (Haregeweyn et al., 2015; Ziadat and Taimeh, 2013). Six decision factors that could potentially affect soil erosion were used in AHP scale-based MCDA modeling to locate soil erosion hotspots. In this work; evaluation, comparison, and prioritization of each soil erosion influencing factors were undertaken to reduce the number of listed factors influencing soil erosion through consultation with senior experts, FGD with local farmers, informal interviews with keynote informants, and based on scientific findings (Assefa et al., 2015; Beshah, 2003; Haregeweyn et al., 2015; Mhiret et al., 2019; Nyssen et al., 2005, 2006). In addition to the previously given facts, a comprehensive survey was performed throughout the watershed to determine what causes were worsening soil erosion on the ground. Finally, the six soil erosion affecting factors of land cover, slope, soil, topographic wetness index (TWI), rainfall, and the presence of gullies in the Girana watershed were considered and decided to use for the identification of soil erosion hotspots. Based on four soil erosion sensitivity class was adopted (Table 4).

Table 4

Factors sensitivity class.

Sensitivity class and notation	Explanation
S1- highly sensitive	Factors significantly accelerate erosion
S2- moderately sensitive	Factors sensitive but has an opportunity to reduce erosion (not extreme or excessive)
S3- marginally sensitive	Factors significantly reduce erosion (accelerate only a small extent)
S4- currently not sensitive	Factors that cannot support erosion

Land use land cover (LULC)

LULC regulates soil particle detachability and transport, as well as the infiltration of soil water. Based on soil erosion suitability analysis, the LULC criterion map was created by reclassifying the LU map of 1988 and the current (2019).

Rainfall

Rainfall's effect on soil erosion is based on its erosivity and kinetic energy considerations of falling rain, and it represents a measure of the erosive power, intensity, and contribution of rain to runoff in a typical year (Morgan, 2005). An erratic nature but a high-intensity rainfall pattern has a significant impact on susceptibility to soil erosion.

Rainfall criterion raster map of Girana watershed was prepared using the most preferable method, ordinary kriging surface interpolation technique, in spatial analysis tool of ArcGIS 10.4 environment (Benhamrouche et al., 2015; Goovaerts, 2000) using precipitation concentration index (PCI) data as an input. The point rainfall data from eight stations; seven rainfall stations surrounding of Girana watershed and the Mersa station found in the watershed (Table 5).

Table 5

Summary of rainfall stations used for rainfall raster map generation.

No.	Stations	Location		Time resolution (Years)
		Latitude	Longitude
1	Hayek	11.305	39.680	30
2	Gishen	11.516	39.354	15
3	Mersa	11.664	39.666	28
4	Sirinka	11.751	39.614	26
5	Woldia	11.827	39.592	28
6	Hara	12.670	41.420	13
7	Goshmeda	11.540	39.245	30
8	Wogaltena	11.590	39.221	28
9	Marye	11.317	39.734	15

The missed data was filled using the arthmatic mean method (De Silva et al., 2007) before using it for the analysis. Long-term mean monthly rainfall calculated from daily data was used to compute PCI. The reason for articulating rainfall in PCI was mainly because it is a powerful statistical descriptor to check the rainfall variability in the watershed. PCI indicates the temporal distribution of precipitation; as the value increases; the more concentrated the precipitation which is essential for soil erosion analysis (De Luis et al., 2011; Li et al., 2010; Martin-Vide, 2004; Oliver, 1980). The PCI values were calculated as given by Oliver (1980) (equation 8).Where Pi is the rainfall amount of the ith month; and Σ = summation over the 12 months.

PCI values of less than 10 indicate the uniform monthly distribution of rainfall, values between 11 and 20 indicate high concentration, and values above 21 indicate very high concentration (Oliver, 1980).

Soil type

Based on their physical and chemical qualities, which are determinants of soil susceptibility to erosion, the type of soil plays an essential role in influencing soil erosion processes. Soil detachability, particle mobility, soil hydraulic conductivity and infiltration rate, soil water storage capacity, and water flow paths are all controlled by soil types (Assefa et al., 2015; Berhanu et al., 2013; and Lilly, 2010). Soil type map of Girana watershed was extracted and the soil criterion map has been generated based on soil type suitability classes to soil erosion from soil type shapefile of Awash basin obtained from the Ministry of Water Resource Irrigation and Energy.

Slope

The slope is one of the most important topographical features affecting soil erosion (Guerra et al., 2017; Srinivasan and Engel, 1991). The erosivity of runoff increases with the velocity of runoff water and the volume of runoff that accumulates, and it is particularly high on steep slopes and slopes with longer slope lengths. A slope criteria map based on the slope suitability class was produced using a DEM with a resolution of 30 m.

Topographic wetness index (TWI)

A TWI is a physically-based indication of the effect of local topography on runoff flow direction and accumulation, and it describes how a terrain profile controls the distribution of water and the areas subject to water accumulation. It reflects the spatial distribution of surface saturation and runoff, which is a key factor in simulating soil erosion (Easton et al., 2008). TWI was extracted from DEM using the formula developed by Beven and Kirkby (1979), in an ArcGIS 10.4 raster calculator.Where, α is 'specific' catchment (runoff contributing) area (i.e. the upslope inflowing area normalized for a measure of contour length) and β, is the slope gradient in a degree at the grid cell calculated from the DEM. The criterion map of the TWI map was reclassified based on its susceptibility to soil erosion.

Potential location of gullies

Gullies are large and deep erosion depressions or channels that accelerate soil erosion significantly by delivering sediment downstream and can be the main erosion process in some regions (Gutierrez et al., 2009; Hughes and Prosser, 2012). To anticipate the vulnerability of a particular field of the watershed to gully formation, the threshold concept was used. When the contributing area combined with the local slope surpasses a specified threshold, a gully incision is expected. The method proposed by Govers (1994), has been used to predict the potential location and spatial pattern of gullies in the study area. Upslope contributing area or flow accumulation (α) and local slope (β) were generated from the DEM using GIS. Stream power index and TWI were determined Using α and β as input.

Stream Power Index (SPI): is the rate of energy of flowing water that is exerted on the bed and bank of a channel. It is the product of the flow contributing area and slope (Equation 11) in which the higher SPI indicates a possible source of soil erosion by concentrated flow detachment risk (Assefa et al., 2015).

TWI gives an idea of the spatial distribution sources for runoff generation in the watershed. However, the possible locations of gully have been predicted depending on the TWI and SPI and the criterion map was prepared based on soil erosion suitability analysis (Mhiret et al., 2019).

AHP scale based MCDA modeling of soil erosion susptablity

AHP scale-based MCDA is the most widely used complex decision-making and highly preferable tool in various fields involving both quantitative and qualitative factors (De Steiguer et al., 2003; Forman and Gass, 2001; Mardani et al., 2015). AHP, developed first by Thomas L. Saaty in the 1970s, is a structured technique for organizing and analyzing decisions. The key rationale for using this technique include: (1) it represents an accurate approach to quantifying the weights of decision factor, (2) provides a comprehensive and rational framework for structuring a decision problem, for representing and quantifying its elements, for relating those elements to overall goals, and for evaluating alternative solutions (FEMA, 2010), (3) in the AHP, the decision problem first decomposes into a hierarchy of subproblems, in which the decision-maker can evaluate the relative importance of its various elements by pairwise comparisons and converts these evaluations to numerical values/priorities/for each of the decision factors that can be processed and compared to one another rationally and consistently over the entire range of the problem to represent the alternatives' relative ability to achieve the decision goal (Li et al., 2019; Saaty, 2008).

In applying AHP based MCDA, Modelling the objective as a hierarchy containing the decision goal, the alternatives for reaching it, and the criteria for evaluating the alternatives; establish priorities among the deciding factors of the hierarchy by making a series of judgments based on pairwise comparisons; synthesize these judgments to yield a set of overall priorities for the hierarchy; check the consistency of the judgments to evaluate the extent to which the selection of decision factors was consistent in the identification of soil erosion hotspot areas; and a final decision based on the results of this process was the general procedures that the researchers were followed (De Steiguer et al., 2003; Forman and Gass, 2001; Saaty, 2008).

Information collected from FGD and keynote informants; the researcher's detailed field survey and experience in the Girana watershed were used as input in the weighting and ranking of decision factors. Also, General insights have been collected from scientific findings in the Ethiopian highlands (Assefa et al., 2015; Mhiret et al., 2019; Setegn et al., 2009). Pairwise comparison of the AHP method (Saaty, 1977) based on the solution of the Eigenvalue problem was used to assign a weight of each of the decision factors. Each of the six decision factors (LULC, TWI, soil type, rainfall, gully locations, and slope) was matched head-to-head with each other and a comparison matrix has been prepared to express the relative importance of each decision factor (Abeyou Wale et al., 2013) (equation 11).

As stated in Table 6, a scale of importance was broken down from 1 to 9 and the comparison was done using the scales of importance 1/9, 1/8...8, and 9 (Saaty, 1977), as shown in Table 6. The highest value denotes extreme importance, whereas the lowest value denotes absolute triviality) (Saaty, 1977, 1990).Where; aij gives the relative importance of the decision factors/alternative i and j and w1, w2 … wn, are the weights. The matrix for pairwise comparison of alternatives Ai = [aij] represents the preference between individual pairs of decision factors. The weights of the decision factors were estimated and prioritized from the normalized score table using matrix analysis after normalizing the corresponding eigenvector by the cumulative eigenvector (Assefa et al., 2015; Mhiret et al., 2019).

Table 6

Rating the relative importance of decision factors based on the AHP scale.

1/9	1/7	1/5	1/3	1	3	5	7	9
extremely	very strongly	strongly	moderately	equally	moderately	strongly	very strongly	extremely
less important					more important
2, 4, 6, and 8 are intermediate values between the consecutive odd numbers.
Description
Extremely important		the evidence favoring one criterion over another is of the highest possible order of affirmation
Very strongly important		one criterion is favored very strongly over another, its dominance is demonstrated in practice
Strongly important		experience and judgment strongly favor one criterion over another
Moderately important		experience and judgment slightly favor one criterion over another
Equally important		two criteria's contribute equally to the objective

All decision factor layers from MCDA factor generation and reclassification were multiplied by appropriate weights resulting from a pairwise comparision of decision factors. The final output map shows watershed-wise composite erosion index that relates erosion intensity of the area under a relative contribution of the given criteria. The procedure has been done for both years (1988 and 2019 G.C) which enables the revealing of the effect of LULC change on the magnitude and extent of soil erosion hotspot areas in the study area.

The degree of consistency for the judgment was tested statistically in pair-wise comparisons by using consistency ratio (CR) and the consistency index (CI) (Saaty, 1977), which have been used to obtain the rating to verify in what measure the judgment supplied is consistent (Eqs. (12) and (13)).Where RI is the random index (the mean of consistency index of 500 randomly filled matrices), ʎmax is the average maximal Eigenvalue, and ‘n’ represents the matrix order. The average value of ʎmax should be consistent with the regression of the random indexes formula (Alonso and Lamata, 2006), given by Eq. (14) below.

If CR is less than 10% (the inconsistency is less than 10% of 500 randomly filled matrices) then the matrix is of an acceptable consistency. Saaty (1977), calculated the random indices based on the number of decision criteria as shown in Table 7 below.

Table 7

Value of RI for the corresponding number of decision factors/alternatives.

n	3	4	5	6	7	8	9	10
RI	0.58	0.9	1.12	1.24	1.32	1.41	1.45	1.49

The criterion map of the six major decision factors was reclassified and generated in raster format for the weighted overlay model (Figure 2) to identify soil erosion hotspot areas for two given referenced years using weighted overlay under spatial analysis tool in the ArcGIS environment.

Figure 2

A model toolset showing a Workflow of the weighting overly analysis for soil erosion hotspot area identification of (a) for 1989 and (b) for 2019 using MCDA in arc GIS.

Validation

The aim of the proposed research was to examine and map soil hotspot locations (areas with high erosion risk) of the years 1989 and 2019. Extensive field visits and physical observation based on the location of hotspot locations as well as prospective gully sights verified the 2019 year findings. Watershed management priority areas have been determined based on the analysis results, field observation, and local knowledge consultation via focus groups and interviews with key informants.

Conceptual framework of the research

The general framework of the research has given in the figure below.

Results and discussion

This section consists of two parts: the first one, the results of image analysis on the state of LULC changes over the last 30 years (1989 and 2019) and its effect on the extent and magnitude of soil erosion hotspot area, while the second section contains the result of MCDA modeling to locate the erosion hotspot areas of Girana watershed.

The effect of rainfall on soil erosion

Rainfall erosion index maps are critical for evaluating soil erosion since rainfall has a significant impact in soil erosion. Rainfall station densities and time series data play an essential role in the preparation of a reliable rainfall erosivity index map (Bewket and Conway, 2007; Panagos et al., 2015). The PCI was calculated for the nine sites in the watershed to check for rainfall heterogeneity (Table 8 and Figure 3), and the results show that rainfall is not uniform between these stations. Six of the stations have PCI values between 15 and 20, implying a high concentration, while the other three have PCI values greater than 20, indicating a very high concentration.

Table 8

Long term mean annual rainfall (Ramm) and PCI of all stations.

Stations	Hara	Hayik	Mersa	Sirinka	Woldia	Gishen	Goshmeda	Marye	Wegaltena
Ramm	852.2	1216.6	1022.1	1090.2	1071.9	812.3	613	1328.6	800.5
PCI	17.1	15.7	16.0	15.3	17.6	23.2	24.7	18.1	21.9

Figure 3

Precipitation concentration index of all stations.

The PCI raster criterion map was created in ArcGIS 10.4's spatial analysis tool using a linear semivariogram model under ordinary kriging interpolation technique, and this raster map was reclassified using natural break statistics in the reclassify tool. Figure 4a presented a rainfall index map of the Girana watershed. According to the reclassified PCI map (Figure 4b and Table 9), 12.9% of the watershed area receives erosive rainfall, whereas the other 87.1 % receives rainfall that is moderately susceptible to soil erosion. The highest PCI values imply that a few rainy days contribute a large percentage of yearly precipitation, whereas the lowest PCI value indicates that a large number of rainy days contribute a large percentage of annual precipitation.

Figure 4

(a) rainfall raster map (b) re-classified rainfall map.

Table 9

PCI range and erosion sensitivity class.

PCI (%)	Class	Area (ha)	%
15.45–16.22	S1	3226.8	12.9
16.22–17.06	S2	21742.1	87.1

The Wollo highlands, where the research area is located, have bimodal rainfall, with annual rainfall being lower than in the Amhara region's mono-modal rainfall zones (Bewket and Conway, 2007). Rainfall variability is also higher in places with low annual rainfall (Merasha, 1999, 2003).

The effect of soil type on soil erosion

To create a soil criteria map based on soil type suitability classes, a soil type map of the area was taken from the Awash basin soil type shapefile. Five major soil types were identified, and were further divided into two erosion sensitivity classes based on their available water content and drainage class (Table 10). Eutric Cambisols (40.8%) dominate the watershed, followed by Eutric Regosols (21.6%), which have moderately well drainage classes. It is Leptosols, which covers almost 14.3% of the watershed and has an imperfect drainage class. The soil types in the Girana Watershed were depicted in Figure 5(a). According to the reclassified soil map (Figure 5b), 14.3 % of the soil is highly to erosion, whereas the remaining 85.7 % is marginally sensitive.

Table 10

Soil type and description of Girana watershed (Batjes, 1997).

No.	Soil type	Area (%)	Drainage class	AWC (mm)	Texture class	sensitivity class
1	Dystric Nitisols	13.5	moderately well	150	loam	S3
2	Eutric Cambisols	40.8	moderately well	150	loam	S3
3	Eutric Regosols	21.6	moderately well	150	Sandy loam	S3
4	Leptosols	14.3	imperfect	15	clay loam	S1
5	Vertic Cambisols	9.8	moderately well	150	clay (light)	S3
Total		100.0

Figure 5

(a) soil type map (b) re-classified soil map.

The effect of slope on soil erosion

Slope is most important topographical element that influences erosion and degradation (Kosmas et al., 2000a; Ziadat et al., 2013). The most important topographical element that influences erosion and degradation is the slope (Kosmas et al., 2000a; Ziadat et al., 2013). The Girana watershed is a leaf-shaped watershed in which hilly areas and high slope ranges dominate the east and west, while low slope ranges dominate the centre and downstream parts of the study area. The higher slope value represents the steeper slope, highly affected by soil erosion while the lower slopes value represents gentler (Table 11) (Nearing et al., 1991; Ziadat et al., 2013). Figure 6(a) shows the slope map of the Girana watershed, it ranges from 0 to 73°. The reclassified slope map (Figure 6b) indicated that from the entire area 7% is highly sensitive, 17.5% moderately sensitive, 23% marginally sensitive and 51.5% currently not sensitive to soil erosion.

Table 11

Slope range and erosion sensitivity class.

Slope (degree)	Erosion sensitivity group	Area (km2)	Percent
0–9	S4	128.58	51.5
9.001–20	S3	57.44	23
20.001–35	S2	43.7	17.5
35.001–73	S1	19.95	8

Figure 6

(a) slope map and (b) the reclassified slope map.

The effect of LULC on soil erosion

Spatial data on the surface cover types enables in assessing the resistance of terrain units to erosion as a result of surface protection (Ziadat et al., 2013). Based on the analysis of Landsat image with verification of field data, six main LULC classes of the Girana watershed were identified and the rate of LULC change has been presented in Table 12. The result revealed that there were again and losses of LULC in the last 30 years using 1989 and 2019 as a reference line in the Girana watershed (Table 12 and Figure 7(a) and (b)). Maps generated from ERDAS imagine 2014 (Figure 7a) were used as a criterion map after the land-use suitability classes have been done in a GIS environment.

Table 12

Summary of land use land cover change between 1989 and 2019.

Land use class	Suitability class	1989		2019		Rate of change (%)	Rate of change (ha/yr.)
		Area (ha)	%	Area (ha)	%	1989–2019
Bare land	S1	570.2	2.3	497.2	2.0	-12.8	-2.4
Residential area	S1	145.4	0.6	374.7	1.5	157.7	7.6
Grassland	S2	699.3	2.8	569.2	2.3	-18.6	-4.3
Cultivated land	S2	15147.0	60.4	18817.8	75.1	24.2	122.4
Forest area	S4	5194.8	20.7	1256.7	5.0	-75.8	-131.3
Shrub land	S3	3313.4	13.2	3555.3	14.2	7.3	8.1

Figure 7

Land use/land cover map of (a) 1989 and (b) 2019 Girana watershed.

According to the reclassified land use map (Figure 8a), in 1989, 2.9 % of the watershed was highly sensitive to soil erosion, 63.2 % was moderately sensitive, 13.2 % was marginally sensitive, and 20.7 % was not sensitive. In 2019, the reclassified land use map (Figure 8(b)) showed that 3.5 % of the watershed is highly sensitive to soil erosion, 77.4% is moderately sensitive, 14.2% is marginally sensitive, and 5% is not sensitive to soil erosion.

Figure 8

The reclassified land use map of (a) 1989 and (b) 2019 of Girana watershed.

The effect of topographic wetness index (TWI) on soil erosion

TWI is one of the most important surface parameters for modeling soil erosion hotspot areas as it shows the influence of topography on the saturation excess runoff generation process (Beven and Kirkby, 1979; Cheng et al., 2014; Hallema et al., 2016; Uuemaa et al., 2018). Figure 9(a) presented the TWI map of the Girana watershed. In the resulted TWI map, a high value of TWI indicates that the larger the contributing area and the lower slope gradient, and increased soil moisture, which has a higher correlation with soil erosion processes and vice versa (Conoscenti et al., 2013). The erosion sensitivity class for TWI is presented in Table 13. The reclassified TWI map (Figures 3 and 7b) was based on Ballerine (2017), and indicated that from the entire area 0.65% is highly sensitive; 3.84% moderately sensitive; 95.51% marginally sensitive to erosion.

Figure 9

(a) topographic wetness index map and (b) the reclassified erosion sensitivity class-map.

Table 13

TWI sensitivity class.

TWI	Erosion sensitivity group	Area (km2)	Percent
up to 11.5	S3	238.5	95.51
11.5–16.5	S2	9.58	3.84
16.5 to high	S1	1.62	0.65

The effect of gully on soil erosion

According to Zegeye et al. (2016), the rate of gully extension in the Ethiopian Highlands is high every year, and it is a key cause of increased soil erosion and land degradation. It has a negative impact on reservoir longevity, as well as the collapse and removal of failed material by runoff from gully banks. It happens when surface runoff is confined in a narrow concentrated stream and begins to erode ground surface channels. The threshold concept of two topographic parameters, i.e., SPI and TWI, was applied to predict the sensitivity of a particular field to gully formation (Lulseged and Vlek, 2005) in the Girana watershed. The potential locations of gullies were predicted based on a threshold value of > 18 for SPI and > 6.8 for TWI. From the figure, small gullies (plot-level) were not captured by the threshold it was because that gully formation follows a stream route. The Gully location map of the study area has been presented in Figure 10(a). As described in Figure 10(b), gully locations were given a high sensitive class while no gully location was grouped as a marginally sensitive class (Assefa et al., 2015). Gully's location was verified with an intensive field visit whereas Google Earth was used for verification of the actuality of gullies in inaccessible areas of the watershed.

Figure 10

(a) gully location map and (b) of erosion sensitivity class.

Soil erosion hotspot area mapping using AHP based MCDA modelling

Using six primary decision criteria, AHP-based MCDA modeling was used to identify soil erosion hotspot locations in the Girana watershed. The Pairwise comparison matrix was constructed by comparing factors one to one using a Pairwise comparison scale that ranges from 1 to 9. In the Girana watershed (Table 15), LULC change (on the left) is much more important than soil type (on the top), then a value of 9 was assigned at their intersection of LULC change row and soil column. Conversely, soil (on the left) of Table 14, is much less important than LULC change (on the top) therefore the reciprocal was assigned (i.e., 1/9). The highest value assigned for LCT could be because farmers dynamically changed their land use every year. This caused the initiation of soil erosion if it is not accompanied by proper land management practice. Gully location and TWI were ranked as the second most important criteria. Slope, rainfall, and soil types are the third, fourth, and fifth important criteria's in the identification of soil hotspot areas of Girana watershed respectively (Table 14 and Figure 11).

Table 15

Multi-criteria based weighted overlay analysis result of soil erosion hotspot area.

scale	Erosion sensitivity class	1989		2019		Rate of Change (%)	Rate of change (ha/year)
		Area (ha)	%	Area (ha)	%
1	S4	1969.9	7.92	406.4	1.63	-79.37	-52.1
2	S3	20907.1	84.06	20117.5	80.88	-3.78	-26.3
3	S2	1991.2	8.01	4333.8	17.42	117.65	78.1
4	S1	3.6	0.01	14.0	0.06	285.71	0.3
		24871.8	100	24871.8	100

Table 14

Pairwise comparison matrix.

	LU/LC	Gully location	TWI	slope	rainfall	soil type	Criteria Weight
LU/LC	1	5	5	7	7	9	0.499
Gully location	0.20	1	1	3	4	5	0.159
TWI	0.20	1.00	1	3	4	5	0.159
slope	0.14	0.33	0.33	1	3	4	0.088
rainfall	0.14	0.25	0.25	0.33	1	5	0.065
soil type	0.11	0.20	0.20	0.25	0.20	1	0.029

Figure 11

The weight of each decision criteria in soil erosion sensitivity.

The AHP analysis i.e. the consistency ratio, the consistency index, and weighting were calculated using IDRISI software. AHP based MCDA modeling is based on the solution of the Eigenvalue problem, i.e., the weights of factors were computed by normalizing the Eigenvector by its cumulative and multiplied by 100%.

The consistency ratio (9.1%) and ʎmax (6.57) of the comparison matrix computed results revealed the reliabilities of weights that were deemed consistent. To assign values for each sensitivity class (S1, S2, S3, S4), the percent value of each factor was divided into four using natural break statistics. Each map was then summed up using the weighted overlay method of spatial analysis tool in ArcGIS 10.4. Finally, the total raster value was re-classified equally into four sensitivity regions: S1, S2, S3, and S4 (Figure 12).

Figure 12

Soil erosion hotspot area map of (a) 1989 and (b) 2019.

For the last 30 years, parts of the Girana watershed areas which are marginally sensitive to soil erosion (S3) were noted as the dominant erosion form; diminished by 26.3 ha per year (ha/yr) on average. It was S4 (watershed areas that are not currently sensitive to soil erosion) which showed the highest rate of diminish for the given time references by an average of 52.1 ha/yr. The S2 (parts of the watershed areas which are moderately sensitive for soil erosion) was noted as the second most dominant soil erosion form for the last 30 years in the Girana watershed; showed the highest rate of expansion on an average of 78.1 ha/yr (Table 15). S1 (Areas that are highly sensitive to soil erosion) have also shown gradual increment for the last thirty years; 0.3 ha/yr on average.

In addition to steep sloppy areas, the distributions of watershed areas that are relatively sensitive to soil erosion were high in mild and flat slope areas. This could be owing to high TWI and gully densities in the middle and downstream parts of flatter areas, which have been identified as the second most influential factor for soil erosion hotspot areas in the Girana watershed. This was in line with other studies (Bayabil et al., 2010; Beven and Kirkby, 1979; Buchanan et al., 2014; Chapi et al., 2015; Guzman et al., 2013; Soulis et al., 2009; Panjabi, 2015). In addition to the effect of slope, the existence of bare lands and cultivation of steep slope areas play a vital role in the existence of the second soil erosion form (moderately sensitive for soil erosion) dominantly in the upper part of the Girana watershed (Mhiret et al., 2019; Kakembo and Rowntree, 2003).

Dynamics of land use/land cover change and soil erosion hotspot area

The susceptibility of LULC types varies to soil erosion based on the change from one type to another. Considering the degree of susceptibility to soil erosion of each land-use type, its area declining and increasing has a direct effect on the rate of expansion of soil erosion in the watershed (Bahadur, 2009; Zare et al., 2017).

In 1989 cultivation land followed by forest areas, were the top two dominant LU types in the Girana watershed whereas cultivation land, followed by shrublands was the top two dominant LU types in 2019 of the Girana watershed. As shown by Figure 13, Maximum loss was observed in forest areas between the two reference years, while a maximum gain/expansion was observed in residential areas from 1989 to 2016.

Figure 13

Land use land cover of Girana watershed for 1989 and 2019.

A high rate of LULC change was observed in the Girana watershed for the last thirty years (Table 12 and Figure 14) and the change consists of both the loss and the gain in all types of LULCs. The highest change was observed in residential areas, which increased by 157.7%, and forest areas, which decreased by 75.8%, in 30 years period (1986–2016).

Figure 14

The net gain or loss of (a) individual land use/cover categories and (b) individual soil erosion sensitivity class from 1989 to 2019.

As described in Figure 11, it is LULC that plays a great role in the expansion of soil erosion in the case of the Girana watershed (covers a weight of about 49.9%). The residential area and cultivated land were the two most expanded land use types at a high rate (157.7% and 24.2% respectively) from 1989 to 2019 (Figure 14a). Together with other decision factors, the two LU types described above have played for the expansion of S1 & S2 soil erosion class in the study area by 285.71% and 117.65% respectively (Figure 14b). In contrast, forest areas have shown a high rate of reduction (75.8%) and its contribution was high for the reduction of S4 soil erosion class by 79.37% in the Girana watershed for the last 30 years.

Conclusion

The result of the AHP scale-based MCDA model analysis has shown that a historical fluctuation of Soil erosion hotspot areas is attributable to the considerable LU/LC changes in the Girana watershed. Multi-temporal remote sensing data (1989–2019) analysis shows the interrelationship between human demand, the dynamics of LC in the watershed ecosystems, and the areal extent of soil erosion that occurs in the study area. The quantitative evidence of LULC change shows that shrubland, cultivated land, and residential areas showed an increment in areal coverage between 1989 and 2019, respectively. High population growth, urbanization, and a collective demand for food and related products may be the causes for the increment of cultivated land and residential areas. On the other hand, forestland, grassland, and bare lands showed a diminution in their areal extent. This was due to the transformation of forestland and grassland into other land use/land cover types. As revealed from the field survey and confirmed by GIS and RS analysis of satellite images the LU changed dominantly from forestland to cultivated land and residential area due to the driving force of population pressure and urbanization. This resulted in the reduction of the regulatory capacity of the land, which in turn, aggravates soil erosion-sensitive areas.

The soil erosion sensitivity class of S1 and S2 areas of the Girana watershed increased by 285.71% and 117.65% respectively, in the last 30 years. In contrast, S3 and S4 classes decreased by 3.78 % and 79.37 % respectively, in the given time resolution. This shows that soil erosion sensitivity is directly proportional to LULC dynamicity.

The unique feature of this investigation can be best described as identifying the hotspot areas of soil erosion by incorporating rainfall as one of the determinant factors for soil erosion in which most studies done in Ethiopian highlands have been ignored. Moreover, the spatial distribution of soil erosion hotspots was observed in the Girana watershed both at steep slope areas and in flatter areas of the watershed. Along with soil erosion decision factors considered for the analysis, for the existence of soil erosion source areas in flatter parts of the watershed, the contribution of saturated areas may be high. Therefore in the study area, future works shall focus on hydrological monitoring at both sub-watershed scale and at the main outlet for the analysis of such runoff contributing areas which are the main source of sediment load and gully expansion in the watershed level. Further, this concept may be extended for the estimation of soil loss in different time resolutions to evaluate the effectiveness of afforestation programs and soil and water conservation practices done at a watershed level.

The finding of this study can be used as a basis for decision making; can provide useful information for designing land use planning to regulate the effect of LC change and other soil erosion decision factors, and for sustainable development of the Girana watershed. A hotspot area of soil erosion is increased temporally in the Girana watershed, even though soil and water conservation measures are under implementation for the last three decades. Therefore, an integrated watershed management policy of the Girana watershed shall reconsider the hotspot areas of soil erosion. Finally, it could be concluded that soil erosion hotspot area maps can be effectively used to formulate appropriate management strategies and for the effective implementation of soil erosion conservation measures.

Declarations

Author contribution statement

Belachew Beyene: Conceived and designed the experiments; Performed the experiments; Analyzed and interpreted the data; Contributed reagents, materials, analysis tools or data; Wrote the paper.

Funding statement

This work was supported by Woldia University.

Data availability statement

Data will be made available on request.

Declaration of interests statement

The authors declare no conflict of interest.

Additional information

No additional information is available for this paper.