Next Article in Journal
Spatial Variations and Distribution Patterns of Soil Salinity at the Canal Scale in the Hetao Irrigation District
Previous Article in Journal
Analysis of the Role of Precipitation and Land Use on the Size of the Source Area of Shallow Landslides
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

A Mathematical Method for Estimating the Critical Slope Angle of Sheet Erosion

1
College of Environmental & Resource Sciences, Zhejiang University, Hangzhou 310058, China
2
Zhejiang Provincial Key Laboratory of Agricultural Resource and Environment, Zhejiang University, Hangzhou 310058, China
3
Ministry of Education Key Laboratory of Environment Remediation and Ecological Health, Zhejiang University, Hangzhou 310058, China
*
Author to whom correspondence should be addressed.
Water 2023, 15(19), 3341; https://doi.org/10.3390/w15193341
Submission received: 21 August 2023 / Revised: 14 September 2023 / Accepted: 21 September 2023 / Published: 23 September 2023
(This article belongs to the Section Water Erosion and Sediment Transport)

Abstract

:
Estimating the critical slope angle (CSA) for sheet erosion is important for the precision estimation of sheet erosion and the development of erosion control practices. This study developed mathematical equations considering rainfall intensity and soil infiltration to efficiently estimate both instantaneous (at a given instant during rainfall) and cumulative CSAs, while also providing a valuable explanation for the change in CSA. The mathematical equations were consistent with observations from runoff plots (NSE = −1.01) of loess soils from Zhangjiakou (China) and simulation results (NSE = 0.96) from the Water Erosion Prediction Project model for a loam soil in Montana (USA). Estimated instantaneous CSA determined by the mathematical equations increased as the ratio of rainfall intensity to soil infiltration (I/f) increased, resulting in higher observed cumulative CSA after heavy versus normal rainfall events. Heavy rainfall, compacted soil, and varying rainfall duration affected the CSA by changing the I/f ratio. Maximum instantaneous CSA provided a better prediction of changes in soil erosion dynamics than that obtained from CSAs determined by field observations or experimental simulations. The mathematical equations illustrate the underlying physical mechanisms by which rainfall intensity and soil infiltration affect the CSA through changing the shear stress of overland flow. The results of this study provide critical information for guiding the development of effective soil erosion control strategies.

1. Introduction

Soil erosion is a persistent global environmental threat causing fertilizer loss, soil degradation, and a reduction in soil water/nutrient storage capacity, which in turn leads to a deterioration of water quality, decreased land productivity, and enhanced flooding [1,2,3]. Although many countries have adopted a series of soil erosion control measures, such as returning farmland to forests, building terraces, and increasing vegetation coverage, soil erosion remains a severe problem. Globally, soil erosion in 2012 amounted to an estimated 35.9 billion tons and exhibited an increase of 2.5% compared with 2001, driven primarily by cropland expansion [4]. Increased soil erosion affects sustainable agricultural development and is identified as a major environmental issue worldwide [4,5].
As the primary form of overland flow erosion, sheet erosion is a common phenomenon along hillslopes during rainfall events [6,7,8]. Sheet erosion is affected by many factors, such as vegetation/residue cover, soil texture and infiltration, rainfall intensity, slope length, and slope angle [6,9,10,11]. Unlike other factors, the relationship between slope angle and sheet erosion is not monotonic but increases at first and then decreases at higher slope angles [12,13]. As such, there is a unique slope angle, termed the critical slope angle (CSA), which induces the maximum amount of sheet erosion. The CSA can be attributed to the alteration in overland flow discharge induced by slope angle, which subsequently affects the shear force of overland flow [14,15]. The shear force of overland flow can regulate the transition from sheet erosion to rill erosion, which often leads to a greater erosion hazard [16,17,18]. Therefore, determining the CSA is crucial for investigating the progression of hillslope erosion, ranging from sheet erosion to rill erosion, and for accurately designing measures to control soil erosion.
Estimates of the CSA are typically derived from field observation, experimental simulation, and mathematical methods. Field observation and experimental simulation can provide accurate regional CSA information based on the local climate and soil conditions. These studies have demonstrated that CSA tends to increase with higher rainfall intensity and greater soil bulk density [19,20,21]. However, it is important to note that these studies typically require a considerable field workload over long time periods [22,23]. In contrast, mathematical methods estimate the CSA by analyzing the factors influencing soil erosion with fewer labor and time commitments [14,15]. Mathematical methods to estimate the CSA include soil erosion equation calculations and the shear stress derivation of overland flow. Soil erosion equations for obtaining the CSA are based on calculating soil loss at various slope angles. According to numerical simulation, it was observed that the CSA increased with rainfall duration in the same rainfall intensity [13]. Nevertheless, soil erosion equations do not express the physical mechanisms of CSA change, and several parameters need to be calibrated for their use in different regions. The mathematical method of shear stress analysis to obtain the CSA is based on deducing the shear stress change of overland flow on hillslopes with different slope angles [14,15]. This method assumes that, for a given hillslope, soil loss increases with the shear force of overland flow. The analysis of shear stress requires fewer parameters compared to the optimization of soil erosion equations and provides an explanation that elucidates the physical mechanism of CSA appearance [14,24]. However, previous studies on shear stress analysis assumed a significantly higher magnitude of rainfall intensity compared to soil infiltration intensity in their derivation.
In general, the method of shear stress analysis obtains the CSA more easily and quickly, and the physical mechanism is partially elucidated. However, previous studies on shear stress have had two main weaknesses. First, studies have not fully considered the influence of changes in rainfall intensity and soil infiltration on the CSA. This may produce appreciable differences compared to other methods. For example, for a loess soil, field observations and experimental simulations determined that the CSA varied within a range of 20–30° [12,22,25], a soil erosion equation calculated a CSA of ~25° [24], whereas the CSA estimated by shear stress analysis was 41.5–50° [15]. Second, the meaning of shear stress can also differ between practical studies. For example, rainfall intensity and soil infiltration rates change continually during rainfall events, while shear stress in past studies only reflected the instantaneous state rather than its accumulation during the whole rainfall event. Consequently, previous studies have failed to analyze the underlying physical mechanisms that elucidate the impacts of rainfall intensity, soil compaction, and rainfall on the CSA.
In this work, mathematical equations for estimating instantaneous and cumulative CSAs that consider rainfall intensity, soil infiltration, and different formulations of unit width discharge were derived on a theoretical basis. The major objectives were to (1) verify the efficacy of the mathematical equations for CSA calculation based on the results of previous studies; (2) model changes in instantaneous and cumulative CSAs according to mathematically derived equations; and (3) address the implications of the CSA for applications in soil erosion control. These results are expected to provide a stronger theoretical basis for predicting the CSA, with subsequent application in designing improved erosion control practices.

