Estimation of soil loss using remote sensing and GIS-based universal soil loss equation in northern catchment of Lake Tana Sub-basin, Upper Blue Nile Basin, Northwest Ethiopia

Abstract

Background

Soil erosion is one of the major environmental challenges and has a significant impact on potential land productivity and food security in many highland regions of Ethiopia. Quantifying and identifying the spatial patterns of soil erosion is important for management. The present study aims to estimate soil erosion by water in the Northern catchment of Lake Tana basin in the NW highlands of Ethiopia. The estimations are based on available data through the application of the Universal Soil Loss Equation integrated with Geographic Information System and remote sensing technologies. The study further explored the effects of land use and land cover, topography, soil erodibility, and drainage density on soil erosion rate in the catchment.

Results

The total estimated soil loss in the catchment was 1,705,370 tons per year and the mean erosion rate was 37.89 t ha⁻¹ year⁻¹, with a standard deviation of 59.2 t ha⁻¹ year⁻¹. The average annual soil erosion rare for the sub-catchments Derma, Megech, Gumara, Garno, and Gabi Kura were estimated at 46.8, 40.9, 30.9, 30.0, and 29.7 t ha⁻¹ year⁻¹, respectively. Based on estimated erosion rates in the catchment, the grid cells were divided into five different erosion severity classes: very low, low, moderate, high and extreme. The soil erosion severity map showed about 58.9% of the area was in very low erosion potential (0–1 t ha⁻¹ year⁻¹) that contributes only 1.1% of the total soil loss, while 12.4% of the areas (36,617 ha) were in high and extreme erosion potential with erosion rates of 10 t ha⁻¹ year⁻¹ or more that contributed about 82.1% of the total soil loss in the catchment which should be a high priority. Areas with high to extreme erosion severity classes were mostly found in Megech, Gumero and Garno sub-catchments. Results of Multiple linear regression analysis showed a relationship between soil erosion rate (A) and USLE factors that soil erosion rate was most sensitive to the topographic factor (LS) followed by the support practice (P), soil erodibility (K), crop management (C) and rainfall erosivity factor (R). Barenland showed the most severe erosion, followed by croplands and plantation forests in the catchment.

Conclusions

Use of the erosion severity classes coupled with various individual factors can help to understand the primary processes affecting erosion and spatial patterns in the catchment. This could be used for the site-specific implementation of effective soil conservation practices and land use plans targeted in erosion-prone locations to control soil erosion.

Background

Soil erosion is one of the severe land degradation problems in many parts of the world (Kim et al. 2005). It’s a natural process that varies according to natural and anthropogenic factors leading to increased runoff from more impermeable sub surfaces, loss of nutrient-rich topsoil, soil productivity, biodiversity, and also indirect damage to the environment (Pimentel 2006). Estimates showed that about 85% of land attenuation globally is because of soil erosion reducing crop productivity by about 17%, affecting the soil fertility initially and in the long term resulting in land desertion (Singh and Panda 2017). It is commonly regarded as a major environmental problem and has become a serious threat to soils and sustainable agricultural production, which has effects on the country’s economy. Poor soil management, low soil fertility, high population growth, intensified use of stressed resources and expansion of agricultural frontiers to marginal and fragile lands are quite common in these areas that need attention.

Land degradation is a major concern in many African countries that lack adequate data, infrastructure, scientific support, and land concentration to combat the problem. Severely eroded highlands in Ethiopia is a common feature of the landscape, which is a primary environmental problem and has significant effects on the crop yields, the productivity of the agricultural lands and their overall sustainability. Hurni (1993) estimated that the gross soil loss for cropland was at 42 t h⁻¹ year⁻¹, while the average soil loss for all land in the highlands was at 100 t h⁻¹ year^{−1 (}Jan and David 1995). The estimated annual crop production losses to soil erosion to be 1 to 2% (Bezuayehu et al. 2002; Tenaye 2020) and these losses represent an annual decline of 0.3% in the value of global production (Scoones 2001). This poses a major threat to food security, the livelihood of farmers, and poverty in the country since 80% of its population solely depends on agriculture as a source of employment and income. According to Hurni (1985), the soil erosion rates in Ethiopia are extreme, with about 3.5 billion tons of soil eroded each year from the highlands of which a considerable portion (1.5 billion tons) leaves the country through its rivers. FAO reported that about 27 million hectares of the highland area were significantly eroded and over 14 million hectares of the highland area were seriously eroded and over 2 million hectares are described as beyond the point of no return (FAO 1986; Moges et al. 2019).

The major cause of land degradation in highland areas is soil erosion by water which occurs within rills, inter-rills, gullies, stream channels, and active construction sites (Tsegaye, 2017). Soil characteristics, intense rainfall, unique climatic conditions, complex terrain and land use are important elements governing soil erosion by water (Barber, 1984). The areas with intense population and livestock density recorded soil erosion by water varied from 3.4 to 84.5 t h⁻¹ year⁻¹ and mean annual soil loss of 32 t h⁻¹ year⁻¹ (Berry et al. 2003). As per the report of CSA (2019), decades of continuous population growth rates estimated at 2.7% per annum between 1950 and 2019 have induced an increasing demand for crop and livestock products.

This results in over-farming, intensive grazing, agricultural expansion towards steep slopes and marginal lands and deforestation (Negasi et al. 2018). Consequences of this include shallow soils, diminished water holding capacity, secondary salinization, and reduced forest and woody vegetation cover (current estimate at 1.25 ha year⁻¹ and 1.8 ha year⁻¹, respectively) (Asefa et al. 2003; FAO, 2015). Currently, there is less than 3% forest area and the rate of deforestation was 15,000 ha year⁻¹ in the period 2010–2020 (FAO, 2015; Temesgen et al. 2014). As per FAO (2015) over the last five decades, the impact of land degradation on soil productivity has been declining at an alarming level with many habitable areas are transforming into drylands and wastelands as well as high soil erosion rates. Inappropriate land management practices, land-use changes, as well as intense road construction in fragile areas, are among the most well-known causes of land degradation and desertification. Therefore, it is important to estimate the soil loss, the spatial pattern, and the extent of soil erosion risk within the catchment. The estimates will provide a detailed understanding of the processes and factors that affect soil erosion as well as potential soil losses at land-use change and landscape position. This study is essential to improve land management practices for planning and implementation of effective soil conservation measures for sustainable management and thereby reducing the soil erosion risk.

Soil erosion models have been widely used around the world for estimation of soil erosion by runoff. The most widely applied soil erosion models are Universal Soil Loss Equation (USLE) and its improved version the Revised Universal Soil Loss Equation (RUSLE), the European Soil Erosion Model (Euro SEM), the Soil Loss Estimation Model for Southern Africa (SLEMSA), the Water Erosion Prediction Project model (WEPP), the Agricultural Non-Point Source Pollution Model (AGNPS), and the Revised Morgan–Morgan-Finney model (RMMF). These models were developed for estimating soil erosion best suited for a situation, but applying these models to any other location is problematic, which requires validation and specific local data (Smit, 1999; Jim, 2005). Although various erosion models developed globally, USLE (Wischmeier and Smith, 1978) is the most widely used erosion prediction models, which predicts only the soil losses from sheet and rill erosion under specified cropping and management system condition, because of its simple means of application (Smit, 1999). Most of the improvements from the RUSLE model (Renard et al. 1997) are difficult to adapt to the new locations since it requires detailed, continuous data for its calculation (Kim et al. 2005) such as maximum 30-min rainfall intensity, rainfall kinetic energy and soil erodibility. Such essential data is not available on our research site. Despite there are some limitations in applying USLE for a large area and/or the tropic regions as it was developed at the unit plot scale and temperate environments, we believe the USLE model is suitable for estimating long-term average soil erosion in the study area. The combination of available data sources used with USLE and geographic information system (GIS) technology is a viable option to calculate soil erosion, which would allow targeted attention towards a solution to reduce future soil erosion and provide a first-order method for prioritizing areas to be examined and remediated (Pham et al. 2018; Kim et al. 2005). Soil erosion assessment for many micro watersheds in Ethiopia have been done at a regional scale using lumped methods. However, estimation of spatially distributed soil erosion on grid cell basis and micro units has not been adequately addressed. Therefore, in the present study, an attempt has been made to calculate the soil erosion for the Northern catchment of Lake Tana basin in NW Ethiopia using USLE with remote sensing and GIS technologies. Our study is carried out with the following specific objectives (i) to identify and characterize areas with high erosion potential (ii) to identify the probable spatial relationships between erosion rates and geographic factors namely topography, altitude, land use, soil erodibility, drainage density, and land cover (iii) to understand the soil erosion due to land transformations in the catchment (iv) to calculate soil erosion potential during the dry and the wet seasons and (v) to determine the most important USLE factors that control the soil erosion risk in the study area. This is a first-time study that erosion assessment was conducted with a paucity of data and resources in the catchment.

Study area

