Introduction
Soil erosion represents one of the most critical environmental challenges globally, threatening agricultural sustainability, ecosystem health, and rural livelihoods (Montgomery 2007, Panagos et al. 2018, Borrelli et al. 2020). The degradation of soil resources is particularly severe in developing countries characterised by mountainous terrain, high rainfall erosivity, and agricultural dependence (Hurni et al. 2015, Nyssen et al. 2015). In Ethiopia, where agriculture forms the backbone of the economy and supports over 80% of the population, soil erosion poses an existential threat to food security and sustainable development (Haregeweyn et al. 2015, 2017).
The Ethiopian highlands, which cover approximately 45% of the national landmass, are among the most erosion-prone regions globally owing to steep slopes, intense seasonal rainfall, and historical land degradation (Haregeweyn et al. 2015, Nyssen et al. 2015). Annual soil loss rates in these highlands have been estimated to be between 16 and 300 Mg·ha−1·yr−1, significantly exceeding the soil formation rate of 1–2 Mg·ha−1·yr−1 (Hurni 1983, Tadesse, Abebe 2014). The economic impact of soil erosion in Ethiopia is substantial, with estimates suggesting annual costs equivalent to 2–3% of the national GDP through reduced agricultural output, sedimentation of reservoirs, and loss of ecosystem services (Shiferaw, Holden 2001, Tamene, Vlek 2008). Recognising this crisis, Ethiopia has implemented various soil and water conservation measures over recent decades. These include physical conservation structures, such as stone bunds, soil bunds, bench terraces, and check dams, as well as biological interventions including afforestation, area exclosures, and agroforestry systems (Haregeweyn et al. 2015, Adimassu et al. 2017). The sustainable land management programme (SLMP), initiated in 2008 and expanded in subsequent phases, has been particularly influential in coordinating watershed rehabilitation efforts across the country (Wolancho 2015). The productive safety net programme (PSNP) has further contributed through community mobilisation for soil and water conservation works (Adimassu et al. 2017). While isolated studies have documented localised improvements following these interventions (Nyssen et al. 2004, Adimassu et al. 2014), comprehensive assessments of soil erosion dynamics at the watershed scale over extended periods remain limited. Several models have been developed to assess soil erosion, ranging from empirical approaches such as the universal soil loss equation (USLE) and its revised version revised universal soil loss equation (RUSLE) to process-based models including WEPP and EUROSEM (Merritt et al. 2003, de Vente et al. 2013). Among these, RUSLE has been widely adopted owing to its relatively simple structure, moderate data requirements, and adaptability to diverse landscapes (Wischmeier, Smith 1978, Renard et al. 1997). The integration of RUSLE with geographic information systems (GIS) and remote sensing data has further enhanced its utility for spatial soil erosion assessment across large areas (Prasannakumar et al. 2012, Ganasri, Ramesh 2016).
Previous applications of RUSLE in Ethiopia have primarily focused on single time point assessments of erosion risk (Bewket, Teferi 2009, Gelagay, Minale 2016, Tamene et al. 2017). These studies have consistently identified areas of high erosion risk but have not adequately captured the temporal dynamics of soil erosion in response to evolving land cover and management practices. Moreover, many existing studies have not clearly distinguished between the factors that control the spatial pattern of erosion and those that drive temporal changes, leading to potentially circular interpretations when only the cover management factor varies between assessment periods. Understanding these dynamics and their drivers is crucial for evaluating the effectiveness of past conservation investments and guiding future interventions (Nyssen et al. 2004, Haregeweyn et al. 2015). The Muger River Watershed in the central Ethiopian highlands represents an appropriate case study for examining temporal soil erosion dynamics. The watershed has experienced signifcant land use and land cover changes over recent decades and has been a target of various soil and water conservation initiatives (Dibaba et al. 2020). Despite this, no comprehensive assessment of erosion patterns spanning more than two decades has been conducted in this watershed.
This study aims to address this gap by quantifying and analysing soil erosion dynamics in the Muger River Watershed over 24 years (2001–2025) using the RUSLE model integrated with GIS and cloud computing platforms. The specific objectives are to: (1) assess the spatial distribution of soil erosion risk in 2001 and 2025; (2) quantify changes in erosion rates and risk patterns over the study period; (3) analyse the relationship between topographic factors, land cover change, and soil erosion dynamics while distinguishing between spatial and temporal drivers; and (4) evaluate the implications of these fndings for future soil conservation efforts and policy development in the Ethiopian highlands.
Materials and methods
Study area
The Muger River Watershed is situated in the central Ethiopian highlands within the Blue Nile (Abay) basin, covering an area of approximately 7347 km2 (Fig. 1). The watershed boundary was delineated using the HydroBASINS Level 07 dataset (Lehner, Grill 2013), which provides standardised subbasin boundaries derived from high-resolution hydrological data. The watershed is characterised by diverse topography, with elevations ranging from 974 m to 3538 m above sea level and a mean elevation of 2308 m. The terrain is predominantly hilly to mountainous, with a mean slope of 11.0° and approximately 24.5% of the watershed having slopes >15°.
Fig. 1.
Location and topographic characteristics of the Muger River Watershed in the central Ethiopian highlands. The watershed covers approximately 7347 km2 within the Blue Nile (Abay) basin, spanning elevations from 974 to 3538 m above sea level. The map displays the watershed boundary derived from HydroBASINS Level 07, drainage networks, and elevation distribution. Steep terrain (>15-degree slope) comprises approximately 24.5% of the watershed area.
The climate is classified as subhumid tropical with distinct wet (June–September) and dry (October-May) seasons. Mean annual precipitation across the watershed ranges from approximately 1037–1429 mm, with a watershed average of 1205 mm based on Climate Hazards Group InfraRed Precipitation with Station data (CHIRPS) for the period 2001–2025. Mean annual temperature varies between 15°C and 20°C, decreasing with elevation. Rainfall is typically intense and seasonal, concentrating signifcant erosion potential during the wet months (Dibaba et al. 2020).
Soils in the watershed are predominantly Nitisols, Vertisols, and Cambisols under the FAO classifcation system, with varying degrees of degradation (Hurni et al. 2015). The watershed exhibits considerable soil textural diversity, with sand content ranging from 16% to 46% (mean 32.7%), clay content from 26% to 62% (mean 40.3%), silt content from 19% to 39% (mean 27.0%), and soil organic carbon from 1 to 17 g·kg−1 (mean 4.6 g·kg−1), as derived from the OpenLandMap soil database (Hengl et al. 2017). Natural vegetation consists primarily of dry Afromontane Forest and woodland, though much has been converted to agricultural land over recent decades. Agriculture is the dominant economic activity, with smallholder mixed crop-livestock farming systems prevailing throughout the watershed. Major crops include teff, wheat, barley, maize, and pulses. The watershed has been the site of various soil and water conservation initiatives since the early 2000s, including the SLMP and community-led watershed management efforts (Adimassu et al. 2014). Conservation interventions have included construction of stone bunds, soil bunds, bench terraces, and check dams, as well as establishment of area exclosures and tree planting activities, implemented with varying intensity across the watershed (Gebremichael et al. 2005, Haregeweyn et al. 2015). Stone bunds and soil bunds represent the most widely adopted physical conservation structures in the Ethiopian highlands, functioning by reducing slope length and decreasing overland flow velocity (Amsalu, de Graaff 2007, Adimassu et al. 2014). Area exclosures, where degraded lands are protected from grazing and cultivation to allow natural vegetation regeneration, have also been increasingly implemented across the watershed (Descheemaeker et al. 2006).
RUSLE model
The RUSLE was employed to estimate soil erosion rates. RUSLE is an empirical model that predicts long-term average annual soil loss resulting from sheet and rill erosion (Renard et al. 1997). The model calculates soil erosion as a product of five factors, as expressed in Eq. (1):
where:– A is the computed soil loss (Mg·ha−1·yr−1),
– R is the rainfall runoff erosivity factor (MJ·m-m·ha−1·h−1·yr−1),
– K is the soil erodibility factor (Mg·h·M-J−1·mm−1),
– LS is the slope length and steepness factor (dimensionless),
– C is the cover management factor (dimensionless), and
– P is the conservation support practice factor (dimensionless).
Data sources
Multiple datasets were utilised to derive the five RUSLE factors for both 2001 and 2025 (Table 1). Topographic data were obtained from the Shuttle Radar Topography Mission (SRTM) digital elevation model (DEM) at 30 m spatial resolution (Farr et al. 2007). Rainfall data were sourced from the CHIRPS, which provides daily precipitation estimates at 0.05-degree (~5 km) resolution (Funk et al. 2015). Mean annual precipitation was computed from daily CHIRPS data spanning the full study period (2001–2025) to ensure consistency between the rainfall erosivity estimates and the temporal scope of the erosion analysis.
Table 1.
Data sources used for RUSLE model parameterisation.
| Data type | Source | Resolution | Period/year | Used for |
|---|---|---|---|---|
| DEM | Shuttle Radar Topography Mission (SRTM) | 30 m | 2000 | LS factor, slope, P factor |
| Precipitation | CHIRPS | 5 km | 2001–2025 | R factor |
| Soil properties | OpenLandMap | 250 m | Static | K factor |
| Land cover (2001) | MODIS MCD12Q1 | 500 m | 2001 | C factor (2001) |
| Land cover (2025) | Dynamic World | 10 m | 2024–2025 | C factor (2025) |
| Watershed boundary | HydroBASINS Level 07 | Subbasin scale | Static | Study area delineation |
Soil property data, including sand, clay, and organic carbon content, were extracted from the OpenLandMap soil database at 250 m resolution (Hengl et al. 2017). Land cover information for 2001 was derived from the MODIS Land Cover Type product (MCD12Q1 Collection 6.1) at 500 m resolution (Sulla–Menashe, Friedl 2018), while land cover for 2025 was obtained from the Google Dynamic World dataset, which provides near-real-time land cover classification at 10 m resolution derived from Sentinel 2 imagery (Brown et al. 2022). The watershed boundary was delineated using the HydroBASINS Level 07 dataset (Lehner, Grill 2013). All datasets were processed and analysed using the Google Earth Engine (GEE) cloud computing platform for initial computation and exported for detailed analysis using Python with geospatial libraries.
RUSLE factors computation
Rainfall erosivity factor (R)
The rainfall erosivity factor represents the erosive power of rainfall and its associated runoff. The R factor was calculated using the equation proposed by Morgan et al. (1984), which has been widely applied in tropical regions:
where:– R is the rainfall erosivity factor (MJ·mm·ha−1·h−1·yr−1),
– P is the mean annual precipitation (mm).
Mean annual precipitation was derived from CHIRPS daily rainfall data by computing annual totals for each year from 2001 to 2025 and subsequently averaging across the 25-year record. This approach ensures that the R factor reflects the full range of rainfall conditions experienced during the study period, rather than a subset that might not be representative of the complete temporal scope of the erosion analysis.
Soil erodibility factor (K)
The soil erodibility factor represents the susceptibility of soil to erosion based on its intrinsic physical and chemical properties. The K factor was calculated using the equation developed by Williams et al. (1984):
where:– Fcsand is the coarse sand factor,
– Fsilcl is the silt and clay factor,
– Forgc is the organic carbon factor,
– Fhisand is the high sand content factor.
These subfactors were calculated as follows:
where:– SAN, SIL, and CLA are the sand, silt, and clay contents (%), respectively,
– OC is the organic carbon content (%), and SN1 = (1 - SAN/100).
These soil properties were obtained from the OpenLandMap soil database (Hengl et al. 2017).
Topographic factor (LS)
The topographic factor accounts for the combined effect of slope length and steepness on erosion. In this study, the LS factor was estimated using a simplified slope-based approach derived entirely from the SRTM DEM at 30 m resolution. A quadratic relationship between slope gradient and the LS factor was applied:
where S is the slope gradient in degrees.This simplified formulation captures the non-linear increase in erosion potential with slope steepness, consistent with the principles described by Moore and Burch (1986) and the slope steepness relationships of McCool et al. (1987). The LS values were capped at a maximum of 50 to prevent extreme values in areas of very steep terrain. This slope-based approach was adopted as a computationally efficient alternative to flow accumulation-based methods, avoiding potential inconsistencies that can arise when deriving flow accumulation from different hydrological datasets at the watershed scale.
Cover management factor (C)
The cover management factor represents the effect of land cover and cropping management practices on soil erosion. C factor values were assigned based on land cover classifications for both 2001 and 2025. Since different land cover products were used for the two time periods, separate classification schemes were applied (Table 2).
Table 2.
C factor values assigned to land cover classes for 2001 (MODIS) and 2025 (Dynamic World).
| Panel A: MODIS MCD12Q1 classes (2001) | |||
|---|---|---|---|
| MODIS class | Land cover type | C value | Source |
| 1–5 | Forest (all types) | 0.03 | Hurni (1985) |
| 6 | Closed shrubland | 0.03 | Bewket and Teferi (2009) |
| 7 | Open shrubland | 0.05 | Bewket and Teferi (2009) |
| 8 | Woody savannas | 0.05 | Bewket and Teferi (2009) |
| 9 | Savannas | 0.08 | Renard et al. (1997) |
| 10 | Grasslands | 0.01 | Renard et al. (1997) |
| 11 | Permanent wetlands | 0.01 | Renard et al. (1997) |
| 12 | Croplands | 0.21 | Hurni (1985) |
| 13 | Urban and built up | 0.00 | Gelagay and Minale (2016) |
| 14 | Cropland/natural mosaic | 0.15 | Hurni (1985) |
| 16 | Barren or sparsely vegetated | 0.45 | Renard et al. (1997) |
| 17 | Water bodies | 0.00 | N/A |
| Panel B: Dynamic world classes (2025) | |||
|---|---|---|---|
| DW class | Land cover type | C value | Source |
| 0 | Water | 0.00 | N/A |
| 1 | Trees | 0.03 | Hurni (1985) |
| 2 | Grass | 0.01 | Renard et al. (1997) |
| 3 | Flooded vegetation | 0.01 | Renard et al. (1997) |
| 4 | Crops | 0.21 | Hurni (1985) |
| 5 | Shrub and scrub | 0.03 | Bewket and Teferi (2009) |
| 6 | Built up | 0.00 | Gelagay and Minale (2016) |
| 7 | Bare | 0.45 | Renard et al. (1997) |
| 8 | Snow and ice | 0.00 | N/A |
For the 2001 assessment, MODIS MCD12Q1 Land Cover Type 1 (IGBP classification) classes were reclassifed to corresponding C values based on published literature (Hurni 1985, Renard et al. 1997, Bewket, Teferi 2009). For the 2025 assessment, Dynamic World land cover classes (Brown et al. 2022) were assigned C values from the same literature sources, ensuring consistent erosion sensitivity assignments for equivalent land cover types across both periods. The C values assigned to each land cover type are properties of the vegetation and management condition, not of the assessment period, and therefore remain constant for equivalent cover types regardless of year.
Support practice factor (P)
The conservation support practice factor reflects the effect of soil conservation measures on erosion rates. Conservation practices, such as contouring, strip cropping, terracing, and construction of stone bunds, reduce erosion by modifying flow patterns and reducing the effective slope length. Since detailed spatial information on the distribution and condition of conservation structures was not available at the watershed scale, P-values were assigned based on slope classes following the approach of Wischmeier and Smith (1978) and adapted for Ethiopian highland conditions by Hurni (1985), as presented in Table 3. This slope-based assignment serves as a proxy under the assumption that steeper slopes are less amenable to effective conservation practice implementation.
Table 3.
P factor values assigned to slope classes.
| Slope class [%] | P-value | Description |
|---|---|---|
| 0–7 | 0.55 | Gentle slopes with moderate conservation potential |
| 7–11.3 | 0.60 | Moderate slopes with some conservation measures |
| 11.3–17.6 | 0.80 | Moderately steep slopes with limited conservation |
| 17.6–26.8 | 0.95 | Steep slopes with minimal conservation |
| >26.8 | 1.00 | Very steep slopes with no effective conservation measures |
Data processing and analysis
All geospatial data were processed and analysed using a combination of GEE for initial computation and factor derivation, and Python with geospatial libraries (including NumPy, Rasterio, GeoPandas, and SciPy) for detailed statistical analysis and visualisation. Factor layers generated in GEE were exported as GeoTIFF files and resampled to a common spatial resolution of 100 m for pixel level analysis.
Soil erosion rates were calculated for 2001 and 2025 by applying Eq. (1) with the respective C factor maps, while R, K, LS, and P factors were held constant across both time periods. It is important to note that because the C factor is the only temporally variable input in the model, any observed differences in erosion between the two assessment periods are, by construction, a consequence of changes in land cover. This methodological constraint is acknowledged in the interpretation of results.
Based on the calculated erosion rates, the watershed was classified into six erosion risk categories following Thapa (2020):
– Low [0–5 Mg·ha−1·yr−1],
– Moderate [5–10 Mg·ha−1·yr−1],
– High [10–20 Mg·ha−1·yr−1],
– Very High [20–40 Mg·ha−1·yr−1],
– Severe [40–80 Mg·ha−1·yr−1], and
– Very Severe [>80 Mg·ha−1·yr−1].
To analyse the relationship between topography and erosion dynamics, the watershed was divided into four slope classes: 0–5°(gentle), 5–15°(moderate), 15–25°(steep), and >25°(very steep). Mean erosion rates were calculated for each slope class for both 2001 and 2025.
Statistical analyses included calculation of descriptive statistics (mean, median, standard deviation, minimum, and maximum) for erosion rates, computation of areas and proportions under different erosion risk categories, and quantification of absolute and percentage changes. Paired t tests were applied to compare pixel-level erosion rates between 2001 and 2025 for each slope class, and effect sizes were quantified using Cohen’s d. Spatial correlation analysis between C factor change and the erosion change was conducted using Pearson’s correlation coefficient. Sensitivity analysis was performed to identify which input factors most strongly influenced the spatial variation of model outputs. Standardised regression coeffcients were computed by regressing standardised RUSLE factor values (R, K, LS, C, and P) against standardised erosion values at 50,000 randomly sampled pixel locations within the watershed. To quantify the uncertainty associated with literature-derived factor values, a Monte Carlo simulation (n = 1000) was conducted by independently varying C and P factor values for each assessment year within plus or minus 20% of assigned values (Panagos et al. 2014). The simulated distribution of erosion reduction percentages was used to derive 95% confidence intervals around the observed estimate.
Climate variability assessment
To ensure that observed erosion changes were not confounded by rainfall variability, CHIRPS daily precipitation data were analysed for the full study period (2001–2025). Annual precipitation totals were computed for each year, and the study period was divided into an early phase (2001–2013) and a late phase (2014–2025). A two sample Welch’s t test was applied to compare the mean annual precipitation between the two phases. The coefficient of variation (CV) was calculated to assess interannual rainfall stability. This analysis allows assessment of whether any systematic shift in rainfall could have contributed to observed changes in erosion rates independently of land management changes.
Results
Spatial distribution of RUSLE factors
Rainfall erosivity factor (R)
The R factor across the Muger watershed ranged from 401.5 to 538.7 MJ·mm·ha−1·h−1·yr−1, with a mean of 459.1 MJ·mm·ha−1·h−1·yr−1 and a standard deviation of 27.1 (Fig. 2A). Higher R factor values were observed in the northern and western portions of the watershed, corresponding to areas of higher mean annual precipitation. The spatial pattern exhibited a gradual decrease from northwest to southeast, following the general precipitation gradient of the central Ethiopian highlands.
Fig. 2
Spatial distribution of RUSLE model factors across the Muger River Watershed: (a) – rainfall erosivity (R), (b) – soil erodibility (K), (c) – topographic factor (LS), (d) – cover management factor (C) for 2001, (e) – cover management factor (C) for 2025, and (f) – conservation support practice factor (P).
Soil erodibility factor (K)
The K factor ranged from 0.028 to 0.043 Mg·h·MJ−1·mm−1, with a mean of 0.036 Mg·h·MJ−1·mm−1 and a standard deviation of 0.002 (Fig. 2B). The spatial distribution of K factor values was relatively homogeneous compared to other RUSLE factors, with slightly higher erodibility observed in the central portion of the watershed. The predominant soil types in the watershed (Nitisols and Vertisols) exhibited moderate erodibility values, consistent with clay-loam to clay-textured soils.
Topographic factor (LS)
The LS factor exhibited the highest spatial variability among all RUSLE factors, ranging from 0.065 to 44.4 with a mean value of 1.89 and a standard deviation of 2.74 (Fig. 2C). Higher LS values were strongly associated with steeper slopes and topographically diverse areas. The northeastern and southwestern portions of the watershed, characterised by rugged terrain and steep slopes, displayed the highest LS values (>10), while the central plains and valley bottoms showed substantially lower values (<1).
Cover management factor (C)
The C factor exhibited notable differences between 2001 and 2025 (Figs 2D and 2E). In 2001, C factor values ranged from 0.01 to 0.21 with a mean of 0.183, reflecting the predominantly agricultural character of the watershed at that time. By 2025, the range expanded from 0.00 to 0.45 while the mean decreased to 0.119, representing an approximately 35% improvement in the average level of vegetative protection against erosion. The spatial pattern of C factor changes revealed substantial shifts in land cover composition, with the most notable improvements occurring in previously cultivated areas across the central and western portions of the watershed.
Support practice factor (P)
The P factor ranged from 0.55 to 1.00 across the watershed, with a mean value of 0.79 (Fig. 2F). Since the P factor was assigned on the basis of slope categories, its spatial pattern closely followed the topography of the watershed. Lower P values (0.55–0.60), indicating greater conservation effectiveness, were associated with gentler slopes in the central plains, while higher values (0.95–1.00) characterised the steeper terrain in the northeastern and southwestern regions.
Soil properties
The spatial distribution of soil properties within the watershed is presented in Figure 6 and summarised in Table 4. Sand content ranged from 16% to 46% (mean 32.7%, SD 4.6%), clay content from 26% to 62% (mean 40.3%, SD 4.0%), silt content from 19% to 39% (mean 27.0%, SD 2.9%), and soil organic carbon from 1 to 17 g·kg−1 (mean 4.6 g·kg−1, SD 1.1 g·kg−1). The high mean clay content is consistent with the dominance of Nitisols and Vertisols in the watershed, while the relatively low organic carbon content reflects the advanced state of soil degradation in the intensively cultivated central Ethiopian highlands. The moderate spatial variability in textural properties, as indicated by the standard deviations, contributes to the limited range of K factor values observed across the watershed.
Soil erosion rates and spatial distribution
Erosion rates in 2001
The estimated soil erosion rates for 2001 ranged from 0.005 to 117.6 Mg·ha−1·yr−1, with a mean of 4.80 Mg·ha−1·yr−1 and a median of 2.00 Mg·ha−1·yr−1 (Table 5). The standard deviation of 7.65 Mg·ha−1·yr−1 indicated considerable spatial variability in erosion rates across the watershed. Approximately 70.6% of the watershed experienced low erosion rates (0–5 Mg·ha−1·yr−1), while 15.7% showed moderate erosion (5–10 Mg·ha−1·yr−1). Areas with erosion rates exceeding 10 Mg·ha−1·yr−1 constituted 13.7% of the watershed, with the highest rates concentrated on steep slopes in the northeastern and southwestern regions (Fig. 3A).
Fig. 3.
Spatial distribution of soil erosion risk in the Muger River Watershed: (a) – 2001, (b) – 2025, and (c) – net change in erosion rates (2025 minus 2001). Negative values in (c) indicate erosion reduction.
Erosion rates in 2025
By 2025, soil erosion rates ranged from 0 to 156.6 Mg·ha−1·yr−1, with a mean of 2.18 Mg·ha−1·yr−1 and a median of 0.83 Mg·ha−1·yr−1 (Table 5). The standard deviation decreased to 3.90 Mg·ha−1·yr−1, indicating reduced spatial variability relative to 2001. Low erosion areas expanded to cover 88.0% of the watershed, while moderate erosion areas decreased to 8.1%. The proportion of the watershed experiencing erosion rates exceeding 10 Mg·ha−1·yr−1 decreased to 3.8% (Fig. 3B).
Changes in erosion patterns (from 2001 to 2025)
The comparison of erosion rates between 2001 and 2025 revealed a significant overall decrease in soil loss. Mean erosion decreased by 2.63 Mg·ha−1·yr−1 (54.7% reduction), while the median decreased by 1.17 Mg·ha−1·yr−1 (58.7% reduction). The spatial change map (Fig. 3C) shows that erosion reduction was not uniform across the watershed but exhibited spatial patterns closely related to topography and the distribution of land cover change.
The areal extent of different erosion risk classes changed substantially between 2001 and 2025 (Table 6, Fig. 4). The area under low erosion risk increased by 1282.1 km2 (24.7%), while all higher risk categories showed decreases. The most substantial proportional reductions occurred in the Severe category (88.6% decrease) and the Very High category (81.1% decrease).
Fig. 4.
Comparison of watershed area under each erosion risk class for 2001 and 2025: (a) – area in km2 and (b) – percentage of watershed area.
Table 6.
Area distribution by erosion risk class in 2001 and 2025.
| Erosion Risk Class | 2001 area | 2001 | 2025 area | 2025 | Change in area | Change |
|---|---|---|---|---|---|---|
| [km2] | [%] | [km2] | [%] | [km2] | [%] | |
| Low (0 to 5) | 5,185.5 | 70.6 | 6,467.6 | 88.0 | +1,282.1 | +24.7 |
| Moderate (5 to 10) | 1,151.8 | 15.7 | 595.8 | 8.1 | –556.0 | –48.3 |
| High (10 to 20) | 660.1 | 9.0 | 221.6 | 3.0 | –438.5 | –66.4 |
| Very high (20 to 40) | 287.7 | 3.9 | 54.4 | 0.7 | –233.3 | –81.1 |
| Severe (40 to 80) | 60.0 | 0.8 | 6.8 | 0.1 | –53.2 | –88.6 |
| Very severe (>80) | 1.4 | <0.1 | 0.3 | <0.1 | –1.1 | –76.4 |
Soil erosion by slope class
Analysis of erosion rates by slope class revealed a strong and consistent relationship between topography and erosion dynamics (Table 7, Fig. 5). In 2001, mean erosion rates increased progressively with slope gradient, from 0.48 Mg·ha−1·yr−1 on gentle slopes (0–5°) to 23.49 Mg·ha−1·yr−1 on very steep slopes (>25°). By 2025, this general pattern persisted but the magnitude of erosion was substantially reduced across all slope classes, with the largest absolute and proportional reductions observed on steeper terrain.
Fig. 5.
Mean annual soil erosion rates by slope class for 2001 and 2025, with error bars representing standard deviation and percentage change annotated above each pair.
Fig. 6.
Spatial distribution of soil properties across the Muger River Watershed: (a) – sand content (%), (b) – clay content (%), (c) – silt content (%), and (d) – organic carbon content (g·kg−1). Data derived from OpenLandMap at 250 m resolution.
Table 7.
Mean erosion rates by slope class in 2001 and 2025.
| Slope Class | Area | Erosion 2001 | Std dev 2001 | Erosion 2025 | Std dev 2025 | Absolute change | Percent change |
|---|---|---|---|---|---|---|---|
| [°] | [cell grids] | [Mg·ha−1·yr−1] | [Mg·ha−1·yr−1] | [%] | |||
| 0–5 | 243.745 | 0.48 | 0.69 | 0.39 | 0.57 | –0.10 | –20.2 |
| 5–15 | 335.893 | 3.26 | 2.54 | 1.97 | 2.30 | –1.29 | –39.6 |
| 15–25 | 114.306 | 10.13 | 5.26 | 4.72 | 5.65 | –5.41 | –53.4 |
| >25 | 51.575 | 23.49 | 14.95 | 6.39 | 8.09 | –17.10 | –72.8 |
The most pronounced improvement occurred in the steepest slope class (>25°), where mean erosion rates decreased by 17.10 Mg·ha−1·yr−1, representing a 72.8% reduction. The 15–25-degree class showed a 53.4% decrease, while the 5–15-degree class decreased by 39.6%. Gentle slopes (0–5°) exhibited a more modest 20.2% reduction, reflecting the already low baseline erosion rates in these areas.
Paired t tests confirmed that erosion reductions were statistically significant across all slope classes (p < 0.001 in all cases; Table 8). Effect sizes, quantified using Cohen’s d, increased systematically with slope gradient: negligible for gentle slopes (d = 0.18), medium for moderate slopes (d = 0.54), large for steep slopes (d = 0.85), and large for very steep slopes (d = 1.12). This pattern indicates that the magnitude of erosion improvement increased substantially with slope steepness.
Land cover change and C factor dynamics
Land cover analysis revealed substantial compositional changes between 2001 and 2025. In 2001, the watershed was overwhelmingly dominated by cropland (86.3% of the watershed area within the valid analysis domain), with grassland covering 13.0% and forest and shrubland together accounting for <1%. By 2025, cropland had decreased to 48.6%, while tree cover expanded to 15.7%, and shrub and scrub cover increased to 23.2%. Built-up areas emerged at 4.5%, and bare land occupied 1.1% of the watershed.
These land cover changes directly influenced the C factor distribution. In 2001, approximately 85% of all watershed pixels had a C value of 0.21, corresponding to the cropland class, while approximately 12% had a value of 0.01, corresponding to grassland. By 2025, the distribution was more heterogeneous: 48.6% of pixels retained the cropland value of 0.21, while 38.9% shifted to a C value of 0.03 (corresponding to trees and shrub cover), 6.7% had a value of 0.01 (grass), and 4.7% had a value of 0.00 (water and built-up areas). This redistribution resulted in the observed decrease in mean C factor from 0.183 to 0.119.
It is noted that the land cover comparison involves two different classification products at different spatial resolutions (MODIS at 500 m for 2001 and Dynamic World at 10 m for 2025), which may introduce some inconsistencies in the absolute area estimates for individual land cover classes. The differences in spatial resolution mean that the coarser MODIS product is more likely to classify mixed pixels under the dominant class (cropland), potentially overestimating cropland extent in 2001. Accordingly, the emphasis of this analysis is placed on the net change in the C factor and its relationship to erosion, rather than on direct class-to-class transition mapping.
Spatial correlation analysis between C factor change and erosion change yielded a Pearson’s correlation coefficient of 0.42 (p < 0.001), indicating that approximately 17% of the spatial variability in erosion reduction could be explained by changes in the C factor alone. The remaining variability reflects the multiplicative interaction between C factor change and the spatially varying topographic (LS) and other factors in the RUSLE equation: a given change in C value produces larger erosion reductions where LS values are high (steep slopes) than where they are low (gentle terrain).
Sensitivity and uncertainty analysis
Sensitivity analysis based on standardised regression coefficients at 50,000 randomly sampled pixels revealed that the spatial variation of erosion across the watershed is primarily governed by the topographic factor. The LS factor yielded the highest standardised coefficient (0.800), followed by the C factor (0.133), P factor (0.101), K factor (0.089), and R factor (0.006). The overall regression model explained 75.8% of the variance in spatial erosion distribution (R2 = 0.758). The dominance of the LS factor reflects the strong topographic control on erosion patterns within the watershed, consistent with the large range of LS values (0.065–44.4) relative to the more limited spatial variability of other factors.
It is important to distinguish this spatial sensitivity from the temporal attribution of erosion change. Because the R, K, LS, and P factors were held constant between the two assessment periods, any observed difference in erosion between 2001 and 2025 is, by model construction, entirely attributable to changes in the C factor. The spatial sensitivity analysis, therefore, characterises which factors control where erosion is highest or lowest, while the temporal change in erosion reflects where and how much land cover has changed, amplified by the local values of the other factors.
Monte Carlo simulation (n = 1000) with independent random variation of C and P factors within plus or minus 20% of assigned values for each assessment year produced a mean simulated erosion reduction of 54.4%, with a 95% confidence interval of 38.7–67.3%. The observed reduction of 54.7% falls near the centre of this distribution. The standard deviation of simulated reductions was 7.6%. These results indicate that even under the most conservative uncertainty assumptions, the watershed experienced substantial erosion reduction over the study period.
Climate variability analysis
Analysis of CHIRPS precipitation data for the full study period (2001–2025) revealed no significant temporal trend in mean annual precipitation. Mean annual precipitation was 1183.7 mm for the early period (2001–2013, n = 13 years) and 1221.2 mm for the late period (2014–2025, n = 12 years), representing a non-significant 3.2% increase (Welch’s t = – 1.151, df = 23.0, p = 0.262). The CV of 6.8% indicated low interannual variability, characteristic of relatively stable subhumid tropical climates. These findings confirm that the observed 54.7% reduction in mean soil erosion cannot be attributed to systematic changes in rainfall patterns and is instead consistent with the documented changes in land cover and conservation practices.
Discussion
Drivers of erosion reduction and the role of land cover change
The 54.7% reduction in mean soil erosion observed in the Muger River Watershed between 2001 and 2025 represents a substantial improvement in the erosion status of this highland landscape. This reduction is reflected in the decrease in mean C factor from 0.183 to 0.119, corresponding to a 35% improvement in the average vegetative protection provided by the land surface. The expansion of tree cover from <1% to 15.7% and shrub cover from <1% to 23.2%, alongside a decline in cropland from 86.3% to 48.6%, are consistent with the well documented trajectory of land rehabilitation in the Ethiopian highlands driven by conservation programmes and changing land management practices (Haregeweyn et al. 2015, Nyssen et al. 2015, Gashaw et al. 2019).
An important methodological consideration, raised during peer review of an earlier version of this manuscript, concerns the interpretation of factor importance in a temporal comparison where only one factor varies. Because R, K, LS, and P are held constant between the 2001 and 2025 assessments, the observed differences in erosion are, by the structure of the RUSLE equation, entirely a function of changes in C. This does not diminish the validity of the finding; rather, it means that the erosion reduction quantified here reflects the cumulative effect of land cover improvements as mediated by the local topographic and soil conditions. The observation that the greatest absolute reductions occurred on the steepest slopes (17.10 Mg·ha−1·yr−1 on slopes >25°) illustrates this interaction: the same proportional decrease in C produces a far larger erosion reduction where the LS factor amplifies the erosive potential. The spatial sensitivity analysis provides complementary insight. Across the watershed at any single point in time, erosion rates are over whelmingly controlled by topography (LS coefficient = 0.800), with C contributing modestly to spatial variability (coefficient = 0.133). This distinction between factors that control the spatial pattern of erosion and those that drive its temporal change is essential for accurate interpretation and should be recognised in erosion modelling studies that compare different time periods using empirical models such as RUSLE.
Role of specific conservation measures
The erosion reduction documented in this study reflects the cumulative impact of multiple conservation interventions implemented across the watershed over more than two decades. These interventions can be broadly grouped into physical conservation structures that modify overland flow and reduce effective slope length, and biological measures that increase vegetative cover and reduce the impact of rainfall on the soil surface. Physical structures, particularly stone bunds and soil bunds, have been the most widely implemented conservation measures in the Ethiopian highlands. Stone bunds reduce erosion by intercepting surface runoff, reducing flow velocity, and trapping eroded sediment, with documented effectiveness of 50–68% erosion reduction relative to untreated slopes in the Tigray highlands (Gebremichael et al. 2005, Nyssen et al. 2009). Soil bunds serve a similar function and have been shown to reduce runoff by 20–30% and soil loss by 40–60% in the central Ethiopian highlands (Adimassu et al. 2014, Belayneh et al. 2019). Bench terraces, though more labour-intensive to construct and maintain, are particularly effective on steep slopes where they substantially reduce the effective slope length and surface runoff (Amsalu, de Graaff 2007). Check dams installed in gullies and ephemeral drainage channels further contribute to sediment retention and reduction of channel incision (Grum et al. 2017). In the context of the RUSLE framework, these physical structures primarily influence erosion through their effect on the conservation support practice factor (P), though their influence cannot be fully captured by the slope-based P factor assignment used in this study.
Biological interventions directly reduce the cover management factor (C) by increasing vegetation density and ground cover. Area exclosures, in which degraded lands are protected from grazing and cultivation to allow natural vegetation regeneration, have been widely established across the Ethiopian highlands since the early 2000s (Descheemaeker et al. 2006). These exclosures promote recovery of native woody species and grasses, transforming degraded cropland or grazing land (C = 0.21 or higher) into shrubland or woodland (C = 0.03), representing a reduction in erosion susceptibility of approximately 85% for that factor alone. Active reforestation with native and exotic species complements this passive regeneration and has contributed to the expansion of tree cover observed in the 2025 land cover data (Mhiret et al. 2019). The combined effect of these biological measures is reflected in the substantial decrease in mean C factor documented in this study.
The synergistic combination of physical and biological conservation measures, whereby stone bunds reduce surface runoff and create moisture conditions favourable for vegetation establishment, is widely recognised as more effective than either approach alone (Nyssen et al. 2009, Tewodros et al. 2020). The comprehensive nature of Ethiopia’s SLMP and the PSNP, both of which promote integrated packages of physical and biological measures at the watershed scale, is consistent with the magnitude and spatial extent of erosion reduction observed in the present study.
Comparison with other studies
The baseline erosion rates observed in this study for 2001 (mean 4.80 Mg·ha−1·yr−1) are lower than the 42 Mg·ha−1·yr−1 reported by Hurni (1983) for cultivated fields in the Ethiopian highlands. This difference can be attributed to several factors: (1) Hurni’s estimates were based on plot scale measurements specifically targeting cultivated fields, whereas the present study provides watershed averaged values that include low erosion areas such as grasslands and valley bottoms; (2) the 100 m analysis resolution smooths localised erosion peaks compared to field measurements; (3) approximately 33% of the Muger watershed lies on gentle slopes (<5°) that contribute minimal erosion; and (4) some conservation activities predating the 2001 baseline may already have reduced erosion rates. On steep slopes (>25°), the 2001 mean of 23.49 Mg·ha−1·yr−1 more closely approaches Hurni’s estimates, suggesting that spatial averaging and methodological scale differences largely account for the apparent discrepancy.
The 54.7% erosion reduction observed across the watershed over 24 years exceeds most documented intervention impacts in the Ethiopian literature. Descheemaeker et al. (2006) reported approximately 40% erosion reduction from exclosures in the Tigray highlands after 7 years. Nyssen et al. (2009) documented 68% sediment reduction from stone bunds in small northern Ethiopian catchments. At larger scales, Haregeweyn et al. (2015) reported 38% reduction in the Enabered watershed, and Tamene et al. (2017) observed 41% reduction in the Debre Mawi watershed over shorter assessment periods. The larger reduction observed in the present study likely reflects the longer assessment period (24 years versus 5–10 years in most cited studies), allowing cumulative effects of conservation measures to develop, and the synergistic effect of multiple intervention types implemented across the watershed.
In a broader context, comparable multitemporal RUSLE studies in China (Wei et al. 2007), India (Mandal, Sharda 2013), and Thailand (Ongsomwang, Thinley 2015) have reported erosion reductions of 15–35% over similar time spans, though direct comparisons are limited by differences in environmental conditions, conservation intensity, and assessment methodologies.
Implications for sustainable land management
The findings of this study carry several implications for watershed management in the Ethiopian highlands. First, the results provide empirical evidence that sustained conservation investment over a period of more than two decades can achieve meaningful reversal of land degradation trends at the watershed scale. The Monte Carlo analysis, yielding a 95% confidence interval of 38.7–67.3% for the erosion reduction, confirms the robustness of this conclusion even under substantial parameter uncertainty.
Second, the disproportionate improvement on steep slopes highlights the importance of targeting conservation interventions according to landscape vulnerability. Slopes exceeding 25°, which constitute approximately 7% of the watershed, exhibited a 72.8% reduction in erosion and the largest effect size (Cohen’s d = 1.12). This suggests that investments in conservation on steep terrain yield particularly large returns. Future efforts should continue to prioritise remaining areas of high erosion risk, which in 2025 still account for 3.8% of the watershed, predominantly on steep slopes in the northeastern sector.
Third, the integrated application of physical structures and biological measures appears more effective than isolated interventions, as evidenced by the concurrent reduction in both effective slope length (through terraces and bunds) and surface vulnerability (through vegetation recovery). Policies that promote reforestation, establishment of area exclosures, and maintenance of existing conservation structures should therefore remain central to land management strategy.
Fourth, the decrease in areas classified as high to very severe erosion risk from 13.7–3.8% of the watershed represents a significant achievement in reducing sediment delivery to downstream water bodies, with implications for reservoir longevity and water quality in the Blue Nile basin (Wolancho 2015, Lemma et al. 2019, Assefa et al. 2021). Finally, the study demonstrates the practical value of multitemporal erosion assessment for evaluating conservation programme effectiveness and guiding adaptive management. Regular reassessment of erosion risk using updated land cover data can identify areas where conservation gains are being maintained, as well as areas where degradation may be resuming, enabling responsive allocation of limited resources.
Limitations and uncertainty
Several limitations should be acknowledged when interpreting the results. First, the RUSLE model estimates sheet and rill erosion only and does not capture gully erosion, mass movements, or bank erosion, all of which can be locally significant in the Ethiopian highlands (Nyssen et al. 2004). Second, the use of global datasets for soil properties, while necessary given the absence of detailed local surveys, may not fully capture fine-scale soil variability. Third, the assignment of P factor values based solely on slope is a simplification that does not reflect the actual spatial distribution and maintenance status of conservation structures (Hurni 1985). Field verification would improve the accuracy of this factor. Fourth, the comparison of land cover between 2001 and 2025 utilises products derived from different sensor systems at different spatial resolutions (MODIS at 500 m versus Dynamic World at 10 m). The coarser MODIS product is likely to overestimate the extent of the dominant class (cropland) in landscapes with fine-scale heterogeneity, while the finer Dynamic World product can resolve smaller patches of forest and shrubland. This resolution disparity may contribute to the apparent magnitude of land cover change and, consequently, to the estimated erosion reduction. The analysis therefore emphasises the integrated C factor change and its relationship to erosion dynamics, rather than relying on detailed class-to-class transition mapping. Fifth, the simplified slope-based LS factor formulation used in this study, while computationally straightforward and reproducible, does not incorporate flow accumulation and may therefore underestimate the LS factor in areas of concentrated flow convergence. However, this limitation applies equally to both assessment years and does not affect the temporal comparison.
Despite these limitations, the convergence of evidence from multiple analytical approaches (erosion rate comparison, risk class analysis, slope stratified analysis, climate assessment, and uncertainty analysis) supports the principal finding of substantial, statistically significant erosion reduction over the 24-year study period.
Conclusion
Soil erosion dynamics in the Muger River Watershed were quantified over 24 years (2001–2025) using the RUSLE model integrated with multisource geospatial data. Mean annual soil erosion decreased from 4.80 to 2.18 Mg·ha−1·yr−1, representing a 54.7% reduction. The proportion of the watershed under low erosion risk (0–5 Mg·ha−1·yr−1) expanded from 70.6% to 88.0%, while the proportion exceeding 10 Mg·ha−1·yr−1 decreased from 13.7% to 3.8%. Erosion reductions were statistically significant across all slope classes (p < 0.001), with the largest effect sizes on slopes exceeding 25° (Cohen’s d = 1.12, 72.8% reduction). The mean cover management factor decreased from 0.183 to 0.119, reflecting a shift from cropland-dominated to more heterogeneous land cover. No significant trend in mean annual precipitation was detected over the study period (p = 0.262), confirming that the erosion reduction is not attributable to changes in rainfall. Monte Carlo simulation indicated that the finding of substantial erosion reduction is robust to parameter uncertainty (95% confidence interval: 38.7-67.3%). Sensitivity analysis demonstrated that topography governs the spatial distribution of erosion, while temporal changes in erosion rates are a function of land cover improvement.
Data availability statement
The datasets generated and analysed during this study are available from the corresponding author on reasonable request. Source datasets include SRTM DEM (30 m), CHIRPS precipitation data (2001 to 2025), OpenLandMap soil properties, MODIS MCD12Q1 land cover (2001), and Google Dynamic World land cover (2024 to 2025). The Muger watershed boundary was delineated from HydroBASINS Level 07. All source datasets are publicly available from their respective providers as cited in the manuscript. Processing was conducted using Google Earth Engine and Python.
Acknowledgements
The author thanks the anonymous reviewers whose constructive comments significantly improved this manuscript. The author acknowledges the data providers: NASA/USGS for the SRTM Digital Elevation Model, the Climate Hazards Group (University of California, Santa Barbara) for CHIRPS precipitation data, OpenLandMap/ISRIC for soil property datasets, NASA for the MODIS land cover products, Google for the Dynamic World land cover dataset, and WWF/McGill University for the HydroBASINS dataset.