Strength Reduction Method in the Stability Assessment of Vegetated Slopes

Slope stability remains a significant and challenging problem in geotechnical engineering that applies to both natural and man-made slopes. Furthermore, the problem becomes more challenging when the slope is partially saturated and covered with vegetation. Estimating the stability of such a slope requires the application of advanced constitutive models and the incorporation of various factors that may influence the final result.

The presence of vegetation on the slope provides additional reinforcement by roots but also contributes to the development of negative pore pressures which are governed by the evapotranspiration process resulting from water uptake by plant roots. In these conditions, the hydro-mechanical soil behaviour should ideally be captured by a transient and coupled constitutive model, to more accurately reflect the response of the soil-root composite to changes in external conditions.

The most popular and widely used methods for slope stability assessment are those based on limit equilibrium (LE). The basic assumption in these methods is the pre-definition of the failure surface, which is established in the preliminary calculations phase. The sliding mass constrained by that failure surface is further divided into a set of parallel slices, and it is assumed that the limit equilibrium condition is obeyed at the base of each slice. Slope stability is determined in terms of the factor of safety (FOS), which is calculated as the ratio between the total resistance force and the total gravitational driving force. The classical and established LE methods, including those elaborated by Fellenius [1], Janbu [2], Morgenstern-Price [3], and Spencer [4], are widely used in engineering practice because they enable the relatively fast and efficient assessment of FOS, especially when slope geometry and soil conditions are straightforward.

However, when conditions and geometry are complex, and when it is difficult to pre-define the slip surface, the strength reduction method can be considered to determine slope stability [5,6,7]. This technique is often implemented in geotechnical engineering when the calculations are performed with finite element methods, especially together with the Mohr-Coulomb constitutive model. The soil strength parameters, namely effective cohesion and effective friction angle, are gradually decreased until failure occurs. The corresponding value of the strength reduction factor (SRF) is then used as an alternative to the factor of safety (FOS) calculated with LE methods.

The strength reduction method can be applied to slopes modelled with elasto-plastic soil constitutive models [8]. Practically, however, it is mostly used with the Mohr-Coulomb model [5]. Furthermore, Griffiths and Lu [7] also presented the possibility of applying the strength reduction method to slopes for which the stability assessment requires the consideration of partial saturation (infiltration, evaporation).

Schneider-Muntau et al. [9] proposed the application of the strength reduction method for the assessment of the stability of a slope built of the barodetic material, applying the constitutive model for clays [10, 11]. That concept can be utilized for slopes formed from overconsolidated clays and is based on the increase of the void ratio (e) and reduction of the friction angle, corresponding to the critical state (φ’_c).

This paper presents the numerical model of the vegetated slope simulated using the Modified Cam-Clay (MCC) constitutive model, which has been extended for unsaturated conditions by Tamagnini [12] and for root-reinforced soils by Świtała et al. [13], Świtała and Wu [14], Świtała et al. [15], and Świtała [16]. Therefore, this fully coupled hydro-mechanical model, implemented in finite element software Abaqus [15,17], enables the consideration of both mechanical root-reinforcement and evapotranspiration. To perform the stability analysis of the slope, the strength reduction technique has been implemented in the model. However, in the case of the Modified Cam-Clay-based model, the classical strength reduction method cannot be applied as it is done for the Mohr-Coulomb constitutive model, where high values of preconsolidation pressure are reflected in increased cohesion. In the model used in this paper, cohesion is not accounted for, whereas preconsolidation pressure is directly given and, together with the void ratio, plays a crucial role. Therefore, in the considered case, the strength reduction is based on the increase of the void ratio and reduction of the critical friction angle.

