Potential Soil Loss Estimation and Erosion-Prone Area Prioritization Using RUSLE,GIS,and Remote Sensing in Chereti Watershed,Northeastern Ethiopia

Abstract

Soil erosion by water is the major form of land degradation in Chereti watershed, Northeastern Ethiopia. This problem is exacerbated by high rainfall after a long period of dry seasons, undulating topography, intensive cultivation, and lack of proper soil and water conservation measures. Hence, this study aimed to estimate the 23 years (1995-2018) average soil erosion rate of the watershed and to identify and prioritize erosion-vulnerable subwatersheds for conservation planning. The integration of the revised universal soil loss equation (RUSLE), geographic information system, and remote sensing was applied to estimate the long-term soil loss of the watershed. The RUSLE factors such as rainfall erosivity (R), soil erodibility (K), topography (LS), cover and management (C), and support and conservation practices (P) factors were computed and overlayed to estimate the soil loss. The result showed that the annual soil loss rate of the watershed ranged up to 187.47 t ha⁻¹ year⁻¹ in steep slope areas with a mean annual soil loss of 38.7 t ha⁻¹ year⁻¹, and the entire watershed lost a total of about 487 057.7 tons of soil annually. About 57.9% of the annual watershed soil loss was generated from 5 subwatersheds which need prior intervention for the planning and implementation of soil conservation measures. The integrated use of RUSLE with GIS and remote sensing was found to be indispensable, less costly, and effective for the estimation of soil erosion, and prioritization of vulnerable subwatersheds for conservation planning.

Keywords

Soil erosion RUSLE GIS remote sensing Chereti watershed

Introduction

Soil erosion by water is the detachment and removal of fertile topsoil as a result of rainfall and runoff causing the soil to deteriorate.¹ It is the global environmental problem that negatively affects the productivity of the natural ecosystem and agriculture.^2
-4

Physical land degradation, especially soil erosion, depletes soil fertility, removes organic matter, and leads to loss of the topsoil that provides water and nutrient holding capacity.⁵ The overall effect is to threaten the ecosystem services, and the sustainable development goals (SDGs) adopted by the United Nations (UN) in 2015.^6,7 In reversing the challenges, land degradation neutrality (LDN) and land restoration, which is part of SDGs could be an alternative to stop net erosion or bringing the erosion rate in equilibrium with the soil formation rate.⁸

In Ethiopia, where agriculture is the mainstay of the economy, soil erosion is one of the major and continuous forms of land degradation that affect the sustainability of agricultural production in the country.^9
-12 The problem is particularly severe in the highlands of the country owing to the high erosive power of rainfall, intensive agricultural practices, and high population density.^3,9,12,13

Different recent studies^3,9
-12 were undertaken to estimate soil loss in different parts of Ethiopian highlands, and all reported that soil erosion by water is the major problem that impends agricultural productivity and sustainability. Soil erosion in the highlands is among the many factors that affect rural livelihoods by reducing land productivity and aggravating poverty and food insecurity.^9,12

Chereti watershed is a part of Ethiopian highlands where soil erosion by water is the major problem exacerbated by different factors such as high rainfall after a long period of dry seasons, undulating terrain, intensive cultivation, and lack of proper soil and water conservation measures. The problem is particularly severe in croplands and steep slope areas. Unlike other areas in the highlands, no research has been done so far on soil loss estimation in the watershed.

The estimation of soil erosion and mapping its spatial variation at the watershed level is crucial for decision-makers and land resource managers who are involved in reducing soil resource degradation. Different methods were developed to estimate the rate of soil erosion and erodibility at field and watershed levels. According to Pham et al,¹⁴ universal soil loss equation (USLE), soil and water assessment tool (SWAT), and water erosion prediction project (WEPP), the erosion productivity impact calculator (EPIC), the agricultural nonpoint source (AGNPS) are some of the tools developed by researchers to estimate soil loss. Previous studies employed geospatial techniques such as geographic information systems (GIS) and remote sensing to investigate soil erodibility and soil erosion in different regions of the world either integrated with these and other models or independent applications. Puno et al¹⁵ applied GIS and Geospatial-interface Water Erosion Prediction Project (GeoWEPP) model to predict soil erosion and sediment yield by preparing 4 input files corresponding to climate, slope, land management, and soil properties in Taganibong watershed, Philippines. Panagos et al¹⁶ developed a high-resolution soil erodibility map of Europe through integrating soil texture and organic carbon from Land Use/Cover Area Frame Survey (LUCAS) data with spatial data such as latitude, longitude, remotely sensed, and terrain features. Alexakis et al¹⁷ used satellite remote sensing observations (Sentinel-2 and Landsat 8), field spectroscopy, artificial neural networks (ANN), soil textural and chemical analysis, and GIS to estimate soil organic matter, calcium carbonate equivalent (CaCO₃), and soil erodibility (K) factor to be used as the main data input in the Revised Universal Soil Loss Equation (RUSLE) model to estimate the soil erosion risk in the Akrotiri cape in Chania, Crete, Greece. Shi et al¹⁸ applied visible-near infrared (Vis-NIR) spectroscopic assessment to predict the soil aggregate stability (resistance of the soil to external erosive forces) in the Belgian Loam Belt, Belgium, and found that Vis-NIR spectroscopy is a promising technique that enables fast and efficient large-scale assessment of soil’s resistance to external erosive forces. Other authors^19
-21 also employed the integration of GIS and remote sensing with SWAT model to estimate soil erosion.

Universal soil loss equation developed by Wischmeier and Smith in 1965 and later revised into RUSLE by Renard et al in 1997 is the most commonly used empirical model to estimate soil erosion in areas where there is a lack of measured data.^14,22
-26 The RUSLE model predicts the long-term average annual rate of erosion based on rainfall patterns, soil type, topography, cropping system, and management practices.^2,14,27 It predicts annual soil loss as a product of factors such as rainfall erosivity (R), soil erodibility (K), topography (LS), cover and management (C), and support practice (P).² Angima et al²⁸ employed the RUSLE model to predict and map soil erosion risk of Kianjuki Catchment in the central Kenya using GIS.

Different previous studies^{3,9
-13,29,30} in Ethiopia have also used the combined application of GIS and remote sensing with the RUSLE model in predicting the rate of soil erosion. However, many of these studies have employed single-date satellite image acquired during the dry season for both cover and management (C) factor and support and conservation practice (P) factor estimation of RUSLE, with less consideration on the effect of high vegetation and crop cover, and high management practices on soil erosion during the summer season. Also, these previous studies adopted manually assigning C-factor value from the values suggested and used in different published studies after the classification of the satellite image, with less emphasis on the application of Normalized Difference Vegetation Index (NDVI) to estimate C-factor during the prediction of soil loss using RUSLE.

This study was conducted in Chereti watershed using the integration of RUSLE, GIS, and remote sensing techniques with the aims (1) to estimate the long-term average soil loss rate of the watershed and (2) to identify and prioritize erosion vulnerable subwatersheds for planning and implementation of conservation measures. To achieve our research objectives, unlike other similar studies in Ethiopian highlands, we have used 2 Landsat 8 OLI images acquired during winter and summer seasons to estimate P-factor and C-factor, respectively. We have also employed NDVI to compute the C-factor value of RUSLE using the Landsat 8 OLI image data acquired during the summer season.

Materials and Methods

Study area

Chereti watershed is located in the North Wollo Zone of Amhara National Regional State, Northeastern Ethiopia (Figure 1). It lies between 11°42’02”N and 11°48’41”N latitude, and 39°31’41”E and 39°42’43”E longitude. The watershed is a part of Northern highlands with an area coverage of 125.82 km². The landscape of the watershed is quite diversified. The topography of the watershed is characterized by gently sloping flat terrains, steep to very steep hillside slope, and deeply incised V-shaped valleys with an altitude range from 1569 m.a.s.l at the outlet of the watershed to 3446 m.a.s.l in the southwestern border of the watershed. The watershed is dissected by several small tributaries. Chereti is the largest and longest perennial river in the watershed that drains throughout the year. It flows from the southwest part of the watershed to the northeast direction. There are also several small tributaries in the watershed flowing toward this largest river.

Figure 1.

Location map of the study area.

Woyina Daga (cool humid highlands), Dega (temperate, subhumid highlands), and Wurch (cold highlands) are the 3 agroclimatic zones found in the watershed. The majority of the watershed area (81.55%) is characterized by Woyina Daga agroclimatic zone (cool and humid highlands). The remaining 18.05% and 0.4% of the watershed area are characterized by the Dega agroclimatic zone (temperate, subhumid highlands) and Wurch (cold highlands) agroclimatic zone, respectively. The Dega agroclimatic zone is found in northern and southwestern parts of the watershed, and the Wurch (cold highlands) agroclimatic zone is found in the southwestern border of the watershed.

