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

: 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 inﬁltration to efﬁciently 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 inﬁltration ( 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 ﬁeld observations or experimental simulations. The mathematical equations illustrate the underlying physical mechanisms by which rainfall intensity and soil inﬁltration affect the CSA through changing the shear stress of overland ﬂow. The results of this study provide critical information for guiding the development of effective soil erosion control strategies.


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.

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]: 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): 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.
Water 2023, 15, x FOR PEER REVIEW 3 of 16

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]: 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): 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 = ρgn 0.6 X 0.6 I  cosα f t 0.6 sin 0.7 α 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 = According to Equation (2), τ A and τ A x (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: To facilitate the calculation of differential τ A with respect to α, Equation (2) was rearranged as: Assuming a constant Manning coefficient (n), differentiating Equation (7) with respect to α leads to: By equating the left-hand side to zero, the expression becomes: According to Equations ( 8) and ( 9), the τ A 5/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 b 1 according to Equation ( 9).Accordingly, in Equation ( 8) the relationship between τ A 5/3 /dα and 0 is as follows: Therefore, τ A 5/3 has a maximum value when the slope angle is b 1 .Because τ A and τ A x have the same trend, τ A also has a maximum value when the slope angle equals b 1 .
As for the SIL, using the same method as above, by virtue of dτ A /dα = 0 (according to Equation ( 3)) it can be expressed as:

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 (t 1 and t 2 ) are different at different slope angles (α).Therefore, cumulative CSA should be analyzed according to the rainfall and soil infiltration processes.

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 (I 10 ) and soil stable infiltration rate [22].
Water 2023, 15, 3341 5 of 15 A previous equation for the CSA [15] was used to contrast with Equation (9): where α is the CSA, N 0 is the friction coefficient (nondimensional), d is representative grain size of soil (mm), r s 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.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: 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).

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: 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.

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 • (I 10 = 10.9 mm min −1 ), while the minimum discrepancy was 2 • (I 10 = 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.

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

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 Water 2023, 15, 3341 7 of 15 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.
direct comparison between Equation ( 9) and observed values.

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.

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).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).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.

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 t 1 , overland flow did not occur, and the cumulative CSA was 0 • .With increasing rainfall progression in the t 2 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 (t 2 ), before leveling off as the I/f ratio decreased the end of the phase.During the t 3 phase, there was no runoff and hence the cumulative CSA did not change.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].

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 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 t 4 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 t 5 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, I 5 (maximum 5 min rainfall intensity) is recommended as a more conservative value to estimate the maximum instantaneous CSA in practice.

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.

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.

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.

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.

Figure 1 .
Figure 1.Two derivations for unit width discharge change: (a) same horizontal projective length and (b) same incline length.

Figure 1 .
Figure 1.Two derivations for unit width discharge change: (a) same horizontal projective length and (b) same incline length.

Water 2023 ,
15, x FOR PEER REVIEW 7 of 16

Figure 2 .
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 .
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 3 .
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 4 .
Figure 4. Relationship between ratio of rainfall intensity and soil infiltration rate (I/f) and instantaneous CSA (critical slope angle).

Figure 4 .
Figure 4. Relationship between ratio of rainfall intensity and soil infiltration rate (I/f) and instantaneous CSA (critical slope angle).

Figure 4 .
Figure 4. Relationship between ratio of rainfall intensity and soil infiltration rate (I/f) and instantaneous CSA (critical slope angle).

Figure 5 .
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 6 .
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).
is cumulative shear stress at point A, T A 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, t 1 is time of overland flow initiation, and t 2 is time of overland flow cessation.

Table 1 .
Rainfall and soil parameters used for the CSA calculation.
2.3.2.Validation by Water Erosion Prediction Project (WEPP) Model SimulationsWEPP model results were used to validate the estimated cumulative CSA calculated by the mathematical equations (Equation (

Table 2 .
Parameters used in the WEPP (Water Erosion Prediction Project) model and Equation (4).