Impacts of land use/land cover and climate changes on soil erosion in Muga watershed, Upper Blue Nile basin (Abay), Ethiopia

Soil erosion is one of the major threats in the Ethiopian highlands. In this study, soil erosion in the Muga watershed of the Upper Blue Nile Basin (Abay) under historical and future climate and land use/land cover (LULC) change was assessed. Future LULC was predicted based on LULC map of 1985, 2002, and 2017. LULC maps of the historical periods were delineated from Landsat images, and future LULC was predicted using the CA–Markov chain model. Precipitation for the future period was projected from six regional circulation models. The RUSLE model was used to estimate the current and future soil erosion rate in Muga watershed. The average annual rate of soil erosion in the study area was increased from about 15 t ha−1 year−1 in 1985 to 19 t ha−1 year−1 in 2002, and 19.7 t ha−1 year−1 in 2017. Expansion of crop cultivation and loss of vegetation caused an increase in soil erosion. Unless proper measure is taken against the LULC changes, the rate of soil loss is expected to increase and reach about 20.7 t ha−1 year−1 in 2033. In the 2050s, soil loss is projected to increase by 9.6% and 11.3% under RCP4.5 and RCP8.5, respectively, compared with the baseline period. Thus, the soil loss rate is expected to increase under both scenarios due to the higher erosive power of the future intense rainfall. When both LULC and climate changes act together, the mean annual soil loss rate shows a rise of 13.2% and 15.7% in the future under RCP4.5 and RCP8.5, respectively, which is due to synergistic effects. The results of this study can be useful for formulating proper land use planning and investments to mitigate the adverse effect of LULC on soil loss. Furthermore, climate change will exacerbate the existing soil erosion problem and would need for vigorous proper conservation policies and investments to mitigate the negative impacts of climate change on soil loss.