The Northern catchment of Lake Tana, a sub-basin of the Blue Nile basin (known as the Abbay river basin), is located in the north-western highlands of Ethiopia, between longitudes 37°00′ and 37°48′ E, and latitudes 12°08′ and 12°46′ N. It covers an area of 292,230 ha or approximately 22.1% of the lake’s total catchment area is at an elevation between 1751 − 3014 m a.m.s.l (Fig. 1). The catchment area is drained by the Gabi Kura, Derma, Megech, Gumero, and Garno Rivers include minor drainage networks of streams. They rise from the slopes of the Ethiopian highlands in the North Gondar zone and flow south into Lake Tana. Most of these rivers are perennial but highly seasonal in their flow, which also acts as tributaries to Lake. The Megech sub-catchment occupied the largest area (27.6%) in the total study catchment area, followed by the Derma (20.6%), Gabi Kura (19.3%), Gumero (18.5%), and Garno (14%). These are bordered by low hills and ridges and drained by many small streams and gullies. The total annual runoff from the northern catchment of Lake Tana is 689.2 × 106 m³ and it is about 9.7% of the total inflow to the Lake Tana sub-basin (Dessie et al. 2015).

Fig. 1

Location of the study area

The study area has a complex terrain that was formed by large-scale tectonic and volcanic activity (Mohr 1971) (Fig. 2). It is composed of rugged mountains, hills with steep slopes and valleys in the north; mostly flat to gently sloping plains flanked by some low hills in south and west; steep slope hills, flat to dissected, undulating flood plains in east and southeast; along the eastern and western margins mostly ridges. The catchment gradient is strongly influenced by the hillslope processes, leading to a major influence on the physical attributes of stream ecosystems, morphology and drainage density. About 11.2% of the catchment area on very steep (> 45% gradient) and steep (20–45% gradient) slopes, leading to severe floods and interruption to stream-banks. Conversely, about 88.8% of the catchment on moderately steep (10-20% gradient), strong (3-10% gradient) and gentle (0-3% gradient) slopes leading to diverse channel gradients. The study area was divided into three major altitudinal zones based on topography and geology: the lowland zone (1751–2000 m), the middle zone (2000–2400 m) and the highland zone (2400–3010 m) cover 61.7%, 24%, 14.2% of the catchment, respectively. These altitudinal zones have a range of slope gradients. The lowland zone shows gentle slopes, strong slopes, hilly slopes and steep slopes 16.8%, 66.1%, 15.4% and 1.7% respectively; middle zone shows 4.8%, 36.5%, 36.5% and 22.1% respectively; and highland zones 3.1%, 27%, 36%, and 34% respectively. A typical dendritic drainage pattern was observed in the study area. However, the curious river course follows the original drainage pattern radiating from volcanic ridges in the mountainous part.

Fig. 2

Terrain and flow direction in the northern catchment of Lake Tana

The climate in the study area is distinctive of the region showing tropical highland monsoon characteristics with different climate zones: Humid subtropical (Cwa) and Subtropical highland oceanic (Cwb) in the highland regions, and Tropical wet and dry or savanna climate (Aw) in the lowland regions (Köppen and Geiger, 1930). Seasonal rainfall in the study is mainly driven by the migration of the Inter-Tropical Convergence Zone (ITCZ) with an average annual rainfall of 1453 mm and reaches its peak in July and August. Most precipitation occurs in a single wet season (called ‘Kiremt’) between June to September, when the ITCZ reaches its northernmost position over northern Ethiopia. It is followed by the dry season (lesser rainfall in October to December called ‘Bega’ and considerably lesser rainfall in January to early May called ‘Belg’). We observed no significant trend in the mean rainfall in any season/month during 1997–2018. The precipitation increases normally with elevation ranging from 618 to 3259 mm in the study area however, it varies over short distances due to fast-changing topographic conditions. The mean annual temperature of 22 years was 19 °C, ranging from 15 to 22 °C in high altitudes and 24–30 °C in low altitudes while the mean monthly temperatures ranging from 15.1 °C in August to 32.5 °C in April.

The major geologic units in the study area have occurred in four distinct episodes of Afro–Arabian Large Igneous Province volcanism (34.05–23.75 Ma) onto the timeline of the Eocene–Oligocene transition during the Cenozoic (Prave et al. 2015) and the dominant rock types are Eocene–Oligocene basalts (lower mafic volcanic unit); Oligocene felsites (felsitic volcanic unit); Oligocene sedimentary rocks, aphanitic and amygdal basalt (sedimentary-basalt association unit); and Oligocene–Miocene basalts (upper mafic volcanic unit). The physiographic setting, drainage characteristics, unique climate conditions, depth of soils, and bi-modal mineral compositions cause various types of soils in the study area. There are six major soil types: Eutric vertisols (VRe, 38.1%), Lithic leptosols (LPq, 28.5%), Haplic luvisols (LVh, 20.5%), Chromic luvisols (LVx, 7.4%), Humic nitisols (NTu, 5.1%), Eutric leptosols (LPe, 0.5%) characterized by loam, clay loam, sandy clay loam, and clay texture (medium to fine granular). Most of the soils occur within 100 cm of the soil surface in the catchment except the Lithic and Eutric leptosols as they occur within 10 cm and 30 cm of the soil surface, respectively. The larger part (67%) of the study area is characterized as imperfect to poor drainage class while the remaining area is moderately well. The land-use types in the study area are namely, cropland, grassland, shrubland, plantation forest, forest, barren land, and the built-up area. Agriculture land-use was widespread in the catchment area, mostly dependent on rainfall, and it is the main economic activity and source of livelihood.

Materials and methods

The USLE model (Wischmeier and Smith, 1978) was applied for calculating soil erosion in the study catchment using the following empirical equation as a product of five major erosion governing factors:

$$ {\text{A}} = {\text{R}} \times {\text{K}} \times {\text{LS}} \times {\text{C}} \times {\text{P}} $$

(1)

where A is the average annual soil loss (tons ha⁻¹ year⁻¹), R is the rainfall erosivity (MJ mm ha⁻¹ h⁻¹ year⁻¹), K is the soil erodibility factor (tons ha h ha⁻¹ MJ⁻¹ mm⁻¹), LS is the topographic factor (dimensionless), C is the land cover management factor (dimensionless), and P is the support practice factor (dimensionless).

In order to apply the USLE model in the study catchment, the input data were collected from various available data sources. We also used improved methods developed by various researchers to calculate the USLE factors where it practically applicable because the required data were not available. The USLE model was implemented in a raster-based GIS environment on a cell-by-cell basis, which could be brought under thorough investigation.

Rainfall erosivity factor (R)

Rainfall erosivity is the erosive power caused by rainfall that causes soil loss, and it can be determined by the product (EI₃₀) of the total kinetic energy (E) of the storm times the maximum 30-min intensity (I₃₀). R factor is also an index of rainfall erosion which is the average annual total of the storm El values in a particular locality. In this study, we calculated the annual rainfall erosivity using the empirical equation proposed by Hurni (1985), because the R factor of USLE requires long-term rainfall intensity data (I₃₀) which are not available for the study area. The equation used to calculate the annual rainfall erosivity is:

$$ {\text{R}} = 0. 5 5*{\text{P}} - 4. 7 $$

(2)

where P is the average annual precipitation (mm).

Rainfall data from 35 rainfall stations within and around the study area over 21 years (1997–2018) were used to estimate the rainfall erosivity (R) in the catchment. The spatial distribution of mean annual rainfalls at all stations of the study area is shown in Fig. 3. The correlation between rainfall, altitude, and topography was assessed. The R-factor was calculated based on Eq. 2 for each rainfall station and the results interpolated in a GIS using the ordinary Kriging method. Many authors have established various relationships to calculate the R-factor using annual and monthly rainfall data (Renard and Freimund, 1994; Diodato and Bellocchi, 2007). We also calculated the R-factor using the other mostly accepted erosion index methods such as Fournier Index (FI), Modified Fournier Index (MFI), and Precipitation Concentration Index (PCI) to study the effect of rainfall length on erosivity are:

$$ {\text{Fournier Index}}\left( {\text{FI}} \right) = \frac{{p^{2} }}{\text{p}} $$

where p is the precipitation in the wettest month and P is the total annual rainfall.

$$ {\text{Modified Fournier Index }}\left( {\text{MFI }} \right) = \frac{{\mathop \sum \nolimits_{{{\text{i}} = 1}}^{12} p_{i}^{2} }}{\text{p}} $$

where $p_{i}$ = the monthly rainfall depth (mm) in i month, and p = the annual rainfall (mm).

Fig. 3

Spatial distribution of annual rainfall in the northern catchment of Lake Tana

$$ {\text{Precipitation Concentration Index }}\left( {\text{PCI}} \right) = 100 \mathop \sum \limits_{I = 1}^{12} p_{i}^{2} /{\text{p}}^{2} $$

where pi is the monthly rainfall and P is the total annual rainfall.

Soil erodibility factor (K)

The soil erodibility factor (K) is one of the key factors of soil erosion which expresses the susceptibility of a soil type to erosion, and is usually regarded as the rate of soil loss per erosion index unit (Wischmeier and Smith 1978). The major soil properties affecting the K-factor are soil texture (sand, silt and clay composition), organic matter content, soil structure index, and the soil permeability index which are used in soil erodibility estimation. Therefore, K is one of the most challenging factors, requiring substantial time, cost, and resources for field surveys and analyses (Bahrami et al. 2005). However, some researchers found a relationship between soil organic matter content, soil texture, and the K factor (Nguyen and Thai 1999). In our study, the harmonized world soil database (HWSD) version 1.21 (FAO/IIASA/ISRIC/ESBN/ISS-CAS/JRC 2013) was used because the above field survey data of soil profiles are not available for the study area. It composes the GIS raster image file linked to an attribute database that provides information about the composition of each soil mapping unit and standardized soil parameters for top and subsoil. The K values for soil types in the study area were estimated based on percent sand, silt and clay composition (soil texture) relative to percent organic matter (i.e., organic carbon multiplied by a factor value, 1.72) at < 2 and ≥ 2 as per the USDA technical manual (USDA 2008). The dominant soil types and associated soils of the catchment are shown in Fig. 4.

