Assessment of soil erosion risks in response to land-use and land-cover changes in Coka watershed, Southern Ethiopia

ABSTRACT

Soil erosion is a serious problem arising from agricultural intensification and other anthropogenic activities which can cause deterioration of soil fertility and productivity. This study aims to estimate total amounts of soil loss using the RUSLE and NDVI models and identify major areas of concern for management interventions through site specific in the Coka watershed, Southern Ethiopia. The RUSLE model required the integration of thematic factors such as rainfall, length and steepness of the slope, land-use and land-cover, soil erodibility, and control practices. These factors were calculated using remote sensing data and GIS. The results indicated that total area coverage from study of watershed under the degree of severity of erosion and priority was <5(very low), 5 − 10 (low), 10 − 18 (moderate), 18 − 50 (high), 50 − 100 (very high) and >100 t ha^–1yr^–1 (severe). Soil loss rates ranged between 2.33 in bare land and 237.16 t ha^–1yr^–1 in forest. Based on the findings, erosion risk having estimated soil loss (≥ 18 t ha^–1yr^–1) which is above tolerable range for the country containing 73% (2723 ha) from the total area coverage of watershed. Thus, appropriate soil and water conservation practices should be adopted to restore degraded lands and environment by giving attention for erosion hotspot areas.

KEYWORDS:

• Environment
• erosion risk
• land-use and land-cover
• RUSLE
• soil productivity

1. Introduction

Soil erosion is a global major environmental problem leading to loss of nutrient-rich top surface soil, increased runoff from the more impermeable subsoil, and decreased water availability to plants (Fu et al., 2006; Ganasri & Ramesh, 2015; Pal et al., 2021). Modifications in land-use and management primarily drive soil erosion. The FAO-led Global Soil Partnership reports that 75 billion tons of soil are eroded every year from arable lands worldwide, which equates to an estimated financial loss of US$400 billion per year (GSP (Global Soil Partnership), 2017). The average soil erosion rates throughout the world are estimated between 12 and 15 t ha⁻¹ y⁻¹ (Biggelaar et al., 2003). It means that every year the land surface loses about 0.96–1.2 mm of soil (FAO, 2019), which is results in estimated financial loss of US$ 64 − 80 billion per year (FAO, 2017). In addition, soil erosion threatens sustainable agricultural production and water quality in tropical and semi-arid regions (Prasannakumar et al., 2012). Thus, soil erosion is a global economic, environmental, and ecological problem, mainly affecting natural resources and agricultural production (Chakrabortty et al., 2020a; Paveen & Kumar, 2012; Ighodaro et al., 2013; Pal et al., 2022).

Several studies have proven the effectiveness of remote sensing data in monitoring soil erosion, and land-use and land-cover (LULC) changes in Sub-Saharan Africa and elsewhere (Atoma et al., 2020; Gebreegziabher et al., 2021; Tadesse et al., 2017). Seutloali and Beckedahl (2015) reported that LULC change is the major causal factor driving soil loss in upland watersheds, particularly in Sub-Saharan Africa. It is the main approach through which human beings influence the ecosystem and contains plentiful information of human activities (Kishtawal et al., 2010; Pielke et al., 2011). As in the other parts of the world, the fastest change in LULC was due to increasing population pressure on forest resources, and agricultural land expansion to meet the graduated food demand, increasing demand for fuel wood and raw construction material, and the other anthropogenic influences (Gebreegziabher et al., 2021; Samuel et al., 2016; Uwemeye et al., 2020).

In Ethiopia, water erosion is the leading cause of degradation of the soil capital and the environment. It affects, with varying intensities, about 1.9 billion tons of topsoil from the highlands by soil erosion and losses an estimated USD$ 1 billion y⁻¹ of land in Ethiopia (Addis & Klik, 2015; Atoma et al., 2020; FAO, 1984; Sonneveld, 2002). This amount of soil could have added about 1.5 million tons of grain to the country’s annual harvest (Lulseged & Vlek, 2008). Soil erosion is a severe threat to food security and directly impacts livelihoods among rural communities in Ethiopia (Berhan et al., 2014; Samuel et al., 2016). There is wide range of estimates on the annual rate of soil erosion by water is the dominant degradation process and occurs particularly on cropland, with annual soil loss rates on average of 42 ton ha⁻¹ y⁻¹ for croplands, and up to 300 ton ha⁻¹ y⁻¹ in extreme cases (Tamrie, 1995; Tesfaye et al., 2013). Bewket and Teferi (2009) estimated annual soil loss at ton ha⁻¹ y⁻¹ for the entire Chemago watershed based on field assessment of rill and inter-rill erosion. Haile and Fetene (2012) estimated that about 97.04% of Kilie catchment, East Shoa, has 0–10 ton ha⁻¹ y⁻¹ erosion rate. Temporal patterns of LULC changes, especially in areas susceptible to accelerated soil erosion, provide further reason to reevaluate former qualitative approaches, considering the worldwide increase in croplands by 279 million hectares during 1985 − 2013 (Gebreegziabher et al., 2021). Therefore, in recent years, local authorities increasingly concerned about the issue of soil fertility have funded several research projects aimed at the conservation of water resources and limitation of soil erosion.