Introduction
Soil is an indispensable resource; however, it has been affected by human beings since the beginning of agriculture (Amundson et al. 2015). Land use/land cover (LULC) and climate are expected to change in the future as a result of human activities (Field and Barros 2014) and expected to influence soil (Mengistu et al. 2015;de Hipt et al. 2019;Anache et al. 2018;Berberoglu et al. 2020;Hu and Gao 2020). Soil erosion has become a significant threat to terrestrial ecosystems, and it denotes a vital environmental risk (Sun et al. 2013).
In Ethiopia, soil erosion has been recognized as a serious environmental problem (Tamene and Vlek 2008;Haregeweyn et al. 2015Haregeweyn et al. , 2017. The loss of top fertile soil by water erosion creates severe limitations to sustainable agricultural land use, which leads to reduced soil productivity and food insecurity (Hurni et al. 2015).
Soil erosion jeopardizes the sustainability of agriculture and leads to siltation of streams, lakes, dams, and reservoirs (Haregeweyn et al. 2017;Yaekob et al. 2020) and downstream ecosystems impacts (Fenta et al. 2016). As a result, it can have serious implications for irrigation agriculture and related investments in different parts of the country. Thus, soil erosion caused by water is a severe problem of the country both in its onsite (e.g., productivity loss, soil fertility loss) and off-site effects (e.g., siltation) (Demessie 2015;Wagena et al. 2016;Berihun et al. 2019;Chimdessa et al. 2019;Yaekob et al. 2020;Aneseyee et al. 2020;Belihu et al. 2020). The risk varies across different areas, depending on the watershed landscape, local climatic conditions, soil characteristics, land use, and soil and water conservation and management practices. Therefore, it is a potential threat to the national food supply.
Numerous studies in Ethiopia (e.g., Mengistu et al. 2015;Taye et al. 2018;Tadesse et al. 2017;Haregeweyn et al. 2017;Gelagay and Minale 2016;Molla and Sisheber 2017;Hassen and Assen 2018;Miheretu and Yimer 2018;Anache et al. 2018;Ebabu et al. 2019;Berihun et al. 2019;Kidane et al. 2019;Aneseyee et al. 2020) reported that changes in climate and LULC could significantly affect the intensity of soil erosion at various spatiotemporal scales. For instance, Mengistu et al. (2015), reported the influence of temperature and precipitation on biomass production and soil organic carbon, thereby increasing soil erosion in the Abay river basin. An increase in precipitation levels and intensity will lead to intensified erosion (Field and Barros 2014).
Similarly, studies on the soil erosion caused by climate and LULC change were conducted in the Upper Blue Nile river basin of Ethiopia (e.g., Steenhuis et al. 2013;Demessie 2015;Mengistu et al. 2015;Wubie et al. 2016;Haregeweyn et al. 2017;Moges and Bhat 2017;Tadesse et al. 2017;Aneseyee et al. 2020). For instance, a study by Moges and Bhat (2017) in Rib watershed, and Tadesse et al. (2017) in Yezat watershed of the Upper Blue Nile basin, Ethiopia, showed that change in LULC of an area led to soil erosion. Mengistu et al. (2015) in the Abay river basin reported that climate change leads to soil erosion. These indicate that LULC change and changing precipitation patterns have an impact on soil erosion. Thus, to better understand the effects of LULC and climate change on soil erosion at a watershed level, a holistic and multidisciplinary approach is required. Modeling future soil loss rates due to climate and LULC change is an essential step (Mullan et al. 2012).
The effects of change in LULC and climate on soil erosion can be investigated by associating field measurements of soil erosion variables with different LULC classes and climate data (Adugna and Abegaz 2016;Ebabu et al. 2019). However, measurement of soil erosion is often not possible over the required temporal and spatial scale. Thus, field studies have to be complemented by soil erosion model simulation (Giri et al. 2015;Woznicki et al. 2016). Furthermore, results from soil erosion models have been used to predict impacts of climate and LULC change on soil erosion and as a scientific basis for soil erosion control and management at different scales (Addis and Klik 2015;Serpa et al. 2015;Woldesenbet et al. 2018;Anache et al. 2018).
Some researchers have developed and used different LULC change models depending on their study objectives and backgrounds. For instance, Araya and Cabral (2010) in Setúbal and Sesimbra, Portugal; de Oliveira Barros et al. (2018) in Montes Claros, Brazil; Omrani et al. (2017) in Luxembourg; Arsanjani et al. (2013) in Tehran; Xie et al. (2007) in China. Currently, the most widely used models in LULC change monitoring and prediction in Ethiopia are cellular and Markov chain models in Ethiopia (e.g., Gidey et al. 2017;Gashaw et al. 2018;Hishe et al. 2020;Kura and Beyene 2020;Fitawok et al. 2020;Mohamed and Worku 2020), and these models are verified in the Ethiopian context. Accordingly, the incorporation of GIS, CA-Markov model, climate model, and soil erosion models have been used to assess the effects of change in LULC and climate on soil erosion.
Studies including Maeda et al. (2010) in Kenya; Plangoen et al. (2013) in Thailand; Routschek et al. (2014) in Saxony, German; Ferreira et al. (2015) in the south of Portugal; Mengistu et al. (2015) in the Upper Blue Nile basin of Ethiopia; and Perović et al. (2019) in Vranjska Valley of Serbia were conducted to evaluate the potential impacts of climate change on soil erosion using soil erosion models with different scenarios. On the other hand, Sharma et al. (2011) in India; Plangoen et al. (2013) in Thailand; Ferreira et al. (2015) in the south of Portugal;de Hipt et al. (2019) in Burkina Faso, applied soil erosion model to simulate the potential effects of LULC change on soil erosion. The combined and separate effect of the current LULC and climate change on soil erosion were evaluated by several studies (e.g., Ferreira et al. 2015;Mengistu et al. 2015;de Hipt et al. 2019;Perović et al. 2019). Most of the studies focused on the separate effects of climate and LULC changes on soil erosion in different parts of the world. Although some researches were conducted on the combined and individual impacts of climate and land-use changes on soil erosion, the context is not yet well understood in the Upper Blue Nile basin. In addition, the environmental impact of future climate and LULC are still contentious issues and unresolved problems and require further research (Simane et al. 2013;Demessie 2015). Hence, predicting the impact of future climate and LULC changes on soil erosion is very important to design appropriate land use planning and adaptation and mitigation measures under local and regional scales.
The purpose of this study was to analyze the impact of LULC and climate change on soil erosion in the Muga watershed. This study would help to plan and use appropriate land use planning, soil and water conservation practices, decision making and policy development, thereby reduce soil erosion in the study area and watersheds with similar settings. Moreover, applying the CA-Markov chain model, climate model, and RUSLE model at the watershed level is the main contribution of this study in Ethiopia.

Study area
This study was carried out at Muga watershed, located within 10° 05ʹ 00˝ N to 10° 43ʹ 48˝ N and 37° 49ʹ 12ʺ E to 38° 8ʹ 56ʺ E, in the Upper Blue Nile basin, Ethiopia (Fig. 1). The area coverage is about 423 km 2 from the top of mount Choke to the watershed's gauging station. The watershed's altitude varies from slightly over 2384 m above mean sea level (a.s.l) in the southern part to 4088 m a.s.l. The study site has two distinctive seasons: a wet season from May to October and a dry season that extends from November to April.
Based on the records from 32 years  at nearby meteorological stations, the annual rainfall depth ranges from 1020.7 to 1165.52 mm. More than 85% of the rains fall during the wet season. The rainfall variability has significant impacts on agricultural production, hydrological processes, and soil erosion (Simane et al. 2013). The mean monthly minimum and the maximum temperature are about 9.3 °C and 23.7 °C, respectively (Belay and Mengistu 2019).
The recent LULC map of the study area published in Belay and Mengistu (2019) shows that 74% of the watershed is covered by cultivated land, followed by grassland (11%) and shrub-bushland (9.8%), while forest land (2.3%) and urban areas (2.3%) covered a very small proportion of the study area. The cultivated land and urban areas showed an increasing trend between 1985 and 2017 (Belay and Mengistu 2019). Forest and shrub-bushlands are the most dominant in the upper escarpment of the watershed, while most of the cultivable lands and very small urban areas are found in the lower escarpment of the watershed.
According to the data obtained from the GIS department of the Ethiopian Ministry of Water and Energy, the geology of the study area consists of tertiary extrusive and intrusive deposits. The soil of the watershed is dominated by eutric vertisols (50.2%) and eutric cambisols (26.5%), followed by Haplic Luvisols (17.3%) and Haplic Nitisols (6%) (BCEOM 1998). Eutric Vertisols predominantly occurs in the lower parts of the study watershed.

Data sources and processing
The modeling of changing soil erosion requires past, present, and future climate and land use data as model input. Thus, LULC and climate change scenarios were simulated. The RUSLE model was used for this study to obtain the scenarios of soil erosion. The main components of the data set involved in the study were described in detail below.

Image processing and preparation of spatial drivers
To assess the impacts of LULC changes on soil erosion, LULC maps were prepared. Landsat Thematic Mapper (TM) of 1985, Landsat Enhanced Thematic Mapper Plus (ETM + ) of 2002, and Landsat Operational Landsat Imageries (OLI) of 2017 with path 169 and row 053 images were among the spatial data sets used and processed to generate LULC maps. For this study, the acquisition years were selected based on a 15-year interval to easily visualize changes in spatiotemporal LULC patterns.
Furthermore, Landsat images were selected for years that align with anticipated significant LULC changes in the study area. Accordingly, the 1985 image is indicative of conditions after the collectivization of land resources with the promotion of Agricultural Producer Cooperatives in Ethiopia including the study watershed during the Dergue regime (Crewett et al. 2008). The year 2002 represents the period aftermath of the Dergue regime, and the year was selected to evaluate the redistribution of land to farmers. Finally, to include recent changes and the study area's current biophysical status, the 2017 Landsat image was used. However, due to Landsat images' accessibility and quality from USGS archives for the study watershed, a certain discrepancy in the time interval (+ 2 year) was taken into account. The images were downloaded from the USGS website (http:// earth explo rer. usgs. gov/) ( Table 1).
The satellite imageries used in the study area were cloud-free and taken at the same season. The satellite images have gone through several pre-processing stages, such as layer stacking, mosaic, geometric correction, radiometric enhancement, and sub-setting using the Earth Resource Data Analysis System (ERDAS) imagine ® 2015 and ArcGIS ® 10.4 software. All the satellite data sets were projected to the Universal Transverse Mercator map projection system zone 37N and datum of World Geodetic System 84 (WGS84), which ensure consistency between data sets during analysis.
In this study, unsupervised and supervised image classification techniques were used to classify the Landsat images used in the study. The unsupervised classification was initially applied prior to the field survey using the visual interpretation method to differentiate various land use/cover types in the studied watershed. Unsupervised classifications were carried out using Iterative Self-Organizing Data Analysis (ISODATA) clustering algorithm. The LULC maps of the study area were produced using the pixel-based supervised image classification with the maximum likelihood classification algorithm (Congalton and Green 2019). Supervised image classification is a recommended classification approach to achieve good results when adequate training data are available for the study area (Lillesand et al. 2015;Congalton and Green 2019). First, about five LULC classes (grassland, cultivated land, shrub-bush land, forest, and urban area) were identified from the images. Then, about 300 samples (60 samples per LULC class) were randomly selected from five LULC types as training points, whereas 300 samples were used for accuracy assessment based on Lillesand et al. (2015) and Congalton and Green (2019). Besides, high-resolution data (SPOT images and Google Earth), topographic maps, and data from field observations and knowledge of the elderly were used. The accuracy assessments of all the classified images were conducted using ground control points collected from: (1) a topographic map, (2) Woody Biomass Inventory and Strategic Planning Project data, (3) knowledge of elderly people, field data, and (4) Google Earth data sets on ERDAS IMAGINE environment. Ground control points were collected from each LULC type and used as a reference for the accuracy assessment. The classification's overall accuracy for the years 1985, 2002, and 2017 is 85.7%, 87.9%, and 89%, respectively, with kappa values of 0.82, 0.86, and 0.86. The accuracy statistics and kappa coefficient values are well above the recommended values (Congalton and Green 2019).
In this study, the CA-Markov chain model was used to simulate and predict the future land cover map of 2033. Different sets of spatial and non-spatial data were used as input to the CA-Markov model. The simulation and prediction processes were performed using LULC maps of 1985, 2002, and 2017, roads and towns, distance to a river, slope, and elevation maps. The LULC maps of 1985, 2002 were mainly used to generate transitional matrix using the Markov chain process. Besides maps of a slope, road, settlement, and elevations maps were used to create transitional potential maps. The data sets were used in combination to predict future change in the LULC using the CA-Markov chain model. Potential land-use change drivers were identified through literature reviews (Bewket and Teferi 2009;Teferi et al. 2013;Gidey et al. 2017;Hishe et al. 2020) and field interviews with farmers, local farming experts, regional land bureau officials, and through spatial correlations. A digital elevation model (DEM) of 30 m resolution from the Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) was used to prepare the watershed's slope, elevation, and stream network of the study watershed. The watershed road network was also downloaded from OpenStreetMap (https:// www. opens treet map. org/), and its geometric consistency was verified in QGIS software. The population data set for the study area was obtained from the Ethiopian Central Statistical Agency (CSA). Consequently, these data sets were produced using ArcGIS software and then exported to IDRISI Selva to run the potential transition maps. The list and sources of data sets used to analyze future land use map and RUSLE model in the Muga watershed is presented in Table 2.

Climate data
Two climate data sets: baseline and modeled climate data, were used in this study. The baseline climate data were used to estimate the baseline soil loss and validate the predicted climate scenarios. The rainfall data were compiled from seven meteorological stations operating within and around the study watershed from 1985 to 2017; we obtained the data from the National Metrological Agency (EMA) of Ethiopia.
Climate models were used to quantify the relative change in the current and future climate, often used as an input to the soil erosion models. The present study used six regional climate models (RCMs) with the driving model ICHEC-EC-EARTH under CORDEX (Coordinated Regional Climate Downscaling Experiment)-Africa (Fick and Hijmans 2017;Shamir et al. 2015). COR-DEX-Africa provides projected climate outputs at a relatively higher spatial resolution (50 km). The models were: HadGEM2-ES, CSIRO-Mk3.6.0, GFDL-ESM2M, CanESM2, MIROC, and NorESM1-M, because these models are appropriate based on earlier research conducted in other parts of Ethiopia (Alemseged and Tom 2015;Teklesadik et al. 2017;Worku et al. 2018).
The representative concentration pathways (RCPs) scenarios RCP4.5 and RCP8.5 were considered for this study to drive the RUSLE model (Alemseged and Tom 2015). Policymakers usually focus on events occurring in the 2050s compared to the events in the far future (Weber 2006). Thus, the authors select climate change data from 2018 to 2050 rather than the end of the century.
The CORDEX grid points were extracted using MAT-LAB software. It is recognized that climate model output data contain systematic errors and cannot be used immediately in soil and hydrological simulations (Christensen et al. 2008;Seibert 2010, 2012;Chen et al. 2016). Thus, the RCM climate data outputs in CORDEX-Africa under emission scenarios RCP4.5 and RCP8.5 were bias-corrected for this study.

Land use land cover prediction using CA-Markov chain model
Several LULC change models were developed (e.g., Veldkamp and Lambin 2001;Parker et al. 2003;Verburg et al. 2004;Koomen and Stillwell 2007). Some of the most widely used LULC models include statistical models (regression) (Yalew et al. 2016 (Omrani et al. 2017), and multi-agent-based models (Xie et al. 2007;Arsanjani et al. 2013). These approaches are often combined to create an integrated model that determines the probabilities of LULC changes. Frequently, the choice of models depends on the study's objective and the level of complexity required (Nainggolan et al. 2012). The Markov Chain model of LULC change has been widely used to predict LU changes (Sang et al. 2011) and usually integrated with a cellular automata (CA) model with an agent-based model to express the human interaction on the landscape (Xie et al. 2007;Ralha et al. 2013). This is due to its reliability and compatibility with many geospatial technologies (Halmy et al. 2015). The reliability and compatibilities of the CA and Markov chain model was checked in the Ethiopian context by Gidey et al. (2017), Gashaw et al. (2018) and Hishe et al. (2020). Thus, this study used the CA-Markov chain model integrated with the multi-criteria evaluation to predict LULC of 2033 in the Muga watershed. The Markov chain model provides the probability and extent of LULC change, whereas the CA model operates on neighborhood interactions and spatial distribution.
The complex interrelationships between the physical and human factors lead to the conversion from one land use type to another (Eastman 2012). Thus, with the CA-Markov model, a multi-criteria evaluation (MCE) technique was used to support land allocation decision process using different land use criteria (Al-Sharif and Pradhan 2014; Deng et al. 2015). The MCE is useful for studying the land use suitability for the possible conversion of certain land use to another, revealing each criterion's influence and importance (Yalew et al. 2016;Gidey et al. 2017).

Preparation of suitability maps
LULC can be driven by a multitude of socio-economic and biophysical factors. The future change in LULC can be determined by the inherent change and external factors/spatial variables, such as the proximity to towns, distance from a river and city, elevation, slope, and areas suitable for each change in each class. A slope was considered as a constraint for cultivation and urban areas because steep slopes prohibit both cultivation and urban expansion.
An integrated evaluation procedure was used to generate potential transition maps based on biophysical and socioeconomic indicators. In this regard, more than five biophysical and socio-economic variables including elevation, slope, distance to towns, and distance to roads were considered (Additional file 1: Fig. S2). The selected transition potential maps are in different units. Therefore, the maps were converted into a uniform measurement scale through standardization techniques for weighted overlay analysis (WOA) (Reshmidevi et al. 2009;Zabihi et al. 2015). Once all the maps were standardized, the weight for each criterion was calculated using the analytical hierarchy process (AHP) (Yalew et al. 2016;Gidey et al. 2017).
Consequently, specific weights were assigned to each factor (Hishe et al. 2020) and used to compute suitable maps in IDRISI software. The relative weights for a group of factors were defined based on the authors' knowledge and experience about the studied landscape, a review of the scientific literature, and the farmers' and local experts' opinions. The highest weight value is the most influential factor, while the lowest is a less important LULC change factor. After determining each criterion layer's relative importance through a pairwise comparison matrix, these values were entered using IDRISI software to produce associated weights and consistency ratio value. Table 3 shows inputs to the pairwise comparison for the AHP analysis to determine weights. The weights produced from the AHP procedure using the inputs in Table 3 are between 0 and 1, where 0 denotes the no probability and 1 the high probability. The consistency ratio of the pairwise comparisons for the computation of criteria weights is shown in Table 3, which is an acceptable range (Eastman 2012). Consequently, a weighted overlay analysis was performed using IDRISI software.
As shown in Table 4, the 1985-2002 transition matrix and the 2002 transition potential maps were integrated to simulate LULC map of 2017. The data were integrated using the CA model to simulate LULC of 2017. The same procedure was used to predict the LULC of 2033. The transition matrix from 2002 to 2017 and transition potential map of 2017 were integrated to run the prediction for 2033.

Validation of LULC prediction model
Model calibration and validation are essential in predicting future decadal changes, where no data sets are available for predicted data accuracy (Srivastava et al. 2014;Singh et al. 2015). To validate the model, the actual LULC map of the year 2017 was compared with the 2017 map simulated by CA-Markov and based on the kappa statistics and a comparison of each simulated LULC class with the real class has been used. According to the CA-Markov model's required preparation, a calibration map for 2017 was prepared. The number of repetitions in the model was set equal to the number of years between the reference map and the map predicted by the model (15 years). The actual 2017 LULC map was then used as a reference map to compare with the 2017 simulated LULC map results. Accordingly, the analysis result has shown consistency between the reference and the simulated LULC maps of 2017.
Furthermore, an attempt was also made to examine its accuracy using the kappa coefficient. The Kappa variations used to validate the CA-Markov model for LULC change predictions in this study were computing the kappa coefficient using a validated tool in IDRIS Selva software. Kappa statistics are used for testing accuracy are the traditional Kappa (K standard ) or Kappa for no information/ability (K no ), Kappa for location (K location ), and quantity of correct cells (Kappa for quantity). Kappa variations have been strongly recommended and widely used to validate LULC change predictions (Pontius Jr 2002;Singh et al. 2015).
The Kappa coefficient for quality and location was computed. The statistics show that K no is 0.8803, K location is 0.8270, and K standard is 0.8141 (overall). The K values are greater than 0.8, showing the model accuracy to predict future LULC change. In other words, the kappa coefficient result indicates the model's ability to specify grid cell level location of future change is nearly perfect (Mondal et al. 2016). The result indicates that K values of 0.80 and above are considered strong and reasonable to make plausible future projections (Yang et al. 2014;Gidey et al. 2017). The LULC scenario was prepared for 2033, accounting for the implications of future land use on soil erosion due to the variation in C factor. Figure 2 shows  the process for simulating the future LULC using CA-Markov chain model.

Climate projection and bias correction
Climate projection Climate scenarios have served as an essential tool in climate change research in the past and will likely continue to do so in the future. There are four greenhouse gas scenarios recently published in the fifth IPCC Assessment Report (AR5), called the representative concentration pathways named based on their possible range of radiative forcing values, such as (RCPs) 2.5, 4.5, 6, and 8.5. In the present study, we selected two climate change scenario data (Representative Concentration Pathways: RCP4.5 and RCP8.5) that cover the entire range of radiative forcing from the newly available Coupled Model Intercomparison Project Phase 5 (CMIP5). The output of the RCM ensemble of the Coordinated Regional Climate Downscaling Experiment (CORDEX) for African domain projection was used as an input to the hydrological model. RCP4.5 scenario is medium-term emission scenario at the stabilization level of about 4.5 w/m 2 (about 650 ppm CO 2 equivalent) not exceeding this value by the year 2100 van Vuuren et al. 2011), supposing that all countries around the world undertake emission mitigation policies . RCP8.5 scenario is the worst-case scenario in terms of greenhouse gas emissions, with no clear climate policy. The RCP8.5 scenario shows a radiative forcing pathway leading to 8.5 w/m 2 (greater than 1370 ppm CO 2 -equivalent) in 2100 (van Vuuren et al. 2011). Climate scenario data sets consist of a historical run  and projection (2018-2050) with a spatial resolution of 0.44° based on the emission scenarios of RCP4.5 and RCP8.5. In this study, climate data up to 2050 (2018-2050) was used rather than considering the end of the century, as small and large dams (with a lifespan of 30 to 50 years) in the upper Blue Nile basin of Ethiopia, including the study watershed, has been studied, designed, and started to implement for irrigation and hydropower purposes, which will be affected by a changing climate in the short period of time; and policymakers often focus on what is happening sooner than what is happening in the distant future (Weber 2006). The CMIP5 data for two specific periods (historical and future periods) were downloaded from https:// esgf-node. llnl. gov/ proje cts/ esgf-llnl/. The grid points of the regional climate scenarios (daily precipitation and maximum and minimum air temperature data) were extracted using MATLAB software. The software was also used to convert NetCDF files into text format, taking into account the spatial references of a particular weather station. Seven stations were used to extract the downloaded data: Bichena, Debre Work, Felege Birhan, Kuy, Rob Gebeya, Yetemn, and Yetnora. The models were: HadGEM2-ES, CSIRO-Mk3.6.0, GFDL-ESM2M, CanESM2, MIROC, and NorESM1-M, as these models are appropriate based on earlier research conducted in other parts of Ethiopia (Alemseged and Tom 2015;Teklesadik et al. 2017;Worku et al. 2018) (Table 5). In this study, the six models' ensemble mean was computed and used as input for the RUSLE model.

Bias correction of regional climate models simulation
The relatively low spatial resolution of GCMs and RCMs outputs systematic errors, such as estimation of climate variables (over or under), inaccurate estimates of seasonal variations of precipitation (Christensen et al. 2008), and more wet-less days compared to observed data (Ines and Hansen 2006). As a result, the data from RCPs were bias-corrected to prevent over-or-under estimation and to ensure a realistic representation of the future climate. Because climate model output cannot be used as direct input data for hydrological simulations (Christensen et al. 2008;Seibert 2010, 2012;Chen et al. 2016), the model performance in simulating the observed precipitation and maximum and minimum temperature were assessed, and those with unacceptable time seriesbased metrics results were bias-corrected. The historical RCMs and observed data were used for bias correction. Therefore, RCM outputs (precipitation) are typically adjusted to eliminate any bias (Maraun 2012;Teutschbein and Seibert 2012).
Several bias correction techniques were developed to correct the bias, ranging from very simple to more complex methods, such as quantile mapping (QM), general quantile mapping (GQM), power transformation (PT), and linear scaling (LS) (Teutschbein and Seibert 2012). They can be classified according to the degree of their complexity and simpler methods, such as scaling factors, and more sophisticated techniques, such as probability mapping. Although the RCM climate variables' bias correction significantly improves hydrological modeling, there is a major drawback. All bias correction methods are based on the assumption of static model errors (Teutschbein and Seibert 2012). This means that the correlation algorithm for current climate conditions and its parameters are considered valid even under the changing climate conditions. According to Teutschbein and Seibert (2012) and Sisay et al. (2017), linear scaling is a simple technique and provides better efficiency in correcting the ensemble RCMs climatological biases. In this study, the linear scaling bias correction method was applied to adjust the raw ensemble of climate scenarios (rainfall) output simulation data.
The linear scaling bias correction method was selected following a review of the study by Teutschbein and Seibert (2010), which evaluated five bias correction methods for precipitation and showed that linear scaling is suitable for precipitation. The linear scaling method uses a correction factor for each month based on long-term historical RCM data and the observed data. Correction factors were applied for future climate scenarios of the RCP4.5 and RCP8.5 emissions. The bias correction was used to estimate and remove the bias correction in future RCMs output. The linear scaling bias correction method (V.1.0) Microsoft Excel file described by Shrestha (2015) was used in this study to adjust the climate model's average value to observations appropriately.
To verify the bias correction method and the improvement obtained after bias correction, the average daily observed and simulated data (before and after bias correction) were again compared with the observed data by model evaluation statistical methods, such as the correction coefficient (r) and root mean square error (RMSE). Statistical results of the bias correction show that the bias correction significantly improves the simulated data as the RMSE values decrease, the SD values closer to  Sisay et al. 2017). The two RCP scenarios (RCP4.5 and RCP8.5) were taken into account to estimate the rainfall erosivity factor.

