Numerical Analysis of an Explicit Smoothed Particle Finite Element Method on Shallow Vegetated Slope Stability with Different Root Architectures

Planting vegetation is an environmentally friendly method for reducing landslides. Current vegetated slope analysis fails to consider the influence of different root architectures, and the accuracy and effectiveness of the numerical simulations need to be improved. In this study, an explicit smoothed particle finite element method (eSPFEM) was used to evaluate slope stability under the influence of vegetation roots. The Mohr–Coulomb constitutive model was extended by incorporating apparent root cohesion into the shear strength of the soil. The slope factors of safety (FOS) of four root architectures (uniform, triangular, parabolic, and exponential) for various planting distances, root depths, slope angles, and planting locations were calculated using the shear strength reduction technique with a kinetic energy-based criterion. The results indicated that the higher the planting density, the stronger the reinforcement effect of the roots on the slope. With increasing root depth, the FOS value first decreased and then increased. The FOS value decreased with an increase in slope angle. Planting on the entire ground surface had the best improvement effect on the slope stability, followed by planting vegetation with a uniform root architecture in the upper slope region or planting vegetation with triangular or exponential root architecture on the slope’s toe. Our findings are expected to deepen our understanding of the contributions of different root architectures to vegetated slope protection and guide the selection of vegetation species and planting locations.


Introduction
Landslides and debris flows are common natural disasters that cause environmental damage, human casualties, and economic losses [1,2]. To tackle this problem, measures, such as nailing, vegetation, ground improvement, geosynthetic reinforcement, and improved drainage have been taken. Among these, vegetation is an economically, sustainable, and environmentally friendly bio-remediation technique [3,4]. Further studies on the effects of vegetation on slope stability are essential.
It is now understood that vegetation contributes to the stability of civil infrastructure, including shallow slopes (i.e., slopes with vertical depths less than 2 m), road and railway slopes, dams, embankments, and dykes [5][6][7][8][9]. Vegetation roots (e.g., from grasses, shrubs, and trees) are believed to stabilise slopes and slow the large-scale movement of landslides by strengthening soils with increased cohesion [10][11][12][13]. Planting trees on slopes can reduce the occurrence of shallow landslides by up to 95% compared with similar areas without trees [14,15].
To date, a large body of literature has documented studies that have focused on quantifying the contribution of roots to soil shear strength. These studies included in situ The PFEM is based on an updated Lagrangian (UL) fashion for modelling the motion of a continuum medium. The continuum medium is discretised into a set of Lagrangian nodes (particles) that contain and transmit all of the information present. The computational mesh was built using the Delaunay triangulation technique, and the boundary of the computational mesh was identified using the alpha-shape approach. A mesh was then used to solve the governing equations. However, in the PFEM, excessive mesh distortion is avoided by frequent remeshing, and the state variables (i.e., stress and strain) are mapped from old Gauss points to new ones. This mapping procedure inevitably introduces errors, which increase the complexity of the calculation process [30,32,42].
The SPFEM [50,51] uses a strain smoothing technique for nodal integration based on the PFEM to achieve the balance of governing equations at nodes or particles. In the SPFEM, all of the field variables are calculated at particles instead of Gauss points in the PFEM, which can avoid information transfer between Gauss points and particles, thus reducing the calculation error. In addition, the SPFEM possesses the upper-bound property and provides a conservative estimate for problems in geomechanics. Finally, low-order triangle elements can be used directly without volumetric locking [30,32,42]. The SPFEM can consider the entire dynamic failure process of a slope in slope-stability analysis and therefore simulate the large deformation and post-failure of soil to obtain a more reliable FOS value [32].
In SPFEM, the computation domain is discretized into strain smoothing cells associated with nodes. The physical volume, Ω, is correspondingly discretized into particles. As shown in Figure 1, The smoothing cell associated with the particle k is created by connecting sequentially the mid edge point to the central points of the surrounding triangular elements of particle k. A strain smoothing operation is then performed based on the set of smoothing cells that are created based on the triangulation mesh [32]. In the smoothed strain technique, the area, A k , of the smoothing cell, Ω k , associated with node k is calculated by: where N k is the number of elements related to particle k and A j is the area of the jth element. The smoothed strain matrix is calculated by: where B j represents the strain gradient matrix used in the standard FEM for the jth element. More details regarding the strain smoothing technique are available in Liu et al. [51] and Chen et al. [52].
The discretization of the computational formulations of eSPFEM is now briefly presented. The motion of a continuum can be described as: where, ρ is the material density, a is the acceleration, σ is the Cauchy stress tensor, and b is the specific body force density. Considering the principle of virtual displacement and the divergence theorem, the weak form is expressed as: where u is the displacement vector, Ω represents the configuration domain, S represents the boundary, and τ S is the prescribed traction. After the node-based discretization, the above equation reads: where T is the total number of nodes in the computation domain. The above discrete form can then be written in the vector or matrix form as: in which: where F ext and F int are termed as the external and internal forces, respectively, and M is the diagonal mass matrix. A typical computational cycle of eSPFEM is shown as follows [43]: (1) Generate mesh using Delaunay triangulation and alpha shape technique; (2) Acquire basic data of elements and nodes; (3) Calculate smoothed strains of nodes; (4) Update stresses of nodes through constitutive integration; (5) Calculate internal forces of nodes; (6) Update velocities and position.
Compared with the implicit SPFEM [30], the eSPFEM adopts an explicit time-integration scheme, which has a more concise formulation, lower computational cost, and wider applications. We used a self-developed code that aims to develop a GPU-accelerated SPFEM for large deformation analysis in geomechanics based on CUDA [43], which is released by NVIDIA to perform high-performance computing and has gained popularity rapidly in geomechanics recently [53,54]. Details of the eSPFEM theory can be found in Yuan et al. [32,42] and Zhang et al. [43,55,56].