2. Data and Methods

2.1. Basic Formulas

When a homogeneous and isotropic hillslope is assumed, the equation for shear stress includes the Manning equation and 1D kinematic wave approximation to the St Venant equations. Correspondingly, the equation for shear stress can be written as [14,15,24]:
τ = ρ g n 0.6 IL cos α fL 0.6 sin 0.7 α ,  
where τ is shear stress of overland flow, ρ is density of overland flow, g is acceleration of gravity, n is the Manning coefficient, I is rainfall intensity, f is soil infiltration rate, L is the incline length of slope, α is slope angle, and (ILcosαfL) is unit width discharge.
The same horizontal projective length (SHPL) (Figure 1a) and same incline length (SIL) (Figure 1b) were used to change unit width discharge in this derivation. For any point A (A′) on the slope (Figure 1), the shear stress of overland flow at point A (A′) can be expressed as (Appendix A.1):
τ A = ρ g n 0.6 IX f X cos α 0.6 sin 0.7 α ,  
τ A = ρ g n 0.6 IX cos α fX 0.6 sin 0.7 α ,  
where τA is shear stress at point A, τA is shear stress at point A′, X is the horizontal projective length from A to the top of the slope (Figure 1a), X′ is the incline length from point A′ to the top of the slope (Figure 1b), and α is slope angle.
According to Equations (2) and (3), cumulative shear stress (the sum of shear stress at point A reflecting the soil erosion amount at the same point during the entire rainfall process) can be written as:
T A = t 1 t 2 ρ g n 0.6 X   0.6 I t f t cos α 0.6 sin 0.7 α dt ,
T A = t 1 t 2 ρ g n 0.6 X 0.6 I t cos α f t 0.6 sin 0.7 α dt ,  
where TA is cumulative shear stress at point A, TA is cumulative shear stress at point A′, I(t) is change in rainfall intensity with time, f(t) is change in soil infiltration rate with time, t1 is time of overland flow initiation, and t2 is time of overland flow cessation.

2.2. Mathematical Equation Derivations

2.2.1. Derivation of Instantaneous CSA Estimation Equation

According to Equation (2), τA and τAx (x > 0) have the same trend with a change in slope angle α (0° < α < 90°) because τA ≥ 0, hence the relationship between I and f can be expressed as:
I f   1 cos α   τ A     0 .
To facilitate the calculation of differential τA with respect to α, Equation (2) was rearranged as:
τ A 5 3 = ρ g 5 3 nX I f 1 cos α sin 0.7 α 5 3 .
Assuming a constant Manning coefficient (n), differentiating Equation (7) with respect to α leads to:
d τ A 5 3 d α =     1 6 ρ g 5 3 nX sin 1 6 α I 7 cos α f 6 + cos 2 α .          
By equating the left-hand side to zero, the expression becomes:
I f = 6 + cos 2 α 7 cos 3 α > 1 cos α   0 ° <   α <   90 ° .  
According to Equations (8) and (9), the τA5/3 term has an extremum (i.e., a maximum point). It is assumed that the arbitrary ratio of I/f is c (the value of c is greater than 1 because overland flow occurs only when rainfall intensity exceeds soil infiltration), and the corresponding slope angle (α) is b1 according to Equation (9). Accordingly, in Equation (8) the relationship between τA5/3/ and 0 is as follows:
d τ A 5 3 d α = 1 6 ρ g 5 3 n X f sin 1 6 α c   7 cos α 6 + cos 2 α > 0 ,         α <   b 1   d τ A 5 3 d α =   1 6 ρ g 5 3 n X f sin 1 6 α c   7 cos α 6 + cos 2 α < 0 ,         α > b 1   .  
Therefore, τA5/3 has a maximum value when the slope angle is b1. Because τA and τAx have the same trend, τA also has a maximum value when the slope angle equals b1.
As for the SIL, using the same method as above, by virtue of A′/dα = 0 (according to Equation (3)) it can be expressed as:
I f = 7 cos α 7     13 sin 2 α > 1 cos α   0 ° <   α   <   90 ° .

2.2.2. Derivation of Cumulative CSA Estimation Equation

According to Equations (4) and (5), cumulative shear stress is related to the rainfall and soil infiltration processes. Even though the process of rainfall is the same, the initiation and cessation times (t1 and t2) are different at different slope angles (α). Therefore, cumulative CSA should be analyzed according to the rainfall and soil infiltration processes.

2.3. Validation of Mathematical Equation Method

2.3.1. Validation by Field Observations

Field observations were used to validate the equations for instantaneous CSA (Equations (9) and (11)). The validation data came from the runoff plots of 6 natural rainfall events in He et al. [22]. The runoff plots (40°47′ N, 114°50′ E) were located in Zhangjiakou city, Hebei province (China). Equation (9) was used to estimate the CSA because the runoff plots had the same horizontal projective length. The rainfall and infiltration parameters in Equation (9) were the maximum 10 min intensity (I10) and soil stable infiltration rate [22].
A previous equation for the CSA [15] was used to contrast with Equation (9):
0.6 sin 2 α 0.7   cos 2 α sin 0.3 α cos 0.4 α ( N 0 sin α + cos α ) = d r s r r nL I f   0.6 ,
where α is the CSA, N0 is the friction coefficient (nondimensional), d is representative grain size of soil (mm), rs and r are soil dry bulk density and water density, respectively (N m−3), n is the Manning coefficient, L is slope incline length (m), I is rainfall intensity (mm−1), and f is soil infiltration rate (mm min−1). The characteristic parameters for the loess soil experimental trials obtained from Liu et al. [24] and Zhu et al. [26] are shown in Table 1.

2.3.2. Validation by Water Erosion Prediction Project (WEPP) Model Simulations

