Integrated use of GIS, RS and USLE model for LULC change analysis and soil erosion risk mapping in the Lake Hawassa Watershed, Southern Ethiopia

Abstract

Soil erosion is a worldwide problem leading to loss of nutrient rich top soil resulting in reduction of agricultural productivity that threatens the livelihood of agrarian societies. Assessment of soil loss due to water-induced erosion using Models is essential to quantify and map erosion prone areas for possible application of soil conservation measures. In the current study, GIS, RS, and USLE model were integrated with the aim to quantify LULC classes, estimate soil loss rates, and identify priority areas for conservation in the Lake Hawassa Watershed for the two temporal instants, 1975 and 2020. The use of thematic recoding and model maker during LULC classification process has improved the accuracy of thematic classes, producing more than 89% for both classified images. The LULC change analysis showed that agricultural land and built-up area expanded by 0.44% and 15.26% per annum, respectively over the last 45 years. The declines in vegetation cover coupled with steep slopes have contributed to the increase in soil loss from 1975 to 2020 by 65%. The model estimated annual soil erosion rate in 1975 that ranged from nil to 158.77 t ha⁻¹yr⁻¹ with a mean annual soil loss rate of 25.16 t ha⁻¹yr⁻¹. Likewise, the annual soil loss estimated for the year 2020 varied from zero to 262.48 t ha⁻¹yr⁻¹ with a mean annual soil loss rate of 28.20 t ha⁻¹yr⁻¹. Generally, a total of 23,580.59 and 27,783.35 tons of soil have been lost per annum for the year 1975 and 2020, respectively, which is a direct reflection of vegetation clearance. The identified highest three severity classes combined accounted for 50.54% in 1975 and 55.40% of the entire soil loss in 2020 and these losses occurred on areas less than 2% in both 1975 and 2020 periods. Besides, 23.48% and 31.04% of soil loss rates exceeded the established tolerable limit of soil loss for the study area. The integrated use of USLE model, GIS and RS have identified erosion risk areas that could help apply conservation measures.

Keywords:

• Soil erosion risk map
• Lake Hawassa Watershed
• LULC change
• USLE model
• GIS
• RS

1. Introduction

Water erosion is the main cause of soil degradation (Efe et al. 2008; Bekele and Gemi 2021) and about 80% of the world’s agricultural land suffers from moderate to severe erosion (Ritchie et al. 2005) and it is one of the biggest threats to agricultural productivity in Ethiopian highlands (ibid.). Of the total global agricultural soil degradation, water erosion accounts for 56% of the observed damage (Gashaw et al. 2017). Soil erosion by water has been challenging and continuous problem in Ethiopia for decades (Hurni 1988) and land degradation in this country accounts for 8% of the global total, the most series problem being the removal of fertile topsoil by water (Ayalew and Selassie 2015). Water and soil are vital natural resources to maintain environment and livelihoods of humankind, particularly for countries with agrarian economy like that of Ethiopia. Mainly, water is at the heart of most efforts to address economic growth, energy production, human health, food security, and poverty reduction (Berhanu et al. 2013). On the other hand, when the raindrops, which are the major sources of water, fall on a soil devoid of vegetation, the water becomes trouble. Hence, land with no vegetation cover to control soil erosion by reducing the direct impacts of raindrops and the velocity of runoff will always be vulnerable to erosion (Qiang et al. 2016; Negese 2021). According to the Ethiopian highland reclamation study report, 27 million ha or almost 50% of the highland areas were eroded, 14 million ha seriously and over 2 million ha has become beyond reclamation level (Assefa et al. 2015; Bekele and Gemi 2021).

As in most tropical developing nations, the unsustainable and exploitative land use practices have accelerated soil erosion in many parts of Ethiopia (Shiferaw 2011; Tessema et al. 2020) reducing the protective plant cover and exposing the topsoil to high-intensity rainfall. Besides, the dependency of majority of population in Ethiopia on rain fed farming coupled with population pressure, overgrazing, and deforestation have led to accelerated erosion and land degradation. In the highlands of Ethiopia, which carry the highest density of population, soil loss rates on single fields are up to 300 t/ha/year, where rainfall erosivities are highest (Hurni 1993). Soil loss rates are highest from cultivated land, the country average being 42 t/ha/year (Hurni 1993; Shiferaw 2011) and the resulting annual production loss due to soil loss is estimated to be in the order of 2% (ibid.).

The empirical analysis of previous studies confirmed that the mean annual erosion rate of the study watershed was high, 37 t ha⁻¹yr⁻¹ (Degife et al. 2021), however, the soil erosion rate could be much higher where wide and deep gullies are formed. The change in land use and land cover in Ethiopia also played a substantial role in changing hydrological processes such as an increase in surface runoff volumes, reduction of infiltration, and reduction of groundwater recharge in the country (Tesfaye et al. 2017; Negese et al. 2021). Furthermore, LULC change is one of the factors responsible for the sedimentation or siltation of lakes in the country at large and in the Lake Hawassa Watershed in particular. The extinction of Lake Cheleleka with an area of about 12 km² in 1972, estimated average depth of 5 m and a storage volume of ca. 60 × 10⁶ m³ (Ayenew 2004) due to runoff and sedimentation is a manifestation of LULC changes.

To solve the soil degradation problems, wider conservation programs were launched in Ethiopia, after 1970s and 1980s (Bewket 2007; Tessema et al. 2020). Consequently, vast areas have been covered with terraces and millions of trees have been planted in different parts of the country. However, this has not solved the problem of land degradation and the success rate has been minimal (Desta et al. 2005). The reasons for failure are owed to poor watershed management approaches that were especially focussed on physical soil and water conservation practices, which are not based on scientific reasoning like constructing structures without determining the amount of runoff expected (Adimassu et al. 2014). To minimize soil degradation, we need to identify areas with high risk of soil erosion using soil erosion models and thereby use the results to plan and implement the best-suited soil conservation measures.

Currently, there are many models around the world used to estimate soil erosion by runoff, which incorporate the topographic features with land use and soil characteristics with Geographic Information System (GIS) and Remote Sensing (RS) (Star et al. 1997). Among the various erosion models developed globally, Universal Soil Loss Equation (USLE) (Wischmeier and Smith 1978; Balabathina et al. 2020) and its Revised Universal Soil Loss Equation (RUSLE) are the most widely used erosion prediction models (Masha et al. 2021). The USLE model has been widely used as it combines the impact of climate, soil, topography and land use on soil loss with limited data (Wischmeier and Smith 1978; Bekele and Gemi 2021). USLE in its original form does not predict spatial erosion risk status, which was the major limitation of this model. However, integrating the model with GIS had solved this limitation (Fistikoglu and Harmancioglu 2002) and it became a viable option to estimate soil erosion. Considering the effectiveness of integrating erosion models with GIS and RS, many researchers quantified the rate and severity of soil erosion worldwide and in Southern Ethiopia (Oliveira et al. 2015; Wolka et al. 2015; Marondedze and Schütt 2020).

Soil erosion studies have been widely conducted in Ethiopia, however, estimation of spatially distributed soil erosion on grid cell basis and micro units has not been adequately addressed (Balabathina et al. 2020). Similarly, in spite of the presence of erosion in the present study area, the rate of soil loss and its spatial patterns are poorly understood. The study conducted in the Watershed by Wolka et al. (2015) covers only the Lake Cheleleka wetland, which is about 5% of the current study area extent. Therefore, the present study has been initiated to estimate the soil loss rate in the Lake Hawassa Watershed by integrating USLE Model, GIS and RS techniques. The study objectives are to (i) classify remote sensing data from two points in time and map LULC classes, (ii) estimate soil loss rate in the watershed using USLE factors, (iii) develop the spatial pattern of erosion severity map, and (iv) recommend possible management measures for optimal soil conservation intervention. We believe that the soil loss estimates and erosion prone area maps provided a clear understanding of the erosion processes at landscape level and the results will guide methodological decision-making and enable conservation planners to protect the watershed from further damage.

2. Materials and methods

2.1. Description of the study area

