Study area description

The study was conducted at the Agewmariam watershed situated in the Wag Himra Administration zone in Amhara Regional State, Northern Ethiopia. The study area is located at 38° 55’ 10” to 38° 56’ 10” E longitudes to 12° 31’ 40” to 12° 32’30” N latitudes (Figure 1). The elevation of the watershed ranges from 2075 to 2393 m above the sea level. The study area is characterized by a unimodal rainfall pattern; the main rain extends from late June to early September. The mean annual rainfall of the area varied from 689.3 to 1087.9 mm, and the mean minimum and maximum annual temperatures are 12.8°C and 28°C, respectively. The climate of the study area is dry semi-arid lowland.³² The dominant soil types are eutric cambisols (64.8%) and eutric regosols (35.2%).³³ The landholding size in the watershed is characterized by a small and fragmented plot with less than 0.75 ha per household.³⁴ The major crops grown are sorghum (Sorghum bicolor L.), wheat (Triticum aestivum L.), barley (Hordeum vulgare L.), teff (Eragrostis tef [Zucc.] Trotter), chickpea (Cicer arietinum L.), and faba bean (Vicia faba L.). The dominant vegetation species are Acacia etbaica, Acacia tortilis, Acacia seyal, Albizia amara, Olea europaea, and Cordia africana. According to the Bureau of Agriculture, the numbers of heads of the households and the total population in the study area were 259 and 1113, respectively.

Figure 1.

Location map of the study area.

Data sources and methods

Soil loss modeling

The USLE²⁸ was used for soil loss modeling. The soil loss was calculated based on grid cell analyses by multiplying the respective USLE factor values (R, K, LS, C, and P) using the following equation

(1)

where A is the annual soil loss (t ha⁻¹ year⁻¹), R is the rainfall erosivity factor (MJ mm h⁻¹ ha⁻¹ year⁻¹), K is the soil erodibility factor (Mg ha⁻¹ MJ⁻¹ mm⁻¹), LS is the slope length factor (dimensionless), C is the management factor (dimensionless), and P is the conservation practice factor (dimensionless).

Data sources and analysis

Erosivity (R) factors

The R factor was calculated using the equation given by Hurni³⁵ adapted for the Ethiopian condition which has been derived from spatial regression analysis³⁶ using equation (2).

The R factor is given by a regression equation as Hurni³⁵

(2)

where R is the erosivity factor and P is the mean annual rainfall (mm year⁻¹). The source of precipitation data for the current and future climate is described in climate modeling below.

Soil erodibility (K) factor

Soil samples were collected from the study watershed by a grid of 100 × 100 m, and the sample for analysis was taken from the center of each grid. The 155 soil samples were collected to estimate the physical and chemical properties of the soil. The analyzed soil properties are organic matter content, soil texture, soil structure, and soil permeability. The soil sample analyses were performed based on standard laboratory procedures. Particle size distribution was analyzed using the hydrometer method,³⁷ whereas organic matter was analyzed by the wet combustion method of Walkley and Black as outlined by Nelson and Sommers.³⁸

Soil structure was identified under field conditions with the help of a soil structure assessment kit to determine soil structural class code, shape, and size as adopted from the USLE nomograph.²⁸ The permeability class code was obtained from soil textural classes,³⁹ which are encoded from the textural triangle based on the observed soil texture.

Finally, the K factor was calculated from estimated soil properties, namely texture, organic matter, structure, and permeability.^40,41 Using the following equation (equation (3))

(3)

(4)

where C is % of clay (<0.002 mm), L is % of silt (0.002-0.05 mm), Armf is % of very fine sand (0.05-0.1 mm), OM is the organic matter content (%), p is a code indicating the class of permeability, and s is a code for structure size.

Soil sample points of erodibility factors were converted to surface data by interpolation techniques using ordinary kriging in the ArcGIS10.3.1 environment and the Gaussian model.⁴²

Slope length and gradient (LS) factor

