Estimating the Best Exponent and the Best Combination of the Exponent and Topographic Factor of the Modified Universal Soil Loss Equation under the Hydro-Climatic Conditions of Ethiopia

Abstract

The effect of the topographic factor of the Modified Universal Soil Equation (MUSLE) on soil erosion and sediment yield is not clear. Except for the coefficient, soil erodibility, cover, and conservation practice factors of the MUSLE, an individual effect of the exponents and topographic factors of the MUSLE on soil erosion and sediment yield can be seen by applying the model at different watersheds. A primary objective of this paper is to estimate the best exponents and topographic factors of the MUSLE under the hydro-climatic conditions of Ethiopia. For the sake of the calibration procedure, the main factors of the MUSLE that directly affect the soil erosion process, such as cover, conservation practice, soil erodibility, and topographic factors, are estimated based on past experiences from the literature and comparative approaches, whereas the parameters that do not directly affect the erosion process or that have no direct physical meaning (i.e., coefficient a and exponent b) are estimated through calibration. We verified that the best exponent of the MUSLE is 1 irrespective of the topographic factor, which results in the maximum performance of the MUSLE (i.e., approximately 100%). The best exponent that corresponds to the best equation of the topographic factor is 0.57; in this case, the performance of the model is greater than or equal to 80% for all watersheds under our consideration. We expect the same for other watersheds of Ethiopia, while for other exponents and topographic factors, the performance of the model decreases. Therefore, for the conditions of Ethiopia, the original exponent of the MUSLE is changed from 0.56 to 0.57, and the best equations of the topographic factor are provided in this paper.

1. Introduction

Williams(1975) developed the MUSLE using 778 storm-runoff events collected from 18 small watersheds [1,2], with areas varying from 15 to 1500 ha, slopes from 0.9 to 5.9%, and slope lengths of 78.64 to 173.74 m (Hann et al.1994) as cited in [3]. The MUSLE is given by

$$y=a{\left(Qq\right)}^{b}KLSCP$$

where y is the sediment yield in metric tons, a is the coefficient and b is the exponent ($a=11.8$ and $b=0.56$ for USA, where the MUSLE was originally developed), Q is the runoff volume in $m^3$, q is the peak runoff rate in ${\mathrm{m}}^3{\mathrm{s}}^{-1}$, K is the soil erodibility factor in $0.01\ast tons\ast acre\ast hour\ast acr{e}^{-1}\ast yea{r}^{-1}\ast foo{t}^{-1}\ast ton{s}^{-1}\ast inc{h}^{-1}$, L is the slope length factor, S is the slope steepness factor, C is the cover factor, P is the soil conservation practice factor. Essentially, the MUSLE was developed for a small agricultural watershed, where the extent of erosion is from sheet to rill erosion.

However, we cannot exactly tell whether it considers gully erosion or not. To apply the MUSLE for a large watershed (in our case, the Hombole Watershed is 762,281 ha, Mojo Watershed is 150,282 ha, Gumera Watershed is 127,805 ha, and Gilgel Gibe 1 Watershed is 292,809 ha), one approach that was proposed is using the MUSLE in Soil and Water Assessment Tool (SWAT) environment.

This may be because the sediment yield can be more accurately estimated if the large watershed is divided into subwatersheds (area < 2590 ha) to compensate for nonuniformly distributed sediment sources; the effect of watershed hydraulics and sediment particle size can be included by routing the sediment yield from subwatersheds to the large watershed [1]. As part of the evaluation of the model, we considered the specific behavior of the MUSLE, experiences of other researchers, the physical connection between factors of the MUSLE, and the suitability of the model toward a specific location.

If we consider the specific behavior of the MUSLE, we found that the MUSLE showed better performance in the case of directly measured flow data rather than indirectly obtained flow data using indirect methods [3]. The model also provides appropriate estimates at a watershed rather than an experimental plot as was reviewed and reported by [3]. In this connection, if we consider SWAT, SWAT uses an indirect method (such as the Soil Conservation Service Curve Number) to generate runoff, then it uses the MUSLE to estimate the soil loss from a hydrologic response unit (which is similar to a plot scale), then SWAT routes sediment output in channels to the outlet of a large watershed. However, this also leads to accumulative error at the end due to uncertainty in the definition of a channel, channel depth, and width in the SWAT environment.

If we consider the experiences of some other researchers, the MUSLE is unsuitable for the prediction of the sediment yield for small storms [4]. However, the slight variation in hydrological response of a watershed in terms of the sediment yield might be changing in the antecedent hydrological conditions, the spatial and temporal distribution of rainfall, availability of eroded sediment throughout the watershed, which is not taken into account by the MUSLE as for many other lumped models [4].

If we consider the physical connection between factors of the MUSLE, as far as the runoff energy for soil detachment and sediment transport; the physical connection between the runoff, soil erodibility, topographic, cover, and soil conservation practice factors is convincing; however, further refining the physical connection between the factors may become necessary. For instance, the cover and soil conservation factors play a role to break runoff energy so as to protect soil loss due to runoff.

As the slope length becomes increasingly larger, there is a possibility that erosion from the upper part of the slope becomes deposited at the lower part of the slope (for instance, if we consider the last runoff from the slope-field after the end of rainfall). This is because, depending on the magnitude of the runoff and its sediment transport capacity, the runoff takes up more soil particles and becomes concentrated on its way to the bottom of the slope. In other words, the energy of the runoff decreases as resistance against flow increases along the length of the slope, and its shear force decreases.

If we consider the suitability of the model toward a specific location, the MUSLE has been observed to give good results in various applications in some parts of tropical Africa (Ndomba, 2007) as cited in [5], and it has been successfully demonstrated in sub-Saharan Africa [5]. As per the experimental plot result of sheet erosion at Enerta study site in Ethiopia, the MUSLE was better at estimating soil loss from a cultivated field than the USLE [6].

Therefore, based on the above limitations and advantages, the MUSLE should be tested at a watershed scale rather than a plot scale (which is similar to hru) under the hydro-climatic conditions of Ethiopia, using directly measured flow data. In this regard, we have daily average flow records but we do not have daily peak flow records in our database. As we explained above, the effect of the slope length factor on soil erosion and sediment transport should also be investigated. Therefore, based on these problems, we consider the following cases for the evaluation of the model.

If we consider the simulation time step, the daily sediment yield may not reflect daily watershed information, such as land cover, soil erodibility, and conservation activities. The reason for this can be soil erosion, sediment transport, deposition, consolidation, and re-suspension are quite complex processes, which depend on physical, biological, mechanical, and chemical activities within a large heterogeneous watershed.

Due to these complex processes, the soil that was eroded at an unknown last time can be transported, deposited, consolidated, re-suspended, and can reach outlets at different times. Therefore, the measured sediment at the outlet at the current time may not reflect the current information about the watershed; it rather reflect the unknown last time. This may be because sediment that was deposited along the length and the bottom of the slope by small runoff energy at a previous time can be transported by high runoff energy at the current time.

In the original development of the USLE, the annual soil erodibility factor was taken to compute the annual soil loss from the unit plot. Based on a previous [7] formulation, we can conclude that the annual soil erodibility is the average of soil erodibility ranging from loose to compacted soil due to rainfall impact. As the soil erodibility factor of the USLE and MUSLE is the same, the annual time step is preferred over the daily time step (in the case of SWAT).

