Assessment of the combined effects of land use/land cover and climate change on soil erosion in the Sile watershed, Ethiopian Rift Valley Lakes Basin

Abstract

The combined impacts of land use/land cover (LULC) changes and climate change are significantly affecting soil and water ecologies, calling for proper natural resource management. The aims of this study were to assess the LULC dynamics between 1987 and 2050 and estimate the impact of LULC and climate changes (under two climate scenarios (RCP4.5 & RCP8.5)) on soil erosion rate in Sile watershed, Ethiopian rift valley lakes basin (ERVLB). A supervised method with maximum likelihood classification was applied in the ERDAS Imagine software for assessing the LULC of 1987, 2003, and 2020. A classified past LULC map serves as a baseline to predict 2050 LULC by considering different driver variables using Multi-Layer Perceptron Neural Network and Cellular Automata-Markov Chain Model integrated with TerrSet software. A revised universal soil loss equation model (RUSLE) has been used to compute an input parameter in ArcGIS to estimate the rate of soil loss. The input parameters were erosivity (R-factor), soil erodibility (K-factor), topography (LS-factor), conservation practice (P-factor), and cover (C-factor). The findings revealed that with the business-as-usual scenario, the substantial LULC alteration will be prolonged over the coming several years. The average soil loss under the RCP4.5 and RCP8.5 scenarios with a future 2050 LULC was predicted to be 56.48 t-ha⁻¹ year⁻¹ and 57.11 t-ha⁻¹ year⁻¹ and then it will be expected to increase 36.40% and 38.19%, respectively. In order to sustainably manage the natural resources, interventions in prioritized sub-watersheds have been suggested at present and for the future. Additionally, it needs ecological sound conservation policies and strategies to mitigate the adverse impacts of climate change on soil erosion.

Keywords:

• LULC
• soil loss
• RUSLE
• climate change
• CA-Marcov

1.

PUBLIC INTEREST STATEMENT

Climate change is a global concern affecting many landscape characteristics. Land use changes as an important variable in soil erosion are expected to change under projecting climatic scenarios. Evaluation of future soil erosion is crucial for devising appropriate interventions. Moreover, research results from such studies contribute to planning, decision-making, and policy development which could make possible sustainability of natural resources.

1. Introduction

Land is an essential resource and principal basis for human well-being and livelihood as it provides us food, shelter, and multiple other ecosystem services (Mahata & Sharma, 2021; Prăvălie, 2021). As the global population continues to grow, so does the demand for land. However, land is a finite resource, and its availability is limited (Mueller et al., 2014). Therefore, the competition among these uses puts immense pressure on available land resources (UNCCD, 2019). Consequently, the lands’ positive multidimensional role has been heavily altered over the past decades, amid large-scale land degradation that has affected numerous ecosystems across the globe (Prăvălie, 2021).

Land degradation is a pressing environmental issue that poses a substantial threat to the ecosystems, food security, and socio-economic stability worldwide (AbdelRahman, 2023; Roy et al., 2022). These impacts tend to be more severe for communities whose livelihoods directly depend on land resources (Pricope et al., 2022). Land degradation is influenced by a combination of social (e.g., population pressure, land tenure system, and cultural practices), economic (e.g., market demand and economic policies), ecological (e.g., climate change, loss of biodiversity, and soil erosion), and physical (e.g., geomorphology, hydrology, and natural disasters) drivers. These drivers can vary in significance depending on the region and context, but they often interact and exacerbate each other (Mahata & Sharma, 2021).

Soil erosion is one of the most important forms of land degradation (Prăvălie, 2021; Vanmaercke et al., 2021). Soil erosion results in a multitude of negative effects, such as reducing the water-holding capacity and infiltration rates, decreasing organic matter and nutrients, degrading soil structure, and reducing soil depth (Poesen, 2018; Sartori et al., 2019). Soil erosion under climate change is most directly affected by the changes in extreme precipitation (Eekhout & de Vente, 2022; Roy et al., 2022; Teng et al., 2018). According to the Fifth Assessment Report (AR6) of the IPCC (Intergovernmental Panel on Climate Change), global mean precipitation and surface temperature have changed significantly, and these changes are very likely to continue during the twenty-first century (IPCC, 2022). Recent projections indicate that the rate of soil erosion will increase globally (Borrelli et al., 2020). It is also well thought that land use land cover (LULC) change has a significant impact on soil erosion as well, and interacts with climate change (Bogale, 2020; Eekhout & de Vente, 2022).

Assessing the impact of LULC and climate change on soil erosion is crucial to understand its environmental challenges, identify erosion-prone areas and develop effective soil and water conservation (SWC) measures (Kanito & Feyissa, 2021). Estimation of soil erosion involves a diverse array of methodologies and tools (Borrelli et al., 2021; Guo et al., 2019). Field-based methods like sediment sampling and erosion plot studies provide more detailed and direct measurement (Hsieh, 2022). Empirical and process-based models consider factors such as climate, soil properties, slope, and land cover to predict erosion rates (Borrelli et al., 2021). Among them, the Revised Universal Soil Loss Equation (RUSLE) has been extensively used to estimate the long-term effects of climate and LULC changes on soil erosion rates (Borrelli et al., 2021; Teng et al., 2018). Moreover, remote sensing and Geographic Information Systems (GIS) play a pivotal role in soil erosion assessment, offering powerful tools to monitor and combat this critical environmental challenge (Chaminé et al., 2021). The integration of RUSLE with GIS can provide a comprehensive erosion quantification framework (Girma & Gebre, 2020; Moges & Bhat, 2017; Phinzi & Ngetar, 2019).