The study area is part of the Great Ethiopian Rift Valley and it is located 275 km south of the Capital City, Addis Ababa, along the high way to Moyale, Kenya. The Lake Hawassa watershed is geographically located between 6° 49′ and 7° 14′ N latitude and 38° 17′ and 38° 44′ E longitude and it covers an area of 1302.94 km². It constitutes seven districts, of which four districts (19.17%) fall under Oromia Region while the remaining three districts (80.83%) fall under the newly established Sidama Regional State.

Ethiopia is located in the region where the main rainy season is from June through September. However, Lake Hawassa watershed (Figure 1) has even more extended period of wet season (March to October) with mean monthly rainfall varying from 85 to 133 mm (Gebreegziabher 2004). The mean annual precipitation based on 15 to 56 years of recorded meteorological data amounts to 955.85 mm.

Figure 1. Map showing location of the study area.

Topography of the area exhibits flat low-lying plains, steep hills, and dissected rolling land features with an altitude ranging from 1668 to 2984 m above mean sea level. The mean annual temperature of the area varies from 12.5 to 26.01 °C as computed from the meteorological data recorded for ten years in Hawassa and Wondo Genet stations (Wondrade et al. 2014).

The major LULC of this area is agriculture, which constitutes nearly 60% of the entire watershed with the projected human population of 1,103,507 (Wondrade et al. 2016). The vegetation cover is mainly characterized by natural and plantation forest, at Abaro Mountain, the scrubland located at the northwestern part of the watershed serving as a sanctuary to shelter Swayne’s Hartebeest, and woody vegetation mostly found in the eastern part. Settlement is heavily expanding all over the study area including in the wetland. Currently there are more than 1000 thatched houses illegally built in the wetland indicating the expansion of settlement in the study watershed. Natural forest on the Wondo Genet escarpments have always been under continued threat due to illegal logging, settlement and agricultural expansion, demand for firewood and construction materials. Particularly, northern and northwestern part of the watershed represents the border between Oromia and Sidama Region, where vegetation cover has been depleted leading to severe erosion, which requires establishment of inter-regional collaboration to thwart further threats to ecosystem and thereby reduce erosion risk.

An overlapping pair of calderas namely, Hawassa and Corbeti containing Lake Cheleleka and Hawassa geologically characterizes the study area. The Caldera complex forms a giant elliptical depression 30–40 km wide on the rift floor (Gebreegziabher 2004). The watershed was formed during the tertiary uplifting following the regression of the Mesozoic Sea towards the southeast. The uplifting was associated with the intrusion of basaltic magma through fissures. Later tectonics and geomorphic processes resulted with the present shape of the area (Eshete 2009). The majority of the watershed, the flat caldera floor, is composed of lacustrine sediments of Pleistocene age and the north-west part of the study areas are composed of pumice, un-welded tuffs, obsidian and pitchstone, while the hills and the steep escarpment to the north of Lake Cheleleka are rhyolitic and trachytic lava flows (Belete et al. 2013). The hills of eastern escarpment are composed of Nazret silicic volcanic comprising ignimbrites, un-welded tuffs, ash flows, rhyolites and trachyte, while hilly lands in the east are underlain by rocks of the Dino formation comprising ignimbrites, tuffs, water lain pyroclastic and occasional lacustrine beds (ibid.).

2.2. Data types and sources

The process of soil erosion by water is a complicated system controlled by a multitude of factors comprising LULC, slope of the terrain, climate condition, soil characteristics, conservation practices and interaction between them (Yesuph and Dagnew 2019). Hence, the estimation of soil loss due to erosion requires several data related to these factors. The data required for this study were collected from different sources. Shuttle Radar Topography Mission (SRTM) DEM with a grid resolution of 30 m × 30 m, Landsat 1 satellite image with MSS sensor, and Sentinel 2 A-MSI were downloaded from Earth Explorer, US Geological Survey website, http://www.earthexplorer.usgs.gov. Long-term rainfall data from 11 rain gauge stations found within and in close proximity to the study watershed that ranges from 15 to 56 years were acquired from the National Meteorological Agency of Ethiopia (Table 1).

Table 1. Rainfall stations and their corresponding mean annual precipitation.

Station	Longitude (^o)	Latitude (^o)	Altitude (m)	No. of years	MAP^a (mm)
Hawassa	38^o28′59.2′′	7^o03′53.8′′	1694	51	968.24
Wondo Genet	38^o36′56.1′′	7^o04′55.7′′	1742	56	1043.35
Yirba Dubancho	38^o01′48.1′′	6^o55′45.7′′	2023	47	947.30
Haisawita	38^o33′42.1′′	6^o54′07.4′′	2267	51	1082.30
Shashemene	38^o36′01.5′′	7^o12′01.5′′	1927	46	708.55
Tula	38^o28′29.7′′	6^o57′28.9′′	1873	15	957.38
Shamena	38^o17′25.8′′	7^o04′52.5′′	2072	16	665.93
Belela	38^o25′07.4′′	6^o56′30.8′′	1866	25	887.35
Morocho	38^o24′32.9′′	6^o52′20.3′′	1831	43	1207.81
Aje	38^o21′14.1′′	7^o17′44.9′′	1855	25	858.61
Kofele	38^o46′59.4′′	7^o04′59.2′′	2660	27	1187.59

^a MAP = Mean Annual Precipitation.

To identify the major soil types covering the Lake Hawassa watershed and compute K-Factor, digital soil map of the world published by Food and Agriculture Organization-United Nations Educational, Scientific and Cultural Organization (FAO-UNESCO 1974) was used. The soil map was obtained from Oromia Agricultural Bureau.

2.3. Image data sets and pre-processing

Image datasets were selected based on the available images covering our study area and the amount of funding at our disposal. Accordingly, two satellite images were used in the current study: The first image was Landsat1 acquired on 27 October 1975 and the second was Sentinel-2A acquired on 19 November 2020. The latter is a product of European Space Agency (ESA), but both images were downloaded from Earth Explorer, US Geological Survey website http://www.earthexplorer.usgs.gov. Such images have been used for LULC change analysis in similar scenario (Aneseyee et al. 2019; Eisenberg and Muvundja 2020). Both images that cover our area were cloud cover free and acquired during the same vegetation season not to affect the classification result (Congalton and Green 2009). The two images were checked for geometric correctness with the scanned and mosaicked reference topographic map prepared for spatial compatibility and both images were projected to Universal Transverse Mercator (UTM), Zone 37 North. As part of image pre-processing, the reflective spectral bands of Landsat1-MSS: 4(green), 5(red), 6(NIR), and 7(NIR) with a pixel size of 60 m × 60 m and the reflective bands of Sentinel-2A-MSI: 2(blue), 3(Green), 4(red), 8(NIR), and 12(SWIR) with a pixel size of 10 m × 10 m were stacked into one multi-band image. The SWIR band of Sentinel-2A image was resampled from 20 m to 10 m to match the other bands during stacking process. Part of each raster image data covering an area including DEM was then extracted by importing shapefile of the study area in Erdas Imagine Software Version 2014 to make it readily available for classification. However, before classification, an appropriate classification scheme has to be devised. Such classification scheme should be mutually exclusive and given due attention (Congalton and Green 1993; Wondrade et al. 2014). Since there is no one ideal classification of LULC and considering the heterogeneous character of the landscape, LULC types identifiable by remote sensing are organized on the basis of modified USGS Level I classification scheme (Jensen 2005). Based on the defined classification scheme, the area was described into nine LULC classes (Table 2).

Table 2. LULC classification scheme.