Fig. 4

Dominant soil type and the associated soils in the northern catchment of Lake Tana

Topographic Factor (LS)

Topographic Factor (LS) is the slope length-gradient factor that represents the effect of topography on soil erosion rates and it 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). Of the five factors used in the USLE, the LS factor is the most problematic (Renard et al. 1997). The USLE was originally developed on gentle slopes and the LS factor was one-dimensional in calculating the average annual sheet and rill erosion per unit area at the watershed or even larger scales. However, since the topography was shown as two-dimensional, the LS factor is more difficult to estimate than the other terms in the equation (Ligonja and Shrestha 2015; Van Remortela et al. 2004). However, many researchers agreed that the amount of soil loss depended on the three-dimensional distribution of terrain (Mitasova et al. 1996; Moore and Wilson 1992). Numerous methods were proposed to calculate the combined LS factor for complex terrain over the past few decades (McCool et al. 1987; Desmet and Govers 1996; Zhang et al. 2013). A simplified equation suggested by Moore and Wilson (1992) using a unit contributing area to calculate the LS for three-dimensional terrain as follows:

$$ LS = \left( {\frac{{A{\text{s}}}}{22.13}} \right)^{m} \left( {\frac{\sin \beta }{0.0896}} \right)^{n} $$

(3)

where A_s is the up-slope contributing area that could characterize the effect of converging as well as diverging terrains on soil erosion unlike the horizontal slope length (λ) which only applies to 2-dimensional and uniform hill slopes; m (0.4–0.6) and n (1.2–1.3) are the slope length and steepness exponents, respectively; $\beta$ is the slope angle (radian); the values 22.13 m (72.6 ft.) and 0.0896 rad (5.14°) are the length and slope of the standard USLE plot.

For predicting LS factor values on a grid cell basis, the Eq. 3 should be multiplied by (m + l) as proposed by Griffin et al. (1988). Pham et al. (2018) reported that A_s can also be calculated based on the multiple-flow direction algorithm using the digital elevation model (DEM) along with the ArcGIS, as suggested by Freeman (1991). DEM data for this study was obtained from the Advanced Space-borne Thermal Emission and Reflection Radiometer Global Digital Elevation Model (ASTER GDEM 2011) with a spatial resolution of 30 m. In this study, we computed the LS factor using the following equation suggested by Mitasova and Mitas (1999):

$$ \begin{aligned} LS & = \left( {m + 1} \right)*\left( {\frac{{FA*{\text{cell}}\;{\text{size}}}}{22.13}} \right)^{m} \\ & \quad *\left( {\frac{{{\text{sin }}\left( {{\text{slope}}\;{\text{angle}}*0.01745} \right)}}{0.0896}} \right)^{n} \end{aligned} $$

(4)

where FA is the flow accumulation, cell size is the size of DEM data (30 × 30 m), slope angle in radians, and m = 0.5 (0.4–0.7) and n = 1.3 (1.0–1.4) are the exponent values assigned as recommended by Mitasova et al. (1996) and Liu et al. (2000) because the study area has a very complicated terrain with dense stream network resulting in dominating rill erosion and also gully erosion. Determination of the flow direction followed by calculating a raster layer of accumulated flow to each cell was performed using Arc Hydro tools.

Land cover management factor (C)

The C-factor represents the combined effect of all the interrelated cover, crops, and crop management variables on soil erosion rate. It is the most important factor required for land-use policy decisions as it represents conditions that can be most easily managed to reduce erosion. The C-factor is the ratio of soil loss from land cultivated under specific conditions to the corresponding loss from clean-tilled and continuous fallow lands (Wischmeier and Smith 1978). The C-factor values range between 1 and 0 where C equals to 1 shows lack of cover and the C near-zero shows a very strong cover. Assigning the C values for the various land covers by choosing representative values from Tables given by (Wischmeier and Smith 1965) is widely applied in some studies in Ethiopia. By adapting the USLE to the Ethiopian highlands, Hurni also suggested in 1985, the C-factor values for various land cover types. However, for this study, we used satellite images to compute the Normalized Difference of Vegetation Index (NDVI), which is an index of vegetation abundance, for calculating the C values. We chose time-series 30 m resolution of Landsat 8 OLI images: path 170, row 51 from January 11, April 17, June 20, September 24, and December 13 of the year 2018. Remote sensing technology can provide a lot of information about the land surface through the NDVI, which is positively correlated with the amount of green biomass and shows differences in green vegetation coverage (Van der Knijff et al. 2000). We computed the NDVI for each satellite image using the following equation:

$$ NDVI = \left( {\frac{NIR - RED}{RED + NIR}} \right). $$

(5)

NIR and RED are the surface spectral reflectances in the Near-infrared (0.85–0.88 µm) and Red bands (0.64–0.67 µm). The spectral responses depend on many factors that may be affected by various surface conditions. We extracted the spectral reflectances from all 8 bands of the selected time series Landsat-8 imageries. Several studies have established linear relationships between NDVI and USLE C-factor, but the correlations are still quite low (De Jong 1994; Van der Knijff et al. 1999). Many researchers calculated the C-factor with NDVI values by developing different equations (Durigon et al. 2014; Ahmet 2010). However, in this study, the C-factor was calculated using Van der Knijff et al. (2000) equation who found the relationship between them in an exponential function that is more realistic than the linear equation. We generated the mean NDVI from the time series Landsat 8 images data for 5 months: January, April, June, September, and December of the year 2018 and used to create the C-factor map with the following equation:

$$ C = \exp \left[ { - \alpha \cdot \frac{NDVI}{{\left( {\beta - NDVI} \right)}}} \right] $$

(6)

where $\alpha$ = 2 and $\beta$ = 1.

An estimate of C_i for each pixel was also obtained from NDVI of ith time-series map using Eq. 7 to study the trend of soil erosion risk according to the seasonal variation of the C-factor so that variability during dry and wet months can be estimated.

Support practice factor (P)

The support practice factor (P) is defined as the impact of land use or farming system on soil erosion (Pham et al. 2018). It represents erosion prevention practices such as contouring, strip-cropping, and terracing to reduce the amount of soil erosion by the runoff. P-factor is the ratio of soil loss by a particular support practice to that of straight-row farming up and down the slope. P-factor is considered the most uncertain value due to difficulties in its estimation of specific land-use types and farming systems of the specific land plot (Morgan and Nearing 2011). Further, it is time and cost-intensive. P value equal to 1 indicates no erosion control solution. Soil conservation practices are very poor in the study area. Some researchers suggested that the P-value rather depends on the slope inclination (Shin 1999; Lufafa et al. 2003), while others used farming practices on the slope to calculate P values (Stone and Hilborn 2012). We also observed that the soil conservation practices are very poor in the study area, and significant measures are yet to be implemented. Some areas in the lowlands of the catchment have adopted the use of stone walls and eucalyptus trees, but they are poorly maintained. Moreover, the map of conserved areas was not available for the study area and hence this study adopted the P-factor values as suggested by Shin (1999), which includes cultivation method, land-use types, and land slopes. In this, the P-factor is calculated based on the relationship between terracing and slope in the agricultural field areas and is estimated according to the relation of both the farming type and the slope. Slope gradient ranges (in degrees) and types of land-use for this study were generated from Landsat-8, OLI imagery of December 2018 (pan merged 15 m), and DEM of the study area, respectively. To generate a LULC map of the study area, the imagery was first pre-processed to convert raw DN values to top of atmosphere (TOA) reflectance that minimizes spectral differences caused by different factors. The unsupervised classification was performed on NDVI imagery as well as all-original bands of the Landsat data to classify the satellite image into 15 land-use classes and in turn, these classes were classified into 11 classes using supervised classification. Supervised-Maximum Likelihood classification was used for this analysis using ERDAS IMAGINE 2014 software. LULC classes identified in the catchment were: dense forest, plantation forest, shrubland, cultivated land, cash crops, perennial crops, grassland, bare land, bare rock, urban built-up area, and water. After classification, the accuracy of the LULC map was assessed by computing the error matrix. The classification had an overall accuracy of 0.905 (90.5%) with a kappa coefficient of 0.897 and this indicates a strong agreement with the reference data according to Landis and Koch (1977) classification. To assign the P-factor value, we combined the LULC type and slope range maps in raster using spatial analyst tools in Arc GIS and assigned P values to each combination by choosing representative values proposed by Shin (1999).

Estimation of soil erosion risk

All USLE factors (R, K, LS, C, and P) were derived as raster layers of 100 m grid-cell after processing the thematic layers of the factors. These raster layers were then overlaid and multiplied together according to the USLE model for estimating the spatial distribution of average annual soil erosion in the catchment. GIS technologies were used to simulate this model. The total soil erosion from the study catchment and its sub-catchments were also estimated separately by integrating grid-cells with estimated erosion in tons/year. The grid cells were further classified into five erosion severity levels, according to their estimated erosion rates, for identification of critical erosion-prone areas in the catchment: very low erosion (0–1 t ha⁻¹ year⁻¹), low erosion (1–5 t ha⁻¹ year⁻¹), moderate erosion (5–10 t ha⁻¹ year⁻¹), high erosion (10–50 t ha⁻¹ year⁻¹), and extreme erosion (≥ 50 t ha⁻¹ year⁻¹) as proposed by Pham et al. (2018).

