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−1 year−1 in steep slope areas with a mean annual soil loss of 38.7 t ha−1 year−1, 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.
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.1 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.5 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.8
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 studies3,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,14 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 al15 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 al16 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 al17 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 (CaCO3), 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 al18 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 authors19–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).2 Angima et al28 employed the RUSLE model to predict and map soil erosion risk of Kianjuki Catchment in the central Kenya using GIS.
Different previous studies3,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 km2. 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.
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.31 Eutric Regosols, Lithosols (now are Leptosols),31 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.
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))24 given below
where, A is estimated annual soil loss (t ha−1 year−1); R is rainfall erosivity factor (MJ mm ha−1 h−1 year−1); K is soil erodibility factor (t ha−1 MJ−1 mm−1); 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 Smith27 and Yin et al,33 rainfall erosivity (R) factor is the product of kinetic energy of rainstorm and its maximum 30-min intensity (I30). However, due to the lack of rainfall kinetic energy and rainfall intensity (I30) data for the study area, the empirical equation that estimates R-value from total annual rainfall suggested by Hurni,34 for Ethiopian condition was adopted to compute the rainfall erosivity (R) factor of the watershed. This equation was also adopted by similar studies3,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 al33 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))34 given below in Raster Calculator of ArcGIS 10.4 Spatial Analyst Extension.
where R is the rainfall erosivity factor (MJ mm ha−1 h−1) 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.33
Table 2.
Meteorology stations and calculated mean annual rainfall (mm).
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.37 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 Hurni34 for the Ethiopian condition as adopted by several studies3,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,34 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,42 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:
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 researchers3,9,12,32,44,45 assigned published C values after classifying remotely sensed data into different LULC classes. Others authors1,4,23,25,43,46 employed NDVI, one of the vegetation indices, introduced by Rouse et al,47 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.48
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.
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,49 for tropical regions and adopted by de Carvalho et al,50 was employed to estimate the C-factor value of the watershed from the calculated NDVI. According to Phinzi and Ngetar51 and Almagro et al,44 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.
where CrA 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,3 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 Smith27 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 by3,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,27 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.
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−1 year−1(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−1 MJ−1 mm−1 in areas with Vertic Cambisols to 0.20 t ha−1 MJ−1 mm−1 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−1 MJ−1 mm−1 mainly dominated by croplands with Eutric Cambisols.
Table 3.
Soil types of the watershed and their respective erodibility (K) factor value.
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.
According to Haque and Basak,54 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,55 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).
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 Smith27 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).
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−1 year−1 in flat areas to 187.47 t ha−1 year−1 in upper undulating areas of the watershed (Figure 6) with a mean annual soil loss of 38.7 t ha−1 year−1. 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 km2.
The soil loss severity class adopted by Haregeweyn et al,10 Belayneh et al,3 Zerihun et al,30 and Yesuph and Dagnew,12 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−1 year−1) 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.
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.
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−1 year−1 generated from SWS2 subwatershed to 25.0 t ha−1 year−1 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−1 year−1. 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.
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.
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−1 year−1, 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−1 year−1) 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.
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 studies3,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−1 year−1 with a mean annual soil loss of 38.71 t ha−1 year−1. 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,12 which estimated the mean annual soil loss rate of 37 t ha−1 year−1for the Beshillo watershed in the Blue Nile Basin of Ethiopia. Alemu and Melesse37 were estimated a mean annual soil loss of 37 t ha−1 year−1 for the Anjeni watershed in Northwestern highlands Ethiopia. Belayneh et al3 also found the mean annual soil loss of 42.67 t ha−1 year−1 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−1 year−1.36,56,57A study made by Wolka et al46 in Chaleleka wetland watershed of the central rift valley of Ethiopia estimated a mean annual soil loss rate of 45 t ha−1 year−1. Another study made by Zerihun et al30 found a mean annual soil loss of 49 t ha−1 year−1for Demcha district in Northwestern Ethiopia. Gelagay and Minale36 also reported a mean annual soil loss rate of 47 t ha−1 year−1 in the Koga watershed of Northwestern Ethiopia. Woldemariam et al,29 estimated a mean erosion rate of 51.04 and 34.26 t ha−1 year−1 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 al39 reported a mean soil loss of 31 t ha−1 year−1 for Welmel Watershed in Genale Dawa Basin of Ethiopia. The mean annual soil loss of this study (38.71 t ha−1 year−1) was found to be 2 times beyond the maximum soil loss tolerance value of 18 t ha−1 year−1 for Ethiopian highlands as determined by Hurni,34 Gashaw et al,9 reported a mean soil loss rate of 23.7 t ha−1 year−1 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−1 year−1 estimated by Sisay,58 for the Wondo Genet watershed. A study conducted in Chemoga watershed in the northwestern highlands of Ethiopia by Bewket and Teferi,59 estimated an average soil erosion rate of 93 t ha−1 year−1 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 Dagnew12 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 al30 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 Melesse37 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,3 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 Dagnew12 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 al3 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 studies3,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 studies3,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 studies3,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, CaCO3, 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−1 year−1 with a mean annual soil loss of 38.7 t ha−1 year−1. 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.
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.
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.