LULC class	Description
Agricultural land	The land occupied by agricultural activities. This category includes areas currently covered by crop, under preparation for cropping, limited areas temporarily left as fallow, and rural farmsteads with sparse bushes and trees.
Forest	Areas dominated by natural and plantation trees, which are greater than 5 m tall. The forest is with an overall mean stand density of >700 stems/ha and a canopy cover of greater than 30%.
Woody vegetation	Land on which trees, bushes, shrubs and grasses are mixed with different proportions. It also contains patches of land covered with perennial crops, such as Coffee, Khat, Enset, and Sugarcane.
Scrub land	Area regarded as Swayne’s Hartebeest sanctuary and covered with dense grasses, bushes, shrubs and small-scattered tress.
Grass land	Areas dominated by herbaceous vegetation. It includes privately and communally owned greeneries.
Water body	Land covered by water body, such as lakes, ponds, rivers, reservoirs
Swamp	Areas where water table is above, near or below the surface of the earth. The area accommodates grassland, few tress and limited thatched houses.
Built-up	Areas characterized by impervious layers, such as buildings of all kinds, paved roads, and playgrounds.
Bare land	Land surface devoid of vegetation, rock outcrops, and quarry sites.

2.4. Image classification and accuracy assessment

Perhaps the most basic form of LULC analysis within the field of remote sensing is LULC classification, which involves the association of features within remotely sensed imagery with specific LULC maps (Aplin 2004; Wondrade et al. 2014). Classifying most African landscapes including the current study area is particularly challenging because of the heterogeneous biophysical features captured by remote sensing instrument and fragmented land-holding system. Classification of the images was performed based on the classes defined in the classification scheme using Supervised (Maximum Likelihood) classification algorithm available in Erdas Imagine software version 2014. Besides, Thematic Recode and Model Maker image correction tools in Erdas Imagine were used to minimize misclassification errors in each thematic class. Supervised classification method demands the need to use training areas to train the computer algorithm the brightness values that will represent one category of LULC in each band of the digital image (Wondrade et al. 2014) having in mind a set of desired classes. Training samples were collected from ancillary data such as mosaic of topographical map sheets from 1976, 1978 and 1979, and aerial photographs from 1972/73 for the Landsat image acquired in 1975. For the image acquired in 2020, training areas were selected from google earth image of 2020. To ensure representativeness of the spectral signature of the LULC class to be discerned, each training area ranges in size from few to large number of pixels. Accordingly, 385 and 1186 training areas were used for the image in 1975 and 2020, respectively.

The overall classification of remotely sensed data cannot be 100 percent accurate as errors originate from different sources including the classification algorithm itself (Smits et al. 1999; Okeke and Karnieli 2006). Before using the classified thematic maps for post-processing, one has to know the accuracy of the classification results. Thus, in remote sensing, accuracy assessment is mandatory (Foody 2002; Okeke and Karnieli 2006) and is important for providing information about the quality of the produced classification result. The process of accuracy assessment was performed in Erdas Imagine software version 2014 by collecting 309 and 809 reference points for the image in 1975 and 2020, respectively. The reference points were generated by stratified random sampling method using topographical maps, aerial photo, and google earth image as a reference data. The collected reference points were compared with the classified image in an error matrix (Foody 2002; Congalton and Green 2009) to check the level of classification accuracy. This will produce accuracy measures that will enable us to accept or reject the classification results.

Besides, the average temporal and annual rates of LULC change dynamicity were estimated for the purpose of comparison. The annual rates of change were computed by dividing the temporal rates of change by the number of years in the interval (Zhang et al. 2010; Wondrade et al. 2016), which is calculated based on the following formula: $R_a=R_t/\mathrm{\Delta }Y,$ and $R_t={(A}_c-A_p)/A_p*100,$ where $R_{a }\mathit{\text{and }}{R}_{t }$ are the annual and temporal change rates of the LULC class in question, respectively. While $A_c-A_p$ are the area of the target LULC class at the current and previous study points in time, respectively, and $\mathrm{\Delta }Y$ is the length of the study period measured in the unit of years.

2.5. Soil erosion model and estimation of soil loss

There are different models used to estimate soil erosion in watershed and such models vary based on the ease with which they are applied and the input data they demand to predict soil loss due to erosion. The Universal Soil Loss Equation (USLE) model proposed by Wischmeier and Smith (1978) is the most widely used model in different countries. The reason for the popularity of USLE model is that it is a good compromise between applicability in terms of required input data and relatively good reliability of obtainable soil loss estimates (Alewell et al. 2019). The model was developed ‘to compute long time average soil losses from sheet and rill erosion’, while considering numerous physical and management variables (Wischmeier and Smith 1978; Eisenberg and Muvundja 2020). The model, integrated with GIS and RS, was applied for the estimation of soil erosion in the study watershed. It estimates the annual soil loss rate in ton/hectare/year by multiplying the respective USLE factor values (R, K, LS, C, and P) as per the recommendations of Wischmeier and Smith (1978), Hurni (1985) and Belay et al. (2020) (Eq. (1)(1) $$A=R \mathrm{*} K \mathrm{*} LS \mathrm{*} C \mathrm{*} P$$(1) ). Each factor is a numerical estimate, which indicates how an area is affected by soil erosion. (1) $$A=R \mathrm{*} K \mathrm{*} LS \mathrm{*} C \mathrm{*} P$$(1) where A is mean annual soil loss in tons per hectare per year (t ha⁻¹yr^–1), R is the rain fall runoff erosivity factor in MegaJoule Millimeter per hectare per hour per year (MJ mm ha⁻¹ h⁻¹ yr⁻¹), K is the soil erodibility factor in tons hour per MegaJoule per Millimeter (t*h MJ⁻¹ mm⁻¹), LS is slope length and steepness factor (unitless), C is the cover and management factor (unitless), and P is the support practice factor (unitless).

The USLE model was implemented on a cell-by-cell basis using Spatial Analyst Tool, Map Algebra, and Raster Calculator in ArcGIS 10.1 environment to investigate soil loss (Hurni 1985; Yesuph and Dagnew 2019). The flow chart of the methods used to estimate soil erosion rate is depicted in Figure 2.

Figure 2. Flow chart of the methodology.

2.5.1. Rainfall runoff erosivity factor (R)

Rainfall-runoff erosivity is the primary factor causing soil erosion and accounts for about 80% of the soil loss (Gezahegn and Arus 2020). Erosivity of rain is the product of storm kinetic energy (KE) and maximum 30-min intensity (EI30) (Wischmeier and Smith 1978; Renard et al. 1997). However, the data on rainfall kinetic energy and rainfall intensity are not available in our study area to compute rainfall erosivity. In this study, we calculated the annual rainfall erosivity using the empirical equation (2)(2) $$\left[R=0.362 \mathrm{*} P+47.6\right],\mathrm{ }$$(2) utilized by Bagegnehu et al. (2019) and applied by other researchers (Bagegnehu et al. 2019; Bekele and Gemi 2021) in Ethiopia. The R factor value was calculated in GIS Environment based on mean yearly precipitation for the period 1965–2020. (2) $$\left[R=0.362 \mathrm{*} P+47.6\right],\mathrm{ }$$(2) where P is mean annual precipitation

Rainfall data from 11 rain gauge stations were used to estimate the rainfall erosivity (R) in the watershed. The missing rainfall data for all areas were filled in using Inverse Distance Weighting (IDW) interpolation method using the tool available in ArcGIS 10.1.

2.5.2. Soil erodibility (K)

