Spatiotemporal modeling of soil loss and sediment yield due to land use and land cover changes: implications for watershed management

ASTER DEM

ASTER digital elevation model

GIS

geographic information systems

GPS

global position system

KII

key informant interview

LULC

land use land cover

MoWE

Ministry of Water and Energy

NMSA

National Meteorological Service Agency

RUSLE

revised universal soil loss equation

RS

remote sensing

SRTM

Shuttle Radar Topography Mission

SDR

sediment delivery ratio

SY

sediment yield

Soil erosion is a pervasive environmental challenge, threatening global sustainability and exacerbating food insecurity through the loss of fertile soil, nutrients, and organic matter. This challenge is amplified by increasing land-use modifications and subsequent land degradation, directly contributing to escalating rates of annual soil loss and diminished crop productivity (Abiye et al. 2023). Notably, water erosion stands as a primary driver of land degradation, heavily influenced by human activities and their impact on land utilization (Yirgu 2022). Globally, the scale of this issue is immense, with over 24 billion tons of fertile soil lost annually from cropland alone (AbdelRahman et al. 2023). Further compounding the problem, various disturbances affect soil in 33% of cropland, and 20% of irrigated land suffers from varying degrees of salinization. Indeed, water erosion is the leading cause of soil degradation worldwide, with approximately 80% of global agricultural land experiencing moderate to severe erosion (Dutta 2016).

Globally, growing human activity has had an enormous effect on natural ecosystems and the services they provide. Estimating the effects of land use and land cover (LULC) on ecosystem service values (ESVs) is crucial for determining how land use changes affect human well-being (Mersha et al. 2025). One of the biggest environmental issues facing the world today is soil erosion, which has both immediate and long-term implications (Moisa et al. 2021). Because nations with lower incomes are less equipped to handle soil erosion and replace lost nutrients, the economic effects are more severe. These nations are also experiencing rapid population growth, which encourages the extension of production to marginal and unstable territories and the increased usage of already limited resources. These activities worsen erosion and lower production, resulting in a population–poverty–land degradation cycle (Vlek et al. 2008).

In regions such as Africa, Asia, and Latin America, where the majority of the populace relies on agricultural activities, the issue of soil erosion presents a notable obstacle to the process of food cultivation (Gashaw et al. 2017). In contrast to the situations observed, soil erosion has primarily affected Africa, with over half of the one billion affected individuals globally living on the continent (Mandal et al. 2023).

Soil erosion in Ethiopia, particularly in response to LULC dynamics, has been the subject of numerous studies (Negese 2021; Shifaw et al. 2024; Reta Roba et al. 2025). A central concern within this body of research was the expansion of cultivated land at the expense of natural vegetation, and the consequential impact on hydrological processes and soil stability.

A study in Ethiopia demonstrated that the progressive conversion of forest land, shrub land, and grassland to cultivated areas over the past four decades has resulted in a significant increase in soil erosion rates, sediment yield, surface runoff, and mean wet monthly and annual stream flows (Negese 2021). Conversely, this LULC transformation has led to a reduction in dry average monthly flow, groundwater recharge and flow, and evapotranspiration (ET), indicating a disruption of the natural hydrological balance.

Supporting these findings, Reta Roba et al. (2025) documented a dramatic shift in LULC patterns, characterized by a substantial increase in cultivated land, from 62.3% to 77.0%, coupled with a sharp decline in forest cover, from 13.3% to 3.8%. This land cover transformation directly correlated with a significant rise in maximum soil loss, escalating from 726.7 tons/ha/year in 1993 to 937.8 t ha⁻¹ year⁻¹ in 2023. These findings reinforce the notion that agricultural expansion is a major driver of accelerated soil erosion in the region.

A consistent theme across these studies, as well as the research conducted by several authors (Abiye et al. 2023) in the Maybar watershed, emphasizes the disproportionate contribution of specific LULC types to soil loss. Bare land, cultivated land, and grazing land consistently emerge as significant sources of erosion, while forest land, shrub land, and water bodies exhibit minimal soil loss. This reinforces the understanding that land management practices and land cover composition are critical determinants of soil erosion rates in Ethiopian watersheds.

The relationship between LULC change and soil erosion severity is consistently emphasized in the literature. Studies by other authors (Gashaw et al. 2019; Kerbe et al. 2023) underscore the critical role of LULC change as a predictor of past, present, and future erosion patterns across varying spatial scales. These findings suggest that understanding and monitoring LULC dynamics is essential for effective soil erosion management and mitigation strategies.