The annual simulation time step enables the consideration of gully erosion (gully erosion is usually estimated on an annual basis [8]; it is important to note here that gully erosion is a common problem in Ethiopia(e.g., [9,10,11,12,13])); to consider the gradual soil erosion process and gradual changing activities, such as the cyclic behavior of agricultural activities, conservation practice, flood protection activities, plant growth, and harvest with respect to the rainfall pattern and extreme events in a one-year full cycle.

If we consider the hydrologic response unit (hru) in the SWAT environment, as the number of hrus becomes increasingly larger, we should consider the spatial variability of land use, soil, and slope over the entire watershed. To test the MUSLE at a watershed scale, sediment and flow routing in stream channels of SWAT are not considered (it is important to note here that there is uncertainty in the definition of a channel, channel width, and depth in the SWAT environment). Therefore, we only considered hrus to calculate the areal weighted average to capture the spatial variation of soil, cover, conservation practice, and topography.

If we consider the calibration parameters, all parameters ($a,b,Q,q,K,L,S,C,P$) of the MUSLE can potentially be used for calibration and validation [14]. Researchers [15] conducted global sensitivity analysis (Monte Carlo sampling) of the parameters of the MUSLE using the extended Fourier amplitude sensitivity test (EFAST) method. Accordingly, the exponent b is the most sensitive parameter to predict the amount of soil loss, followed by P, a, LS, C and q, and k’s influencing variables, such as organic matter, soil structure class, and soil permeability class.

In addition, researchers [16] used Sobol’s sensitivity analysis and found that the coefficient a and the exponent b were the most sensitive parameters of the MUSLE model contributing about 66% of the variability in the output sediment yield, at upper Malewa catchment in Kenya. On a storm event basis, researchers [5] estimated the location parameters ($a=12.4$ and $b=0.51$) of the MUSLE for Ofuloko watershed in Nigeria.

In some studies, only the exponent of the model was calibrated, which is logically more acceptable as was reviewed and reported by [3]. The calibrated sediment does not reflect the actual soil erodibility and conservation practice factors on the ground unless they are otherwise measured.

To accept our calibration, we should also check the calibrated value of the soil erodibility and conservation practice factors against the actual ones on the ground. This is because their product effect is reflected in the MUSLE rather than their individual effect during the calibration of the sediment yield. Unless otherwise, we cannot reach a certain conclusion that these factors are truly affecting the soil erosion process. For a given uniform watershed, the temporal variation of the soil erodibility, cover, and conservation practice factors is expected.

As the temporal variation of these factors is difficult to measure in a large watershed, we may estimate them through calibration. However, it is highly preferable if these factors are measured and studied at a temporal and spatial scale to understand their effect on soil erosion in a particular field. Any change in these factors affects the coefficient of the MUSLE; this is because only a product effect of the coefficient and these factors is reflected in the MUSLE rather than their individual effect during the calibration of the sediment yield.

As compared to the other parameters of the MUSLE, the individual effect of the exponent of the MUSLE is reflected during the calibration of the sediment yield. Therefore, estimating the exponent of the MUSLE through calibration is more feasible than other parameters of the MUSLE. For a given uniform watershed, the topographic factor does not change with time (i.e., it has a constant effect), and the effect of the topographic factor can be seen when the MUSLE is applied at different watersheds.

From this explanation, the independent effect of the exponents and topographic factors of the MUSLE can be seen by applying the model at different watersheds. In general, runoff and sediment data reflect the hydro-climatic conditions of a particular area, which independently affect the overall calibration process. Therefore, our main objective is to estimate the best exponent of the MUSLE and to estimate the best combination of the exponents and topographic factors of the MUSLE by applying the model at different watersheds of Ethiopia.

For the sake of the calibration procedure, the main factors of the MUSLE that directly affect the soil erosion process, such as cover, conservation practice, soil erodibility, and topographic factors, are estimated based on the past experiences from the literature and comparative approaches, whereas the parameters that do not directly affect the soil erosion process or that have no direct physical meaning (i.e., coefficient a and exponent b) are estimated through calibration.

2. Materials and Methods

2.1. Description of Study Areas

To begin our work, we considered four watersheds, the Gumera watershed in the Abbay River Basin, Gilgel Gibe 1 watershed (at Assendabo) in the Omo-Gibe River Basin, and Hombole and Mojo watersheds in the Upper Awash River Basin, in Ethiopia. We describe the hydro-climate, land use, and soil of the study areas based on the data that were prepared or obtained from different sources. Records of the climate, sediment and flow data of each watershed are given in the

Table A1. Data types and record period.

Watershed	Data Type	Record Period
Hombole and Mojo	Climate data	1986–2020
	Flow data	1990–2016
	Sediment data	1989–2015
Gumera	Climate data	1986–2020
	Flow data	2000-2017
	Sediment data	1990–2017
Gilgel Gibe 1	Climate data	1986–2020
	Flow data	2000–2015
	Sediment data	1990–2017

. The reference coordinate system EPSG:4326-WGS 84 is used to describe the geographic location of each watershed.

2.1.1. Upper Awash River Basin

Upper Awash River Basin drains into the Koka hydroelectric power reservoir. The basin comprises two main gauged watersheds: Hombole and Mojo watersheds, which cover 65.26% and 12.87% of the total area of the basin, respectively, and the basin also includes an ungauged watershed, which covers 21.87% of the total area of the basin. The total drainage area of the basin is estimated to be 11,680.25 km${}^2$.

In the basin, there are active socio-economic activities, such as agricultural, industrial, and commercial activities. On the other angle, the basin experienced catastrophic flooding, and land degradation problems due to severe gully erosion. The gully erosion assessment in the basin was reported by [11].

For the Hombole and Mojo watersheds (see Figure 1), the monthly average rainfall and monthly average outflow discharge at the main outlet point of the watershed are given in Figure 2 and Figure 3, respectively. The annual average maximum and minimum temperatures are 25.56 and 10.06 °C, respectively. The dominant soil types of the watersheds are given in the Figure 6. Land-use changes have been observed in the watersheds as shown in the Figures 9–12.

2.1.2. Gumera Watershed

Gumera watershed drains into Lake Tana. Its geographic boundary is given in the Figure 4. The total drainage area of the watershed is estimated to be 1278.05 km${}^2$.

The monthly average rainfall, and monthly average outflow discharge at the main outlet point of the watershed are given in Figure 2 and Figure 3, respectively. The annual average maximum and minimum temperatures are 25.38 and 10.02 °C, respectively. The dominant soil types of the watershed are given in Figure 7. Land-use changes have been observed in the watershed as shown in the Figures 13 and 14.

2.1.3. Gilgel Gibe 1 Watershed

Gilgel Gibe 1 watershed drains into the Gilgel Gibe 1 hydroelectric power reservoir. Its geographic boundary is given in Figure 5. The total drainage area of the watershed is estimated to be 2928.09 km${}^2$.