WEPP model results were used to validate the estimated cumulative CSA calculated by the mathematical equations (Equation (4)). The WEPP model results came from Wu et al. [13]. The simulation validation was conducted on a loam soil from Montana (USA) [13,27]. The characteristic parameters for the loam soil are shown in Table 2. Equation (4) was used to estimate the cumulative CSA because the simulation results from the WEPP model were based on the same horizontal projective length. The WEPP model has been demonstrated to accurately model hillslope erosion processes and it has been well tested by many studies (e.g., Cochrane and Flanagan [28]; Abaci and Papanicolaou [29]). A modified Green–Ampt (GA) slope infiltration model [30] was used to simulate soil infiltration processes in the WEPP model. Since only soil horizontal surface infiltration rate was required in Equation (4), the GA model [31] was used to analyze changes in soil infiltration with time, which was different from Wu et al. [13]. The GA model can be expressed as:
f = k 1 + SM / h p   ,
h p = k t + SM   ln 1 + h p / SM ,  
where f is soil infiltration rate (mm min−1), k is effective saturated hydraulic conductivity (mm min−1), M is effective porosity, S is the wetting front soil suction head (m), and hp is cumulative infiltration depth (mm).

2.4. Validation of Mathematical Equation Method

Ratios of rainfall intensity to soil infiltration from 0 to 10 were used to model the change in instantaneous CSA based on Equations (9) and (11). In contrast to instantaneous CSA, cumulative CSA is influenced by different rainfall patterns. Thus, three rainfall patterns were utilized in the analysis of cumulative CSA based on Equations (4) and (5). According to the relationship between rainfall intensity and soil infiltration, the three rainfall patterns were divided into light, normal and heavy rainfall regimes. Light rain was defined as rainfall intensities lower than the soil infiltration rate, thereby resulting in no overland flow. As a result, only the normal and heavy rainfall regimes were analyzed with respect to runoff generation. The total amount of rainfall was assumed to be the same for both the normal and heavy rainfall regimes. Equations for rainfall intensity and soil horizontal surface infiltration rate with time were assumed as:
I t = 6 t + 1 1 , ( h e a v y   r a i n ) I t = 0.276 t 2 + 1.38 t , ( n o r m a l   r a i n )                             f t = 0.5 t + 0.75 , ( s o i l   i n f i l t r a t i o n )                                 .
where t is time (nondimensional), I(t) is variation in rainfall intensity with time for heavy rainfall regime, I′(t) is variation in rainfall intensity with time for normal rainfall regime, and f(t) is variation in soil infiltration with time.

3. Results and Discussion

3.1. Validation of Mathematical Equations

3.1.1. Comparison to Field Plot Observations

The instantaneous CSAs estimated by Equation (9) displayed a good fit with those observed from field runoff plots (NSE = −1.01, RMSE = 5.28) of loess soils from Zhangjiakou city (China). The estimated CSA value obtained from Equation (9) surpassed that obtained from Equation (14) (RMSE = 25.27). Previous theoretical studies did not fully consider soil infiltration as a factor [14,15], resulting in rainfall intensity having no effect on CSAs (Figure 2a). The maximum discrepancy between the estimated values obtained from Equation (9) and the observed CSAs was 10° (I10 = 10.9 mm min−1), while the minimum discrepancy was 2° (I10 = 7.5 mm min−1). The variability in soil moisture content prior to rainfall leads to inconsistent actual soil infiltration rates during each rainfall event. The CSAs estimated by Equation (9) employed a stable infiltration rate instead of the actual infiltration rate, which could potentially account for the observed discrepancy between the estimated and observed CSAs. Overall, Equation (9) effectively estimated the CSA when a more accurate soil infiltration rate was utilized.

3.1.2. Comparison to WEPP Simulation Results

Cumulative CSAs estimated by Equation (4) were in good agreement with simulation results from the WEPP model (NSE = 0.96, RMSE = 1.73) for a loam soil in Montana (USA) (Figure 3b). The estimated and simulated cumulative CSAs displayed the same changes in trends (divided into three periods) and varied within a range of about 0–33.6° and 0–32.5°, respectively, with a maximum difference of only 2.6° (t = 4) (Figure 3a). The small differences may have been caused by the soil infiltration calculation and/or soil erodibility. Compared to the GA model (used in Equation (4)), the modified GA model (used in the WEPP) considered the influence of slope angle on soil infiltration [30]. However, there was little difference between the GA model and the modified GA model when rainfall intensity was small (i.e., the ratio of rainfall intensity to soil infiltration rate was less than 2) [30]. On the other hand, the simulated CSA calculated by the WEPP model considered soil erodibility, which is affected by slope angle [32], whereas Equation (4) only considered the shear stress of overland flow. Equation (4) needed two parameters (rainfall intensity and soil infiltration rate) to estimate the CSA, whereas the WEPP model required seven parameters. Thus, the mathematical equations derived in this study were simpler and more efficient to apply.

3.2. Characteristics of Instantaneous and Cumulative CSAs

3.2.1. Modeling the Change in Instantaneous CSA

The modeled instantaneous CSA increases with increasing the I/f ratio from 0 to 10 based on Equations (9) and (11) (Figure 4). This indicates that the instantaneous CSA increases with rainfall intensity in the same soil, or conversely that a soil with a high infiltration capacity has a lower instantaneous CSA in a given rainfall intensity. In physical terms, the I/f ratio denotes the net unit width discharge over a horizontal surface. Therefore, Equations (9) and (11) delineate the slope angle at which the maximum shear stress occurs on a hillslope under a certain net unit width discharge. The shear stress of overland flow increases with increases in slope angle and overland flow depth (Figure 1). However, the overland flow depth and slope angle have an inverse relationship [14,15]. Therefore, shear stress has a maximum value at a specific slope angle, which means that the CSA is evident. As the I/f ratio increases, the effect of the slope on overland flow depth decreases over a range of slope angles (Figure 5a). This leads to an increase in the instantaneous CSA as the I/f ratio increases (Figure 5b).
There were some differences in determining instantaneous CSAs based on the same horizontal projective length (SHPL) versus the same incline length (SIL). Instantaneous CSAs based on SHPL were greater than those based on SIL. The difference can be divided into two components. When the I/f ratio was less than 2.1, the difference was less than 5°, but then the difference increased with increasing I/f (Figure 4). On the other hand, the limit values for the two instantaneous CSAs were not the same. The limit value based on SHPL was 90°, whereas the limit value based on SIL was 47.5°. This is due to differences in overland flow depth caused by the unit width discharge characterization. Under SHPL, an increase in infiltration led to a decrease in unit width discharge as the slope angle increased (Figure 1a). In contrast, under SIL, the decrease in unit width discharge was due to a decrease in rainfall (Figure 1b). When the ratio was small (i.e., less than 2.1), the overland flow depth (determined by unit width discharge) exhibited a similar change under both conditions, resulting in a similar CSA (Figure 5). However, as rainfall had a greater effect on overland flow depth than soil infiltration, overland flow depth varied less under SHPL than SIL for a range of slopes when the I/f ratio varied. This led to a wider difference in shear stress between the two conditions, so that the discrepancy in the instantaneous CSA between the two conditions increased as the I/f ratio increased (Figure 5). These characteristics imply that, in an area with frequent rainstorms, control measures for soil erosion utilizing instantaneous CSAs should specifically state whether they are based on SHPL or SIL conditions.