Furthermore, the implications of soil erosion extend beyond environmental degradation, directly impacting agricultural productivity. The low average crop-yields in Ethiopia, when compared with international standards, are attributed, in part, to soil fertility decline resulting from topsoil removal through erosion (Shiferaw 2011; Gashaw et al. 2014). This connection between soil erosion and agricultural output highlights the socio-economic consequences of land degradation in the Ethiopian context.

The drivers of LULC change, such as rapid population growth, cultivation on steep slopes, clearing of vegetation, and overgrazing, are the main factors that accelerate soil erosion in Ethiopia (Vlek et al. 2008). As a result, the annual rate of soil erosion surpasses the annual rate of soil formation within the country (Belay & Mengistu 2021). Ethiopia experiences an annual loss exceeding 1.5 billion tons of topsoil from its highlands due to erosion, an amount that could have potentially contributed around 1.5 million tons of grain to the nation's harvest. This serves as a clear indication of the severe conditions by which soil erosion causes food insecurity, which makes urgent the intervention of effective soil and water management (Vlek et al. 2008).

The spatiotemporal assessment of soil erosion and sediment yield, particularly in relation to LULC changes, is essential for developing targeted conservation strategies and prioritizing critical sub-watersheds (Belay & Mengistu 2021). Research has consistently demonstrated the efficacy of integrating geospatial methodologies with the Revised Universal Soil Loss Equation (RUSLE) to quantify the impacts of LULC changes on soil erosion and sediment yield (Luvai et al. 2022; Gitima et al. 2023; Epple et al. 2025). For instance, Gashaw et al. (2019), in a study of the Andassa watershed within the upper Blue Nile basin, documented a significant expansion of cultivated land and built-up areas between 1985 and 2015, concurrent with a decline in natural vegetation cover. This LULC transformation resulted in a substantial increase in annual soil erosion rates, from 35.5 to 55 t ha⁻¹ year⁻¹, and sediment yield, from 14.8 to 22.1 t ha⁻¹ year⁻¹.

Understanding soil erosion is essential for the strategic development and prioritization of strategies within a watershed, as well as for comprehending the intricacies of the erosion phenomenon and its interconnected processes. The soil erosion assessment and mapping of areas susceptible to soil loss contribute significantly to the implementation of proper soil and water conservation practices and the effective management of the entire ecosystem within watersheds (Gelagay & Minale 2016). The mean yearly soil-loss information per unit of land area can be ascertained using the Universal Soil Loss Equation (USLE) and RUSLE (Kidane et al. 2019; Srinivasan et al. 2019).

The RUSLE involves fixed parameters such as topographical features (LS-factor) and soil characteristics (K-factor). These can remain constant over a long period. Controlled parameters such as rainfall (R-factor), cover management (C-factor), and land management (P-factor), can control the severity of soil erosion (Renard et al. 1997). However, the most-controlled component is cover management (C-factor) because it reflects the influence of vegetation type, density, cropping, and management practices on soil erosion (Kulimushi et al. 2021; Mariye et al. 2022).

Several studies used different methods to estimate soil loss and sediment movement in response to LULC dynamics. For instance, in Shekar & Mathew (2024), the empirical models utilized in their investigation were RUSLE to estimate soil loss, and variables such as LULC, rainfall, DEM, and soil were used in this study. The RUSLE model incorporates several key parameters that have been identified as significant contributors to soil loss, like the land-use management factor (C-factor), the length and steepness factor (LS-factor), the agricultural practice factor (P-factor), the soil erodibility factor (K-factor: units in tons ha h/MJ/ha/mm), and the rainfall erosivity component (R-factor: units in MJ mm/ha/h/year). The RUSLE equation was used also by Arias-Muñoz et al. (2023) to quantify soil erosion rates, and CA-Markov chain analysis was implemented to forecast the LULC change for the 2030 year. The InVEST model has been applied in several studies conducted in Ethiopia to analyze the impacts of LULC changes on soil erosion and sediment delivery in the Winike watershed, Omo Gibe Basin (Aneseyee et al. 2020).