The monthly average rainfall, and monthly average outflow discharge at the main outlet point of the watershed are given in Figure 2 and Figure 3, respectively. The annual average maximum and minimum temperatures are 25.36 and 11.7 °C, respectively. The dominant soil types of the watershed are given in Figure 8. Land-use changes have been observed in the watershed as shown in the Figures 15 and 16.

2.2. Preparation of Soil Maps

Soil data is required to estimate the soil erodibility factor of the MUSLE. The necessities of preparing soil maps are to assign a specific type of soil from a general category of the soil and to maintain the spatial variability of soil. For all our watersheds, national soil maps of Ethiopia, which we obtained recently from the River Basin Authority of Ethiopia, show the general category of soil. To assign a specific type of soil, we locate the shapefile of each watershed on a harmonized world soil data map, and we clip the harmonized world soil data map to the size of our watersheds in the QGIS environment.

Then, we compare the national soil maps of Ethiopia, the harmonized world soil map, and the field observation report from the International Soil Reference and Information Centre on QGIS. Particularly for the Upper Awash River Basin, we have two soil maps that were prepared at different times from the River Basin Authority of Ethiopia. Based on these two soil maps, we maintain the spatial variability of soil right after the specific type of soil was assigned by locating an areal coverage of the specific type of soil on the old map that completely lies inside a large area of another specific type of soil on the current soil map. Therefore, the soil maps of each watershed, which are finally prepared, are given in Figure 6, Figure 7 and Figure 8.

2.3. Preparation of Land Use Maps

Preparation of land use data is necessary to estimate the cover and conservation practice factors of the MUSLE. Based on our assessment of land use and land cover by the support of Google Earth Pro, Planet Explorer, literature review (e.g., [17,18,19]), and land use maps from the River Basin Authority of Ethiopia, land-use change has been observed in the study areas. As the basis of classification of land use maps, dominant land use classes are categorized at 30 m spatial resolution. This is an acceptable level of spatial dimension to consider the spatial variability of land use at a tolerable level of accuracy. As a result, land use maps of each watershed are prepared based on a comparative approach and logical sequence.

To prepare land use maps by the comparative approach, sample geographic coordinate points with defined land use classes are collected from the Global land service map; forest and agricultural land on historical imagery in the Google Earth Pro at different acquisition dates are digitized. A time demarcation of the land-use change classification depends on a number of available baseline land use maps per watershed, the time boundary of the Global land service maps, and historical imagery in the Google Earth Pro. As a result, the time demarcations of land-use changes for the Hombole and Mojo watersheds are 1989–2000, 2001–2008, 2009–2012, and 2013–2015 and for the Gumera and Gilgel Gibe 1 watersheds are 1989–2009 and 2010–2015.

During the comparison of the above land use data files with the baseline national land use maps of Ethiopia on QGIS and Google Earth Pro; the vector data files are converted from the shapefile to the Keyhole Markup Language (KML) and vice versa. To prepare a land use map by the logical sequence, we check whether a change in land-use from one class to another is possible or not (for example, is the change from urban to agriculture possible?) such as the comparison of different land use data files that were prepared or acquired from different sources at the specified time demarcation.

Particularly for the Upper Awash River Basin, land use classes, found on the previous baseline map but not on the latter map, are included on the latter map based on the logical sequence, and vice versa. Furthermore, also missing land use classes, such as water bodies, are added on either of the maps during the comparison of the maps with other sources, such as historical imagery in Google Earth Pro, while following these procedures, the land use maps that are finally prepared for each watershed are given in Figure 9, Figure 10, Figure 11, Figure 12, Figure 13, Figure 14, Figure 15 and Figure 16.

2.4. Sediment Rating Curves

A sediment rating curve is required to generate sediment data from corresponding flow data. The linear regression equation and nonlinear regression equations, such as power function, the second and third-order polynomial function can be used to model the sediment rating curve (e.g., [20]). Different authors indicate that the power function is a commonly used nonlinear regression approach to model the sediment rating curve (e.g., [21,22,23]). The power function is given by:

$$C=a{Q}^b$$

where C is the suspended sediment load or concentration, Q is the discharge, a is the coefficient, and b is the exponent. Different authors reviewed physical meanings associated with the coefficient a, and the exponent b (e.g., [21,24,25]). Accordingly, the coefficient a represents an index of soil erodibility, whereas the exponent b is considered as an index of erosivity and transport capacity of a river. Thus, the power function can be derived by interpreting or deducting the Modified Universal Soil Loss Equation (MUSLE), where its topographic, soil erodibility, cover, and conservation practice factors describe a site-specific condition of a given watershed, and these factors affect the coefficient a of the power function at defined hydro-climatic conditions.

For the sake of simplicity of regression analysis, the nonlinear regression Equation (in our case, the power function) can be transformed to the simple linear regression equation by log-transform of both sides of the nonlinear equation. Accordingly

$$log\left(C\right)=log\left(a\right)+blog\left(Q\right)$$

If

$$y=log\left(C\right),d=log\left(a\right)\mathrm{and}x=log\left(Q\right)then,y=bx+d$$

The Least Squares, Reduced Major Axis Line (R.M.A.L) or other regression methods can be applied to find the best-fit regression line on logarithms of the suspended sediment load or concentration and discharge data, and back transform of the linear equation results in the power function. Despite that there are no generally accepted procedures to model the sediment rating curve, we proceed with the Least Squares regression method, which is based on the minimum sum of squared errors to estimate the coefficient b and the constant d of the best-fit linear regression equation on logarithms of suspended sediment concentration and discharge data.

$$b=\frac{{\sum }_{i=1}^n(x_i-\overline{x})(y_i-\overline{y})}{{\sum }_{i=1}^n{(x_i-\overline{x})}^2}$$

$$d=\overline{y}-b\overline{x}$$

Apart from choosing sediment load–discharge [26], logged mean loads within discharge classes [25] or sediment concentration–discharge [20] approaches, correction factors ($y=CF\ast a{Q}^b$) (e.g., [22,25]) and power function with some additive constant can be used [22,27] to improve the sediment rating curve. Furthermore, to improve the sediment rating curve, we may use the data consistency or homogeneity test in order to determine the data classes at specific hydro-climatic conditions.

While considering the above advantages and limitations to model the sediment rating curve, the relationship between discharge and the suspended sediment concentration rate is checked against land-use changes, seasonal weather variations or rainfall patterns, and period of land tillage. Accordingly, the sediment rating curve that is drawn while considering the rainfall and discharge relationship for the Gilgel Gibe 1 watershed, shows some improvement provided that one extreme discharge 319.65 ${\mathrm{m}}^3{\mathrm{s}}^{-1}$ on 23 August 2009 (no similar record in the daily average discharge from 1990 to 2015), which corresponds to the suspended sediment concentration 0.53 kgm${}^{-3}$, is removed from the records as part of the data quality check.

Furthermore, some data replication is possible to improve the sediment rating curve, due to the assumption that two measurements that are taken at very small time differences are almost the same as we only consider a pattern of record rather than a period of record, and also data record does not show watershed information. Accordingly, the sediment rating curve is drawn for the Gumera watershed, showing some improvement (the change in the coefficient of determination is from $R^2=0.324$ to $R^2=0.5091$) if it is a significant improvement.