Slope length and gradient factors were analyzed using a Shuttle Radar Topography Mission digital elevation model with a resolution of 30 m × 30 m by ArcGIS 10.3.1. which is available by the United States Geological Survey (USGS) Earth explorer. The flow accumulation, slope steepness, and slope gradient were estimated. The LS map was generated using equation (5) developed by Wischmeier and Smith²⁸ following similar approaches conducted by other researchers^22,25,43

(5)

where LS is the slope length and steepness factor and X is the slope length (m). The cell value is 30 m, m is a variable slope length exponent, and S is a slope gradient (%).

Climate change modeling

Both the baseline and future climate data were used in this study. The baseline daily rainfall data were acquired from Koziba weather station (38° 55’ 19” E, 12° 29’ 3” N) and the Ethiopia Meteorological Agency for the period 1989 to 2018. The projected daily and monthly precipitation data were downloaded from CMIP5 of the Royal Netherlands Meteorological Institute (KNMI) Climate Explorer between 2011 and 2100 under 2 Representative Concentration Pathways, RCP 4.5 (medium stabilization scenario) and RCP 8.5 (high baseline emission scenario). A detailed description of RCP can be found in Van Vuuren et al.⁴⁷ The data of 20 models (ACCESS1-0, bcc-csm1-1, BNU-ESM, CanESM2, CCSM4, CESM1-BGC, CSIRO-Mk3-6-0, GFDL-ESM2G, GFDL-ESM2M, HadGEM2-CC, HadGEM2-ES, inmcm4, IPSL-CM5A-LR, IPSL-CM5A-MR, MIROC5, MIROC-ESM, MPI-ESM-LR, MPI-ESM-MR, MRI-CGCM3, and NorESM1-M) from the CMIP5 ensemble and embedded in “Monthly CMIP5 scenario runs” were provided by KNMI Climate Explorer’s online File Depot. The RCP 4.5 and 8.5 scenarios were generated under the time scales of the near future (2011-2040), middle future (2041-2070), and far future (2071-2100). Hence, the future climate scenarios were simulated based on the observed climate scenarios.

A statistical bias correction method

Climate model output data are regularly given biased compared with observed time series, and they are making correction procedures necessary.⁴⁸ Bias correction methods are correct for the current climate conditions and also valid for future conditions. There are several methods of statistical bias correction analyzing ways.⁴⁹ Among the statically bias correction method the delta change approach was used in this study area. The future period of the study area was in 2011 to 2040, 2041 to 2070, and 2071 to 2100 simulated to the baseline period from 1989 to 2018. Future climate time series were constructed using the delta change approach. This involves observed climate time series by mean changes (differences or ratios of changes) simulated with global climate models (GCMs). The changes were determined as monthly precipitation changes (in %) values from the base period during 1989 to 2018. The formula of the delta change bias correction method is as follows⁴⁹

(6)

(7)

P is the precipitation in mm, µ is the mean, M is a month, cont is the simulate baseline period, obs is the observed and scen is the simulated scenario future period RCP 4.5 and 8.5.

Statistical techniques and model performance assessment

Descriptive statistics such as standard deviation, coefficient variance, and relative error were analyzed using IBM SPSS software version 20. The performance of model fit related to precipitation projection and observed data were assessed based on the coefficient of determination (R²), root mean square error (RMSE), and Nash-Sutcliffe efficiency (NSE). Statistical bias correction was employed to prevent overfitting and underfitting of the data. The future period of 2020s, 2050s, and 2080s was simulated from the CMIP5 scenarios of RCP 4.5 and RCP 8.5 to produce the future estimate of the erosivity factors. Finally, estimates of the future period’s soil loss were derived using predicted erosivity factors and the other erosion factors which were affected indirectly by climate change in the study area. Figure 2 shows the detailed process of the methodology.

Figure 2.

Flowchart of current and future climate change impact on soil erosion.

DEM indicates digital elevation model; GCM, global climate model; LULC, land use land cover; RCM, regional climate model; RCP, Representative Concentration Pathways.

The current rate of soil loss by water

The annual precipitation of the study area was 890.23 mm. The result showed that the average erosivity (R) factor value in the study area was 499.2 MJ mm ha⁻¹ h⁻¹ year⁻¹.

