Spatial modeling of erosion hotspots using GIS-RUSLE interface in Omo-Gibe river basin, Southern Ethiopia: implication for soil and water conservation planning

Soil degradation due to soil erosion is one of the major environmental threats in developing countries. In resource limited conditions, computing the spatial distribution of soil erosion risk has become an essential and practical mechanism to implement soil conservation measures. This study aimed to assess the spatial distribution of soil loss in Omo-Gibe river basin using the integration of computer-based RUSLE and ArcGIS 10.7.1 to identify areas that require erosion prevention priority. Once raster layer of the input parameters was created, overlay analysis was carried to assess the spatial distribution of soil loss. The estimated annual soil loss varies from 0–279 t ha−1 yr−1 with a mean annual soil loss of 69 t ha−1 yr−1. The empirical analysis also confirmed that the basin losses a total of about 89.6 Mt of soil annually. Out of the total area; 7% was in very sever class, 4.8% was found in the sever and 8.7% was categorized in very high range. The remaining area were ranging from low to high erosion risk class. The influence of the combined LS factor for soil loss is significant. It was observed that small area of the Omo-Gibe basin contributed for the significant amount of soil loss. The finding of this study is in a good agreement with previous studies. Compared to the country permissible soil loss rate, 26% of the entire basin significantly exceeds the country threshold value (TSL = 18 t ha−1 yr−1). As a result, precedence and immediate attention should be given to those erosion prone areas. The study output could deliver watershed management experts and policy makers for better management implementation and resource allocation based on the local context.


Background
The world-wide adverse influence of soil erosion has been considered as the most critical issues resulting in both on-site and offsite effects (Zhou et al. 2014;Aiello et al. 2015;Zhou and Wu, 2008). In developing countries like Ethiopia, soil erosion by water and the resulting land degradation (Alexandridis et al. 2015;Panditharathne et al. 2019) are an alarming problem that led to a significant economic crisis and negative environmental threat (Shiferaw et al. 2009;Fazzini et al. 2015). Several scholars confirmed that anthropogenic factors (Balabathina et al. 2019) mainly rapid population growth, cultivation on steep slopes and rapid land use changes due to intensive agricultural practices aggravate soil erosion in Ethiopia (Abate, 2011;Tesfaye, 2015). Soil erosion and nutrient depletion is one of the major bottlenecks affecting the sustainability of agricultural production (Adugna et al. 2013;Wolka et al. 2015;Beshir and Awdenegest, 2015).
Modeling the spatial distribution of erosion risk can be used as predictive tools (Popp et al. 2000) and has become essential for land managers and policy makers to Girma and Gebre Environ Syst Res (2020) 9:19 implement sustainable intervention measures (Fernandez et al. 2003;Lu et al. 2004). Quantitative estimation of soil loss to identify erosion prone areas provides useful information to implement suitable intervention measures (Wischmeier and Smith 1978;Shi et al. 2004;Haregeweyn et al. 2013). In this regard, the integration of Revised Universal Soil Loss Equation (RUSLE) with geographic information system (GIS) provides a rather simple and yet comprehensive erosion quantification framework (Gelagay and Minale, 2016). In 2018, Gashaw (2018) and his colleagues has been used the RUSLE empirical erosion model as it is clear and relatively requires simple computational input data.
The advancement in GIS and remote sensing technologies assists to attain the RUSLE input parameters at limited costs and with reasonable accuracies (Phinzi and Ngetar 2019). The RUSLE model is the most widely-used model for erosion assessment and conservation planning (Meshesha et al. 2012 andWolka et al. 2015). The model has been tested in different part of Ethiopia by modifying some of the factors and found valid (Meshesha et al. 2012;Belayneh et al. 2019). In this background, this assesses the spatial distribution of soil loss and identify areas that require prior soil conservation measures using RUSLE integrated with GIS at Omo-Gibe river basin.

Study area description
Omo-Gibe River Basin is situated in the South-West part of Ethiopia, between 4°30′ and 9°30′ N latitude, and between 35° and 38° E longitude (Fig. 1). It flows from the northern highlands through the lowland zone to discharge into Lake Turkana at the Ethiopia and Kenya in the south. It encompasses parts of two National Regional States; Oromia which occupies the north-eastern part of the basin and the rest of the basin, the study area, is situated in the Southern People's Regional States. It is drained by two major rivers; Gibe and Gojeb river. The northern part of the basin has several tributaries which the largest are the Walga and Wabe rivers. The Tuljo and Gilgel Gibe rivers drains to the Gibe (Water Works Design Supervision Enterprise (WWDSE) and South Design and Construction Supervision Enterprise (SDCSE) 2013).