Some recent studies that have applied the RUSLE model to investigate the soil loss through physical field assessments and indicated that the high rates of soil erosion in Ethiopia are mainly caused by the combined effects of deforestation, overgrazing, detrimental cultivation practices with an emphasis on small-seed crops that require fine tillage, poverty, land fragmentation, expansion of cropland onto steep slopes, and unsustainable use of natural resources (Amsalu et al., 2007; Billi, 2015; Chakrabortty et al., 2022; Descheemaeker et al., 2006; Haile & Fetene, 2012; Samuel et al., 2016; Sonneveld & Keyzer, 2003). Erosion could also generate deposition of soil materials in the reservoirs, irrigation schemes, and waterways downstream (Cerda & Doerr, 2008). To curb such problems, assessment, planning, and implementation of locally acceptable and effective soil conservation practices are urgently needed (Sonneveld & Keyzer, 2003).

A wide range of models were developed to assess the rate of soil erosion. Model simulations are the basic methods for soil erosion prediction, usually, there are two types of prediction models such as physically-based models, and empirical models. Physically-based models theoretically have a greater transferability than empirical models, are more likely to achieve sound predictions (Schmidt et al., 1999; Wainwright & Mulligan, 2013). The prediction of rates of soil erosion typically involves the use of empirical models (Leh et al., 2013), thus the empirical models are still widely used for soil erosion assessment both on the watershed and regional scale. Of all the empirical models, the Universal Soil Loss Equation (USLE; Wischmeier & Smith, 1978) and the Revised Universal Soil Loss Equation (RUSLE) are used probably the most widely used empirical models for water erosion assessment and water conservation plans worldwide (Kinnel, 2010; Renard et al., 1997; Oliveira et al., 2011). The RUSLE model reveals that soil erosion is greatest on cultivated land than other land-use types (Gimenez-Morera et al., 2010; Hurni, 1993; Pal & Chakrabortty, 2019). It has been widely used in predicting soil erosion loss because of its simplicity, limited data requirement, and compatibility with Geographical Information Systems (GIS) (Angima et al., 2003; Bewket & Teferi, 2009; Lu et al., 2004; Hurni, 1985; Milward & Mersey, 1999). It has been recommended that where data for validation of model are scarce, with careful consideration of transfer functions and parameters, the model can be used for cost-effective conservation planning purposes (Angima et al., 2003; Chowdary et al., 2013; Haile & Fetene, 2012; Milward & Mersey, 1999; Garg & Jothiprakash, 2012; Prasannakumar et al., 2012; Wischmeier & Smith, 1978). Risse et al. (1993) reported that among six RUSLE factors topography factors (LS) and the C factor exert the greatest effect on the global model efficiency. The Normalized Difference Vegetation Index (NDVI) values have been developed to estimate the RUSLE C factor (De Asis & Kenji Omasa, 2007; Durigon et al., 2014; Van der Knijff et al., 1999). Therefore, a new method for calculating the RUSLE C factor, based on NDVI rescaling, was proposed by Durigon et al. (2014). It is an important model to produce realistic C values for RUSLE and is suggested to be used in tropical regions.

Recently, due to climate change which causes climate variation in some parts of the world, particularly rainfall patterns, this will increase and enhance soil erosion, especially in areas where habitat changes occur. Rainfall erosivity is the potential ability of rain to cause erosion (Gebreegziabher et al., 2021). Thus, Coka watershed faced food insecurity problem, climate change (uncommon rainfall, hotness and coldness), occurrence of drought, lack of crop productivity and dramatic population growth for last 30 years till present (CSA, 2021; WAO (Woreda Agricultural office), 2022), which is a 56.3% increment of population growth with requirement of construction, crops for food security, firewood and charcoal for fuel or income, and grazing. These all anthropogenic activities have significant impacts on the environment and ecology of natural resource management and assurance of food security. Food security depends on soil productivity whereas soil fertility and productivity depend on topsoil nutrient content, which is a reservoir of essential nutrients and plays an important role in cycling essential nutrients. As a result, government and non-governmental organizations were taking action to satisfy the food need of people in the study area for long periods of time but the problem is continued till present with the need of solution. The long period food insecurity and draught problems of study area may be due to loss of topsoil by water erosion in response of land use and land cover change. Therefore, spatial correlation between climate change, land-cover change, and soil erosion using global models, such as RUSLE, can effectively assist in the planning and prioritization of soil and water conservation. Thus, the purpose of this study was to quantify and analyze soil loss due to LULC change by water erosion with the help of RUSLE and NDVI to identify major areas of concern for management interventions in the Coka watershed, Southern Ethiopia.

2. Materials and method

2.1. Study area

The study area is located at about 195 km west of Hawassa city in the southern Ethiopia, 330 km south west of Addis Ababa. Geographically, the Coka watershed is bounded by latitude 7°12ʹ10”−7°18ʹ20”N and longitude 37°31ʹ0”−37°34ʹ25”E covering total area is about 3731 ha (Figure 1). Elevation of Coka watershed ranges from 771 − 2524 m a.s.l. Most of the area in this watershed is under high mountains with very steep slopes. Average annual temperature is 30.0°C. Annual rainfall ranges between 932.7 and 1727.8 mm. The main rainy season is between June and September while the short rains fall between February and April. Coka watershed has bimodal rainfall distribution small (belg) and major (Kiremt). Rains in kiremt season are intensive which uses for crop production and causes soil erosion during this time at different parts in the study watershed. Geologically, the area is covered by trap series volcanics of the tertiary period that are largely acidic rocks such as rhyolites, ignimbrites, and trachites covering the basement complex of the Precambrian rocks (Elias, 2016). Moderate to deeper soils are found in the moderate to gentle slopes and shallower soils are found on steeper slopes in Coka watershed. According to Alemu and Tadele (2018), the soil type of surrounding the study watershed is Nitisol.