Soil erodibility (K) factor in the central and eastern parts of the watershed ranged from 0.079 to 0.129 Mg h MJ⁻¹ mm⁻¹ and in northern and southern parts of the watershed, it ranged from 0.13 to 0.173 Mg h MJ⁻¹ mm⁻¹, which is shown in Figure 3B. The K value of the study area ranged from 0.079 to 0.173. The current result also agrees with similar findings reports in the K values of tropical soils that range from 0.06 to 0.48,⁵⁰ and most Ethiopian soils have the K values ranged from 0.05 to 0.6.⁵¹ In contrary to our result, values ranging from 0.018 to 0.033 from Gumra-Maksegnit watershed in the northwestern Amhara⁴² and values ranging from 0.0008 to 0.0073 in western Iran⁵² were reported. These lower values may be recorded because of high organic matter and surface aggregates. In our study, the mean organic matter was 1.3% in a range from 0.86% to 2.59%, which low organic matter.⁵³ The high K value indicates high erodible and vulnerability to soil erosion which could be attributed to low clay and organic matter contents as less aggregation of soil colloids.⁵⁴

Figure 3.

Model parameter map to estimate soil loss in the study area.

The slope length and gradient (LS) factor values ranged from 0 to 311 as illustrated in Figure 3A. The average topographic factor value of the study watershed was 16.33 with a standard division of 26.7. The result showed that the northern and southern parts of the study area, the hillsides, and along the gullies had higher LS values. It represents higher susceptibility to erosion; perhaps there will be a greater accumulation of runoff and high velocity.^{22,23,25,28,55}

The crop management (C) factor for different land uses was derived from satellite images based on land use and land cover maps and its attribute data analysis. The average C value in the Agewmariam watershed was 0.053. The cultivated land had a maximum cover factor which could indicate higher erosion. It covers the largest part of the study area indicated in Table 1 and Figure 3C.

Table 1.

Land use, area coverage, and cover management factor for the study area.

Erosion control practice (P) factor values were assigned 0.8 and 1, regardless of their slope for shrub and forest. However, P value for agricultural land was given accounting for its slope. Hence, the agricultural land is additionally subdivided into 6 classes based on the slope percentage. Different P values were assigned for each slope class (0%-5%, 5%-10%, 10%-20%, 20%-30%, and >50%) as illustrated in Table 2 and Figure 3D. High P values are obtained from agricultural land on slope classes greater than 30%.

Table 2.

Conservation practice (P) factor values.

The average annual soil loss was 25 t ha⁻¹ year⁻¹, but the soil loss varies from 0 t ha⁻¹ year⁻¹ in the plain area to 897 t ha⁻¹ year⁻¹ in the hilly landscapes of the study area. The maximum soil loss happens on the hilly terrains and the mainstream because of the high length and steepness factor value (50-311) and greater than the 30% slope value.⁵⁷ The erosion risk was categorized into 5 classes as shown in Figure 4. Accordingly, 67.2% of the total area has a slight erosion rate of 0 to 11 t ha⁻¹ year⁻¹ and characterized by low-risk areas. The remaining areas (32.8%) were categorized as follows: moderate erosion risk area 5.4% (11-18 t ha⁻¹ year⁻¹), high erosion risk area 5.9% (18-30 t ha⁻¹ year⁻¹), very high erosion risk area 3.3% (30-50 t ha⁻¹ year⁻¹), and severely affected erosion risk area 18.4% (50-897 t ha⁻¹ year⁻¹).

Figure 4.

The current (1989-2018) annual soil loss of Agewmariam watershed.

Statistical bias correction of global climate data

The GCM output data underestimate in the summer season (July and August) and overestimate in the winter season. The corresponding model RCP and simulated RCP scenarios from observed mean monthly precipitation were analyzed as indicated in Table 3. The coefficient of determination (R²), NSE, RMSE, and percent bias (PBAIS) of the precipitation data are shown in Table 3. The result shows that close up an agreement to the observed value illustrated in Table 3.

Table 3.

Statistical measures of monthly precipitation model output simulated and observed data sets.