Descriptive statistics

The log-linear regression analysis was performed using SPSS 20 to determine the relative sensitivity of each input factor of the USLE affecting soil erosion in the catchment. In this technique, we selected over 100 locations from soil erosion susceptible areas on the derived A-factor map and those values were superimposed with all other independent factors, R, K, LS, C, and P. The Logarithmic transformation of the Eq. (1) was:

$$ \begin{aligned} { \ln }\left( {\text{A}} \right) & = { \ln }\left( {{\text{R }}*{\text{K }}*{\text{LS }}*{\text{C }}*{\text{P}}} \right) \\ & = { \ln }\left( {\text{R}} \right) \, + { \ln }\left( {\text{K}} \right) \, + { \ln }\left( {\text{LS}} \right) \, + { \ln }\left( {\text{C}} \right) \, + { \ln }\left( {\text{P}} \right) \\ \end{aligned} $$

(7)

where ‘ln’ is the natural logarithm.

The Eq. 7 can be expressed as:

$$ \begin{aligned} { \ln }\left( {\text{A}} \right) & = \, \beta \, + \, \beta_{\text{i}} *{ \ln }\left( {{\text{X}}_{\text{i}} } \right) \, + \, \beta_{\text{j}} *{ \ln }\left( {{\text{X}}_{\text{j}} } \right) \, + \, \beta_{\text{k}} *{ \ln }\left( {{\text{X}}_{\text{k}} } \right) \\ & \quad + \, \beta_{\text{l}} *{ \ln }\left( {{\text{X}}_{\text{l}} } \right) \, + \, \beta_{\text{m}} *{ \ln }\left( {{\text{X}}_{\text{m}} } \right) \end{aligned}$$

(8)

where β is a standardized coefficient value of different units of the input factors; β_i is the estimated regression coefficient which quantifies the association between the factor (X_i) and A; ln (X_i) is the natural logarithm of ith input factor.

Drainage density (Dd)

Drainage density (Dd) was defined as the ratio of the total length of streams in a watershed over its contributing area (Horton 1945). The drainage density has been analyzed with many parameters like slope gradient, soil type, elevation, soil erosion, soil erodibility and sediment yield (Collins and Bras 2010; Gregory and Walling 1968). Total drainage density for the entire catchment was calculated with the following equation proposed by Horton (1945):

$$ Dd = \frac{{L_{u} }}{A} $$

(9)

where, Dd = Drainage density; $L_{u}$ = total stream length in a catchment; A = total contributing area of the catchment.

The levels of Dd were derived as well as the dry and wet drainage densities of the catchment were calculated using focal statistics in ArcGIS since the terrain was complex. Based on field observations, we considered the Dd values less than 0.86 km/sq.km as low drainage densities while the Dd values were greater than 0.86 km/sq.km as high drainage densities. The effects of Dd together with other geographical parameters such as elevation, slope, and soil conditions on soil erosion were estimated.

Stream Power Index (SPI)

Although the USLE model accounts for rill and inter-rill erosion, it does not account for soil loss from gullies or mass wasting events such as landslides (Rubianca et al. 2018). Stream power index (SPI) is a measure of the erosive power of flowing water and it appropriates the locations where the gullies are likely to form on the landscape. There are several models to predict gully erosion by considering various factors like topographical, hydrological, geological, and environmental conditions. DEM and GIS techniques were used to map the gully-erosion susceptibility in the study area. Drainage density, slope, soils, and rainfall are critical factors promoting gullies (Arabameri et al. 2019). The effects of these factors on gully erosion were analyzed. The SPI was calculated using the following equation:

