Spatial Modeling of Soil Erosion Risk and Its Implication for Conservation Planning: The Case of the Gobele Watershed, East Hararghe Zone, Ethiopia

: Soil erosion by water has accelerated over recent decades due to non-sustainable land use practices resulting in substantial land degradation processes. Spatially explicit information on soil erosion is critical for the development and implementation of appropriate Soil and Water Conservation (SWC) measures.The objectives of this study were to estimate the magnitude of soil loss rate, assess the change of erosion risk, and elucidate their implication for SWC planning in the Gobele Watershed, East Hararghe Zone, Ethiopia. Applying remote sensing data, the study ﬁrst derived the Revised Universal Soil Loss Equation (RUSLE) model parameters in an ArcGIS environment and estimated the soil loss rates. The estimated total soil loss in the watershed was 1,390,130.48 tons in 2000 and 1,022,445.09 tons in 2016 with a mean erosion rate of 51.04 t ha − 1 y − 1 and 34.26 t ha − 1 y − 1 , respectively. The study area was divided into eight erosion risk classes ranging from very low to extremely high. We established a change detection matrix of the soil erosion risk classes between 2000 and 2016. The change analysis results have revealed that about 70.80% of the soil erosion risk areas remained unchanged, 19.67% increased in total area, and 9.53% decreased, showing an overall worsening of the situation. We identiﬁed and mapped areas with a higher and increasing erosion risk as SWC priority areas using a Multi-criteria Decision Rules (MCDR) method. The top three priority levels marked for the emergency SWC measures account for about 0.04%, 0.49%, and 0.83%, respectively. These priority levels are situated along the steep slope areas in the north, northwest, south, and southeast of the Gobele Watershed. It is, thus, very critical to undertake proper intervention measures in upslope areas based on the priority levels to establish sustainable watershed management in the study area.


Introduction
Soil erosion is a naturally occurring environmental process by which soil materials are displaced, transported, and deposited in downstream areas by wind, water, or gravitational forces [1][2][3]. In the context of water-caused soil erosion, removal of soil particles is the result of raindrops, while surface runoff carried out the transportation process [4]. Though soil erosions are the result of the interplay between soil erodibility and rainfall erosivity factors, inappropriate human practices such as cultivation in upslope areas, deforestation, an extension of urban areas and roads, and uncontrolled and overgrazing aggravate the problem [5][6][7][8]. In connection to this, Knapen et al. [9], reported that Israel [70] reported the mean annual soil loss rate of 58.30 t ha −1 y −1 and recommended the implementation of conservation measures to reduce the on-site and off-site effects of soil erosion in the Dire Dam Watershed. In a related study conducted on the recently dried Lake Haramaya catchment in the East Hararghe Zone, Senti et al. [71] suggested the importance of an integrated physical soil erosion control and conservation measures to tackle the on-site and off-site effects of soil erosion. However, none of these addressed the spatial changes among different soil erosion risk grades over time. Therefore, this study was designed to (i) estimate the magnitude of soil loss rates in 2000 and 2016; (ii) assess the spatial changes among soil erosion risk classes between 2000 and 2016; and (iii) identify priority areas for SWC in the Gobele Watershed, East Hararghe Zone, Ethiopia.

Description of the Study Area
The Gobele Watershed is located in the east Hararghe Zone, Oromia Regional State, Ethiopia. The astronomic location of the watershed extends from 8°50'10" N to 9°20'30" N latitude and from 41°41'10" E to 42°11'30" E longitude, with elevation ranging between 974 and 3,264 meters above mean sea level ( Figure 1). The Gobele Watershed, with a surface area of 237,786.44 hectares, is one hydrological watershed within the Wabi Shebelle Basin. Topographically, about 77.51%, 21.58%, and 0.91% of the total study area have a slope gradient ranging from 0% to 10%, 10% to 30%, and 30% to 100%, respectively. The mean annual rainfall of the watershed is 820.01 mm, with August (152.31 mm) and April (1.16 mm) being the wettest and the driest months, respectively [72].

Data Sources
In this study, several geospatial datasets were collected and processed in a raster format to suit the RUSLE model for estimating the soil loss. We obtained the mean annual rainfall data for a period of sixteen years (1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015) with twelve rain gauge stations from the National Meteorological  (Table S1) [72]. We used two cloud-free satellite images of Landsat-7 ETM+ (Enhanced Thematic Mapper Plus) image (path/row 168/054) captured on 26 March 2000, and Landsat-8 OLI (Operational Land Imager) image (path/row 168/054) captured on 14 March 2016. The images were downloaded from the United States Geological Survey's (USGS) Earth Explorer website [73]. In addition, a digital elevation model (DEM) with a spatial resolution of 30 m was accessed from Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) Global Digital Elevation Model (GDEM) website [74]. The soil data covering the study area was accessed from the FAO website in Environmental System Research Institute (ESRI) shapefile format [75]. Furthermore, field observation was made between September 2016 and November 2016 to collect reference data of land use/land cover (LULC) classes. A total of 150 ground truth data were collected from the field stratified proportional to each LULC class to support image classification and accuracy verification. Handheld Global Positioning System (GPS) was used to collect the required data. Figure 2 shows sampled LULC types in the study area. downloaded from the United States Geological Survey's (USGS) Earth Explorer website [73]. In addition, a digital elevation model (DEM) with a spatial resolution of 30 m was accessed from Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) Global Digital Elevation Model (GDEM) website [74]. The soil data covering the study area was accessed from the FAO website in Environmental System Research Institute (ESRI) shapefile format [75]. Furthermore, field observation was made between September 2016 and November 2016 to collect reference data of land use/land cover (LULC) classes. A total of 150 ground truth data were collected from the field stratified proportional to each LULC class to support image classification and accuracy verification. Handheld Global Positioning System (GPS) was used to collect the required data. Figure 2 shows sampled LULC types in the study area.

Satellite Images Preprocessing
Before LULC classification and change detection, one should preprocess the distorted and degraded images to ensure the results of adequate quality with a more correct and faithful representation of the real ground scene. In the context of the current study, this involved removing or diminishing any undesirable image characteristics that occurred during the acquisition process [76,77]. This was applied to all raw Landsat image data on board ETM+ and OLI sensors. The Landsat images were preprocessed at each band level using the ENVI software version 5.1 (Exelis Visual Information Solutions, Inc. Boulder, CO, USA). The first step involved the re-projection of the Landsat images into World Geodetic System (WGS 84) spheroid, a spatial reference system of the Universal Transverse Mercator (UTM) with Datum Zone 37 N. The preprocessing includes a series of sequential operations. This includes an atmospheric and radiometric correction to diminish the effects of clouds and the sun elevation angle of satellite images taken during different periods and from different sensors. Geometric rectification accurately links the Landsat imagery with a ground truth data and other ancillary data sets, and masking operations [78,79]. Then, single-band images were combined

Satellite Images Preprocessing
Before LULC classification and change detection, one should preprocess the distorted and degraded images to ensure the results of adequate quality with a more correct and faithful representation of the real ground scene. In the context of the current study, this involved removing or diminishing any undesirable image characteristics that occurred during the acquisition process [76,77]. This was applied to all raw Landsat image data on board ETM+ and OLI sensors. The Landsat images were preprocessed at each band level using the ENVI software version 5.1 (Exelis Visual Information Solutions, Inc. Boulder, CO, USA). The first step involved the re-projection of the Landsat images into World Geodetic System (WGS 84) spheroid, a spatial reference system of the Universal Transverse Mercator (UTM) with Datum Zone 37 N. The preprocessing includes a series of sequential operations. This includes an atmospheric and radiometric correction to diminish the effects of clouds and the sun elevation angle of satellite images taken during different periods and from different sensors. Geometric rectification accurately links the Landsat imagery with a ground truth data and other ancillary data sets, and masking operations [78,79]. Then, single-band images were combined to get multi-band composite images [80]. Additionally, image enhancement operations, including density slicing, contrast adjustment, edge enhancement and color composite were employed to enhance the interpretability of image data [81]. Since the entire scene of ETM+ and OLI images respectively covers 170 km by 185 km and a 190 km by 180 km surface area [73], the area of interest (AOI) that covers the study area (237,786.44 ha) was extracted using vector polygon layer and the subset tool available in ERDAS IMAGINE software version 10 (ERDAS, Inc., Norcross, GA, USA).

LULC Classification
We used a computer-aided digital image classification procedure to classify satellite images to generate thematic LULC maps of the study area based on known features on the ground [82]. For this purpose, the training signatures were collected for each LULC class from satellite images aided by sampled field reference data, and knowledge of experts and people of the locality about the earlier state of land cover in the watershed. A supervised classification method based on the Maximum Likelihood Classifier (MLC) was applied to classify the Landsat image of 2000 and 2016 separately [83,84]. We applied a 3 × 3 moving window majority filtering operation for neighboring cells in classified LULC images to minimize salt-and-pepper effects. A total of seven LULC classes were identified in the study area for the years 2000 and 2016 (Table 1). This includes bare land, cultivated land, settlements, forest, grazing land, shrubland, and water body, with a share of each class in 2000 contributes 67,021.03 ha (28.19% of the total study area), 64,159.60 ha (26.98% of the total study area), 1199.16 ha (0.50% of the total study area), 7728.48 ha (3.25% of the total study area), 299.16 ha (0.13% of the total study area), 96,822.90 ha (40.72% of the total study area), and 556.11 ha (0.23% of the total study area), respectively ( Figure 3a). As shown in Table 1, each LULC class in 2016 accounted for 49,932.80 ha (21% of the total study area), 135,972.81 ha (57.18% of the total study area), 13,320.50 ha (5.60% of the total study area), 4794.84 ha (2.02% of the total study area), 1863.97 ha (0.78% of the total study area), 31,704.40 ha (13.33% of the total study area), and 197.12 ha (0.08% of the total study area), respectively ( Figure 3b). During the study period, shrubland, bare land, forest, water body have decreased from 96,822.90 ha (40.72% of the total study area), 67,021.03 ha (28.19% of the total study area), 7728.48 ha (3.25% of the total study area), and 556.11 ha (0.23% of the total study area) to 31,704.40 ha (13.33% of the total study area), 49,932.80 ha (21% of the total study area), 4794.84 ha (2.02% of the total study area), and 197.12 ha (0.08% of the total study area), respectively. On the contrary, areas covered by cultivated land, settlements, and grazing lands have increased from 64,159.60 ha (26.98% of the total study area), 1199.16 ha (0.50% of the total study area), and 299.16 ha (0.13% of the total study area) to 135,972.81 ha (57.18% of the total study area), 13,320.50 ha (5.60% of the total study area), and 1863.97 ha (0.78% of the total study area), respectively.

Accuracy Assessment
The thematic layers of the classified LULC images for 2000 and 2016 were validated using ground truth data. Out of the total sampled ground truth data collected from the field stratified proportionally to each LULC classes, we used 50 points as a reference data for image classification. The remaining 100 points were used to examine the classification accuracy of the LULC images. A lack of field reference data on the state of historical LULC of the study area necessitated using knowledge of the local elders and experts to minimize the misclassification of remote sensing images. Thus, the ground truth data collected between September 2016 and November 2016 were used as a reference for image classification and validation of classification accuracy in 2000 and 2016. We calculated the classification accuracy in terms of producers' and users' accuracy, overall accuracy, and Kappa Statistics [85]. The results revealed that the overall classification precision obtained per LULC map in 2000 and 2016 was 84.26% and 94.44%, respectively ( Table 2). The Overall Kappa Statistics (K^) calculated for each LULC images in 2000 and 2016 were 0.815% and 0.934%, respectively. The highest producers/users and K^ is found in the water bodies and cultivated land in classified images for 2000 and 2016.

Accuracy Assessment
The thematic layers of the classified LULC images for 2000 and 2016 were validated using ground truth data. Out of the total sampled ground truth data collected from the field stratified proportionally to each LULC classes, we used 50 points as a reference data for image classification. The remaining 100 points were used to examine the classification accuracy of the LULC images. A lack of field reference data on the state of historical LULC of the study area necessitated using knowledge of the local elders and experts to minimize the misclassification of remote sensing images. Thus, the ground truth data collected between September 2016 and November 2016 were used as a reference for image classification and validation of classification accuracy in 2000 and 2016. We calculated the classification accuracy in terms of producers' and users' accuracy, overall accuracy, and Kappa Statistics [85]. The results revealed that the overall classification precision obtained per LULC map in 2000 and 2016 was 84.26% and 94.44%, respectively ( Table 2). The Overall Kappa Statistics (Kˆ) calculated for each LULC images in 2000 and 2016 were 0.815% and 0.934%, respectively. The highest producers/users and Kˆis found in the water bodies and cultivated land in classified images for 2000 and 2016.