For the Hombole and Mojo watersheds, the sediment rating curves are drawn without any pre-conditions. This is because the above pre-conditions do not work for these two watersheds. For the Mojo watershed, two inconsistent records of the rainfall (extremely large and small), flow and sediment on 7 August 1996 and 6 August 2003, are removed from the records as part of the data quality check. The sediment rating curves of all watersheds are given in Figure 17.

2.5. Estimating Factors of The MUSLE

The original factors of the USLE represent the average value to estimate the annual sediment yield. The unit plot [7] represents the worst case for the maximum soil erosion at a given rainfall event. It is practically impossible to directly measure each field slope, slope length, the temporal variation of soil erodibility, instantaneous runoff, cover change, and conservation practice for a large watershed. In the actual field, the field slope and length are not uniform, which means they are irregular.

The topographic, soil erodibility, and cover and conservation practice factors depend on the spatial resolution of the digital elevation model (DEM), soil, and land use maps, respectively. Therefore, in the actual sediment modeling, average field slope length [28] and slope steepness or simply topographic factor [14], average runoff, average soil erodibility factor [29], and average cover and conservation practice factors are taken.

2.5.1. Estimationof Runoff Factor

In the MUSLE, the runoff factor is the product of the total runoff volume and peak runoff rate. Researchers [14] reviewed that the runoff factor represents the energy used in transporting as well as in detaching sediment, which acts as the best indicator for predicting sediment yield for the individual storm event. To estimate the runoff factor, the peak runoff rate and/or volume of runoff can be obtained by direct measurement of the runoff on a storm-event basis, as well as using indirect methods, such as the Soil Conservation Service Curve Number (SCS CN) method, Rational method, flood routing, unit hydrograph, etc. In our case, we used the daily average discharge to estimate the annual total runoff volume and yearly peak runoff rate for the annual sediment yield estimation. The reasons for why we use directly measured flow data and why we estimate the annual sediment yield are addressed in Section 1.

2.5.2. Estimation of Soil Erodibility Factor (K-Factor)

Researchers [7] defined soil erodibility factor as the soil loss rate per erosion index unit for a specified soil as measured on a unit plot; the unit plot is defined as a 72.6-ft length of uniform 9% slope continuously in clean-tilled fallow; it is the continuous fallow tilled up and down the slope. The soil erodibility factor is given by [7]:

$$K=\frac{{\sum }_{n=1}^N{\left(A\right)}_n}{{\sum }_{n=1}^N{\left({EI}_{30}\right)}_n}$$

where A is the event soil loss from the unit plot in tons/acre/year, E is the storm kinetic energy in 100 × foot-tons/acre, $I_{30}$ is the maximum 30 min intensity in inch per hour, and K is the soil erodibility factor in $0.01\ast tons\ast acre\ast hour\ast acr{e}^{-1}\ast yea{r}^{-1}\ast foo{t}^{-1}\ast ton{s}^{-1}\ast inc{h}^{-1}$. It is important to note that the soil erodibility factor represents the worst or the maximum possible erosion from the unit plot with the specified field slope and length.

At the same rainfall impact pressure, less soil erosion condition that is different from the worst condition considers the soil cover and conservation practices on the same field slope and length. On the unit plot, or any unit plot for that matter, the temporal and spatial variation of the soil erodibility depend on the types of soil; the quite complex interaction of physical, biological, chemical, and mechanical processes.

From the soil erodibility table and Equations (see Figure 18), we can reveal that the soil erodibility factor varies from 0 to 1, where 0 indicates the soil that is hard to erode, whereas 1 represents easily erodible soil by the same rainfall impact pressure under otherwise similar soil erosion conditions. From this range of the soil erodibility factor, we can conclude that soil erodibility refers to the degree of being easy to erode a given soil.

The soil erodibility factor (K-factor) can be estimated by direct field measurement or by using different empirical equations or a soil erodibility nomograph.

(1)

The K-factor that was originally developed for the soil conditions of the USA [7]:

$$K=\frac{\left\{\left[2.1\ast M^{1.14}\ast \left({10}^{-4}\right)\ast \left(12-a\right)\right]+3.25\ast \left(b-2\right)+2.5\ast \left(c-3\right)\right\}}{100}$$

where K = the soil erodibility in $0.01\ast tons\ast acre\ast hour/acre\ast year\ast foot\ast tons\ast inch$; $M=(\%silt+\%veryfinesand)\ast (100-\%clay)$, $M=$ Particle-size parameter, silt (%) = percentage of silt, % very fine sand = percentage of very fine sand (0.1 to 0.05 mm), clay (%) = percentage of clay, $a=$ percentage of organic matter, $b=$ soil structure code used in soil classification, and $c=$ profile permeability class. For soils containing less than 70 percent silt and very fine sand, the nomograph [7] is used to solve the above equation.

Some comments on this equation: we do not have a percentage of very fine sand in our database to test the equation. Our source of data is the harmonized world soil data, which includes the texture, reference soil depth, drainage class, available water capacity, sand, silt and clay fraction, bulk density, gravel content, organic carbon content, pH, cation exchange capacity, base saturation, total exchangeable bases, calcium carbonate content, gypsum content, sodicity, and salinity content. As land tillage and mechanical compaction (due to rainfall impact) change the structure of the soil; the structure of tilled, bare, or compacted soil varies at temporal and spatial scales. As soil permeability depends on soil texture and organic matter, their relationship should be explicitly shown. Unrealistic values were obtained for tropical soils from the equation’s erodibility nomograph (Mulengera and Payton, 1999; Ndomba, 2007) as cited in [5].

(2)

The K-factor (Williams and Renard, 1983) as cited in [30] and similar equation in [31,32].

$$K=(0.2+0.3\ast exp(-0.0256\ast S_a\ast (1-\frac{s_i}{100})))\ast {\left(\frac{S_i}{c_L+s_i}\right)}^{0.3}\ast (1-\frac{0.25c}{c+exp(3.72-2.95c)})\ast (1-\frac{0.7S_N}{S_N+exp(-5.51+22.9S_N)})$$

where $S_a$ = sand (%); $S_i$ = silt (%); $C_L$ = clay (%); $S_N=1-(S_a/100)$; C = organic carbon.

(3)

The K-factor that was tested in the soil conditions of the Philippines [33]:

$$K=\left[0.043\ast pH+\frac{0.62}{OM}+0.0082\ast S-0.0062\ast C\right]\ast Si$$

where $pH$ = pH of the soil, $OM$ = organic matter (%), S = sand content (%), C = clay ratio = % clay/(% sand + % silt), and $Si$ = silt content = % silt/100.

(4)

The K-factor that was originally developed for the volcanic soil of Hawaii, USA (El-Swaify and Dangler, 1976) as cited in [34]:

$$K=-0.03970+0.00311\ast x_1+0.00043x_2+0.00185x_3+0.00258x_4-0.00823x_5$$

where $x_1$ = unstable aggregate size fraction (<0.250 mm)(%), $x_2$ = modified silt (0.002–0.1 mm) (%) × modified sand (0.1–2 mm) (%), $x_3$: % base saturation, $x_4$ = silt fraction (0.002–0.050 mm) (%), and $x_5$ = modified sand fraction (0.1–2 mm) (%).