Figure 1. Location map of the study area.

In the Coka watershed, the natural vegetation has been damaged by human agricultural activities. Thus, the current vegetation is primarily the secondary vegetation types that restored since the national soil and water conservation project began in 2019. The farming system is dominated by mixed crop-livestock system which consists of intensive and continuous crop cultivation and animal husbandry. The main cultivated crops are maize (Zea mays L.), teff (Eragrostisteff Z.), sorghum (Sorghum bicolor L.), wheat (Triticum durum Desf.), and potato (Solanum tuberosum L.), barley and legumes such as haricot bean, fava bean and pea. The homestead areas are largely planted in enset (Ensetteventricossum), coffee (Coffee Arabican) with undergrowth of root crops (sweet potato, cassava, ginger and yams).

2.2. Data

Time series Landsat (TM, ETM⁺, and OLI-TIRS) satellite data of 1988, 1998, 2008, and 2018 with path 169 and row 055 were used for developing land-use/land-cover maps of the study area and to determine C and P-factors used in the RUSLE model. Landsat data were downloaded from United States Geological Survey (USGS) Center for Earth Resources Observation and Science (EROS) (https://earthexplorer.usgs.gov/). Shuttle Radar Topographic Mission (SRTM) digital elevation data of 30 m resolution were used in this study. Digital Image data files were downloaded in zipped files from the United States Geological Survey (USGS). The Universal Transverse Mercator (UTM) map projection system zone 37 N and datum of World Geodetic System 84 (WGS84) were used for all the datasets to ensuring consistency between datasets during analysis. Field observations were also conducted to fix up training sites, to check ground truth, and verify the final output of the maps. All satellite imageries were geometrically rectified with the help of a topographic map (1:50,000) of the study area obtained from the Ethiopian Mapping Authority. The soil map of the study area was obtained from the soil database compiled by Food and Agricultural Organization, collected from the Ministry of Agriculture, Federal Government of Ethiopia. Remotely sensed data combined with further spatial information in a GIS environment to assess the extent and rate of annual soil loss. Data required for mapping the soil erosion and their sources are given in Figure 2 and Table 1.

Figure 2. Processing flowchart for modeling soil erosion in the catchment.

Table 1. Observation data set and source.

Data type	Format	Source
Rainfall data	Table (xls)	National Meteorological Agency (NMA), Ethiopia
Satellite image	Raster (tif)	United States Geological Survey (http://earthexplorer.usgs.gov/)
Soil properties	Vector image (shf) +Tables (xls)	Food and Agricultural Organization (FAO), Ministry of Agriculture, Federal Government of Ethiopia
Topographic data	Raster (tif)	Shuttle Radar Topographic Mission (SRTM), (http://earthexplorer.usgs.gov/)

2.3. Model selection and parameterization

RUSLE model was used to estimate and map soil loss rates in the study watershed, however, it includes numerous improvements, such as factors and incorporation of the influence of surface. Improved empirical equations were used for the computation of rainfall erosivity (R), eroded factor (K), and crop management (C) factor. Moreover, some parameters of the model have been adjusted for Ethiopian situations (Hurni, 1985). Today computed in a raster-based GIS, the model provides information on potential water erosion on a cell by cell basis. Annual soil loss is a function of five factors described by the expression (Renard et al., 1997):

(1) $\mathrm{A}=\mathrm{R}\times \mathrm{K}\times \mathrm{LS}\times \mathrm{C}\times \mathrm{P}$(1)

where R is the rainfall-runoff erosivity factor (MJ mm ha⁻¹ h⁻¹ y⁻¹), K is soil erodibility (t ha h ha⁻¹ MJ⁻¹ mm⁻¹), LS is slope length-steepness, C is land-use and land-cover; and P is the erosion control (conservation support) practices factor (dimensionless, ranging between 0 and 1).

2.3.1. Rainfall and runoff erosivity factor (R)

Rainfall-runoff erosivity is the primary factor causing soil erosion and accounts for about 80% of the soil loss (Wischmeier & Smith, 1978). The R factor governs the forces of rainfall to disconnect and move soil particles in a specified area. It is computed by multiplying Kinetic energy (E) of rainfall and maximum intensity of rain in 30 min (Renard et al., 1997; Wischmeier & Smith, 1978). In areas where rainfall intensity data is not commonly available, total annual rainfall has been adjusted to estimate R-factor for different regions (Bewket & Teferi, 2009; Gebreyesus & Kirubel, 2009; Le Roux et al., 2008; Xin et al., 2011). As reported by Gitas et al. (2009) and Paveen and Kumar (2012) it can be also analyzed by using annual rainfall. The R factor value was calculated based on mean yearly precipitation for the period of 30 years, which was collected from 23 meteorological stations nearby the Coka watershed, using the erosivity computation formula of Hurni (1985):

(2) $\mathrm{R}=\mathit{\hspace{0.17em}}-8.12+0.562\mathit{\hspace{0.17em}}\mathrm{P}\mathit{\hspace{0.17em}}$(2)

where R is rainfall erosivity factor, and P is annual average rainfall (mm).

2.3.2. Soil erodibility factor (K)

The K-factor determines the property of soil particles that govern their resistance to detachment and transport by rainfall. It is strongly related to the soil properties, plays an important role in soil conservation strategies (Shabani et al., 2014), and reveals the soil erosion per rainfall erosivity index (Paveen & Kumar, 2012). The most important soil properties that affect soil erosion are soil organic matter content, soil texture, drainage ratio, and soil structure (Renard et al., 1997; Robert & Hilborn, 2000; Wischmeier et al., 1971). Soil data in Ethiopia often do not have detailed information about these soil parameters e.g., structure and permeability. Textures are among the principal factors affecting K-values. Soil erosion depends on the relative easiness of the soil for detaching and transporting forces. As stated by Sheikh et al. (2011) and Paveen and Kumar (2012), there is a relationship between organic matter, texture, and K-factor.

In this study, the K factor was estimated from the textural class and soil organic matter due to the high costs of obtaining the data of all four soil parameters (Randolph, 2004). Besides, some investigators have estimated soil erodibility only from soil color (Bewket & Teferi, 2009; Gebreyesus & Kirubel, 2009). For this study, 20 composite soil samples from the top 0 − 30 cm were collected from land-use and land-cover classes and were analyzed for soil texture, and organic matter in Hawassa University’s soil laboratory (Table 2). Particle size distribution was determined by the Bouyoucos hydrometer method (Van Reeuwijk, 1992), and textural classes were made following the USDA system of textural classification. Soil organic carbon was estimated using the Walkley and Black (1934) wet digestion method, and the content of organic matter was obtained by multiplying the content of organic carbon by a factor of 1.724. The mentioned soil parameters were used to compute the K-factor based on the following Equation:

(3) $\mathrm{K}=2.1\times 10-6\times \mathrm{M}1.14\times \left(12-\mathrm{OM}\right)+0.325\times \left(\mathrm{P}-2\right)+0.025\times \left(\mathrm{S}-3\right)$(3)

Table 2. (K) factor values corresponding to the soil texture and organic matter content Randolph (2004).

Textural class	Organic Matter (%)
	≤2	2–4	≥4
Fine sand	0.16	0.14	0.10
Very fine sand	0.42	0.36	0.28
Loamy sand	0.12	0.10	0.08
Loamy very fine sand	0.44	0.38	0.30
Sandy loam	0.27	0.24	0.19
Very fine sandy loam	0.47	0.41	0.33
Silt loam	0.48	0.42	0.33
Clay loam	0.28	0.25	0.21
Silty clay loam	0.37	0.32	0.26
Silty clay	0.25	0.23	0.19

where M is (percentage silt + percentage very fine sand) (100% clay), OM is the percentage of organic matter, P is permeability class, and S is structure class.

2.3.3. Slope gradient factor (LS)

The dimensionless slope length and steepness factor, LS, represent the effect of slope gradient on soil loss, can be determined as a product of the slope length (L) and slope steepness (S; Panagos et al., 2015). This factor was calculated from the Shuttle Radar Topography Mission (SRTM), digital elevation model (DEM) of the watershed at a 30 × 30 m spatial resolution obtained from the United States Geological Survey (USGS) Earth Explorer. Therefore, the LS factor maps were created using ArcGIS (Carlos et al., 2016; Lulseged et al., 2017; Zhang et al., 2013). For a large basin, the grid resolution is always an important factor that influences the model output (Wang et al., 2012). The L is influenced by grid size and steepness while S is only by the steepness so the change of grid size affects the steepness values which accordingly affect the L and S factor values (A.H. Oliveira et al., 2013). Thus, the LS factor was computed using the following equation as recommended by Simms et al. (2003).

(4) $\mathrm{LS}\mathit{\hspace{0.17em}}={\left(\frac{\mathit{FA}\ast \mathit{cellsize}}{22.13}\right)}^{\mathit{m}}\ast {\left(\frac{sin\left(\mathit{\beta }\right)/0.01745}{0.09}\right)}^{\mathit{n}}$(4)

where FA is Flow Accumulation, cell size is DEM size (30 × 30 m), $\mathit{\beta }\mathit{\hspace{0.17em}}$ is slope gradient (in radians). The mean LS value for the study area is m = 0.5, and n = 1.2 were recommended by Mitasova and Mitas (1999) and Liu et al. (2000).

2.3.4. Estimation of cover management factor (C)

The cover and management factor (C) is expressed as the soil loss ratio from an area with a certain cover and management. It might be the most important factor to predict the past, actual, and future soil erosion severity because of its dynamism in short term (Van der Knijff et al., 2000; Karaburun, 2010). This factor was prepared from the corresponding land-use and land-cover maps and fixed values were assigned to every LULC class. The C-factor is dimensionless and ranges from near 1 and 0. Where C is 1 means the surface is treated as barren land, C is 0 means very strong cover effects and well-protected soil. On the contrary, the removal of vegetation for cropping and/or grazing leads to the depletion of Soil Organic Matter (SOM) and increased runoff and soil erosion. To account for such differences, the land-use and land-cover of the study area were derived from a supervised classification of the Landsat satellite image and verified with field information collected from 200 samples taken from training sites. The seven land-use and land-covers were identified overall classification accuracy, using ERDAS imagine 9.2. Hence, four C-factor raster grids were derived for 1988, 1998, 2008, and 2018 to assess the impact of LULC change on soil erosion. Descriptions of the land-cover categories of the study watershed are shown in Table 3.

Table 3. Description of land-use and land-cover classes in the Coka watershed.

Land-use types	Description
Cultivated land	This land-use type includes areas currently under crop, fallow land under preparation and rural settlement housing.
Bare lands	Land use type no vegetation cover, or with rocky land.
Settlement	Urban area (MudullaTwon extension)
Bush/shrub land	Land covered by small trees, bushes/shrubs.
Grasslands	Land covered by grass and small shrubs dominated by grass
Forest lands	Area dominated by high and dense natural forests and dense plantation trees including eucalyptus.
Water body	Back flow water of Gilgel gibe 3 dam.

Normalized Difference Vegetation Index (NDVI) was used in this study for gaining information about the seasonal growth of vegetation condition, vegetation dynamics and as one of the input parameters for estimating the potential of erosion using RUSLE model. The NDVI, which is positively correlated with the amount of green biomass and gives an indication of differences in green vegetation coverage (Van der Knijff et al., 2000). Many researchers calculated the C factor with different equations (Durigon et al., 2014; Van der Knijff et al., 2000; Karaburun, 2010), and the sum of these measurements using the formula Eq. (5)(5) $\mathit{C}=\frac{\mathit{\hspace{0.17em}}\left(-\mathit{NDVI}+1\right)}{2}$(5) :

(5) $\mathit{C}=\frac{\mathit{\hspace{0.17em}}\left(-\mathit{NDVI}+1\right)}{2}$(5)

where C is land-use and land-cover factor, and NDVI is Normalized Difference Vegetation Index.

(6) $\mathit{NDVI}=\frac{\left(\mathit{NIR}-\mathit{RED}\right)}{\mathit{\hspace{0.17em}}\left(\mathit{RIR}+\mathit{RED}\right)}$(6)

where NIR is the reflectance value in near-infrared band, red is the reflectance value in the visible red band.

2.3.5. Estimation of support practice factor (P)

The P-factor represents the role of conservation practices in reducing erosion. In this perspective, soil erosion can be reduced by adjusting the flow pattern and direction of runoff. The type of conservation measures implemented and require conserved areas to be mapped can determine the value of P (Daniel et al., 2015). As observed during field visits, contour plowing is the main soil erosion control practiced among the farmers in cultivated lands improved by the construction of soil and stone bunds in some parts in the study area. The P-factor modifies the potential of erosion by contouring, strip cropping, and terraced contour farming (Kim & Julien, 2006; Kuok et al., 2013; Wischmeier & Smith, 1978). According to Wischmeier and Smith (1978), the support practices (P) factors can be computed for agricultural lands of different slope classes but non-agricultural lands are aggregated together and the P-values range from 0 to 1 depending on the slope classes of arable lands and non-arable lands. P-values can be derived from either image classification using remote sensing data, or even expert knowledge (Panagos et al., 2015). Some studies recommended that the P-value is highly varied due to slope inclination (Karamage et al., 2016; Lufafa et al., 2003; Wenner, 1980; Wischmeier & Smith, 1978); it is calculated by farming practices (Stone & Hilborn, 2012). On the other hand, the P-values suggested by Bewket and Teferi (2009) and Asmamaw and Mohammed (2019) by considering two types of land-uses agricultural and non-agricultural land uses and their corresponding slope. For this study, the P factor was computed by the cultivation and slope based on land-use map and DEM database as used by Shin (1999). To further understand the process behind the LULC dynamics, major causes of the changes were explored using key informant interviews and focus group discussions, including from local communities, government officials, and NGO workers in the study site.

3. Results

3.1. Rainfall and runoff erosivity factor map

The various RUSLE factors identified in the study area are shown in Table 5 and Figure 3. The R-factor values varied within ranges 653.8 − 680.5 (low), 680.5 − 711.3 (moderate), and 711.3 − 760.1 (high) MJ ha⁻¹ y⁻¹ over area coverage of 23%, 45%, and 32%, respectively (Figure 3). Table 4 shows the average, maximum, and minimum annual R factor values for each rain gauge were reported. Moreover, the monthly R factor distribution was reported. The effect of rainfall on soil erosion is high (760.1 MJ ha⁻¹ y⁻¹) in the northern part and gradually decreases at the southern part due to low elevation and presence of forestland, bush-land, and grassland, reaching a minimum erosivity value of 653.8 MJ ha⁻¹ y⁻¹. Moreover, it has been observed two dominated rainfall periods (June to July and August to September) recorded about 80% of the total mean rainfall. Despite the spatial variation of the erosivity factor, the entire watershed is highly erosive as the R-factor is by far greater than the global average R-factor of 760.1 MJ ha⁻¹ y⁻¹.

Figure 3. Thematic maps of R-factor (a), K-factor (b), LS-factor (c), C-factor (d), and P-factor (e) of the Coka watershed.

Table 4. Area distribution of RUSLE factors values in Coka watershed.

R Factor (Mj ha^−1 yr^−1)			Area
			(ha)	(%)
653.8 − 680.5			858	23
680.5 − 711.3			1679	45
711.3 − 760.1			1194	32
K factor	Rank
0.10 − 0.15	Low		410.4	11
0.15 − 025	Moderate		2537.1	68
0.25 − 0.37	High		783.5	21
Ls factor
0 − 2			1007.4	27
2 − 5			8-Dec	19
5 − 10			895.4	24
10 − 15			671.6	18
15 − 30			447.7	12
LULC types	C-factor	Soil loss (t ha^−1yr^−1)
Forest land	0.001	2.3	859	23
Bushland	0.01	23.3	685	18.3
Grassland	0.01	23.3	631	16.9
Urban area/settlement	0.05	21.9	23	0.6
Bare land	0.6	236.2	66	1.8
Cultivated	0.15	65.8	1435	38.5
Water body	0	0	32	0.9
LULC types	Slope (%)	P-factor
Agricultural land use	0 − 5	0.1	20	1.4
	5 − 10	0.12	58	4
	10 − 20	0.14	152	10.6
	20 − 30	0.19	609	42.4
	30 − 50	0.25	489	34.1
	50 − 100	0.33	107	7.5
			1435	100
Non-agricultural land-use	All	1	2296

Table 5. Soil erodibility (K) factor values of different textures and organic matter content.

Sample site	LULC types	Northing	Easting	OM %	Soil Texture class	K-value
TC1	FL	340298	798994	2.21	Loamy Sand(LS)	0.10
TC2	GL	340160	798622	2.98	Sandy Loam(SL)	0.24
TC3	FL	34095	799779	3.20	Loamy Sand(LS)	0.1
TC4	BuL	340002	798665	3.20	Sandy Loam(SL)	0.24
TC5	BaL	341107	802536	2.48	Sandy Loam(SL)	0.24
TC6	CL	341062	802861	1.43	Sandy Loam(SL)	0.27
TC7	FL	340238	805701	3.07	Loamy Sand(LS)	0.1
TC8	FL	340963	802826	2.80	Sandy Loam(SL)	0.24
TC9	CL	341361	804929	1.52	Silty clay loam(SiCL)	0.37
TC10	GL	341301	806789	2.01	Silty clay(SiC)	0.25
TC11	CL	341224	806876	1.58	Silty clay loam(SiCL)	0.37
TC12	CL	341323	803765	1.91	Sandy loam(SL)	0.27
TC13	BuL	341312	801547	2.63	Silty clay loam(SiCL)	0.32
TC14	FL	340293	802985	3.63	Clay loam(CL)	0.25
TC15	BaL	340959	799240	1.76	SandyClay Loam(SCL)	0.28
TC16	GL	340954	799779	3.07	Sandy Clay Loam(SCL)	0.25
TC17	FL	337357	802524	4.67	Clay(C)	0.19
TC18	BUL	339897	802548	2.64	Clay(C)	0.23
TC19	GL	341122	802275	2.22	Clay(C)	0.23
TC20	FL	340859	801022	3.17	Clay(C)	0.23

T-Tembaro, C-Coka watershed, CL- cultivated land, FL-forest land, BUL- bush land, BaL-bareland, GL- Grassland.

3.2. Soil erodibility factor map

The soil erodibility map in the study area is reported in Figure 4. Soil erodibility ranged from 0.10 − 0.15 (low), 0.15 − 0.25 (moderate) and 0.25 − 0.37 (high) within area coverage of 11%, 68% and 21%, respectively (Table 4). K-values show an increasing trend in order such as loam sandy < clay < sandy loam and sandy clay loam < silty clay loam (Table 5). The result noticed that lowest K-value (0.10) is showed high organic matter and loamy sand soil texture in forestland, bush-land, and foot slopes (Figure 4). The highest K-value (0.37) has resulted from low organic matter and silty clay loam soil texture in cultivated land, bare-land and grassland, and middle and upper slopes (Figure 3).

Figure 4. Spatial distribution of soil loss within the study area.

3.3. Slope gradient factor map

The LS-factor represents the overall contribution of topographic features such as length and steepness of the slope on soil erosion. The LS factor is expected in the catchment due to the high value of the contributing inputs. The topography accelerates progressively the runoff downslope and increases the erosivity of runoff. The result indicated that ranges from 0 to 30 and considered 0 − 2 (lowest) and 15 − 30 (highest) represented by 27% and 12%, respectively (Figure 3 and Table 4). About 54% of total areas were found in LS classes between 5 − 30. The result revealed that high LS-factor predominates the southern, and middle parts to eastern highlands and those areas are characterized by very steep slopes.

3.4. Cover and management practice factor map

The dimensionless C-factor values ranged from 0.60 in bare-land to 0.001 in dense forestland with area coverage from total watershed 1.8% and 23.0%, respectively. The observed considerable change in LULC is likely to impact and expose the watershed to severe soil erosion due to the drastic increase of high C-factor in space and time. The highest area coverage is by cultivated-land covering 38.5% with C-value of 0.15 (Figure 4, Table 4). This showed that the risk in the order such as forestland < bush-land < grassland < built-up-area < cultivated-land < bare-land which was with minimum (2.33) in forestland and maximum (236.16) t ha⁻¹ y⁻¹ in bare-land (Table 4), Which fact can increase the surface runoff due to the absence of vegetation while the latter plays a role in dissipating the rainfall energy before it impacts the soil. The proper conservation measures are the alternative to counter the high occurrence of soil erosion. P-values ranged from 0.1 − 0.33 for agricultural land-use based on the slope gradient classes and 1 for all non-agricultural land-use (Figure 3). Area coverage of the highest (42.4%) and lowest 20 ha (1.4%) represented by P-values of 0.19 and 0.1 from agricultural land-use areas (Table 4). This signifies inexistent conservation measures and a high risk of soil erosion.

3.5. Distribution of soil erosion

The soil erosion rate and the total soil loss are presented in Figure 4. Amounts of annual soil loss and priority of land management were used to classify the degree of severity of erosion and priority into six classes as <5 (very low), 5 − 10 (low), 10–18 (moderate), 18 − 50 (high), 50 − 100 (very high) and >100 (severe) soil erosion risk with coverage of the study area 8%, 6%, 13%, 11%, 20%, and 42%, respectively (Table 6). The actual erosion rate is high (18 − 50) to severe (>100) t ha⁻¹ y⁻¹) soil erosion distributed on steep slope landscape, cultivated-land, bare-land, and high rainfall area of northern parts and east-west middle and upper parts of Coka watershed. About, 2723 ha (73%) of the total study area was found to have soil loss rates above 18 t ha⁻¹ y⁻¹. The soil loss exceeding 100 t ha⁻¹ y⁻¹ covered 1567 ha (42%) of the total Coka watershed which was cultivated and poorly vegetated steep slopes area. The mean soil loss rate in the area was 62.15 t ha⁻¹ y⁻¹ which the soil depth may be lowering by 4.97 mm y⁻¹.

Table 6. Estimation of soil erosion rates.

Soil loss (t ha^−1yr^−1)	Erosion risk	Priority area	Area
			(ha)	(%)
<5	Very low	V	299	8
5–10	Low	VI	224	6
10–18	Moderate	III	485	13
18–50	High	IV	410	11
50–100	Very high	II	746	20
> 100	Severe	I	1567	42

4. Discussion

Habitat cover is an important factor affecting soil erosion. As known, the vegetation has an impact effect on soil erosion by protecting the surface soil either from the raindrop impact or runoff detachment and transport, and this effect for soil conservation usually increases as the vegetation cover increases, which is named as the C factor in RULSE. The spatial distribution and pattern of vegetation can greatly influence the flow process and thus greatly correlated with soil erosion (Atoma et al., 2020; Chakrabortty et al., 2020b; Gebreegziabher et al., 2021; Pal & Shit, 2017; Tadesse et al., 2017). Land-use and land-cover changes in the Coka watershed for last 30 years (during 1988–2018) shows that most of the area is cultivated-land and rural-settlement, bare land, built-up-area and water body have increased at an annual rate of 23.1, 2.2, 0.8 and 1.1 ha/year but forestland, bush-land and grassland have decreased at an annual rate of 14.4, 4.1 and 8.7 ha/year, respectively (Tadele et al., 2022). In contrast of that noted by Tadele et al. (2022), projected that the increase in cultivated-land, bare-land, built-up-area, and waterbody but decreases the forestland, grassland and bush-land from 2018 to 2060 years. This indicates that LULC change was occurred in the Coka watershed during specified periods. In the Coka watershed, the RUSLE model has revealed that areal distribution of mean annual soil loss quantity varied from <5 to >100 t ha⁻¹ y⁻¹ which highly affected north, northern-east, and northern-west parts.

Present results are in conformity with Asmamaw and Mohammed (2019) demonstrated that it varied from less than 10 to over 80 t ha⁻¹ y⁻¹ in Ethiopia. Soil loss rates between 18 − 50, 50 − 100, and >100 t ha⁻¹ y⁻¹ are above 18 t ha⁻¹ y⁻¹which is above the soil loss tolerable limit (loss tolerance of 18 t ha⁻¹ yr⁻¹ and formation rate of 2 − 22 t ha⁻¹ y⁻¹) for Ethiopia developed by Table 5, (Hurni, 1983a, Hurni, 1983b). Similarly, FAO (2019) showed that the threshold soil loss value between 4.5 − 11.2 t ha⁻¹ y⁻¹ established by USDA and a range for threshold soil loss between 1 t ha⁻¹ y⁻¹ for shallow sandy soils and 5 t ha⁻¹ y⁻¹on deeper and well-developed soils set by the European Environment Agency. The average annual soil loss rate from the study area was 62.15 t ha⁻¹ y⁻¹which may be lowering the soil depth by 4.97 mm. This study finding shows that about 2723 ha (73%) of the total area was found to have soil loss rates above 18 t ha⁻¹ y⁻¹. This finding is in line with Asmamaw and Mohammed (2019) above 75% of Gerado catchment and Abate (2011) over 50% of Borena district that the rate of soil loss exceeds the tolerable soil loss rate estimated for Ethiopia. These results agree with Addis and Klik (2015) reported that soil erosion is greater than rate of soil formation which reduces the existing production systems with decreased soil productivity. Similarly, Panagos et al. (2016) and FAO (2019) reported that annual crop yield reduction was estimated between 0.1 − 0.4% due to soil erosion. The loss of topsoil leads to a decrease in nutrient-holding and nutrient-supplying power of soil (FAO, 2019).

Concerning the deviations (%) between the analytically and empirically estimated R factors, there is no clear remark on the prevalence of a certain method, since more than one display similar behavior by Arnoldus (1977) equation seems to perform best). The R-values exceed from 50 MJ ha^–1yr⁻¹ which the amount and degree of erosion depend on the rainfall precipitation variability in the study area as opined by Meliho et al. (2016). R-values exceeding 50 MJ ha^–1yr⁻¹indicated that the entire watershed is subject to a high climatic aggressiveness. Management factors that improve the soil carbon, cover and conservation efforts could reduce the effective K factor (Gebreegziabher et al., 2021). The highest K-value (0.37) has resulted from low organic matter and silty clay loam soil texture which is least resistant and highly susceptible to soil erosion in the study area. According to (Das et al., 2020) silt particles tend to be more susceptible to soil erosion than sand or clay particles. These findings are in line with Asmamaw and Mohammed (2019) reported that the K-values ranged between 0.19 − 0.42 in the Gerado catchment, Ethiopia. Similar to the K-values of most soils of Ethiopia vary between 0.05 − 0.6 (FAO, 1984), tropical soils 0.06 − 0.48 (El-Swaify et al., 1992).

Soil erosion is primarily driven by modifications in land-use and management. The mean annual soil loss LULC changes ranged between 2.33 in forestland and 236.16 t ha⁻¹ y⁻¹ in bare-land with average soil loss 62.15 t ha⁻¹ y⁻¹ which may be lowering the soil depth by 4.97 mm in the study area (Table 4). All LULC types expect forestland observed as soil erosion hotspot areas with requirement of management priority. The present results agree with Asmamaw and Mohammed (2019) reported that land-use and land-cover and topographic factors as the major agents of soil erosion in Ethiopia. Bui et al. (2010) and FAO (2019) demonstrated that vulnerable soils soil loss of 1 t ha⁻¹ y⁻¹ equivalent to a value of soil depth lowering of 0.08 mm.

The length and steepness of a slope affect the area and are accounted by the LS-factor in RUSLE model. In addition to steepness and length, other factors such as compaction, consolidation and disturbance of the soil were also considered while generating the LS-factor. Erosion increases with slope steepness but, in contrast to the L-factor representing the effects of slope length, the RUSLE makes no differentiation between rill and inter-rill erosion in the S-factor that computes the effect of slope steepness on soil loss (Carlos et al., 2016; Lulseged et al., 2017; Zhang et al., 2013). Approximately, 54% of total study areas were found in LS-classes between 5 − 30 which are largest LS-value range and highest area coverage (Table 4). Most parts of steep slope study area faced from high (18–50) to severe (>100) t ha⁻¹ y⁻¹ erosion hazard which was bare-land, cultivated and poorly vegetated steep slopes area. This finding shows that the severity of erosion hazard is governed by the combined interaction of environmental (topographic) and anthropogenic (LULC change) factors. According to Zhou et al. (2013), when the slope steepness was < 10° and > 40°, the mean soil loss varied from 11.6 and 59.7 t ha⁻¹ y⁻¹, respectively. The same author states that high soil loss from steep slope areas was due to high runoff from steep gradients. This is supported by FAO (2019), erosion increases in a linear relationship with slope increases, and the velocity of the runoff increases as the slope gradient increases.

Soil erosion effect might become up with the consequence of land degradation and lowering soil productivity in the north, northern-east, and northern-west parts in the study area. Management interventions might be reduced soil erosion and improve soil productivity by increasing the soil macro and micronutrients availability for plants. Daniel et al. (2015) reported that the reduction of soil erosion following improvement in the soil organic matter content. Similarly, Ndah et al. (2015) proposed deduction of livestock as a conservation device that possibly reduces erosion by reducing runoff in the grassland. According to (FAO, 2019), the loss of nutrient-supplying power can be replaced by increasing use of fertilizer while there are economic costs to increased fertilizer use as well as environmental costs to minimize agrochemical contamination of surface waters. Additionally, Bizoza (2014) reported that soil amendments through animal wastes reduce bulk density and compaction which ultimately reduce runoff and soil erosion by increasing pore spaces and infiltration capacity. Therefore, the values obtained from RUSLE modeling were regarded as indicators of soil erosion and associated long-term environmental change rather than the tool for precise estimation of net soil loss.

5. Conclusion

A quantitative assessment of soil erosion in the Coka watershed using the RUSLE equation in a GIS interface successfully provided a set of spatially different management approaches to stop erosion in the study area. The process of calculating RUSLE factors, especially the up scaling methods for R-factor and C-factor, may be helpful for the relative studies. The study findings have identified the presence of high to severe soil erosion in most parts of the Coka watershed of north, north-eastern, north-western, Ethiopia. The study divided the Coka watershed into five erosion severity categories and areas coverage from 299 ha (8%) to 1567 ha (42%) of the study area were observed by the soil loss rates from very low (<5) to severe (>100 t ha⁻¹ y⁻¹), respectively. This indicates that large amount of study area covered by potentially severe erosion. From this, spatial variability in soil loss severity shows the necessity to prioritize land management interventions to the hotspots of soil erosion in the study area. The predicted amount of soil loss and its spatial distribution can provide a basis for comprehensive management and sustainable land-use for the watershed. The areas with high and severe soil erosion hotspots warrant special priority for the implementation of control measures and ultimately improve the livelihood of the local communities. Thus, the government and local community should make mechanisms to protect loss of topsoil by erosion and mitigate climate change which can improve soil fertility and productivity to ensure food security of local community in the study area.

Data availability

All data sets are available on request from the corresponding author (Tadele). The data are not publicly available because they contain information that could compromise the privacy of research participants.

Acknowledgments

We are thankful to the Center of Environmental science at Addis Ababa University and Wachemo University for providing academic and financial support. Mr. Tadele also thanks Dr. Surya, for constant encouragement. We also wish to thank the two anonymous reviewers, and Dr. Daren Jones, Associate Editor, whose suggestions were helpful in significantly improving this manuscript.

Disclosure statement

No potential conflict of interest was reported by the author(s).