The study demonstrated how to use geographic mapping techniques to identify flood impacts produced by river catchments. GIS technology in combination with remote sensing data can be used to study how floods affect river catchments in the Niger Delta. According to their study findings, the Niger and Forcardos river catchments rank ‘high,’ while the Orashi River catchment is the most prone to flooding due to its high runoff parameters and short concentration period and the Ikoli, Nun, and Bomadi River catchments, which indicate ‘medium.’ It is therefore advised that the Orashi, Niger, and Forcados river catchments be prioritized in order to reduce the impact of floods in the Niger Delta research region (Oborie & Rowland 2023).

The Brokena watershed, which is the primary focus area of this study, is situated in the lower regions of the Awash River basin and is known for its high agricultural productivity. However, the watershed has experienced human-induced land degradation, primarily attributed to agricultural expansion and unregulated LULC modifications. Therefore, in order to establish watershed management strategies and rehabilitate the degraded watershed, it is essential to evaluate the effects of LULC changes on soil erosion and sediment yield in this watershed.

While previous studies have employed RUSLE and InVEST SDR in Ethiopian watersheds, this research distinguishes itself by uniquely examining the Borkena watershed, a critical agricultural area experiencing rapid LULC changes. Furthermore, this study offers a more comprehensive understanding of the watershed's dynamics by focusing on the combined impact of cropland expansion and urban growth on soil erosion and sediment yield over a 30-year period. This extended temporal analysis, spanning three decades, provides a crucial long-term perspective on soil erosion trends within the watershed, extending beyond the limitations of studies focused solely on current or recent conditions. This study distinguishes itself from previous research on Ethiopian watersheds by focusing on the Borkena watershed, a region characterized by a unique confluence of rapid urbanization and intensive agricultural practices. Unlike many studies that primarily focus on the impact of agricultural expansion on soil erosion, this research examines the synergistic effects of both urban growth and agricultural intensification. The Borkena watershed is experiencing a particularly rapid expansion of built-up areas, coupled with the intensification of cropland cultivation on steep slopes, which presents a distinct challenge to soil stability (Mersha et al. 2024). Furthermore, the watershed's varied topography and specific soil types, which exhibit a higher susceptibility to erosion than in other regions, contribute to a unique environmental context. By quantifying the spatiotemporal dynamics of soil loss and sediment yield in response to these combined pressures, this study provides novel insights into the complex interactions between land use changes and soil degradation in a rapidly evolving agricultural landscape. This analysis is crucial for developing agroclimatic-specific watershed management strategies that address the specific challenges posed by the Borkena watershed's unique characteristics.

Therefore, the establishment of strategies for managing watersheds and the restoration of degraded watersheds rely heavily on a comprehensive understanding of the effects of LULC modifications on soil erosion and sediment yield in the target watershed. As a result, this study examines the phenomenon of soil erosion and sediment yield, considering both spatial and temporal dimensions resulting from LULC changes. The specific objectives of this study were as follows: (1) to map the trends of LULC change across time; (2) to estimate the spatiotemporal trend of soil erosion risk in association with the LULC change; and (3) to estimate the spatiotemporal trend of sediment yield at the watershed linked with LULC change.

The Borkena River originates in Kutaber Woreda, located at the confluence of the Abay and Awash River basins. From its source, the river flows southeastward, converging with the Jara River and traversing the Cheffa Swamp, before ultimately joining the Awash River. The river's headwaters are situated within the steep chains and escarpments of the northern plateau, proximate to the Afar Rift. A significant tributary, the Berberie River, arises from the hilly terrain northeast of Kombolcha town, meandering through the town before its confluence with the Borkena River upstream of the Addis Ababa–Dessie highway bridge.

The watershed's topography varies significantly, transitioning from flat plains to steep slopes. An analysis of the digital elevation model (DEM), categorized according to FAO slope classes, revealed the following distribution: 34% of the watershed exhibits gentle slopes (0°–7°), 14% exhibits slightly undulating slopes (7°–17°), and 17% exhibits moderately steep slopes (17.6°–26.8°). The remaining 15% of the land area is classified as steep or very steep, with slopes exceeding 26° (Table 1).

