Estimating the critical shear stress for incipient particle motion of a cohesive soil slope

The critical shear stress is a vital reference indicator for soil erosion. Soil erosion will occur when soil slope suffers from a exceed shear stress, and then causing soil loss and destruction of soil structure. In this work, an equation was proposed based on the force equilibrium of a single particle to estimate the critical shear stress for incipient particle motion of a cohesive soil slope. This formula is characterized by its physical significance, and the critical shear stress for incipient slope soil motion can be easily calculated when the soil properties and the slope angle are known. Moreover, the seepage-runoff coupled model and the excess shear stress equation are introduced in this paper. Two parameters, namely the weight of rushed soil particles and the discharge of water, must be measured in the scouring tests. Through linear regression, the tested τc-values were obtained to validate the τc-values calculated by the formula derived from the critical shear stress. In addition, two other formulas were compared with the derived formulas, which considered more parameters with physical significance. Finally, the influence of all parameters on the critical shear stress was analyzed: the porosity of the soil, the specific gravity of the soil and the slope gradient had less influence on the critical shear stress; the critical shear stress was negatively influenced by the particle diameter and positively influenced by the internal friction angle of the soil.

where the shear stress with general transport was defined as the critical shear stress 19 . Wilcock and Southard considered the critical shear stress to be approximately equal to the shear stress that produces a small transportation rate 20 . Choi and Kwak analyzed three different incipient modes of particles: rolling, sliding, and lifting 15 . Shvidchenko et al. took the fractional transport rate as an indicator for judging the incipient conditions 21 . Singh et al. used visual observations along with the quantitative measurement of transport rates as an incipient criterion 16 . The incipient standards described above are unquantifiable and diverse, which results in measured critical shear stresses that may not turn out to reflect the actual conditions of the critical shear stress.
Direct measurement of the critical shear stress requires complex experimental equipment. Normally, the critical shear stress can be indirectly measured through scouring tests by using the excess shear stress equation. The results are empirical and lack of physical significance. The main contents in this paper include: (1) proposing a theoretical formula to estimate the critical shear stress of soil slope, which is based on the soil properties (i.e., soil unit weight, porosity, particle size, specific gravity, and internal friction angle) and the slope angle; (2) introducing the seepage-runoff coupling model and the excess shear stress equation to validate the derived formula of the critical shear stress; (3) conducting the deep analysis to explain how are the soil properties and the slope angle influence the critical shear stress.