Revised Universal Soil Loss Equation model
In the present study, we used the RUSLE model for erosion modeling. RUSLE is the revised form of the USLE model (Wischmeier and Smith 1978), with substantial improvements in predicting the yearly amount of soil loss; and it was revised by Nyssen et al. (2009) to the Ethiopian conditions. The annual soil loss (A) of the study watershed was calculated by overlaying five raster layers using the following equation: where A is the average annual soil loss per unit area (t ha −1 year −1 ), R is rainfall erosivity (MJ mm ha −1 h −1 year −1 ), K is soil erodibility (t ha −1 MJ −1 mm −1 ), LS is the slope length and steepness factor, C is a cover and management factor, and P is the conservation practice factor. The procedures and techniques we used for these factors are presented in the following sections. Rainfall erosivity factor (R-factor) indicates the input, which represents the effect of rainfall intensity on soil erosion and requires in-depth, unceasing precipitation data for its calculation (Renard and Freimund 1994;Angima et al. 2003). The rainfall erosivity is calculated as the total storm energy multiplied by the maximum 30 min intensity (Renard 1997). However, such data is not available for the study area, limiting the R factor spatiotemporal application. Empirical relationships between measured rainfall amount and R-factor values widely used as an alternative approach (Hurni 1985;Nigussie et al. 2014;Mengistu et al. 2015). This method is widely used by previous researchers in Ethiopia and other countries (e.g., Bewket  (1) A = R * K * LS * C * P In the present study, the annual rainfall erosivity was calculated according to the method of Hurni (1985) to the Ethiopian highlands as the following equation: where R is the rainfall erosivity factor (in MJ mm ha −1 h −1 year −1 ), and P is the mean annual rainfall (mm).
After accounting the R factor of the baseline (1985-2017) and future period (2018-2050) using the mathematical equation, the R factor raster was created using the Inverse Distance Weighted (IDW) interpolation method in the Geostatistical Analysis extension of Arc-GIS 10.4 software.
Soil erodibility factor (K-factor) represents soil vulnerability to erosion (Wischmeier and Smith 1978). According to the soil difference, the K value of an area is also different depending on the parameters of soil texture, organic matter content, and permeability (Renard 1997). The soil erodibility was determined from the available soil data of the Blue Nile river basin (scale 1:250,000), which was obtained from the Ministry of Water and Energy. In this study, we used the nomograph method to calculated the K value as proposed by Wischmeier et al. (1971) where f p is the particle size parameter (unitless), p om is the percent organic matter (unitless), s struc is the soil structure index (unitless), f perm is the profile-permeability class factor (unitless), p silt is the percent silt (unitless) and p clay is the percent clay (unitless).
In Eq. 3 the factor (1.292) is needed to convert K-factor from the imperial to the international system units (i.e., SI metric units) (Streile et al. 1996). The soil structure index, S struc , is 4 for blocky, platy, or massive soil, 3 for medium or coarse granular soil, 2 for fine granular soil, and 1 for very fine granular soil (Mengistu et al. 2015). The profile-permeability class factor, f perm , is 1 for very slow infiltration, 2 for slow infiltration, 3 for slow to moderate infiltration, 4 for moderate infiltration, 5 for moderate to rapid infiltration, and 6 for rapid infiltration (Streile et al. 1996). In general, Eq. 3 can help capture relative differences in erodibility between soil types and help appropriate the resistance to erosion of different soils under consideration (Ganasri and Ramesh 2016). The soil type and soil erodibility map of the watershed is presented in Fig. 3.
Topographic factor (LS-factor) refers to the effect of topography on soil erosion. The slope length (L) and slope gradient (S) factors are joined in a single index, (2) R = −8.12 + 0.562P in which f p = p silt (100 − p clay ) Belay and Mengistu Ecol Process (2021) (2033) A (t ha -1 yr -1 ) = R*K*LS*C*P Soil loss estimation LS-factor, to describe the topographic factor for soil loss. The slope length refers to the distance from the point of origin of overland flow to the point, where either the slope gradient declines enough in which sedimentation starts or the runoff water enters a well-defined channel (Renard 1997). The LS factor represents (Fig. 4) the ratio of soil loss per unit area on-site to the corresponding loss from a 22.13 m long experimental plot with a 9% slope (Renard 1997). Higher slope lengths may increase the overland flow and lead to more land surface soil erosion (Moore and Burch 1986a). Moreover, higher slope gradients also promote the runoff rate and bring out more soil erosion.
In this study, the L factor was computed following Eq. (4) proposed by Moore and Burch (1986a;. The algorithm (Eq. 5) recommended by Desmet and Govers (1996) was used to calculate the S factor. It is based on flow accumulation and slope steepness. The slope steepness (S) factor was calculated for high (> 9%) and low slope land (< 9%) from the slope angle. Finally, the L and S factor were multiplied to derive the LS factor for the study watershed using the Spatial Analysis Tool Map Algebra Calculator in the ArcGIS 10.4 environment and the spatial variability of the slope length steepness factors (Renard 1997): where L and S indicate slope length and steepness factor (dimensionless); λ is the slope length (in meter), cell size is the size of the grid cell (for this study 30 m); m is an adjustable slope length exponent to β; while β represents the rill to interrill erosion ratio; and θ is the slope angle (in degrees). (4) The study area's topography is steep with non-uniform terrain on the upper parts of the watershed, whereas the area's lower escarpment is gentle and uniform. The combined steeper and longer slopes resulted in larger cumulative runoff volumes with high velocity and erosive power.
Cover management factor (C-factor) represents the effect of vegetation and the cover on the amount of soil erosion (Bewket and Teferi 2009;Haregeweyn et al. 2017). The individual values of C vary between 0 for a completely non-erodible condition, and 1 implies conditions more erodible than those normally experienced under unit plot conditions, which can occur for conditions with very extensive tillage, and C is strongly related to land use.
In Ethiopia, satellite-derived images are a good proxy for land cover on relatively large basins and watershed levels and were applied in the Upper Blue Nile basin of Ethiopia (Bewket and Teferi 2009;Mengistu et al. 2015;Molla and Sisheber 2017;Taye et al. 2018). In this study, C values were determined based on the most recent (2017) and predicted future (2033) land use data, as suggested in the literature (Table 6).
Conservation practice (P-factor) reflects the effects of support practices, such as terracing and contour tillage, to reduce the rate of soil erosion (Renard 1997). The higher the P factor, the less effectively the practice facilities deposition to take place close to the source. It is primarily used to assess the effects of conservation measures implemented on soil loss (Mengistu et al. 2015). Our field visits indicated that soil and water conservation activities are not widely practiced in the study area. However, farmers commonly use stone bunds and contour farming to protect the soil from erosion. Still, they are poorly maintained as implementation was carried out in a topdown approach, and there was no map of conserved areas in the sub-watershed. In areas where contour plowing is widely practiced, Hurni (1985) suggested using P values of 0.9 and 0.8 for agricultural and non-agricultural lands, respectively. Thus, we used a similar method in the present study to assign P values to the study area for the baseline period (2017) and future scenario (2033).

Table 6 C-factor and P-factor values for the respective land use/land cover classes of Muga watershed, Upper Blue Nile Basin, Ethiopia
Land use/land cover classes C factor P factor Grassland 0.05 (Hurni 1985;Bewket and Teferi 2009) 0.8 (Hurni 1985;Adugna et al. 2015) Cultivated land 0.15 (Hurni 1985;Adugna et al. 2015) 0.9 (Bewket and Teferi 2009;Adugna et al. 2015;Taye et al. 2018) Shrub-bushland 0.02 (Mengistu et al. 2015) 0.8 (Hurni 1985) Forest 0.05 (Molla and Sisheber 2017) 0.8 (Hurni 1985) Urban area 0.05 (Moges and Bhat 2017) 0.8 (Hurni 1985) Finally, the annual soil loss was estimated on a cellby-cell basis of multiplying the five RUSLE factors using Eq. (1). As Landsat images and ASTER-DEM used in this study had 30 m spatial resolution, all the raster maps were resampled to 30 × 30 m cell size and re-projected to UTM Zone 37° N, WGS 1984 datum. Table 6 shows the C values and the P values derived from LULC information for the baseline and future period. The workflow of the methodology used in this study is shown in Fig. 5

Land use/land cover change in Muga watershed
The produced LULC maps of the Muga watershed for the three reference years (1985, 2002, and 2017) are presented in Fig. 6. The trend analysis made for the two consecutive periods, 1985-2002 and 2002-2017, showed spatiotemporal changes in LULC classes.
The magnitude and extent of LULC change for the study watershed between 1985 and 2017 were obtained from Belay and Mengistu (2019). The LULC change between 2017 and 2033 were calculated (Table 7). The predicted LULC map for 2033 is shown in Additional file 1. As shown in LULC map of 2033, cultivated land remains a dominant land use type, which accounts for about 76.8% of the study watershed. From 2017 to 2033, the areal coverage of grassland decreased by 35.4%, while urban area, shrub-bushland, forest, and cultivation area would increase by 2.9%, 10.4%, 6.8% 3.8%, respectively. It should be noted that the decline of grassland area from 2017 to 2033 is probably due to the conversion of grasslands to cultivated land and shrub-bushland, as shown in Table 7.

Impacts of land use/land cover changes on soil erosion
LULC change in the area may increase or decrease soil erosion. The LULC dynamics of the Muga watershed from 1985 to 2017 was studied by Belay and Mengistu (2019). In the present study, to assess the impact of LULC change on soil erosion, the values of C and R factors were changed, while the other factors (i.e., soil erodibility, conservation practice factor, and slope length and slope steepness) kept constant.
The results of the study show that the average annual rate of soil erosion in the Muga watershed was increased from approximately 15 t ha −1 year −1 in 1985 to 19 t ha −1 year −1 in 2002, and 19.7 t ha −1 year −1 in 2017 (Table 8). The study results also indicate that if the LULC changes are not managed, the soil loss rate is expected to continue in 2033, which is expected to reach 20.7 t ha −1 year −1 in 2033. Areas with soil erosion rates over 15 t ha −1 year −1 during the study years were widely distributed on the upper and steep areas of the watershed, while those with an erosion rate of less than 10 t ha −1 year −1 was mainly concentrated in the gentle areas of the lower and upper part of the watershed (Fig. 7). Increasing some types of LULC, such as cultivable land in steep areas and irrigated agriculture, will accelerate soil erosion by reducing soil cover in the future.
The validation and consistency of the model output was compared with the quantitative outputs of previous  experimental observations and similar empirical studies conducted in Ethiopia, mainly in the northwestern highlands. In addition, selected field observations were carried out. In supporting this process, the color printed model output soil erosion severity map was taken in the field to check the reality on the ground. Consequently, the estimated rate of soil loss and the spatial patterns are generally realistic as compared with the findings in the field and from the results of previous studies (Hurni 1983a(Hurni , b, 1985Bewket and Teferi 2009;Mengistu et al. 2015). The relative contribution of each LULC type to soil erosion was assessed in the study watershed, and the result showed noticeable differences among LULC types in 1985LULC types in , 2002LULC types in , 2017LULC types in , and 2033. The highest mean annual soil erosion rate was predicted on cultivated land (e.g., 18 t ha −1 year −1 in 1985, 19 t ha −1 year −1 in 2002, and 21.7 t ha −1 year −1 2017).  Although the area under cultivation is expected to decrease by 2033, the average annual rate of soil erosion is expected to be higher than 1985, 2002, and 2017. The causes of the higher soil loss estimated from cultivated land may be related to the encroachments of sloping plowing in mountainous areas of the watershed, as the expense of expansion of urban areas into cultivation. The results of this study suggest that the rate of soil loss in shrub-bushland areas will be higher than in grassland and urban area in the future (2033), which may be due to the scanty vegetation and the steepness of the area, where the shrub-bushland area is located. This indicates that in the year 2033, the area of shrub-bushland is expected to have less capability to protect the soil from erosion than in the year 2017. Thus, future LULU changes under the business as usual scenario, with increasing LULC changes in the study area (i.e., expansion of cultivated land) contributing to more soil erosion, which is anticipated to increase further with higher precipitation. It is estimated that in 2033 there will be higher rates of soil erosion in cultivated, followed by grassland and shrubbushland, which are the key areas to prevent soil erosion in the Muga watershed in the future.
As shown in Fig. 7 and Table 9, the percentage areas under low to medium-level soil erosion were decreased between 1985 and 2033. In contrast, the portion of areas under high and very high soil erosion classes was increased during the same periods, but this change is relatively low. The increase in the areas with high and very high soil erosion rates between 1985 and 2033 may be due to the conversion of some parts of the watershed with low and moderate erosion classes into the range of high and very high erosion. Between 1985 and 2017, there was an upward trend in the area occupied by severe and very severe erosion classes, while the severe erosion category is expected to decline in 2033. It was also found that the area classified as very severe erosion is uniform in 2017 and 2033.

Impacts of climate change on soil erosion
Rainfall amount and rainfall erosivity under RCP4.5 and RCP8.5 for the study area showed an increasing trend, which is anticipated to affect soil loss negatively. Increasing rainfall will lead to an increase in annual rainfall erosivity, which increases the annual rate of soil loss. Thus we analyzed the average annual rate of soil erosion from the watershed by modifying climatic data during the study period to assess the relationship between soil erosion and climate change in 2017 and 2050.
In this study, we simulated the future rate of soil erosion due to changes in climatic conditions that may increase the risk of soil and land degradation in the Muga watershed, which can also affect agricultural productivity and livelihoods of the local community. The results are shown in Table 10 and Fig. 8. The results show that the average annual rate of soil loss in the Muga watershed in 2050 under RCP4.5 and RCP8.5 will be 21.6 t ha −1 year −1 and 22.2 t ha −1 year −1 , respectively, which is equivalent to an annual soil loss of 917,900.6 t year −1 and 960,396.0 t year −1 , respectively. Therefore, the rate of soil erosion is predicted to increase by 9.7% and 12.7% under RCP4.5 and RCP8.5, respectively, compared to the reference period. This means that the mean annual soil Table 9 Annual soil erosion risk classes and area coverage under land use/land cover map of 198,2002, 2017, and 2033   erosion rate in the Muga watershed showed an increasing trend in the future compared to the reference period due to climate change.
Recognizing the soil loss response helps to know the rate of soil erosion changes with climate change, which helps to identify priority areas for implementing soil management measures in the study watershed. Increased soil erosion due to climate changes can affect the watershed's local ecosystems and can cause hydrological changes in streams originating from the watershed.
Comparing LULC and climate change, as shown in Table 11, climate changes are expected to have a higher impact on soil erosion than LULC change in the future.

Combined impacts of future land use/land cover and climate change on soil erosion
In addition to assessing the individual impacts of future LULC and climate change, the combined and synergistic effects of these changes on soil erosion were also evaluated in the Muga watershed. By replacing the two factors of rainfall erosivity (R) and cover management factor (C), and keeping all other factors constant in the RUSLE model, soil erosion was predicted using future LULC and climate change scenarios: 2033 LULC map and climate scenarios under RCP4.5 and RCP8.5 in 2050. The magnitudes and rates of mean annual loss of soil under future LULC and climate change scenarios are presented in Table 11.
The highest soil erosion rate for the study watershed is predicted when LULC change is combined with climate change under the RCP8.5 scenario, which is estimated about 22.8 t ha −1 year −1 . Similarly, the average annual rate of soil loss under RCP4.5 is estimated about 21.6 t ha −1 year −1 . Therefore, changes in LULC under RCP4.5 and RCP8.5 scenarios showed a slight variation in annual soil loss rates. It should also be noted that the ranges of relative change in average annual soil loss rates for the future period due to climate change are larger than LULC change. Table 11 shows that soil erosion in the Muga watershed appears to be more sensitive to future climate changes than future LULC changes. The magnitude of the soil loss rates is expected to increase in the future due to the synergy effects of climate and LULC changes. As it is shown in Table 11, the combined effects of LULC and climate change on soil loss rate are expected to increase in the future period.
LULC changes are expected to exacerbate the rate of soil loss by 5%, while climate change is also predicted to increase the rate of soil loss by approximately 9.7% and 12.7% under RCP4.5 and RCP8.5, respectively. When LULC and climate change act together, the average annual soil loss rate increases to 13.2% and 15.7% under RCP4.5 and RCP8.5, respectively, which is higher than the individual effects of LULC and climate change. Thus, the combination of modeled LULC and climate change is expected to have a substantial impact on soil erosion in the future. Increasing average annual soil erosion rates caused by LULC and climate change will significantly impact land use planning and implementation in the Muga watershed. Besides, as the LULC types change, adverse effects may be intensified.

Discussion
Soil erosion prevention patterns varied considerably depending on location and time. Globally, as environmental change accelerates, more frequent and intense changes in LULC associated with a higher frequency of extreme climate events will increase soil degradation and make the ecosystem less resilient to natural disturbance (Cramer et al. 2018;Guerra et al. 2020). The results of the study show a spatial and temporal variability in soil erosion rate due to climate dynamics and changes in LULC. According to Zhou et al. (2008), a soil vegetation cover of more than 78% greatly reduce erosion by water. The results of this study show that less than 5% of the study watershed is covered with forest; hence, the study watershed is sensitive to soil erosion. Furthermore, the expansion of cropland at the expense of forest and shrubland reduced the protection and soil organic matter of the soil and exposed it to the impact of climate change thereby accelerated soil erosion rate (Wijitkosum 2012).
The results of future climate projection showed that the mean annual precipitation in the watershed is expected to increase by 17.15% and 20.27% under RCP4.5 and RCP8.5, respectively, compared to the historical period of 1985-2017. This result agrees with previous studies (e.g., Conway and Schipper 2011; Kassie et al. 2014;Abera et al. 2018;Fentaw et al. 2018) that predicted the increase in precipitation in the other parts of Ethiopia.
Muga watershed also showed significant dynamics in LULC from 1985 to 2017. The results of the study show that LULC changes played a significant role in increasing the soil erosion rate in the study watershed. This finding coincides with previous research conducted in the Upper Blue Nile basin (e.g., Gessesse et al. 2015;Gelagay and Minale 2016;Moges and Bhat 2017;Tadesse et al. 2017), who reported that LULC change was a responsible factor for the increased soil erosion in their respective study areas. Table 8 shows the estimated soil loss in 1985, 2002, 2018, and 2033. The estimated average soil loss was increased due to LULC dynamics that occurred between 1985 and 2017. This shows that LULC has a significant impact on soil loss by water erosion. The results of this study show that the average annual rate of soil loss in the Muga watershed (19.0 t ha −1 year −1 in 2002, 19.7 t ha −1 year −1 in 2017, and 20.7 t ha −1 year −1 in 2033) was substantially higher than the average soil loss of Upper Blue Nile basin (16 t h −1 year −1 ) reported by Mengistu et al. (2015), and Fenta et al. (2021) who reported 16.5 t ha −1 year −1 average soil loss rate for Ethiopia.
This study also shows that about 25% of the study watershed has a soil loss rate of 20 t ha −1 year −1 and above, which is higher than the tolerable soil loss limits estimated for Ethiopia. Tolerable soil loss rates suggested for Ethiopia is 12 t ha −1 year −1 (Hurni 1983b). According to Khosrokhani and Pradhan (2014), the rate of soil formation in tropical areas is generally slow, and soil loss of more than 1 t ha −1 year −1 is regarded as irreversible soil erosion. However, soil loss of 1 t ha −1 year −1 or less is considered as an acceptable soil erosion rate. The predicted average soil loss for the study year exceeds the tolerable soil erosion rates of 12 t ha −1 year −1 and the estimated soil formation rate of 2 t ha −1 year −1 for Ethiopia (Hurni 1983b), which can affect soil productivity.
In the presents study, the average annual rate of soil loss varies with LULC types, with cultivated land being the main contributor (about 18 t ha −1 year −1 in 1985, 19 t ha −1 year −1 in 2002, and 21.7 t ha −1 year −1 in 2017), followed by grasslands and shrub-bushlands (Table 8).
The lowest soil loss rates were observed in forest and   (Table 8). This is due to vegetation cover and the corresponding low C-factor values (Zhou et al. 2008). The highest rates of soil loss on cultivated land indicated that land conversion from natural vegetation (e.g., forest and shrub-bushland) to cultivated land would exacerbate land degradation due to soil erosion. Similarly, Fenta et al. (2021) estimated an average annual soil loss rate of 36.4 t ha -1 year -1 from cultivated land for Ethiopia's highlands. This is twice the overall average rate of soil loss (16.5 t ha -1 year -1 ), they reported for the whole highlands of Ethiopia. Contrary to our research results, Taye et al. (2018) reported the highest soil loss rate of 38.7 t ha -1 year -1 from grassland compared with 7.2 t ha -1 year -1 from cropland.
Studies showed that rainfall is one of the most sensitive factors for soil erosion (Zhang and Nearing 2005;Her et al. 2019;Berberoglu et al. 2020;Borrelli et al. 2020;Ciampalini et al. 2020;Eekhout and De Vente 2020). The present study results show an increase in the future period's annual rainfall erosivity compared to the baseline. The RCP8.5 scenario is expected to have the highest amount of erosivity factor and will have the highest erosion rate, followed by RCP4.5 in 2018-2050. The rate of soil erosion under RCP4.5 and RCP8.5 is expected to increase by 9.7% and 12.7%, respectively, compared to the baseline period.
The spatial distribution of soil erosion in the Muga watershed is shown in Fig. 8. As shown in Table 8 The results of this study are consistent with the results in the Upper Blue Nile basin of Ethiopia (e.g., Bewket and Teferi 2009;Mengistu et al. 2015). However, the rate of soil erosion in the study area is lower than the value previously reported by Yesuph and Dagnaw (2019) for Beshillo catchment of the Blue Nile Basin, Ethiopia, which is 37 t ha −1 year −1 . According to Kouli et al. (2009), soil loss rates above 10 t ha −1 year −1 will not reverse for 50 to 100 years. Considering this threshold, the total area of the study watershed, where the risk of soil erosion exceeds the soil loss tolerance is expected to be 19,348 ha and 21,057 ha in the 2050s under RCP4.5 and RCP8.5, respectively.
The future scenario shows that LULC combined with climate change substantially increases the average soil erosion by 2070s globally (Borrelli et al. 2020). The results of this study also showed that when projected land use is combined with simulated climate change, the mean annual soil erosion rate is expected to increase by 13.2% under LULC map of 2033 and RCP4.5 and 15.7% under LULC map of 2033 and RCP8.5 compared with the baseline.
The results of this study showed that the trend of soil erosion increased during the years of study and is expected to continue in the future due to LULC and climate change. As a result, it affects agricultural productivity and hydrological process in the study watershed. Thus, there should be appropriate land management strategies that take into account the future LULC and climate change to manage these ecological processes, because the LULC and climate change have a profound impact on the ecosystem of the Muga basin.

Conclusion
The study demostrated the impacts of LULC and climate changes on soil erosion using an integrated approach of CA-Markov chain, climate and soil erosion models. The outcomes indicate that soil erosion rate in the Muga watershed shows an increasing trend from 19.7 t ha −1 year −1 in 2017 to 20.7 t ha −1 year −1 in 2033 due to LULC change. Furthermore, by the 2050s, the rainfall erosivity factor may increase, which can lead to more soil erosion rate. Hence, the soil loss rate in Muga watershed is projected to increase to 22.0 t ha −1 year −1 and 22.8 t ha −1 year −1 under RCP4.5 and RCP8.5 scenarios, respectively, due to higher erosive power of the future intense rainfall. When the combined effect of LULC and climate change considered, the average annual soil loss rate was increased by 13.2% and 15.7% under RCP4.5 and RCP8.5, respectively, which is much higher than the individual effects of LULC and climate change. The change in soil erosion rate showed a spatial variation, hence the upper escarpment of the watershed which consists of steep slopes is highly vulnerable to soil erosion hazard.
The outcome of this investigation would be imperative in the decision and implementation of appropriate soil and water conservation methods and local climate adaption strategies within the study area and for deducing the changes in the future. Likewise, the results obtained from this study can be used by local and regional government agencies, developers and policymakers to diminish the rate of soil loss in the study watershed. Integrated use of CA-Markov chain, climate, and soil erosion models have demonstrated to provide relevant information about the impacts of LULC and climate change on soil erosion at a watershed level.