Soil erosion in the Ethiopian Rift Valley lakes basin (ERVLB) is a significant environmental challenge, with widespread and severe consequences (Aga et al., 2018; Taye et al., 2023). It is exacerbated by anthropogenic (e.g., deforestation, unsustainable agricultural practices, and overgrazing) and natural factors (e.g., rainfall and topographic variability). The combination of these factors has led to extensive soil degradation, loss of fertile topsoil, reduced agricultural productivity, biodiversity loss, and sedimentation in rivers and lakes (Ayalew et al., 2022; Moges et al., 2020; Roy et al., 2022). The expansion of agricultural land at the expense of forestlands, grasslands, and shrublands was the most detrimental factor for severe soil erosion risks (Girma et al., 2023; Moisa et al., 2021; Wassie, 2020). This could get worse in the future because of the forecasted increase in population and the occurrence of extreme climate events throughout the 21^st century (Eekhout & de Vente, 2022; Saha et al., 2022).

Given the potential for climate and LULC change to increase soil erosion, assessing climate and LULC change impact on soil erosion is essential to implement proper watershed management interventions and local climate adaption strategies (Belay & Mengistu, 2021; Moisa et al., 2021). Moreover, interventions are local-characteristics specific that could not be obtained from recommendations on the shelf, thus requiring specific studies in areas of concern. In this context, no study has yet been conducted to provide scientific evidence about the impact and trends of past and future climate and LULC dynamics in the Sile river watershed, ERVLB. This knowledge gap is crucial as it allows for a more comprehensive understanding of the environmental changes occurring in the watershed. Therefore, the objectives of this study were to (i) assess the long-term (1987–2050) LULC changes and (ii) quantify the combined impact of LULC and climate changes on soil loss rate in the Sile river watershed.

2. Materials and methods

2.1. Study area

Sile River watershed is located in the ERVLB and geographically located between 5°50″30′–6°4″30′ N latitude and 37°17′30″–37°35′0″ E longitude (Figure 1). Sile River flows from west to east and drains into Lake Chamo. The watershed encompasses two districts, Arbaminch zuria and Bonke districts, and covers an area of 27,610 hectares. The watershed receives mean annual rainfall of 941 mm and average minimum and maximum temperatures were 14°C and 33.5°C, respectively.

Figure 1. Location of the study area.

2.2. Dataset