We do not have unstable aggregate size fraction or modified silt and sand data in our database to test the equation.

(5)

Williams (1995) proposed the following K-factor as cited in [35]:

$$\begin{array}{c}K=f_{csand}\ast f_{cl-si}\ast f_{orgC}\ast f_{hisand}\hfill \\ f_{csand}=0.2+0.3exp[-0.256m_s(1-\frac{m_{silt}}{100})]\hfill \\ f_{cl-si}={\left(\frac{m_{silt}}{m_c-m_{silt}}\right)}^{0.3}\hfill \\ forgC=1-\frac{0.25\ast orgC}{orgC+exp[3.72-2.95orgC]}\hfill \\ f_{hisand}=1-\frac{0.7(1-\frac{m_s}{100})}{1-\frac{m_s}{100}+exp[-5.51+22.9(1-\frac{m_s}{100})]}\hfill \end{array}$$

(6)

Other soil erodibility equations are mentioned in [32,34,35,36,37,38,39,40].

To test the soil erodibility equations on the basis of the original definition of the soil erodibility by [7], the following conditions should be fulfilled. From the MUSLE,

$$K=\frac{y}{a{\left(Qq\right)}^b\ast LSCP}$$

where K represents the worst condition for the maximum erosion case when the slope-field length is 22.13 m and the slope angle is 9%. In this case, no cover and conservation practices are employed in the field to give protection against soil erosion; the land is tilled up and down the slope, and therefore the maximum erosion is expected. In the above equation, K represents the maximum erosion case when the observed sediment yield (y) is due to soil erosion from a field with a specified slope length, slope angle, cover, and conservation practice. If we take $C=P=1$, K represents the maximum erosion from the field with the specified slope length and angle.

However, our observed sediment yield does not represent the worst conditions for the maximum erosion case; we have some magnitude of the cover and soil conservation practice to give protection against soil erosion, and land is not tilled up and down the slope. Therefore, in this case, K represents the minimum value as compared to the actual value that will be obtained from the soil erodibility equation for the worst conditions for the maximum soil erosion case ($K_{eq}$)

$$K_{min}=\frac{y}{LSa{Q}^b}<<k_{eq}.$$

For our watersheds, the minimum K value is calculated by replacing the annual sediment load, runoff volume, and topographic factor (the reasons why we use the annual erodibility factor are given in Section 1). However, at this point, the coefficient and exponent of the model are not known for our watersheds; therefore, the minimum reference value ($K_{min}$) is set. Based on the soil data we have, the actual soil erodibility factor is calculated using the soil erodibility equations that were proposed by [33,35], and Williams and Renard (1983) as cited in [30]. Accordingly, the graphs of the K-factor are shown in Figure 18.

From the graph, a reasonable actual erodibility graph for our watersheds lies between the minimum K-factor graph and the calculated K-factor using Williams’s (1995) equation as cited in [35]. To proceed with Williams’s (1995) equation as cited in [35], Williams’s sub-K-factors are calculated and compared based on the silt and sand content, clay and silt content, and organic carbon content of our soil data. Comparatively speaking, the soil erodibility increases if teh silt content increases, and the sand and clay content decreases. This is because the interaction between soil particles ranges from the loose interaction for silt soil to the strong interaction for clay soil. Humus, manure, organic matter, and the organic carbon content decreases soil erodibility as it binds the soil particles together, or it provides protective cover for the soil particles.

From Figure 19, the soil erodibility factors conform to the general comparisons stated above. Therefore, we use Williams’s (1995) equation as cited in [35] to calculate the soil erodibility factor using soil data of the watersheds under our consideration and watersheds of Ethiopia in general.

2.5.3. Estimation of the Slope Steepness and Slope Length Factors

The slope steepness factor (S) is the ratio of soil loss from a field slope gradient to soil loss from the 9% slope under otherwise identical conditions [41]. A high rate of soil loss is associated with steep slopes [29,42], and soil-loss prediction is more sensitive to the slope steepness than slope length [43].

Slope length is defined as the distance from the point of origin of overland flow to the point where either the slope gradient decreases enough that deposition begins or the runoff water enters a well-defined channel that may be part of a drainage network or a constructed channel [7]. It is important to note that the definition of the slope length relies on the conditions at which the unit plot was constructed by [7]; the unit plot represents the worst condition for the maximum soil erosion case.

Therefore, for the worst condition for the maximum erosion case, the slope length is the shortest distance from the origin of overland flow to the point where deposition takes place or enters stream channels. The slope lengths would rarely have a constant gradient along their entire length, and the slope irregularities would affect the amount of soil movement to the foot of the slope [7]. The slope length factor is given by [7]:

$${L=\left(\frac{\lambda }{{\lambda }_o}\right)}^m$$

where $\lambda$ is the slope length, and ${\lambda }_0$ is the unit plot length = 72.6 ft = 22.13 m. ${\lambda }_0$ is also defined as the horizontal projection of the slope length (e.g., [42,44,45,46]). In one term, the slope steepness factor (S) and slope length factor (L) together are called the topographic factor ($LS$-factor). The topographic factor is the ratio of soil loss per unit area from a field slope length and gradient to that from the 22.1 m length of uniform 9% slope under otherwise identical conditions [7]. Different equations have been suggested at different locations to estimate the topographic factor while taking into account site-specific conditions.

• The topographic factor that was proposed at the topographic condition of USA [7]:

  $$LS={\left(\frac{\lambda }{72.6}\right)}^m(65.41si{n}^2\theta +4.56sin\theta +0.065)$$

  where $\lambda =$ slope length in feet, $\theta =$ angle of the slope, and m = dependent on the slope (0.5 if the slope > 5%, 0.4 if the slope is between 3.5% and 4.5%, 0.3 if the slope is between 1% and 3%, and 0.2 if the slope is less than 1%).
• McCool et al. (1987) improved the $LS$-factor from classic USLE for use in terrain with steeper slopes as cited in [14] for use in RUSLE [42]:

  $$\begin{array}{c}L={\left(\frac{\lambda }{22.13}\right)}^m\hfill \\ m=\frac{sin\theta }{sin\theta +0.269{\left(sin\theta \right)}^{0.8}+0.05}\hfill \\ S=3.0{\left(sin\theta \right)}^{0.8}+0.56for\lambda <4m\hfill \\ S=10.8\ast sin\theta +0.03for\lambda >4mands<9\%\hfill \\ S=16.8sin\theta -0.50for\lambda >4mands>9\%\hfill \end{array}$$

  where $\lambda$ is the slope length in meters, m is the dimensionless parameter, $\theta$ is the angle of field slope in degrees = tan${}^{-1}$ (s/100), and S is the field slope as a percentage.