Landsat imagery with a 30 m resolution and minimal cloud cover was acquired from the United States Geological Survey (USGS) EarthExplorer. Specifically, Landsat 5 Thematic Mapper (TM) data were utilized for the years 1993 and 2003, while Landsat 8 Operational Land Imager/Thermal Infrared Sensor (OLI/TIRS) data were employed for 2013 and 2023. Although the Landsat 8 OLI/TIRS data has a 15 m panchromatic band, for consistency, the 30 m resolution was used for all datasets. All imagery was selected from the same seasonal period, specifically March or April, across the respective study years. Detailed specifications of the satellite data are presented in Table 2. Ground validation was conducted using reference samples derived from Google Earth Pro for the years 1993 to 2023, supplemented by a field survey in 2023 to ensure accuracy.

Image preprocessing was conducted using ArcGIS 10.8 software. The catchment boundary was utilized to delineate the study area. A comprehensive field survey was undertaken between May and July 2023, employing topographic maps, Google Earth, and preliminary LULC classification maps derived from satellite imagery for the respective reference years. This field validation, providing direct observational data of the study site, was integrated with a supervised digital image classification technique.

Multi-temporal Landsat satellite imagery, spanning four distinct dates, was analyzed to assess LULC dynamics using the supervised classification method. A combined methodological approach was developed to facilitate the interpretation, analysis, mapping, and quantification of the collected data.

LULC datasets for the study watershed, covering the period from 1993 to 2023, along with their associated drivers, were generated. Supervised classification, as detailed in Table 2, was employed to derive LULC maps from Landsat 5 TM (1993), Landsat 7 ETM + (2003), and Landsat 8 OLI/TIRS (2013 and 2023) imagery.

For the supervised classification using the maximum likelihood classification (MLC) method, approximately 450 ground-truth data points were collected for each of the LULC classes: cropland, forest, shrub land, water body, bare land, and built-up area. To classify the Landsat imagery from 1993 to 2023, reference data were derived from Google Earth imagery corresponding to the respective time-periods. Additionally, focus-group discussions were conducted with local elders to supplement the reference data and provide contextual information. The classified imagery was subsequently evaluated against the Google Earth imagery to ensure accuracy. ArcGIS 10.8 software was used for all mapping procedures. Monthly rainfall data, covering the period from 1990 to 2020, were obtained from the Ethiopian Meteorological Institute (Table 3).

The quantitative estimation of soil erosion was computed four times in 1993, 2003, 2013, and 2023. The K, LS, and P factors were kept constant, whereas the R and C factors were periodically changed as the LULC changed (Kulimushi et al. 2021). Five classes of severity were used to reclassify erosion hazard maps: very low (0–5 t/ha/year), low (5–15 t/ha/year), moderate (15–30 t/ha/year), high (30–50 t/ha/year), and very high (>50 t/ha/year), adapted from Kulimushi et al. (2021).

The study also considered topsoil organic carbon (%), soil types (adapted from the revised FAO-UNESCO-ISRIC (1990) soil map legend), and soil texture classes. The soil erodibility factor (K-factor) is expressed in tons ha⁻¹ MJ⁻¹ mm⁻¹. The K-factor quantifies the inherent susceptibility of soil particles to detachment and transport by rainfall (Buraka et al. 2022). This factor is strongly influenced by soil properties and plays a critical role in the development of effective soil-conservation strategies (Shabani et al. 2014) representing soil erosion per unit rainfall erosivity. In this study, a soil map was acquired from the Ministry of Water and Energy (MoWE), and corresponding K-factor values were assigned based on established literature (Table 4). Field surveys were also conducted to verify soil colors and supplement the soil map data.

The K-factor was estimated using the USDA soil texture classification and organic matter percentage (%SOM). The observed K-factor value given in US units is converted to International System units (SI units) (tons ha h ha⁻¹ MJ⁻¹ mm⁻¹) by multiplying by 0.132.

The cover management factor (C-factor) quantifies the influence of cropping and crop management practices on soil erosion rates in agricultural lands (Tsegaye & Bharti 2021; Dapin & Ella 2023; Kerbe et al. 2023). Due to its inherent temporal variability, the C-factor is often considered a critical predictor of past, present, and future soil erosion severity (Gashaw et al. 2017). This dimensionless factor ranges from approximately 0, representing well-protected land cover, to 1, representing barren soil prior to vegetation establishment (Kulimushi et al. 2021).

The six land-uses and land covers were identified with overall classification accuracy using ArcGIS Imagine 10.8. Hence, four C-factor raster grids were derived for 1993, 2003, 2013, and 2023 to assess the impact of LULC change on soil erosion. Descriptions of the land cover categories of the study watershed are shown in Table 5.