Development of the RUSLE Model Factors
The RUSLE is a practical tool for predicting the long-term average annual soil loss attributed to raindrop splash and runoff [63]. According to Jiang et al. [16], "to build the quantification model, as many as possible of the criteria that influence soil erosion should be taken into consideration." In the present study, we used the empirical prediction model of RUSEL, which is easily applicable in different scales with Geographic Information system (GIS) tools [86,87]. The application of RUSLE model require five factors: rainfall erosivity (R) factor, soil erodibility (K) factor, slope length and steepness (LS) factor, cover management (C) factor, and conservation practice (P) factor determined in an ArcGIS environment ( Figure 4) and multiplied together to estimate the amount of the soil loss rates (Equation (1)) [63]. Since the input model layers were acquired from various sources, and at varying scales, resampling procedures provided in digital analysis tools need to apply for the input datasets to be compatible with each other [86][87][88][89].
where A is average annual soil loss (t ha −1 y −1 ); R is the rainfall-runoff erosivity factor; K is a soil erodibility factor; LS is a slope length-steepness factor (dimensionless); C is a cover management factor (dimensionless), and P is a support practice factor (dimensionless). The RUSLE model was run to estimate the soil loss rate for the year 2000 and 2016 separately. Erosion risk area was categorized into eight classes from very low to extremely high, following Uddin et al. [90] for Koshi Basin. Classification of soil erosion risk area was based on the estimated rate of mean annual soil loss (Table 3). By overlaying thematic soil erosion risk map for 2000 and 2016 using a transformation matrix (GIS analysis), the study obtained information about spatial changes among the erosion risk class between 2000 and 2016. Moreover, conservation priority levels were obtained to support the spatial planning and implementation of SWC measures in the study area. Prioritization was based on soil loss assessment results and change in erosion risk classes between 2000 and 2016 with the help of MCDR method [90][91][92]. We set the highest value for areas with a high and increasing risk of soil erosion.
where A is average annual soil loss (t ha −1 y −1 ); R is the rainfall-runoff erosivity factor; K is a soil erodibility factor; LS is a slope length-steepness factor (dimensionless); C is a cover management factor (dimensionless), and P is a support practice factor (dimensionless). The RUSLE model was run to estimate the soil loss rate for the year 2000 and 2016 separately. Erosion risk area was categorized into eight classes from very low to extremely high, following Uddin et al. [90] for Koshi Basin. Classification of soil erosion risk area was based on the estimated rate of mean annual soil loss (Table 3). By overlaying thematic soil erosion risk map for 2000 and 2016 using a transformation matrix (GIS analysis), the study obtained information about spatial changes among  Rainfall Erosivity (R) Factor The R factor reflects the ability of rainfall-runoff to erode the soil particles due to the joint effect of rainfall kinetic energy, duration, and potential [93]. We calculated a mean annual rainfall data using the mean yearly rainfall data covering the period from 1999 to 2015. It follows that the rainfall spatial distribution of the study area was mapped using the spline interpolation method. Finally, spatially distributed R factor value was calculated by Equation (2) (Figure 5b and Table S1) [94,95].
where P is an annual rainfall (mm).   Soil Erodibility (K) Factor The K factor corresponds to the rate of soil loss per rainfall erosion index unit measured on a standard plot by taking into consideration the inherent soil properties [96]. It is the average prolonged effects of soil profile characteristics, soil properties (e.g., soil texture, size, thickness, organic matter, clay type, and permeability) and human activities on soil loss [22,61,63,[97][98][99]. The erodibility of a soil caused by rainfall-runoff will increase proportionally with an increase in the amount of fine sand and silt contents [97][98][99][100]. For instance, the finer and the richer the soil texture in a clay ratio is, a more resistant is soil to particle detachment and the lower the soil erodibility factor is and vice versa. Moreover, the content of the organic matter is a key factor that decides erodibility of soil layers. It contributes to the increment of particle aggregation (due to the presence of chelating agents) and water infiltration [100,101]. When the K factor value of a specific soil class gets higher, more erosion occurs as the soils are exposed to the erosive force of rainfall, splash, or surface flow [102]. According to Yahya et al. [103] the soil erodibility value mostly ranges between 0 and 1, where 0 shows the soil class's sensitivity to erosion while 1 is the high susceptibility of soil class to water erosion (Figure 5c,d). The K factor value was estimated based on a formula adapted from published literature [104][105][106][107] with the FAO harmonized digital soil map [75] as follows: where f csand is a factor that lowers the K indicator in soils with a high proportion of coarse-sand content and higher for soils with little sand; f cl−si gives low soil erodibility factors for soils with a high clay-to-silt ratio; f orgC reduces the K values in soils with a high organic carbon content while f hisand reduce the K value of soil classes with high sand contents (Table S2). The f csand , f cl−si , f orgC and f hisand was calculated using Equations (4)-(7) [101,104]: The highest soil erodibility (K) factor value was found in the Eutric Nitosols (0.42 t h MJ −1 mm −1 ), while the lowest is in Eutric Camisoles (0.33 t h MJ −1 mm −1 ) along the southeastern and northwestern part of the watershed (Figure 5d).

Slope Length and Steepness (LS) Factor
The slope length (L) and steepness (S) factor reflect the effect of terrain and topography on soil erosion [61]. The increase in slope length (L) and slope steepness (S) can result in the higher overland flow speed and higher erosion [1,63,103]. The specific effects of topography on soil erosion are estimated by the dimensionless LS factor as the product of the slope length (L) and slope steepness (S) constituents converge into a point of interest [5,108]. The LS factor is a ratio of soil loss under a given condition of a site for a slope length of 22.13 meters and slope steepness of 9%, free of vegetation and left in a seedbed condition [5,61]. The LS factor (Figure 5f where, L is a slope length factor, λ is the field slope length in meters, and m is a variable slope length exponent related the value of the slope gradient (Figure 5e): 0.5 for slopes steeper than 4.5%; 0.4 for slopes from 3% to 4.5%; 0.3 for slopes from 1% to 3%, and 0.2 on slopes less than 1%. Given the uneven distribution of slope gradients over the study area, the slope steepness (S) was sub-divided into segments using Equation (9) [110,111], as follows: where A i,j−in is the contributing area at the inlet of the grid cell (i, j) is measured in m 2 ; D is the grid cell size (meters); x i,j is sin a i, j + cos a i,j ; A i,j is the aspect direction of the grid cell (i, j); m is the slope length exponent related to the ratio of β of rill erosion and interill erosion Equations (10) and (11) [107,109]: where, θ is slope steepness angle in degrees was calculated by Equation (12)  The C factor is typically associated with the effects of cropping and management practices on soil erosion [14,108]. It is the ratio of soil loss from land with specific vegetation to the corresponding soil loss under clean tilled continuous fallow or management systems to reduce erosion [54,63]. The C factor is a dimensionless factor that ranges from 0 for a completely non-erodible condition in areas with high plant cover to 1 which corresponds to the greater magnitude of soil loss due to very extensive tillage, leaving a very smooth surface that produces much runoff and makes the soil susceptible to erosion [22,63]. In this study, the C factor was derived from the Normalized Difference Vegetation Index (NDVI) interpreted using image data from ETM+ and OLI sensors (Figure 5g,h). The equation to calculate the spectral indices of NDVI is as follows [113,114]: where NIR is the reflectance of near infrared bands, Band 4 of ETM+ and Band 5 of OLI, while Red is the reflectance of visible red bands, which are Band 3 of ETM + and Band 4 of OLI imagery. The spatial distribution of the C factor was calculated using Equation (14) [115]: Land 2018, 7, 25 13 of 25

Support Practice (P) Factor
The P factor reflects the effect of soil conservation practices that reduce the amount and rate of erosion by rainfall-runoff [116]. The P factor is mechanical practices such as the effects of contouring, strip cropping, or terracing and the resultant average annual soil loss rate due to water erosion [103]. The P factor ranges between 0 and 1, where the value closer to 0 shows good conservation practices and the values close to 1 is showing little conservation practices [47]. Several studies have tried to generate the P factor value considering the slope gradient and conservation practices [54,87,117,118]. Due to a lack of field data concerning the conservation practices that have been implemented in the study area, the P factor value corresponds to LULC classes was adopted from published literature [67,119]. The values of the P factor were set to 1.00, 0.40, 0.10, 0.3, 1.00, 1.00, and 0.00 for the bare land, cultivated land, forest, grazing land, settlement, shrubland, and water bodies, respectively. Figure 5i,j display the spatial distribution of the P factor in the study area in 2000 and 2016.

The Soil Erosion Risk in the Gobele Watershed
The soil erosion risk maps are shown in Figure   cultivated land, forest, grazing land, settlement, shrubland, and water bodies, respectively. Figure  5i,j display the spatial distribution of the P factor in the study area in 2000 and 2016.

The Soil Erosion Risk in the Gobele Watershed
The soil erosion risk maps are shown in Figure 6a in 2000 and Figure 6b in 2016 while the statistical details of the soil loss rates and associated erosion risk classes are presented in Table 4     The high, very high, and extremely high areas have also declined from 952.93 ha (0.401% of the total study area), 608.25 ha (0.256% of the total study area), and 141.93 ha (0.06% of the total study area) to 607.86 ha (0.256% of the total study area), 244.76 ha (0.103% of the total study area), and 126.67 ha (0.053% of the total study area), respectively. In the present study, the areas with erosion risk higher than the class value of very low were defined as eroded areas. Thus, the area of eroded classes decreased from 53,465.44 ha (22.49% of the total study area) in 2000 to 30,876.40 ha (12.99% of the total study area) in 2016. This suggests that the magnitude of soil erosion rates has decreased over the study period. This is partly attributed to some conservation measures taken by the local people over the recent years. Table 5 presents the percentage of each erosion risk class transformation between 2000 and 2016.

Spatial Changes in Soil Erosion Risk in the Gobele Watershed
The diagonal numbers showed in italic down the change matrix represents the proportion of each class that stayed unchanged for the total soil erosion risk classes in the study area. The above diagonal elements show the percentage decrease in the erosion risk area while an increase in erosion risk is below the diagonal. The change detection matrix reveals that about 70.80% of the total soil erosion risk classes covered in 2000 (i.e., the sum of the diagonal elements in Table 5) showed no change in 2016. Overall, the erosion risk areas increased by 19.67% of the total study area, and decreased by 9.53%, showing that the soil erosion situation is getting worsening in the study area.
The proportion of the area at very low risk of erosion was the largest unchanged class during the study period, while the area of extremely high-risk was the lowest persistent class. Out of the 77.51%, the very low class covered in 2000 about 68.69% of the total area stayed unchanged, and only about 8.82% were converted to other erosion risk classes. The highest percentage gain was also found in erosion risk in the category of very low. It accounts for 18.47% of the total area due to area converted primarily from low (11.63% of the total area), low-medium (4.38% of the total area), and medium (1.29% of the total area), respectively. Although the soil erosion at risk of very low gained area converted from the low, low-medium, and medium, it lost about 8.77% of the total area to other classes. In addition, the proportion of losses were comparatively higher in the area at low, very low, and low medium risk of soil erosion. It accounts for about 12.19%, 8.77%, and 5.11% of the total area, respectively. On the contrary, the lowest percentage losses were figured out in areas at extremely high, very high and high erosion risk by 0.06%, 0.26%, 0.40% of the study area, respectively.

Identification of Conservation Priorities
Proper identification of areas that are highly vulnerable to soil loss is a critical factor for designing and implementing appropriate SWC measures. Prioritization was done at watershed scales considering areas with a higher soil loss and increases in erosion risk. Thus, we set a higher value for areas with increasing mean annual soil loss. Accordingly, we classified the study area into eight conservation priority levels (Figure 7). As a result, about 104.78 ha (0.04% of the total study area), 1164.27 ha (0.49% of the total study area), 1963.74 ha (0.83% of the total study area) was identified and mapped as the top three conservation priority areas (Table 6). Among the topmost three priority levels, the majority of the first conservation levels have the slope gradient ranges from 30% to 50%, whereas, the second and the third levels were found within slope gradient greater than 50%. Moreover, about 2565.27 ha (79.35% of the total study area) was situated within the Kersa, Kurfa Chele and Girawa districts which are located in the north, northwest, south, and south-west of the watershed. and low medium risk of soil erosion. It accounts for about 12.19%, 8.77%, and 5.11% of the total area, respectively. On the contrary, the lowest percentage losses were figured out in areas at extremely high, very high and high erosion risk by 0.06%, 0.26%, 0.40% of the study area, respectively.

Identification of Conservation Priorities
Proper identification of areas that are highly vulnerable to soil loss is a critical factor for designing and implementing appropriate SWC measures. Prioritization was done at watershed scales considering areas with a higher soil loss and increases in erosion risk. Thus, we set a higher value for areas with increasing mean annual soil loss. Accordingly, we classified the study area into eight conservation priority levels ( Figure 7). As a result, about 104.78 ha (0.04% of the total study area), 1164.27 ha (0.49% of the total study area), 1963.74 ha (0.83% of the total study area) was identified and mapped as the top three conservation priority areas (Table 6). Among the topmost three priority levels, the majority of the first conservation levels have the slope gradient ranges from 30% to 50%, whereas, the second and the third levels were found within slope gradient greater than 50%. Moreover, about 2565.27 ha (79.35% of the total study area) was situated within the Kersa, Kurfa Chele and Girawa districts which are located in the north, northwest, south, and south-west of the watershed.    The remaining conservation priority areas within these levels account for 667.52 ha (20.65% of the total study area) and are confined within Haramaya and Fedis districts located in the eastern and southeastern part of the watershed. The three priority areas are characterized by higher soil loss rates. The areas require urgent intervention in SWC measures. The fourth, fifth and sixth conservation levels cover about 21,104.37 ha (8.88% of the total study area) and currently need of minor SWC measures. The last two conservation priority areas covered 213,449.27 ha (89.77% of the total study area) currently may not need SWC measures. The last two conservation levels representing areas with lower soil loss and erosion risk account for 213,449.27 ha (89.76% of the total area).

Discussion
In Ethiopia, land degradation is the most severe problem that affects agricultural productivity and negatively affects food security [31,35,120,121]. According to the World Bank statistical estimate, land degradation cost to annual agricultural GDP ranges from 2% to 6.75% [122]. To reduce the problem, the Ethiopian government adopted intensive SWC at the national level [35,42,43]. However, the effectiveness of the government's efforts to deal with soil degradation needs up-to-date quantitative information on the extent of soil erosion risk, and its geographical distribution [123]. This study estimated the magnitude of soil loss rate, assessed the transformation of erosion risk, and identified priority areas for conservation. The estimated total soil loss in the Gobele Watershed area was 1,390,130.48 tons in 2000 and 1,022,445.09 tons in 2016. The findings of this study are within the range of the early findings that estimated the annual soil loss rate in the highland areas of Ethiopia from 1248 to 23,400 million tons [40]. Significant variations in annual soil loss rates were also reported in various parts of Ethiopia. Based on an assessment of soil erosion in the Koga watershed, Northwestern Ethiopia, Gelagay and Minale [124] reported that the total annual soil loss accounted for 255,283 tons. Ayalew [33] estimated soil loss and sediment in the Zingin watershed of the Ethiopia highlands for conservation planning and reported a total annual soil loss potential accounts for 57,750.15 tons.
The present study found that the mean annual soil loss accounts for 51.04 t ha −1 y −1 , 34.26 t ha −1 y −1 in 2000 and 2016, respectively. The estimated mean rate of annual soil loss in the study area is much greater than that of the maximum tolerable soil loss estimate at a national scale (18 t ha −1 y −1 ) [33,125], and to the normal soil loss tolerances (from 5 t ha −1 y −1 to 11 t ha −1 y −1 ) [54,126]. It is also higher than the findings of earlier researchers in other areas in Ethiopia [27,33,127], and elsewhere in the world [128,129]. On the contrary, the present study results are lower than the findings across the Northwestern Ethiopia highland [130], and other similar studies' findings [131][132][133]. Several previous studies suggested that vulnerability to a higher rate of water-induced soil erosion occurring in eastern Ethiopia's highland is associated with the adverse effects of LULC changes, an unsustainable land management, and less emphasis is given to SWC practices [7,134,135]. These events resulted in an accelerated on-site soil nutrient loss and off-site sediment accumulation in downstream areas [135]. A complete drying up of Lake Alemaya in the East Hararghe Ethiopian highland over the recent years is one piece of compelling evidence that justifies the negative consequences of soil erosion, in combination with an unsustainable land use management practice [136][137][138].
Based on the estimated rates of mean annual soil loss, erosion risk was grouped into eight classes ranging from the very low to extremely high. We derived thematic maps to show the spatial distribution of erosion risk classes ( Figure 6). The very low dominated major parts of the watershed area in 2000 and 2016 (Table 4) account for 77.52%, 87.02% of the total study area, respectively. The high, very high, and extremely high erosion risks have declined by 0.23%, 0.21%, and 0.05% of the total study area, respectively. Overall, the magnitude of soil loss rates across the Gobele Watershed marginally decreased over the period between 2000 and 2016. The decline in soil loss rates is probably due to some conservation measures were taken by the local people over the recent years. This implies that the local communities within the study area have well recognized the possible environmental threats posed by the water-induced soil erosion and because of this, putting efforts into SWC practices in erosion-prone areas. Supporting the current study's findings, Tadesse et al. [139] noted considerable improvement in the soil erosion and rehabilitation of degraded lands in the Yezat Watershed of North Western Ethiopia following implementation of integrated watershed development programs. Validating these findings, Akale et al. [140] also reported that implementation of upland conservation measures significantly reduced the surface runoff and increased base flow in the Guale and Tikur-Wuha watersheds of the Ethiopian highlands.
In the present study area, contour bounds, terracing, and check-dams are some of the promising soil erosion control practices adopted by the local communities as observed during the field visit. We have also seen that plantation of trees in up-slope areas, especially Khat (Catha edulis) as a perennial cash crop and raw intercropping with maize and sorghum is largely practiced in various parts in the study area, and substantially contributed to the diminishing of soil losses by water erosion (Figure 8). This agrees with Lemessa [141], who highlighted the role of Khat as the most important erosion control practices by minimizing rainfall-runoff velocity in most of the Eastern and Western Hararghe highlands of Ethiopia. Confirming this, Li et al. [67] and Efe et al. [131] pointed out that the existence of vegetation covers reduces the adverse effect of slope gradients on soil erosion intensity. Moreover, the recent study by Cerdà et al. [50] revealed that application of the straw mulch as erosion control practices has been shown to reduce the rainfall-runoff velocity and connectivity of the flows, subsequently decreasing the soil loss rates and sediment fluxes. Keesstra et al. [51] emphasized the nature-based solutions (NBSs), including the soil solutions and landscape solutions as a cost-effective and sustainable conservation strategies for mitigating and restoration of degraded land and improve ecosystem services. Though the SWC implemented over the recent years might have reduced the magnitude of soil loss rate, the eroded areas are widely observed in the study area (Figure 8a-c).
The change analysis results of the present study showed that about 70.80% of erosion risk areas occupied in 2000 continued under the same classes in 2016. The results show that the erosion risk area increased by 19.67% of the total study area, and decreased by 9.53%, which shows that overall soil loss situation is worsening in the study area. These findings agree with those of the recent study by Uddin et al. [90] who reported that the state of soil erosion risk in the Koshi Basin has been worsening following increase in the proportion of the eroded areas over 9.0% of the total basin area between 1990 and 2010. On the contrary, Wang et al. [92] found improvement in the state of erosion risk in the Danjiangkou reservoir area, China, where the eroded areas have declined from 32.1% in 2004 to 25.43% in the 2010 study period. Moreover, Jiang et al. [16] reported that the eroded area has decreased by 61% between 2000 and 2012, despite increases in the intensity of soil erosion by 39% in some areas of the Mount Elgon region, Uganda.
practices has been shown to reduce the rainfall-runoff velocity and connectivity of the flows, subsequently decreasing the soil loss rates and sediment fluxes. Keesstra et al. [51] emphasized the nature-based solutions (NBSs), including the soil solutions and landscape solutions as a cost-effective and sustainable conservation strategies for mitigating and restoration of degraded land and improve ecosystem services. Though the SWC implemented over the recent years might have reduced the magnitude of soil loss rate, the eroded areas are widely observed in the study area (Figure 8a-c).  19.67% of the total study area, and decreased by 9.53%, which shows that overall soil loss situation is worsening in the study area. These findings agree with those of the recent study by Uddin et al. [90] who reported that the state of soil erosion risk in the Koshi Basin has been worsening following increase in the proportion of the eroded areas over 9.0% of the total basin area between 1990 and 2010. On the contrary, Wang et al. [92] found improvement in the state of erosion risk in the Danjiangkou reservoir area, China, where the eroded areas have declined from 32.1% in 2004 to 25.43% in the 2010 study period. Moreover, Jiang et al. [16] reported that the eroded area has Therefore, since there is geospatial variation in soil erosion risk distribution across the watershed area, identification of priority area is the key factor for planning and implementing appropriate SWC [10][11][12][13][14][15][16][17][18][19]51,90,142]. Accordingly, we prioritized areas with a higher and increasing soil loss rate as SWC priority areas (Table 6). This can play a vital role in supporting decision makers and conservation planners to design proper SWC measures based on the severity levels of soil loss in the Gobele Watershed. The first-, second-and third-priority levels account for about 3232.79 ha (1.36% of the watershed area) that need urgent SWC measures. Most of the top three priority areas are spatially situated in Kersa and Kurfa Chele in the northwest, Girawa in the southwest, Haramaya in north, and Fedis in the south and southeast. Moreover, these districts are situated in the steep slope area with slope gradients greater than 30 meters. In line with the findings of our study, Abate [5] reported that the steep slopes are attributed to high soil loss rate and aggravated soil erosion risk in the Borena Woreda of South Wollo Highlands of Ethiopia. A recent report by Karamage et al. [133] for the Nyabarongo River Catchment, Rwanda, also agrees with the present study's findings, showing a higher soil erosion in the very steep slope areas contributing about 73.5% of the total soil loss. They strongly suggested construction of erosion control structures and rehabilitation of vegetation covers as a solution to reduce the soil loss in the catchment area [133]. The fourth, fifth and sixth conservation levels accounting for 21,104.37 ha (8.88% of the total study area) need of minor SWC measures. The last two conservation levels represent areas with lower soil loss and erosion risk accounting for 213,449.27 ha (89.76% of the total area), which currently may not need the emergency SWC measures.

Conclusions
This study estimated annual soil loss rate, assessed the spatial change of erosion risk and mapped priority areas for SWC measures by taking into consideration the severity levels of soil loss across the Gobele Watershed, East Hararghe Zone, Ethiopia. We used the RUSEL model developed in an ArcGIS software environment. The model results offered a reliable quantitative estimate of water-induced soil loss rates and spatial distribution of erosion risk in the study area. The estimated total soil loss at the  51.04 t ha −1 y −1 and 34.26 t ha −1 y −1 , respectively. Thus, the status of soil loss across the Gobele Watershed has slightly improved over the study period partly due to some conservation practices. Eight erosion risk classes were mapped based on the mean soil loss rates showing that a larger portion of the study area are found under the very low, low and low medium classes and only a small part found within the high, very high and extremely high. Despite a reduction in the magnitude of soil loss rates, the spatial dimension of erosion risk is getting worse in the study area. The change analysis results showed that overall, about 70.80% of the erosion risk areas covered in 2000 continued under the same erosion risk classes in 2016 but 19.67% increased and 9.53% decreased in the total study area. This reveals that the soil erosion risk situation is worsening in the study area. We identified and mapped areas with a high and increasing soil erosion risk as conservation priority areas based on the MCDR method. Thus, conservation priorities identified in the present study can serve as a spatial decision support tool, and as input for decision makers and conservation planners for future intervention measures in highly affected areas. It also supports efforts to minimize environmental and economic impacts of soil erosion and offer insights on policy implications on what should be done to establish sustainable watershed management practices in the study area. The emperical RUSEL model and MCDR method outlined in ArcGIS environment for estimating the soil loss rates, assess the spatial change of erosion risk, and prioritization of SWC in the Gobele Watershed can be applied to other similar watersheds following proper validation.
Supplementary Materials: The following are available online at http://www.mdpi.com/xxx/s1. Table S1: Mean annual Rainfall (mm) obtained from National Meteorological Agency (NMA) from 1999 to 2015 and Rainfall erosivity (R) factor. Table S2: Attributes of Soil units and calculated soil erodibility (K) factor in the Gobele Watershed, East Hararghe Zone, Ethiopia.