Soil erodibility factor (K) is the average soil loss for a specified soil and it is a measure of the susceptibility of soil particles to detachment and transport by raindrop impact and runoff (Wischmeier and Smith 1978; Chadli 2016). All equations for determining K factor are associated to soil organic matter, soil texture, and percentage of sand, silt, and clay in soil (Basri et al. 2019). The observed data of the local soil properties in the Lake Hawassa Watershed are not available. Hence, the major soil types for soil erodibility analysis were extracted from digital soil map in a vector format. The data on soil properties were analyzed following the methods used by several authors (Chadli 2016; Chuenchum et al. 2019; Eisenberg and Muvundja 2020), in which the soil erodibility factors of USLE were calculated using percentages of silt, clay, sand and organic carbon fractions (Eqs. (3)–(7)). (3) $$K_{\mathit{USLE}}=K_{ws}=f_{\mathit{csand}} \mathrm{*} f_{cl-si} \mathrm{*} f_{\mathit{orgc}} \mathrm{*} f_{\mathit{hisand}}$$(3) (4) $$f_{\mathit{csand}}=\{0.2+\left.0.3 \mathrm{*} \exp \left[-0.256 \mathrm{*} m_s\left(1-\frac{m_{\mathit{silt}}}{100}\right)\right]\right\}$$(4) (5) $$f_{cl-si}={\left(\frac{m_{\mathit{silt}}}{m_c+m_{\mathit{silt}}}\right)}^{0.3}$$(5) (6) $$f_{\mathit{orgC}}=\left\{1-\left.\frac{0.25 \mathrm{*} \mathit{org}\mathrm{C}}{\mathit{orgC}+\exp \mathrm{*} \left[3.72-2.95 \mathrm{*} \mathit{orgC}\right]}\right\}\right.$$(6) (7) $${\mathrm{ }f}_{\mathit{hisand}}=\left\{1-\left.\frac{0.7\mathrm{*}\left(\frac{m_s}{100}\right)}{\left(1-\frac{m_s}{100}+\exp \mathrm{*} \left[-5.51+22.9 \mathrm{*} \left(1-\frac{m_s}{100}\right)\right]\right)}\right\}\right.$$(7) where f_csand is the function of high-coarse sand content in soil; f_cl-si is the function of clay and silt in soil; f_orgC is the function of organic carbon content in soil; f_hisand is the function of high sand content in soil; m_s is the percentage of sand fraction content (particles with diameters from 0.05 to 2mm)%; m_silt is the percentage of silt fraction content (particles with diameters from 0.002 to 0.05 mm)%; m_c is the percentage of clay fraction content (particles with diameters of <0.002)%; and orgC is the percentage of organic carbon content of the layer (%).

K-Factor values were assigned for each soil type (Table 3), and the soil map with K values was converted to raster format in ArcGIS 10.1 software package for further processing. The areas of each soil units were calculated and shown in Table 3. Here it has to be noted that the only difference of area in 1975 and 2020 was in water body and swamp classes. The areas of water body and swamp written in brackets are for the year 2020 and this happened for the reason that the area, which was water in 1975, has been converted into swamp (Gh) in 2020.

Table 3. The major soil types identified in the Lake Hawassa Watershed.

Soil unit name	Soil unit symbol	DOM-SU^a	Area (km^2)	Area (%)	K-Factor
Vitric Andosols	Tv	ANz	13.18	1.01	0.146754234
Chromic Vertisols	Vc	VRe	25.74	1.98	0.140277284
Mollic Andosols	Tm	ANm	56.04	4.30	0.128627684
Lithosols	I	LPq	120.75	9.27	0.139088159
Orthic Luvisolls	Lo	LVh	215.42	16.54	0.086107144
Eutric Fluvisols	Je	FLe	478.45	36.73	0.118939786
Chromic Luvisols	Lc	LVx	229.38	17.61	0.137361812
Umbric Gleysols	Gh	GLu	67.05(78.99)	5.15	0.122632491
Water Body	WBD	WBD	96.55 (84.61)	7.41	–
Total			1302.56	100.00	–

^a DOM-SU represents Dominant soil unit.

2.5.3. Slope length and steepness (LS) factor

Slope length and slope steepness are topographic factors that represent the effect of both slope length (L) and slope steepness (S) on soil erosion. LS is defined as the estimated ratio of soil loss per unit area from a field slope to soil loss from a 22.13 m length of uniform 9% slope (Wischmeier and Smith 1978; Balabathina et al. 2020). The slope length L is defined as the distance between the upslope starting point of a slope segment to the downslope point where deposition begins (Panagos et al. 2015; Marondedze and Schütt 2020). Increasing slope length and steepness per unit area results in increased runoff and flow velocity and consequently in increased soil loss (Desmet and Govers 1996). Numerous methods were proposed to calculate the combined LS factor for complex terrain over the past few decades (Desmet and Govers 1996; Zhang et al. 2013). However, a simplified equation suggested by several authors (Wischmeier and Smith 1978; Gashaw et al. 2017; Girmay et al. 2020; Tessema et al. 2020) using a unit contributing area to calculate the LS was used in the current study (Eq. (8)(8) $$LS={\left(\frac{X}{22.1}\right)}^m\left(0.065+0.045 \mathrm{*} S+0.0065 \mathrm{*} S^2\right)$$(8) ). (8) $$LS={\left(\frac{X}{22.1}\right)}^m\left(0.065+0.045 \mathrm{*} S+0.0065 \mathrm{*} S^2\right)$$(8) where $X$ is slope length (m), and it is the product of flow accumulation and cell value, $m$ is a variable slope-length exponent that depends on slope steepness, $S$ is slope gradient in percentage, cell value is the resolution of the DEM used. $m=0.5$ if the % slope is 5 or more, 0.4 on slopes of 3.5 to 4.5%, 0.3 on slopes of 1 to 3%, and 0.2 on uniform gradients of less than 1%.

In this study, the SRTM-DEM with grid cell size of 30 m × 30 m, that was obtained from usgs.gov website was used. The flow accumulation and LS factor were calculated and mapped as a raster layer from DEM using hydrology toolset found in ArcGIS Spatial Analyst Tool box.

2.5.4. Cover and management (C) factor

The cover and management (C) factor represents the ratio of soil loss from land with specific vegetation to the corresponding soil loss from clean-tilled, continuous fallow conditions (Wischmeier and Smith 1978; Morgan 2005; Gelagay and Minale 2016). The cover management is a dimensionless factor that ranges from 0 for a completely non-erodible condition in areas with high vegetation cover to 1 which corresponds to the greater magnitude of soil loss due to very extensive tillage that produces much runoff (Renard et al. 1997; Woldemariam et al. 2018). The vegetation type, stage of growth, and canopy percentage conditions mainly affect the value of C factor (Hagos 2020).

Due to spatial and temporal variations, many studies used remote sensing data to classify LULC classes for quantification of C values (Beskow et al. 2009; Molla and Sisheber 2017). The procedure for classification and analysis of satellite images to generate LULC classes is given in Sec. 2.4. The values of C factor for each LULC classes were obtained from published tables (Table 4). Then, the corresponding C factor value of each LULC class was assigned based on the literature recommendation in the highlands of Ethiopia (Yesuph and Dagnew 2019).

Table 4. LULC classes and values of C factor adopted from published sources.

LULC classes	C value	Sources
Agricultural land	0.15	(Hurni 1985; Tessema et al. 2020)
Water	0.00	(Ganasri and Ramesh 2016; Girma and Gebre 2020; Teshome et al. 2021)
Swamp	0.01	(Bewket and Teferi 2009; Teh 2011)
Built-up	0.09	(Ganasri and Ramesh 2016; Gashaw et al. 2017; Teshome et al. 2021)
Bare land	0.60	(Bewket and Teferi 2009; Tessema et al. 2020)
Forest	0.003	(Wischmeier and Smith 1978; Morgan 2005; Ganasri and Ramesh 2016)
Woody vegetation	0.06	(Bewket and Teferi 2009; Degife et al. 2021)
Scrub land	0.04	(Jain and Kothyari 2000; Eshete 2009)
Grass land	0.05	(Hurni 1985; Yesuph and Dagnew 2019)

2.5.5. Support practice (P) factor

The support practice (P) factor is the ratio of soil loss with a specific conservation practice like contouring, strip-cropping, or terracing measures to the corresponding loss with straight-row farming up and down the slope (Wischmeier and Smith 1978; Gashaw et al. 2017). The P factor ranges between 0 and 1, where the value closer to 0 shows good conservation practices and the values close to 1 shows little conservation practices (Ganasri and Ramesh 2016; Woldemariam et al. 2018).