3.2.2. Modeling the Change in Cumulative CSA

Substituting Equation (13) into Equation (4) yielded the modeled cumulative CSA under normal and heavy rainfall regimes. The modeled cumulative CSA increased with an increasing I/f ratio but did not continuously decrease as rainfall intensity decreased (Figure 6). The normal rain pattern was divided into three phases (Figure 6a). At t1, overland flow did not occur, and the cumulative CSA was 0°. With increasing rainfall progression in the t2 phase, the cumulative CSA increased from 0° to 33°. The cumulative CSA increased quickly with an increasing I/f ratio at the beginning of this phase (t2), before leveling off as the I/f ratio decreased at the end of the phase. During the t3 phase, there was no runoff and hence the cumulative CSA did not change.
In contrast to the normal rainfall regime, the heavy rainfall pattern had only two phases (Figure 6b). A change in the cumulative CSA only occurred in the t4 phase, in which it decreased from 47° to 35° as the I/f ratio decreased. Similar to the normal rainfall regime, the cumulative CSA was a constant angle in the t5 phase. These characteristics are attributed to the mechanism regulating the change in instantaneous CSA with a changing I/f ratio. When the I/f ratio increases, the slope angle with the maximum erosion amount increases, leading to a cumulative CSA increase. However, the cumulative erosion amount did not decrease with increasing rainfall intensity for the hillslope during the rainfall event, hence the cumulative CSA was at a constant value when rainfall intensity was less than the soil infiltration rate.
The heavy rainfall regime had a larger observed CSA (35°) (cumulative CSA after rainfall) than the normal rainfall regime (31°) (Figure 6). This indicates that the same amount of rainfall and soil infiltration could result in different observed CSA values, with the heavy rainfall regime having a larger CSA. Although rainfall events usually induce various types of soil erosion (e.g., splash, sheet, rill, gully) under natural conditions, the amount of sheet erosion in a heavy rainfall regime tends to be larger than that in a normal rainfall regime for the same slope [13]. Therefore, a significantly higher observed CSA implies a larger amount of sheet erosion for the same area.
The modeled characteristics of instantaneous CSA (Figure 4) and cumulative CSA (Figure 3a) across a range of I/f values demonstrate the effects of rainfall intensity [19,20], soil bulk density [21], and rainfall duration [13] on the CSA. For example, rainfall intensity affects the I/f ratio directly (Figure 4), whereas rainfall duration and soil bulk density indirectly affect the ratio by changing the soil infiltration rate. The soil infiltration rate decreases with increasing rainfall duration, resulting in an increasing I/f under steady rainfall conditions (Figure 3a). Moreover, soil compaction can reduce the soil infiltration rate [33,34], which increases the I/f ratio under similar rainfall conditions (Figure 4). Therefore, the CSA increases with rainfall duration and soil compaction. As a result, soil erosion control strategies should be made in accordance with the CSA, such as increasing soil infiltration (Figure 4) by increasing vegetation cover, improving soil structure, or converting steep cropland by terracing it [35,36,37].

3.2.3. Modeling Instantaneous CSA and Cumulative CSA in Different Rainfall Conditions

Substituting Equation (15) into Equation (9) yielded the modeled instantaneous CSA under normal and heavy rainfall regimes. The modeled instantaneous CSA varied with the I/f ratio, and the maximum instantaneous CSA occurred at the peak of the ratio during the rainfall process (Figure 6). Specifically, the maximum instantaneous CSA was 36° during normal rainfall, while it reached 47° in heavy rainfall. Notably, the cumulative CSA after rainfall exhibited similar values (31° and 35°) in the two different rainfall patterns. This indicates that the disparity between the two commonly used CSAs (maximum instantaneous and cumulative CSA after rainfall) is more pronounced in heavy rainfall events. Moreover, even if the observed CSAs exhibited similarity (with a difference of only 4°) in different rainfall events, there may still have been significant variations in the maximum instantaneous CSAs (with a difference of up to 11°) due to disparities in rainfall processes. Therefore, the differentiation of erosion mechanisms among various rainfall occurrences can be more effectively elucidated by instantaneous CSA rather than cumulative CSA.
Therefore, soil erosion control measures should be developed based on the maximum instantaneous CSA instead of the CSA obtained from field observations or experimental simulations, especially for heavy rainfall events. Given its ability to better characterize rainfall intensity, especially during heavy rain, and its associated erosion potential, I5 (maximum 5 min rainfall intensity) is recommended as a more conservative value to estimate the maximum instantaneous CSA in practice.

3.3. Effect of Manning Coefficient (n) on the CSA