Modelling the Mechanical Effect of Roots
Plant roots extend into the soil matrix and form a soil-root composite material with a mechanical effect that enhances the shear strength of the soil [57]. This effect is often considered as additional soil cohesion (known as apparent root cohesion) [23,58,59]. Different theoretical models have been created to estimate the mechanical effect of roots on slope stability, such as the Wu model [17,60] and the fibre bundle model (FBM) [61]. The Wu model assumes that all roots are moved and broken simultaneously, which is not true in reality. However, it is the most commonly used model owing to its simplicity. The Wu model is a perpendicular root-strengthening model established using two variables: root area ratio (RAR) and root tensile strength [48].
The Wu model is used to describe the increase in soil shear strength caused by the mechanical effect of the roots, as follows: where c r is the additional soil cohesion and ζ represents the correction factor, which takes into account the influence of roots breaking progressively in reality on soil shear strength; in the present study, ζ equals 0.4 [62]. T r is the root tensile strength. For simplicity, the effects of the root diameter distributions were ignored, and a constant root tensile strength was considered [63]. Available current data on tree roots indicate that T r ranges from 5 to 60 MPa [11,64]; in the present study, T r equals 20 MPa.
R f stands for the root orientation factor, which is defined as follows: where θ represents the angle between the root and failure surface when the root breaks, and ϕ is the effective angle of internal friction. It should be noted that roots do not always grow perpendicular to the slope's surface, as their growth is influenced by ambient conditions (e.g., gravity and nutrition) [57]. In many cases, θ is between 48 • and 72 • at failure, so the range of sin θ + cos θ tan ϕ is narrow, i.e., approaching 1.2 [11]. In the present study, R f equals 1.2.
The RAR (in Equation (1)), which is defined as the proportion of the cross-sectional area of the soil occupied by roots, is determined as follows: where A r and A represent the cross-sectional area of root and soil, respectively. In the present study, RAR equals 0.45% [64]. Many researchers have estimated the value of c r for various vegetation species growing in different environments, and its typical values vary from 1.0 to 94.3 kPa [29]. According to Equation (1), the original c r score in this study is 43.2 kPa.