Data collection and processing
In a given area, the magnitude of soil erosion rate is varied, both temporally and spatially, due to the existing Fig. 1 Boundary of the study area local condition mainly biophysical and land management variables (Wischmeier and Smith 1978;Renard et al. 1997 andMorgan 2005). The central idea is-the collection of spatial data is crucial (Lulseged et al. 2017). To this end, the following datasets were collected from different sources and processed using the conventional methods in ArcGIS 10.7.1 environment.

Soil erodibility factor (K)
Soil erodibility (K), the susceptibility of soil towards erosion, is highly dependent on the inherent properties of the soil (Wischmeier and Smith, 1978;McCool et al. 1995). K is the rate of soil loss per rainfall erosion index (R) for a standard condition of bare soil, recently tilled up-and-down slope with no conservation practice and on a slope of 5° and 22 m length (Renard et al. 1997 andMorgan, 2005). Soil map, in vector format, collected from Ethiopian Ministry of Water, Irrigation and Energy (MoWIE) was converted in to raster map using the Feature to Raster tool in ArcGIS. The K factor value was then estimated based on a formula (Eq. 2) adapted from Williams (1995) as follows in raster calculator (Gezahegn et al. 2018): where f csand is a factor that gives low soil erodibility factors for soils with high coarse-sand contents and high values for soils with little sand, f cl-si is a factor that gives low soil erodibility factors for soils with high clay to silt ratios, f orgC is a factor that reduces soil erodibility for soils with high organic carbon content, and f hisand is a factor that (1) R = −8.12 + (0.562 * P) (2) K RUSLE = f csand × f cl−si × f orgC × f hisand reduces soil erodibility for soils with extremely high sand contents. The factors are calculated as be (Eq. 3-6) (Neitsch et al. 2002): where m s is the percent sand content (0.05-2.00 mm diameter particles), m silt is the percent silt content (0.002-0.05 mm diameter particles), m c is the percent clay content (< 0.002 mm diameter particles), and orgC is the percent organic carbon content of the layer (%).

Slope length and steepness factor (LS)
In a particular are, the effect of topography on soil erosion is represented by its slope length and steepness condition. According to Wischmeier and Smith (1978) and Schmidt et al (2019), considering the two as a single topographic factor, LS, is more convenient. LS was generated using freely available 30 m*30 m resolution digital elevation model (DEM) from ASTER DEM. L factor was calculated by Eq. 7 using raster calculator (Wischmeier and Smith, 1978;Desmet and Govers, 1996). The equation considers the flow accumulation and adding a ratio of rill to interrill erosion (Schmidt et al, 2019).
where L i,j is the slope length factor for grid cell (i,j), A i,j-in is contributing area (flow accumulation) in m 2 at the inlet of grid cell with coordinates (i,j), D is the grid cell size in m (20 m in this study). X i,j = sinα i;j + cosα i,j, α i;j is the aspect direction of the grid cell with coordinate (i,j), m is a variable slope-length exponent is related to the ratio β of rill and interrill erosion ( Fig. 8) and β-value was calculated by Eq. 9 (McCool et al. 1989). (3) where θ is the slope angle in degree. According to McCool et al. (1989), soil loss increases more rapidly with slope steepness than it does with slope length and the slope steepness (S) factor was computed by Eqs. 10 and 11 using raster calculator.

Cover and management factor (C)
The C-factor represents conditions that can most easily be managed to reduce erosion (McCool et al. 1995). For this study, cloud free Landsat 8 Imagery from https :// glovi s.usgs.gov of 2018 was attained to classify the LULC on the basis of spectral signatures and terrain characteristics. Prior to classification, image preprocessing was implemented and classification was processed using supervised image classification technique according to the desired decision rule of maximum likelihood algorithm (Abiyot et al. 2018;Suji et al. 2015) using ERDAS IMAGINE 2018. Classification accuracy was assessed by using overall accuracy (OA) and Kappa coefficient (K) based on error matrix. The error or confusion matrix was calculated by comparing the classification results and ground truth data (Hassan et al. 2016). The accuracy was assessed by using 521 reference points. The OA is the total percentage of pixels correctly classified. It was calculated as ratio (Eq. 12) between the number of correctly classified pixels (a) and the total number of pixels (b) used for accuracy assessment (Shao et al. 2016).
Kappa test was performed to measure the extent of classification accuracy as it considers all of the elements of the error matrix (Hassan et al. 2016;Bouaziz et al. 2017), was computed as a statistical measure of interrater reliability, was computed as given by Eq. 13 (Shao et al. 2016): where r is the number of rows in the matrix, x ii is the number of observations in row i and column i, x i+ and x +i are the marginal totals of the row i and column i, respectively, and N is the total number of observations. Hence, the derived thematic LULC raster map was imported in ArcMap to assign the corresponding (9) β = sinθ 0.0896 / 3(sinθ) 0.8 + 0.56 (10) S = 10.8 sin θ + 0.03 for slope in percent < 9% (11) S = 16.8 sin θ + 0.5 for slope in percent ≥ 9% C-factor values adopted from different literatures (Table  4) (Erdogan et al. 2006;Hurni, 1985;Bewket and Teferi, 2009;Tadesse and Abebe, 2014;Abate, 2011;Gelagay and Minale, 2016;Ganasri and Ramesh, 2015;Eweg et al. 1998;Wischmeier and Smith, 1978) and raster map of C-factor was generated.