The Manning coefficient, as a crucial parameter in the calculation of shear stress (Equation (1)), exerts a significant influence on the CSA. The Manning coefficient (n) is significantly influenced by various factors affecting sheet flow, including but not limited to rainfall intensity, slope gradient, slope surface roughness, sediment size, and vegetation [38,39,40]. According to the assumption stated for Equation (1), n only required the consideration of rainfall intensity, slope gradient, and surface roughness in this study. These factors exhibit a direct proportionality to the Manning coefficient [40,41,42].
Rainfall intensity (40 mm/h–120 mm/h) exhibited an approximate 10% increase in n [40], while the slope gradient (0–20°) demonstrated a roughly 5% increment [41]. Therefore, the actual field CSA for this rainfall event exceeded the value calculated by the mathematical equations in this study, as rainfall intensity and slope gradient increased. However, the effect of rainfall duration on soil surface roughness remains uncertain. Several studies have suggested that soil surface roughness may either increase [41] or decrease [43,44] with rainfall duration. Additionally, it had been posited that soil erosion could lead to an increase in the surface roughness of initially smooth surfaces, while decreasing the roughness of initially rough surfaces [45]. These also remains an uncertainty regarding the variability of the Manning coefficient caused by roughness during rainfall. The CSA calculated through mathematical equations in this study may potentially have underestimated its actual value, with a potential underestimation of approximately 5° (I/f < 2.5, slope gradient < 50°). However, this underestimation did not affect the result that the CSA increased with the I/f ratio.

3.4. Improvements of the Mathematical Equations

The mathematical equations developed in this study improved the estimation accuracy of the CSA by fully considering rainfall intensity and soil infiltration. The equations also illustrate the underlying physical mechanisms by which the CSA is regulated by rainfall intensity and soil infiltration through changing the shear stress of overland flow. Compared with previous mathematical equations, the newly developed mathematical equations required only two parameters, making them easier and more convenient to use.
According to the characteristics of CSAs, soil erosion control measures should be based on the maximum instantaneous CSA instead of the CSA observed from field observations after rainfall. Meanwhile, the basic conditions for the same horizontal projective length and same incline length should be demonstrated for areas with frequent rainstorms when making land use policies according to CSAs, such as designating the utilization of sloping land.
Although mathematical equations can provide reliable estimates of CSAs, vegetation and soil detachment resistance were not considered in this study. On the one hand, as an essential landscape element, vegetation can directly influence soil erosion through increasing soil infiltration rates, decreasing soil erodibility, and increasing flow resistance [13,46]. Moreover, vegetation may change the rainfall–runoff dynamics [36] and soil erodibility may change during a rainfall event. Soil detachment resistance decreases as rainfall intensity and soil moisture content increase [47,48]. Furthermore, soil detachment resistance decreases as the slope angle increases [15,24]. Therefore, vegetation and soil detachment resistance may appreciably influence soil erosion dynamics at various scales and are highlighted for further investigation.

3.5. Limitations

The influence of steep slopes (>50°) on soil infiltration and gravity introduces significant uncertainty to the calculation of CSAs [13,49]. Therefore, when applied to steep slopes, the mathematical equations derived from this study should be cautiously utilized and solely used as a reference in real conditions. Simultaneously, the derived mathematical equations do not account for the impact of soil compaction induced by heavy rainfall on CSAs. As the derived equations solely consider the shear force of overland flow and do not take into account the influence of other factors (e.g., the impact of rainfall on soil structure), more field and experimental data are needed to validate the mathematical equations across a much wider range of soil conditions.

4. Conclusions

This study derived a mathematical method for estimating the instantaneous and cumulative CSAs that included the influence of rainfall intensity and soil infiltration. The newly derived mathematical equations predicted the CSA more efficiently than those obtained from natural observations and erosion model simulations. Furthermore, these equations provide valuable insights into the influence of various factors (e.g., soil bulk density, rainfall intensity and duration) on the change in the CSA. The efficacy of the mathematical equations was well evaluated by comparisons to natural runoff–erosion observations and soil erosion model simulations. According to the derived equations, the instantaneous CSA increased as the rainfall intensity to soil infiltration (I/f) ratio increased. Moreover, heavy rain had a higher cumulative CSA than normal rainfall events. Factors influencing the CSA (e.g., soil bulk density, rainfall intensity and duration) through their effects on changing the I/f ratio were well characterized by the mathematical equations. The instantaneous CSA can be used to characterize soil erosion dynamics during rainfall–runoff events. Most notably, soil erosion control strategies should be introduced in accordance with the maximum instantaneous CSA instead of the CSA observed after rainfall, especially for heavy rainfall events. The actual values of CSAs calculated using mathematical equations in this study may be underestimated due to the variability of the Manning coefficient (n). Future studies should include the influence of other soil erosion factors (e.g., vegetation, soil detachment resistance, etc.) on the CSA to better understand soil erosion processes and design more effective erosion control practices.

Author Contributions

Data curation, Y.W.; Investigation, Y.P.; Validation, Z.P.; Writing—original draft, M.W.; Writing—review and editing, D.C. All authors have read and agreed to the published version of the manuscript.

Funding

This work was supported by the National Key Research and Development Program of China (2021YFD1700802) and the National Natural Science Foundation of China (42107393 and 42177352).

Data Availability Statement

The data presented in this study are available in insert article.

Acknowledgments

We thank Randy A. Dahlgren (University of California–Davis) for his preliminary review of the manuscript.

Conflicts of Interest

The authors declare no conflict of interest.

Notation

The following symbols are used in this paper:
τshear stress of overland flow
ρdensity of overland flow
gacceleration of gravity
nManning coefficient
Irainfall intensity
fsoil infiltration rate
Lincline length of slope
αslope angle
τAshear stress at point A (Figure 1a)
τ′Ashear stress at point A′ (Figure 1b)
X horizontal projective length from A to the top of the slope (Figure 1a)
X′incline length from A′ point to the top of the slope (Figure 1b)
TAcumulative shear stress at point A during rainfall (Figure 1a)
T′A cumulative shear stress at point A′ during rainfall (Figure 1b)
I(t)change in rainfall intensity with time during rainfall
f(t)change in soil infiltration rate with time during rainfall
t1time of overland flow initiation
t2time of overland flow cessation.
N0friction coefficient (nondimensional)
dgrain size of soil (mm)
rssoil dry bulk density (N m−3)
rwater density (N m−3)
keffective saturated hydraulic conductivity (mm min−1)
Meffective porosity
Swetting front soil suction head (m)
hpcumulative infiltration depth (mm)
kt transport coefficient (WEPP model)
τcadjusted soil critical shear stress (WEPP model)
fmeanmean of effective soil hydraulic conductivity ((mm h−1))

Appendix A

Appendix A.1. Derivation of the Shear Force Equation for Overland Flow