$$ {\text{SPI}}_{\text{i}} = { \ln }({\text{DA}}_{\text{i}} *{ \tan }\left( {{\text{G}}_{\text{i}} } \right) $$

(10)

where SPI _i is the stream power index at grid cell i; DA_i is the drainage basin area (flow accumulation at grid cell i multiplied by grid cell area), and G_i is the slope at grid cell i in radians. The critical zones were classified into two: gully erosion and no gully erosion by using a threshold value of 7.2. The areas of extreme erosion potential in the catchment were compared with the areas of gully formation.

Land cover change detection

In this study, Spatio-temporal satellite data were used to evaluate the effect of a typical land cover change on soil erosion rate. Some land cover features, like water bodies and vegetation, have very specific spectral reflectance characteristics that facilitate the separation from other features as per spectral indices. However, separation of built-up areas (especially rooftops) from barren lands (bare soils, bare rock, and sands), sandy gravel grounds, and shallow stream beds using a single index are a critical and often overlooked challenge due to similarities in the spectral properties of the land use classes. Using the spectral indices derived from satellite images for land use mapping is an operational approach as it enables the mapping at a higher degree of accuracy which is highly comparable to the above indistinguishable land use areas (Xu 2007; Zha et al. 2003). The characteristics of the reflected energy in different regions of the spectrum for a specific land cover can be utilized to produce various indices that are used as an alternative data source for land cover characterization. In this study, a multiple index approach was constructed using three-band combinations of spectral indices for mapping seven different land cover classes.

The Indices such as Normalized difference tillage Index (NDTI), Soil adjusted vegetation index (SAVI), and Modified normalized difference water index (MNDWI) were selected in the order as components of multi-index data set following the several combinations of other indices such as built-up indices (NDBI, BUI, BAEI, BSI), vegetation index (NDVI), and water index (NDWI) and their performances were tested in the pre-classification. The maximum likelihood supervised classification technique was employed to create the spectral signatures of significant land cover categories (waterlogged, crop fields, barren land, grassland, shrubland, forest, plantation forest, and built-up area). After ensuring accuracy for each classified image, a detailed post-classification land-use change detection analysis was performed to assess the resultant change from 1986 to 2018. The equations for the above multispectral indices used in this study are:

$$ {\text{NDTI}} = \left( {\left( {{\text{SWIR 1}} - {\text{SWIR 2}}} \right)} \right)/\left( {\left( {{\text{SWIR 1}} + {\text{SWIR 2}}} \right)} \right) $$

(11)

$$ {\text{SAVI}} = \left( {\left( {{\text{NIR }} - {\text{ Red}}} \right)} \right)/\left( {\left( {{\text{NIR}} + {\text{ Red }} + \, 0. 5} \right)} \right)* 1. 5 $$

(12)

$$ {\text{MNDWI}} = \left( {\left( {{\text{Green}} - {\text{SWIR 1}}} \right)} \right)/\left( {\left( {{\text{Green}} + {\text{SWIR 1}}} \right)} \right) $$

(13)

In this study, classification accuracy assessment was conducted by comparing the reference satellite image with the classified image using some random sampling points. A total of 320 random sampling points of the reference image were used for accuracy assessment of every classified image. A stratified random sampling of the image was adopted to calculate the classification accuracy of each classified image in this study. Because in this sampling method, each land cover class found an equal probability is to be observed (Haque and Basak (2017). The accuracy assessment result of each classified image (2018, 1986) was quantified by using the error matrix. According to the error matrix reports, these classifications had an overall accuracy of 93.4% with a kappa coefficient of 0.925, and 88.12% with a kappa coefficient of 0.864 for the satellite images 2018 and 1986, respectively.

The methodology described above was implemented in our study area at the northern catchment of Lake Tana to develop substantial data inputs into the database for estimating soil erosion in Ethiopia.

Results

Rainfall erosivity factor (R)

The average annual rainfall at 35 weather stations within and around the study area varies between 643 mm and 3259 mm. Results using the ordinary least square method showed that a weak positive correlation between the rainfall and altitude, with a coefficient of correlation of r² = 0.04 but do not show any significant differences between them (P > 0.05) (Fig. 5a); whereas a nearly moderate positive correlation between the rainfall and topography, with a coefficient of correlation of r² = 0.12 and showed significant differences between them (P < 0.05 level) (Fig. 5b). This inconsistency is due to the fast-changing topography of the catchment which is surrounded by a chain of mountain ridges and low hills, and alternation between valleys and ridges within short distances (Veeranarayana et al. 2019). The estimated R-factor value ranges from 349 to 1788 MJ mm ha⁻¹ h⁻¹ year⁻¹ and showed that the distribution of rainfall was uneven in the study area (Fig. 13a). About 66.6% of the research site is comprised of the R values greater than 1000 MJ mm ha⁻¹ h⁻¹ year⁻¹. The other erosivity indices namely Fournier Erosivity Index (FI), Modified Fournier Erosivity Index (MFI), and Precipitation Concentration Index (PCI %) results showed that the FI values ranged from 20.1 to 225.12 (avg. 112.63) that indicated low to extremely severe erosion risk (Oduro-Afriyie 1996), and the MFI values ranged from 109.4 to 613.7 (avg. 291.5) that indicated moderate to very high risk (CEC 1992) and the PCI % values ranged from 10.9 to 40.66% (avg. 20.97%) that indicated moderately seasonal to highly seasonal (Oliver 1980). Various indices of rainfall erosivity, annual rainfall, and altitude for different stations in the study are shown in Fig. 6a. The trends of monthly average rainfall at various stations are shown in Fig. 6b. The most precipitation events occurred between June and September that confirm the Köppen climatic classification showing no alternation of humid and dry months within the wet season (Köppen and Geiger 1930).

Fig. 5

a Correlation between rainfall and altitude (is not significant at the 0.05 level), and b rainfall and topography in the northern catchment of Lake Tana (is significant at the 0.05 level)

Fig. 6

a R-factor values from various erosion indices and b average monthly rainfall at 35 stations in the northern catchment of Lake Tana during 1997–2018

The mean annual rainfall during the wet season for the study area was 1173 mm and the R-factor value ranges from 272 to 1305 MJ mm ha⁻¹ h⁻¹ year⁻¹ (Fig. 7a). Therefore, soil loss mostly occurs during the rainy season and particularly during major storms, however, the post- and pre-monsoon months that have less intense events with low rainfall also cause some extent of erosion. During the dry season (October–May), the mean annual rainfall for the total catchment was 260 mm (range 135–478 mm) and the R-factor value ranged from 69 to 256 MJ mm h⁻¹ h⁻¹ year⁻¹ (Fig. 7b). The dry months, March, April, and November have light precipitation (21–33 mm) that occurred in two or three non-consecutive days in each of these months whereas, less intense events with low rainfall in May (105 mm) and October (66 mm). During this time, soil, especially from the areas where the agricultural land is left as it is with no soil conservation measures after sowing or harvesting, and where the land is prepared for new plantation forest, is being eroded considerably even by small rainfall intensities. The Megech sub-catchment has the highest R values followed by the Gumero and Derma sub-catchments while the Garno and Gabi Kura sub-catchments have the lowest R values (Table 1).

Fig. 7

a R-factor map of the northern catchment of Lake Tana during the wet season and b dry season

Table 1 Average annual rainfall and rainfall erosivity (R) of sub-catchments of the northern catchment of Lake Tana (1997–2018)

Soil Erodibility Factor (K)

The results also indicate that soil erodibility factor (K) in the study area ranges from 0.20 to 0.34 t ha⁻¹ MJ⁻¹ mm⁻¹ (Fig. 13b). The results presented in Table 2 show that about 66.8% of the study area has a K-factor value of 0.24, 0.33, and 0.34 for the soil types Eutric vertisols (clay) Lithic leptosols (clay loam), and Eutric leptosols (loam) respectively and is considered as high due to having low permeability, organic matter and imperfect drainage, whereas K-factor value of 0.2 and 0.21 for the soil types Chromic-Haplic luvisols and Humic nitisols respectively and is considered as low due to having acceptable soil permeability, moderately well drainage. The high K-factor value soil types are naturally more prone to soil erosion due to their physical structure, texture, permeability, and organic matter content. The organic matter (OM) content of the soils in the study area ranges from 0.67 to 4.2% and the bulk density ranges from 1.18 to 1.51 kg/dm³. The K-factor values for the sloping landscapes are significant. The area about 60% of the high altitudes, 89% of the mid-altitudes, and 55% of the lowlands in the study comprise the high K-factor values that showed more prone to soil erosion. The results also showed that the high K-factor value of 0.33 and above occurred in the sub-catchments Megech, Gumero, and Garno causing severe soil erosion rates (Table 3).

Table 2 Soil types and soil erodibility in the northern catchment of Lake Tana

Table 3 K-factors for the northern catchment and its sub-catchments of Lake Tana

Topographic Factor (LS)

DEM data show that the terrain of the study area is very complex, with 20% of the area having a slope steeper than 25°. The spatial distribution of the combined LS-factor was derived from flow accumulation and slope using the DEM of the study area (Fig. 13c). The LS-factor values vary from 0 to 96.47. The results presented in Table 4 show that nearly 98.4% of the study area showed the LS-factor value 1 and below which indicates the slope is very complex and gentle, and slope lengths are shorter, resulting in a very high flow rate and extreme soil erosion potential. The high LS values occurred only 1.6% of the study area on steep slopes and high flow accumulation areas, particularly in the upper parts of all sub-catchments. It indicates that this area is highly vulnerable to soil erosion. The LS values 1 and below were found in over 98% of all sub-catchment areas, whereas the LS values 10 and below were in over 1% of the Megech, Garno, and Gumero sub-catchments.

Table 4 Distribution of LS-factors in the northern catchment of Lake Tana

Land cover management factor (C)

The 5 months average of time series NDVI varies from − 0.194 to 0.48 (Fig. 8), while the C-factor in the study area varies from 0.20 to 0.91 (Fig. 13d). The results presented in Table 5 show that about 45.3% of the study area comprised higher C values ranging from 0.63 to 0.91, whereas lower C values comprised only 19.3% ranging from 0.20 to 0.53. This indicates the effects of cropping and management practices, vegetation canopy, and ground cover on soil erosion rates are high in half of the study area. The results also show that the trend of C-factor change in the study area was due to the seasonal variation of the various land covers (Fig. 9). The range of C-factor values during the wet (0.037–1.46) and dry season (0.18–1.35) showed that there was a huge variation in cropping management practices, particularly in highlands to low-lying floodplains of the study area (Fig. 10a, b). Consequently, about 46.1% of the lowland area and 41% of the highland area has a higher C value of 0.63 and above. The higher C values indicate high soil erosion risk.

Fig. 8

Average Normalized Difference of Vegetation Index (NDVI) map of the northern catchment of Lake Tana, 2018

Table 5 C-factor values for sub-catchments of the northern catchment of Lake Tana

Fig. 9

Seasonal variation of the C-factor in the northern catchment of Lake Tana in 2018

Fig. 10

a Seasonal C-factor map of the northern catchment of Lake Tana during the wet season and b dry season

Support practices factor (P)

P-factor depends on the erosion control practices at each land use type. However, land-use type affects the rate of soil erosion through farming support practices (P) on the slope and the land cover management (C) (Figs. 11, 12). The P-factor map for the study area is shown in Fig. 13e and the area under different support practice factors is shown in Table 6. The P-factor in the study area varied between 0.003 and 1.0. P-value 1 indicates that there are no conservation practices. About 21% of the study area has a P-value of 1. The results indicate that significant soil conservation measures are to be implemented in the study area for controlling soil erosion and protecting the soil on marginal and steep slopes.

Fig. 11

Slope map of the northern catchment of Lake Tana

Fig. 12

Land-use types in the northern catchment of Lake Tana in the year 2018

Fig. 13

USLE factor maps of the northern catchment of Lake Tana, 2018 for soil loss estimation. a rainfall erosivity (R) (MJ mm ha⁻¹ h⁻¹ years⁻¹); b soil erodibility (K) (tons ha h ha⁻¹ MJ⁻¹ mm⁻¹); c slope length and steepness (LS); d cover management (C) and e support practices (P) factor

Table 6 P-factor values for land use groups on slope gradients (Shin 1999)

Estimated soil losses and area prioritization

The average annual soil loss was estimated on 100 m × 100 m grid-cells by multiplying raster layers of all USLE factors (R, K, LS, C, and P) in ArcGIS. The total estimated soil loss from the study catchment was 1,705,370 tons and the mean estimated erosion rate was 37.89 t ha⁻¹ year⁻¹, with a standard deviation of 59.2 t ha⁻¹ year⁻¹. The cells with estimated erosion rates in the study area were classified into five groups according to their severity level: very low erosion (0–1 t ha⁻¹ year⁻¹), low erosion (1–5 t ha⁻¹ year⁻¹), moderate erosion (5–10 t ha⁻¹ year⁻¹), high erosion (10–50 t ha⁻¹ year⁻¹), and extreme erosion (≥ 50 t ha⁻¹ year⁻¹) (Fig. 14). The risk map shows that the spatial distribution of soil loss in the catchment was variable. The estimated erosion rates and total annual soil losses in the study catchment are shown in Table 7. The soil loss estimated for the catchment ranges from 0 in the flat areas to over 50 t h⁻¹ year⁻¹ in degraded lands, along the steep channel banks, marginal and very steep slopes. About 51.5% of the study area has a soil loss value of 0. This indicates that the soil erosion from these areas is negligible.

Fig. 14

Estimated soil erosion rates and erosion levels in the northern catchment of Lake Tana

Table 7 Estimated soil loss and erosion rates in the northern catchment of Lake Tana

The results presented in Table 7 show that the cells with estimated erosion rates of 50 t ha⁻¹ year⁻¹ or more (extreme erosion potential) comprised 2.9% of the total study area and the soil loss under this class accounts for 47.1% of the total estimated soil loss (809,848 tons). About 87.6% of the study area was classified as being under very low, low, and moderate soil erosion potential and has a soil loss value of less than 10 t ha⁻¹ year⁻¹. The total soil loss from these areas combined accounts for only 17.9% of the total soil loss of the study area (298,136 tons). Areas with high and extreme soil erosion classes (erosion rates of 10 t ha⁻¹ year⁻¹ and more) represent high-priority areas for implementation of soil conservation practices in the study catchment. Although such classes together comprised only 12.4% of the total study area (36,317 ha), they contributed 82.1% of the total estimated soil loss (1,407,234 tons) with a mean soil erosion rate of 38.75 t ha⁻¹ year⁻¹. This mean erosion rate is greater than the tolerable soil limit estimated for Ethiopia (2–18 t ha⁻¹ year⁻¹) by Hurni (1985). This mean rate is also compared with the adjacent catchments as well as other basins in Ethiopia shown in Table 8.

Table 8 Comparison of soil loss estimations with results of various studies conducted in Ethiopia

In this study, we also calculated the total soil erosion for each sub-catchment by adding together the soil erosion for all the cells within each sub-catchment and were classified into five soil erosion rate classes proposed by Pham et al. (2018). The results presented in Table 9 show that the soil erosion rates in each sub-catchment of the study area. The Megech sub-catchment comprised 27.6% of the total study area. The total annual soil loss from this area was 848,756 tons (sharing 49.9% of the total soil loss), with a mean erosion rate of 40.98 ± 69.91 t ha⁻¹ year⁻¹. The Derma sub-catchment comprised 20.6% of the total study area and the total annual soil loss from this area was 207,933 tons (accounted 12.2% of the total soil loss), with a mean erosion rate of 46.80 ± 71.89 t ha⁻¹ year⁻¹. The Gabi Kura sub-catchment comprised 19.3% of the total study area and the total annual soil loss from this area was 73,392 tons (accounted 4.3% of the total soil loss), with a mean erosion rate of 29.66 ± 25.78 t ha⁻¹ year⁻¹. The Gumero sub-catchment comprised 18.5% of the total study area and the annual soil loss was 325,656 tons (accounted for 19.1% of the total soil loss) with a mean erosion rate of 30.95 ± 36.17 t ha⁻¹ year⁻¹. The Garno sub-catchment comprised 14.1% of the total study area and the annual soil loss was 244,832 tons (accounted for 14.4% of the total soil loss) with a mean erosion rate of 30.01 ± 28.17 t ha⁻¹year ⁻¹.

Table 9 Estimated erosion rates and annual soil loss in sub-catchments of the northern catchment of Lake Tana

The results also indicate that soil erosion rates varied across the five sub-catchments in the study area. Megech sub-catchment has the largest amount of area in the northern catchment of Lake Tana and also has the largest amount of area with high and extreme erosion potential (51.5% and 32.7% of the total high and extreme erosion areas of the study area) followed by the Gumero sub-catchment comprised 18.5% of the study area and the amount (16.8% and 23.3%) of the area with high and extreme erosion potential, The Garno sub-catchment comprised 14.1% of the study area and the amount (16.6% and 19.4%) of the area with high and extreme erosion potential, Gabi Kura sub-catchment comprised 19.3% of the study area and the amount (4.3% and 9.8%) of the area with high and extreme erosion potential and Derma sub-catchment comprised 22.9% of the study area and the amount (10.8% and 14.9%) of the area with high and extreme erosion potential. The areas with high and extreme soil erosion potential in each sub-catchment are the critical areas that require urgent soil and water conservation measures. These differences in reported values could be attributed mainly to the relative strength of influence of erosion governing factors such as topographic, support practices, soil erodibility, land cover, rainfall, and also anthropogenic activities. If measures are not applied to the areas identified as at risk, the agricultural production in these areas will be severely affected, and this will cause food insecurity (Birhan and Assefa 2017).

About 88% of the total soil loss occurred alone during the wet season, with a mean erosion rate of 42.4 ± 77.1 t ha⁻¹ whereas, about 12% of the total soil loss occurred during the dry season, mostly in pre- and post-monsoon months, with a mean rate 34.2 ± 32.7 t ha⁻¹. This indicates that soil erosion is highly seasonal. The annual estimated soil erosion rates during the wet and dry seasons are shown in Fig. 15a and b respectively. The soil conditions affected by the previous rainfall of the wet season played a vital role in the following dry season for land cover, as was seen in the lower parts of the catchment. The early and receding rains also have a major influence on soil erosion in the sub-catchments; Megech has the highest amount (46%) of the total soil loss during the dry season, followed by Gumero (20%), Garno (16%), and Derma (13%).

Fig. 15

a Erosion levels of the northern catchment of Lake Tana during the wet and b dry season

Descriptive statistics

The log-linear regression analysis results showed that the average annual estimated soil erosion rate (A) had a significant correlation and there was no multi-collinearity with each input factor of the USLE model (P < 0.043, VIF < 10). This indicated that the impact of each input factor of the USLE on annual soil erosion rate was significant (Pham et al. 2018). The results presented in Table 10 show that the estimated standardized coefficients, (β) values ranging from 0.151 to 0.563, and were used in Eq. 9 for multiple linear regressions of the average annual estimated soil erosion rate (A) and each input factor of the USLE model for the study area as follows:

$$\begin{aligned} { \ln }\left( {\text{A}} \right) & = 0. 1 5 5*{ \ln }\left( {\text{R}} \right) + 0. 2 1 7*{ \ln }\left( {\text{K}} \right) + 0. 5 6 3*{ \ln }\left( {\text{LS}} \right) \\ & \quad + 0. 1 5 1*{ \ln }\left( {\text{C}} \right) + 0. 4 6 8*{ \ln }\left( {\text{P}} \right) \end{aligned}$$

(14)

The β values in Eq. 14 indicate that the relative influential strength of each input factor on the annual soil erosion rate. The LS-factor had the strongest influence on soil erosion rate (β = 0.563) followed by the other factors, P (β = 0.468), K (β = 0.217), R (β = 0.155), and C (β = 0.151).

Table 10 Standardized coefficients (β) for USLE model-independent factors

Factors influencing soil erosion risk and erosion rates

Land cover

Land use is represented by the values of the land cover factor (C). About 45.3% of the study area has a C-factor value of greater than 0.63. This area comprised 38.6% and 51.7% of the total high and extreme erosion potential area in the study area, respectively. It indicates that nearly all high and extreme erosion areas had barren land use. Only 19.4% of the study area has low C-factor values ranging from 0.2 to 0.53 and it comprised 29.1% and 16.5% of the total high and extreme erosion potential area, respectively (Table 11). About 35.3% of the study area has the C-factor values ranging from 0.53 to 0.63 and this comprised 32.3% and 31.8% of the total high and extreme erosion potential areas in the study respectively. This indicates that nearly all high and extreme erosion areas had agricultural land use.

Table 11 C-factor and estimated erosion in the northern catchment of Lake Tana

Topography

The topography had the strongest influence on soil erosion rates in the study area. The results presented in Table 12 show that various slope gradients and estimated erosion rates in the study area. Strong slopes (3°–10°) comprised 54% of the study area with a mean erosion rate of 25.6 ± 28.8 t ha⁻¹ year⁻¹ and 34.1% of the total estimated soil loss. Moderately steep slopes (10°–20°) comprised 22.8% of the catchment, with a mean erosion rate of 42 ± 66.2 t ha⁻¹ year⁻¹ and 32.2% of the total estimated soil loss. Steep slopes (20°–45°) comprised only 11.1% of the catchment and 30% of the total estimated soil loss with a mean erosion rate of 46.1 ± 37.6 t ha⁻¹ year⁻¹. Very steep (> 45˚) slopes comprised only a small fraction of the entire catchment. Steep and moderately steep lands in the study area comprised over 30% of the high and extreme erosion potential area. There was a high rate of soil erosion from these slopes and it needs more attention to control such high erosion. This could be due to the agricultural land-use on marginal and steep slopes in the study area.

Table 12 Slopes* and estimated erosion in the northern catchment of Lake Tana

Land use

Table 13 shows the land use type and estimated erosion rates in the study area. About 27.6% of the catchment comprised the barren land and 24.8% of the total estimated soil loss, with a mean erosion rate of 35.9 ± 57.2 t ha⁻¹ year⁻¹. This has the highest amount (31.6%) of the area with extreme erosion potential. About 28.7% of the study area comprised agricultural land use and 13.6% of the total estimated soil loss, with a mean erosion rate of 41.1 ± 66.1 t ha⁻¹ year⁻¹, and having 29.7% of the total extreme erosion potential area. About 15.3% of the catchment comprised the plantation forest and 28.7% of the total estimated soil loss, with a mean erosion rate of 34.9 ± 49.3 t ha⁻¹ year⁻¹, and having 21% of the total extreme erosion potential area. About 13.4% of the catchment comprised the shrubland and 12.1% of the total estimated soil loss, with a mean erosion rate of 33 ± 47.5 t ha⁻¹ year⁻¹, and having only 5.9% of the area with extreme erosion potential. 9.7% of the catchment comprised the grassland and 12.4% of the total estimated soil loss, with a mean erosion rate of 41.64 ± 74.2 t ha⁻¹ year⁻¹, and having 9.9% of the area with extreme erosion potential. The mean channel bank erosion rate for the catchment was 27.2 ± 40.1 t ha⁻¹ year⁻¹. Much of this erosion occurs as the weakened stream banks fail and break down due to increased runoff. The results indicate that most of the soil loss occurred from the barren land, cropland, plantation forest, and shrubland. These areas have the largest amount of extreme erosion cells and are contributing to high rates of soil loss.

Table 13 Land-use type and estimated erosion in the northern catchment of Lake Tana

Table 14 shows the land-use types (cropland, plantation forest, barren land) and the respective cells with extreme erosion potential for each sub-catchment in the study area. The extreme erosion cells are found higher in cropland (33.2%), plantation forest (36.3%), and barren land (43.0%) in the Megech sub-catchment than those in all other four sub-catchments. Next to the Megech, Gumero has areas of extreme erosion potential in the cropland (26.6%), plantation forest (22.4%), and barren lands (19.1%). The extreme erosion cells are the least in the three land use (8.7%, 5.9%, and 5.7%) in the Gabi Kura.

Table 14 Land-use type and estimated erosion in sub-catchments of the northern catchment of Lake Tana

Land use and topography

Table 15 shows the slopes versus land-use type and estimated erosion rates. Strong slopes (3°–10°) comprised 54% of the study area, 61% of the agricultural cells, and have the highest amount (39.3%) of areas with extreme erosion potential, whereas Steep slopes (20°–45°) comprised only 8.5% of the agricultural cells, but 27.3% of the cells with extreme erosion potential were located on such slopes. Strong slopes comprised the largest amount (49.3%) of the plantation forest and 27.5% of the cells with extreme erosion potential, whereas, Steep slopes comprised only 21% of the plantation forest and 34.8% of the cells with extreme erosion potential were on such slopes. Strong slopes comprised the largest amount (52.1%) of barren land and 33.5% of the area with extreme erosion potential, whereas, moderately steep slopes comprised only 29.2% of barren land and the highest amount (41%) of cells with extreme erosion potential on such slopes.

Table 15 Slopes, land use types and estimated erosion in the northern catchment of Lake Tana

Drainage density (Dd)

The drainage density (Dd) estimated for the study area varies from 0 to 1.5 km sq km⁻¹ and the overall drainage density was 0.79 km sq km⁻¹. Figure 16 shows the spatial distribution of the Dd of the study area. The overlying map of drainage density and soil loss reveals that about 44.8% of the study area (131,390 ha) has a high drainage density value of greater than 0.86 km sq km⁻¹. The soil loss of the area under this class accounts for 45.1% of the total soil loss of the study area and has 41.3% of the cells with extreme erosion potential. Whereas, about 55.2% of the study area (160,985 ha) has a low drainage density (Dd < 0.86 km sq km⁻¹) and the soil loss area under this class accounts for 54.9% of the total soil loss and has 58.7% of the cells with extreme erosion potential (Table 16). About 29.5% of the catchment is characterized by wet (active) drainage density (86,297 ha), whereas 70.5% of the catchment is characterized by dry drainage density (206,297 ha). The altitude and soils of the study area also affect the rate of soil loss from the Drainage density areas (Sumantra and Padmini 2015). In the high altitudes (upstream section) of the study area, low Dd and soil loss are high. In this study, the C factor and R factor is the reason for high soil loss. In the low altitudes (downstream section), high Dd and soil losses are very high due to several fingertip streams and the shallow depth of the soils promotes the soil loss. The soils, Lithic leptosols, and Eutric vertisols with high erodibility factors occur in both high and low Dd areas. These areas have the largest amount (65%) of extreme erosion cells.

Fig. 16

Drainage density map of the northern catchment of Lake Tana

Table 16 Drainage densities and Soil erosion in the northern catchment of Lake Tana

Soil erodibility and topography

Strong slopes comprised 68% of the Vertisols (K value 0.24) as well as 64.3% of the Luvisols (K value 0.20), and 57.8% and 48.6% of cells with extreme erosion potential, respectively. About 37.8% of the Leptosols (K value 0.33) occurred on moderately steep slopes and 41.2% of cells with extreme erosion potential were located on such slopes. The steep slopes comprised 32.3% of Nitisols (K value 0.21) and 45.2% of cells with extreme erosion potential thereon.

Gully erosion susceptibility

Drainage density, land-use, topography, soil erodibility, and stream power were considered as influencing factors to gully in our study. The Stream Power Index (SPI) values ranged from -1.71 to 16.29 in the catchment. These susceptible areas were then classified into the gully and no gully erosion by using a threshold value of 7.71 (Fig. 17). The gullies occurred in 15.52% of the study area (45,375 ha). About 44.1% of the gully area was identified in the high Dd areas with the presence of 37.6% of cells with extreme erosion potential (3969 ha), whereas the remaining 55.9% of the gully area was in low Dd areas with the occurrence of 62.4% of cells with extreme erosion potential (4606 ha). The gullies in low Dd areas have higher extreme erosion potential than those in high Dd areas in all sub-catchments. The larger areas of gullies were identified in areas of K values 0.24 (Vertisols), 0.33 (Leptosols), and 0.20 (Luvisols) respectively. The gully areas in plantation forest and barren lands were larger than compared with those in the crop-, shrub- and grasslands. 50.5%, 22.2%, 16.5%, and 10.8% of gullies occurred on strong slopes, hilly slopes, steep slopes, and gentle slopes respectively.

Fig. 17

Critical areas of gully and no gullies in the northern catchment of Lake Tana

Land use cover change (LUCC) detection analysis

The output of thematic change detection analysis (Fig. 18) and the corresponding statistics of the land use cover change from 1986 to 2018 (Table 17) showed that the highest amount of land use cover conversion had taken place from barren land to cultivated land (4.29%) and plantation forest (3.69%). A considerable amount of land cover changed from the shrub-, grass- and natural forest lands to cultivated land and plantation forests (1.01–1.30%). Other land changes or conversions were below 1% and not exempt from extreme erosion potential because much of the cells with extreme erosion potential were located even in such smaller land cover changed areas. Results of land use cover change (LUCC) statistics for 1986 and 2018 showed the typical behavior of each land cover type considering the change dynamics with remarkable precision. In summary, forest, barren land, grassland, and shrublands decreased by 64%, 23.4%, 19.0%, and 17.5% of the total land-use area respectively during the 32 years (1986–2018). Conversely, plantation forests and cultivated land increased (56.6% and 32.4% respectively). Artificial surface and waterlogged areas also increased (101.6% and 37.1% respectively) Table 18). The calculated accuracy assessment result of each classified image (2018, 1986) is shown in Table 19. The results in the table were summarized and quantified by using the error matrix. The accuracy results indicate that the classifications in this study had a strong agreement (kappa statistics > 0.8) with the reference data.

Fig. 18

Land use cover change in the northern catchment of Lake Tana during 1986–2018

Table 17 Land use cover change in the northern catchment of Lake Tana from 1986 to 2018

Table 18 Net land-use cover change in the northern catchment of Lake Tana from 1986 to 2018

Table 19 Accuracy assessment result of classified images

Discussion

The results revealed the high spatial variability of soil erosion potential in the catchment. The average annual soil erosion risk map indicates that 12.4% (36,617 ha) of the study catchment is critical due to the presence of the high and extreme erosion-prone grid cells that contributes about 82.1% (1,400,109 tons) of the total estimated soil loss in the catchment. The highest concentration of the grid cells with extreme erosion potential comprised only 2.9% of the total catchment, but it contributes 47.1% of the total estimated soil loss with a mean erosion rate of 156.8 ± 285 SD t ha⁻¹ year⁻¹; and cells with high erosion potential comprised 9.5% of the total area that contributes 35% of the total sol loss, with a mean erosion rate of 20.8 ± 9.8 SD t ha⁻¹ year⁻¹. The mean estimated soil erosion rate from the catchment was 37.89 t ha⁻¹ year⁻¹, with a standard deviation of 59.2 t ha⁻¹ year⁻¹. This average erosion rate for the catchment is greater than the 18 t ha⁻¹ year⁻¹ soil loss tolerance limit estimated for the Ethiopian highlands by Hurni (1985). The extreme soil erosion rate (156.8 ± 285 t ha⁻¹ year⁻¹) for 2.9% of the catchment far exceeds the maximum tolerable soil erosion limit. The mean estimated soil erosion rate from the catchment is comparable with the similar studies that have been conducted for other catchments in Lake Tana basin, Blue Nile basin, and Ethiopia. For example, soil losses were estimated to be 39.8 t ha⁻¹ year^{−1 (}Estifanos 2014), 42.67 t ha⁻¹ year^{−1 (}Mengie et al. 2019), 49 t ha⁻¹ year⁻¹ (Mengesha et al. 2018), 42–47.4 t ha⁻¹ year⁻¹ (Tegegne and Biniam 2017; Gelagay and Minale 2016), and 37.9 t ha⁻¹ year⁻¹ (Asnake and Amare 2019); 27.5 t ha⁻¹ year⁻¹ (Haregeweyn et al. 2017) from catchments Ribb, Gumara, Dembecha, Koga, Beshillo and Upper Blue Nile basin, respectively. However, the mean rate is comparatively low to the estimated soil loss of 100 t ha⁻¹ y⁻¹ from all highlands of Ethiopia (FAO 1986). The comparison study illustrates that our results are in line with the previous works that have been carried out in neighboring and other watersheds of Ethiopia.

Of the five sub-catchments in the present study, Megech has the highest amount (49.9%) of the total estimated soil erosion and highest amount (32.7%) of the cells with extreme erosion potential followed by Gumero (19.1% and 23.3%), Garno (14.5% and19.4%), Derma (12.2% and 14.9%) and Gabi Kura (4.3% and 9.8%). It was also observed that the extent and magnitude of soil erosion were spatially varying. The northern and eastern parts of the catchment as well as the areas of the river banks have extreme erosion potential values whereas, the western parts of the catchment do not show such high values. The climate, altitude, topography, drainage density, soil conditions, soil types, geology, sediment delivery, sizes of sub-catchments, and land uses are markedly different in these parts. The differences in these parameters lead to the spatial and temporal variability of soil erosion potential. The main USLE factors that influence the soil erosion risk in the catchment are the slope length and steepness factor (LS-factor), the support practice factor (P-factor), as well as the land cover factor (C-factor). The land cover depends on climatic factors (e.g., rainfall) and soil conditions (texture, organic matter, structure, permeability, and topography). There is a close relationship between the soil erodibility and topography indicates the susceptibility of the catchment to erosion. The susceptibility also depends on the land cover, sensitivity, rainfall erosivity and seasonal variations. As the K and R factors, that depend on the natural soil types and the rainfall, cannot be altered (Stone and Hilborn, 2012), interventional catchment management strategies to prevent and control soil erosion shall focus on the LS, P, and C factors.

The relationship between rainfall and altitude has no significant difference in the catchment due to the complex and highly varied topography. As the topography of the catchment has an imposing effect on the rainfall distribution, it is important to consider the topography factor in soil erosion studies to evaluate the impact of land management practices in areas sensitive to land degradation. The altitude in the catchment does not directly influence the rainfall erosivity index, but the altitudinal gradient of precipitation does. The highly significant correlation between rainfall and other erosivity indices also shows that the high and extreme erosion risk in the study area was seasonal. The catchment has a unimodal rainfall pattern during June–September that peaks in July and August, and it shows that these critical months are more susceptible to erosion.

The rainfall has a good relationship with land cover. During the dry season (October-May) the rate of vegetation growth and agricultural activities are limited. Seasonal variations, land cover, drainage density, and the slopes played a vital role in soil erosion in the study area. According to the climatic seasonality, June, July, August, and September are the months with the highest risk of soil erosion in scattered barren lands and croplands. In May when usually rainy season starts and in October when the season ends in Ethiopia the soil erosion risk remains low but considerable in grasslands, plantation forests, temporarily abandoned agricultural lands, tilled lands on slopes, and near the river outlets of the catchment.

The high soil erosion risk in the highland area of the catchment is confirmed from a high-resolution satellite image and evidence of colluvial–alluvial formations deposited in the downstream banks of the lowland area of the catchment. These flooded areas with scarce land cover, sandy formations without consistency, and the undulating slopes produce morpho-dynamic conditions in the lowland area. However, once the intensity of the flood is reduced, these areas are being used for cultivation from October onwards. The relationship between rainfall and soil erodibility, drainage density, topography, and land cover can be seen as a key factor to understand the soil erosion in the lower parts of the catchment. Here, though the land cover is dependent primarily on the rainfall, the soil conditions also promote the growth of vegetation. The seasonal variation of the land cover such as natural forest, plantation forest, shrubs, seasonal grasses, and agriculture causes a change in the C factor for a specific month of erosion risk in the catchment. No agricultural activities are conducted in the croplands in the upper parts of the catchment during the dry season (December–April). The land is mostly devoid of any cover during the dry season. The first abundant rainfall that starts in May initiates vegetative growth. However, these initial rains cause severe erosion in highland parts of the catchment due to less land cover, improper tillage methods, and poor practices of conservation and management of land and terraces. Basic physical conditions in Ethiopia, which impact land degradation, include topography and rainfall variability from year to year and place to place, particularly in the drier parts of the highlands. The sequence of drier years with reduced vegetation cover followed by wetter years with heavy rainfall is conducive to high levels of soil loss (Berry et al. 2003). The higher C values in nearly half of the northern catchment indicate the occurrence of less vegetation cover and high estimated erosion. An extensive area of the catchment is on strong slopes; however, the areas with extreme erosion potential are mostly on moderately steep slopes.

Of all the land use types, most of the estimated erosion occurred in barren lands followed by agricultural lands and plantation forests. The relationship between topography and land use is another key factor to understand the soil erosion in the catchment. Very steep slopes were sparse in the catchment, comprising only 0.1% agricultural cells of extreme erosion potential. Moderate and intense steep slopes together comprised 34% of the entire catchment, and 52% of the cells with extreme erosion potential were located on those slopes. Of course, agriculture seemed to favor flat areas with 80% of agriculture occurring on 0-10% slopes and comprised 66% of the catchment. These flatter slopes were not exempt from extreme erosion because 48% of the cells with extreme erosion potential were located on them. Further, these slopes comprise 40% of the plantation forest and 40% of the barren land with 67% and 62% of the cells with extreme erosion potential, respectively.

The wet and dry densities in the study area indicate that the drainage density (Dd) of the catchment depends on the precipitation, infiltration capacity, underlying rock, soil texture, slope, altitude, vegetation, and hydraulic conductivity of the underlying soil. Higher Dd values indicate lower infiltration rates and higher surface flow velocity (Yalcin 2008) and are often related to a high sediment yield transported through the river network, high flood peaks, steep hills, and low suitability for agriculture. The largest area (55.8%) of the catchment is characterized by low Dd and 58.7% of cells with extreme erosion potential, whereas 44.2% of the catchment is characterized by high Dd and 41.3% of cells with extreme erosion potential.

A large percentage of the catchment, concentrated at low altitudes (1751–2000 m a.m.s.l), is characterized by 77.5% of the high Dd, strong slope gradient (3°–10°), and the largest (47.8%) of the area with extreme erosion potential. Such lower parts of the catchment area associated with the presence of several extremely narrow and scattered streams, and shallow depth of the Eutric vertisols and Luvisols. Thus, the fine clay soil and several streams promote soil loss in this area. Very less percent of the catchment occurred in high altitudes (2400–3100 m a.m.s.l) on moderately steep–steep slopes of > 10° and is characterized by high Dd and the least (12.8%) of the area with extreme erosion potential. It is associated with the rough soil surface and very shallow depth of the Leptosols and Nitosols. About 48.6% of the low Dd area in the catchment occurred at low elevations (1751–2000 m a.m.s.l) having only 8.8% of the cells with extreme erosion potential, whereas 51.4% of the low Dd area occurred at high elevations (2000–3014 m a.m.s.l) having 91% of the cells with extreme erosion potential.

The LS, C, and R factors are found in this study as the causative factors for high soil loss. At the river outlets, the high Dd is not related to high soil loss as the elevation and slope of the region is too low. In many areas of the northern catchment, the slopes are irregular and lead to runoffs into small drainage ways and flow at a velocity sufficient to detach and transport soil particles and thus further the hazard of soil erosion. Gully erosion occurs in the catchment due to highly erodible soils such as Vertisols and Leptosols, and misuse of soil and water resources. The plantation forests, barren lands, croplands, moderately steep and steep slopes, and low Dd areas were very susceptible to gully erosion that has been matched with the generated erosion risk map.

The high soil erosion in the catchment is associated with the highland fringe areas where the critical slopes exist. Field observations also showed the degradation and transformation of agricultural lands into wastelands. The presence of low clay content in loamy soils indicates that soil erodibility is high. In the upper part of the catchment, the drainage forms relatively steep narrow gorges that can attribute to small soil depth and high flow permeability that lower the drainage density and increase the surface runoff to lower parts of the catchment. The high soil loss in the upstream parts of the catchment is due to the shallow depth of the leptosols and the deeply weathered rocks, whereas, in the downstream parts of the catchment, the high soil loss results from the poor agricultural practices in croplands and high rainfall erosivity of uplands. The high percentage of barren lands in the hilly terrains, presence of leptosols, fallow lands, and rainfall erosivity of the uplands causes moderate to high amounts of soil loss in the downstream sections of the catchment. The rate of increase in the human population of Ethiopia (from 42 to 109 million between 1986 and 2019) continues to put a great demand for agricultural land. Such demand harms the natural resource conditions and ecosystems of grasslands, shrublands, and natural forests in the highlands as they are being converted to agricultural lands. The thematic change analysis helps to understand such land cover conversions with remarkable precision.

Conclusions

Soil erosion rate is estimated on a 1 ha grid scale for the study catchment and its sub-catchments based on the USLE model. The entire catchment has a moderate to high soil erosion severity and the mean annual soil erosion rate of the catchment is 37.89 ± 59.2 t ha⁻¹ year⁻¹. A major portion (77.1%) of the catchment area was classified under the slight erosion class, whereas the rest falls under the moderate to extreme erosion class. The area under high to extreme erosion potential with 10 to over 50 t ha⁻¹ year⁻¹ erosion severity is about 12.4% of the total area of the catchment (36,317 ha) that contribute about 82.1% of the total soil loss, which should be the highest priority for management efforts. The study revealed that Megech is a critical sub-catchment having the highest number of extreme and high erosion potential grid cells (47.1%), and thereby with the highest priority for implementing conservation practices, followed by Gumero (18.4%), Garno (17.2%), Derma (11.7%), and Gabi Kura (5.6%) sub-catchments. Multiple regression analysis results revealed that soil erosion was most sensitive to the topographic factor followed by the other USLE factors i.e., the support practice, soil erodibility, crop management, and rainfall erosivity. The study identified that spatial relationships between soil erosion rates and other factors such as slopes, land cover, land-use, drainage density, and soil erodibility were inconsistent. The highest intensity of soil erosion from the catchment was principally attributed to the steep slope and land-use changes. High soil erosion risk showed in the catchment during the rainy season only. However, initial rains caused severe erosion in highland parts of the catchment due to poor conservation practices and management of land. Barren land exhibited the highest soil erosion rates, followed by the croplands and plantation forest in the catchment. The results showed that gully erosion in all sub-catchments was attributed to higher steep slopes in land-use. This identification of the erosion-prone areas critical for effective catchment management helps to conserve soils. Spatial distribution of estimated erosion rates on grid cell basis and erosion severity classes coupled with various individual factors can help to understand the primary processes affecting erosion and recommend soil erosion prevention and control measures in extreme erosion-prone areas. This study also suggests that the applicability of the USLE model for the study catchment is to be evaluated with the sediment yield data collected at the outlets of the catchment.