The support practice factor (P-factor) represents the influence of conservation practices, such as contouring, terracing, strip cropping, and tillage, on soil erosion rates, particularly in relation to tillage direction (Kulimushi et al. 2021). The implementation of various conservation methods significantly affects this factor (Belayneh et al. 2019). In this study, a fixed P-factor value of 0.75 was adopted, signifying the presence of unsustainable conservation measures, as suggested by Karamage et al. (2016).

The proportion of gross erosion delivered to the outlet of the watershed area in a specific timeframe is known as the sediment delivery ratio, or SDR (Kidane et al. 2019). The SDR may be estimated using a variety of methods, but they can all lead to the same result. The SDR is dependent on several variables, including watershed size, channel slope, stream density, LULC, sediment supplies, and rainfall–runoff (Kulimushi et al. 2021; Moisa et al. 2021). Topographical features such as stream height ratio and stream length influence sediment transportation within the watershed (Kulimushi et al. 2021).

The InVEST model has been applied in several studies conducted in Ethiopia to analyze the impacts of LULC changes on soil erosion and sediment delivery in the Winike watershed, Omo Gibe Basin (Aneseyee et al. 2020). The InVEST SDR model (Gashaw et al. 2021) (i.e., InVEST 3.9.0), which is a spatial explicit model working at the cell size of the DEM, was applied to estimate annual soil loss and sediment export in the Borkena watershed. InVEST SDR first computes annual soil loss using the soil loss algorithm (Renard et al. 1997).

SY indicates sediment yield (t ha⁻¹ year⁻¹), SDR indicates SDR, and A indicates annual mean soil loss (t ha⁻¹ year⁻¹). The sediment delivery distributed (SEDD) model provides estimates about the distributions of sediment yield in a certain area using two major components such as soil erosion and sediment delivery ratio (SDR), as shown in Equation (6) (Batista et al. 2017). SEDD partitions a certain area into morphological units and determines the SDR for each unit of a watershed (Ferro & Porto 2000; Yan et al. 2018). The SDR expresses the probability that eroded particles on a given upland location will reach the nearest stream channel (Ferro & Minacapilli 1995), and ultimately the outlet of the watershed (Batista et al. 2017). Hence, the SEDD model is imperative to estimate the amount of sediment yield from each unit of a watershed.

The LULC classification for the years 1993, 2003, 2013, and 2023 demonstrated overall accuracies of 87.40%, 93.11%, 97.11%, and 87.01%, respectively, with corresponding Kappa coefficients of 0.79, 0.89, 0.95, and 0.77 (Table 6). The comprehensive accuracy assessment methodology, encompassing a confusion matrix and independent validation points, is detailed in Mersha et al. (2024). According to Monserud (1990), Kappa values within the range of 0.70–0.85 are considered excellent, indicating a substantial agreement between the classified imagery and ground-truth data. Analysis of LULC patterns revealed a consistent expansion of cropland, increasing from 72,666 ha in 1993 to 82,486 ha in 2023 (Table 7). Similarly, the built-up area expanded throughout the study period. Conversely, forest and shrub land areas decreased (Table 7). The expansion of built-up areas primarily occurred at the expense of agricultural land and natural vegetation. Additionally, vegetation cover was transformed into built-up areas, ridges, and agricultural land, as noted by Naikoo et al. (2020).

Analysis of the time-series mean monthly rainfall data revealed a range from 968 to 1,098 mm/year, with an overall mean of 1,034 mm/year. Notably, the trend in mean annual rainfall showed a slight increase over the study period. Specifically, the mean annual rainfall for the selected years was 1,025 mm in 1993, 980.06 mm in 2003, 1043.336 mm in 2013, and 1186.08 mm in 2021 (Figure 6). These values were subsequently used to calculate the R-factor, reflecting the erosive power of rainfall within the watershed.

The mean annual rainfall for the study area was calculated for four distinct periods: 1025.8 mm in 1993, 1089.34 mm in 2003, 977.0 mm in 2013, and 1037.07 mm in 2023. Using the Hans Hurni method, specifically the regression equation developed for Ethiopian conditions as detailed in Hurni (1985), the rainfall erosivity (R-factor) was computed for each corresponding year. The resulting R-factor values were 566.33, 601.94, 539.15, and 572.31 MJ mm ha⁻¹ h⁻¹ year⁻¹, respectively. These values are further illustrated in Figure 4 and summarized in Table 9.

To determine the rainfall erosivity (R-factor) across the study period (1993–2023), the mean annual precipitation for each respective year was utilized. The resulting long-term erosivity values for the watershed ranged from 566 to 602 MJ mm ha⁻¹ h⁻¹ year⁻¹.

Notably, the Kutaber District's highland sub-watershed exhibited the highest recorded erosivity rating. This observation suggests that this specific terrain receives significantly higher rainfall compared with downstream areas. Consequently, the watershed displayed considerable variability in R-factor values. Such fluctuations in rainfall erosivity are expected to influence the estimated soil erosion rates.

Consistent with these findings, Kidane et al. (2019) demonstrated that in Ethiopian watersheds, highland regions, characterized by higher rainfall and steeper slopes, contribute the greatest sediment yield. This reinforces the understanding that areas with elevated erosivity, as observed in the Kutaber highlands, are likely to experience increased soil erosion and sediment transport.

The K-factor, representing soil erodibility (susceptibility to detachment by rainfall), ranged from 0.018 to 0.0203. Specific values for the studied soil types were: chromic vertisols (0.0185), eutric cambisols (0.0203), eutric leptosols (0.0183), and eutric regosols (0.0201). These values, shown in Figures 5(a) and 5(b), fall within expected ranges.

The K-factor, expressed in tons ha⁻¹ MJ⁻¹ mm⁻¹, was also influenced by topsoil organic carbon (%), soil types (adapted from the revised FAO-UNESCO-ISRIC (1990) soil map legend), and the soil nomograph. The soil nomograph is a commonly used method that considers the relative proportions of soil texture, permeability, soil structure, and organic matter content (Wischmeier & Smith 1978; Fu et al. 2006; Addis & Klik 2015).

Several significant soil categories that defined the Borkena watershed were identified. Most mapping units have distinct soil qualities because of the range of landforms within the watershed. There are four distinct types of soil in the research area: chromic cambisols, eutric cambisol, eutric leptosols, and eutric regosols. Among these, eutric cambisol (53%) and eutric regosols (20%) were more dominant in the study area (Table 10).

The C-factor takes into consideration how plant cover prevents erosion from water (Belayneh et al. 2019). There is more soil erosion by water in places with little or no vegetation. Conversely, because vegetation better shields the soil surface from erosion, soil erosion is less common on vegetated ground (Kulimushi et al. 2021). Thus, by changing the LULC types into more vegetated surface coverings, soil erosion may be decreased. As a result, the C component may be the most crucial in lowering soil erosion. Creating C-factor watershed maps for the research (Figure 8) from the corresponding LULC temporal layers, C factors were assigned for each LULC type based on the literature (Table 7).

An approach to combating the widespread incidence of soil erosion is to implement proper conservation measures (Gashaw et al. 2017). The P-factor represents the role of conservation practices in reducing erosion (Adediji et al. 2010). The P-factor has a value between 0 and 1. Generally speaking, places without conservation procedures are given a value of 1, whereas regions with effective conservation practices are given minimum values that are close to 0. Therefore, the lower the P value, the more effective the practice of protecting against erosion (Panagos et al. 2015). There are very few soil conservation practices in the study watershed that are present in the various farm plots. Because there are no major conservation practices followed in the study area, 1 was assigned to the P-factor, as suggested by Panagos et al. (2015), which does not affect the calculation (Figure 9).

The mean annual soil erosion rate increased from approximately 25.5 t ha^{− 1} year⁻¹ in 1993 to 35.5 t ha^{− 1} year⁻¹ in 2003 and to 42.3 t ha^{− 1} year^{− 1} in 2023 as a result of the LULC changes that occurred between 1993 and 2023 (Table 12). After evaluating the soil erosion rates of each kind of LULC, it was found that cropland showed the highest mean annual soil erosion rate, measuring 51.6 t ha^{− 1} year⁻¹ in 1993 and 65.5 t ha^{− 1} year⁻¹ in 2023 (Table 12). On the other hand, the average yearly rate of soil erosion calculated from cropland and the entire watershed is higher than the rate calculated from shrub land, water body, built-up areas, and bare land. As a result, from around 1,896,080 tons in 1993 to 3,232,907 tons in 2013 and 3,469,077 tons in 2023, the watershed's total soil loss increased (Table 13). Cropland accounted for almost 94% of the total soil loss during the periods (Table 13). Less than 10% and 3.4% of the total soil loss was attributed to the shrub land and bare land, respectively, and the contributions of the forest, water body, and built-up area were negligible.