The Manning equation and the unit width discharge formula were used to analyze the shear stress of overland flow in this study. The equations are well applied in runoff and soil erosion, especially in bare land [14,46]. For any point A on the slope (Figure A1), the shear stress of overland flow at point A is:
τ A   = ρ g h A sin α ,
where τA is shear stress, ρ is density of overland flow, g is acceleration of gravity, hA is the depth of overland flow, and α is slope angle. The unit width discharge at point A is:
q A =   v A h A ,
where qA is unit width discharge and vA is overland flow velocity. The relationship between vA and hA can be expressed by the Manning equation as:
v A = 1 n h A 2 / 3 S 1 / 2 ,
where n is the Manning coefficient and S is surface slope angle. According to the 1D kinematic wave model approximation, S is expressed as sinα [15,50], so vA can be expressed as:
v A   = 1 n h A 2 / 3 sin α 1 / 2 .  
Substituting Equations (2) and (4) and into Equation (1), τA can be expressed as:
h A = n 0.6 q A 0.6 sin 0.3 α ,
τ A = ρ g n 0.6 q A 0.6 sin 0.7 α .  
Figure A1. Shear stress of overland flow on hillslope during rainfall.
Figure A1. Shear stress of overland flow on hillslope during rainfall.
Water 15 03341 g0a1
For any point A on the slope (Figure 1a), it is assumed that the horizontal projective length from this point to the top of the slope is X. When the slope angle (α) changes, from point A to the top of the slope, the total rainfall is the same but the total infiltration amount changes, and the unit width discharge in point A can be expressed as:
q A = 0 X dq dx dx = IX f X cos α ,
where qA is unit width discharge at point A, I is rainfall intensity, and f is soil infiltration rate.
It is assumed that incline length from point A′ to the top of the slope is X′ (Figure 1b). When the slope angle (α) changes, the unit width discharge at point A can be expressed as:
q A = 0 X dq dx dx = IX cos α fX .
where qA is unit width discharge at point A′.
Substituting Equations (7) and (8) into Equation (6) yields the instantaneous shear stress τA and τ’A.