• Foster et al. (1977) and McCool et al. (1987, 1989) proposed the following equations for the calculation of the $LS$-factors as cited in [34]
  • $L={\left(\frac{\lambda }{72.6}\right)}^m$
  • $m=\frac{\beta }{1+\beta }$ (Foster et al., 1977) as cited in [34]
  • $\beta =\frac{\frac{sin\theta }{0.0896}}{3.0\ast {\left(sin\theta \right)}^{0.8}+0.56}$ (McCool et al., 1989) as cited in [34]
  • $S=10.8\ast sin\theta +0.03$ if the slope (s) is less than 9% (McCool et al., 1987) as cited in [34]
  • $s=16.8sin\theta -0.5$ if the slope is greater than or equal to 9% ( McCool et al, 1987) as cited in [34]
  • $S=3.0\ast {\left(sin\theta \right)}^{0.8}+0.56$ if the slope length is shorter than 4.6 m (McCool et al., 1987) as cited in [34], for the condition where water drains freely from slope end, and it is assumed that inter-rill erosion is insignificant on slopes shorter than 4.6 m [42], where $\lambda$ is the slope length (ft), $\theta$ is the angle of the slope, and m is the dependent on the slope (0.5 if the slope > 5%, 0.4 if the slope is between 3.5% and 4.5%, 0.3 if the slope is between 1% and 3%, and 0.2 if the slope is less than 1%).
  As a remark, when conditions favor more inter-rill and less rill erosion, as in cases of consolidated soils, such as those found in no-till agriculture, m should be decreased by halving the $\beta$ value, where a low rill to inter-rill erosion ratio is typical of the conditions on rangelands [42]. With thawing, and cultivated soils dominated by surface flow, a constant value of 0.5 should be used (McCool et al., 1989, 1993) as cited in [42]. When freshly tilled soil is thawing, in a weakened state and primarily subjected to surface flow, we use the following (McCool et al., 1993) as cited in [42].

  $$\begin{array}{c}S=10.8sin\theta +0.03s<9\%\hfill \\ S={\left(\frac{sin\theta }{0.0896}\right)}^{0.6}s>9\%\hfill \end{array}$$
• The slope factor that is approximately equal to the $LS$-factor at the topographic condition of the Philippines [33].

  $$S=a+b\ast {S_L}^{4/3}$$

  where S is the slope factor, a = 0.1, and b = 0.21; SL is the slope in percent.
• The $LS$-factor was developed for the the topographic condition of Britain [47]:

  $$LS={\left(\frac{\lambda }{22}\right)}^{0.50}\ast \left(0.065+0.045s+0.0065s^2\right)$$

  where $\lambda$ is the slope length (m), and s is the slope steepness (%).
• Apart from the $LS$-factor of the USLE/RUSLE, the Chinese Soil Loss Equation [48] was developed while taking into consideration the Chinese soil environment and topographic conditions (including the modified equation that can calculate $LS$-factor in >10° conditions) [49]. In the Chinese soil loss equation, the $LS$-factor is calculated by [49].

  $$\begin{array}{c}L={\left(\frac{\lambda }{22.1}\right)}^m\hfill \\ m=0.2for\theta \le 1.7\%\hfill \\ m=0.3for1.7\%<\theta \le 5.2\%\hfill \\ m=0.4for5.2\%<\theta \le 9\%\hfill \\ m=0.5for\theta >9\%\hfill \\ S=10.8sin\theta +0.03for\theta <9\%\hfill \\ S=16.8sin\theta -0.05for9\%\le \theta <17.6\%\hfill \\ S=21.9sin\theta -0.96for\theta \ge 17.6\%\hfill \end{array}$$

  where $\lambda$ is the slope length (m), m is the variable slope-length exponent, and $\theta$ is the slope angle (°).
• Other equations of the slope or slope length factor are mentioned in [14,43,48,49,50,51,52].

To estimate the topographic factor ($LS$-factor) for our watersheds, SWATplus is used to define as many hydrologic response units (hrus) as possible to consider an areal distribution of the slope steepness and slope length. In the TxtInOut folder of the SWATplus, area and topography information of each hru are stored in the hru and topography files, respectively. These files are exported to an excel spreadsheet for analysis.

The area, slope, and slope length of each hru are used to estimate the $LS$-factor for each hru using the above equations and Equations (1) and (2) stated below. The weighted average of the $LS$-factors is taken to represent the watershed (it is important to note that the sediment or flow routing techniques in the SWATplus are not employed in this paper due to one or more reasons are stated in Section 1). The best-fit methods are chosen during the calibration of the annual sediment yield (see the calibration stage below).

2.5.4. Estimation of Cover Factor (C-Factor)

The C-factor is the ratio of soil loss from a field with specified cropping to that from clean-tilled, continuous fallow under otherwise similar conditions. These similar conditions are no soil conservation works (land is tilled up and down the slope), soil, slope steepness, slope length, and the rainfall impact pressure is the same for both cropped field and fallow area. The C-factor is related to the land use and land cover, and it is the reduction factor for soil erosion vulnerability [53].

Therefore, the C-factor lies between 0 and 1, which describes the extent of vegetation cover to protect soil from erosion in a given catchment. Its value closer to 0 indicates dense vegetation cover, whereas its value closer to 1 indicates poor vegetation cover. Essentially, surface cover or canopy protects soil erosion by decreasing rainfall impact energy; however, it may have less importance to protect sediment transport from a field.

To some extent, we can say that surface cover affects soil erosion by reducing the transport capacity of the runoff water (Foster, 1982) and by causing deposition in ponded areas (Laflen, 1983) as cited in [34] and also by decreasing the surface area susceptible to raindrop impact [34]. In addition, the plant-root depth and distribution and porosity increase the infiltration rate of rainfall water into the soil, and thus they play a role in reducing soil loss (Jeong et al. 2012) as cited in [54].

Although the C-factor value can be taken from the literature or determined in situ, an extensive literature review compiling the potential soil loss rates of different crop and forest covers compared to likely soil loss rates of bare soil can be used to determine likely C-factor values of a particular site [51]. The published guidelines [7,42], the revised C-factor (Cai et al., 2000) as cited in [55] and the Normalized Difference Vegetation Index [39,40] can be used to compute the C-factor. In our case, the annual or average annual cover factor for each land use category is adopted based on the assessment of the literature.

Researchers [51] reviewed C-factors for the general types of land use and land cover. For our watersheds, the adopted cover factor for each land use is shown in

Table 1. The assigned cover and conservation practice factors for each land use of the watersheds under our consideration.

Land-Use Category	C-Factor	P-Factor
Acacia	0.01	1
Acacia Bushland/Thicket	0.01	1
Acacia Shrubland/Grassland	0.01	1
Agricultural land	0.525	0.52
Bare Land	1	1
Dispersed Acacia	0.01	1
Dispersed Shrub	0.01	1
Eucalyptus	0.001	1
Fir/Cedar Forest	0.001	1
Forest	0.001	1
Forest; Montane broadleaf	0.001	1
Grassland	0.01	1
Grassland, Herbaceous Wetland	0.01	1
Grassland; unstocked (woody plant)	0.01	1
Herbaceous Wetlands	0.01	1
Montane Broadleaf Evergreen Woodland	0.001	1
Rocky Bare Land	1	1
Secondary Semi-deciduous Forest/Woodland	0.001	1
Semi-Desert Grassland with Shrubland	0.01	1
Shrubland	0.01	1
Tropical Forest	0.001	1
Plantations	0.001	1
Tropical Plantations	0.001	1
Urban	0	1
Water Bodies	0	0
Wetland	0.01	1
Woodland	0.01	1