Mechanical Effect of Root-Soil on Shallow Slope Stability
Most of the results from the direct-shear tests showed that roots increased the cohesion [60], whereas the friction angle remained mostly unchanged [64]. For undrained loading, the shear strength of saturated vegetated soils, s, can commonly be determined by incorporating the additional soil cohesion, c r , into the Mohr-Coulomb failure criterion [64,65], which can be modified as follows: where, c is effective cohesion and σ n is total normal stress.
The shear strength reduction technique was adopted to analyse the slope stability, which is defined as follows: where s f and ϕ f are the reduced shear strength and the reduced effective angle of internal friction, respectively. SRF is the shear strength reduction factor. First, the initial value of the SRF was set to an adequately low value to keep the slope stable under gravitational loading. Then, the value of the SRF increases gradually until the slope becomes unstable and failure occurs. The critical value of the SRF leading to slope instability is considered to be the FOS of the slope [30]. The spatial location of thick roots determines the arrangement of related thin roots; thus, root distribution presents a high degree of complexity [8]. Existing methods for quantifying root architecture include the extraction of roots, the complete washing of soil, and image analysis of roots [66]. The spatial distribution of roots is an important factor in determining the reinforcement behaviour and mechanical properties of roots, and the generalisation of root morphology is essential for evaluating the influence of vegetation on slope stability [9,67].
Based on experimental observations, researchers have summarised four typical patterns of root geometry: uniform distribution [10,68], triangular distribution [68], parabolic distribution [69], and exponential distribution [10]. Table 1 lists the characteristics of different root architectures, typical species, and their growing regions. Figure 2 shows the root architecture of four different real roots.

Root Architectures Characteristics Typical Species Growing Regions
Uniform distribution A root system with a large taproot and large horizontal lateral roots Aleppo pine [10] Mediterranean, Southern Europe, Asia, and North Africa [70] Pulsatilla pratensis [68] Sub-polar areas of Europe, Asia, North America Central, and Eastern Europe [71] Triangular distribution A root taproot system with small lateral roots Trigonella balansae [68] Europe, and Asia [72] Parabolic distribution A concentrated root system Cynodon dactylon [69] North Africa, Asia, Australia, and Southern Europe [73] Exponential distribution A plate-shaped root system Beech and Mature oaks [10] Europe, and North America [74,75]  Root biomass can be expressed by root volume, mass, area, or length, and the most commonly adopted measure is the soil area occupied by the roots [76]. In order to achieve a fair comparison, the symmetric parts of the four root architecture profiles are normalised to the same unit area. Figure 3 describes their normalised function curves. It is assumed that these root architectures have the same root depth and are homogeneously distribute in the root zone.

Root Architecture Functions on Slope
The coordinate system on the slope can be defined as x o y , and the extent between the bottom of the root zone and the slope's surface is defined by the root depth, z r . Figure 4 describes the boundaries of the four different root zones, and the green areas indicate the root zones. The mathematical functions derived from the root zone boundaries are as follows: 1.
Uniform distribution