References

  1. Mandal, D.; Giri, N.; Srivastava, P. The magnitude of erosion-induced carbon (C) flux and C-sequestration potential of eroded lands in India. Eur. J. Soil Sci. 2020, 71, 151–168. [Google Scholar] [CrossRef]
  2. Pimentel, D. Soil erosion: A food and environmental threat. Environ. Dev. Sustain. 2006, 8, 119–137. [Google Scholar] [CrossRef]
  3. von Ruette, J.; Lehmann, P.; Or, D. Effects of rainfall spatial variability and intermittency on shallow landslide triggering patterns at a catchment scale. Water Resour. Res. 2014, 50, 7780–7799. [Google Scholar] [CrossRef]
  4. Borrelli, P.; Robinson, D.A.; Fleischer, L.R.; Lugato, E.; Ballabio, C.; Alewell, C.; Meusburger, K.; Modugno, S.; Schuett, B.; Ferro, V.; et al. An assessment of the global impact of 21st century land use change on soil erosion. Nat. Commun. 2017, 8, 2013. [Google Scholar] [CrossRef] [PubMed]
  5. Montgomery, D.R. Soil erosion and agricultural sustainability. Proc. Natl. Acad. Sci. USA 2007, 104, 13268–13272. [Google Scholar] [CrossRef]
  6. Huang, C.H. Empirical analysis of slope and runoff for sediment delivery from interrill areas. Soil Sci. Soc. Am. J. 1996, 60, 1279–1280. [Google Scholar] [CrossRef]
  7. Parsons, A.J. How reliable are our methods for estimating soil erosion by water? Sci. Total Environ. 2019, 676, 215–221. [Google Scholar] [CrossRef]
  8. Zhang, Q.W.; Wang, Z.L.; Wu, B.; Shen, N.; Liu, J.E. Identifying sediment transport capacity of raindrop-impacted overland flow within transport-limited system of interrill erosion processes on steep loess hillslopes of China. Soil Tillage Res. 2018, 184, 109–117. [Google Scholar] [CrossRef]
  9. Jouquet, P.; Janeau, J.L.; Pisano, A.; Hai Tran, S.; Orange, D.; Luu Thi Nguyet, M.; Valentin, C. Influence of earthworms and termites on runoff and erosion in a tropical steep slope fallow in Vietnam: A rainfall simulation experiment. Appl. Soil Ecol. 2012, 61, 161–168. [Google Scholar] [CrossRef]
  10. Zhang, X.; Yu, G.Q.; Li, Z.B.; Li, P. Experimental study on slope runoff, erosion and sediment under different vegetation types. Water Resour. Manag. 2014, 28, 2415–2433. [Google Scholar] [CrossRef]
  11. Wang, H.; Zhang, G.H. Temporal variation in soil erodibility indices for five typical land use types on the Loess Plateau of China. Geoderma 2021, 381, 114695. [Google Scholar] [CrossRef]
  12. Cheng, Q.J.; Ma, W.J.; Cai, Q.G. The relative importance of soil crust and slope angle in runoff and soil loss: A case study in the hilly areas of the Loess Plateau, North China. GeoJournal 2008, 71, 117–125. [Google Scholar] [CrossRef]
  13. Wu, S.B.; Yu, M.H.; Chen, L. Nonmonotonic and spatial-temporal dynamic slope effects on soil erosion during rainfall-runoff processes. Water Resour. Res. 2017, 53, 1369–1389. [Google Scholar] [CrossRef]
  14. Horton, R.E. Erosional development of streams and their drainage basins-hydro physical approach to quantitative morphology. Geol. Soc. Am. Bull. 1945, 56, 275–370. [Google Scholar] [CrossRef]
  15. Liu, Q.Q.; Chen, L.; Li, J.C. Influences of slope gradient on soil erosion. Appl. Math. Mech. 2001, 22, 510–519. [Google Scholar] [CrossRef]
  16. Giménez, R.; Govers, G. Flow detachment by concentrated flow on smooth and irregular beds. Soil Sci. Soc. Am. J. 2002, 66, 1475–1483. [Google Scholar] [CrossRef]
  17. Govers, G.; Giménez, R.; Govers, R.; Van Oost, K. Rill erosion: Exploring the relationship between experiments, modelling and field observations. Earth-Sci. Rev. 2007, 84, 87–102. [Google Scholar] [CrossRef]
  18. Liu, G.; Zheng, F.L.; Wilson, G.V.; Xu, X.M.; Liu, C. Three decades of ephemeral gully erosion studies. Soil Tillage Res. 2021, 212, 105046. [Google Scholar] [CrossRef]
  19. Liao, Y.S.; Cai, Q.G.; Cheng, Q.J. Critical topographic condition for slope erosion in hilly-gully region of Loess Plateau. Soil Water Conserv. 2008, 2, 32–38. [Google Scholar] [CrossRef]
  20. Chen, X.A.; Cai, Q.G.; Zhang, L.C.; Qi, J.Y.; Zheng, M.G.; Nie, B.B. Research on critical slope of soil erosion in a Hilly Loess Region on the Loess Plateau. J. Mt. Sci. 2010, 28, 415–421. [Google Scholar] [CrossRef]
  21. Foster, R.L.; Martin, G.L. Effect of unit weight and slope on erosion. J. Irrig. Drain. Eng. 1969, 95, 551–561. [Google Scholar] [CrossRef]
  22. He, J.J.; Cai, G.Q.; Liu, S.B. Effects of slope gradient on slope runoff and sediment yield under different single rainfall conditions. J. Appl. Ecol. 2012, 23, 1263–1268. [Google Scholar] [CrossRef]
  23. Fu, S.H.; Liu, B.Y.; Liu, H.P.; Xu, L. The effect of slope on interrill erosion at short slopes. Catena 2011, 84, 29–34. [Google Scholar] [CrossRef]
  24. Liu, Q.Q.; Xiang, H.; Singh, V.P. A simulation model for unified interrill erosion and rill erosion on hillslopes. Hydrol. Process. 2006, 20, 469–486. [Google Scholar] [CrossRef]
  25. Chen, Y.Z. A preliminary analysis of the processes of sediment yield in small catchment on the Loess Plateau. Geogr. Res. 1983, 2, 35–47. [Google Scholar] [CrossRef]
  26. Zhu, P.; Zhang, G.; Wang, H.; Xing, S. Soil infiltration properties affected by typical plant communities on steep gully slopes on the Loess Plateau of China. J. Hydrol. 2020, 590, 125535. [Google Scholar] [CrossRef]
  27. Flanagan, D.; Nearing, M. USDA-Water Erosion Prediction Project: Hillslope Profile and Watershed Model Documentation; NSERL Report No. 10; National Soil Erosion Research Laboratory, USDA-Agricultural Research Service: West Lafayette, IN, USA, 1995.
  28. Cochrane, T.A.; Flanagan, D.C. Representative hillslope methods for applying the WEPP model with DEMS and GIS. Trans. ASAE 2003, 46, 1041–1049. [Google Scholar] [CrossRef]
  29. Abaci, O.; Papanicolaou, A.N.T. Long-term effects of management practices on water-driven soil erosion in an intense agricultural sub-watershed: Monitoring and modelling. Hydrol. Process. 2009, 23, 2818–2837. [Google Scholar] [CrossRef]
  30. Chen, L.; Young, M.H. Green-Ampt infiltration model for sloping surfaces. Water Resour. Res. 2006, 42, W07420. [Google Scholar] [CrossRef]
  31. Green, W.H.; Ampt, G.A. Studies on soil physics. J. Agric. Sci. 1911, 4, 1–24. [Google Scholar] [CrossRef]
  32. Ascough, J.C.; Baffaut, C.; Nearing, M.A.; Liu, B.Y. The WEPP watershed model: I. Hydrology and erosion. Trans. ASAE 1997, 40, 921–933. [Google Scholar] [CrossRef]
  33. Alaoui, A.; Helbling, A. Evaluation of soil compaction using hydrodynamic water content variation: Comparison between compacted and non-compacted soil. Geoderma 2006, 134, 97–108. [Google Scholar] [CrossRef]
  34. Fullen, M. Compaction, hydrological processes and soil erosion on loamy sands in east Shropshire, England. Soil Tillage Res. 1985, 6, 17–29. [Google Scholar] [CrossRef]
  35. Yu, M.; Zhang, L.; Xu, X.; Feger, K.-H.; Wang, Y.; Liu, W.; Schwärzel, K. Impact of land-use changes on soil hydraulic properties of Calcaric Regosols on the Loess Plateau, NW China. J. Soil Sci. Plant. Nut. 2015, 178, 486–498. [Google Scholar] [CrossRef]
  36. Chen, L.; Sela, S.; Svoray, T.; Assouline, S. Scale dependence of Hortonian rainfall-runoff processes in a semiarid environment. Water Resour. Res. 2016, 52, 5149–5166. [Google Scholar] [CrossRef]
  37. Cui, Z.; Wu, G.L.; Huang, Z.; Liu, Y. Fine roots determine soil infiltration potential than soil water content in semiarid grassland soils. J. Hydrol. 2019, 578, 124023. [Google Scholar] [CrossRef]
  38. Zhang, G.H.; Luo, R.T.; Cao, Y.; Shen, R.C.; Zhang, X.C. Impacts of sediment load on Manning coefficient in supercritical shallow flow on steep slopes. Hydrol. Process. 2010, 24, 3909–3914. [Google Scholar] [CrossRef]
  39. Vastila, K.; Jarvela, J. Characterizing natural riparian vegetation for modeling of flow and suspended sediment transport. J. Soil. Sediment. 2018, 18, 3114–3130. [Google Scholar] [CrossRef]
  40. Shen, E.; Liu, G.; Dan, C.; Chen, X.; Ye, S.; Li, R.; Li, H.; Zhang, Q.; Zhang, Y.; Guo, Z. Estimating Manning’s coefficient n for sheet flow during rainstorms. Catena 2023, 226, 107093. [Google Scholar] [CrossRef]
  41. Ding, W.F.; Li, Y.L.; Wang, Y.F.; Cheng, D.B.; Zhang, P.C.J. Study on runoff hydrodynamics of purple soil slope under the rainfall simulation experiment. Soil Water Conserv. 2010, 24, 66–69. [Google Scholar] [CrossRef]
  42. Giménez, R.; Govers, G. Interaction between bed roughness and flowhydraulics in eroding rills. Water Resour. Res. 2001, 37, 791–799. [Google Scholar] [CrossRef]
  43. Zhao, L.; Liang, X.; Wu, F. Soil surface roughness change and its effect on runoff and erosion on the Loess Plateau of China. J. Arid Land 2014, 6, 400–409. [Google Scholar] [CrossRef]
  44. Zheng, Z.C.; He, S.Q.; Wu, F.Q. Changes of soil surface roughness under water erosion process. Hydrol. Process. 2014, 28, 3919–3929. [Google Scholar] [CrossRef]
  45. Bahddou, S.; Otten, W.; Whalley, W.R.; Shin, H.C.; El Gharous, M.; Rickson, R.J. Changes in soil surface properties under simulated rainfall and the effect of surface roughness on runoff, infiltration and soil loss. Geoderma 2023, 431, 116341. [Google Scholar] [CrossRef]
  46. Fu, S.; Mu, H.; Liu, B.; Yu, X.; Liu, Y. Effect of plant basal cover on velocity of shallow overland flow. J. Hydrol. 2019, 577. [Google Scholar] [CrossRef]
  47. Fox, D.M.; Bryan, R.B. The relationship of soil loss by interrill erosion to slope gradient. Catena 2000, 38, 211–222. [Google Scholar] [CrossRef]
  48. Wang, C.F.; Wang, B.; Wang, Y.Q.; Wang, Y.J.; Zhang, W.L.; Yan, Y.K. Impact of near-surface hydraulic gradient on the interrill erosion process. Eur. J. Soil Sci. 2020, 71, 598–614. [Google Scholar] [CrossRef]
  49. Gabet, E.J.; Dunne, T. A stochastic sediment delivery model for a steep Mediterranean landscape. Water Resour. Res. 2003, 39, 1237. [Google Scholar] [CrossRef]
  50. Woolhiser, D.A.; Liggett, A. Unsteady, one-dimensional flow over a plane—The rising hydrograph. Water Resour. Res. 1967, 3, 753–771. [Google Scholar] [CrossRef]