Support and conservation practices factor (P)
The erosion control practice factor (P) reflects the effects of measures to reduce the amount and rate of water runoff and thus soil erosion. According to McCool et al. (1995), the P-factor mainly represents how surface conditions affect flow paths and flow hydraulics. For this study, P-value was adopted from Wischmeier and Smith (1978) where the slope of the area was correlated with the land use types as management activities are highly dependent on slope of the area. It was implemented in Ethiopia by different scholars (Yesuph and Dagnew, 2019;Gelagay and Minale, 2016;Beshir and Awdenegest, 2015;Gizachew, 2015;Gashaw et al. 2018;Tiruneh and Ayalew 2015). Procedurally, the identified LULC and slope in percent of the study area has been reclassified based on the desired range (Table1). The LULC was grouped in to cultivated land and other land uses. The category under "other land uses" were given the P-value of 1 regardless of their slope class whereas the cultivated land was classified into six slope classes. Spatial analyst Boolean-And operation was implemented on the reclassified cell values of the two input raster (LULC and slope) in ArcMap. Then the outputted LULC and slope combinate raster was used to assign the corresponding p values with the help of ArcMap editor tool and P factor value raster map was produced.

Quantitative estimation of soil loss
Annual soil loss rate was estimated by Eq. 14 superimposing and multiplying the respective preceding RUSLE factor values interactively using Raster Calculator in Arc-GIS 10.7.1 (McCool et al. 1995). where A is the amount of soil erosion (t ha −1 yr −1 ), R is the rainfall erosivity factor (MJ mm ha −1 h −1 yr −1 ), K is the soil erodibility factor (t hr MJ −1 mm −1 ), LS is a slope length and steepness factor (dimensionless), C is a cover factor that accounts for the land use and cover (LULC) class (dimensionless), and P is a conservation/management support practice factor (dimensionless).

Results and discussion
This study was intended to assess the spatial distribution of soil erosion using GIS-RUSLE interface model (Fig. 2). In this regard, each input parameters were derived from different data sources and the results were discussed as follows:

Rainfall erosivity (R-value)
According to the IDW result (Table 2), the average annual rainfall of the study area ranges Fig. 2 The GIS-RUSLE interface framework implemented in this study between 546 and 1243 mm. The R-factor analyzed using Eq. 2 indicated that, the value ranges from 299 to 690 MJ mm ha −1 h −1 yr −1 with a mean of 563 MJ mm ha −1 h −1 yr −1 . The spatial distribution (Fig. 3) confirms that there is a variation of rainfall erosivity value in the entire basin. The upper and central part of the basin is relatively dominated by high R value. A study conducted by Gerawork (2014) on Gibe-III dam catchment (central part of Omo-Gibe river basin) also revealed the rainfall variability and its potential impact on erosion. At the lower part and near to the basin outlet, R value accounts for low range (Turmi and Dimeka stations). These trends were mainly reliant on the spatial distribution of the average annual rainfall. Hence, spatio-temporal variations in soil loss are highly connected with variations in the R factor (Shina et al. 2019) and this revealed the importance of rainfall erosivity to assess erosion hotspot areas.

Soil erodibility (K-value)
Physico-chemical properties of the soil affect its resistance to detachment and transportation. The K value in Omo-Gibe river basin ranged from 0 to 0.22 t hr MJ −1 mm −1 (Fig. 4). Most of the lower part of the basin was dominated by loamy soil and characterized with high K value ranging from 0.17 to 0.22 t hr MJ −1 mm −1 ; hence these soils are highly affected by erosion. On the other hand, the upper and central part is more of clay dominant. Clay loam and sandy nature of soil was also found. These soils have a moderate K value ranging from 0.14 to 0.16 t hr MJ −1 mm −1 while low K value (< 0.14 t hr MJ −1 mm −1 ) was observed in the lower border and outlet of the basin that revealed the area more resists the impact of rainfall kinetic energy. In 2019, Melku (2019) and his colleagues were confirmed the central part (Geshy subcatchment) is dominated by clays with slow infiltration rate. Clay and sandy dominated soils have low K value because of resistant to detachment and high infiltration rates respectively. Silty loam soils have moderate to high K values, as the soil particles are moderately to easily  . 3 Mean annual rainfall and R-factor value detachable, infiltration is moderate to low (Gizaw and Degifie, 2018). Wall et al (1987), concluded that soils with high infiltration rates, higher levels of OM and improved structure have a greater resistance to erosion. Nonetheless, soil erosion is a function of many factors. Gizaw and Degifie (2018) calculated the mean K value for central part of the Omo-Gibe basin (Gilgel Gibe-I catchment) as 0.358 t hr MJ −1 mm −1 .

Topographic (LS-value)
The slope angle and aspect direction of the study basin is ranging 0-64 0 and − 1 to 359.8 0 respectively. The β-value derived from Eq. 9 is found between 0 and 3.02. The m-value is ranging from 0-0.75. In the central part of the basin, the gradient was characterized with steep slope ranging up to 64 0 (Fig. 5). In the southern part, the area has been dominated by lower elevation and flat slope (< 3 0 ). The combined LS factor varies from 0 in the lower part of the basin to 46.58 in the steepest central and upper part of the basin (Fig. 5). As noted by Shiferaw et al. (2016), steep slopes with dissected hills characterize the highland part while the low lands are characterized by relatively gentle and undulating slopes. Areas with higher LS values are generally located in the steep slope terrain.

Cover and management (C-value)
Based on the 2018 Landsat image analysis, seven different LULC (Table 3) classes were identified where 54.4% of the study basin was dominated by cultivated land. Forest cover and grassland was accounted for 16.2% and 14.8% respectively. The smallest class was settlement/ built-up area. A study conducted in Gigel Gibe-I catchment (central part of Omo-Gibe) also found the same trends of LULC (Gizaw and Degifie, 2018). The OA (86.94%) and K statistical test (0.8) conducted using 521 sampled reference observations has been found within the range of the acceptable limits (Anderson et al. 1976) (Table 3). For each LULC types, the corresponding C-values were assigned in Table 4 and Fig. 6 showed the spatial distribution of C values. Most of the southern part was covered by shrubs and forest cover dominated the north western part of the study basin; hence this part of the basin has assigned lower C value. Though cultivated land Fig. 4 Major soil type and the computed K factor value was found in the entire area, mostly dominants the north and north eastern part. So that this part including bare ground was accounted for the highest C value. In such area it is unlikely to reduce the direct impact of rainfall as compared to perennial and forest cover (Gelagay and Minale, 2016).

Erosion management (P-value)
The entire LULC was grouped in to cultivated land where further classified into six slope class and given P-values and land uses classified under forest, water body, grass, shrub, bare and built-up were grouped as other land uses given the P-value of 1. The P value obtained from the correlation between the slope in percent and LULC of the study area was ranged from 0.1 to 1. As can be seen in Fig. 7, P value of 1 is observed in the southern and some of the north western central parts of the basin. On the other hand, lower P-value (0.1) is dominant on the upper (north east) part of the basin. The higher the P value, the more the area is dominated by grass cover and shrub land where erosion management practice were not implemented (Gelagay and Minale, 2016).

Soil loss estimation (A)
With the help of ArcGIS 10.7.1, the RUSLE input parameters (R, K, LS, C and P) were created as a raster layer in a grid format. These thematic layers were then processed in raster calculator tool to generate map that shows the annual soil loss of the study area. The GIS-RUSLE based estimation showed the annual soil loss ranging from 0 in flat terrain to 279 t ha −1 yr −1 in the steep slope central area and extended to the upper part of the basin (Fig. 8).
The mean annual soil loss is. 69 t ha −1 yr −1 and the entire basin losses a total of about 89.6 Mt of soil annually. Compared with the tolerable soil loss limit (TSL), 26% (1,494,066.6 ha) of the entire basin area is by far higher than the maximum limit (18 t ha −1 yr −1 ) determined by Hurni (1985). Gizaw and Degifie (2018), reported mean annual soil loss of 62.98 t ha −1 yr −1 in the central part of Omo-Gibe (Gilgel Gibe-I catchment). For the year 2013, 60.9 t ha −1 yr −1 mean soil loss was recorded in Jimma Zone ranged from 1.6 to 232.4 t ha −1 yr −1 (Beshir and  Awdenegest, 2015). In the Ethiopian highlands, soil erosion ranges from 16-300 t ha −1 yr −1 (Hurni, 1988). According to Tesfaye and Bogale (2019), high amount of soil loss rate was recorded in the upper catchment of Omo-Gibe river basin (Gilgel Gibe -III watershed) due to deforestation, spares land cover, shallow soil depth and high rainfall intensity. Most of the central parts of the Omo-Gibe river basin is characterized by steeply sloping terrain; hence higher soil loss was estimated in this area. Out of the total soil loss, 44% (highest amount) was contributed from Weyibe Zigna Zege subbasin where 35% of its slope exceed 15 0 and the lowest amount (2.9%) was drown from Hamerkake Omo subbasin where more than 95% of the sub-basin has a slope lower than 15 0 .
This denotes there was spatial soil loss variability and the influence of the combined LS factor for soil loss is significant in the central and upper part of Omo-Gibe river basin than the lower part. As noted by Adugna et al. (2015), the amount of soil loss in Ethiopia mainly relay on the degree of slope gradient, characteristics of rainfall intensities and type and/or intensity of land cover. Though there are different factors, several studies were also indicated the impact of LS factor on soil erosion in different part of Ethiopia. According to Gashaw et al. 2018, high soil erosion rate (237 t ha −1 yr −1 ) was recorded in the hilly terrain of the Geleda watershed, northern Ethiopia. Gezahegn et al. 2018, found a mean annual soil loss of 51.04 t ha −1 yr −1 for the year 2000 in the eastern part of the country. In 2014, Amsalu and Mengaw also reported an extremely high amount of soil loss (504.6 t ha −1 yr −1 due to steep terrain) much greater than the soil loss estimation of the Ethiopian highlands. Gebreyesus and Kirubel (2009), observed the highest soil loss from steep slopes in Medego watershed. Generally, the finding of this study is in a good agreement and realistic compared with previous studies.

Prioritization for soil conservation planning
In limited resource condition, prioritizing erosion hot spot area is important for land managers and policy makers to emphasize on effective planning and implementation of appropriate intervention measures (Gashaw et al. 2018;Lu et al. 2004;Shi et al. 2004 andHaregeweyn et al. 2013). To prioritize for intervention planning, the study area was subdivided into eight major sub-basins and further categorized in to six erosion severity classes (Table 5): about 53.29% of the area was categorized under low erosion risk which extends from 0-7 t ha −1 yr −1 that contributes 5.5% of total soil loss; 17.61% is characterized in moderate class ranging between 7-15 t ha −1 yr −1 and accounted 11.6% of soil loss; 8.56% of the area is found in the high erosion risk class (15-25 t ha −1 yr −1 sharing 10.8% of soil loss); 16.3% of the annual soil loss was derived from 8.7% of the entire area which was categorized under very high erosion risk (25-45 t ha −1 yr −1 ), while the remaining 4.82% and 7.01% area was categorized as sever (45-60 t ha −1 yr −1 accounted 16.5% of total soil loss) and very sever (> 60 t ha −1 yr −1 where 39.3% of the total soil loss was contributed from) respectively. As confirmed by several studies (Amare et al. 2014;Gashaw et al. 2018;Abate, 2011), it was observed that small area of the Omo-Gibe basin contributed for the significant amount of soil loss. Areas characterized under very sever erosion risks were given the first priority and the viceversa for soil conservation planning.

Conclusions
Soil erosion by water is considered as the major cause degradation processes. Understanding the extent and its spatial distribution is essential to make sustainable land management more effective with limited resources especially in developing countries like Ethiopia. The empirical analysis of this study indicated that the mean annual soil loss from the entire area was 69 t ha −1 yr −1 which significantly exceed country's threshold value. Though there were several factors, the influence of topographic factor was significant on soil erosion variability. The overall result explicitly revealed that more than half of the area was categorized between low to moderate soil loss range, whereas 26% was identified as erosion hotspots accounted for high to very sever risk class. It was observed that small area of the Omo-Gibe basin contributed for the significant amount of soil loss. Out of the total soil loss, 44% (highest amount) was contributed from Weyibe Zigna Zege sub-basin where 35% of its slope exceed 15 0 . As a result, urgency should be given to those erosion prone areas for watershed management planning and implementation. Further studies should also be conducted to understand the major causes of soil erosion in the study basin.

Abbreviations
DEM: Digital elevation model; GIS: Geographic information system; IDW: Inverse distance weighted; LULC: Land use land cover; RUSLE: Revised universal soil loss equation; TSL: Tolerable soil loss.