1. Introduction
The commonly used soil erosion models are the Universal Soil Loss Equation (USLE), Revised Universal Soil Loss Equation (RUSLE), and Modified Universal Soil Loss Equation (MUSLE). Physically speaking, the MUSLE is more appropriate than the USLE or RUSLE [
1]. Williams (1975) developed the MUSLE using 778 storm runoff events collected from 18 small watersheds [
2,
3], with areas varying from 15 to 1500 ha, slopes from 0.9 to 5.9%, and slope lengths of 78.64–173.74 m (Hann et al. 1994), as cited in [
4]. The MUSLE is given by
where
y is the sediment yield in metric tons,
a is the coefficient,
b is the exponent (
and
for USA, where the MUSLE was originally developed),
Q is the runoff volume in m
,
q is the peak runoff rate in m
s
,
K is the soil erodibility factor in
,
L is the slope length factor,
S is the slope steepness factor,
C is the cover factor, and
P is the soil conservation practice factor. For the applications of the MUSLE at different locations, different equations of the topographic and soil erodibility factors have been suggested. To apply the MUSLE to a large watershed (for our case, the Hombole Watershed is 762,281 ha, the Mojo Watershed is 150,282 ha, the Gumera Watershed is 127,805 ha, and the Gilgel Gibe 1 Watershed is 292,809 ha), one approach that was proposed was using the MUSLE in the Soil and Water Assessment Tool (SWAT) environment. This may be because sediment yield can be more accurately estimated if the large watershed is divided into subwatersheds (area < 2590 ha) to compensate for non-uniformly distributed sediment sources. The effect of the watershed hydraulics and sediment particle size can be included by routing the sediment yield from the subwatersheds to the large watershed [
2]. As part of evaluation of the model, we considered the specific behavior of the MUSLE, the experiences of some 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, it was found that the MUSLE showed better performance in the case of directly measured flow data rather than indirectly obtained flow data using indirect methods [
4], and the model also provided appropriate estimates at a watershed rather than an experimental plot, as was reviewed and reported by [
4]. In this connection, if we consider the SWAT, it uses an indirect method (like the Soil Conservation Service Curve Number) to generate the runoff and then uses the MUSLE to estimate the soil loss from a hydrologic response unit (which is similar to a plot scale). Then, the SWAT routes the 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 the channel, channel depth, and width in the SWAT environment.
If we consider the experiences of some other researchers, the MUSLE is unsuitable for prediction of the sediment yield for small storms [
5]. However, the slight variation in the hydrological response of a watershed in terms of sediment yield might be changing in the antecedent hydrological conditions, the spatial and temporal distribution of rainfall, and the availability of eroded sediment throughout the watershed, which is not taken into account by the MUSLE as it is for many other lumped models [
5].
If we consider the physical connection between the factors of the MUSLE, as far as runoff energy for soil detachment and sediment transport are concerned, the physical connection between the runoff, soil erodibility, topographic factor, cover, and soil conservation practice factors is convincing, but further refining the physical connection between the factors may become necessary. For instance, the cover and soil conservation factors play a role in breaking runoff energy so as to protect soil loss due to runoff. As the slope length becomes larger and larger, there is a possibility that erosion from the upper part of the slope gets 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 gets concentrated on its way to the bottom of the slope. In other words, the energy of the runoff decreases as resistance against the 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 [
6], and it has been successfully demonstrated in Sub-Saharan Africa [
6]. As per the experimental plot result of sheet erosion at the Enerta study site in Ethiopia, the MUSLE was better at estimating soil loss from a cultivated field than the USLE [
7].
Therefore, based on the above limitations and advantages, the MUSLE should be tested at a watershed scale rather than a plot scale under the hydroclimatic 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. The MUSLE considers sheet to rill erosion, but we cannot exactly tell whether it considers gully erosion or not. It follows that improving the MUSLE becomes important for possible application in data-scarce areas. The underlying physical assumption to improve the MUSLE is that the amount of potential energy in the runoff is proportional to the shear stress for sediment transport from a slope field and the kinetic energy of the runoff at the bottom of the slope field for gully formation. In the improved MUSLE, the peak runoff rate is eliminated from the variables of the original MUSLE. For the evaluation of the improved MUSLE, we consider the following cases.
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 that soil erosion, sediment transport, deposition, consolidation, and re-suspension are quite complex processes which depend on the physical, biological, mechanical, and chemical activities within a large heterogeneous watershed [
1]. Due to these complex processes, the soil that was eroded at an unknown last time can be transported, deposited, consolidated, resuspended, and reach an outlet at a different time. Therefore, the measured sediment at the outlet at the current time may not reflect the current information about the watershed; it rather reflects 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 the formulation in [
8], we can conclude that the annual soil erodibility is the average soil erodibility ranging from loose to compacted soil due to rainfall impact. As the soil erodibility factors of the USLE, MUSLE, and improved MUSLE are the same, the annual time step is preferred over the daily time step. The annual simulation time step enables taking into account gully erosion (Gully erosion is usually estimated on an annual basis [
9], and it is important to note here is that gully erosion is a common problem in Ethiopia (e.g., [
10,
11,
12,
13,
14])) to take into account gradual soil erosion processes and gradual changing activities like 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 1-year full cycle.
If we consider the hydrologic response unit (hru) in the SWAT environment, as the number of hru becomes larger and larger, we better take into account the spatial variability of land use, soil, and slope all over the watershed. To test the improved MUSLE at a watershed scale, the sediment or flow routing in the stream channels of the SWAT is 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 the soil, cover, conservation practice, and topography.
If we consider the calibration parameters, choosing the appropriate calibration parameters is essential for modeling. For a given uniform watershed, 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 improved MUSLE, and this is because only a product effect of the coefficient and these factors is reflected in the improved MUSLE rather than their individual effect during the calibration of the sediment yield. Compared with the other parameters of the improved MUSLE, the individual effect of the exponent of the improved MUSLE is reflected during calibration of the sediment yield. Therefore, estimating the exponent of the improved MUSLE through calibration is more feasible than through other parameters of the improved 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 improved MUSLE is applied at different watersheds. From this explanation, the independent effect of the exponent and topographic factor of the improved MUSLE can be seen by applying the model at different watersheds. In general, runoff and the sediment data reflect the hydroclimatic conditions of a particular area, which independently affect the overall calibration process. Therefore, our main task is to estimate the best exponent of the improved MUSLE and the best combination of the exponent and topographic factor of the improved MUSLE by applying the model at different watersheds of Ethiopia. For the sake of the calibration procedure, the main factors of the improved MUSLE which directly affect 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 which do not directly affect the soil erosion process or which have no direct physical meaning (i.e., coefficient a and exponent b) are estimated through calibration.
3. Improving MUSLE
The assumption of the USLE and RUSLE is that the rainfall intensity is uniformly distributed over a slope field area. Likewise, in the case of the MUSLE, the runoff is uniformly distributed over the slope field area. The total runoff volume is considered for sediment yield modeling. In both cases, the slope length represents the worst condition for the maximum erosion case, and therefore, the slope length is the shortest distance from the origin of the runoff to the bottom of the slope. The underlying physical assumption to improve the MUSLE is that the amount of potential energy of the runoff is proportional to the shear stress for sediment transport from a slope field and the kinetic energy of the runoff at the bottom of the slope field for gully formation. The potential energy of an object due to its position is given by
where
is the potential energy,
m is the mass of the object,
g is the acceleration due to gravity,
h is the height referring to the position of the object from a reference point,
is the density of the object, and
v is the volume of the object.
In the process of all of the rainfall in the slope field area, the total potential energy of the runoff is given by [
23]
where
is the total potential energy of the runoff (
J),
is the density of water (kg m
−3),
g is the acceleration due to gravity (m s
−2),
B is the slope width (m),
L is the slope length (m),
is the slope angle, and
Q is the runoff volume (m).
Equation (
1) can be applied when the runoff volume flows down from the top of the slope field to the bottom of the slope field, and a temporal variation of the runoff volume is considered. However, this equation does not consider every other runoff volume that will begin at every slope height and flows down simultaneously to the bottom of the slope (i.e., the runoff that will start from any point on the entire slope surface due to the uniform distribution of rainfall on the entire slope field), and it also does not consider a temporal variation of every other runoff volume. Therefore, let us assume that every runoff volume due to raindrops on the entire slope length flows down to the bottom of the slope. Let us say the position of first runoff volume, as well as those of the second, third, and so on along the length of the slope are at the slope heights
h,
,
, and so on from the bottom of the slope, respectively. One runoff volume takes over the position of another runoff volume as it flows down to the bottom of the slope.
The total potential energy of the first runoff volume due to its changes in position as it flows down from the height
h to the bottom of the slope is equal to
:
The total potential energy of the second runoff volume due to its changes in position as it flows down from the height
(let us say just immediately after the first runoff volume) to the bottom of the slope is equal to
:
The total potential energy of the third runoff volume due to its changes in position as it flows down from the height
(let us say just immediately after the second runoff volume) to the bottom of the slope is equal to
, and so on:
Therefore, the total potential energy of all runoff volumes due to their change in position along the length of the slope to the bottom of the slope is equal to
:
Then, we can substitute Equations (
2)–(
4) into Equation (
5):
where
and
are the heights of the slope just immediately after heights
h and
, respectively, and so on:
Each runoff volume (
) in a very small area (
) changes in time (
) with the rainfall intensity (
) and soil infiltration rate (
):
where
is the runoff depth:
Then, we can substitute Equation (
10) into Equation (
9):
It is noteworthy that only the depth of the runoff volume changed with a small change in time. If it rained continuously for some duration of time, then the runoff volume increased at every point of the slope length. Therefore, the total runoff volume at a particular point of the slope field is a function of the runoff depth changing with time, where the bottom area of the runoff volume does not change with time:
Substituting Equation (
12) into Equation (
8) yields
For our case, the total potential energy per a unit area (
) is important parameter for soil erosion and sediment transport along the length of the slope, and therefore
Substituting Equation (
13) into Equation (
14) yields
The bottom area of each runoff volume is constant, and therefore
where
Q is the runoff volume in m.
Substituting Equation (
16) into Equation (
15) yields
We can evaluate the integral as follows:
The trigonometric relationship between the slope length, slope angle, and height is given by
Substituting Equation (
19) into Equation (
18) yields
We should note that as every runoff volume flows down simultaneously to the bottom of the slope, the runoff concentrates at the bottom of the slope, and it is proportional to the rainfall intensity. The runoff concentration leads to gully formation at the bottom of the slope. The total potential energy of the runoff is proportional to the rainfall intensity.
We can compare the total potential energy of the runoff per unit area (Equation (
20)) with the MUSLE:
From the relationship in Equation (
21), we can reveal the following.
In the MUSLE, is the term which contributes to the energy of the runoff, whereas K, C, and P are the coefficients which contribute to the energy dissipation (i.e., K, C and P can be taken as friction resistances against the flow, and the values of K, C, and P vary from 0–1).
Therefore, we can use
in place of
in the MUSLE:
For a given watershed, the topography does not change over time (i.e., the slope length (
L) and slope angle (
) are constant). Therefore, the energy of the runoff is proportional to the sediment yield by defining the proportionality constants:
If the topographic factor of the MUSLE is considered, then
where
is the constant.
Therefore, we can use
in place of
in Equation (
23):
Therefore, we call Equation (
25) the improved MUSLE.
3.1. Estimating the Theoretical Exponent of the Improved MUSLE
As we discussed in the Introduction, 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 improved MUSLE or not. This approach led us to estimate the theoretical exponent of the improved 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 affect the soil erosion process. For a given field, no change in the cover, conservation practice, or soil erodibility factors of the improved MUSLE will be expected in the small simulation period. At the end of the simulation period, only variation of the coefficient and the exponent of the improved MUSLE with the simulation time step affects the sediment yield output (see the proof steps below for a change in the runoff volume (
Q) and the peak runoff rate (
q)). If the variations of the coefficient and the exponent of the improved 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 the start, let us consider 1 and 2 unit simulation time steps and a 2-unit simulation period. No change in the factors of the improved MUSLE will be expected for about a 2-unit simulation period. Therefore, the soil loss from a field at the 1-unit simulation time step for about a 2-unit simulation period is equal to the sum of the soil loss at the end of the first and next 1-unit time:
where suffixes 1 and 2 correspond to the runoff volume (
Q), indicating the first and second simulation at the 1-unit simulation time step or interval. It is noteworthy that
, and
P are the same for the 2-unit simulation period. The coefficient and the exponent are the same at the 1-unit simulation time step. □
The soil loss from the field at the 2-unit simulation time step for about a 2-unit simulation period is expressed as
where
and
indicate a value of the coefficient (
a) and exponent (
b) at the 2-unit simulation time step. It is noteworthy that the total runoff volume (
Q) at the end of the 2-unit simulation period is equal to the sum of the runoff volumes at the end of the first and next 1-unit time.
In either case, the sediment yield is the same, and therefore
If there is no variation in the coefficient or exponent with small variation in the simulation time step, then
Equation (
26) is false if
or
. In this case, the coefficient and the exponent of the improved MUSLE change with a change in the simulation time step for a given total simulation period. Equation (
26) holds true when
. This implies that we consider the total runoff volume per storm event (i.e., from the beginning of runoff to the end of the runoff from a slope field) to estimate the total sediment load. The objective of the improved MUSLE is to estimate the total sediment load transported from the beginning of runoff to the end of the runoff from the slope field. Therefore, the theoretical exponent of the improved MUSLE is one. It is a theoretical exponent because the left and right sides of Equation (
26) represent the theoretical linked expressions without knowledge of the observed sediment. With knowledge of the observed sediment, the actual exponent of the improved MUSLE can be obtained by applying the model at the selected watersheds of Ethiopia.
3.2. Estimating the Factors of the Improved MUSLE
The improved MUSLE shares the same factors with the USLE, RUSLE, and MUSLE. The original factors of the USLE represent the average value for estimating the annual sediment yield. The unit plot [
8] represents the worst case for the maximum soil erosion in 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, cover, and conservation practice factors depend on the spatial resolution of the digital elevation model (DEM) and soil and land use maps, respectively. Therefore, in the actual sediment modeling, the average field slope length [
24] and slope steepness, or simply the topographic factor [
25], average runoff, average soil erodibility factor [
26], and average cover and conservation practice factors, are taken.
3.2.1. Estimation of Runoff Factor
In the MUSLE, the runoff factor is the product of the total runoff volume and peak runoff rate. In the improved MUSLE, the peak runoff rate is eliminated. To estimate the runoff factor of the improved MUSLE, the volume of the runoff can be obtained by direct measurement of the runoff on a storm event basis and also by using indirect methods such as the Soil Conservation Service Curve Number (SCS CN) method, rational method, flood routing, or unit hydrograph. For our case, we used the daily average discharge to estimate the annual total runoff volume for the annual sediment yield estimation. The reasons for why we used directly measured flow data and why we estimated the annual sediment yield are addressed in the Introduction.
3.2.2. Estimation of Soil Erodibility Factor (K-Factor)
The authors in [
8] defined the 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 feet in length uniform 9% slope continuous in a clean-tilled fallow. It is the continuous fallow that is tilled up and down the slope. The soil erodibility factor is given by [
8]
where
A is the event soil loss from the unit plot in
,
E is the storm kinetic energy in 100 ∗ foot-tons per acre,
is the maximum 30-min intensity in inches per hour, and
K is the soil erodibility factor in
. 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, a lesser soil erosion condition that is different from the worst condition takes into account the soil cover and conservation practice 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 type of soil and the quite complex interaction of physical, biological, chemical, and mechanical processes. From the soil erodibility table or equations, 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 under the same rainfall impact pressure in otherwise similar soil erosion conditions. From this range of the soil erodibility factor, we can conclude that soil erodibility refers to the degree of ease in eroding 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:
The procedure to test the above soil erodibility equations on the basis of the original definition of the soil erodibility are given in [
1]. According to the procedure, we used Williams’s (1995) equation as cited in [
32] to calculate the soil erodibility factor using the soil data of the watersheds under our consideration: the watersheds of Ethiopia in general.
3.2.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 the soil loss from the 9% slope under otherwise identical conditions [
38]. A high rate of soil loss is associated with steep slopes [
26,
39], and soil loss prediction is more sensitive to the slope steepness than slope length [
40].
The 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 [
8]. It is important to note that the definition of the slope length relies on the conditions at which the unit plot was constructed according to [
8]. 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 the overland flow to the point where deposition takes place or enters the 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 [
8]. The slope length factor is given by [
8]:
where
is the slope length,
is the unit plot length of 72.6 ft = 22.13 m, and
is also defined as the horizontal projection of the slope length (e.g., [
39,
41,
42,
43]). In one term, the slope steepness factor (S) and slope length factor (L) together is 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 a uniform 9% slope under otherwise identical conditions [
8]. 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 for the topographic conditions in the USA [
8]:
where
the slope length in feet,
the angle of the slope, and
m is dependent on the slope (0.5 if slope >5%, 0.4 if slope is between 3.5% and 4.5%, 0.3 if slope is between 1% and 3%, and 0.2 if slope is less than 1%).
McCool et al. (1987) improved the LS factor from the classic USLE for use in terrain with steeper slopes, as cited in [
25], for use in the RUSLE [
39]:
where
is the slope length in meters,
m is the dimensionless parameter,
is the angle of the 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 [
31]:
;
(Foster et al., 1977), as cited in [
31];
(McCool et al., 1989), as cited in [
31];
if the slope (s) is less than 9% (McCool et al., 1987), as cited in [
31];
if the slope is greater than or equal to 9% ( McCool et al., 1987), as cited in [
31].
if the slope length is shorter than 4.6 m (McCool et al., 1987), as cited in [
31], for the condition where water drains freely from the slope end, and it is assumed that inter-rill erosion is insignificant on slopes shorter than 4.6 m [
39]. Here,
is the slope length (ft),
is the angle of the slope, and
m is dependent on the slope (0.5 if 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 the conditions favor more inter-rill and less rill erosion, as in the cases of consolidated soils like those found in no-till agriculture,
m should be decreased by halving the
value, where a low rill to inter-rill erosion ratio is typical of the conditions in rangelands [
39]. 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 [
39]. When freshly tilled soil was thawing, in a weakened state, and primarily subjected to surface flow, we used the following (McCool et al., 1993), as cited in [
39]:
The slope factor which is approximately equal to the LS factor in the topographic conditions of the Philippines [
30]:
where
S is the slope factor,
a = 0.1,
b = 0.21, and
is the slope in percent.
The LS factor was developed in the topographic conditions of Britain [
44]:
where
is the slope length (m) and
s is the slope steepness (%).
Apart from the LS factor of the USLE and RUSLE, the Chinese Soil Loss Equation [
45] was developed while taking into consideration the Chinese soil environment and topographic conditions (including the modified equation that can calculate the LS factor in >10° conditions) [
46]. In the Chinese soil loss equation, the LS factor is calculated by [
46]
where
is the slope length (m), m is the variable slope length exponent, and
is the slope angle (°).
Other equations for the slope or slope length factor are mentioned in [
25,
40,
45,
46,
47,
48,
49].
To estimate the topographic factor (LS) for our watersheds, SWATplus was 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 SWATplus, the area and topography information of each hru were stored in the hru and topography files, respectively. These files were exported to an Excel spreadsheet for analysis. The area, slope, and slope length of each hru were used to estimate the LS factor for each hru by using the above equations and Equations (
27) and (
28), stated below. The weighted average of the LS factors was taken to represent the watershed. The best fit methods were chosen during calibration of the annual sediment yield (see the calibration stage below).
3.2.4. Estimation of the Cover Factor (C Factor)
This is the ratio of soil loss from a field with specified cropping to that from a clean-tilled continuous fallow under otherwise similar conditions. These similar conditions are no soil conservation works (land is tilled up and down the slope) and the soil, slope steepness, slope length, and rainfall impact pressure being the same for both the cropped field and fallow area. The C factor is related to land use and land cover, and it is the reduction factor to soil erosion vulnerability [
50]. Therefore, the C factor lies between 0 and 1, which describes the extent of vegetation cover to protect the soil from erosion in a given catchment. Its value being closer to zero indicates dense vegetation cover, whereas its value being closer to one indicates poor vegetation cover. Essentially, the surface cover or canopy protects from soil erosion by decreasing the rainfall impact energy, but it may have less importance in protecting sediment transport from a field. To some extent, we can say that the 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 [
31], and also by decreasing the surface area susceptible to raindrop impacts [
31]. In addition, the plant root depth and distribution as well as the 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 [
51].
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 the likely soil loss rates of bare soil can be used to determine the likely C factor values of a particular site [
48]. The published guidelines [
8,
39], the revised C factor (Cai et al., 2000) as cited in [
52], and the Normalized Difference Vegetation Index [
36,
37] can be used to compute the C factor. For our case, the annual or average annual cover factor for each land use category was adopted based on the assessment of the literature. The authors in [
48] reviewed the 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. To estimate an areal weighted average of the cover factor for our watersheds, SWATplus was 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 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 A4,
Figure A5,
Figure A6,
Figure A7,
Figure A8,
Figure A9,
Figure A10 and
Figure A11) to estimate the areal coverage of each land use class in QGIS, and then the corresponding C and P factors can be assigned.
3.2.5. Estimation of Soil Conservation or Erosion Control Practice Factor (P-Factor)
This is the ratio of soil loss associated with a specific support practice in the corresponding soil loss when cultivation is performed up and down the slope [
31] 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 [
50]. These conservation practices change the direction and speed of runoff [
39]. It mainly reduces the transport of soil particles by blocking runoff and breaking its speed, but it does not reduce rainfall impact energy to reduce soil erosion. Therefore, the P factor ranges from 0 to 1, where 0 represents strong conservation practice (no soil loss from a field is expected) and 1 represents the worst conditions for maximum erosion due to a lack of conservation practice and when the land is tilled up and down the slope such that 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 this by giving their P factor a value of 1.0 [
48]. Some P factors can be ignored if some C factors already account for the presence of a support factor such as intercropping or contouring [
48]. All non-agricultural lands were also assigned a value of one if no feasible conservation measures were applied [
26,
51,
52]. 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 [
48]. The authors in [
4] revealed 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 [
8], using land use and land cover maps [
26,
51,
52], and through field measurement (see the literature review report in [
4]). For our case, the annual soil conservation practice factor for each land use category is adopted based on the assessment of the literature. The authors in [
48] reviewed the P factors for general types of land use and land cover. The adopted P factor for the land use and land cover category of each watershed is shown in
Table 1. The areal weighted average of the P factor was found in the same manner as the cover factor.
3.2.6. Estimation of Coefficient a and Exponent b through Calibration
For a chosen value of the exponent
b, the best fit corresponding value of coefficient
a was estimated through calibration. The selection of the best exponent and the best equation among those listed above and below (see Equations (
27) and (
28)) for the topographic factor was performed after calibration of the observed and simulated sediment (i.e., the improved MUSLE was used to estimate the sediment load).
Figure 7 shows sample graphs of the calibrated sediment when the topographic factor was calculated using the equation that was proposed in [
8].
During calibration, the Nash–Sutcliffe efficiency corresponds to each LS factor, exponent
b and coefficient
a are evaluated, and graphs of exponent
b versus the Nash–Sutcliffe efficiency and coefficient
a versus exponent
b are drawn for each watershed, as shown in
Figure 8,
Figure 9,
Figure 10,
Figure 11,
Figure 12,
Figure 13 and
Figure 14. For a chosen value of
b, we tested seven different equations of the topographic factor for each watershed. Therefore, we could have as many graphs as possible.
where
J is the slope in percent and
is the slope length. For a 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 took △
y/22.1 as the field slope length.