Root Architecture Functions after Coordinate Transformation
The coordinate transformation formula from the coordinate system x o y to the rectangular coordinate system xoy is assembled as: where α is the complementary angle of the slope angle and β, C 1 , and C 2 are parameters related to the planting position. Equation (10) can be expressed as: According to Equations (15) and (20), the scope boundaries of the uniform distribution were obtained.
According to Equations (16) and (20), the scope boundaries of the triangular distribution were obtained. when z r 2 − 2 tan α > 0, According to Equations (17) and (20), the scope boundaries of the parabolic distribution were obtained.
According to Equations (18) and (20), the scope boundaries of the exponential distribution were obtained. (21)-(25) are discussed under three different conditions: vegetation growth on the slope's surface, upper slope region, and lower slope region. Figure 5 shows the geometric parameters of the vegetated slope (taking a uniform distribution as an example), green areas indicate the root zones. When vegetation is planted on the slope's surface: where a is the horizontal interval between root centres on the slope's surface, d a is the horizontal distance between the root centre of the highest vegetation and the near edge of the upper slope region, d a = a/2, i tree is the amount of vegetation on the slope's surface, 1 ≤ i tree ≤ l 2 /a, and l 2 /a is an integer.
When vegetation grows on the upper slope region: where b is the horizontal interval between root centres on the upper slope region, d b is the horizontal distance between the root centre of the vegetation closest to the slope's surface and the top edge of the slope's surface, d b = b/2, j tree is the number of vegetation on the upper slope region, 1 ≤ j tree ≤ l 1 /b, and l 1 /b is an integer.
When vegetation is planted on the lower slope region: where g is the horizontal interval between the root centres in the lower slope region, d g is the horizontal distance between the root centre of the vegetation closest to the slope's surface and the bottom edge of the slope's surface, d g = g/2, k tree is the amount of vegetation in the lower slope region, 1 ≤ k tree ≤ l 3 /g, and l 3 /g is an integer.

Numerical Implementation
The geometry and boundary conditions of the slope-stability problem are shown in Figure 6. The boundary conditions were set as rollers along the left and right vertical boundaries and were fully fixed at the base. The soil behaviour was modelled by an elastic-perfectly plastic Mohr-Coulomb material, and the material properties are listed in Table 2. The material parameters were obtained by referring to [30,77]. The slope was divided into lower bedrock and upper soil. The root zone is the soil elements influenced by vegetation [29]. It is assumed that two adjacent root zones may overlap without affecting root distribution. Figure 7 shows root zones and overlapping root zones of different root architectures in vegetated slopes (slope angle of 45 • , root depth of 1 m, root zone area of 2 m 2 , planting distances of 2.5 m, and planting location on the slope's surface). To guarantee both the numerical accuracy and the computational efficiency, a nonuniform initial particle distribution is assumed. A very fine mesh discretisation is used for the soil and root zones near the slope's surface, and a coarser mesh discretisation is used for the bedrock, which is situated away from the slope's surface. The meshes adopt six-node triangular elements. Table 3    In this study, the eSPFEM is utilised to simulate the failure process of the entire vegetated slope and predict the large deformation behaviour and final deposit of the slope failure. A kinetic energy-based criterion [32] is adopted to analyse the stability of vegetated slopes. This approach is based on the relation of the kinetic energy of the slope with the simulation time [32]. The peak value of kinetic energy curve can be considered as an indicator of the critical state, because the failure of the slope is related to large deformation, and the kinetic energy increases significantly after failure occurs. The value of SRF at the peak of the kinetic energy curve is interpreted as the FOS of the slope, and kinetic energy reaches a steady state in a short period after the peak.
The simulation is divided into two stages: first, the displacements of all of the particles are fixed, and gravity is applied to all of the particles to achieve the initial stress field; then, the particles are allowed to move, and after the initial shear failure, the unstable soil moves and reaches an equilibrium state at a new slope configuration. The horizontal displacement of point P, which is located at the slope's toe, is monitored (see Figure 5).
Slope failure was calculated using the shear strength reduction technique by gradually increasing the SRF. For different SRFs, the simulation was implemented with eSPFEM, and a physical time of 6 s [32] was considered for each SRF to obtain a new steady state for the slope after failure. To evaluate the effect of roots on slope stability, five variables (root architecture, planting distance, root depth, slope angle, and planting location) were systematically varied to calculate FOS values under different conditions.

Comparison of the Instability of the Vegetated Slope and Bare Slope
Case 1 compares the failure mechanisms of the bare and vegetated slopes. The eSPFEM was used to analyse the stability of bare and vegetated slopes with a slope angle of 45 • (planting location is the slope's surface, root depth of 1 m, root zone area of 2 m 2 , planting distance of 2.5 m, and a uniform root distribution) for various SRFs. Figure 8 shows the evolution of the kinetic energy with time. Figure 9 shows the variation in the maximum horizontal displacement at the slope's toe with time. Figure 10 shows the equivalent plastic strain and final configurations for various SRFs, and the grey zones indicate the initial root areas.  For the bare slope (Figure 8a), when SRF ≤ 1.3, the slope is stable, and no obvious kinetic energy occurs. When SRF ≥ 1.4, the kinetic energy-time curve exhibits a peak value, and then the kinetic energy reaches a steady state within a short time. Therefore, the critical SRF value was 1.3, which was regarded as the FOS of the slope. Displacement at the slope's toe (Figure 9a) also occurred when SRF ≥ 1.4. In addition, from the results of the equivalent plastic strain (Figure 10a), it can be seen that when SRF = 1.2, a plastic strain develops at the weak band between bedrock and soil as well as at the slope's toe. When SRF = 1.3, a narrow local band of plastic strain was detected, and the plastic zone subsequently extended from the slope's toe to the slope top. The plastic strain band was dark near the top of the slope, indicating that the slope was not completely damaged. When SRF = 1.4, a continuous band of plastic strain localisation was observed, and slope failure occurred, accompanied by a large increment in displacement. The deformation of the slope was mainly concentrated on the slope's surface, and the maximum deformation occurred at the centre of the slope's surface. For the vegetated slope (Figure 8b), no obvious kinetic energy appears when SRF ≤ 1.5, whereas the kinetic energy-time curve shows a peak value when SRF ≥ 1.6. The FOS of the slope was 1.5. A large deformation displacement at the slope's toe (Figure 9b) also occurs during FOS ≥ 1.6. Moreover, according to the equivalent plastic strain diagram (Figure 10b), when SRF = 1.6, the slope produces a continuous band of plastic strain, and the upper soil mass creates a displacement. The depth of the shear band corresponded to the root zone depth.

Effects of the Planting Distance on the Stability of the Vegetated Slopes
The research object of case 2 was a vegetated slopes with a slope angle of 45 • , root depth of 1 m, root zone area of 2 m 2 , and planting location on the slope's surface. The effects of the root architecture on slope stability are discussed when the planting distances, a, are set to 5 m, 2.5 m, and 1.25 m, respectively. Figures 11 and 12 show the maximum horizontal displacement at the slope's toe after the failure of the vegetated slopes with the increase in SRF for different planting distances.  As shown in Figures 11 and 12, when the planting distance on the slope's surface was 2.5 m or 1.5 m, the slope stability was significantly improved compared with the bare slope (Figure 11b,c). However, the effect was not evident when the distance was 5 m (Figure 11a). When the planting distance is 5 m, the slope FOS is 1.3 for the four root architectures, which is equal to that of the bare slope. When the distance is 2.5 m, the FOS of the uniform and exponential root architectures is 1.5, and that of the triangular and parabolic root architectures is 1.4. The FOS of the four root architectures was 1.7 when the distance was 1.25 m.

Role of the Root Depth in the Stability of the Vegetated Slopes
The study object of case 3 is a vegetated slope with a slope angle of 45 • , planting distances of 2.5 m, root zone area of 2 m 2 , and planting location on the slope's surface. The effects of root architecture on the slope stability are discussed when the root depths, z r , are set to 0.5 m, 0.75 m, 1.0 m, 1.25 m, and 1.5 m, respectively. Figures 13 and 14 show the maximum horizontal displacement at the slope's toe after the failure of the vegetated slopes with an increase in SRF for various root depths. As shown in Figure 13, after slope failure, the horizontal displacement at the slope's toe with root depths of 0.5 m and 1.5 m is smaller than that with root depths of 0.75 m, 1.0 m, and 1.25 m. When the root depth is 1.5 m, the FOS of the triangular and exponential root architectures is 1.6, whereas that of the uniform and parabolic root architectures is greater than 1.8. The FOS of the four root architectures was 1.6 when the root depth was 0.5 m.

Effects of the Slope Angle on the Stability of the Vegetated Slopes
The research object of case 4 is a vegetated slope with a planting distance of 2.5 m, root depth of 1 m, root zone area of 2 m 2 , and planting location is the slope's surface. The influence of root architecture on slope stability is discussed when the slope angles, β, are set to 40 • , 45 • , and 50 • , respectively. Figures 15 and 16 show the maximum horizontal displacement at the slope's toe after the failure of vegetated slopes with an increase in the SRF for different slope angles.
As shown in Figures 15 and 16, when the slope angles are 40 • , 45 • , and 50 • , the FOS are 1.3, 1.5, and 1.7, respectively, for the uniform and exponential root architectures; for triangular root architecture, the FOS are 1.2, 1.4, and 1.6, and for parabolic root architecture, the FOS are 1.3, 1.4, and 1.6, respectively.

Influence of the Planting Location on the Stability of the Vegetated Slopes
The study object of case 5 is a vegetated slope with a slope angle of 45 • , planting distance of 2.5 m, root depth of 1 m, and root zone area of 2 m 2 . The effects of the root architecture on slope stability when the location of the root zone changes are discussed. We referred to Chok et al. [29] for the planting locations. Figure 17 shows the slopes with the planting positions of the slope's surface, slope's toe, the slope's surface and toe, upper slope region, lower slope region, upper and lower slope regions, and entire ground surface, as well as their final equivalent plastic strains during FOS = 1.6. The grey area represents the root zone (taking uniform root architecture as an example). Figures 18 and 19 show the maximum horizontal displacement at the slope's toe after slope instability with an increase in SRF for different planting locations.  Figure 17 shows that different planting positions have different influences on the shear band of the slope failure. The improvement in slope stability was the most obvious in the entire ground surface planting. Planting at the slope's toe can effectively reduce the sliding displacement of the slope soil.
As shown in Figure 16, vegetation on the entire ground surface (Figure 18g) had the best effect on the slope stability (FOS > 1.8). Planting on the lower slope region (Figure 18e) has little impact on the FOS, which is equal to that of the bare slope (FOS = 1.3). Planting on the slope's surface (Figure 18a) was better than that on the lower slope region, and the influence of uniform and exponential root architectures (FOS = 1.5) was better than that of triangular and parabolic root architectures (FOS = 1.4). The effects of planting on the slope's toe (Figure 18b) and the slope's surface and toe (Figure 18c) were both similar. Both were better than the slope's surface, and the impacts of the triangular and exponential root architectures (FOS = 1.7) were better than those of uniform and parabolic root architectures (FOS = 1.5). The effects of planting on the upper slope region (Figure 18d) and upper and lower slope regions (Figure 18f) were similar, and the influence of uniform architecture (FOS > 1.8) was better than that of the other three root architectures (FOS = 1.5). In Figure 19, the planting of uniform root architecture (Figure 19a) on the upper and lower slope regions, upper slope region, and entire ground surface (FOS > 1.8) are better than those on the slope's surface and toe, as well as on the slope's surface (FOS = 1.5). The planting effects of parabolic root architecture (Figure 19c) on the upper and lower slope regions, slope's surface, and upper slope region were weaker than those of the uniform root architecture. The effects of triangular ( Figure 19b) and exponential root architectures (Figure 19d) are similar, and they have a greater advantage when they are present at the slope's toe, as well as the slope's surface and toe (FOS = 1.7).

Discussion
In this study, the slip surface of the bare slope was mainly distributed in the weak layer at the interface of the rock and soil, whereas the slip surface of the vegetated slope was distributed at the bottom of the root zone, and the upper layer of soil entrains the root system to slip. The tensile strength and adhesion properties of the roots reinforced the soil. Plant roots with high tensile strength increase the confining stress of the soil through their compact root matrix system [29]. The results of the finite element analysis showed that reinforcement of the root can improve the stability of the slope, reduce the displacement of the landslide, and increase the FOS value.
The FOS value of the vegetated slope increased with decreasing planting distance. The higher the planting density, the stronger the root-strengthening effect. It should be noted that if the planting distance exceeds a certain range, then the slope stability may not improve. In this study, when the planting distance on the slope was 5 m, the effect was similar to that on the bare slope. The stability of a slope with a uniform root architecture is most sensitive to the planting distance.
For many slopes, the root depth is usually limited by bedrock, which is usually shallow and less than 2 m [58]. The failure depth of most slopes is between 0.5 m and 1 m, and the root zone plays a role in mechanical stability only when the root depth reaches the deeper soil layer [78]. With an increase in root depth, the FOS first decreased and then increased. When the root depth is 0.5 m, the slope's surface is covered with overlapping roots, which is similar to the geosynthetic reinforcement on the slope's surface. When the root depth was 1.5 m, it was similar to installing anti-slide piles or anchors. In addition, for uniform and parabolic root architectures, a greater root depth is more beneficial for slope stability. For triangular and exponential root architectures, shallow roots were more conducive. In addition, the stability of a slope with a uniform and parabolic root architecture is more sensitive to its root depth.
The slope FOS of the four root architectures decreased with an increase in slope angle. Regardless of whether the slope is steep or has a slight incline, uniform and exponential root architectures are more effective in improving the slope stability.
The position of the plants on the slope also affects their contribution to stability. Except for planting in the lower slope region, the stability of the vegetated slope in the other planting positions was better than that of the bare slope. If there is no restriction on planting location and vegetation quantity, planting on the entire ground surface of a slope has the best effect on improving slope stability, which is similar to the conclusion of [29]. Moreover, it is better to plant vegetation with a uniform root architecture in the upper slope region or plant vegetation with a triangular or exponential root architecture on the slope's toe. When the vegetated slope is unstable, the depth of the plastic strain at the top and toe of the slope is shallower; therefore, it is more conducive for roots to play a role. Tensile cracks may occur at the top of a slope, and the roots provide traction and support soil to prevent landslides [79]. The roots bear pressure at the toe of the slope, which can act as a support and inhibit soil sliding [80].

Conclusions
In this study, the eSPFEM was utilised to simulate the instability of vegetated slopes with large deformations. This method can reasonably predict the deformation process of the slope structure and the final deposition, avoiding the difficulty of numerical calculations and loss of calculation accuracy.
The Mohr-Coulomb constitutive model was extended by introducing the additional soil cohesion, c r , generated by roots and therefore increasing the shear strength of the soil. The boundary functions of the four root architectures (uniform, triangular, parabolic, and exponential) on the slope were derived. The FOS values of the four root architectures for various planting distances, root depths, slope angles, and planting locations were calculated using the shear strength reduction technique with a kinetic energy-based criterion, and the effects of root architecture on slope stability were evaluated.
Their results showed that roots can effectively improve slope stability and reduce landslide displacement. The higher the planting density, the stronger the root-strengthening effect. With an increase in root depth, the FOS first decreased and then increased. For uniform and parabolic root architectures, deeper roots are more beneficial to slope stability, whereas for triangular and exponential root architectures, shallower roots are more conducive. The FOS decreases with the increase in slope angle; uniform and exponential root architectures are more effective in improving slope stability, no matter whether the slope is steep or has a slight incline. Vegetation at the slope's toe can effectively reduce the sliding displacement of the slope soil. Planting on the entire ground surface had the best effect on slope stability improvement, followed by planting vegetation with uniform root architecture in the upper slope region or planting vegetation with triangular or exponential root architectures on the slope's toe.
This study provides valuable information on the contribution of different root architectures to slope stability and can guide the selection of vegetation species and planting locations, which can contribute to improving slope stability and optimising the management of mountain shelter forests. Nonetheless, the limitation of this model lies in the need for accurate parameterisation and a large amount of calculation, so it is currently only applicable to small slopes.

Data Availability Statement:
The data presented in this study are available upon request from the corresponding author.

Conflicts of Interest:
The authors declare no conflict of interest.