The 23 years (1995-2018) rainfall data of 4 stations obtained from the National Meteorological Agency of Ethiopia shows that the mean annual rainfall of the watershed varies from 946.17 to 1371.29 mm. The watershed receives its maximum rainfall during the summer season (locally known as kiremt) which extends between June to September. The annual mean minimum and mean maximum temperatures of the watershed are 13.94°C and 26.76°C, respectively, with maximum temperature in the watershed occurs from February to May.

The digital soil map obtained from Ethiopia Ministry of Water, Irrigation, and Energy shows that Eutric Regosols, Eutric Cambisols, Lithosols, and Vertic Cambisols are the soil types that existed in the watershed (). Eutric Cambisols is the dominant soil type which covers an area of about 6242.7 ha (49.6%) from the total area of the watershed. Relatively flat areas and cultivated lands of the watershed are largely covered by Eutric Cambisols.³¹ Eutric Regosols, Lithosols (now are Leptosols),³¹ and Vertic Cambisols are found in the hilly and upstream parts of the watershed.

From the field observation and later supported by the classification of Landsat 8 satellite image, it was found that forest land, cropland, sparse shrubland, dense shrubland, built-up area, and river beds are the land use and land cover (LULC) types in the watershed (Table 4). Cropland shares the greatest size from all other land-use types which account for an area of 7003.3 ha (55.7%) followed by dense shrubland (area dominantly occupied by densely grown short trees) which accounts for about 3906.2 ha (31.1%) of the total area of the watershed. Climatically, the watershed is conducive for the production sorghum (Sorghum bicolour) and Ethiopian Teff (Eragrostis tef) which are the 2 major crops grown in the watershed during the rainy season. During the field visit, it was also observed that farmers practice small-scale irrigation during the dry season through diverting water from Chereti river.

Data types and sources

Soil erosion by water is affected by topography, vegetation cover, climate, soil condition, and other management practices. Hence, the assessment of soil erosion requires data related to these factors. The required data for this study were collected from different sources. Landsat 8 Operational Land Imager (OLI) satellite images with a spatial resolution of 30 m acquired on February 1, 2019, and July 11, 2019, were downloaded from US Geologic Survey website. The 2 Landsat 8 OLI images have been downloaded with L1TP correction level or Level-1 precision and terrain corrected product level (provides radiometric and geodetic accuracy by incorporating ground control points [GCPs]). However, to further enhance the quality of the Landsat images, atmospheric corrections such as haze reduction, noise reduction, and histogram equalization were performed for the 2 Landsat 8 OLI images after layer stacking the spectral bands. The first Landsat image acquired on February 1, 2019, was used to prepare the LULC map for P-factor estimation. The acquisition date of this Landsat image was purposively selected to easily distinguish the LULC types during the classification process since the season is dry with the lowest percent or zero monthly cloud cover. The second Landsat image acquired on July 11, 2019, was used to compute the NDVI map for C-factor estimation. To better understand the effect of crop and management (C) factor on the rate of soil erosion, it is important to consider the season when vegetation and crop cover is high. Hence, the Landsat image acquired on July 11, 2019, was purposively considered because it is a relatively cloud-free Landsat image acquired during the summer season when vegetation and crop cover in the watershed is high.

STRM GDEM Digital Elevation Model (DEM) with 30-m spatial resolution was also downloaded from US Geologic Survey website (Table 1) to prepare the slope map for LS and P-factor estimation. The 23 years (1995-2018) precipitation data of 4 meteorology stations were obtained from the National Meteorological Agency of Ethiopia. The digital soil map of the watershed was obtained from the Ethiopian Ministry of Water, Irrigation, and Energy. For supervised classification of Landsat 8 satellite image and accuracy assessment, ground truth data of different LULC types were collected from the field using Global Positioning System (GPS GARMIN 60×) receiver. Reference data from inaccessible areas within the watershed were digitized from Google Earth.

Table 1.

Data types and sources used in the study.

