A critical evaluation of predictive models for rooted soil strength with application to predicting the seismic deformation of rooted slopes

This paper presents a comparative study of three different classes of model for estimating the reinforcing effect of plant roots in soil, namely (i) fibre pull-out model, (ii) fibre break models (including Wu and Waldron’s Model (WWM) and the Fibre Bundle Model (FBM)) and (iii) beam bending or p-y models (specifically Beam on a Non-linear Winkler-Foundation (BNWF) models). Firstly, the prediction model of root reinforcement based on pull-out being the dominant mechanism for different potential slip plane depths was proposed. The resulting root reinforcement calculated were then compared with those derived from the other two types of models. The estimated rooted soil strength distributions were then incorporated within a fully dynamic, plane-strain continuum finite element model to assess the consequences of the selection of rooted soil strength model on the global seismic stability of a vegetated slope (assessed via accumulated slip during earthquake shaking). For the particular case considered in this paper (no roots were observed to have broken after shearing), root cohesion predicted by the pull-out model is much closer to that the BNWF model, but is largely over-predicted by the family of fibre break models. In terms of the effects on the stability of vegetated slopes, there exists a threshold value beyond which the position of the critical slip plane would bypass the rooted zones, rather than passing through them. Further increase of root cohesion beyond this value has minimal effect on the global slope behaviour. This implies that significantly over-predicted root cohesion from fibre break models when used to model roots with non-negligible bending stiffness may still provide a reasonable prediction of overall behaviour, so long as the critical failure mechanism is already bypassing the root-reinforced zones.