The coupled hydro-mechanical model for partially saturated soils reinforced by plant roots [13, 14] is based on the concept of Modified Cam-Clay [18], which has been extended for unsaturated conditions by Tamagnini [12]. The yield surface shape corresponds to that of a classical Modified Cam-Clay [18]. (1) f=q2−M2p′(pc′−p′=0 f = {q^2} - {M^2}\left[ {p^\prime(p_c^\prime - p^\prime} \right] = 0 in which p’ and q represent stress i.e., mean effective pressure and deviatoric stress, respectively. M denotes the slope of the critical state line and p’_c represents the effective preconsolidation pressure.

The partial saturation of the soil and mechanical root-reinforcement are accounted for in the model through the multi-hardening mechanism. This mechanism monitors the increment of effective preconsolidation pressure, which comprises three components corresponding to specific phenomena [13, 16]: (2) p˙c′=p˙cSAT′+p˙cUNSAT′+p˙cROOT′ \dot p_c^\prime = \dot p_{c\left( {SAT} \right)}^\prime + \dot p_{c\left( {UNSAT} \right)}^\prime + \dot p_{c\left( {ROOT} \right)}^\prime

The first component is responsible for the hardening resulting from the increment in volumetric plastic strains, denoted as ε˙vp \dot \varepsilon _v^p (following the original MCC model [18]): (3) p˙cSAT′=1+eλ−κpc′ε˙vp \dot p_{c\left( {SAT} \right)}^\prime = {{1 + e} \over {\lambda - \kappa }}p_c^\prime\dot \varepsilon _v^p where λ and κ are slopes of the normal compression and swelling line, respectively.

The second component reflects the hardening/softening induced by the changes in the soil degree of saturation S˙r \left( {{{{\rm{\dot S}}}_r}} \right) , and follows the solution proposed by Tamagnini [12]: (4) p˙cUNSAT′=−bpc′S˙r \dot p_{c\left( {UNSAT} \right)}^\prime = - bp_c^\prime{\dot S_r}

In the above equation, the parameter b, which determines the extent to which changes in the degree of saturation affect the preconsolidation pressure increment and the initial value of the preconsolidation pressure (as shown in Eq. 6), is a constitutive parameter that varies with soil type.

Finally, the last component describes the mechanical soil hardening induced by the root-reinforcement. In this case, the preconsolidation pressure increment is proportional to the increment of the, so-called, root activation strain ε˙r {\dot \varepsilon _r} , which is defined as a sum of volumetric and deviatoric strains, [13]. (5) p˙cROOT′=RPpc′eε˙rmrini \dot p_{c\left( {ROOT} \right)}^\prime = {R_P}p_c^\prime e\dot \varepsilon_r m_r^{ini}

The incorporation of the activation strain in the above equation allows for the determination of the level of root strength activation in response to the actual strain level in the soil. In the last component, two constitutive parameters are applied: R_P, which is dependent on the vegetation type, and mrini m_r^{ini} , defining the percentage of the mass of roots in the total mass of the considered root zone. Further details on the model's development, validation, calibration, implementation and application can be found in the papers [13,14,15,16].

In this paper, the strength reduction procedure is based on the one presented in the study by Schneider-Muntau et al. [9], which was implemented for the Barodesy model for overconsolidated clays. The procedure involves gradually increasing the void ratio e, which leads to a reduction in peak soil strength and corresponds to a reduction in preconsolidation pressure. This is justified by Fig. 1 (adapted from Atkinson [19]) which follows the principles of critical state soil mechanics. Simultaneously, the effective friction angle (φ’_c), corresponding to the critical state, is gradually decreased, leading to a reduction in strength at the critical state.

Figure 1.

In the presented model, only changes in preconsolidation pressure increment are accounted for in unsaturated soil conditions and root reinforcement. Therefore, for the partially saturated and vegetated slope being considered, the initial preconsolidation pressure will be enhanced. This evolution is given by the following formula [13, 15]: (6) pc,ini=pc0expb1−Sini+Rpmrini {p_{c,ini}} = {p_{c0}}exp\left[ {b\left( {1 - {S_{ini}}} \right) + {R_p}m_r^{ini}} \right] where p_c0 is the initial value of the preconsolidation pressure resulting from loading history, and S_ini is the initial degree of saturation.

The strength reduction technique implemented in this model considers an increase in the void ratio and a decrease in the friction angle corresponding to the critical state, which in the presented model is accounted for in terms of the slope of the critical state line M. The following relation is used [20]: (7) M=6sinϕc′3−sinϕc′ M = {{6sin\left( {\phi _c^\prime} \right)} \over {3 - sin\left( {\phi _c^\prime} \right)}}

The friction angle corresponding to the critical state is reduced as follows [9]: (8) tanϕc,r′=tanϕc′RF {\tan}\left( {\phi _{c,r}^\prime} \right) = {{{\rm{\;tan}}\left( {\phi _c^\prime} \right)} \over {RF}} where RF is the reduction factor.

In the considered procedure, the void ratio (e) is increased, leading to the decrease in the preconsolidation pressure. The maximum void ratio that can be obtained corresponds to the maximum void ratio at normal consolidation. The void ratio is increased as follows: (9) er=e−emaxRF+emax {e_r} = {{e - {e_{max}}} \over {RF}} + {e_{max}}

The strength reduction is implemented in Abaqus, [17] in the user material subroutine, along with the considered constitutive model. The reduction factor is treated as a field variable that increases with the increase of step time. The void ratio is a solution-dependent variable, updated according to the current stress state and the value of the reduction factor. The critical friction angle, reflected in terms of the slope of the critical state line, is a material constant.

The strength reduction implemented in the model described in the previous sections is applied in the numerical simulations, which aim to assess the stability of the vegetated slope. The model slope is covered with a uniform vegetation layer. The model assumes that the root zone is homogeneous, with enhanced properties, and thus, individual roots are not considered separately. Initially, the slope is unsaturated and the value of the negative pore pressure is set to −10 kPa. However, after the geostatic equilibrium step, the pore pressures are distributed in the slope and range from about −8 kPa at the base of the model to about −12 kPa near the slope surface. The unsaturated conditions are reflected in the initial value of the preconsolidation pressure, according to Eq. 6. The geometry of the slope, finite element discretization and boundary conditions are presented in Fig. 2. The parameters required to perform the calculations, i.e., constitutive model parameters and hydraulic characteristics of the slope soil, are listed in Table 1. The values of these parameters are assumed based on experimental investigations performed by Askarinejad et al. [21], Askarinejad [22] and Casini et al. [23], on silty sand.

Figure 2.

Furthermore, the value of the initial preconsolidation pressure, based on the past loading history, is set to 15 kPa, to assure the initial stability of the considered model slope. The parameter for partial saturation, which defines the level of the influence of the degree of saturation on the preconoslidation pressure increment, is assumed to be equal to 3.0.

The changes in the soil degree of saturation with increasing suction are considered in van Genuchten's equation, given by the following formula [24]: (10) Sr=Sres+1+uw−uaa1/1−nSsat−Sres {S_r} = {S_{res}} + \left[ {1 + {{\left( {{{{u_w} - {u_a}} \over a}} \right)}^{1/\left( {1 - n} \right)}}} \right]\left( {{S_{sat}} - {S_{res}}} \right) with the residual value of the degree of saturation S_res =0.412 and the degree of saturation for the fully saturated state S_sat=1.0 [22]. u_w-u_a is the difference between pore water and atmospheric pressure. Furthermore, the unsaturated hydraulic conductivity dependent on the current value of the degree of saturation can be calculated according to the formula proposed by Averjanov [25], with the fitting parameter β=1.7 for the considered soil [22].

The numerical simulations are performed in Abaqus [17] finite element software, with the model discretized using CPE8RP finite elements (eight-node plane strain quadrilateral with biquadratic displacement, bilinear pore pressure and reduced integration points). In the first step, the geostatic equilibrium is achieved, followed by the application of rainfall on the slope surface (as shown in Fig. 2) in the second step. The rainfall intensity is set to 11 mm/h, which is classified as heavy rainfall [26]. Immediately after the rainfall step, the strength reduction procedure is conducted, as described in the previous section.

In the case of slopes composed of material that has undergone high preconsolidation pressures, the slip surface is usually deep-seated, and thus, the enhanced strength of root-reinforced superficial zones does not significantly contribute to slope stability. On the other hand, for rainfall-induced landslides, which remain a serious threat in many regions worldwide, shallow, superficial soil layers are usually affected, resulting in translational landslides. These landslides are particularly dangerous in urban areas, and therefore, more efficient mitigation measures need to be considered to prevent them. Vegetation plays a significant role in reinforcing and stabilizing shallow slope soil layers [26,27,28,29,30] and has been used for decades as the simplest, most ecological, and most economically attractive measure for slope stabilization.

4.2.

Transient analyses

An additional set of calculations has been performed to assess the changes in slope stability affected by rainfall infiltration, considering various durations of the rainfall step. Time series calculations have also been conducted for the vegetated slope. The resulting values of the FOS for the reference slope and vegetated slope are presented in Fig. 3, for the considered durations of rainfall steps. The factor of safety decreases with increasing rainfall duration, which is observed for both the bare and vegetated slopes. In this analysis, the rainfall step duration varies between about 11 and 14 h. The intense rainfall affects the preconsolidation pressure significantly, contributing to the reduction of the strength in the superficial soil layer. However, the presence of vegetation may compensate to some extent for this reduction. The highest influence of mechanical root-reinforcement is observed when the highest strain increments are reported, which takes place in the first phase of the landslide initiation. The influence of vegetation reduces with longer rainfall periods.

Figure 3.

Fig. 4 shows the contour of the square root of the second invariant of the deviator strain at the end of calculations for the bare slope, presented on a deformed mesh. These results are similar for all the considered cases. On the other hand, Fig. 5 presents the corresponding results for the vegetated case (R_p = 4 and mrini=1.5% m_r^{ini} = 1.5\% ). In that case, the failure surface is not yet clearly detectable. However, the FE solution convergence is no longer achieved, which means that the force equilibrium is disrupted.

Figure 4.

Figure 5.

Until now, the strength reduction technique has mainly been applied in calculations using the Mohr-Coulomb constitutive model. The technique proposed for other models is described in the paper by Schneider-Mundau et al. [9], who implemented it for barodetic material. To reach the residual values of strength parameters, the void ratio is increased, and the critical friction angle is decreased until slope failure occurs (i.e., the analysis ceases to converge). The final value of the reduction factor is treated as the factor of safety. Increasing the void ratio results in a decrease in the preconsolidation pressure, while reducing the friction angle leads to a reduction in critical strength. This paper proposes the similar solution for the calculation of the factor of safety of the slope modelled with the Modified Cam-Clay model. Due to the fact that the extended version of the model enables modelling unsaturated soil conditions and taking into account root reinforcement and their effects on the preconsolidation pressure, the proposed technique can be used to assess the stability of vegetated slopes subjected to rainfall infiltration. The applied strength reduction technique is valid only for overconsolidated soils. Therefore, after the rainfall step and prior to the reduction application, the value of the preconsolidation pressure should be checked.

The results of numerical simulations demonstrate that the model can estimate the factor of safety of slopes where stress-strain relations are modelled using the Modified Cam-Clay model. Additionally, the stability of vegetated slopes has been assessed following a rainfall episode of a certain duration. Although the model is still in testing, initial reasonable results can be reported.

The model is based on several assumptions, including a homogeneous vegetation layer with enhanced properties as described in the ‘Model description’ section. The void ratio and saturated hydraulic conductivity of that layer are assumed to be equal to the corresponding parameters of the bare soil. However, to accurately reflect the nature of this layer, the values and variability of these parameters with respect to different conditions, should be investigated in more detail. Furthermore, to obtain root constitutive parameters for a certain vegetation species, the model calibration based on results from triaxial or direct shear tests performed on root-permeated soil samples is necessary. [13].

The implemented strength reduction technique allows assessment of the stability of vegetated slopes. However, to observe the enhancement caused by the root-reinforcement, the root strength must be mobilized, i.e. strain develops in the root zone. If the slope is built from cohesive materials and is not affected by any disruptive factor, the strength reduction results in the development of deep-seated failure surfaces that do not intersect the root zone. Consequently, activation strains do not develop in the superficial soil layer containing the roots, and the factor of safety remains unaffected.

The proposed idea for strength reduction should be extensively tested in future research. For instance, the results from the finite element method could be compared with those from limit equilibrium methods, or with strength reduction using appropriate parameters for Mohr-Coulomb. Additionally, different rainfall intensities and durations can be investigated, and the results compared. Furthermore, the effects of different shapes and depths of root zones, as well as varying root parameters, could also be explored.

Prevention against rainfall-induced landslides and early warning systems are crucial to avoid their dangerous consequences. The performed numerical simulations have confirmed the well-known fact that the factor of safety of a slope subjected to rainfall decreases with increasing rainfall duration. The presence of a homogeneous vegetation layer may contribute to increasing the stability of the superficial slope zone and, thus, a higher factor of safety of the slope subjected to rainfall. This translates into additional time for warning and evacuation, making responsible and wise slope stabilisation with vegetation a measure with positive outcomes.