Methods
Theory. Expression derivation of the critical shear stress. In this section, we derive a new formula for the critical shear stress based on the force equilibrium of a single particle (Fig. 1).
In addition to gravity, the normal force between particles and cohesion, the particle is also subjected to forces induced by water flow, including buoyant force, seepage force, drag force and uplift force (Fig. 1).
The expressions for effective gravity F G , uplift force F L and drag force F D are given as follows 12,22 : where d is the particle diameter (m), γ' s is the unit weight of the soil particle (kN/m 3 ), C L and C D are the coefficients (dimensionless) of uplift force and drag force, respectively, ρ w is the density of water (kg/m 3 ), and u b is the velocity (m/s) of the water flow near the slope surface. The seepage force can be expressed as: where n is the porosity of soil (dimensionless). For cohesive particles, it can be considered that cohesion is generated by the contact of film water around the particles, and the expression is written as 23 : where q 0 is the bonding force per unit area, with q 0 = 1.3 × 10 10 N/m 2 ; δ 0 is the diameter of a water molecule, with δ 0 = 3 × 10 −10 m; δ 1 is the thickness of the film water, with δ 1 = 4 × 10 −7 m; and m is half of the closest gap distance between the particles, which can be calculated by Eq. (7) 23 .
(2) Figure 1. Force analysis for a single particle on the surface of a soil slope. where φ is the friction angle of the soil (°). In this context, sliding failure is adopted to describe the incipient slope soil motion. When the sliding force is equal to the resistance, the soil particles are in a critical state (see Eq. 10).
We substitute Eqs. (2), (4), (5), and (9) into Eq. (10), and then the expression for the critical surface shear velocity u bc can be written as: where G s is the specific gravity of the soil particle, with G s = γ' s /γ w ; γ w is the unit weight of water (kN/m 3 ); α is an integrated parameter (m 3 /s 2 ), with α = q 0 δ 3 0 (1/m 2 − 1/δ 2 1 )/ρ w . The friction velocity u * is empirically related to the surface shear velocity u b and their relation is written as 24 : where β is an empirical parameter whose value varies between 5.6 and 8.51. Herein, β = 5.6.
Furthermore, the surface shear stress τ a is related to the friction velocity u * , which is expressed as 25 : Thus, the parameters in the critical state can be expressed as follows: where u *c is the critical friction velocity (m/s) and τ c is the critical surface shear stress (Pa). Combining Eq. (11) with Eq. (14), the critical surface shear stress τ c can be given by the following equation: This formula can be employed to estimate the critical shear stress for the incipient motion of the soil.
Excess shear stress equation. According to Eq. (1), the excess shear stress equation can be rewritten as follows: In this formula, the erosion rate ε and the surface shear stress τ a are the parameters that need to be measured. In this way, the erodibility coefficient k d and the critical shear stress τ c of the soil can be obtained based on linear regression 26 .
In addition, the erosion rate ε is easy to measure, as it is defined by: where M is the weight (kg) of sediment washed away, t is the scouring time (s), and S is the area of the eroded zone (m 2 ). Moreover, it is vital to determine the surface shear stress τ a and have it measured. Herein, we modeled a simplified soil slope (Fig. 2a) and considered runoff-seepage coupling 18 . The soil slope has a slope gradient of θ, and the soil properties include porosity n and permeability K. In this model, the motions of runoff and seepage are described by the Navier-Stokes equation and the Brinkman-extended Darcy equation, respectively. The runoff and seepage velocities are denoted as u and v, respectively, and b denotes the thickness of the soil layer. In addition, there is impermeable bedrock under the soil layer.  18 derived the velocity distribution in the y direction, which is expressed as Eq. (18), and the velocity distribution is shown in Fig. 2b.
where u and v are the velocities of runoff and seepage (m/s), respectively; θ is the slope gradient (°); h is the depth of runoff (m); b is the thickness of the soil layer (m); n and K are the porosity and permeability of the soil (m 2 ), respectively; and J is the hydraulic gradient, with J = sin θ.
Then, according to Newton's law of friction, the expression of surface shear stress can be derived as: In general, the average velocity V of the runoff cross-section is used to describe the incipient sediment motion, which is given by Manning's equation as follows: where V is the average velocity (m/s) of the runoff cross-section; λ is Manning's coefficient, which is used to indicate the slope roughness; and R is the hydraulic radius (m) of the runoff in the testing flume, which is approximately equal to the runoff depth h.
According to Eq. (20), the discharge per unit width of runoff can be written as: where q is the discharge per unit width (m 2 /s). Combining Eq. (19) with Eq. (21), we can obtain an alternative expression for the surface shear stress: Therefore, in laboratory tests, the surface shear stress τ a can be measured indirectly when the discharge per unit width is known. According to a series of scouring tests, a fitted line can be obtained, resulting in the k d and τ c -values. Furthermore, the obtained τ c -values lack physical significance.
Experiment. Experimental device. The scouring device (Fig. 3) mainly consists of a water tank, a guiding flume, a testing flume, and a collection barrel.
τ a = 2γ w  flume is used to place the testing material and observe the scouring process. A bolted connection exists between the guiding flume and the testing flume, which allows the slope gradient to be adjusted by changing the height of the jack to rotate the testing flume. (d) Collection barrel: This barrel is used to collect the discharged water and the sediment washed away. The collected water and sediment are weighed in a certain period to obtain the discharge per unit width and the scouring rate.
Materials and parameters. The testing material was taken from a highway slope site on Chengdu Third Ring Road, China. The soil material object after drying and smashing is shown in Fig. 4a, and the grading curve is shown in Fig. 4b. In addition, the physical parameters of the material and other related parameters can be found in Table 1.
Testing process. First, the flow rate of the valve should be calibrated by marking the valve position for a certain flow rate (0.4 L/s, 0.5 L/s and 0.6 L/s). Second, the soil material should be filled in the testing flume with the flat slope surface parallel to the flume surface. Then, the desired slope gradient can be adjusted through the length of the jack. Before the scouring experiment begins, a low-discharge flow should be ensured and infiltrates into the soil layer, which will stop when runoff from the slope surface merely begins to form. The scouring process will last for 10 min, but the experiment will end if the slope fails. During the scouring process, the collection barrel is replaced every 2 min. Finally, when the scouring process is finished, the scoured sediment collected in the barrel is dried and then weighed. Using the weight of the sediment M, the scouring rate can be calculated by Eq.   Table 2.
According to Eq. (17), the erosion rate ε can be calculated for different times t (Fig. 5). For each curve, the ε-values varied slightly with time, indicating that the slope was steady. Thus, we adopted the average ε-value of the curves to represent the erosion rate of the scouring process.
The four fitted lines (Fig. 6) were obtained through a series of ε − τ a points categorized by the slope gradient. These lines contain the coefficient values of the excess shear stress equation, as derived in Table 3.
By comparison, theoretical values of the critical shear stress are calculated from Eq. (15), which are also given in Table 3. The fitted results for the critical shear stress have small relative errors (< 15%) with the theoretical results for θ = 26.6°, 29.7° and 33.7°, while the relative error when θ = 38.7° is much larger than 20%.     critical shear stress. We avoided complex equipment and unquantifiable criteria for the critical shear stress and found this simple method to determine the critical shear stress by simply measuring the weight of soil particles being washed away and the water discharge. This method is discussed as follows.    (16), was applied to certify the accuracy of the simple measurement. Here, the erodibility coefficient k d adopts the fitted values for the corresponding slope gradient, while the critical shear stress τ c adopts the theoretical results of Eq. (15). The calculated results were compared with the testing results in Fig. 7.
The results indicate that the slope gradient and the flux have a positive effect on the erosion rate ε. Taking the calculated results as an example, when the slope gradient is 26.6°, the result with a flux of 0.5 L/s is 1.36 times greater than that with a flux of 0.4 L/s; when the flux is 0.5 L/s, the result with a slope gradient of 29.7° is 1.08 times greater than that with a slope gradient of 26.6°. In addition, it was found that the theoretical results obtained from Eq. (16) are larger than the testing results.
In the analysis of the results, most of the error results were found to be less than 15%, indicating a reliable estimation through this measurement. However, when θ = 38.7°, there is a relatively large difference between the fitted and theoretical τ c -values, with a relative error of more than 20%, and the ε-values have a similar error result. This may be due to different flow patterns, as laminar flow is assumed in the runoff-seepage coupled model, but the flow motion is actually turbulent. However, the erodibility coefficient k d in erosion rate estimation needs to be measured by scouring tests.
Regarding the comparison with published formulas. To validate whether Eq. (15) can be widely used in other references, test data from the reported articles were collected (see Table 4). The experimental data from these references are shown in Fig. 8.
Through linear regression, we obtained the critical shear stress and the erodibility coefficient of the five samples in Table 5. In addition, the τ c -values calculated by Eq. (15) are also given in Table 5.
By comparison, Eq. (15) also shows advantages in accurately estimating the critical shear stress of soils with a maximum relative error value of 31.89% and a maximum absolute error value of 0.1326 Pa.
Moreover, we adopted the following two empirical formulas for comparison purposes 31 : where d 50 is the mean diameter of the soil (m), and P c is the percentage of clay by weight (%). Furthermore, we employed the test data obtained by Lei et al. 32 . The soil sample was a typical silt-clay soil from the Loess Plateau of China. The clay particles, whose size is less than 0.01 mm, account for approximately 56% of the weight. Table 6 presents the τ c -values for different slope gradients from Lei et al. 32 and the calculated results of Eqs. (15), (23) and (24).
The τ c -values calculated from Eq. (15) were found to be the most consistent with the experimental τ c -values among the three equations. In addition, the calculated τ c -values from Eq. (15) vary with the slope gradient, while   www.nature.com/scientificreports/ the results calculated by the other two equations do not. On the one hand, the last two formulas are incomplete, as they focus on only one parameter. On the other hand, both formulas are empirical and lack physical significance. Moreover, Eq. (15) is derived based on the force equilibrium of a single particle. In the derivation process, multiple parameters are considered, namely particle diameter d, porosity n, specific gravity G s , slope gradient θ, and internal friction angle φ. Therefore, it can be concluded that the critical formula (Eq. 15) derived in this paper is applicable to estimate the critical shear stress for the incipient particle motion of cohesive soil slopes.
Regarding sensitivity. According to Eq. (15), the critical shear stress τ c is related to the particle diameter d, porosity n, soil unit weight γ s , slope gradient θ, internal friction angle φ, and integrated parameter of bonding force α. In addition, α is influenced by d and γ s from Eqs. (6) and (7). Here, these parameters are discussed to determine their influences on the estimation of the critical shear stress τ c , and the parameter values are given in Table 7.
Due to the influence of d and γ s on the bonding force, the parameter α was first analyzed with the increase of d. The slope gradient θ = 10° and the values of other parameters can be found in Table 1. The results of α varying with d and γ s can be seen in Fig. 9. Generally, α is negatively correlated with d and positively correlated with γ s . After the above discussion, the variation in the critical shear stress τ c with d and γ s can be further studied. Correspondingly, the α-curves should be used properly in the calculation process of τ c with similar soil unit weights.
The analyzed results are plotted in Fig. 10a. τ c is positively correlated with γ s , and it first decreases and then increases with increasing d. When d is less than 1 mm, the difference between the four curves is obvious. The critical shear stress τ c decreases from 11.54 to 1.20 Pa with d = 0.01 mm when γ s decreases from 18.0 kN/m 3 to 17.0 kN/m 3 , and it decreases from 0.40 to 0.07 Pa with d = 0.1 mm. The decreasing gap of τ c is enlarged with decreasing of soil particle diameter. Nevertheless, there is almost no difference between the four curves when d Table 5. Erodibility coefficient and critical shear stress obtained by linear regression and comparison between the fitting critical shear stress values and the calculation results from Eq. (15).  The sensitivity analysis of τ c for the other three parameters (porosity n, slope gradient θ, and internal friction angle φ) is also plotted in Fig. 10b-d. The α-values refer to the α-curve in Fig. 9 with γ s = 18.0 kN/m 3 .
In Fig. 10b, τ c shows a slight difference with diverse porosity values, and the difference can be observed when d ≥ 1 mm. In addition, τ c is negatively correlated with n. When n increases from 0.30 to 0.35, τ c decreases from 3.88 to 3.70 Pa with d = 10 mm and decreases from 0.42 to 0.40 Pa with d = 1 mm. The variation difference is enlarged as d increases.
In Fig. 10c, τ c can be obviously influenced by θ when d ≥ 0.1 mm, and τ c is negatively correlated with θ. When θ increases from 8° to 12°, τ c decreases from 0.41 to 0.38 Pa with d = 0.1 mm and from 0.48 to 0.26 Pa with d = 1 mm. Similarly, the variation difference is enlarged as d increases.
In Fig. 10d, τ c can be obviously influenced by φ when d ranges from 0.01 to 10 mm, and τ c is positively correlated with φ. When φ increases from 30° to 35°, τ c increases from 11.30 to 13.35 Pa with d = 0.01 mm, from 0.23 to 0.29 Pa with d = 0.2 mm, and from 3.16 to 4.85 Pa with d = 10 mm. The variation difference first decreases and then increases with d.
According to the sensitivity analysis, the critical shear stress τ c of slope cohesive soil is mainly influenced by the soil particle diameter d, internal friction angle φ, and soil unit weight γ s while that of slope cohesionless soil is mainly influenced by the soil particle diameter d, internal friction angle φ, porosity n and slope gradient θ.

Conclusions
In this paper, a formula has been derived for estimating the critical shear stress for the incipient particle motion of a cohesive soil slope. This formula is based on the force equilibrium of a single particle, which makes the derived formula physically significant.
To validate the derived formula, a series of scouring tests were conducted, and the excess shear stress equation was introduced to estimate the tested values of the critical shear stress. The results show that the calculated results from the derived formula are in good agreement with the tested results. In addition, two other formulas were employed for comparison with the derived formula. The derived formula shows great advantages in considering multiple parameters and has obvious physical significance.
The parameters influencing the critical shear stress were analyzed. The soil porosity, soil specific gravity and slope gradient have less influence on the critical shear stress. In addition, the critical shear stress is negatively influenced by the particle diameter and positively influenced by the internal friction angle.