Soil erosion within the watershed varied dramatically. Steep terrains experienced an average annual loss of 144 t ha⁻¹ year⁻¹, while plains showed negligible erosion (0 t ha⁻¹ year⁻¹) (Figure 10). The watershed's total annual soil loss was estimated at 0.097 million tons. In 1993, 78% of the watershed was categorized as low risk for soil erosion (0–11 t ha⁻¹ year⁻¹), with the remaining areas classified as moderate (11%), high (8%), very high (2%), and severely affected (1%). By 2023, the low-risk area decreased to 62%, while moderate, high, very high, and severely affected areas increased to 7%, 21%, 7%, and 3%, respectively (Table 11). The watershed's average soil loss rate increased from 25.5 t ha⁻¹ year⁻¹ in 1993 to 42.3 t ha⁻¹ year⁻¹ in 2023, with intermediate rates of 33.5 and 35.5 t ha⁻¹ year⁻¹ in 2003 and 2013 (Table 12). These rates are comparable to the 18 t ha⁻¹ year⁻¹ reported by Hurni (1998) for the Ethiopian highlands.

The current result also agrees with similar findings. For instance, The central and northern highlands had an average soil erosion rate of 35 t ha⁻¹ year⁻¹, according to FAO (1986), and Bewket & Teferi (2009) revealed that the Chemoga watershed in the northwest highlands of the Blue Nile basin had an average soil erosion rate of 93 t ha⁻¹ year⁻¹. In the Koga watershed of the Blue Nile basin, a recent study by Gelagay & Minale (2016) found that the average rate of soil erosion was 47.4 t ha⁻¹ year⁻¹. A study conducted in the Central Rift Valley watershed by Mosello et al. (2015) revealed a comparatively high rate of soil erosion, exceeding 45 t ha⁻¹ year⁻¹. Given that 44% of the watershed is made up of undulating plains with a slope of less than 12%, the topography of the area may be the reason for the comparatively low average soil erosion rate reported in the currently examined watershed.

The Borkena watershed was classified into five soil erosion priority classes (Table 11, Figure 10). Notably, 22% of the watershed (27,834 ha), exceeding the tolerable erosion limit of 11 t ha⁻¹ year⁻¹, contributed 63% of the total soil loss (78,270 ha). It can be observed from the assigned classes that the different priority areas contributed differently to the total erosion rate. For instance, the priority class severe covers only 9% of the entire watershed, which amounts to 12,787 ha. Priority classes II and III combined cover approximately 29% of the watershed, amounting to 35,198 ha. Similar to the findings of this study, other studies have indicated that small areas of the watershed contributed to a significant amount of soil loss. For example, in the Borena watershed, places with very high to extremely severe soil loss accounted for 29.85% of the total predicted soil loss, but the same areas also provided 60.03% of the total soil loss (Abate 2011). In the Wondo Genet catchment, 23.5% of the watershed was responsible for 54.54% of the soil loss (Gelagay & Minale 2016). On the other hand, there is a little higher quantity of soil loss in priority class V. Despite accounting for 70.65% of the yearly total soil loss, the rate of soil loss in these locations is still within the allowable limit. Therefore, prioritizing soil conservation efforts in the higher priority classes, as recommended by Abate & Singh (2011), is crucial for maximizing impact.

The continuous increase of the mean SDR in the study area is generally due to the expansion of non-vegetated LULC types and the dwindling of the vegetated class. On the other hand, the higher SDR obtained in cultivated compared with vegetated LULC covers in the watershed study is reasonable because the probability of eroded materials reaching the nearest stream channel from the non-vegetated covers is higher than that from the vegetated LULC (Tesfaye et al. 2014).

The annual sediment yield was approximately 801,760 tons in 1993, 1,444,795 tons in 2013, and 1,448,964 tons in 2023, in accordance with the rise in soil erosion and SDR. During all of these times, over 95% of the total sediment output came from cropland. The shrub land (less than 5%) is the second factor contributing to these increases. However, Table 17 indicates that the contributions of built-up areas, water bodies, and forests are not very significant.