The P factor could be mapped by collecting data from frequent field observation (Bewket and Teferi 2009; Amsalu and Mengaw 2014). However, in the Lake Hawassa Watershed, there were few attempts to construct conservation measures at limited locations, for which no complete data was found. Besides, the map of conserved areas was not available and hence the P factor was assessed using major LULC types and slope interaction as suggested by Wischmeier and Smith (1978) and similar studies have been applied for Ethiopian conditions (Bewket and Teferi 2009; Gelagay and Minale 2016; Gashaw et al. 2017; Yesuph and Dagnew 2019; Negese et al. 2021). This method categorizes the study area into agricultural land and other lands. The agricultural lands are classified into six slope groups, converted to vector format and P-factor values were assigned, while P-Factor value of 1 was allocated to all non-agricultural lands (Gashaw et al. 2017) including agricultural lands with slope greater than hundred percent, if it exists (Table 5). The slope gradient ranges in percentage were generated from SRTM-DEM in ArcGIS environment, while the types of LULC classes were produced by classifying Landsat-MSS and Sentinel-2A-MSI satellite images. At the end, P factor was assigned by combining the LULC types and slope maps in a raster using Spatial Analyst Tools in ArcGIS 10.1.

Table 5. Slope classes and support practice (P) factor, (Wischmeier and Smith 1978).

LULC classes	Slope (%)	Area (1975)		Area (2020)		P value
		(ha)	%	(ha)	%
Agricultural class/land	0-5	41,739.63	32.18	43,770.57	34.12	0.10
	5-10	10,820.17	8.34	13,823.12	10.78	0.12
	10-20	4240.68	3.27	8539.09	6.66	0.14
	20-30	804.30	0.62	2282.65	1.78	0.19
	30-50	360.27	0.28	154.04	0.12	0.25
	50-100	81.17	0.06	998.86	0.78	0.33
Other LULC classes	All	71,661.60	55.25	58,705.40	45.77	1.00

3. Results and discussion

3.1. Accuracy assessment

Before using image classification results for further analysis, it is mandatory to assess their accuracy to evaluate the percentage of pixels classified correctly and incorrectly per LULC cluster (Foody 2002; Okeke and Karnieli 2006; Wondrade et al. 2014). One of the most common ways of representing accuracy assessment information is in the form of an error matrix or contingency table (Congalton 1991; Congalton and Green 1993; Smits et al. 1999; Congalton 2001; Okeke and Karnieli 2006). Error matrix provides a detailed assessment of the agreement between the sample reference data and classification data at specific locations, together with a complete description of the misclassifications registered for each category (ibid.). According to Congalton (1991) and Congalton and Green (2009), the use of classification measures such as overall accuracy, Kappa statistics, producer’s accuracy, and user’s accuracy is quite common and all accuracy measures are acceptable if greater than 85%.

The result of accuracy assessment for the image in 1973 and 2020 indicated an overall accuracy of 89.97 and 92.83%, respectively, while corresponding kappa statistics of 88.47 and 91.64% for Landsat-MSS and Sentinel-2A-MSI, respectively. The accuracy measures generated by Erdas Imagine software 2014 are presented in Table 6.

Table 6. Accuracy assessment results of the classified images.

No.	LULC class	1975		2020
		PA^a (%)	UA^a (%)	PA^a (%)	UA^a (%)
1	Agricultural land	88.89	82.76	94.51	90.53
2	Water	95.00	95.00	96.15	92.59
3	Swamp	96.15	89.29	90.32	93.33
4	Built-up	72.70	94.12	96.04	96.04
5	Bare land	91.30	93.33	93.28	94.87
6	Forest	94.34	96.15	92.86	90.70
7	Woody vegetation	90.70	88.64	88.37	91.94
8	Scrub land	92.00	92.00	92.94	94.05
9	Grass land	80.00	80.00	90.38	94.00
-	Overall accuracy (%)		89.97	-	92.83
-	Overall kappa statistics (%)		88.47	-	91.64

^a PA and UA stands for producer’s accuracy and user’s accuracy, respectively.

The producer’s accuracies (PA) in the classified maps of 1975 and 2020 were 72–96%, and 88–96%, respectively, while the user’s accuracies (UA) were greater than 80% for a Landsat image and greater than 90% for the Sentinel image. Overall kappa coefficients were all greater than 0.88 showing a strong agreement between the remotely sensed classification and reference data. Though the achieved accuracy measures were reasonable, there were some misclassifications as expected. Such misclassifications were observed in built-up and agricultural land, which could be attributed to difficulties to discern pixels with mixed spectral signatures. The classification error in the thematic map generated from Landsat could also be due to course spatial resolution. Spectral confusions were also noted in woody vegetation and grassland classes.

3.2. Spatial pattern and dynamics of LULC changes

The spatial pattern and extent of LULC classes play a significant role in reducing raindrop impacts on soil particles (Belay et al. 2020). Tropical regions with poor surface cover are highly susceptible to soil erosion and consequent soil loss following rainfall. As the result of classification, nine LULC classes were identified using supervised maximum likelihood classification algorithm (Figure 3) augmented by field observation and ancillary data. The spatial pattern and magnitude of LULC classes extracted from classified images (Table 7) revealed that agricultural land was the most dominant in the study watershed in both 1975 and 2020. Agricultural land covered 44.79% (583.86 km²) and 53.75% (700.36 km²) area in 1975 and 2020, respectively. For agricultural domain in Africa, Brink and Eva (Brink and Eva 2009) analyzed a 57% increase from the year 1975 to 2000. Similar studies conducted in the highlands of Ethiopia indicated an increase in agricultural land from 43–57% (Ali et al. 2011) and 47.9–64.4% (Gezahegn and Arus 2020). Study from the same watershed conducted between 1973 and 2011 also recorded an increase in agricultural land from 43.6 to 56.4% (Wondrade et al. 2014). This shows that agricultural expansion is at a decelerating rate as other land to be converted into agriculture is becoming scarce and agricultural land is losing its area too. The annual rate of change in agricultural land for the current study between 1975 and 2020 was estimated at 0.44%. The main reasons for land resource degradation and particularly agricultural expansion can be attributed to factors such as population pressure (currently growing at a rate of 2.5%), migration and resettlement, environmental pollution, the change in climate, and rapid change in the lifestyle of people among others (Wassie 2020). The population density of the area was predicted to be 870 persons per square kilometre in 2013, which is much higher than the regional and country average (141 and 84.47 persons per square kilometer) for the same year, respectively (Wolka et al. 2015).

Figure 3. LULC maps of the study watershed.

Table 7. Spatial coverage and proportion of LULC classes and change pattern.

LULC classes	1975		2020		LULC Change 2020-1975	Temporal rate of change 1975-2020	Annual rate of change 1975-2020
	Area		Area
	(km^2)	(%)	(km^2)	(%)	(km^2)	(%)	(%)
Agriculture	583.86	44.79	700.36	53.75	116.5	19.95	0.44
Water	98.60	7.56	100.42	7.71	1.82	1.85	0.04
Swamp	79.69	6.11	63.24	4.85	−16.45	−20.64	−0.46
Built-up	5.27	0.40	41.46	3.18	36.19	686.72	15.26
Bare land	51.23	3.93	39.58	3.04	−11.65	−22.74	−0.51
Forest	137.21	10.53	54.80	4.21	−82.41	−60.06	−1.33
Woody veg.	237.34	18.21	207.27	15.91	−30.07	−12.67	−0.28
Scrub land	80.25	6.16	75.02	5.76	−5.23	−6.52	−0.14
Grass land	30.10	2.31	20.79	1.60	−9.31	−30.93	−0.69
Total	1303.55	100.00	1302.94	100.00	–	–	–