Introduction
Understanding and quantifying the mechanical effect of vegetation on steep slopes began approximately 50 years ago with direct shear tests performed on soil blocks containing roots (Wu 1976(Wu , 2013Stokes et al. 2014). Since then, various approximate models for predicting root reinforcement of soils have been introduced. Generally, these models can be classified into two types: (i) continuum approaches, which consider the root-soil matrix as a homogenous material of increased strength Δτ (also described as root cohesion, c′ r ) or (ii) soil-root interaction approaches which consider roots as a structural element embedded in the soil.
Continuum approaches involve laboratory tests of or numerical simulations of representative elements of rooted soil, with the strength being represented as a Mohr-Coulomb failure envelope or yield surface, as conventionally used for non-vegetated soil, like fibre-reinforced sands (e.g. Michalowski and Čermák 2003;Zaimoglu and Yetimoglu 2012;Wood et al. 2016). This approach is convenient where the dimensions and spacing of the reinforcement are small and behaviour can be homogenised statistically; otherwise, tests are difficult to perform, time consuming and expensive. As for the latter approach, the soil-root interaction properties can be estimated from axial root properties which can be determined from axial tension or pull-out tests of the roots (e.g. Van Beek et al. 2005;Docker and Hubble 2008;Fan and Su 2008;Mickovski et al. 2009;Sonnenberg et al. 2010;Loades et al. 2010;Comino et al. 2010;Schwarz et al. 2011;Boldrin et al. 2017;). The additional resistance within the soil due to the presence of roots may then be introduced into stability calculations either as boundary forces (Greenwood et al. 2004;Greenwood 2006) or used to evaluate c′ r for use in the Mohr-Coulomb failure envelope equation (Waldron 1977;Wu et al. 1979;Pollen and Simon 2005).
For the latter approach, the most widely used soil-root interaction models are the family of fibre break models, commonly known as Wu and Waldron's Model (WWM, Wu 1976;Waldron 1977;Wu et al. 1979) and the Fibre Bundle Model (FBM, Pollen and Simon 2005). Both models assume that roots are highly flexible with negligible bending stiffness and will break (structurally) in tension during soil shear deformation, so that the additional strength provided by the roots is a function of root properties only (i.e. tensile strength of roots, root density and root orientations); however, an indirect effect of the soil properties is incorporated in the way these influence the growth of the roots and therefore the aforementioned root properties. The major difference between the two models lies in the ability of FBM to model progressive failure as the weakest roots within the root system break first (Thomas and Pollen-Bankhead 2010), with load shared between different diameters of root by either: (i) equal load applied to individual roots regardless of root dimension, (ii) load apportioned by root diameter or (iii) load apportioned by root cross-sectional area.
For plants with larger structural roots where root bending, rather than axial breakage, may be more dominant, considering the roots as flexible cables (fine roots) or bending beams (coarse/structural roots) subject to lateral loadings provides an alternative means of estimating root reinforcement, e.g. using py models, as reported by Duckett (2013), Liang et al. (2015) and Meijer et al. (2019). Such models use a set of transverse forcedisplacement (p-y) springs, which may be highly non-linear, to model the root-soil interaction in bending. They are Landslides Original Paper computationally efficient (at least compared to continuumbased finite element simulations) and implicitly incorporate the effects of soil properties as well (even where these may vary along the length of a root); however, further development would be required to generalise such analyses into analytical or finite difference-based models which are simple to use in practice.
For plants with shallower, fine and fibrous root systems, axial fibre breakage models may not always work (e.g. Loades et al. 2010) as roots are pulled out of the soil before breaking. Waldron and Dakessian (1981) reported a pull-out-based method for estimating root reinforcement; however, this model has not been widely adopted due to its dependence on root strain, which is relatively difficult to estimate in practice. Schwarz et al. (2010) proposed a more complicated pull-out-based model, the Root Bundle Model (RBM), which incorporated some features of the geometry (e.g. root length, root diameter, root branching pattern and root tortuosity) and mechanics (e.g. maximum tensile strength, Young's modulus, root-soil interfacial friction). Such a model contributes to understanding the pull-out behaviour of roots; however, it is seldom used in engineering practice or implemented in numerical codes since developed due to its complexity in input parameters (Wu et al. 2015).
In contrast, Wu (2013) proposed that the WWM (and by extension, FBM) could be adapted to use the pull-out capacity rather than the breakage tensile strength (whichever is lowest), to capture potential root pull-out. Further studies are, however, required to confirm this as little information is available regarding direct comparison between breakage strength and pull-out strength for different root species in different soil media (Schwarz et al. 2011;Sonnenberg et al. 2011;Kamchoom et al. 2014). There is some preliminary evidence that analytical models developed in the field of piling engineering may also be an efficient way estimating the pull-out capacity of roots .
This study will develop further the estimation of root reinforcement from pull-out capacity based on a series of laboratory pull-out tests of root analogue segments in sandy soil, in which the confining stresses are varied to simulate segments of root up to 1.5 m deep below ground level. The individual roots will be modelled as straight vertical analogue elements as typically made in prediction models (e.g. WWM and FEM). The root shear strength contributions so determined will then be compared with those derived from the different classes of fibre break model and a beam bending (p-y) model to evaluate the implications of selecting a particular method for quantifying root reinforcement for slope stability calculations. This will be achieved through incorporating the root reinforcements suggested by the different methods within a fully dynamic non-linear finite element (FE) model of a slope subject to earthquake-induced instability. Seismically induced slip is a convenient way of assessing the stability of a slope as instants during the earthquake where the factor of safety is less than 1.0 result in deformation which can be measured accurately, in contrast to trying to directly 'measure' the factor of safety of a slope under static conditions. The FE simulations are validated against a physical model test of the slope conducted in a geotechnical centrifuge and previously reported by Liang et al. (2015).

Laboratory testing of root analogue pull-out
Soil A uniformly graded fine sand (HST 95 Congleton silica sand) was used throughout this study, as this was previously used in the centrifuge slope test. Cohesionless soil (sand) is used as it is possible to pluviate this (dry) around the analogues, while also being an analogue for coarse-grained field soils that we have observed previously in field applications (e.g. Meijer et al. 2018). It is a specific fraction of the sand extracted at Bent farm, Congleton, Cheshire. The sand had a mean particle size D 50 = 0.16 mm, minimum dry bulk density ρ min = 1462 kg/m 3 , and maximum dry bulk density ρ max = 1795 kg/m 3 . This sand has been widely used in previous geotechnical research at the University of Dundee (e.g. Al-Defae et al. 2013;. For all tests, the sand was air pluviated dry, resulting in a density of ρ = 1636 ± 8 kg/m 3 and a relative density of 55-60%. The critical state friction angle of the sand is φ 0 cs =32°across a range of relative densities (9-93%) and effective confining stresses (5-200 kPa), as measured in direct shear tests (Al-Defae et al. 2013).

Model root analogues
Root analogues have previously often been made of either rubber or wood as contrasting analogue materials (e.g. Mickovski et al. 2007Mickovski et al. , 2010Sonnenberg et al. 2011) with material properties (strength and stiffness) which bracket typical mechanical properties of plant roots. Neither of these materials is ideal; however, Liang et al. (2015 pioneered the use of 3D printing to fabricate root analogues using Acrylonitrile Butadiene Styrene (ABS) plastic, which could exhibit representative mechanical behaviour to real roots of woody root species such as trees and shrubs, as shown in Fig. 1 (after Liang et al. 2015), based on uniaxial tensile testing within an Instron 4204 loading frame. The 3D printed analogues could be considered as a stack of fibres aligned unidirectionally, and such structure is very similar to the cellular structure of real roots, with overlying layers of tissue. Among them, the xylem layers, which consists of long, cylindrical cells that are joined from end to end and provide unidirectional fibre orientation, play a significant role in mechanical behaviour, driving the characterisation of tensile strength (Karam 2005). As a result, the analogues model tensile strength particularly well, including a dependence on diameter (power function, as commonly used to fit measured data of field roots), and while they are stiffer than field roots, they are a closer representation than either wood or rubber analogues. In Fig. 1, the 'real root' data was collated from the literature (Mora et al. 2009;Warren et al. 2009;Mickovski et al. 2009;Mao et al. 2012).
Throughout this study, straight vertical elements with 150-mm anchorage length were used to simulate individual roots. These rods, with diameters of 1.6 mm, 3 mm and 12 mm, represent roots 1.5 m long in the 1:10 scale centrifuge test summarised later. A steel hook was attached to the top of each root analogue used in the pull-out tests (either using epoxy-resin adhesive for smaller analogues or a screw in the larger ones) so that the Instron 5985 loading frame used to perform the pull-out tests could 'pick-up' the root with minimal disturbance using a horizontal bar.

Model preparation and test procedure
Vertical high-density polyethylene (HDPE) plastic tubes, of 150mm inner diameter and 500 mm in length, were used as model containers. Such dimensions of the tubes (D Pipe > 8D root ) ensures that any boundary effects of the container on the pull-out resistance are minimised according to previous suggestions for piles (Randolph 1981(Randolph , 2003. The soil was initially pluviated to a depth of 140 mm, at which point a single analogue was pushed vertically into the soil by 20 mm. Pluviation was then continued until the analogue was completely surrounded by sand. This was identical to the procedure used in the centrifuge test. In the centrifuge test, the root analogues represented roots 1.5 m deep, at which depth the confining vertical effective stress would be approximately 24 kPa. To simulate higher confining stresses in the 1:10 scale tests, slotted circular surcharge weights which just fitted within the tubes were added on the soil surface. The slot allowed both easy placement around the analogue after pluviation and also root pull-out through the gap while minimising stress non-uniformity over the surface area. Four levels of confining stresses (q = 0 kPa, 4 kPa, 8 kPa, 12 kPa) were considered in this study, which represented 0 m, 0.25 m, 0.50 m and 0.75 m below ground level. The maximum stress applied was Landslides limited by the number of weights which could be stacked within the tube.
The root analogues were pulled out at a speed of 10 mm/min using the aforementioned load frame (Instron 5985L7706, Instron Inc., UK). The capacity of the load cell used was 30 kN with an accuracy of 1mN.

Interpretation using the beta (β) method
Assuming a root analogue acts as a miniature pile with capacity provided from interface friction along the root 'shaft', the uplift capacity (F p ) of a segment of length (L) and constant diameter (D) will be given by: where K is the coefficient of lateral earth pressure, δ' is the friction angle mobilised at the root-soil interface, z is the depth, and σ 0 V z ð Þ is the vertical effective stress at depth z. In the tests, at the soil where γ is the unit weight of the soil. The beta method takes its name from the combination of the coefficient of lateral earth pressure and interface friction angle into a single parameter (K tan δ ′ = β) multiplying the vertical effective stress. From the pull-out tests, F p was measured and used to backcalculate the value of β representing that of a root segment at a depth of z = q/γ + L/2 ≈ q/γ. The overall resistance of a 1.5-m-long prototype root in the centrifuge could then be obtained using Eq.
(1) and the distribution of β with depth as measured from the tests.

Finite element modelling
General assumptions Finite element analyses were performed using the geotechnical engineering FE code PLAXIS 2D 2015. Two-dimensional planestrain dynamic analyses were conducted in the time domain in order to model the subsequent seismic response of a vegetated slope. The constitutive model used was the 'Hardening Soil model with small-strain stiffness' (HS small; Schanz et al. 1999), which can simulate the non-linear stress and strain dependent behaviour of geo-materials and also the limiting stiffness shown at very small strains (G 0 ). This specific constative model has previously been verified to be effective at simulating the dynamic behaviour of the HST95 sand used in this study (see Al-Defae et al. 2013). A summary of the adopted values of the model parameters for this sand at the density used in both the pull-out tests and the centrifuge test are listed in Table 1.
The FE mesh employed is shown in Fig. 2. The slope is 2.4 m high from toe to crest, with a further 0.8 m of soil underneath, which represents one side of a small height embankment, such as might support road or rail infrastructure at the crest. The slope angle is 27°, representing a 1:2 slope. The soil is modelled with 15noded triangular elements with 12 Gaussian points.
The boundary conditions (see Fig. 2b) were modelled as an extension of both the left and right boundaries of centrifuge model (see Fig. 2a) to represent a semi-infinite soil condition provided by the Equivalent Shear Beam (ESB) container within centrifuge tests (See Zeng and Schofield 1996). The boundary conditions along the bottom boundary of the mesh are fixity in the vertical direction and prescribed values of acceleration as a function of time in the horizontal direction, while viscous boundaries, which allow lateral deformation in reaction to normal stress and incorporate nonreflecting elements, were applied along the vertical sides of the model in both directions based on the method described by Lysmer and Kuhlemeyer (1969).
The input accelerations consisted of eight successive earthquake motions, from three different historical events with distinct peak ground acceleration (PGA), duration and frequency content, as simulated in the centrifuge test. The sequence of motions is summarised in Table 2. Further details about these motions can be found in Liang et al. (2015).

Root-soil matrix modelling
Two types of physical models were used to represent three types of natural tree root systems in the centrifuge tests (see Liang et al. 2015. Among them, plate/heart root system was modelled as an idealised group of straight vertical rods (see Fig. 3b), and this may be a reasonable simplification for such two types of root system, where the main vertical roots which cross the shear plane grow downwards from the main horizontal lateral roots. This idealisation is used in existing analytical models described earlier. However, this idealisation may not be suitable to simulate a tap root system where lateral roots are interlocked by the main tap roots. For this specific root system, a 3D root cluster model (see Fig. 3a) was also used. It should be noted that the straight vertical rods group model (see Fig. 3b) was arranged to have the same cross-sectional distribution at the level of the middle of the 3D root cluster model, which was designed to identify the corresponding root morphology effect (after ). In the following section, the idealised group of straight vertical rods will be used to assess the suitability of different models for root reinforcement and the 3D root cluster model will be shown as a reference for comparison. Figure 4 shows a comparison of measured root reinforcement (in terms of additional soil shear strength due to roots, c′ r ) riginal Paper Landslides provided by the 3D printed root analogues used in this study with some in situ Direct Shear Apparatus (DSA) test data on 14 young trees and 1 shrub (data collected by . In this figure, the data of the additional shear strength provided by both the straight vertical rods group model and the 3D root cluster model were obtained from large laboratory DSA tests across different shear planes and confining normal stresses (Liang et al. 2015. As can be seen in Fig. 4, most of the in situ test data available were concentrated on the top 0.5 m due to the limitations of the shear apparatus. In addition, measured root reinforcement data deeper than 1 m for mature trees (S. Mexicana, E. camaldulensis and M. ericifolia, after Shields and Gray 1992; Abernethy and Rutherfurd 2001) is also shown in Fig. 4 for comparison. These data points were obtained by using a root reinforcement model (WWM) according to the field data of root distribution and root tensile strength.
The comparison clearly indicates that the printed root analogues have root reinforcement highly comparable with field roots. It should be noted here that the straight vertical rods group model demonstrated generally higher magnitude of root reinforcement with depth, compared to the 3D root cluster model. This is not unexpected as the individual roots of the straight vertical rods group model had higher anchorage length in the lower layer of soil.
The spatial distribution of root groupings used in this study is shown in Fig. 2, representing the areas rooted in the centrifuge test (see later). The root-soil matrix was modelled using a composite set of soil blocks (see Fig. 2) with a distinct additional soil shear strength due to roots c′ r , added to the HST95 soil properties in these zones (e.g. Li et al. 2016, Temgoua et al. 2016). This parameter was determined for the case of various fibre break and bundle models (WWM and FBM), fibre pull-out model using the results of the previously described laboratory tests, and a recently published beam bending model (Liang et al. 2015). It should be mentioned here that the determination of the size of the zone of the root group influence (the width of the blocks in Fig. 2) is essential for accurately predicting the global slope performance. Liang et al. (2015) suggested to use the  Landslides actual extreme boundary circumscribing the root analogues (also known as the critical rooted zone) and it was employed in the modelling (Fig. 3b). Roots can be stretched by 10-20% of their length before failure while most soils fail at strains around 2% (Tobin et al. 2007). When a root-soil system is subjected to shear loading, the soil will typically reach peak strength ahead of the roots, with the roots then providing enhanced shear strength until they themselves fail. The value of c′ r used for the WWM simulations described here assumes that all of the roots' strengths are mobilised and broken simultaneously with the contribution of the roots to shear strength being: where T rn is the tensile strength of a root of size (n), RAR n is the root area ratio (root cross-sectional area as a fraction of the crosssectional area of the critical rooted zone) of roots of this size, and R θ is a root orientation factor, calculated using where φ ' is the effective stress friction angle of the soil, generally taken as the critical friction angle (see Table 1) in practice (e.g. Wu 1976;Wu et al. 1979;Docker and Hubble 2008) and θ is the angle of shear distortion of the (vertical) root, estimated by where x is the shear displacement at failure (peak shear resistance), and Z is the thickness of the shear zone (Fig. 5a). Wu (1976) found that R θ was fairly insensitive to normal variation in θ and φ ' (40-70°and 25-40°, respectively), with values ranging from 0.92 to 1.31. Hence, a constant value of 1.2 was used in practice to replace R θ .
In the FBM, c r ' is also a function of RAR and tensile strength, but the roots are considered to break progressively rather than simultaneously. When some roots break, the total shear force is redistributed among the remaining roots with this being apportioned to each root in a bundle according to one of three potential assumptions (after Pollen and Simon 2005): where n is the root number ordered from strongest to weakest, n ∈ [1, N]; N is the total number of roots; j is the weakest root removed at each simulation step, j ∈ [1, N]; and T rj is the tensile strength of the weakest remaining root. For the above three assumptions, the breaking order of each root can be evaluated by T rj , T rj D j and T rj D j 2 , as a function of root CSA, root diameter and root number, respectively (Mao et al. 2012).
The WWM and FBM assume that root breakage represents the maximum tension a given root can carry; however, in reality, some roots will slip through the soil before breakage if the interface shear strength between the root and soil is low. A root may either pull up from the stable subsoil below the shear plane, or pull downwards through the slipping soil mass, depending on the depth of slip surface (i.e. the length of root in each part of the soil). The root is considered as a series of segments such that for any given element (Fig. 5(b)) below the shear plane, the tensile force and shaft resistance are in equilibrium: where T i-1 is the tensile stress generated at the top section of ith element, T i is the tensile stress generated at the bottom section of ith element, l i is the length of ith element and F i is the shaft resistance provided to the ith element by the surrounding soil. The maximum tensile stress T up is generated at the slip plane, while tensile stress reduces to zero at the tip of root. Integration of the stresses on these elements gives: where T up is the maximum tensile stress within the root. In the same way, the generated maximum tensile stress can also be estimated through the integration of root elements above the slip surface: Equation (2) can then be modified for determination of the root cohesion: where T n is the minimum tensile stress generated within the roots: If T n = T up , the part of the root below the slip plane will be pulled up from the soil below; if T n = T down , the part of the root within the slipping mass will be pulled down from the above the slip plane; if T n = T r , the pull-out strength is high enough that root will break.

Results and discussion
Pull-out of model root analogues Figure 6 shows selected pull-out resistance-displacement relationships measured for the model root analogues at varied confining stress. The pull-out resistance experienced a rapid increase initially and reached a peak at ≈ 2 mm displacement. After mobilisation of peak resistance, the pull-out resistance reduced steadily due to a combination of strain softening at the soil-root interface and the reduction of interface area as the root is pulled out, as suggested by Mickovski et al. (2010). The behaviour is very similar to wood and rubber analogues pulled out of similar soil reported by Mickovski et al. (2007), and it was also very similar to that of an unbranched primary root of Pea (Pisum sativum L.) (Hamza et al. 2004) and a tap root of Sunflower (Helianthus annuus L.) (Ennos 1989). Figure 7a shows that the maximum pull-out resistance of roots increased with diameter for all confining stresses. The maximum pull-out resistance was not linearly proportional to the diameter (or surface area). This indicates that the root-soil interface shear Landslides strength varied with diameter. Mickovski et al. (2010) treated this as a type of scale effect between root dimeter and soil particle size based on the standard axial pile and soil nailing theory. As suggested by Meyerhof (1983) and Foray et al. (1998), the interface friction between piles and soil depends on the ratio of the interface shear band thickness to the pile diameter while the shear band thickness is related to the soil particle size and asperity height on the root surface. For a given confining stress and root material, the shear band thickness is constant for all diameters, such that roots with smaller diameter experienced a higher interface friction. Figure 7 also indicates that the maximum pull-out resistance of the root analogues increased with confining stress, which was true for all diameters tested. The variation of peak pull-out resistance against confining stress is shown in Fig. 7b and also exhibited nonlinear behaviour. This will be further discussed in the following section.
The back-calculated β value for pull-out tests using Eq. (1) is shown in Fig. 8. For a given diameter, two distinct values were presented depending on the effective embedded depth of root (corresponding to the confining stress level during the tests). For stresses representing shallow embedment (up to 0.15 m below the ground surface), the calculated β value was approximately 3 times larger than those for stresses representing elements of the root deeper than 0.15 m. This may be attributable to high dilation at very low stress. The interface friction angle δ' is a function both of root roughness and soil properties. API codes for piling (API 2000) recommend to estimate it as a function of soil peak friction angle ϕ', that is where k is a dimensionless coefficient to account for root roughness and soil particle size. For a given diameter in this study, the k value should be constant. Bolton (1986) proposed a model to quantify the effects of dilation at low confining stress in terms of an increase in the soil friction angle ϕ' pk , as shown below: where A is the dimensionless factor to account for strain type: A = 3 for triaxial strain; A = 5 for plane strain. I R is given by where I D is the relative density of the sand, σ'(z) is the mean confining stress in the soil at a known depth z; Q and R are fitting parameters that depend on the intrinsic sand characteristics, which can be simplified to 10 and 1 respectively when 0 < I R < 4 (Bolton 1986), while at very low confining stress level (I R > 4), Chakraborty and Salgado (2010) suggested to use7.1 + 0.75 ln p ′ (for plane strain) and 1, respectively (after Liang and Knappett 2017b). The lateral earth pressure coefficient is also incorporated within β. In this study, the root-soil model was prepared using air pluviation around the root analogues so that the lateral earth pressure coefficient K is initially approximately equal to the undistributed lateral earth pressure coefficient K 0 : During uprooting, the movement of the root causes the surrounding soil to dilate and also experience lateral compression strain; as a result, K increases from K 0 until it reaches the limiting passive state, K p (see Fig. 9, after Knappett and Craig 2012):  Fig. 6 Results of pull-out tests on root segments with varied confining stress: a 1.6-mm-diameter roots; (b) 3-mm-diameter roots; (c) 12-mm-diameter roots The soil lateral earth pressure coefficient should fall between K p and K 0 , and an upper bound of K p and lower bound of K 0 were considered.
The β value expected was calculated as a function of depth using Eqs. (14)-(18) and root roughness coefficient k between 0 and 1 as shown in Fig. 10. When k is close to 1 (fully rough interface), β experiences a significant increase near the ground surface and is relatively constant at depth below 0.25 m, which explains the measured beta values in Fig. 8.

Prediction of rooted soil strength
The predicted pull-out resistance (POR) of the 1.5-m-long root analogues given any possible slip plane position with depth was determined using the measured β values from the previous section (1), and compared with the breakage force so that c r ' for the group of roots allowing pull-out could subsequently be found. As shown in Fig. 11, all three types of roots will be directly pulled out without any breakage regardless of any potential slip plane depths. However, for natural roots, any branches or tortuosity may increase the pull-out resistance of single root (e.g. Mickovski et al. 2007, Schwarz et al. 2010. Once the pullout resistance is high enough to exceed the break force, roots will break before pull-out. The derived root cohesion considering the shear plane at different depths (with change of confining stress) was then determined using Eqs. (12) and (13) and the data in Fig. 11 and validated  Landslides against the measured values. This is shown in Fig. 12 with a comparison to the calculated root cohesions determined using fibre break models (WWM and FBM) and a previously proposed Beam on-non-linear-Winkler-foundation (BNWF) model presented by Liang et al. (2015). The BNWF model appears to demonstrate the best match with the limited DSA test data points. In contrast, using the root cohesion distribution assuming pull out of roots under-predicts the additional resistance contributed by the root analogues, while fibre break models (WWM and FBM) significantly over-predict the additional resistance. The corresponding failure mechanisms of these three models are shown schematically in Fig. 13. It should be noted that following shearing, no roots were observed to fail when roots are subjected to a shear displacement of around 50 mm (see Liang et al. 2015 and the measured angles of shear distortion θ, ranged from 30°to 60°, with resulting R θ of 1.04-1.18; hence, the value of 1.2 recommend by Wu et al. (1979) could provide a reasonable representation of R θ and was used in the calculation. The root with a diameter > 10 mm (central root in Fig. 3b) was not included in the calculation for WWM or FBM considering their rotational behaviour during soil slip (Genet et al. 2008;Stokes et al. 2009;Mao et al. 2012). Despite not including 12-mm-diameter root, WWM and FBM models appear to significantly over-estimate the additional strength contributed by the root analogues. This overestimation is not surprising as no roots were observed to break. The implications of this wide range of predicted root strength contributions was assessed by simulating each of the distributions shown in Fig. 12 within the FE simulations.
Finite element simulations Figure 14 shows a comparison of measured (from centrifuge tests) and simulated (from numerical simulation) settlement at the crest of the model slope under the eight successive earthquake motions. The presence of root analogues caused a significant reduction in permanent slope movement compared to the fallow case, regardless of root morphology. Interestingly, both the centrifuge tests and numerical simulation results consistently suggest that slopes reinforced by the 3D root cluster model induced a greater reduction of the slope crest settlement than the one reinforced by a group of straight vertical rods (see Fig. 14). This is not consistent with the variation of root cohesion profile shown in Fig. 4, for which the straight root case is generally higher. It may therefore not be suitable to relate the root reinforcing effects solely to root cohesion. Apart from root cohesion, another key parameter that may affect the root reinforcing effect is the lateral extent of the root system, namely, the diameter of critical rooted zone. As shown in Fig. 3, the diameter of the critical rooted zone of the straight root group (0.4 m) is much smaller than that of the 3D root cluster model (1.0 m). Thus, more attention should be paid on the importance of the lateral extent of root systems on the overall performance of a vegetated slope in engineering practice. Further parametric studies are needed to more fully quantify the influence of the diameter of critical rooted zone.
The implication of using commonly made model simplifications is examined in Fig. 14 with a direct comparison with the test result for slope reinforced with straight root group. Using the root cohesion distribution determined from the bending (BNWF) model presented by Liang et al. (2015) provides the best match to the measured response in the centrifuge test. In contrast, using the root cohesion distribution assuming uplifting of roots overpredicts the crest settlement compared with the measured value, as this model does not derive any additional resistance from the bending of the root analogues. Interestingly, when using the significantly over-predicted strength derived from WWM (see Fig.  12), a good match is obtained. This indicates that once the root contribution to shear strength is sufficiently high, continuing increases in the magnitude of root reinforcement may not have a significant effect on subsequent response. However, the key issue then becomes the determination of this critical 'minimum' shear strength, for which it is essential to be able to model root-soil interaction correctly. Investigating the failure mechanism of the slope, Fig. 15 shows the mobilised shear strain inside the slope after the earthquake sequence. The presence of roots causes the slip plane to move depending on the root contribution to shear strength. At low root cohesion (Fig. 15b), the increased cohesion is small enough that the critical mechanism passes through the rooted zone (as is conventionally assumed), though the additional shear strength within the rooted zone decreases the total shear

Original Paper
Landslides strain mobilised. When the rooted soil strength is high enough to buttress the movement of the sliding block, the slope is separated into several small sliding blocks, as shown in Fig. 15c and this improves the overall performance of the slope. Further increase of the root cohesion has minimal effect on slope deformation (Fig.  15d) as the soil is already failing in the unreinforced area between the zones reinforced with root analogues.

Implications for engineering application
A key finding of this study is that where roots are clustered on slopes (e.g. beneath individual shrubs or trees), there exists a threshold shear strength distribution beyond which increasing the strength of the rooted zone (e.g. by stronger or more roots) will not provide any further benefit to stability as the critical failure mechanism will already have moved to bypass the stronger zones and fail through the weaker unreinforced zones. It is therefore important that root strength contributions are not simply added to the resistance of the fallow shear plane, with the implication that indefinitely increasing c r ' will result in ever greater stability. Such a finding is important for the selection of suitable species to protect slopes against natural hazards.

Original Paper
Landslides experiments presented here, use of the BNWF model to determine rooted soil strength appeared to result in the best predictions of the slope behaviour. The pull-out model under-predicted reinforcement as the root must deform in shear to generate the axial strains to reach pull-out strength. However, the model is then neglecting the resistance due to bending of the roots, which results in an under-prediction of the reinforcing effect on the slope. The WWM and FBM models over-predicted reinforcement as they similarly rely on the roots to deform to generate breakage tensile stresses; however, because these stresses are very high in these structurally competent root analogues, the bending resistance of the analogues means that the soil will fail in shear around the analogues before this sufficient tensile stress can be achieved (a feature which is captured in the BNWF model via the limiting strength of the soil springs). For woody root species such as trees and shrubs, structural and coarse roots occupy the majority of the total root mass; for example, a study by Parr and Cameron (2004) observed that a spruce tree had a total of 82,500 roots. Among them, coarse roots (> 5 mm) comprised 62%.
FBM/pull-out models may continue to be suitable for modelling fibrous root systems where the roots are of such small diameter that they have negligible resistance in bending. The smallest analogues used in this study were 1.6-mm diameter and none of these were observed to break within the centrifuge test under substantial shear deformation. These analogues are notably much finer than the critical diameter (10 mm) used by Mao et al. (2012) to predict the root cohesion of natural diverse mountain forests via a fibre break model. Further study is required to determine the combination of root diameter and soil conditions below which the bending resistance is negligible and fibre break models can be employed.
It was also observed in this study that the lateral extent of the root system also plays an important role, sometimes even more important than the magnitude of root cohesion, on the improvement of the overall performance of vegetated slopes (see Fig. 14 and detailed discussion in the "Finite element simulations" section). This has important implications for the engineering use of vegetation to protect natural or man-made slopes against landslides, and it may be advisable to select plant species for their propensity of lateral spread and deep rooting, rather than species with the strongest possible roots, so as to maximise the reinforcing effects given by the vegetation. For the field investigation of vegetation sites, it would be desirable to use hand-held devices for rapid testing of root strength properties in multiple locations quickly (e.g. Meijer et al. 2016Meijer et al. , 2018 to better quantify the distribution of rooted soil strength with position, rather than only a limited number of highly detailed yet slow and expensive in situ tests (e.g. large-size direct shear box tests).

Conclusions
The family of fibre break models (i.e. both WWM and FBM) predicted much higher root cohesion than the pull-out models proposed in this study. In all of the experiments, root analogues with diameters ranging from 1.6 to 12 mm did not show any Fig. 15 Comparison of failure mechanism between fallow slope and root-reinforced slope with cohesion derived from different analytical models: a fallow slope; b rooted slope with cohesion from fibre pull-out model; c rooted slope with cohesion from beam bending model; d rooted slope with cohesion from fibre break model (WWM) Landslides breakage upon vertical pull out under confining pressures ranging between 0 and 12 kPa. However, the pull-out models underpredicted resistance compared to a BNWF (root bending) model, which demonstrated the best fit to the test results. When the root cohesion estimated by the different predictive models was applied in a dynamic slope analysis, there appears to exist a threshold value of enhanced rooted soil shear strength in concentrated zones around the root analogues where the position of the critical slip plane would bypass the rooted zones, rather than passing through them. Interestingly, further increase of root cohesion beyond this threshold value introduced very limited influence on the global slope behaviour. This implies that the significantly over-predicted root cohesion made by the family of fibre break models for roots that have non-negligible bending stiffness may still provide a reasonable prediction of the overall behaviour, so long as the critical failure mechanism is already bypassing the rootreinforced zones. However, they may potentially over-estimate the stability improvement if the actual rooted soil strength is low.

Acknowledgements
The James Hutton Institute receives funding from the Scottish Government (Rural & Environmental Services & Analytical Services Division). The authors would also like to express their sincere gratitude to Dr. Gary Callon at the University of Dundee for his assistance in printing the model root analogues and undertaking the test programme. The first author would like to acknowledge the financial support of the China Scholarship Council.

Notation
A, dimensionless factor accounting for strain type; c ′ , cohesion of soil; c' r , additional soil shear strength due to roots; D, diameter; D j , diameter of the weakest remaining root; D n , diameter of root size n; D pipe , inner diameter of plastic tube; D root , diameter of root; D 50 , particle diameter at which 50% is smaller; E ref 50 , triaxial secant stiffness (at 50% of deviatoric failure stress in drained triaxial compression; E ref oed , oedometric tangent stiffness (in compression); E ref ur , unloading-reloading stiffness; F down , pull down force; F i , shaft resistance provided to the ith element; F p , uplift capacity; F up , pull up force; g, acceleration due to gravity (= 9.81 m/s 2 ); G 0 , small strain shear modulus; E ref 0 , reference small strain shear modulus; i, number of element; I D , relative density; I R , relative dilation index; j, weakest root removed at each simulation; k, dimensionless coefficient accounting for root roughness and soil particle size; K, coefficient of earth press at rest; K 0 , coefficient of undistributed earth press at rest; K a , coefficient of active earth press at rest; K p , coefficient of passive earth press at rest; l i , length of ith element; L, length of pile; m ′ , power-law index for stress level; n, root number ordered from strongest to weakest; N, total root number; p, reaction from soil due to the deflection of pile; q, bearing pressure; Q, fitting parameter for relative dilation index; R, fitting parameter for relative dilation index; R f , ratio of deviatoric failure stress to asymptotic limiting deviator stress; R θ , root orientation factor; RAR, root area ratio; RAR n , root area ratio of root size n; T down , tensile stress below the shear plane within the root; T i-1 , tensile stress generated at the top section of i th element; T i , tensile stress generated at the bottom section of i th element; T up , tensile stress above the shear plane within the root; T r , ultimate tensile strength; T rj , tensile strength of weakest remaining root; T rn , tensile strength of a root of size n; T n , maximum tensile stress within the root; x, shear displacement at failure; y, deflection; z, depth of soil; Z, thickness of the shear zone; β, drained interface strength parameter; δ'', root-soil interface friction angle; ε s, 0.7 , shear strain; ρ max , maximum dry bulk density, ρ min , minimum dry bulk density; γ, unit weight; γ unsat , unsaturated unit weight; γ sat , saturated unit weight; θ, angle of shear distortion; Δτ, additional soil shear strength due to reinforcement; σ′, effective confining stress; σ' v , vertical effective stress; ϕ ′ , effective angle of friction; ϕ' cs , critical state angle of friction; ϕ' pk , (secant) peak angel of friction.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http:// creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.