Studies in the upper Blue Nile basin, such as those in the Andassa watershed (Gashaw et al. 2021), have consistently shown cultivated land to be the primary contributor to maximum sediment yield. This aligns with the findings of the current study, where cultivated land, covering over 65% of the watershed, exhibits the highest sediment yield. This dominance is attributed to its extensive coverage and its inherent vulnerability to erosion (Gashaw et al. 2019). Furthermore, cultivated land demonstrates a distinct pattern of extremely high severity soil erosion (exceeding 50 t ha⁻¹yr⁻¹) with an increasing trend, setting it apart from other land use/land cover (LULC) classes (Kulimushi et al. 2021). Because cultivated lands are subject to greater erosion and flow velocity, which in turn contribute to higher sediment yield, the highest predicted sediment yield from cultivated land and the lowest from forest area in the study watershed are generally reasonable.

Similar results were observed in the Maybar watershed, Awash basin, Ethiopia (Abiye et al. 2023) where an analysis of soil erosion response to LULC changes revealed that bare land, cultivated land, and grazing land exhibited the highest annual soil loss, while forest land, shrub land, and water bodies showed the lowest. In the current study area, the most severely soil-erosion-affected regions were located in steep slope areas, whereas gentle slopes experienced minimal erosion. This disparity in sediment yield intensity, relative to soil erosion, can be explained by the dominance of deposition over erosion in gentle slope areas, leading to measurable sediment yield, albeit at very low intensities. Conversely, steep slope areas, where erosion surpasses deposition, exhibit a reduction in sediment yield intensity compared with soil erosion, particularly at the most severe levels.

This study utilized rainfall, soil, LULC and DEM data. The RUSLE was implemented within a GIS environment to estimate soil loss. This involved calculating the R, K, C, P and LS factors. The analysis revealed substantial LULC changes over the past three decades, leading to significant soil erosion, particularly in the steep areas of the watershed. Soil loss reached a maximum of 144 t ha⁻¹ year⁻¹, with cultivated land being the primary contributor. Due to resource constraints, immediate implementation of comprehensive soil conservation measures was not feasible. Therefore, the watershed was divided into five priority areas to guide targeted interventions. The mean annual soil erosion rate increased significantly from 25.5 t ha⁻¹ year⁻¹ in 1993 to 42.3 t ha⁻¹ year⁻¹ in 2023. Concurrently, the SDR and sediment yield showed similar upward trends, rising from 0.41 t ha⁻¹ year⁻¹ to 0.49 and 18.3 to 29.4 t ha⁻¹ year⁻¹, respectively, indicating a high erosion risk. The expansion of cultivated land was identified as the principal driver of these increases.

The anticipated rise in soil erosion poses a significant threat to land productivity, directly impacting the future livelihoods of local communities. Furthermore, the expected increase in sediment yield will likely exacerbate sedimentation at small-scale irrigation sites within the watershed. To mitigate these adverse effects, diversifying livelihoods is crucial to reduce agricultural expansion into marginal lands, thereby minimizing subsequent increases in soil loss and sediment yield entering water bodies (Renard et al. 1997; Aneseyee et al. 2020; Shakir 2022; Abiye et al. 2023). Future research should focus on detecting and mapping the potential regions for afforestation of the high-risk areas.

Based on the findings of this study, we recommend the following: (1) Implement targeted afforestation programs in high-erosion risk areas. (2) Promote sustainable agricultural practices, such as terracing and contour plowing. (3) Enforce stricter regulations on land use changes, particularly regarding cropland expansion and urban development. (4) Conduct regular monitoring of soil erosion and sediment yield to assess the effectiveness of conservation measures. (5) Increase community involvement in watershed management through education and awareness programs.

The authors would like to thank the Ethiopian National Meteorological Institute and the Ethiopian Ministry of Water and Energy for providing us with the meteorological and daily rainfall data, respectively.

E.M. conceptualized the whole article, S.D. developed the methodology, arranged the software, rendered support in data collection; W.M. did formal analysis and validated the data, and E.M. wrote the original draft. S.D. supervised the article, investigated the data, visualized the process, M.A. reviewed the article. All authors have read and agreed to the published version of the manuscript.

The authors received no financial support for the research authorship and/or publication of this article

Data cannot be made publicly available; readers should contact the corresponding author for details.

The authors declare there is no conflict.