A close examination of other lands also revealed that woody vegetation was with the second highest coverage in both classified images, followed by forest in 1975 and water body in 2020. Forest cover declined from 10.53% in 1975 to 4.21% in 2020, while woody vegetation declined at a rate of −1.33% between the two temporal instants. The decline in forest class could have been much more, if it was not for the farmers who are growing Eucalyptus and other trees on their farmsteads to overcome the shortage of construction materials and fuelwood. Of the three LULC classes that increased throughout the study period, viz., built-up, agricultural land, and water body, built-up area increased by 686.72% between 1975 and 2020, which is the highest temporal growth rate recorded between the two dates. Driven by a constantly accelerating urban population in recent decades, urbanization has become one of the most dynamic processes in the context of global land cover transformations and land resource degradation (Bhandari 2010; Wondrade et al. 2016). Such urban expansion processes are occurring at the expense of green spaces. The promotion of Hawassa City as a political centre of the Southern Nations, Nationalities, and People’s Regional State (SNNPRS) and encouraging investment endeavours demonstrated by the establishment of huge industrial parks coupled with the recent upgrading of Sidama Zone to Regional State have contributed to the expansion of built-up area at an annual rate of 15.26%. As presented in the LULC maps of the study area, water body has shown a positive change in area (98.60 km² from 1975 to 100.43 km² in 2020), which is accounted for the surface area expansion of Lake Hawassa despite the conversion of Lake Cheleleka into mudflat swamp area. Lake Cheleleka, which serves as natural trap, has been progressively silted-up over the past nearly four decades at a rate of about 2.25% per annum (Dadi 2020), showing the presence of erosion in the area. The descent in swamp area has shown not only the effect of intensive anthropogenic activities, but also a serious threat to the sustainability of the wetland ecosystem.

Bare land cover declined at an annual rate of −0.51% from 1975 to 2020 due to population pressure and expansion of rain fed agriculture, which elevates land demand and leads to the encroachment to marginal lands that are not suitable for farming. Onsite observations also showed that large extents of bare lands are more prevalent in the western part of the study area and these areas were exposed to gully formations due to erosion. Uncontrolled quarrying of construction materials as a source of instant cash have also led to a wide spread clearance of vegetation cover and soil degradation. Scrubland was the area restricted as Swayne’s Hartebeest Sanctuary, which is located in the northwestern part of the study watershed. However, the area declined from 80.25 km² in 1975 to 75.02 km² in 2020 because of the disturbance imposed by individual intruders and encroaching settlers. Today it is common to see the Swayne’s Hartebeests and domestic animals grazing in unison suppressing the vegetation growth and affecting the wildlife. Grassland occupied the smallest area in 2020 and the second smallest in 1975 next to built-up area. The area of grassland cover has declined at an annual rate of −0.69% between 1975 and 2020. Grassland, which was found in the northern part of the study area in 1975 have visibly been converted into farmland in 2020.

3.3. USLE factors and raster maps

3.3.1. Rainfall erosivity (R) factor

To estimate the erosivity of the study watershed, mean annual precipitation data of 6 and 11 stations were used for 1975 and 2020, respectively. The result of IDW indicated that the mean annual rainfall ranged from 706.41 to 1123.74 mm for 1975 and from 667.02 to 1080.18 mm for the 2020 study periods. The result revealed that eastern part of the study area receives the highest rainfall. It was found that the long-term erosivity values ranged from 302 to 465 and 288 to 450 MJmmha⁻¹ h⁻¹yr⁻¹ (Figure 4a and 4b) for the year 1975 and 2020, respectively. The heaviest rainfall for the study area occurs in the eastern part and it gradually decreases towards the west. It was observed that elevated areas are with more rainfall and consequently induces more erosion due to runoff. Although the rainfall erosivity values of the current study watershed is less than the global average, which is 2000 MJmmh⁻¹ha⁻¹yr⁻¹, the amount of soil loss generated by R factor is substantial (Borrelli et al. 2017; Kidane et al. 2019). Other studies (Ganasri and Ramesh 2016; Kidane et al. 2019) also recorded that of the USLE factors, soil erosion rate is more sensitive to rainfall.

Figure 4. Maps of USLE factors in 1975 and 2020: (a, b) rainfall erosivity (R) factors; (c) soil erodibility (K) factor; (d) slope length and steepness (LS) factor; (e, f) cover and management (C) factors; (g, h) support practice (P) factors in the Lake Hawassa Watershed.

3.3.2. Soil erodibility (K) factor

The erodibility of soil is an expression of its inherent resistance to particle detachment and transport by rainfall, which may depend on the presence of vegetation cover and the erodibility factor (K) of the different soil types. K factors of Ethiopia were computed and the values ranged from 0 to 0.18 t h MJ⁻¹mm⁻¹ (Berhanu et al. 2013). The estimated K-factor values of the study watershed varied between 0.086 and 0.146 t h MJ⁻¹mm⁻¹, the average being 0.127 t h MJ⁻¹mm⁻¹ showing erodibility of the study area is fairly low. This in turn indicates that the influence of soil factor in estimation of soil loss is relatively lower than other factors. Generally, eight types of soil units and water body were identified in the study watershed (Table 3) and the dominant soil type was Eutric Fluvisols followed by Chromic Luvisols and Orthic Luvisols covering an area of 478.45, 229.38, and 215.42 km², respectively. The erodibility map (Figure 4c) shows that Vitric Andosols, Chromic Vertisols, and Lithosols are more vulnerable to soil erosion, with K-factor values of 0.147, 0.140, and 0.139, respectively. These types of soil units were found in Shashemene and Siraro districts, where vegetation cover is highly affected. Here, it has to be noted that K-factor map for the 1975 and 2020 are the same, except that the Lake Cheleleka has been converted into Swamp in the year 2020 and consequently the soil unit is merged with Umbric Gleyosols. The result also revealed that Eutric Fluvisols covered about 37% of the study area, which is predominantly agricultural land.

3.3.3. Slope length and steepness (LS) factor

The LS factor map and flow accumulation were generated using void filled SRTM-DEM with 30mx30m grid cells, while the LS factor was computed using Map Algebra-Raster Calculator of Spatial Analyst Tool in ArcGIS 10.1 environment. The result revealed that the LS factor values ranged from 0 in areas with flat terrain including waterbody to 18.28 in steep slope areas (Figure 4d). Findings of the current study also corroborates with other previous Ethiopian highland study reports (Tessema et al. 2020; Negese et al. 2021). The LS factor was higher in the highland areas, predominantly in the eastern part and along hillsides of the watershed and it is low in the central part due to the flat topography of the terrain. The distribution of LS factor showed that about 99.9% of the watershed corresponds to a LS factor value of less than three. About 50% of the study area (652.36 km²) lies with slope ≤5% (which ranges from flat to gently sloping according to FAO (FAO-UN 2006) slope classification). A relatively small variation of slope length and steepness factor could be accounted for the 50% flat and nearly flat terrain of the watershed.

3.3.4. Cover and management (C) factor

The LULC map of the watershed was generated from satellite images of 1975 and 2020 using Erdas Imagine 2014 Software. The classified LULC change trend indicated that the study watershed experienced substantial changes, which affects the cover and management (C) factor. The values of C factor in the watershed ranges from 0.003 to 0.6 (Table 4) on terrestrial body, the lowest value being assigned to forest and the highest to bare land. The high value of C factor is associated with the area void of vegetation and prone to more soil erosion. Although agricultural land has the greatest coverage in the 1975 and 2020 classified images (Figure 3), vegetation cover and density in 1975 were better, compared to the coverage in 2020, which slowed down erosion process and consequently reduced soil loss rate in 1975. Compared to other LULC classes, bare land followed by agricultural land has higher C factor values and are more vulnerable to erosion risk. The C factor map (Figure 4e and 4f) indicated that the western and most central part of the watershed are areas with high C factor values as these regions were occupied by agricultural and bare lands. These two LULC classes account for about 49% (635.09 km²) and 57% (739.94 km²) of the total study area in 1975 and 2020, respectively.

3.3.5. Support practice (P) factor

Discerned LULC classes and slope in percent were used to generate agricultural and non-agricultural/other lands in ArcGIS 10.1 software package. Since there were no available data on erosion control structures in the study area, P factor values suggested by Wischmeier and Smith (1978) and Gelagay and Minale (2016) were used in each classified agricultural lands and non-agricultural lands. LULC classes such as water, swamp, built-up, bare land, forest, scrub, woody vegetation, and grassland were considered as other lands. Other lands in vector format were assigned P factor value of 1 irrespective of their slope category, while agricultural land was classified into six slope classes and assigned their corresponding P factor values ranging from 0.1 to 0.33 (Table 5). Finally, maps in vector formats were converted into raster to produce P factor maps (Figure 4g and 4h). The result indicated that agricultural land and other lands coverage in 1975 was 44.75 and 55.25%, respectively, while in 2020 agricultural land covered 54.23%, and other lands covered 45.77% of the entire study area. Most part of agricultural land falls under slope class of 0–5% covering an area of 437.71 km². As illustrated by the P factor maps from 1975 and 2020, north eastern, central, and eastern part of the study watershed emerged as areas where support practices have not been in place. Such areas were mainly characterized by forest, woody vegetation; scrubland, grassland, and built-up area cover classes. Except waterbody and swamp areas, most of the areas with high P factor values are also with the highest slope class significantly contributing to soil erosion.