. To estimate an areal weighted average of the cover factor for our watersheds, SWATplus is used to define as many hydrologic response units (hru) as possible to consider the areal distributions of land use and land cover. In the TxtInOut folder of the SWATplus, the area of each hru is stored in the hru file, and the hru’s land use data files are stored in the hru-data file.

These files are exported to an excel spreadsheet for analysis and calculation of the areal weighted average. We can use the shapefile of each land use map (see Figure 9, Figure 10, Figure 11, Figure 12, Figure 13, Figure 14, Figure 15 and Figure 16) to estimate the areal coverage of each land use classes in QGIS, and then the corresponding C and P-factors can be assigned.

2.5.5. Estimation of Soil Conservation/Erosion Control Practice Factor (P-Factor)

The P-factor is the ratio of soil loss associated with a specific support practice to the corresponding soil loss when cultivation is done up and down the slope [34] under otherwise similar conditions. The P-factor describes the effects of practices, such as contouring, strip cropping, concave slopes, terraces, grass hedges, silt fences, straw bales, and subsurface drainage [53].

These conservation practices change the direction and speed of runoff [42]; it mainly reduces the transport of soil particles by blocking runoff and breaking its speed; however, it does not reduce rainfall impact energy to reduce soil erosion. Therefore, the P-factor ranges from 0 to 1, where 0 represents the strong conservation practice (no soil loss from a field is expected), whereas 1 represents the worst condition for the maximum erosion due to lack of conservation practice and when land is tilled up and down the slope, and runoff takes the shortest well-defined channel or route in the field.

The difficulty of accurately mapping support practice factors or not observing support practices leads to many studies ignoring it by giving their P-factor a value of 1.0 [51]. Some P-factors can be ignored if some C- factors already account for the presence of a support factor, such as intercropping or contouring [51]. All non-agricultural lands were also assigned a value of 1 if no feasible conservation measures were applied [29,54,55].

At suitably detailed scales and with enough knowledge of farming practices, using the P-factor may lead to a more accurate estimation of soil loss [51]. Researchers [3] reviewed that considering the temporal variation of the P-factor could significantly improve the performance of the MUSLE, although it has been rarely taken into account.