The efficiency of predicted precipitation from the climate model was evaluated based on the different model efficiency criteria. The coefficient of determination (R²) value for raw data of the model and simulated or corrected data of model comparison to observed precipitation data were 0.82 and 0.99, respectively. The result confirms that there was a good relationship performance between the observed and simulated precipitation. The result agrees with the Werri watershed study conducted by Gebremeskel⁵⁸ in the Tekeze river basin with an efficiency of 0.92. The statistical tests clearly describe that the raw model has a bias that indicted an overall overestimation of the observed mean monthly precipitation. The bias-corrected value indicated a better fit and shows a positive and small PBIAS value of 0.05, RSME of 1.96, and a strong NSE value of 0.96 as indicated in Table 3. The result agrees with the national-level study conducted by Berhanu et al⁵⁹ in whole Ethiopia with the efficiency of RSME (0.87-8.37), NSE (0.95-1), and PBAIS (−1.88 to 2.5).

Future climate scenarios of precipitation and rainfall erosivity

The average annual precipitation was increased by 0.1%, 1.6%, and 1.8% for RCP 4.5 and 2.9%, 5.9%, and 11.8% for RCP 8.5 during the 2020s, 2050s, and 2080s periods, respectively, as indicated in Table 4. Our result concurs with the previous studies conducted at the major Ethiopian river basin including Tekeze,⁶⁰ current and future assessments of soil erosion on the Tibetan Plateau in Asia,¹⁴ and Keleta watershed in the Upper Awash basin of Ethiopia⁶¹ using current and future period under both scenarios (Figure 5).

Table 4.

Change in average annual rainfall erosivity (R) for future period from the base period for RCP 4.5 and RCP 8.5 scenarios.

Figure 5.

Change in average monthly precipitation for the different future period and scenarios: (A) RCP 4.5 and (B) RCP 8.5 concerning the base period.

RCP indicates Representative Concentration Pathways.

Average annual rainfall erosivity was calculated by using equation (2) for the current period (1989-2018) as well as for the future period (2011-2100) under RCP 4.5 and 8.5 scenarios. Similarly, the average annual erosivity was also calculated for the 3 successive periods, 2011 to 2040, 2041 to 2070, and 2071 to 2100 under both scenarios (Figure 6).

Figure 6.

Temporal variation of precipitation and rainfall erosivity in the study area.

RCP indicates Representative Concentration Pathways.

The erosivity value was found to be 492.6, 499.9, and 500.8 MJ mm ha⁻¹ h⁻¹ year⁻¹ for the RCP 4.5 scenario and 506.7, 521.6, and 550.2 MJ mm ha⁻¹ h⁻¹ year⁻¹ for the RCP 8.5 scenario, respectively (Table 4). The rainfall erosivity showed an increasing trend by 0.1%, 1.6%, and 1.8% for the 2020s, 2050s, and 2080s, respectively, under the RCP 4.5 scenario and 2.9%, 5.9%, and 11.8% under the RCP 8.5 scenario as illustrated in Table 4.

Projected future soil loss of the study area

The result showed that under the RCP 4.5 scenario, the average annual soil loss increased by 0.76 t ha⁻¹ year⁻¹ (3.0%), 1.18 t ha⁻¹ year⁻¹ (4.7%), and 1.29 t ha⁻¹ year⁻¹ (5.2%) in the period of 2020s, 2050s, and 2080s, respectively, as illustrated in Figure 7B to D and Table 5. Likewise, under the RCP 8.5 scenario, the average annual soil loss increased by 1.5 t ha⁻¹ year⁻¹ (6.0%), 2.38 t ha⁻¹ year⁻¹ (9.5%), and 3.58 t ha⁻¹ year⁻¹ (14.3%) during the period of 2020s, 2050s, and 2080s, respectively, as shown in Figure 8B to D and Table 5 compared with the base period (1989-2018) as indicated in Figures 7A and 8A.

Figure 7.

Soil erosion under RCP 4.5 for the (B) 2020s, (C) 2050s, and (D) 2080s periods, respectively, concerning the (A) base period.

RCP indicates Representative Concentration Pathways.

Table 5.