3.4. Soil erosion rate estimation and risk mapping

All USLE factors (R, K, LS, C, and P) were derived as raster layers with spatial resolution of 30x30m after processing the thematic layers of the factors (Balabathina et al. 2020). These raster layers were then overlaid and multiplied using Raster Calculator toolset in ArcGIS 10.1. The USLE model estimates the long-term average soil loss per hectare per year to determine the degree of erosion risk in the lake Hawassa Watershed. The mean soil loss of the watershed was calculated using Zonal Statistics as Table toolset found in Spatial Analyst Tool of ArcGIS 10.1.

The result revealed that the annual soil erosion rate due to rill and sheet erosion in 1975 ranged from 0 in areas with flat terrain to 158.77 t ha⁻¹yr⁻¹ (Figure 5a and 5b) in steep slope areas with a mean annual soil loss rate of 25.16 t ha⁻¹yr⁻¹. The total estimated annual soil loss from the study area for the year 1975 was 23,580.59 t yr⁻¹, which is a direct reflection of the wide spread clearance of vegetation cover. Similarly, the annual average soil loss estimated for the year 2020 ranged from zero in flat areas to 262.48 t ha⁻¹yr⁻¹ in areas with steeper slopes and hillsides with the mean annual soil loss rate of 28.20 t ha⁻¹yr⁻¹. A total of 27,783.35 tons of soil have been lost per annum for the year 2020 due to rill and sheet erosion (Table 8). However, the soil erosion rates in the study area could be much higher where wide and deep gullies are formed. The soil loss rates from this study are comparable to a mean soil loss rate of 29.9 t ha⁻¹yr⁻¹ recently reported on soil loss due to rill and sheet erosion at national level based on 25 observation sites in different parts of the country (Haregeweyn et al. 2017; Molla and Sisheber 2017; Kidane et al. 2019). The estimated soil loss rates of the current study also corroborates with the findings of other similar studies, conducted in the highlands of Ethiopia, by Ayalew and Selassie (2015); Kidane et al. (2019) as well as Masha et al. (2021), who estimated mean annual soil loss rates of 25.8 and 28.7, 24.95, and 29.62 t ha⁻¹yr⁻¹, respectively. However, the mean annual soil loss of our study exceeded the maximum soil loss tolerance (18 t ha⁻¹yr⁻¹) that permits economically sustainable crop productivity as stated for Ethiopian agricultural lands (Hurni 1985; Tamene et al. 2006; Degife et al. 2021; Tamene et al. 2022). In a similar pattern the soil loss rate estimated for our study watershed is higher than the soil loss tolerance (5–11 t ha⁻¹yr⁻¹) reported for the Ethiopian highlands (Shiferaw 2011; Ayalew 2015; Bekele 2021). Soil erosion is a severe environmental and economic problem not only in developing countries, but also in several developed countries. Consequently, the most alarming scenario as reported by Kouli et al. (2009) and Samela et al. (2022) is that in areas with a slow rate of soil formation, any soil loss of above 1 t ha⁻¹yr⁻¹ can be considered as irreversible with in a time span of 50–100 years in European condition. However, we are not able to locate studies of irreversible limit of soil loss rate for Ethiopian condition. Soil formation rates for Ethiopia ranges from 2 to 22 t ha⁻¹yr⁻¹, while the soil formation rate established particularly for the Central Rift Valley Zone where our study area is located vary from 6–10 t ha⁻¹yr⁻¹ (Hurni 1983).

Figure 5. Soil erosion risk map of the study watershed.

Table 8. Erosion severity class, area coverage, and proportion of soil loss.

Soil loss range (t/ha/yr)	Soil erosion severity class	Priority class	1975		2020
			Area, ha (%)	Estimated soil loss, t/yr (%)	Area, ha (%)	Estimated soil loss, t/yr (%)
0-0.054	Nil to Slight	V	102,795.10 (78.896)	599.707 (2.543)	101,971.40 (78.258)	581.951 (2.095)
0.054-3.40	Moderate	IV	26,049.96 (19.993)	11,063.991 (46.920)	26,734.32 (20.517)	11,809.084 (42.504)
3.4-11.80	High	III	1222.56 (0.938)	7214.129 (30.594)	1310.40 (1.006)	7613.976 (27.405)
11.8-83.60	Very High	II	222.48 (0.171)	4452.187 (18.881)	276.48 (0.212)	6745.435 (24.279)
>83.60	Severe	I	2.52 (0.002)	250.573 (1.063)	9.36 (0.007)	1032.904 (3.718)
Total	–	–	130,292.62 (100)	23,580.587 (100)	130,301.96 (100)	27,783.349 (100)

The estimated results indicated that 23.48% and 31.04% of the soil loss rates exceeded the presumed tolerable limit of soil loss rate of 10 t ha⁻¹yr⁻¹ contributing total soil loss of 5536.33 tons/year in 1975 and 8624.69 tons/year in 2020. Areas about 3.02 km² in 1975 and 3.64 km² in 2020 absorbed soil loss rates above the tolerable limit. The result also showed that 8.95 and 18.38% of the soil loss rates were greater than the maximum soil formation rate estimated for the country (22 tha⁻¹yr⁻¹) in 1975 and 2020, respectively.

On the other hand, multiple studies demonstrated that the estimated annual soil loss rates of the current study were lower than the rates of local scale studies that estimated soil loss rate of 279 t ha⁻¹yr⁻¹ in Omo-Gibe river basin, southern Ethiopia with a mean annual soil loss rate of 69 t ha⁻¹yr⁻¹ (Girma and Gebre 2020). Soil loss rate of 716 t ha⁻¹yr⁻¹ in the Koga watershed in the highlands of Ethiopia and 935 t ha⁻¹yr⁻¹ in Beshillo catchment of the Blue Nile basin, Ethiopia were also much higher than our findings (Molla and Sisheber 2017; Yesuph and Dagnew 2019). The mean annual soil loss of Koga watershed and Beshillo catchment was 42 and 51 t ha⁻¹yr⁻¹, respectively (ibid.). Evidently, soil loss profile in the country depicts a strong variation spatially and temporally and generally, there seems disparity between the estimated soil erosion rates in the current study and the figures reported in literature across the country. Such variation is expected and may be accounted for several reasons including the use of models and data inputs, variations in biophysical environment, area extent, and management practices. Despite all the variations, the estimated soil loss rates and their spatial distribution across the field are comparable to what was observed in the watershed. The main contributing factors of soil loss are still inappropriate agricultural practices, deforestation, overgrazing, construction activities, population expansion, and climate change (Montanarella 2004; Wassie 2020). Therefore, mitigation strategies have to be developed to those areas with high, very high, and severe soil loss rates to halt soil erosion and allow soil regeneration.

3.5. Classification and prioritization of critical erosion prone areas for possible conservation planning

Based on the estimated annual soil loss rates, the study watershed was classified into five erosion severity and priority classes (Table 8). The erosion risk level classification was performed following the soil loss rate class proposed by Kadam et al. (2018) with minor modification to suit the features of Lake Hawassa Watershed. The severity classes were assigned mean soil loss values from zero to 83.6 t ha⁻¹yr⁻¹ and verbal labels of ‘slight’ to ‘severe’ erosion classes. Under each year, areas covered by soil loss per severity class were given in ha, whereas the proportions were written in brackets. In areas with scarce data, such prioritization of erosion prone areas is indispensable for practitioners and policy makers to develop and implement interventions to mitigate erosion hazards (Haregeweyn et al. 2013; Gashaw et al. 2017).