The soil conservation or erosion control practice factors can be estimated with the help of available tables [7], using land use and land cover maps [29,54,55], and through field measurement (see literature review report in [3]. In our case, the annual soil conservation practice factor for each land use category is adopted based on the assessment of the literature.

Researchers [51] reviewed the P-factors for general types of land use and land cover. The adopted P-factor for land use and land cover category of each watershed is shown in Table 1. The areal weighted average of the P-factor is done in the same way as for the cover factor.

2.5.6. Estimation of Coefficient a and Exponent b through Calibration

For a chosen value of the exponent b, the best-fit corresponding value of the coefficient a is estimated through calibration. The selection of the best exponent and the best equation among listed above and below (see Equations (1) and (2)) for the topographic factor is done after calibration of observed and simulated sediment (i.e., the MUSLE is used to estimate sediment load). Figure 20 shows sample graphs of the calibrated sediment when the topographic factor is calculated using the equation that was proposed by [7].

During calibration, the Nash–Sutcliffe efficiency corresponds to each $LS$-factor, the exponent b and the coefficient a is evaluated, and graphs of the exponent b versus the Nash–Sutcliffe efficiency, and graphs of the coefficient a versus exponent b are drawn for each watershed, as shown in Figure 21, Figure 22, Figure 23, Figure 24, Figure 25, Figure 26 and Figure 27. For a chosen value of b, we test seven different equations of the topographic factor for each watershed. Therefore, we can have as many graphs as possible.

$$LS=(0.02222J^{1.5}+0.03231J+0.1004)\ast 0.2901△y^{0.4002}forJ<5\%$$

(1)

$$LS=(0.02222J^{1.5}+0.03231J+0.1004)\ast 0.2105△y^{0.5004}forJ>5\%$$

(2)

where J is the slope in %, and $△y$ is the slope length. For further description, readers are encouraged to watch the video at https://www.youtube.com/watch?v=w6w8jxbTJfo (accessed on 25 February 2021). For the case of the watersheds under our consideration, we take △y/22.1 as the field slope length.

2.6. Verifying the Best Exponent of the Modified Universal Soil Loss Equation

As we discussed in Section 1, the yearly simulation time step is preferred to address the gradual processes of soil erosion and sediment transport. It is important to prove whether a change in the simulation time step changes the coefficient and the exponent of the MUSLE or not. This approach leads us to find the best exponent of the MUSLE.

Proof. 

If we consider the small simulation time step and the small simulation period, we can maintain the temporal variation of the factors, which directly affects the soil erosion process. For a given field, no change in the cover, conservation practice, and soil erodibility factors of the MUSLE will be expected at the small simulation period. At the end of the simulation period, only in variation of the coefficient and the exponent of the MUSLE with the simulation time step affect sediment yield output (see the proof steps below for a change in runoff and the peak runoff rate).

If the variations of the coefficient and the exponent of the MUSLE with a small change in the simulation time step are detected, then the variations of the coefficient and the exponent with any other simulation time step are confirmed. For the sake of starting, let us consider one and two unit simulation time steps and two unit simulation period; no change in the factors of the MUSLE will be expected for about two unit simulation period. Therefore, soil loss from a field at the one unit simulation time step for about two unit simulation period, is equal to the sum of soil loss at the end of the first and next one unit time;

$$a_1{\left(Q_1{q}_1\right)}^{b_1}KLSCP+a_1{\left(Q_2{q}_2\right)}^{b_1}KLSCP$$

where suffixes 1 and 2 correspond to the runoff volume (Q), and the peak runoff rate (q) indicates the first and second simulation at the one unit simulation time step or interval. We note that K, L, S, C, and P are the same for the two unit simulation period, and the coefficient and the exponent are the same at the one unit simulation time step. Soil loss from the field at the two unit simulation time step for about a two unit simulation period;

$$a_2{\left((Q_1+Q_2)q_1\right)}^{b_2}KLSCP\mathrm{if}\mathrm{the}\mathrm{peak}\mathrm{runoff}\mathrm{rate}\mathrm{is}{q}_1$$

$$a_2{\left((Q_1+Q_2)q_2\right)}^{b_2}KLSCP\mathrm{if}\mathrm{the}\mathrm{peak}\mathrm{runoff}\mathrm{rate}\mathrm{is}{q}_2$$

where $a_2$ and $b_2$ indicate a value of the coefficient (a) and exponent (b) at the two unit simulation time step. We note that the total runoff volume (Q) at the end of the two unit simulation period, is equal to the sum of the runoff volumes at the end of the first and next one unit time; the peak runoff rate will be expected before one unit time or between 1 and 2 unit time.

In either case, sediment yield is the same. Therefore,

$a_1{\left(Q_1{q}_1\right)}^{b_1}KLSCP+a_1{\left(Q_2{q}_2\right)}^{b_1}KLSCP=a_2{\left((Q_1+Q_2)q_1\right)}^{b_2}KLSCP$ if the peak runoff rate is $q_1$

$$a_1{\left(Q_1{q}_1\right)}^{b_1}+a_1{\left(Q_2{q}_2\right)}^{b_1}=a_2{\left((Q_1+Q_2)q_1\right)}^{b_2}$$

If there is no variation of the coefficient and exponent with small variation in simulation time step, then $a_1=a_2=a$ and $b_1=b_2=b$

$${\left(Q_1{q}_1\right)}^b+{\left(Q_2{q}_2\right)}^b={\left((Q_1+Q_2)q_1\right)}^b$$

(3)

In the same way,

$a_1{\left(Q_1{q}_1\right)}^{b_1}KLSCP+a_1{\left(Q_2{q}_2\right)}^{b_1}KLSCP=a_2{\left((Q_1+Q_2)q_2\right)}^{b_2}KLSCP$ if the peak runoff rate is $q_2$

$$a_1{\left(Q_1{q}_1\right)}^{b_1}+a_1{\left(Q_2{q}_2\right)}^{b_1}=a_2{\left((Q_1+Q_2)q_2\right)}^{b_2}$$

If there is no variation of the coefficient with small variation in simulation time step, then $a_1=a_2=a$ & $b_1=b_2=b$

$${\left(Q_1{q}_1\right)}^b+{\left(Q_2{q}_2\right)}^b={\left((Q_1+Q_2)q_2\right)}^b$$

(4)

Equations (3) and (4) are false for a given value of the exponent b. In this case, the coefficient and the exponent of the MUSLE change as a change in simulation time step for a given total simulation period. Equations (3) and (4) hold true when $b=1$ and $q_1=q_2$, and for other values of the exponent b and $q_1=q_2$, it is false. This implies that only one peak runoff rate is possible per storm event (i.e., from the beginning of runoff to the end of the runoff from a slope-field). This means that sediment is transported from the beginning of runoff to the end of the runoff; the objective of the MUSLE is to estimate the total sediment load transported from the beginning of runoff to the end of the runoff.

Therefore, the best theoretical exponent of the MUSLE is 1. This is a theoretical exponent because the left and right sides of Equations (3) and (4) represent the theoretical linked expressions without knowledge of observed sediment. The actual exponent of the MUSLE is estimated by applying the model at selected watersheds. From all graphs (see Figure 21, Figure 22, Figure 23, Figure 24, Figure 25, Figure 26 and Figure 27), the best actual exponent of the MUSLE is 1, which results in a Nash–Sutcliffe efficiency of approximately 1 irrespective of the topographic factor and the three watershed sizes (Hombole, Mojo, and Gumera watersheds). Therefore, the best exponent of the MUSLE is 1. □

3. Results

We confirmed that the best exponent of the MUSLE is 1 irrespective of the topographic factor, which results in the maximum performance of the MUSLE (i.e., approximately 100%). From all graphs (see Figure 21, Figure 22, Figure 23, Figure 24, Figure 25, Figure 26 and Figure 27), if we consider one watershed, we take the exponents and topographic factors, which result in the maximum Nash–Sutcliffe efficiency; however, if we consider two or more watersheds, we take the exponents and topographic factors, which result in the minimum Nash–Sutcliffe efficiency.

Accordingly, the best exponent of the MUSLE is 0.57, which results in a Nash–Sutcliffe efficiency of approximately 0.8 if the topographic factor is calculated using Equations (1) and (2). Therefore, this is the best combination of the exponents and topographic factors of the MUSLE under the hydro-climatic conditions of all watersheds under our consideration.

To determine the best combination of the exponents and topographic factors, the important relationships between the coefficient a and exponent b, the exponent b and the Nash–Sutcliffe efficiency are drawn for the future evaluation of the MUSLE at any watershed. As we can see from the graphs (see Figure 21, Figure 22, Figure 23, Figure 24, Figure 25, Figure 26 and Figure 27), for observed and simulated sediment, as the relationship between the coefficient a and exponent b approaches to power or logarithmic function; the relationship between the exponent b and the Nash–Sutcliffe efficiency approaches to a quadratic function. This relationship can be used to find the best performance of the MUSLE during the calibration of the model.

4. Discussion

Based on our evaluation of the soil erodibility equations, we found that the best equation to estimate the soil erodibility factor was the Williams (1995) equation as cited in [35]. We considered land use maps to assign a value for the cover and conservation practice factors from the past experiences from the literature, and the coefficient a was estimated through calibration. Since only a product effect of the coefficient, soil erodibility, cover, and conservation practice factors are reflected in the MUSLE rather than their individual effect during the calibration of the sediment yield, any change in these factors affects the coefficient of the MUSLE.

Therefore, we do not like to suggest strict procedures to estimate these factors. It is highly preferable if these factors are measured and studied at a temporal and spatial scale to understand their effect on soil erosion in a particular field. This is because the soil erodibility, cover, and conservation practice factors of the MUSLE reflect site-specific conditions. For example, we can talk about the density and pattern of land cover, nature and extent of soil conservation and flood protection work, and the temporal variation of soil properties.

For all watersheds under our consideration, the best exponent of the MUSLE was 0.57, which resulted in a Nash–Sutcliffe efficiency of 0.8 if the topographic factor was calculated using Equations (1) and (2). In this case, the proposed exponent of the model is different from its original exponent (0.56), while for other exponents and topographic factors, the performance of the model decreases. For example, the best exponent of the model is 0.56, which results in a Nash–Sutcliffe efficiency of 0.78 if the topographic factor is calculated using the equation that was proposed by McCool et al. (1987) as cited in [14]. In this case, the proposed exponent is the same as the original exponent of the MUSLE (0.56); however, the performance of the MUSLE decreases as compared to the previous one. Therefore, the performance of the MUSLE is very good for the previous case.

The performance of the MUSLE was tested at a watershed scale using directly measured flow data; it showed good performance (i.e., the performance of the MUSLE is greater than or equal to 80%) for all four watersheds under our consideration provided that the exponents and topographic factors of the original MUSLE were changed. This result supports the literature review report that the model shows better performance at a watershed scale than a plot scale and if it is applied using directly measured runoff data [3].

This also supports the conclusions of some authors as the MUSLE has been observed to give good results in various applications in some parts of tropical Africa (Ndomba, 2007) as cited in [5]. The MUSLE has also been successfully demonstrated in sub-Saharan Africa [5]. In addition, it also supports the experimental plot result of sheet erosion at the Enerta study site in Ethiopia, where the MUSLE was better at estimating soil loss from a cultivated field than the USLE [6].

5. Conclusions

We verified that the best exponent of the MUSLE was 1, which resulted in the maximum performance of the MUSLE. The performance of the MUSLE was greater than or equal to 80% for all four watersheds under our consideration. We expect the same for other watersheds of Ethiopia provided that the exponent of the model is 0.57 and that its topographic factor is calculated using Equations (1) and (2). This can be taken as the best combination of the exponents and topographic factors under the hydro-climatic conditions of Ethiopia. We recommend further investigation of the best combination of the exponents and topographic factors by applying the MUSLE at different watersheds of Ethiopia.

We suggested the best equations of the topographic factor for the conditions of Ethiopia. In the MUSLE, the topographic factor is directly proportional to the soil erosion and sediment yield. However, as the slope length becomes increasingly larger, there is a possibility that soil erosion from the upper part of the slope becomes deposited at the lower part of the slope. Therefore, more research works are required to understand the effect of the slope length on soil erosion and sediment transport.