Data types	Data sources
Digital soil map of the watershed (1:250 000 scale)	Ethiopian Ministry of Water, Irrigation, and Energy
23 years monthly rainfall data (1995-2018)	National Meteorological Agency of Ethiopia
STRM DEM digital elevation model (30-m spatial resolution)	Downloaded from US Geologic Survey (http://earthexplorer.usgs.gov/)
Landsat 8 satellite image (30-m spatial resolution) (acquired on February 1, 2019)
Landsat 8 satellite image (30-m spatial resolution) (acquired on July 11, 2019)
Ground control points (GCP) from LULC types	Collected from the field using GARMIN 60× GPS receiver (February 18-21, 2019)
Reference data of LULC types from inaccessible areas of the watershed	Digitized from Google Earth (February 25, 2019)

Methods

The integration of RUSLE (equation (1)) with GIS and remote sensing was used to estimate the soil loss of the watershed. Spatial data layers of the 5 factors of RUSLE such as rainfall erosivity (R) factor, soil erodibility (K) factor, topographic (LS) factor, cover and management (C) factor, and support and conservations practices (P) factor were prepared using GIS and remote sensing techniques from the data collected from different sources. ArcGIS 10.4 and ERDAS IMAGINE 2014 software were exhaustively used to process, create, and overlay digital data layers of each factor. Average annual soil loss of the watershed was computed by multiplying the 5 factors using the RUSLE equation (equation (1))²⁴ given below

A = R * K * L S * C * P

(1)

where, A is estimated annual soil loss (t ha⁻¹ year⁻¹); R is rainfall erosivity factor (MJ mm ha⁻¹ h⁻¹ year⁻¹); K is soil erodibility factor (t ha⁻¹ MJ⁻¹ mm⁻¹); L is slope length and S is slope steepness factor (dimensionless); C is cover and management factor (dimensionless), and P is support and conservation practice factor (dimensionless).

Methods of RUSLE factors estimation

Rainfall erosivity (R) factor

The rainfall erosivity factor (R) describes the erosivity of rainfall at a particular location based on the rainfall amount and intensity and reflects the effect of rainfall intensity on soil erosion.^23,27,32 According to Wischmeier and Smith²⁷ and Yin et al,³³ rainfall erosivity (R) factor is the product of kinetic energy of rainstorm and its maximum 30-min intensity (I₃₀). However, due to the lack of rainfall kinetic energy and rainfall intensity (I₃₀) data for the study area, the empirical equation that estimates R-value from total annual rainfall suggested by Hurni,³⁴ for Ethiopian condition was adopted to compute the rainfall erosivity (R) factor of the watershed. This equation was also adopted by similar studies^3,13,35,36 in Ethiopia. The mean annual rainfall of 4 meteorology stations (Sirinka, Woldia, Hara, and Mersa) was calculated (Table 2) from raw monthly precipitation records of 23 years (1995-2018) obtained from the National Meteorological Agency of Ethiopia. Since only one meteorology station (Sirinka) is available within the watershed area, the rainfall record of the other 3 nearby meteorology stations to the watershed such as Woldia, Hara, and Mersa was used for R-factor estimation. The mean annual rainfall of the 4 stations was interpolated by inverse distance weighted (IDW) method in ArcGIS10.4 environment with a cell size of 30 m and clipped using the watershed boundary shapefile to generate a continuous rainfall map of the watershed. Yin et al³³ reviewed as about one-quarter of the studies from the Scopus and Google Scholar databases used Inverse Distance Weighting (IDW) spatial interpolation technique to generate rainfall erosivity maps based on the analysis of 51 studies from the United States. The R factor was calculated from the interpolated continuous rainfall raster map by applying the equation (equation (2))³⁴ given below in Raster Calculator of ArcGIS 10.4 Spatial Analyst Extension.

R = - 8.12 + (0.562 * P)

(2)

where R is the rainfall erosivity factor (MJ mm ha⁻¹ h⁻¹) and P is the mean annual rainfall (mm). A similar equation but for United States condition was reviewed as a method for R-factor determination using annual precipitation.³³

Table 2.

Meteorology stations and calculated mean annual rainfall (mm).

Station name	Latitude	Longitude	Mean annual rainfall (mm) (1995-2018)
Sirinka	11.75	39.61	1072.91
Woldia	11.83	39.59	1417.83
Hara	11.84	39.74	824.72
Mersa	11.66	39.67	1033.81

Soil erodibility (K) factor

The soil erodibility (K) factor reflects the effect of soil properties and profile characteristics on soil loss.^11,14,24 It is the measure of the potential susceptibility of soil to the detachment and transport by rainfall and runoff.³⁷ The main soil properties influencing the K factor are soil texture, organic matter, soil structure, and permeability of the soil profile.^{4,10,11,13,25,32,38} Since determining theses soil properties from the soil map without soil sampling and analysis is practically difficult, the K-factor of the watershed for this study was computed based on soil type and soil color method suggested by Hurni³⁴ for the Ethiopian condition as adopted by several studies^{3,30,36,39,40} in Ethiopian highlands. Primarily, the soil types of the watershed were identified through extraction using ArcGIS 10.4 software from Digital Soil Map of Ethiopia (1:250 000 scale) obtained from the Ethiopian Ministry of Water, Irrigation and Energy (EMWIE). Then, the K-value for each soil type was assigned in the extracted vector soil map of the watershed depending on the type and color of soil as suggested by Hurni,³⁴ and the vector soil map was converted into raster format with 30-m cell size using its K value for final estimation of soil loss.

Topographic (LS) factor

The Topographic (LS) factor reflects the effect of surface topography or slope length (meter) and slope steepness on the rate of soil erosion by water.^14,25,32,41 For the computation of the topographic (LS) factor of the watershed, the LS factor estimation equation (equation (3)) suggested by Moore and Burch,⁴² was adopted. This equation is also adopted in different studies.^{13,23,35
-37,43} for LS factor calculation. The flow accumulation map was generated from the DEM (30-m spatial resolution) after filling the DEM and computing flow direction using ArcHydro tools in ArcGIS 10.4 software. The slope (degree) map of the watershed was derived from the DEM map using the Spatial Analyst extension in ArcGIS 10.4 environment (Figure 4A).

Finally, the flow accumulation map and slope map were integrated to prepare LS factor map using the equation (equation (3)) given below in the raster calculator of ArcGIS Spatial Analyst Extension:

L S = {(\frac{F l o w A c c u m u l a t i o n * C e l l s i z e}{22.12})}^{0.4} * {(\frac{S i n θ}{0.0896})}^{1.3}

(3)

where LS is slope length and slope steepness factor, flow accumulation denotes the accumulated upslope contributing area for a given cell, cell size is the pixel size of the DEM (30 m for this study), $θ$ is slope steepness in degree, and 22.12 is the RUSLE unit plot length in meter.

Crop and management (C) factor

Crop and management (C) factor reflect the effect of LULC, cropping, and management practices on the rate of soil erosion,^12,14,25 and it is the ratio of soil loss from land covered by vegetation to the corresponding loss from continuous fallow.^9,23 The C-factor value ranges between 0 implies very strong cover effects and 1 indicates no cover present and the surface treated as barren land.^11,14

Different researchers^{3,9,12,32,44,45} assigned published C values after classifying remotely sensed data into different LULC classes. Others authors^{1,4,23,25,43,46} employed NDVI, one of the vegetation indices, introduced by Rouse et al,⁴⁷ that measures the amount of green vegetation from satellite images based on the characteristic reflectance patterns of green vegetation. Normalized Difference Vegetation Index is the ratio of the difference between Near Infrared and RED band reflectance to the sum of Near Infrared and RED band reflectance (equation (4)); and the values of NDVI ranges from −1 to 1 where higher values are for green vegetation, closest to 0 for bare soil, and negative values for water bodies.⁴⁸

For this study, the NDVI was used to estimate and map the crop and management (C) factor value of the watershed. Normalize Difference Vegetation Index was calculated from the Landsat 8 satellite image by applying the equation (equation (4)) given below using ERDAS IMAGINE 2014 software.

NDVI = \frac{(NIR - RED)}{(NIR + RED)}

(4)

where NDVI is Normalized Difference Vegetation Index, NIR is surface spectral reflectance in the near-infrared band, and RED is surface spectral reflectance in the red band.

After calculating the NDVI of the watershed, the following equation (equation (5)) proposed by Durigon et al,⁴⁹ for tropical regions and adopted by de Carvalho et al,⁵⁰ was employed to estimate the C-factor value of the watershed from the calculated NDVI. According to Phinzi and Ngetar⁵¹ and Almagro et al,⁴⁴ this method is more realistic to determine C-factor values in tropical regions. Raster Calculator of Spatial Analyst Extension in ArcGIS 10.4 software was used to compute and generate the C-factor map.

C_{r A} = 0.1 (\frac{- NDVI + 1}{2})

(5)

where C_rA is the crop and management (C) factor.

The Landsat image with 30-m spatial resolution acquired during rainy season on July 11, 2019, was used to calculate the NDVI and to generate the C-factor value because soil erosion by water and vegetation cover is high during this season.

Support and conservation practice (P) factor

Support and conservation practice (P) factor refers to the effects of conservation practices in reducing the quantity and rate of runoff and the amount of soil erosion.^13,27,32 The P value can be determined by the type of conservation measure implemented,³ and the value is always between 0 and 1, where the value 0 indicates a good erosion-resistant facility made by man, and the value 1 indicates the absence of an erosion-resistant facility.^25,32 However, it is very difficult to estimate the P-factor value when there is a lack of data and no permanent support and conservation practices. The method suggested by Wischmeier and Smith²⁷ was used to estimate the support and conservation practice (P) factor value of the study watershed because there is little information on support and conservation practices in the watershed. This method is an alternative method that uses the combination of LULC data and land slope to estimate the P-factor value when there is a lack of data and no permanent management practices in the study area.^3,13,27 Similarly, the method was adopted by^{3,13,27,36,40} to determine P-factor value.

The LULC map of the watershed was derived through maximum likelihood supervised image classification using ERDAS IMAGINE 2014 software from the Landsat 8 OLI satellite image with 30-m spatial resolution acquired on February 1, 2019. Before image classification, field observation was conducted from February 18, 2019, to February 21, 2019, to have a clear understanding of the main categories of LULC types and to collect representative GCPs (Figure 5A) for the image classification process. Image preprocessing such as layer stacking the spectral bands, and image enhancement technique (atmospheric correction) were performed to improve the quality of the satellite image before classification.

About 116 representative reference points of LULC types were used for the accuracy assessment of the classified image. From the total 116 reference points, about 95 GCPs were collected using a GPS receiver with GPS accuracy of 5 to 7 m, and the remaining 21 reference points of LULC types were digitized from Google Earth on February 25, 2019. More specifically, 19 reference samples from forest land (11 samples digitized from Google Earth), 14 samples from dense shrubland (8 samples digitized from Google Earth), 13 samples from sparse shrubland (2 samples digitized from Google Earth), 45 samples from cropland, 11 samples from the river bed, and 14 samples from the built-up area were collected, and Maximum Likelihood supervised image classification and was done using these reference points. Subsequently, accuracy assessment was done and the confusion matrix (Table 6) was developed using an accuracy assessment tool in ERDAS IMAGINE 2014 software to measure the reliability of the LULC classification.

After the processes of Landsat 8 satellite image classification and, the LULC map of the watershed were broadly categorized into agricultural and other land uses, and a P value of 1 was assigned for other land uses. As it was suggested by Wischmeier and Smith,²⁷ the agricultural land use was reclassified into 6 classes based on the slope (%) of the land, and the respective P value for each class was assigned (Table 6). The classified agricultural land use map based on slope and other land use maps were overlayed after converting into vector format and assigning respective P value. Finally, the overlayed map was converted into a raster format with 30-m pixel size using its P value to make it suitable for pixel-by-pixel overlay analysis to estimate soil erosion (Figure 5D).

After the estimation of all RUSLE factors and preparation of each factor in raster format, the raster layers were overlayed together (Figure 2) using the Raster Calculator of Spatial Analyst Extension in the ArcGIS environment to drive the final soil erosion estimation map. Zonal Statistics as table tool of Spatial Analyst Extension in ArcGIS 10.4 software was used to calculate the average soil loss value of the watershed.

Figure 2.

Flow chart of the methods used to estimate soil erosion.

Approach for validation of model results

Due to poorly available data in the study area, the validity and consistency of the model output was compared with the output of previously published study with similar approaches in Ethiopian highland as used by previous studies.^3,9,12,36,52 Field observation was also conducted to check the model outputs in the study area.

Methods of subwatersheds erosion vulnerability mapping and prioritization

The whole watershed was classified into 18 subwatersheds which were delineated using ArcSWAT Spatial Analyst Extension of ArcGIS 10.4 software, and the mean soil loss value of each subwatershed was estimated using Zonal Statistics as table of Spatial Analyst Extension in ArcGIS 10.4 software after extraction from the soil loss map of the whole watershed. The total annual soil loss of each subwatershed and their contribution to the total annual soil loss of the whole watershed was calculated, and priority level was assigned for each subwatershed based on their mean soil loss value.

Result

RUSLE factors estimation and map preparation

Rainfall erosivity (R) factor estimation

The result of IDW interpolation using precipitation data of 4 meteorology stations showed that the mean annual rainfall of the watershed is ranged from 946.17 to 1371.29 mm with the highest value in the northern part of the watershed (Figure 3A). About 55.4% of the watershed area receives a mean annual rainfall from 1106.22 to 1371.29 mm, and the remaining 44.6% of the watershed area receives 946.17 to 1106.22 mm with the lowest value near the outlet of the watershed.

Figure 3.

Mean annual rainfall map of the watershed (A), rainfall erosivity (R) factor map (B), map of watershed soil types (EMWIE, 2019) (C), and soil erodibility (K) factor map (D).

The rainfall erosivity (R) factor map was derived by applying equation (2) in Raster Calculator of Spatial Analyst Extension in ArcGIS 10.4 software from the interpolated rainfall map, and the result indicated that the R-factor value of the watershed ranged from 523.63 to 762.55 MJ mm ha⁻¹ year⁻¹(Figure 3B) with higher values occurring in the northern part of the watershed, and the potential of rainfall to erode soil gradually decreases toward the south and southeast part of the watershed.

Soil erodibility (K) factor estimation

The major soil types in the watershed are Eutric Cambisols, Eutric Regosols, Vertic Cambisols, and Lithosols with the largest area (49.6%) occupied by Eutric Cambisols soil type. The result in Table 3 and Figure 3D showed that soil erodibility (K) factor value of the watershed varied from 0.15 t ha⁻¹ MJ⁻¹ mm⁻¹ in areas with Vertic Cambisols to 0.20 t ha⁻¹ MJ⁻¹ mm⁻¹ in areas with Eutric Cambisols, Eutric Regosols, and Lithosols soil types. The result indicated that about 91.9% of the watershed area has a K-value of 0.20 t ha⁻¹ MJ⁻¹ mm⁻¹ mainly dominated by croplands with Eutric Cambisols.

Table 3.

Soil types of the watershed and their respective erodibility (K) factor value.

Watershed soil type	Area		Soil color	K-factor	References
Watershed soil type	ha	%	Soil color	K-factor	References
Eutric Cambisols	6242.7	49.6	Brown	0.20	Hurni³⁴; Miheretu and Yimer⁵²; Zerihun et al,³⁰; Tilahun et al,⁵³
Eutric Regosols	2850.2	22.7	Brown	0.20	Hurni³⁴; Miheretu and Yimer⁵²
Vertic Cambisols	1019.2	8.1	Black	0.15	Gelagay and Minale³⁶; Hurni³⁴
Lithosols	2470.2	19.6	Brown	0.20	Hurni³⁴; Miheretu and Yimer⁵²

Topographic (LS) factor estimation

The slope steepness map (Figure 4A) was directly generated from the STRM DEM (30-m spatial resolution) map of the watershed. The flow accumulation map was also derived from the STRM DEM map after filling small imperfections and processing flow direction. The topographic (slope length and slope steepness) factor was computed based on equation (8) using the Raster Calculator of Spatial Analyst Extension in the ArcGIS environment.

Figure 4.

Slope (degree) map of the watershed (A), topographic (LS) factor map (B), NDVI map of the watershed (C), and crop and management (C) factor map (D).

The result in Figure 4B indicated that the slope length and slope steepness (LS) factor value of the watershed varies from 0 in flat areas to 11.81 in stream bank and steep slope areas of the watershed with a mean LS value of 5.09 and standard deviation of 3.31.

Crop and management (C) factor estimation

Normalized Difference Vegetation Index map of the watershed was computed using equation (4) in ERDAS IMAGINE 2014 software from the satellite image acquired on July 11, 2019, and the result showed that the NDVI value of the watershed ranged from −0.015 to 0.595 with the highest values in the forest and dense shrubland areas (Figure 4C).

The C-factor value of the study area was calculated from the NDVI of the watershed using equation (5) in Raster Calculator of ArcGIS Spatial Analyst extension, and the result showed that the C-factor value of the watershed ranged from 0.02 to 0.05 with the lowest value in forest areas of the watershed (Figure 4D) which implies higher cover effect in the forest and dense shrubland areas.

Support and conservation practice (P) factor estimation

Image classification and classification accuracy assessment

Maximum likelihood supervised image classification was performed using a total of 30 training sites also known as area of interest (AOI) and 5 training sites for LULC class which have been created based on the collected reference data from the field using GPS and Google Earth. The classification result in Figure 5B and Table 4 showed that forest land, cropland, dense shrubland, sparse shrubland, built-up area, and river bed are the major LULC types in the watershed. Cropland shares the greatest area coverage with 7019.29 ha (55.66%) from all other types of LULC types and followed by dense shrubland which covers about 3906.20 hectares (31.05%) of the study watershed.

Figure 5.

Reference sample points for image classification and accuracy assessment (A), LULC map of the watershed (B), slope (percent) map of the watershed (C), and support and conservation practice (P) factor map (D).

Table 4.

Watershed LULC types and area coverage.

LULC types	Description	Area
LULC types	Description	ha	Percent
Dense shrubland	Land area dominantly occupied by densely grown short trees	3906.2	31.1
Sparse shrubland	The land area with very sparsely grown short trees and scrubs	1078.7	8.6
Forest area	An area dominantly occupied by dense natural and plantation forest, and riparian forest	375.2	3.0
Cropland	Land areas used for subsistence rainfed agriculture and irrigation	7003.3	55.7
Built-up area	Permanent or temporary houses/buildings, villages	55.0	0.4
River bed	The channel bottom of a river or the ground a river flow over	163.9	1.3

According to Haque and Basak,⁵⁴ satellite image classification accuracy should be done by ground-truthing or by the physical appearance in the study area. In this study, a total of 116 reference points collected from the field and Google Earth were used as reference data (User-defined Points) to check the overall image classification accuracy and develop an error matrix. An error matrix (confusion matrix) results in Table 5 below showed that the classification had an overall classification accuracy of 87.07% with overall Kappa Statistics of 0.83, which indicates a good agreement exists between the reference data and the classified image. According to Foody,⁵⁵ classifications of remotely sensed data with overall classification accuracy >85% has been widely used and accepted for further analysis.

Table 5.

Error matrix (confusion matrix) of the classified image (February 1, 2019).

Classified data	Reference data								Conditional Kappa for each category
Classified data	Dense shrubland	River bed	Built-up area	Sparse shrubland	Forest area	Crop land	Total	User’s accuracy (%)	Conditional Kappa for each category
Dense shrubland	13	0	0	0	1	1	15	86.67	0.8484
River bed	1	9	0	0	0	0	10	90.00	0.8895
Built-up area	0	2	11	0	0	0	13	84.62	0.8250
Sparse shrubland	0	0	3	11	0	0	14	78.57	0.7587
Forest area	0	0	0	2	13	0	15	86.67	0.8405
Crop land	0	0	0	0	5	44	49	89.80	0.8333
Total	14	11	14	13	19	45	116	86.67	0.8484
Producer’s accuracy (%)	92.86	81.82	78.57	84.62	68.42	87.78
Overall classification accuracy = 87.07%
Overall Kappa statistics = 0.8308

The P-factor map was generated by integrating the classified LULC map (Figure 5B) with the classified slope map (Figure 5C) based on Wischmeier and Smith²⁷ method (Table 6). The P-factor value of the watershed varies from 0.1 in croplands to 1 in nonagricultural land uses (Figure 5D). The P-factor map shows that the P value is high in forest and dense shrubland areas. This is due to the fact that forests and dense shrubs have a high contribution in reducing the quantity and rate of runoff and the amount of soil erosion.

Table 6.

Support and management practice (P) value adopted from Wischmeier and Smith (1978).

Land use	Slope (%)	Area		P factor
Land use	Slope (%)	Ha	Percent	P factor
Agricultural land use	0-5	550.5	4.4	0.10
	5-10	1212.6	9.6	0.12
	10-20	2429.6	19.3	0.14
	20-30	1317.1	10.5	0.19
	30-50	1050.7	8.4	0.25
	50-100	442.3	3.5	0.33
Other land uses	0-100	5579.4	44.3	1.00

Potential soil loss of Chereti watershed

The spatial data layers of 5 factors of RUSLE were combined using a raster calculator of Spatial Analyst extension in ArcGIS 10.4 software, and Zonal Statistics as table in Spatial Analyst Extension of ArcGIS 10.4 software was used to calculate mean soil loss of the watershed. The result showed that the annual soil erosion rate of the watershed ranged from 0 t ha⁻¹ year⁻¹ in flat areas to 187.47 t ha⁻¹ year⁻¹ in upper undulating areas of the watershed (Figure 6) with a mean annual soil loss of 38.7 t ha⁻¹ year⁻¹. The result indicated that the entire watershed lost a total of about 487 057.7 tons of soil annually from an area of 125.82 km².

Figure 6.

Soil erosion map of Chereti watershed.

The soil loss severity class adopted by Haregeweyn et al,¹⁰ Belayneh et al,³ Zerihun et al,³⁰ and Yesuph and Dagnew,¹² for Ethiopian highlands was used to reclassify and develop soil loss severity map from the soil loss map of the watershed. The result showed that about 56.3% (7080.0 ha) of the watershed experienced a very slight (0-5 t ha⁻¹ year⁻¹) rate of soil erosion hazard which contributed about 2.0% to the annual soil loss of the whole watershed, and about 26.9% (3379.4 ha) of the watershed is classified under moderate rate of soil erosion hazard which contributed about 76.5% to the annual soil loss of the whole watershed. The remaining 12.9% (1619.6 ha), 3.9% (492.4 ha), and 0.1% (10.8 ha) are characterized by slight, severe, and very severe rates of soil erosion, respectively (Table 7 and Figure 7A). A very slight risk of soil erosion was found in gentle slope areas, and severe and very severe risk of soil erosion was found in steep slope areas of the watershed. During the field visit, high soil erosion was observed in steep and very steep slope areas with a lack of soil conservation practices.

Table 7.

Watershed soil loss severity levels and area coverage.

Soil loss (t/ha/year)	Severity levels	Area		Annual soil loss (ton/year)	Contribution to annual watershed soil loss (%)
Soil loss (t/ha/year)	Severity levels	Ha	Percent	Annual soil loss (ton/year)	Contribution to annual watershed soil loss (%)
0-5	Very slight	7080.0	56.3	9510.7	2.0
5-15	Slight	1619.6	12.9	34 726.6	7.1
15-30	Moderate	3379.4	26.9	372 820.9	76.5
30-50	Severe	492.4	3.9	67 558.5	13.9
>50	Very severe	10.8	0.1	2480.6	0.5

Figure 7.

Soil loss severity class map (A), and soil erosion susceptibility map of subwatersheds in the study area (B).

Subwatershed’s erosion vulnerability mapping and prioritization

Soil erosion vulnerability area identification and prioritization for conservation planning is imperative because it is difficult to implement conservation measures at the same time in all parts of the watershed. A total of 18 subwatersheds (Figure 7B and Table 8) were delineated for soil conservation prioritization.

Table 8.

Mean and annual soil loss of sub-watersheds and their priority level.

Subwatersheds	Area		Mean soil loss (ton/ha/year)	Annual soil loss (ton/year)	Contribution to total annual soil loss (%)	Priority level
Subwatersheds	ha	Percent	Mean soil loss (ton/ha/year)	Annual soil loss (ton/year)	Contribution to total annual soil loss (%)	Priority level
SWS1	647.0	5.1	13.4	22 379.6	4.6	7
SWS2	518.9	4.1	6.2	8364.4	1.7	18
SWS3	176.0	1.4	7.3	3296.3	0.7	16
SWS4	428.9	3.4	14.0	15 522.1	3.2	6
SWS5	595.1	4.7	16.1	24 675.1	5.1	5
SWS6	185.5	1.5	13.0	6230.7	1.3	11
SWS7	501.7	4.0	8.0	10 367.4	2.1	15
SWS8	883.4	7.0	10.3	23 597.6	4.8	14
SWS9	369.6	2.9	12.3	11 742.9	2.4	12
SWS10	470.4	3.7	17.3	20 960.0	4.3	3
SWS11	675.1	5.4	13.2	23 036.8	4.7	10
SWS12	784.6	6.2	11.8	23 855.6	4.9	13
SWS13	543.0	4.3	6.3	8850.5	1.8	17
SWS14	890.0	7.1	13.3	30 601.4	6.3	8
SWS15	943.2	7.5	18.1	44 102.6	9.1	2
SWS16	505.1	4.0	13.3	17 289.6	3.6	9
SWS17	1,549.1	12.3	17.1	68 511.7	14.1	4
SWS18	1,915.8	15.2	25.0	123 776.5	25.4	1

The subwatershed soil erosion vulnerability map in Figure 7B and Table 8 and showed that the mean annual soil erosion rate of subwatersheds varied from 6.2 t ha⁻¹ year⁻¹ generated from SWS2 subwatershed to 25.0 t ha⁻¹ year⁻¹ generated from SWS18 subwatershed.

Among the 18 delineated subwatersheds, SWS18, SWS15, SWS10, SWS17, and SWS5 subwatersheds were identified as subwatersheds with a high rate mean annual rate of soil loss (Table 8) and were given a rank of priority level of first, second, third, forth, and fifth, respectively, for conservation planning. These subwatersheds cover a total area of about 5473.6 ha (43.5%) from the whole watershed and contributes about 57.9% to the annual watershed soil loss.

Soil erosion across LULC type, soil type, and slope

The spatial distribution soil loss in the watershed varied across LULC types, slope, and soil type. The result in Table 9 and Figure 8B showed that the highest mean soil loss rate was generated from sparse shrubland followed by river beds. The likely cause for the highest mean soil loss in sparse shrubland is that majority of the sparse shrubland area is dominated by steep slope areas with very sparse shrubs (nearly bare land). However, about 73.9% of total annual soil loss of the watershed was generated from croplands with a mean soil loss rate of 19.9 t ha⁻¹ year⁻¹. The smallest mean annual soil loss was generated from the forest and dense shrubland area due to the role of vegetation cover in reducing the energy of raindrops and the amount of runoff.

Table 9.

Soil loss across LULC type.

LULC type	Area		Mean soil loss (ton/ha/year)	Annual soil loss (ton/year)	Contribution to total annual soil loss (%)
LULC type	ha	Percent	Mean soil loss (ton/ha/year)	Annual soil loss (ton/year)	Contribution to total annual soil loss (%)
Dense shrubland	3906.2	31.1	3.3	33 199.2	6.8
Sparse shrubland	1078.7	8.6	27.8	77 494.4	15.9
Forest area	375.2	3.0	2.9	2781.9	0.6
Cropland	7003.3	55.7	19.9	359 751.6	73.9
Built-up area	55.0	0.4	21.5	3051.7	0.6
River bed	163.9	1.3	25.5	10 787.7	2.2

Figure 8.

Mean soil loss of sub-watersheds (A), and soil loss across LULC type, slope, and soil type (B).

The topography is one of the factors which significantly affect the spatial distribution of soil loss in the watershed. The result indicated that areas with steep slope topography have greater soil loss rates than areas with gentle slopes. As shown in Table 10 and Figure 8B, the soil loss rates increase as the slope gradient increases across the watershed. About 88.7% of total annual soil loss of the watershed was generated from areas characterized by moderately steep to very steep slope topography. Hence, special emphasis should be given to implementing soil and water conservation measures in steep slope areas of the watershed.

Table 10.

Soil loss across different slope degrees.

Slope description	Slope class (%)	Area		Mean soil loss (ton/ha/year)	Annual soil loss (ton/year)	Contribution to total annual soil loss (%)
Slope description	Slope class (%)	ha	Percent (%)	Mean soil loss (ton/ha/year)	Annual soil loss (ton/year)	Contribution to total annual soil loss (%)
Flat to very gentle sloping	0-2	113.4	0.9	1.2	360.2	0.1
Gently sloping	2-5	567.1	4.5	2.5	3691.5	0.8
Sloping	5-10	1567.2	12.5	5.4	21 983.3	4.5
Strongly sloping	10-15	1588.4	12.7	7.1	29 134.1	6.0
Moderately steep	15-30	3695.0	29.4	12.8	122 560.9	25.2
Steep	30-60	3423.1	27.3	17.1	151 304.0	31.1
Very steep	> 60	1604.5	12.8	38.1	158 003.6	32.4

As shown in Table 11 and Figure 8B, the mean soil loss values also varied across soil types in the watershed. Despite similar erodibility value with Eutric Cambisols soil type, Lithosols and Eutric Regosols soil types were found to be more vulnerable with a mean soil loss rate of 26.1 and 21.7 t ha⁻¹ year⁻¹, respectively, because the soil types are found in steep slope areas of the watershed. Eutric Cambisols has the lowest mean soil loss (7.9 t ha⁻¹ year⁻¹) mainly because most of its area is dominated by a relatively lower slope gradient. This result indicated that the topographic (LS) factor is the primary influential RUSLE parameter to estimate soil loss in the watershed.

Table 11.

Soil loss across soil types.

Soil type	Area		Mean soil loss (ton/ha/year)	Annual soil loss (ton/year)	Contribution to total annual soil loss (%)
Soil type	ha	Percent	Mean soil loss (ton/ha/year)	Annual soil loss (ton/year)	Contribution to total annual soil loss (%)
Eutric Cambisols	6242.7	49.6	7.9	127 224.0	26.1
Eutric Regosols	2850.2	22.7	21.7	159 555.2	32.8
Vertic Cambisols	1019.2	8.1	12.8	33 674.0	6.9
Lithosols	2470.2	19.6	26.1	166 614.3	34.2

Discussion

In this study, multiple data such as 2 Landsat 8 OLI imageries, DEM data, soil data, and rainfall data were integrated with the RUSLE model to estimate potential soil loss of Chereti Watershed. In this study, validation of the model output was challenging due to poorly available input data which is the common problem in developing countries. However, as an alternative, the consistency of the model output of this study was compared with the outputs of previous studies with similar approaches and objectives in Ethiopian highlands. Similarly, previous studies^{3,9,10,12,36,52} in Ethiopia have adopted comparing their result of soil loss with the output of previous studies employing similar approaches to evaluate their RUSLE model results.

The result indicated that the annual soil erosion rate of the watershed ranged from 0 to 187.47 t ha⁻¹ year⁻¹ with a mean annual soil loss of 38.71 t ha⁻¹ year⁻¹. The estimated mean annual soil loss result of this study agreed with the finding of other similar and recent studies made by Yesuf and Dagnew,¹² which estimated the mean annual soil loss rate of 37 t ha⁻¹ year⁻¹for the Beshillo watershed in the Blue Nile Basin of Ethiopia. Alemu and Melesse³⁷ were estimated a mean annual soil loss of 37 t ha⁻¹ year⁻¹ for the Anjeni watershed in Northwestern highlands Ethiopia. Belayneh et al³ also found the mean annual soil loss of 42.67 t ha⁻¹ year⁻¹ for Gumara watershed in the northwestern highlands of Ethiopia which is relatively agreed with the result of this study. The estimated mean annual soil loss of this study was found to be relatively consistent with the overall mean soil loss of the country which is about 42 t ha⁻¹ year⁻¹.^36,56,57A study made by Wolka et al⁴⁶ in Chaleleka wetland watershed of the central rift valley of Ethiopia estimated a mean annual soil loss rate of 45 t ha⁻¹ year⁻¹. Another study made by Zerihun et al³⁰ found a mean annual soil loss of 49 t ha⁻¹ year⁻¹for Demcha district in Northwestern Ethiopia. Gelagay and Minale³⁶ also reported a mean annual soil loss rate of 47 t ha⁻¹ year⁻¹ in the Koga watershed of Northwestern Ethiopia. Woldemariam et al,²⁹ estimated a mean erosion rate of 51.04 and 34.26 t ha⁻¹ year⁻¹ in 2000 and 2016, respectively, for Gobele Watershed, East Hararghe, Ethiopia

On the contrary, other similar studies reported higher and lower values of soil erosion rate than the estimated value of this study. Tessema et al³⁹ reported a mean soil loss of 31 t ha⁻¹ year⁻¹ for Welmel Watershed in Genale Dawa Basin of Ethiopia. The mean annual soil loss of this study (38.71 t ha⁻¹ year⁻¹) was found to be 2 times beyond the maximum soil loss tolerance value of 18 t ha⁻¹ year⁻¹ for Ethiopian highlands as determined by Hurni,³⁴ Gashaw et al,⁹ reported a mean soil loss rate of 23.7 t ha⁻¹ year⁻¹ in Gelada watershed which is much lower than the estimated soil loss of this study. The estimated soil loss of this study was also found to be much higher than 26 t ha⁻¹ year⁻¹ estimated by Sisay,⁵⁸ for the Wondo Genet watershed. A study conducted in Chemoga watershed in the northwestern highlands of Ethiopia by Bewket and Teferi,⁵⁹ estimated an average soil erosion rate of 93 t ha⁻¹ year⁻¹ which is by far greater than the result of this study.

In this study, it was found that about 96% of the watershed area experienced very slight (56.27%), slight (12.87%), and moderate (26.86%) risk of soil erosion, and the remaining 4% of the watershed area dominated by steep and very steep slopes experienced severe and very severe risk of soil erosion. Similarly, Yesuph and Dagnew¹² found that about 85.2% of Beshillo Catchment of the Blue Nile Basin, Ethiopia, experienced very slight (42.9%), slight (31%), and moderate (11.3%) risk of soil erosion. Zerihun et al³⁰ also reported 84% of Dembecha district, Northwestern highlands of Ethiopia, experienced very slight (24%), slight (49%), and moderate (11%) risk of soil erosion.

The result of this study showed that maximum soil loss was estimated in sparse shrubland and river valleys; and low mean soil loss was estimated in the forest and dense shrubland areas. Similarly, Alemu and Melesse³⁷ reported maximum mean annual soil loss in bushlands and river valleys than other land use types in Anjeni watershed, Northwest Ethiopia. Our finding is in agreement with the study made by Belayneh et al,³ for Gumara watershed which reported a low level of mean annual soil loss in forest and shrubland areas. The finding of this study also showed that mean soil loss increases as slope gradient increases. Consistent to this finding, Yesuph and Dagnew¹² reported the spatial distribution of mean soil loss rate increases with an increasing slope gradient in Beshillo Catchment of the Blue Nile Basin, Ethiopia. Belayneh et al³ also reported that mean soil loss increases as the slope gradient increases in Gumara watershed of Northwestern Ethiopia.

What makes our method different from previous researches in Ethiopian highlands?

Several previous studies^{3,9,12,30,36,52} in Ethiopian highlands have used satellite image acquired during the winter season to estimate crop and management (C) factor value of RUSLE model and adopted manually assigning C-factor value from the values suggested and used in different studies after the classification of the satellite image. However, satellite image data acquired during the winter season could not accurately reflect the effect of crop and management (C) factor in the RUSLE model because crop cover and management practices are high during the summer season in Ethiopian highlands. Assigning C-factor value from previously published studies for the LULC type may not also accurately reflect C-factor because cover density varies within each LULC type. In this study, we have used Landsat image acquired during the summer (July) season for C-factor estimation because vegetation cover, crop cover, and management practices are high during this season. We have also employed NDVI to estimate the C-factor value of the RUSLE model, instead of using manually assigning C-factor value from previously published C-values which is the common method to determine C-factor in previous studies in Ethiopian highlands. Previous studies^{3,9,12,29,30,36,53} in Ethiopian highlands also used a single-date Landsat image to determine both the C-factor and P-factor value of the RUSLE model. But, in this study, 2 Landsat images acquired during summer (July) and winter (February) seasons were used to determine the C-factor and P-factor value of RUSLE, respectively.

Limitations of the study

There are various methods to estimate soil erodibility (K) factor based on either field data or secondary data. In this study, we have computed the soil erodibility (K) factor of RUSLE based on soil type and soil color method using soil type map as adopted by several studies^{3,30,36,39,40} in Ethiopian highlands. However, this method may not accurately reflect the effect of soil properties such as soil texture, organic matter, soil structure, and permeability on soil erodibility. Soil color may not also reflect soil erodibility because the color of the soil can be affected by parent material, CaCO₃, and organic matter.

In this study, establishing erosion monitoring sites or representative experimental plots and collecting soil erosion data from the sites or plots was not conducted to validate the estimated soil loss using measured or actual soil loss. Furthermore, a sensitivity analysis was not conducted to identify the most determinant and sensitive RUSLE parameter.

Conclusion and Final Remarks

This study was carried out to estimate and map soil erosion in Chereti watershed using the integration of RUSLE, GIS, and remote sensing techniques. Five RUSLE factors such as rainfall erosivity (R), soil erodibility (K), topographic (LS), crop and cover management (C), and conservation practice (P) factors were computed, mapped, and overlayed together to estimate the potential soil loss of the watershed. In this study, Landsat 8 OLI satellite image acquired during summer (July) season, and Normalized Vegetation Difference Vegetation Index (NDVI) method was used to determine the C-factor value of the RUSLE model.

The result showed that the annual soil loss rate of the watershed ranged from 0 to 187.47 t ha⁻¹ year⁻¹ with a mean annual soil loss of 38.7 t ha⁻¹ year⁻¹. The entire watershed lost a total of about 487 057.7 tons of soil annually. About 95% of the watershed area experienced very slight, slight, and moderate risk of soil erosion, and about 4% of the watershed area is characterized by severe and very severe risk of soil erosion. Soil loss in the watershed varied across slope gradient, LULC type, and soil type. Slight to a very slight risk of soil erosion was found in flat and gentle slope areas, and severe and very severe risk of soil erosion dominated in steep and very slope areas of the watershed. Soil types in steep slope areas of the watershed were found to be more susceptible to soil loss despite similar soil erodibility (K) value with soil types in gentle slope areas. Mean soil loss was also found to be high in river valleys and sparse shrubland areas dominated by steep slopes, and low mean soil loss was estimated in the forest and dense shrubland areas.

From the total 18 subwatersheds in the study area, about 57.9% of the total annual soil loss of the entire watershed was generated from 5 subwatersheds such as SWS18, SWS15, SWS10, SWS17, and SWS5 which need prior intervention for the planning and implementation of soil conservation measures.

Further related research works should pay more attention to the investigation of soil erodibility (K) of the watershed based on field data and other methods that could consider the effect of soil properties such as soil texture, organic matter, soil structure, and permeability on soil erosion.

In this study, the validity of the RUSLE model was evaluated by comparing the model output with the outputs of previous studies with similar approaches in Ethiopian highlands because it was challenging due to poorly available numerical data. Hence, establishing erosion monitoring sites or representative experimental plots and collecting soil erosion data should be conducted to validate the predictive capacity of the RUSLE model in the watershed.

Revised Universal Soil Loss Equation does not consider gully erosion which is now the major problem that seriously threatens and dissects agricultural lands in Chereti watershed. Hence, further study on gully erosion estimation should be conducted in the watershed.

The integrated use of RUSLE with GIS and remote sensing techniques was found to be indispensable, less costly, and effective for the estimation and mapping of soil erosion potential, and identification and prioritization vulnerable subwatersheds for conservation planning and implementation.

Footnotes

Acknowledgements

The authors would like to thank US Geologic Survey for enabling us to download satellite images and the DEM of the studied watershed freely from their website. They would like to acknowledge the National Meteorology Agency of Ethiopia, and the Ethiopian Ministry of Water, Irrigation, and Energy for providing us climate data and digital soil data of the watershed, respectively.

Declaration of conflicting interests:

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding:

The author(s) received no financial support for the research, authorship, and/or publication of this article.

Author Contributions

All authors contributed equally in designing the study, data collection and analysis, interpretation of the results, and writing and revising the manuscript. All authors read and approved the final manuscript.

ORCID iD

Ajanaw Negese

References

El Jazouli

Barakat

Ghafiri

El Moutaki

Ettaqy

Khellouk

. Soil erosion modeled with USLE, GIS, and remote sensing: a case study of Ikkour watershed in Middle Atlas (Morocco). Geoscience Letters. 2017;4:25. doi:10.1186/s40562-017-0091-6.

Ayalew

Deumlich

Šarapatka

Doktor

Quantifying the sensitivity of NDVI-based C factor estimation and potential soil erosion prediction using Spaceborne earth observation data. Remote Sens. 2020;12:1136. doi:10.3390/rs12071136.

Belayneh

Yirgu

Tsegaye

. Potential soil erosion estimation and area prioritization for better conservation planning in Gumara watershed using RUSLE and GIS techniques.’ Environ Syst Res. 2019;8:20. doi:10.1186/s40068-019-0149-x.

Ettazarini

El Jakani

Najoui

Assessment of soil loss risk using integrated remote sensing and geographic information system (GIS) techniques in the Argana Basin, Morocco. Am J Innovat Res Appl Sci. 2017;4:186-194.

Conway

Yost

Kitchen

Sudduth

Veum

KS.

Cropping system, landscape position, and topsoil depth affect soil fertility and nutrient buffering. Soil Sci Soc Am J. 2018;82:382-391.

Keesstra

Mol

de Leeuw

, et al. Soil-related sustainable development goals: four concepts to make land degradation neutrality and restoration work. Land. 2018;7:133.

Lal

Horn

Kosaki

Soil and Sustainable Development Goals; 2018. https://knowledge.unccd.int/sites/default/files/inline-files/Soil_and_Sustainable_Development_Goals_GeoEcology_Essays_ISBN_978_3_510_65425_3_274101900_Ed_Rattan_Lal_Rainer_Horn_Takashi_Kosaki_So.pdf.

Rodríguez-Blanco

Arias

Taboada-Castro

Nunes

Keizer

Taboada-Castro

MT.

Sediment yield at catchment scale using the SWAT (Soil and Water Assessment Tool) Model. Soil Science. 2016;181:326-334.

Gashaw

Tulu

Argaw

Erosion risk assessment for prioritization of conservation measures in Geleda watershed, Blue Nile basin, Ethiopia. Environ Syst Res. 2017;6:1.

10.

Haregeweyn

Tsunekawa

Poesen

, et al. Comprehensive assessment of soil erosion risk for better land use planning in river basins: case study of the Upper Blue Nile River. Sci Total Environ. 2017;574:95-108.

11.

Molla

Sisheber

Estimating soil erosion risk and evaluating erosion control measures for soil conservation planning at Koga watershed in the highlands of Ethiopia. Solid Earth. 2017;8:13-25.

12.

Yesuph

Dagnew

AB.

Soil erosion mapping and severity analysis based on RUSLE model and local perception in the Beshillo Catchment of the Blue Nile Basin, Ethiopia. Environ Syst Res. 2019;8:17. doi:10.1186/s40068-019-0145-1.

13.

Ayalew

Selassie

Soil loss estimation for soil conservation planning using geographic information system in Guang Watershed, Blue Nile Basin. J Environ Earth Sci. 2015;5:126-134.

14.

Pham

Degener

Kappas

Integrated universal soil loss equation (USLE) and Geographical Information System (GIS) for soil erosion estimation in A Sap basin: central Vietnam. Int Soil Water Conserv Res. 2018;6:99-110. doi:10.1016/j.iswcr.2018.01.001.

15.

Puno

Marin

Puno

RCC

Toledo-Bruno

AG.

Geographic information system and process-based modeling of soil erosion and sediment yield in agricultural watershed. Glob J Environ Sci Manage. 2020;7:1-14. doi:10.22034/gjesm.2021.01.01.

16.

Panagos

Meusburger

Ballabio

Borrelli

Alewell

Soil erodibility in Europe: a high-resolution dataset based on LUCAS. Sci Total Environ. 2014;479-480:189-200. doi:10.1016/j.scitotenv.2014.02.010.

17.

Alexakis

Tapoglou

Vozinaki

A-EK

Tsanis

IK.

Integrated use of satellite remote sensing, artificial neural networks, field spectroscopy, and GIS in estimating crucial soil parameters in terms of soil erosion. Remote Sens. 2019;11:1106. doi:10.3390/rs11091106.

18.

Shi

Castaldi

van Wesemael

Van Oost

Vis-NIR spectroscopic assessment of soil aggregate stability and aggregate size distribution in the Belgian Loam Belt. Geoderma. 2020;357:113958. doi:10.1016/j.geoderma.2019.113958.

19.

Aga

Chane

Melesse

Soil erosion modelling and risk assessment in data scarce rift valley lake regions, Ethiopia. Water. 2018;10:1684. doi:10.3390/w10111684.

20.

Diwakar

Kaur

Patel

Hydrologic assessment in a middle Narmada Basin, India using SWAT model. Int J Engineer Techn Comput Res. 2014;2:10-25.

21.

Mosbahi

Benabdallah

Boussema

MR.

Assessment of soil erosion risk using SWAT model. Arab J Geosci. 2012;6:4011-4019. doi:10.1007/s12517-012-0658-7.

22.

Bera

Estimation of soil loss by USLE model using GIS and remote sensing techniques: a case study of Muhuri River Basin, Tripura, India. Eura J Soil Sci. 2017;6:206-215. doi:10.18393/ejss.

23.

Prasannakumar

Vijith

Abinod

Geetha

Estimation of soil erosion risk within a small mountainous sub-watershed in Kerala, India, using Revised Universal Soil Loss Equation (RUSLE) and geo-information technology. Geosci Front. 2012;3:209-215. doi:10.1016/j.gsf.2011.11.003.

24.

Renard

Foster

Weesies

McCool

Yoder

DC.

Predicting Soil Erosion by Water: A Guide to Conservation Planning with the Revised Universal Soil Loss Equation (RUSLE) (Agriculture Handbook Number 703). Washington, DC: United States Department of Agriculture; 1997.

25.

Saha

GIS based soil erosion estimation using RUSLE model: a case study of upper Kangsabati Watershed, West Bengal, India. Int J Environ Sci Nat Res. 2018;13:555871. doi:10.19080/ijesnr.2018.13.555871.

26.

Swarnkar

Malini

Tripathi

Sinha

Assessment of uncertainties in soil erosion and sediment yield estimates at ungauged basins: an application to the Garra River basin, India. Hydrol Earth Syst Sci. 2018;22:2471-2485. doi:10.5194/hess-22-2471-2018.

27.

Wischmeier

Smith

DD.

Predicting Rainfall Erosion Losses: A Guide to Conservation Planning (Agriculture Handbook no. 282). Washington, DC: US Department of Agriculture; 1978.

28.

Angima

Stott

O’Neill

Ong

Weesies

GA.

Soil erosion prediction using RUSLE for central Kenyan highland conditions. Agricul Ecosyst Environ. 2003;97:295-308. doi:10.1016/s0167-8809(03)00011-2.

29.

Woldemariam

Lguala

Tekalign

Reddy

RU.

Spatial modeling of soil erosion risk and its implication for conservation planning: the case of the Gobele Watershed, East Hararghe Zone, Ethiopia. Land. 2018;7:25. doi:10.3390/land7010025.

30.

Zerihun

Mohammedyasin

Sewnet

Adem

Lakew

Assessment of soil erosion using RUSLE, GIS and remote sensing in NW Ethiopia. Geod Reg. 2018;12:83-90. doi:10.1016/j.geodrs.2018.01.002.

31.

IUSS-WRB. World Reference Base for Soil Resources 2014. International Soil Classification System for Naming Soils and Creating Legends for Soil Maps (World Soil Resources Reports No. 106). Rome, Italy: FAO; 2015.

32.

Koirala

Thakuri

Joshi

Chauhan

Estimation of soil erosion in Nepal using a RUSLE modeling and geospatial tool. Geosciences. 2019;9:147. doi:10.3390/geosciences9040147.

33.

Yin

Nearing

Borrelli

Xue

Rainfall Erosivity: an overview of methodologies and applications. Vadose Zone J. 2017;16:1-16. doi:10.2136/vzj2017.06.0131.

34.

Hurni

Erosion—productivity—conservation systems in Ethiopia. IV Int Conf Soil Conser. 1985;1985:654-674.

35.

Tesfaye

Debebe

Fikirie

Soil erosion risk assessment using GIS based USLE model for soil and water conservation planning in Somodo Watershed, South West Ethiopia. Int J Environ Agricul Res. 2018;4:35-43.

36.

Gelagay

Minale

AS.

Soil loss estimation using GIS and Remote sensing techniques: a case of Koga watershed, Northwestern Ethiopia. Int Soil Water Conserv Res. 2016;4:126-136. doi:10.1016/j.iswcr.2016.01.002.

37.

Alemu

Melesse

AM.

Impacts of longterm conservation measures on ecosystem services in Northwest Ethiopia. Int Soil Water Conserv Res. 2020;8:47-55. doi:10.1016/j.iswcr.2019.10.002.

38.

Mohammed

Alsafadi

Talukdar

, et al. Estimation of soil erosion risk in southern part of Syria by using RUSLE integrating geo informatics approach. Remote Sens Appl Soc Environ. 2020;20:100375. doi:10.1016/j.rsase.2020.100375.

39.

Tessema

Jasińska

Yadeta

Świtoniak

Puchałka

Gebregeorgis

EG.

Soil loss estimation for conservation planning in the Welmel Watershed of the Genale Dawa Basin, Ethiopia. Agronomy. 2020;10:777. doi:10.3390/agronomy10060777.

40.

Desalegn

Tezera

Tesfay

Developing GIS-based soil erosion map using RUSLE of Andit Tid Watershed, central highlands of Ethiopia. Journal of Scientific Research Reports. 2018;19:1-13. doi:10.9734/jsrr/2018/40841.

41.

Liu

Tang

Jiang

Zhu

Yang

Song

Regional-scale calculation of the LS factor using parallel processing. Comput Geosci. 2015;78:110-122. doi:10.1016/j.cageo.2015.02.001.

42.

Moore

Burch

GJ.

Physical basis of the length slope factor in the universal soil loss equation. Soil Sci Soc Am. 1986;50:1294-1298.

43.

Belasri

Lakhouili

Estimation of soil erosion risk using the universal soil loss equation (USLE) and geo-information technology in Oued El Makhazine Watershed, Morocco. J Geograph Inform Syst. 2016;8:98-107. doi:10.4236/jgis.2016.81010.

44.

Almagro

Thomé

Colman

, et al. Improving cover and management factor (C-factor) estimation using remote sensing approaches for tropical regions. Int Soil Water Conserv Res. 2019;7:325-334. doi:10.1016/j.iswcr.2019.08.005.

45.

Singh

Panda

RK.

Grid-cell based assessment of soil erosion potential for identification of critical erosion prone areas using USLE, GIS and remote sensing: a case study in the Kapgari watershed, India. Int Soil Water Conserv Res. 2017;5:202-211. doi:10.1016/j.iswcr.2017.05.006.

46.

Wolka

Tadesse

Garedew

Yimer

Soil erosion risk assessment in the Chaleleka wetland watershed, Central Rift Valley of Ethiopia. Environ Syst Res. 2015;4:5. doi:10.1186/s40068-015-0030-5.

47.

Rouse

Jr Haas

Well

Deering

DW.

Monitoring vegetation systems in the great plains with ERTS. ERTS Symp.1973;1:309-317.

48.

Patil

Sharmaand

SK.

Remote Sensing and GIS based modeling of crop/cover management factor (C) of USLE in Shakker river watershed. In: International Conference on Chemical, Agricultural and Medical Sciences (CAMS-2013), December 29-30, 2013, Kuala Lumpur, Malaysia; 2013.

49.

Durigon

Carvalho

Antunes

MAH

Oliveira

PTS

Fernandes

MM.

NDVI time series for monitoring RUSLE cover management factor in a tropical watershed. Int J Remote Sens. 2014;35:441-453. doi:10.1080/01431161.2013.871081.

50.

de Carvalho

Durigon

Antunes

MAH

de Almeida

de Oliveira

PTS.

Predicting soil erosion using Rusle and NDVI time series from TM Landsat 5. Pesquisa Agropecuária Brasileira. 2014;49:215-224. doi:10.1590/s0100-204x2014000300008.

51.

Phinzi

Ngetar

NS.

The assessment of water-borne erosion at catchment level using GIS-based RUSLE and remote sensing: a review. Int Soil Water Conserv Res. 2019;7:27-46. doi:10.1016/j.iswcr.2018.12.002.

52.

Miheretu

Yimer

AA.

Estimating soil loss for sustainable land management planning at the Gelana sub-watershed, northern highlands of Ethiopia. Int J River Basin Manage. 2017;16:41-50. doi:10.1080/15715124.2017.1351978.

53.

Tilahun

Taddesse

Melese

Mebrate

Assessment of spatial soil erosion hazard in Ajema Watershed, North Shewa Zone, Ethiopia. Adv Plants Agricul Res. 2018;8:552-558. doi:10.15406/apar.2018.08.00384.

54.

Haque

Basak

Land cover change detection using GIS and remote sensing techniques: a spatio-temporal study on Tanguar Haor, Sunamganj, Bangladesh. Egypt J Remote Sens Space Sci. 2017;20:251-263. doi:10.1016/j.ejrs.2016.12.003.

55.

Foody

GM.

Harshness in image classification accuracy assessment. Int J Remote Sens. 2008;29:3137-3158. doi:10.1080/01431160701442120.

56.

Chimdessa

Quraishi

Kebede

Alamirew

Effect of land use land cover and climate change on river flow and soil loss in Didessa River Basin, South West Blue Nile, Ethiopia. Hydrology. 2018;6:2. doi:10.3390/hydrology6010002.

57.

Bewket

Abebe

Land-use and land-cover change and its environmental implications in a tropical highland watershed, Ethiopia. Int J Environ Stud. 2013;70:126-139. doi:10.1080/00207233.2012.755765.

58.

Sisay

Landscape scale soil erosion modeling and risk mapping of mountainous areas in eastern escarpment of Wondo Genet Watershed, Ethiopia. Int Res J Arts Soc Sci. 2014;6:1743624. doi:10.14303/irjas.2014.040.

59.

Bewket

Teferi

Assessment of soil erosion hazard and prioritization for treatment at the watershed level: case study in the Chemoga watershed, Blue Nile basin, Ethiopia. Land Degrad Develop. 2009;20:609-622. doi:10.1002/ldr.944.