The result revealed that in 1975, 21.11% and in 2020, 21.74% of the watershed have been experiencing from moderate to severe soil erosion. Combined high, very high, and sever soil loss (50.53% in 1973 and 55.40% in 2020) have been recorded on areas less than 2%. Highland areas are more affected by soil erosion and surface runoff causing the transportation of sediment to the swamp area, which is located at the lowest altitude serving as a natural trap in a closed watershed. This was manifested by siltation of Lake Cheleleka and surface area expansion of Lake Hawassa due to runoff (Gebreegziabher 2004). Such incidences call for an intervention and conservation measures to be implemented focussing on arresting erosion at higher and steep slope areas. By doing so, we will contribute to Goal Number 15 of the United Nations Sustainable Development Goals (SDGs) established by all member countries by the agenda 2030. Particularly, the SDG target 15.3 aimed at restoring degraded land and soil including land affected by desertification, drought, and floods and strive to achieve a land degradation neutral world by 2030 (UN-SDG 2015).

3.6. Limitations

Both modelling and erosion measurements are prone to uncertainty and may require validation using measured data. However, the study area lacks data for model validation and hence, it was not possible to quantitatively validate the model results. Therefore, the USLE model results were based on comparison with the outputs of previously reported studies using similar methods in the highlands of Ethiopia as used by several previous authors (Gelagay and Minale 2016; Gashaw et al. 2017; Yesuph and Dagnew 2019; Negese et al. 2021). We compared the soil loss estimates of Lake Hawassa Watershed to those from previous studies (FAO-UN 1986; Hurni 1988; Sonneveld et al. 1999; Haregeweyn et al. 2017) conducted in Ethiopian highlands. Besides, this study attempted to provide LULC data using coarse resolution images and soil types were derived from small-scale maps. However, the LULC data would have been more accurate if high-resolution images were used and better soil distribution map was obtained if soil units were extracted from reasonably large-scale maps. Though USLE model is a useful tool to estimate soil loss, it was not developed to estimate gully erosion, which involves complex hydrological and soil erosion processes to have a full picture of soil loss in the study watershed. All the limitations mentioned suggest further research to generate comprehensive results of soil loss for better planning and implementation of soil conservation measures.

3.7. Conservation measures

In data scarce areas such as the Ethiopian highlands, USLE model is helpful in identifying erosion hotspot areas to recommend appropriate conservation measures (Gashaw et al. 2020). The soil erosion risk map under Figure 5 and Table 8 portrayed that the watershed with soil loss severity class of moderate and above covering areas of 27,497.52 ha in 1975 and 28,330.56 ha in 2020 require implementation of soil conservation measures at different levels. Particularly, areas with soil loss rate of greater than 10 t ha⁻¹yr⁻¹ (Hurni 1983) need immediate intervention as this rate of soil loss is beyond the tolerable rate for the Ethiopian highlands and for the Rift Valley zone where our study area is situated to ensure sustainability of the watershed. Although safeguarding the natural resources of an area is the responsibility of every citizen, the first step to protect the soil from degradation is creating awareness among the society living in and around the watershed about the wise use of available resources and then identifying the actors that must develop, mobilize, and implement the planned conservation measures.

In mountainous areas, hillsides, and in areas with soil loss above moderate severity classes, the actors should exert concerted effort to cover these areas with vegetation to reduce soil erosion. There was a program in Ethiopia named ‘Food for Work’ aimed to control erosion and afforestation of degraded lands to stimulate production and increase income (Gedamu 2006). Similar program was also reported in China named ‘Grain for Green’, which attempted to increase land cover and it significantly reduced sediment yield and soil erosion (Wolka et al. 2015). Strengthening such practices coupled with the already on-going ‘Green Legacy’ in the country could help to reduce soil erosion and rehabilitate the affected regions.

Improving traditional farming practices, which invite excessive soil erosion by water and controlling population migration to the forested areas, could play a major role in reducing soil loss due to erosion. Restricting overgrazing and trampling of vegetation by freely roaming livestock that leads to the removal of vegetation cover and allowing concentrated flow of water, which enhances the development of rill and gully erosion, could also help to minimize soil erosion.

The current Ethiopian land policy permits ownership of the land to the state while people are entitled to inheritable use right only leaving uncertainties of land tenure. This kind of possession discourages people not to invest on the land they use and such practices require amendment for better management of land resources, as connection exist between land tenure certainty and sustainability.

To identify mechanical, biological, or agroforestry conservation practices as soil conservation measures, erosion prone areas have to be categorized into farmland management, hillside management, and gully rehabilitation as suggested by Gadisa and Midega (2021) and apply the selected conservation measures for better management. The conservation measure should address not only on-site, but also off-site effects of soil degradation as they seriously threaten human health and have substantial economic implications (Montanarella 2004). If the selected conservation measures are applied at the right location and at the right time, it will contribute to the United Nations-Sustainable Development Goal, agenda 2030 pertinent to soil erosion.

4. Conclusion

Natural resources such as land, soil, water, forest, wildlife, and rangeland are considered as resource base for the Ethiopian population (Wassie 2020). Degradation of one or more of these resources will obviously affect the livelihood of people, who depend on them. Soil erosion is particularly a major problem worldwide leading to loss of nutrient rich top soil, increased runoff, and decreased water availability to plants (Ganasri and Ramesh 2016). USLE model was used to estimate soil loss rates at watershed level to estimate soil loss. This estimation was performed with a goal to identify erosion hot spot areas and apply conservation measures. Advancement in remote sensing and GIS technologies in the last few decades coupled with the increased processing capability of large data have augmented the use of more dynamic models. This study estimated annual soil loss rate, total soil loss, and generated spatial distribution of erosion risk areas, and mapped erosion severity class for two temporal instants, 1975 and 2020. The model result indicated that the estimated total soil losses in the Lake Hawassa Watershed were 23,580.59 and 27,783.35 tons per annum for the year 1975 and 2020, respectively. It is also evident that the majority of soil loss was contributed by a small area that means 97.45% of the soil loss was contributed by 21.10% of the area in 1975, while 87.91% of the soil loss in 2020 was contributed by 21.75% of the study area. Following the decline in vegetation cover and expansion of agricultural land, the magnitude and spatial distribution of soil loss have increased over the period (1975–2020) under investigation. The soil loss map also indicated that some of the areas discerned as erosion risk classes in 1975, persistently continued with the same or more erosion risk class in the year 2020. Such areas were affected by soil erosion for a long time and probably have developed into gullies. This shows that soil erosion risk is worsening in some areas in the study watershed. Indubitably, the erosion severity classes identified could help land managers and conservation experts to plan and implement conservation measures to combat erosion hazards. Many researchers focus on commonly used least cost single conservation methods, such as soil bund, stone bund, etc. alone. However, experiences indicated that single and least cost conservation measures have not produced the required results in protecting soil from erosion and we should learn to invest on combined, better, and efficient conservation measures taking into account dry, moist or wet agro-ecological zones to drain excess runoff sideways to the next waterway or river in each applied measure. This is particularly applicable in areas with high rainfall in moist and wet zones. The USLE model provided useful maps and quantitative data about soil loss of the Lake Hawassa Watershed and researchers and practitioners can apply the same method in different watersheds as long as they are aware of the underlying model concepts and parameters and the possible limitations. The increase in the use of USLE model to estimate soil loss by water erosion from watershed level to continental level and its use in all continents make the outputs of the model usable by authors in different parts of the world. We believe that, the result of the current research output will be one of the publication footprints in the area.

Author contributions

The author performed all the conceptualization, methodology design, software application, resource acquisition, data analysis, and manuscript writing tasks.

Acknowledgements

I acknowledge that this research was supported by Hawassa University. I am also indebted to the National Meteorological Agency for providing long-term rainfall data and the anonymous reviewers for their time to review this manuscript.

Data availability statement

The data presented in this research are available upon request.

Disclosure statement

The author reports there are no competing interests to declare.