To estimate the soil loss and identify erosion hotspots in the watershed, various datasets from different sources were collected. The observed daily rainfall data were collected from the Ethiopian Meteorology Institute for the period of 1987 to 2020 (Table 1). Missing rainfall data were filled by using the Normal Ratio method (Sanusi et al., 2017). A total of 102 composite soil samples were collected using stratified random sampling. Soil textural classes were identified using the standard hydrometer method (Gavlak et al., 2005) and soil organic carbon content was determined by Walkley and Black method (Walkley & Black, 1934). Cloud-free satellite imageries were downloaded from the United States Geological Survey (https://earthexplorer.usgs.gov/) (Table 1). The acquisition years (1987, 2003, and 2020) were selected considering the nationwide famine, regime change of the country, and the recent green legacy initiative, respectively (Girma et al., 2022).

Table 1. The metrological stations and satellite data used in this study

Characteristics of the rain gauge stations
Station name	Latitude	Longitude	Elevation (meter)	Mean annual rainfall (mm)
Arbaminch	334487	666294	1,336	888
MirabAbaya	335240	672585	1,223	764
Gidole	309647	667316	2,248	1,257
Gato	327282	643820	1,771	857
Satellite images used in this study for LULC change analysis
Satellite image	Path/Row	Sensor*	Date of acquisition	Spatial resolution
Landsat 4–5 TM	169/056	TM	31 January 1987	30*30
Landsat 7ETM	169/056	ETM+	2 February 2003	30*30
Landsat-8 OLI/TIRS	169/056	OLI/TIRS	27 February 2020	30*30

*TM (Thematic Mapper), ETM+ (Enhanced Thematic Mapper Plus), OLI/TIRS (Operational Land Imager/Thermal Infrared Sensor).

Topographic features were characterized using 30 m resolution digital elevation model (DEM) available at https://search.earthdata.nasa.gov/search. In addition, population density and proximity to road and river were considered as driver variables to predict the future (2050) LULC change (Girma et al., 2022). These data were acquired from OpenStreetMap (https://www.openstreetmap.org/) and Water and Land Resources Center (wlrc-eth.org). To estimate the future soil loss, an ensemble of four RCM models (Table 2) under two representative concentration pathway (RCP) (4.5 and 8.5) scenarios were considered. These RCMs were downloaded from the Coordinated Regional Climate Downscaling Experiment (CORDEX) Africa domain (https://esgf-data.dkrz.de/search/esgf-dkrz/).

Table 2. Description of the driving GCMs used in this study

Regional climate model	Model institute	Country	Driving GCM
Canadian Regional Climate Model version 4 (CanRCM4)	Canadian Centre for Climate Modelling and Analysis (CCCMA)	Canada	CCCma_canESM2
Rossby Center Regional Atmospheric Model (RCA4)	The Swedish meteorological and hydrological institute (SMHI)	Sweden	MOHC-HadGEM2-ES
Rossby Center Regional Atmospheric Model (RCA4)	The Swedish meteorological and hydrological institute (SMHI)	Sweden	CCCma_canESM2
Regional Atmospheric Climate Model Version 2.2 T (RACMO22T)	Koninklijk Nederlands Meteorologisch Instituut (KMNI)	Netherland	ICHEC-EC-EARTH

2.3. Land use land cover analysis

Image pre-processing (layer stacking, haze, and noise reduction), classification using the maximum likelihood classification (MLC) algorithm, and post-processing were performed in ERDAS imagines 2016 (Leta et al., 2021). Based on the analyst’s experience/knowledge of the area, seven LULC classes were identified (Table 3) and statistically validated using equations (Kafy et al., 2021). Ground truthing points (GTPs) collected from Google Earth (485) were used to assess the accuracy of the 1987 and 2003 LULC classification, whereas the GTPs collected during field surveys (281) were used to assess the accuracy of the 2020 LULC classification (Garedew et al., 2009; Leta et al., 2021).

Table 3. Major LULC types and their descriptions (Bewket & Teferi, 2009; Gelagay & Minale, 2016)

LULC cover type	Descriptions
Water body	Covered by water and it may be reservoirs, lakes, large water harvesting ponds, etc.
Natural forest	Area covered by dense and tall trees both natural and plantations
Bare land	Areas that have little or no vegetation cover, mainly with classic gullies and exposed rocks
Shrub land	Areas covered with short shrubs and thorny bushes, usually stony with a very rugged relief
Banana land	Areas covered with banana plantation
Cultivation land	Land covered by annual and perennial crops and fallow lands
Grazing land	Areas with grass cover

The 2050 LULC was predicted using Multi-Layer Perceptron Neural Network (MLP_NN) and Cellular Automata-Markov Chain (CA-MC) model in the TerrSet 18.32 (Girma et al., 2022). Details of MLP_NN and CA-MC model description and procedure used are given in Girma et al. (2022). In this study, slope, elevation, population density, and distance from stream and road were used as driver factors. Kappa statistic index was applied to validate the skill of MLP_NN and CA-MC models using the actual (2020) and predicted (2020) LULC (Girma et al., 2022; Leta et al., 2021). Finally, a long-term LULC change detection analysis on total gain or loss, percent change and rate of change (ha/yr) was performed (Dibaba et al., 2020; Girma et al., 2022).

2.4. Statistical validation of regional climate models

Distribution mapping (DT) bias correction method was employed to adjust the systematic errors in RCMs using Climate Model data for hydrologic modeling (CMhyd) software. Moreover, the performances of the ensemble average RCM model were evaluated using three statistical measures (Equation 1(1) $\mathrm{N}\mathrm{S}\mathrm{E}=1-\frac{{\sum }_{\mathrm{t}=1}^{\mathrm{N}}\hspace{0.17em}{\left({\mathrm{R}}^{\mathrm{O}\mathrm{b}\mathrm{s}\mathrm{e}\mathrm{r}\mathrm{v}}-{\mathrm{R}}^{\mathrm{O}\mathrm{b}\mathrm{s}\mathrm{e}\mathrm{r}\mathrm{v}}\right)}^2}{{\sum }_{\mathrm{t}=1}^{\mathrm{N}}({\mathrm{R}}^{\mathrm{O}\mathrm{b}\mathrm{s}\mathrm{e}\mathrm{r}\mathrm{v}}-\overline{{\mathrm{R}}^{\mathrm{O}\mathrm{b}\mathrm{s}\mathrm{e}\mathrm{r}\mathrm{v}}}{)}^2}$(1) -3(3) $RMSE=\sqrt{\frac{1}{N}{\sum }_{i=1}^N{\left(R^{RCM}-R^{Observ}\right)}^2}$(3) ): correlation coefficient (R²), root-mean-square error (RMSE), and Nash–Sutcliffe measure of efficiency (NSE) (Girmay et al., 2021). The correlation coefficient (R²) is often used to evaluate the linear relationship between mean monthly RCM rainfall and observed rainfall. Values close to 1.0 indicate a better linear relationship between the variables and a value away from 1.0 indicates less agreement. Root Mean Square Error measures how much error there is between two data sets. An RMSE value close to zero indicates better performance (Worku et al., 2018). Nash–Sutcliffe efficiency indicates how well the observed versus simulated data fits the 1:1 line. NSE = 1 corresponds to a perfect match of the model to the observed data (Nash & Sutcliffe, 1970).

(1) $\mathrm{N}\mathrm{S}\mathrm{E}=1-\frac{{\sum }_{\mathrm{t}=1}^{\mathrm{N}}\hspace{0.17em}{\left({\mathrm{R}}^{\mathrm{O}\mathrm{b}\mathrm{s}\mathrm{e}\mathrm{r}\mathrm{v}}-{\mathrm{R}}^{\mathrm{O}\mathrm{b}\mathrm{s}\mathrm{e}\mathrm{r}\mathrm{v}}\right)}^2}{{\sum }_{\mathrm{t}=1}^{\mathrm{N}}({\mathrm{R}}^{\mathrm{O}\mathrm{b}\mathrm{s}\mathrm{e}\mathrm{r}\mathrm{v}}-\overline{{\mathrm{R}}^{\mathrm{O}\mathrm{b}\mathrm{s}\mathrm{e}\mathrm{r}\mathrm{v}}}{)}^2}$(1)

(2) ${\mathrm{R}}^2=\frac{{\sum }_{\mathrm{t}=1}^{\mathrm{N}}{\left({\mathrm{R}}^{\mathrm{R}\mathrm{C}\mathrm{M}}-\overline{{\mathrm{R}}^{\mathrm{R}\mathrm{C}\mathrm{M}}}\right)}^2\ast \left({\mathrm{R}}^{\mathrm{O}\mathrm{b}\mathrm{s}\mathrm{e}\mathrm{r}\mathrm{v}}-\overline{{\mathrm{R}}^{\mathrm{O}\mathrm{b}\mathrm{s}\mathrm{e}\mathrm{r}\mathrm{v}}}\right)}{{\sum }_{\mathrm{t}=1}^{\mathrm{N}}\left({\mathrm{R}}^{\mathrm{R}\mathrm{C}\mathrm{M}}-\overline{{\mathrm{R}}^{\mathrm{R}\mathrm{C}\mathrm{M}}}\right){\sum }_{\mathrm{t}=1}^{\mathrm{N}}\left({\mathrm{R}}^{\mathrm{O}\mathrm{b}\mathrm{s}\mathrm{e}\mathrm{r}\mathrm{v}}-\overline{{\mathrm{R}}^{\mathrm{O}\mathrm{b}\mathrm{s}\mathrm{e}\mathrm{r}\mathrm{v}}}\right)}\hspace{0.17em}$(2)

(3) $RMSE=\sqrt{\frac{1}{N}{\sum }_{i=1}^N{\left(R^{RCM}-R^{Observ}\right)}^2}$(3)

where, $R^{RCM}$ is rainfall of RCMs, ${\mathrm{R}}^{\mathrm{O}\mathrm{b}\mathrm{s}\mathrm{e}\mathrm{r}\mathrm{v}}$ is observed rainfall, the bar over the variables denotes the average throughout 2006–2020, and N represents the analysis period (15 years).

2.5. Soil loss estimation

The study integrated GIS with RUSLE to estimate the soil loss. RUSLE is empirically expressed by Equation 4(4) $A=R\ast K\ast L\ast S\ast C\ast P$(4) (McCool et al., 1995; Meshesha et al., 2012; Wolka et al., 2015) and considered the following input variable: rainfall–runoff erosivity (R-factor), soil characteristics (K-factor), topography (LS-factor), land cover management (C-factor), and support practice (P-factor). In order to estimate the future soil erosion, the values of C, P, and R factors were changed while the others kept constant.

(4) $A=R\ast K\ast L\ast S\ast C\ast P$(4)

where, A is the amount of soil erosion (t-ha-¹ year-¹), R is a rainfall-runoff erosivity factor (MJ mm ha⁻¹ hr ⁻¹year ⁻¹), K is a soil erodibility factor (t-hr ⁻¹ MJ ⁻¹ mm⁻¹), LS is slope length and steepness factor (dimensionless), C is a cover management factor (dimensionless), P is support practice factor (dimensionless);

2.5.1. Rainfall-runoff erosivity factor (R)

The R value was calculated using equation 5(5) $R=-8.12+\left(0.562\ast P\right)$(5) which was developed by Hurni (1985) for Ethiopia and used by several authors (Belay & Mengistu, 2021; Yesuph & Dagnew, 2019). The observed point data was transformed into a spatial dataset using inverse distance weight (IDW) interpolation algorithm in ArcGIS 10.8. The baseline erosivity factor maps were prepared for the three reference years (1987, 2003, and 2020). Therefore, the annual rainfall data of the year 1987 from four sations, the mean annual rainfall data from 1987 to 2003, and the mean annual rainfall data from 1987 to 2020 were used to estimate the R factor value for the years 1987, 2003, and 2020, respectively.

(5) $R=-8.12+\left(0.562\ast P\right)$(5)

where, P is the mean annual rainfall in (mm)

2.5.2. Soil erodibility factor (K)

The K-factor value was estimated using Equation 6 (Williams, 1995) and raster map was generated.

(6) $K=fcsand\times fcl-si\times forgC\times fhisand$(6)

where, fcsand is a factor that gives low soil erodibility factors for soils with high coarse-sand contents and high values for soils with little sand, fcl-si is a factor that gives low soil erodibility factors for soils with high clay to silt ratios, forgC is a factor that reduces soil erodibility for soils with high organic carbon content, and fhisand is a factor that reduces soil erodibility for soils with extremely high sand contents. The different factors indicated in equation 6 are determined using equations 7(7) $fcsand=(0.2+0.3\ast Exp\left[-0.256\ast ms\ast \left(1+\frac{msilt}{100}\right)\right]$(7) -10(10) $fhisand=\left(1-\frac{0.7\ast \left(1-\frac{ms}{100}\right)}{\left(1-\frac{ms}{100}\right)+Exp\left[5.51+22.9\left(1-\frac{ms}{100}\right)\right]}\right)$(10) (Neitsch et al., 2002).

(7) $fcsand=(0.2+0.3\ast Exp\left[-0.256\ast ms\ast \left(1+\frac{msilt}{100}\right)\right]$(7)

(8) $fcl-si={\left(\frac{msilt}{mc+msilt}\right)}^{0.3}$(8)

(9) $forgC=\left(1-\frac{0.256\ast \mathrm{o}\mathrm{r}\mathrm{g}\mathrm{C}}{\mathrm{O}\mathrm{r}\mathrm{g}\mathrm{C}+\mathrm{E}\mathrm{x}\mathrm{p}\left[3.72-2.96\ast \mathrm{o}\mathrm{r}\mathrm{g}\mathrm{C}\right]}\right)$(9)

(10) $fhisand=\left(1-\frac{0.7\ast \left(1-\frac{ms}{100}\right)}{\left(1-\frac{ms}{100}\right)+Exp\left[5.51+22.9\left(1-\frac{ms}{100}\right)\right]}\right)$(10)

where, ms the percent is the sand content (0.05–2.00 mm diameter particles), msilt is the percent silt content (0.002–0.05 mm diameter particles), mc is the percent clay content (<0.002 mm diameter particles), and orgC is the percent organic carbon content of the layer (%).

2.5.3. Topography factor (LS)

The LS factor value was determined by using equation 11(11) $LS={\left(\frac{\beta x}{22.13}\right)}^{0.5}\ast {\left(\frac{sin\left(\theta \right)}{0.0896}\right)}^{1.3}$(11) developed by Moore and Wilson (1992).

(11) $LS={\left(\frac{\beta x}{22.13}\right)}^{0.5}\ast {\left(\frac{sin\left(\theta \right)}{0.0896}\right)}^{1.3}$(11)

where, β is flow accumulation, χ is grid cell size (30), 22.13 is the RUSLE standard plot length, $\theta$ is the slope in degrees (i.e., slope of DEM × 0.01745), and 0.5 is the exponent of slope length.

2.5.4. Cover management factor (C)

The derived thematic LULC raster map was imported in ArcMap to assign the corresponding C-factor values adopted from different literature (Table 4) and to generate a raster map of the C-factor.

Table 4. Land cover classes and assigned cover (C) factor values

LULC class	C-values	References
Banana Land	0.089	Angima et al. (2003)
Barren Land	1	Yesuph and Dagnew (2019); Girma and Gebre (2020)
Cultivation Land	0.15	Yesuph and Dagnew (2019); Girma and Gebre (2020)
Grazing Land	0.05	Yesuph and Dagnew (2019); Girma and Gebre (2020)
Forest Land	0.01	Hurni (1985); Girma and Gebre (2020)
Shrub Land	0.014	Wischmeier and Smith (1978);Girma and Gebre (2020)
Water Body	0.0	Girma and Gebre (2020)

2.5.5. Support practice factor (P)

In this study, we used the P-values suggested by Wischmeier and Smith (1978), which consider only two types of land uses categories, agricultural and non-agricultural, to determine the P-value. Thus, the agricultural lands are classified into six slope categories and assigned P-values, while all non-agricultural lands were assigned a P-value of 1 (Table 5).

Table 5. LULC and slope classes with assigned P-values

Land use type	Slope classes in percent (%)	P-factor
Agricultural land	0–5	0.1
	5–10	0.12
	10–20	0.14
	20–30	0.19
	30–50	0.25
	>50	0.33
Other land use	All	1

Source: Adapted from Wischmeier and Smith (1978).

Finally, the estimated annual soil loss rate was categorized into five classes: very slight (0–5 t ha yr⁻¹), slight (5–11 t ha yr⁻¹), moderate (11–30 t ha yr⁻¹), severe (30–50 t ha yr⁻¹), severe (45–60 t ha yr⁻¹), and very severe (>50 t ha yr⁻¹) to identify priority areas for landscape restoration (Gashaw et al., 2020). Figure 2 shows the overall methodological framework.

Figure 2. Overall study framework used to estimate soil erosion.

3. Result and discussion

3.1. Land use and land cover

The user, producer, and overall accuracies for the different land uses are given in and kappa analysis Table 6 with respect to each year. The result indicates that it is more than the substantial range (≥85%) ensuring a strong level of acceptance (Girma et al., 2022).

Table 6. Classification accuracy assessment result for 1987, 2003, and 2020

Land use type	1987		2003		2020
	PA	UA	PA	UA	PA	UA
Water body	100	100	100	100	100	80
Cultivation	97.6	95.38	98.6	97.33	100	93.89
Barren land	91.7	91.67	93.3	90.32	100	100
Banana land	90.5	90.48	93	91.01	96.29	100
Shrub land	96.2	100	97	97.01	94.74	100
Grazing land	97	98.89	94.9	98.95	90.48	92.68
Forest land	95	96.04	94.6	96.43	88.68	100
Overall accuracy (%)	96.03		96.49		95.73
Kappa Statistics (%)	0.81		0.80		0.94

PA, producer accuracy; UA, user accuracy.

In 1987, the dominant land cover type was shrubland followed by cultivation land and grazing land. The highland part was dominated by forestlands, the midland part of the watershed was dominantly covered by cultivation land, while the low-lying areas were dominated by shrublands Figure 3a. Banana cover was observed on the lower part cloth to the outlet of the watershed. In 2003 and 2020, cultivation land became a major LULC class followed by shrubland and forestland Figure 3(b,c). During these years, bare land was significantly expanded in the middle part of the watershed, while the banana land cover was also increased in the lower part of the watershed. Over the last three decades, Sile River watershed showed massive LULC changes, particularly, the expansion of cultivable lands at the expense of natural vegetation (forestland, shrubland, and grazing land). Moreover, the expansion of banana land cover and bare land was also observed. The result is consistent with previous studies conducted in Ethiopia. For example, Belihu et al. (2020), Dibaba et al. (2020), Regasa et al. (2021), and Girma et al. (2022) indicated that cultivation land and bare land were increased during the last decades in a dramatic manner, while the area covered by forest, grazing land, and shrubland decreased. The ever-increasing human population and the subsequent increased demands on land and water resources are the major driving factors (Dibaba et al., 2020; Regasa et al., 2021; WoldeYohannes et al., 2018).

Figure 3. LULC map (a) 1987, (b) 2003, (c) 2020, and (d) 2050 in Sile watershed.

The 2050 LULC with business as usual scenarios (BAU) were analyzed. The kappa statistics index value (>0.8) indicated the better performance of the MLP_NN and CA-MC in simulating future (2050) LULC (Koko et al., 2020; Leta et al., 2021). As shown in Table 7 and Figure 3d, cultivation land, barren land, and the banana LULC were increased from 12,793 ha to 14,562 ha, 4,133 ha to 4,791 ha, and 1,257 ha to 1,631 ha, respectively. However, shrubland, forestland, and grazing land will decrease by 2050. Similarly, studies also revealed the expectation of these trends in the future (Girma et al., 2022; Regasa et al., 2021).

Table 7. Expected change of future LULC (2050) in study area

LULC category	LULC_2020		LULC_2050		Change (2020–2050)		Annual change
	Area (ha)	%	Area (ha)	(%)	Area (ha)	%	Area	%
Water body	208	1	198	1	−10	−5	−0.35	−0.17
Cultivation	12,793	46	14,562	53	1,769	14	59	0.46
Barren land	4,133	15	4,791	17	658	16	22	0.53
Banana land	1,257	5	1,631	6	629	50	21	1.67
Shrub land	6,539	24	4,711	17	−1,953	−30	−65	−1.00
Grazing land	1,386	5	990	4	−397	−29	−13	−0.95
Forest land	1,295	5	728	3	−697	−54	−23	−1.79
Total	27,610		27,610

3.2. Performance of RCMs

The correlation of RCP 4.5 and RCP 8.5 ensemble average rainfall data with observed rainfall data ranged from 0.699 to 0.794 and 0.50 to 0.802 from different stations, respectively (Table 8). In both the scenarios (RCP 4.5 and RCP 8.5) the ensemble mean of rainfall data was better correlated with the observed stations. The result of Nash–Sutcliffe efficiency (NSE) and root mean square error (RMSE) in both scenarios also confirmed that there was a better relationship in performance between the observed and simulated rainfall. The ensemble average performed better than individual simulated RCM models (Belay & Mengistu, 2021; Worku et al., 2018). Therefore, in this study, four RCM models’ ensemble mean under two scenarios was computed and used to generate an R factor value for 2050.

Table 8. Ensemble RCM performance evaluation

Stations	NSE		RMSE		R^2
	RCP4.5	RCP8.5	RCP4.5	RCP8.5	RCP4.5	RCP8.5
Arbaminch	0.686	0.589	20.684	26.382	0.794	0.50
M/Abaya	0.627	0.756	12.647	12.271	0.723	0.802
Gidole	0.691	0.576	14.245	23.333	0.699	0.646
Gato	0.599	0.745	16.961	14.804	0.764	0.777

3.3. Spatial distribution of the RUSLE factors

It was observed that the highest R value was recorded in the upper western part of the watershed and lowest towards the low-lying areas in the east (Figure 4). The erosivity (R) factor value showed an increasing trend from 1987 to 2050. The result indicated that the spatial distribution of R varied due to the spatial distribution of the long-term mean annual rainfall. The relative higher R value is associated with more rainfall power to erode the soil and thus higher vulnerability to erosion (Kebede et al., 2021).

Figure 4. R factor map (A) 1987, (B) 2003, (C) 2020, (D) 2050 under RCP 4.5, and (E) 2050 under RCP 8.5.

Five soil textural classes were identified, the upper part is dominated by clay, clay loam, and sandy textural class, whereas the lower and middle parts of the watershed were dominated by silty clay loam and sandy clay loam. The soil organic carbon content ranged from 0.128% to 3.9%. The soil erodibility (K) value ranged from 0.1 to 0.18 t-ha⁻¹ MJ⁻¹ mm⁻¹. The highest K-factor value of soil is associated with loamy (very fine sand with silt) particle, thus highly susceptible to detachment and erosion (Kebede et al., 2021). The LS value varied from 0 in the lower part of the watershed to 54 in the steepest and upper part of the watershed (Figure 5). The higher the LS value, the higher the area is susceptible to soil erosion by water (Girma & Gebre, 2020).

Figure 5. K factor (A) and LS factor (B) map of Sile watershed.

The cover management (C) factor value map is generated in Figure 6; the value varied from 0 to 1. The C value is considered as dynamic factors that vary over time (Kavian et al., 2017; Woldesenbet et al., 2020). Therefore, four different C-factor maps were generated for the year 1987, 2003, 2020, and 2050 (Figure 6). Lower C values were observed in the low-lying areas of the watershed. On the contrary, higher C values can be seen in the central and upper parts of the catchment where a large portion of the land is characterized by farmlands and bare soil. The higher C value areas were more prone to water erosion and hence greater soil loss is expected (Olorunfemi et al., 2020).

Figure 6. Cover (C) factor map (A) 1987, (B) 2003, (C) 2020, and (D) 2050.

The conservation practice (P) factor is controlled by the correlation between the slope in percent and LULC of the study area. The highest P-value of the watershed was found in the upstream part of the study area, but the smallest value was indicated in the downstream and midstream parts of the watershed. Four different P-factor maps were prepared for the year 1987, 2003, 2020, and 2050 (Figure 7). The results indicated that the distribution of the P-factor ranges from 0.1 to 1. The lower P value indicates relatively better soil erosion control measures.

Figure 7. Practice (P) factor map for 1987 (A), 2003 (B), 2020 (C), and 2050 (D).

3.4. Soil loss estimation

The average annual soil loss (A) was calculated by overlaying the preceding five maps using RUSLE with the help of ArcGIS. The result indicated that the mean annual soil loss for the year 1987, 2003, and 2020 were 13.05 t ha⁻¹ year⁻¹, 21.04 t ha⁻¹ year⁻¹, and 41.41 t ha⁻¹ year⁻¹, respectively. The total annual soil loss in the study area was 389,456 t year⁻¹, 616,669 t year⁻¹, and 1,184,920 t year⁻¹ for the year 1987, 2003, and 2020, respectively. The findings are consistent with previous studies conducted in different parts of Ethiopia (Gezahegn & Arus, 2020; Kebede et al., 2021; Molla & Sisheber, 2017). The predicted (2050) total annual soil loss under RCP4.5 and RCP8.5 was 1,469,375 t year⁻¹ and 1,485,810 t year⁻¹ with averge estimates of 56.48 t ha⁻¹year⁻¹ and 57.11 t ha⁻¹ year⁻¹, respectively. In order to identify hotspot areas, soil loss values were categorized into five soil erosion severity classes (Figure 8). Most of the upper and central parts of the watershed were characterized by high erosion risk category (>30 t ha⁻¹ year⁻¹) and susceptibility to water erosion. These areas generate considerable runoff and are characterized as areas of high erosion hotspots that require principal soil and water conservation measures.

Figure 8. The mean annual soil loss: (A) 1987, (B) 2003, (C) 2020, (D) 2050-RCP4.5, and (E) 2050-RCP8.5.

The analysis showed the increasing trends of soil loss in the past and will also be extended in the future under both scenarios. For example, the mean annual soil loss from 1987 to 2003 and 2003 to 2020 was augmented by 61% and 97%, respectively. The predicted mean annual soil loss under RCP4.5 and RCP8.5 will increase by 15.07 t ha⁻¹ year⁻¹ (36%) and 15.7 t ha⁻¹ year⁻¹ (38%) as compared to the baseline period (2020), respectively. These results were attributed to the temporal variation of climate variables (i.e. increased in R value) and detrimental changes in LULC (C and P factor values). According to previous studies (Girmay et al., 2021; Negese et al., 2021), LULC alteration and climate change played a significant role in the spatiotemporal variability of the soil erosion rate. Table 9 describes the estimated soil loss per each LULC type. Accordingly, very severe soil loss rates were recorded in bare land and cultivation land, moderate mean annual soil loss was estimated from grazing land, while the lowest average annual soil loss was contributed from banana land cover, forest, and shrubland. Thus, besides its economic feasibility (Ambisa et al., 2019), the banana land cover played a vital role in mitigating the soil erosion risk in the study watershed.

Table 9. Soil erosion under different LULC classes in the study watershed

LULC	Soil loss 1987		Soil loss 2003		Soil loss 2020
	Annual loss (t yr^−1)	Mean (t ha^−1 yr^−1)	Annual loss (t yr^−1)	Mean (t ha^−1 yr^−1)	Annual loss (t yr^−1)	Mean (t ha^−1 yr^−1)
BC	1,028	2	2,626	3.5	2,639	1.9
GL	14,698	15.4	38,067	15.2	38,153	20
BL	31,571	27.9	214,178	35.5	534,664	45.6
CL	257,224	23.3	270,126	29.3	545,637	41.7
SL	51,556	5.1	73,559	7.9	48,466	9.3
FL	33,375	8.6	18,110	4.7	15,361	8.1
Total	389,456		616,669		1,184,920

CL = Cultivation Land, NF = Forest Land, GL = Grazing Land, BC = Banana Cover, BL = Bare Land, SL = Shrub Land.

To identify the spatial variation, extent of soil loss, and prioritize for intervention planning, the study watershed was subdivided into 10 sub-watersheds (Figure 9) and categorized into five erosion severity classes (Table 10). Around 63% (most of the lower areas) of the watershed was found to be between very slight and slight erosion risk. On the contrary, 29% of the Sile watershed was experiencing severe and very severe soil losses. Moreover, the highest amount of soil loss (42%) was contributed from subwatershed numbers 4 and 5. Therefore, these subwatersheds were the first priority for watershed management planning and implementation. Moreover, the results when compared with the tolerable soil loss value in Ethiopia (Hurni, 1985), 16–22 t ha⁻¹ year⁻¹ all sub watershed exceeded the level in different extents. Over 46% area of subwatersheds 1–6 are affected by high soil erosion (>11 t ha⁻¹ year⁻¹) demanding interventions for sustainability of the soil resource.

Figure 9. Estimated soil erosion severity classes based on sub-watershed level in Sile watershed.

Table 10. Soil erosion severity in sub-watershed level and conservation priority class

Sub WS	Severity class (t-ha^−1 yr^−1) and its area coverage (ha)					Total soil loss (t/yr^−1)	%	Priority
	0–5 (Very slight)	5–11 (Slight)	11–30 (Moderate)	30–50 (Severe)	>50 (Very severe)
Sub_1	1,252	103	322	279	976	126,591	11	4
Sub_2	776	53	198	184	508	84,656	7	7
Sub_3	1,020	120	189	188	721	123,159	10	5
Sub_4	1,817	270	500	358	941	247,264	21	2
Sub_5	833	92	213	196	793	258,282	22	1
Sub_6	647	62	78	93	437	75,450	6	8
Sub_7	3,005	176	139	87	323	47,300	4	9
Sub_8	2,634	187	251	253	908	126,988	11	3
Sub_9	2,481	196	110	38	109	9,364	1	10
Sub_10	1,432	156	142	170	594	85,866	7	6
Total (ha)	15,897	1,414	2,143	1,846	6,310	1,184,920
Total (%)	58	5	8	7	22

SW = Subwatershed, (t yr⁻¹) = ton/year, (t ha⁻¹ yr⁻¹) = ton per hectare per year.

5. Conclusion

The study was conducted in the Sile river watershed to assess the combined impact of LULC and climate changes on soil loss rate. The results revealed that the watershed has undergone significant landscape alterations and an alarming increase in soil erosion risk. From 1987 to 2020, there was a substantial shift from natural vegetation (forestland, shrubland, and grazing land) to cultivation land, along with an expansion of banana farming and barren land driven primarily by the increasing demands for land and water resources of a growing human population. Projection for 2050 under business as usual scenarios indicated further expansion of cultivation, banana, and barren lands. The results indicated an increasing trend in soil erosion rates from 1987 to 2050 with substantial annual soil loss, especially in the upper and central parts of the watershed. The findings underscored the correlation between climate change and soil erosion, as evidenced by the rising erosivity factor (R) and the estimated increase in soil loss. Notably, banana land cover appears to play a mitigating role in reducing soil erosion. The sub-watershed analysis underscored the urgent need for targeted watershed management strategies to mitigate the environmental and socio-economic consequences of soil erosion, with priority given to areas experiencing severe soil losses. Overall, addressing the challenges posed by land use and climate changes in the Sile river watershed requires a multidisciplinary and collaborative approach that integrates ecological, social, and policy dimensions to ensure sustainable land management and environmental protection. Additionally, future studies should focus on the effectiveness of different soil conservation strategies and their feasibility in mitigating soil erosion in the identified high-risk areas, taking into account more detailed socioeconomic factors.

Correction

This article has been corrected with minor changes. These changes do not impact the academic content of the article.

Acknowledgment

We are grateful to the Southern Agricultural Research Institute for financial support.

Disclosure statement

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

Additional information

Funding

The work was supported by the Southern Agricultural Research Institute.

Notes on contributors

Birhanu Wolde

The first author, Birhanu Wolde (MSc), is a soil and water conservation researcher at the Agricultural Research Center in Arbaminch, Ethiopia. He has a background in natural resource management and specializes in soil and water conservation engineering.

Awdenegest Moges

Awdenegest Moges (PhD), is an Associate Professor in the Department of Biosystems Engineering, Hawassa University, Ethiopia. His research interests are soil degradation and conservation and watershed management. He has published articles related to soil degradation and conservation, land use dynamics, and land management aspects.

Rediet Grima

Rediet Girma (MSc) is a lecturer and researcher in the Department of Biosystems Engineering, Hawassa University, Ethiopia. His research interests are modeling soil erosion and land degradation, watershed management, morphometric analysis, land use change, GIS and remote sensing, climate change, and hydrology.