Change in soil erosion scenario under RCP 4.5 and 8.5 for the 2020s, 2050s, and 2080s periods, respectively, from the base period.

Figure 8.

Soil erosion under RCP 8.5 for the (B) 2020s, (C) 2050s, and (D) 2080s periods, respectively, concerning the (A) base period.

RCP indicates Representative Concentration Pathways.

The average annual soil loss rate estimated for the entire watershed was 25 t ha⁻¹ year⁻¹ which is comparable to other reports by Hurni³⁵ for the highland of Ethiopia (20 t ha⁻¹ year⁻¹); by Ayalew and Selassie⁶² for the Guang watershed in the Northwestern Ethiopia (24.95 t ha⁻¹ year⁻¹); by Gashaw et al²⁵ for the Geleda watershed in Northwestern Ethiopia (23.7 t ha⁻¹ year⁻¹); and by Amsale and Mengaw²³ for the JabiTehinan watershed in the Northwestern Highlands (30.4 t ha⁻¹ year⁻¹). In contrast to our result, some studies reported a high erosion rate in different parts of the Ethiopian highlands. For instance, Bewket and Teferi²¹ for the Chemoga watershed of the Blue Nile basin in the Northwestern Highlands of Ethiopia (93 t ha⁻¹ year⁻¹) and Gelagay and Minale²⁴ for the Koga watershed in the Blue Nile basin (47.4 t ha⁻¹ year⁻¹). The relatively low average soil erosion rate in our study area could be due to topography with an average slope (27%), and the cultivated areas are located in the gentle slopes of the catchment.

Contrary to our finding, other reported very low average soil erosion rate, for instances, Medago watershed in the Northern Highlands with a rate of 9.63 t ha⁻¹ year⁻¹ by Brhane and Mekonen,⁵⁷ 9.1 t ha⁻¹ year⁻¹ by Gizachew,⁶³ 7.47 t ha⁻¹ year⁻¹ by Pongsai et al,⁴¹ 4.17 t ha⁻¹ year⁻¹ by Eshete,⁶⁴ and 4.81 t ha⁻¹ year⁻¹ by Tiruneh and Ayalew⁷ in the highlands of Ethiopia.

These results may be due to the highly variable topographic nature of the Ethiopian highlands, records of the gentle slope feature, and the land use system. In our study area, a high erosion rate was recorded in the steeper slope area of the watershed ranging from 30% to 83% in the steeper slope of cultivated lands which is similar to other studies.^25,57 The results showed a higher mean soil loss rate under RCP 8.5 than RCP 4.5 as illustrated in Table 5, Figures 4 and 5. Our result also confirmed with other reports by Wang et al⁶⁵ in the Great Lakes region (4.99%-23.2% increased averagely both scenarios), and Simonneaux et al¹³ soil erosion rate increased by 7.9% during the 2080s period in the Rheraya watershed in Morocco.

The rainfall erosivity and soil loss rate were augmented in the period 2020s, 2050s, and 2080s compared with the current period due to high-intensity rainfalls in the study watershed which were also confirmed with other findings.^66–71

This study results showed the linear relationship between the rainfall erosivity and soil loss under both scenarios and 3 successive future periods. The current findings agree with Lanckriet et al⁷² in the May Zeg-zeg catchment in Ethiopia. In contrast to our findings, some studies however reported a decreasing soil loss from the base period. For instance, Zhang et al⁶⁸ revealed that the soil loss decreased from the base period by −4% and −6% for the future period of 2016 to 2045 and 2026 to 2035, respectively. This contradictory result may be due to increased rainfall and intensity in our study area. The precipitation and storm intensity changes can be expected to have a greater impact on the soil loss rate.¹⁰ Accordingly, our results revealed that the future predicted soil loss rate is to increase significantly from the base period for both scenarios due to the higher erosive power of rainfall. The results indicated that the risk of soil loss could be reduced by improving the C factor and the degraded land could be recovered by improving the P factor. Thus, future soil loss can be mitigated through promoting sustainable land use policy, preventing expansion of other land use types into agricultural land use, and taking soil conservation practices.