Figure 1. Two derivations for unit width discharge change: (a) same horizontal projective length and (b) same incline length.
Figure 1. Two derivations for unit width discharge change: (a) same horizontal projective length and (b) same incline length.
Water 15 03341 g001
Figure 2. CSAs (critical slope angles) determined from field observations and mathematical methods: (a) comparison between CSAs estimated by Equations (9) and (14) and field observations; (b) direct comparison between Equation (9) and observed values.
Figure 2. CSAs (critical slope angles) determined from field observations and mathematical methods: (a) comparison between CSAs estimated by Equations (9) and (14) and field observations; (b) direct comparison between Equation (9) and observed values.
Water 15 03341 g002
Figure 3. CSAs (critical slope angles) calculated by different mathematical equations compared to WEPP (Water Erosion Prediction Project) model: (a) comparison of estimated CSAs between Equation (4) and WEPP model; (b) direct comparison between Equation (4) and WEPP model.
Figure 3. CSAs (critical slope angles) calculated by different mathematical equations compared to WEPP (Water Erosion Prediction Project) model: (a) comparison of estimated CSAs between Equation (4) and WEPP model; (b) direct comparison between Equation (4) and WEPP model.
Water 15 03341 g003
Figure 4. Relationship between ratio of rainfall intensity and soil infiltration rate (I/f) and instantaneous CSA (critical slope angle).
Figure 4. Relationship between ratio of rainfall intensity and soil infiltration rate (I/f) and instantaneous CSA (critical slope angle).
Water 15 03341 g004
Figure 5. Changes in overland flow depth and shear stress with slope angle at different I/f ratios: (a) changes in overland flow depth; (b) changes in shear stress (Parameters ρ, g, n, f and L are assumed to be 1 for this analysis).
Figure 5. Changes in overland flow depth and shear stress with slope angle at different I/f ratios: (a) changes in overland flow depth; (b) changes in shear stress (Parameters ρ, g, n, f and L are assumed to be 1 for this analysis).
Water 15 03341 g005
Figure 6. Changes in CSA (critical slope angle) under different rainfall regimes: (a) normal rainfall; (b) heavy rainfall (soil infiltration curves in t3 and t5 indicate restoration of infiltration capacity).
Figure 6. Changes in CSA (critical slope angle) under different rainfall regimes: (a) normal rainfall; (b) heavy rainfall (soil infiltration curves in t3 and t5 indicate restoration of infiltration capacity).
Water 15 03341 g006
Table 1. Rainfall and soil parameters used for the CSA calculation.
Table 1. Rainfall and soil parameters used for the CSA calculation.
ParameterValueParameterValue
I (I10, mm min−1)7.5, 8.9, 7.9, 10.9, 8.2, 8.2rs(N m−3)1300
f, k (mm min−1)0.625r (N m−3)1100
d (mm)0.05L (m)20
N00.047n0.03
M0.45S (mm)60
Table 2. Parameters used in the WEPP (Water Erosion Prediction Project) model and Equation (4).
Table 2. Parameters used in the WEPP (Water Erosion Prediction Project) model and Equation (4).
ParameterValueParameterValue
I (mm h−1)45Kt (m0.5 s2 kg−0.5)0.029
fmean (mm h−1)22.5τc (Pa)0.025
SM (m)0.0173Slope length (m)22.1
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Wang, M.; Chen, D.; Wang, Y.; Pan, Z.; Pan, Y. A Mathematical Method for Estimating the Critical Slope Angle of Sheet Erosion. Water 2023, 15, 3341. https://doi.org/10.3390/w15193341

AMA Style

Wang M, Chen D, Wang Y, Pan Z, Pan Y. A Mathematical Method for Estimating the Critical Slope Angle of Sheet Erosion. Water. 2023; 15(19):3341. https://doi.org/10.3390/w15193341

Chicago/Turabian Style

Wang, Mingfeng, Dingjiang Chen, Yucang Wang, Zheqi Pan, and Yi Pan. 2023. "A Mathematical Method for Estimating the Critical Slope Angle of Sheet Erosion" Water 15, no. 19: 3341. https://doi.org/10.3390/w15193341

APA Style

Wang, M., Chen, D., Wang, Y., Pan, Z., & Pan, Y. (2023). A Mathematical Method for Estimating the Critical Slope Angle of Sheet Erosion. Water, 15(19), 3341. https://doi.org/10.3390/w15193341

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop