Small-scale modelling of plant root systems using 3D printing, with applications to investigate the role of vegetation on earthquake-induced landslides

Abstract

Vegetation has been previously proposed as a method for protecting artificial and natural slopes against shallow landslides (e.g. as may be triggered by an earthquake); however, previous research has concentrated on individual root soil interaction during shear deformation rather than the global slope behaviour due to the extreme expense and difficulty involved in conducting full-scale field tests. Geotechnical centrifuge modelling offers an opportunity to investigate in detail the engineering performance of vegetated slopes, but its application has been restricted due to the lack of availability of suitable root analogues that can repeatably replicate appropriate mechanical properties (stiffness and strength) and realistic 3D geometry. This study employed 3D printing to develop a representative and repeatable 1:10 scale model of a tree root cluster (representing roots up to 1.5 m deep at prototype scale) that can be used within a geotechnical centrifuge to investigate the response of a vegetated slope subject to earthquake ground motion. The printed acrylonitrile butadiene styrene (ABS) plastic root model was identified to be highly representative of the geometry and mechanical behaviour (stiffness and strength) of real woody root systems. A programme of large direct shear tests was also performed to evaluate the additional strength provided by the root analogues within soil that is slipping and investigate the influence of various characteristics (including root area ratio (RAR), soil confining effective stress and root morphology) on this reinforcing effect. Our results show that root reinforcement is not only a function of root mechanical properties but also depends on factors including surrounding effective confining stress (resulting in depth dependency even for the same RAR), depth of the slip plane and root morphology. When subject to shear loading in soil, the tap root appeared to structurally transfer load within the root system, including to smaller and deeper roots which subsequently broke or were pulled out. Finally, the root analogues were added to model slopes subjected to earthquake ground motion in the centrifuge, where it was revealed that vegetation can substantially reduce earthquake-induced slope deformation in the soil conditions tested (76% reduction on crest permanent settlement during slippage). Both the realistic 3D geometry and highly simplified root morphologies, as characterised mechanically by the shear tests, were tested in the centrifuge which, despite exhibiting very different levels of additional strength in the shear tests, resulted in very similar responses of the slopes. This suggests that once a certain minimum level of reinforcement has been reached which will alter the deformation mechanism within the slope, further increases of root contribution (e.g. due to differences in root morphology) do not have a large further effect on improving slope stability.

Introduction

Landslides induced by long and/or intense rainfall or earthquake events have significant effects on lives and infrastructure in many parts of the world (Petley 2012). As an example, in the 2008 Wenchuan earthquake, 69,227 lives were lost and 374,643 people injured, with a further 17,923 listed as missing. During this event, tens of thousands of landslides were triggered over a broad area, some of which buried large sections of some towns and blocked transport links and dammed rivers (Dai et al. 2011). It was estimated that the total losses exceeded £80 billion, and the losses from the earthquake-triggered landslides accounted for over a third of the total earthquake losses (Chen et al. 2008).

Many types of traditional geotechnical methods have been used to improve the slope stability and reduce landslides, such as soil nailing, piles and retaining walls. Compared with those traditional methods, vegetation is an effective and more environmental-friendly approach and has been incorporated into engineering practice (Stokes et al. 2014). To investigate the global performance of rooted slopes and verify the contribution of roots to the behaviour of slopes, some trials (e.g. Smethurst et al. 2006; Leung and Ng 2013) have previously been conducted in the field or on 1:1 scaled slope models. An increased occurrence of shallow landslides has been observed after deforestation of natural forested slopes (Preti 2013; Vergani et al. 2014), such as in southeast Alaska where it was reported that the frequency of landslides increased by 3.8 times after a large-scale decline of yellow cedar (Johnson and Wilcock 2002). However, such large trials are expensive and time consuming, and therefore relatively rare, and have not considered earthquakes as a trigger due to the additional time dependency of a large earthquake occurring while the slope is actively instrumented. An alternative approach which has previously been followed has been to collect information on the root properties (e.g. root tensile strength, root architecture, root cohesion) from stable slopes in situ and perform back calculations of slope behaviour employing existing analytical models or computational models (e.g. Danjon et al. 2008; Mao et al. 2014). This approach highly depends on the accuracy of the analytical models or soil constitutive models selected (Wu 2013). While this has indicated some of the characteristics that impact the global performance of rooted slope, there is a limit to the information that can be gathered from only stable slopes.

Geotechnical centrifuge modelling is an approach which can simulate the global performance of a full-scale soil slope prototype to a high level of fidelity, by achieving similitude of stresses at homologous points within the model and prototype. This can address all of the aforementioned limitations of large trials, being low cost and allowing simulation of external actions (e.g. loading, rainfall, earthquakes) on models with carefully controlled and well-known properties on demand. Despite the popularity of this technique in addressing many problems in geotechnical engineering, only relatively few centrifuge model tests have been performed on vegetation-reinforced slopes (e.g. Sonnenberg et al. 2010, 2011; Ng et al. 2016), and none of these has considered earthquake triggering. In these previous studies, the presence of roots was observed to transform the failure mode of the slopes from progressive block failure to translational failure, which highly improved the stability of the natural slope. The drawback of these studies however lies in the root analogues which were used to model real roots. Straight rods, occasionally with some highly simplified branching patterns following the procedure introduced by Mickovski et al. (2007), were employed to represent highly simplified tap root, herringbone pattern and dichotomous pattern systems. Such simplified structures may mask the effects of more complicated (and realistic) morphologies (Mickovski and Van Beek 2009; Fourcaud et al. 2008a,b; Ghestem et al. 2013). Additionally, the analogues used were typically made of either rubber or wood as materials with low and high stiffness and strength, respectively. While these materials bracket the stiffness and strength of typical roots, neither of their properties is ideal (Liang et al. 2015; Liang and Knappett 2017).

The principal aim of this paper is to present an improved tree root model that can be used repeatably in centrifuge model tests to more realistically represent the geometry, spatial distribution and mechanical properties of tree roots. Root system traits that may influence the mechanical behaviour of tree root systems and the design of the model roots will be reviewed first, before presenting a new modelling procedure that uses 3D printing to repeatably model many of the key properties of roots (in this case, at 1:10 scale). The suitability of the model roots will be identified through performing axial tensile testing on straight rods fabricated using the same method. Large direct shear apparatus (DSA) tests will then be conducted with these root models in fine sand to measure the root contribution to shear strength of the rooted soil composite. Finally, the model roots will be used in centrifuge tests of slopes subjected to earthquake ground motion as a way of inducing soil slip to demonstrate the potential impact of tree roots on seismic slope performance. The results will be discussed in relation to the soil–root interaction behaviour observed in the large DSA tests.

Major concerns in root physical modelling

A body of reviews (e.g. Tobin et al. 2007; Danjon and Reubens 2008) have been conducted on the properties of root systems. Here, the major concerns related to the mechanical behaviour of root systems will be synthesised and updated within the context of developing a new 3D analogue for tree root systems.

Root architecture

Root system architecture (RSA) is a result of a sequence of formation processes including branching, elongation, gravitropic response, thickening and turnover (Thaler and Pagès 1998). Root morphology and root area ratio have been considered to contribute more than root mechanical traits to the additional root cohesion (Docker and Hubble 2008; Fan and Chen 2010; Ghestem et al. 2013).

The lateral extent of root systems (especially of trees) can be considerable, with the majority of them spreading by several meters. Measurements of root distributions around individual isolated trees have revealed that root biomass generally decreases with increasing distance from the stem (Schenk and Jackson 2002a; Göttlicher et al. 2008). It is therefore important to consider root systems as a function of influence zone. Two types of root influence zone have previously been defined, namely critical root zone (CRZ) and zone of rapid taper (ZRT). CRZ extends to approximately one third to half of the zone the roots really occupy (known as the root protection zone (RPZ)). The critical root radius is determined relative to the diameter at breast height (DBH). For each 1 cm of DBH, 18 cm of critical root radius is expected for sensitive, older or unhealthy trees or 12 cm for tolerant, younger, healthy trees (Johnson 1999). The ZRT is more closely related to the mechanical behaviour of roots as this concept defines the zone of dominant structural roots, which have been found to provide more than 80% of the total root mass (Coutts et al. 1999). It is evident that most coniferous trees are supported by 3 to 11 large structural roots (Coutts et al. 1999; Mickovski and Ennos 2003; Tobin et al. 2007). Lateral roots taper rapidly until a maintained diameter of 10–20 mm, by which stage they have lost much of their rigidity and physical strength. Beyond the ZRT, lateral roots extend outwards in a broad zone for many meters, without further decrease in size. The volume of sinker roots is mainly located within twice the DBH radial distance according to Danjon et al. (2005), which provides a quantitative indication for physical modelling.

Plant root systems can penetrate to depths of several meters if unconstrained by soil conditions (Canadell et al. 1996; Jackson et al. 1996; Schenk and Jackson 2002b). Most species have roots at least of 2-m length (Schenk and Jackson 2002b). Despite rooting so deep, the majority of root biomass is concentrated in the top soil and typically shows exponential reduction of root density with increasing depth. An asymptotic equation introduced by Gale and Grigal (1987) has been used by most authors to model the cumulative root fraction with regard to depth for most species:

$$ Y=1-{\beta}^{100 d} $$

(1)

where Y is the cumulative root fraction from the surface to depth d (m) and β is the depth coefficient, which depends upon vegetation type. High values of β correspond to a greater proportion of roots with depth. Jackson et al. (1996) compiled a database of 250 root studies and found that β values for tree, shrub and grass roots across the major biomes are 0.97, 0.978 and 0.952, respectively.

Root mechanical properties

Tree root systems have evolved to an optimum balance of stiffness and strength (Bischetti et al. 2005). The response of structural roots under lateral loading (such as from wind loading or self-weight on slopes) conforms to beam theory. The elastic flexural stiffness of the beam is EI, where E is the Young’s modulus of the material and demonstrates a declining power law trend with diameter, i.e. E ∼ D ^α, where α is typically −0.5 to −2 (Mickovski et al. 2009). I is the second moment of area; if the beam is circular in cross section (of diameter D), I can be calculated using

$$ I=\frac{\pi {D}^4}{64} $$

(2)

The ultimate tensile strength (UTS) of roots is also normally considered to be a function of root diameter

$$ {T}_r={kD}^m $$

(3)

From Genet et al. (2005), the value of k and the exponent m varied within ranges of 23 to 64 and −1.0 to −0.5, respectively, for D in millimetre and T _r in megapascal. There is a trend for angiosperms to have an exponent value near −1.0 and for gymnosperms to have an exponent value near −0.75 (Bischetti et al. 2005); however, it should be noted that many exceptions exist.

Roots can stretch 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.

Development of a 1:10 scale model

Overview of modelling technique

The 3D root models considered in this paper were designed following the procedure shown in Fig. 1. The initial architecture was based on the tap root system of a white oak tree located at the Warnell School for Forestry and Natural Resources, University of Georgia (after Danjon et al. 2008) for which detailed characterisation of the root system was undertaken. It should be noted that ZRT was checked before detailed modelling was conducted to ensure that the root models fitted within the centrifuge container or DSA used for the testing described later and shown in Fig. 2 (DSA) and Fig. 3 (centrifuge) and also within the build volume of the 3D printer. The diameter of structural roots has tapered laterally to less than 5 mm beyond the range of 0.5 m for most species, and these fine components cannot in any case be resolved within the 3D printer (see below) and were therefore not modelled. XYZ coordinates were input manually into the Solidworks 2012 3D computer-aided design (CAD) software, and spline functions were employed to link these coordinates and generate each root segment. The output file (Fig. 4) was directly input to a 3D printer (discussed in more detail below).

Fig. 1

Design process for 3D root model

Fig. 2

Schematic of large direct shear apparatus (DSA)

Fig. 3

Schematic of centrifuge model geometry, showing instrumentation and position of root analogues (dimensions at prototype scale in meters)

Fig. 4

Model tree root cluster. a 3D CAD file with different root diameter classes represented by different colours. b Printed model with temporary scaffold. c Model after dissolution of scaffold. d Straight root analogues

In order to benefit most from using the centrifuge modelling technique to model a prototype slope of larger height at small scale (1:N), without generating unwanted boundary effects from the model container, it is desirable to select a high value of N. However, as suggested by Kutter (1995), for modelling of reinforced slopes, if a structure (e.g. a fine root in this study) is relatively small compared to the median size of the grains (D ₅₀), then the soil may no longer behave as a continuum but more as a set of discrete particles. Therefore, to minimise potential grain size effects on the root–soil interaction, the lowest possible N is preferred. Furthermore, there are limits on the controllable frequency range of all mechanical actuators used to simulate earthquake shaking (this range is 40–400 Hz for the machine used in the later centrifuge testing); hence, selection of a high value of N is again preferred, to avoid losing the lower frequency components of the ground motion which result in most of the slip. A scale of 1:10 (N = 10) was ultimately decided upon as a suitable compromise for model testing given these competing effects. According to typical scaling laws (e.g. Kutter 1995), a model was scaled purely geometrically, with model material properties (e.g. T _r, E) scaled 1:1, i.e. that the analogue roots have similar strength and stiffness compared to real roots. As this was achieved (see later), because of the geometric scaling, the model will not be distorted in shape, and so, the root architecture was the same between model and prototype.

The reduced scale of the model combined with the threshold dimension that could be resolved within the 3D printer (0.75 mm) meant that roots <7.5 mm at prototype scale were not considered. However, roots with diameter in the range of 5–7.5 mm at prototype scale were catalogued into the 7.5-mm range and modelled as such. The fine roots below this size are highly flexible and always break when subject to soil slippage (Stokes et al. 2009) as their area reduces more rapidly than the tensile strength increases as diameter reduces, meaning they have a lower force contribution at break. As a result, the effect of fine roots on slope stability could be investigated separately if desired without considering them to be connected to the structural roots.

The stump

The stump, especially in older and mature trees, integrates various root segments to a coherent whole and plays a specific role in stability. Danjon and Reubens (2008) suggested defining the stump as the portion with a fixed depth of the first-order roots in trees. A value of 0.25 m was employed for real trees planted in slopes by Danjon et al. (2008). The same value was used in this study at prototype scale. The tree component above the ground surface was simplified to be a circular rod. Following this criterion, the diameter of the stump was equal to DBH which was 23 mm at model scale (full details are presented in the following section).

Root diameter and length

All of the roots except the tap root were simplified to circular rods of uniform diameter along their length (though not necessarily straight). The cross-sectional dimensions of these roots were grouped into four classes, as shown in Table 1 (After Watson et al. 1995), with each class being represented by a single diameter which was approximately the midpoint in the class range. These diameters were selected as representative of each class to be similar to those used in previous studies (Sonnenberg et al. 2010), and the use of only a few representative diameters also simplified the subsequent quantification of the mechanical properties. Considering that the shear zone thickness was approximately 20 mm in similar previous DSA tests (see Duckett 2013), if fine roots are held in tension principally (Stokes et al. 2009), roots longer than 20 mm will behave similarly and root length would be expected to have no significant effect on root contribution.

Table 1 Root diameter class for medium and structural roots (not including tap root)

Tap roots are very common across species, and they have been observed in up to 75% of tropical trees (Klinge 1973) and in 73% of Mediterranean woody species (Canadell et al. 1996). In Khuder et al. (2007), the largest vertical root was classified as the tap root. Typical tap root penetrations are within a range of 0.5 to 1.4 m based on limited literature (Crook and Ennos 1998; Danjon et al. 2008). A series of taproot morphologies were recorded by Danjon and Reubens (2008), based on 22 12-year-old Pinus pinaster. These patterns can be regarded as a taper varying with root inclination and branching rate. The pattern of taper will affect the bending moment capacity of the tap root (Goodman 2001). A tap root penetration depth of 1.0 m at prototype scale was selected here. A taper pattern, the diameter of which was 0.12 m (prototype) at the midpoint (reducing from DBH at the ground surface/top), was selected.

The ZRT was determined as a function of DBH, following the procedure introduced by Danjon and Reubens (2008). In Danjon et al. (2005), the radius of ZRT was set to be 2.2 times DBH. This ZRT was approximately achieved within the printer build volume lateral dimension threshold of 50 mm for the model scale DBH of 23 mm (ZRT = 2.17 × DBH).

Upslope/downslope root distribution

Trees on slopes tend to develop a specific asymmetrical architecture in the root system compared with trees growing in level ground. Various factors, that may be related to root growth, e.g. soil strength (Kirby and Bengough 2002), slope angle (Di Iorio et al. 2005; Pérez 2012), vegetation layout (Fan and Lai 2014) and static loading (Genet et al. 2008), have been reported. It is generally believed that roots growing upslope develop to a greater extent than roots growing downslope (see Nicoll et al. 2005; Danjon et al. 2013) despite some exceptions (e.g. Danjon et al. 2005; McIvor et al. 2008; Sonnenberg et al. 2011). This phenomenon appears related to the orientation of the aerial part of the plant which may be vertical or perpendicular to the slope surface (Mao et al. 2015). More roots develop and thicken in the up–down direction of the slope in response to the overturning moments induced by the combination of the inclination and the weight of the stem interacting with the root morphology of a particular species (Danjon et al. 2013). Apart from the inclination and weight of the stem, wind also contributes to the development of root morphology on the slope (Danjon and Reubens 2008). However, wind effects are highly site-dependent and will not be considered further in this paper. An asymmetrical architecture in the root system was here modelled by reducing the ZRT and number of sinkers on one side in the model (see Fig. 4). During the subsequent DSA and centrifuge testing, the model roots were tested with the most heavily rooted side towards the ‘downslope’ direction (e.g. as illustrated for papaya trees with the vertical trunk in Marler and Discekici (1997)).

Fabrication

3D printing techniques (also known as ‘rapid prototyping’ and ‘stereolithography’) were used to fabricate the model roots using a Stratesys Inc. uPrint SE acrylonitrile butadiene styrene (ABS) prototyper. ABS plastic is a low-cost material, ideal for structural applications when impact resistance, strength and stiffness are required. It is widely used for machining pre-production prototypes since it has excellent dimensional stability. The rapid prototyper takes a continuous spool of plastic filament which is melted and injected into place, a process known as fused deposition modelling. In this way, the machine may be thought of as a 3D inkjet printer that uses molten plastic rather than ink. The machine is computer numerically controlled (CNC) and once given an input file containing the 3D geometry (exported directly from the 3D CAD model) will operate entirely automatically. This is useful as the root models described herein took approximately 22 h to fabricate. Given the complex geometry of the model, the rapid prototyper automatically generated a scaffold of a contrast material to support the model roots as they hardened. The scaffold material was subsequently dissolved in a bath of caustic soda. It should be noted here that before the fabrication of root model, preliminary assessments were performed to identify the material properties of the printed ABS plastic via tensile tests of some straight rods (see Fig. 4d) with diameters matching those used in the model root system; the results of this testing will be presented in the following section and compared to actual root material.

Multiple 2D maps showing where roots intersect potential slip planes at different prototype scale depths (Fig. 5) provide an overview of the spatial structure of root reinforcement. Planes at shallow depth showed intersections with the main tap root; below the model scale depth of 100 mm (prototype = 1.00 m), despite the low number of roots (Fig. 6a), the root cross-sectional area (CSA) is high (Fig. 6b) for the following reason: Deeper planes demonstrated abundant medium and fine roots, but the corresponding CSA is relatively low. Root CSA is directly related to root area ratio (RAR) and root biomass, which indicates that the major part of the root model is still concentrated in the upper layers of soil.

Fig. 5

Distribution of roots intersecting four planes at different soil depth for the 3D root model at prototype scale (downslope positive)

Fig. 6

Vertical distribution of model roots. a Number of roots. b Root cross-sectional area (CSA)

Testing methodology

Material testing of individual straight analogue roots

Root analogue mechanical properties (stiffness and strength) were characterised by performing axial tensile tests using an Instron 4204 loading frame. Straight root samples were required for these tests, and these were fabricated using the same method described above. All of the test rods were 200 mm long but were of differing diameters (0.8–16 mm). The UTS of an analogue, T _r , was calculated as the force F in tensile testing required to induce breakage of the root divided by its cross-sectional area a _r . Young’s modulus was determined as the slope of stress–strain curve within the elastic deformation region.

Large direct shear tests at 1g with additional confining stress

Some preliminary direct shear tests on the buried root model using a modified conventional DSA with a free soil surface have previously been reported in Liang et al. (2014). This paper reports a separate and more extensive test programme using a much larger purpose-built DSA (Fig. 2), to simulate a shear plane passing at various depths through the model root system. Compared with the conventional DSA used previously, this apparatus has three main advantages:

• (a) A larger maximum shear displacement can be achieved (75 mm, compared to 20 mm), so that soil containing ductile root analogues could be sufficiently deformed to mobilise the full capacity of the rooted soil and break the roots, if required.

• (b) The top half of the shear box has a depth of 160 mm, which allows the shear plane to be located at different depths along the 150-mm-deep model root cluster by varying the amount of soil fill in the upper half and applying a suitable surcharge to simulate the confining stress at the simulated depth at prototype scale.

• (c) The cross section of the DSA is large compared with the ZRT of the model root clusters allowing low RAR conditions to be modelled and is similar in total area to the slope face in the subsequent centrifuge tests.

The model root cluster was suspended within the centre of the DSA from thin wires attached on the boundary of the DSA. Dry HST95 Congleton silica sand (Bent Farm, Congleton, Cheshire) was then pluviated in air around the model roots to a uniform relative density of I _D  = 55–60%. The basic properties of this sand are shown in Table 2 after Al-Defae and Knappett (2014), and the particle size distribution is shown in Fig. 7. Strength characteristics obtained through direct shear tests yield values of 40°and 32° for peak and residual (or critical state) angles of friction at the target density. Varied confining stress at the shear plane was obtained by changing vertical weights applied on the top of the sample to be able to represent the full range of stress levels expected at different potential shear plane depths in the centrifuge tests.

Table 2 State-independent physical properties of HST95 silica sand (Al-Defae and Knappett 2014)

Fig. 7

Particle size distribution for HST95 silica sand

A total of 14 rooted tests have been performed in this study, as summarised in Table 3. Each rooted test additionally had a reference fallow test with the same testing condition (depth of shear plane and confining stress) such that the root contribution to the shear strength (c _r ′) could be determined. Test series A was designed to investigate the effect of the changing RAR of the 3D root cluster for a constant confining stress at the shear plane of 8 kPa (representing the geostatic stress at a shear plane depth of 0.5 m in the centrifuge prototype slope, and therefore keeping the dilative component of the soil shear strength at the shear plane the same—Bolton 1986). However, the shear plane was made to pass through the root models at different model depths (50, 75, 100 and 125 mm; cross sections on these planes can be seen in Fig. 5) so that the RAR could be varied (3.99, 2.03, 2.86, 1.03%, respectively), as shown in Fig. 12. Here, RAR was defined using the following equation

$$ RAR={A}_r/ A $$

(4)

Table 3 Summary of large DSA testing plan

where A _r is the total root cross-sectional area and A is the area of CRZ (CRZ and ZRT are the same for the 3D root cluster) rather than the total cross-sectional area of the DSA. This was selected based on the previous observation of Liang et al. (2015) that root influence on shearing properties is confined to an area local to the roots.

Test series B utilised a simpler group of straight (vertical) root analogues at the same confining stress (8 kPa), such that the effect of root morphology could be investigated through a comparison with test series A. The straight analogues were 150 mm in length, printed in ABS using the same 3D printer and were installed to provide the same RAR and spatial distribution of roots of different diameter on the shear plane as the 3D root cluster at each of the different depths as shown in Fig. 5.

Test series C utilised the 3D root cluster as in series A; however, the confining stress was varied at the same time as the position of the shear plane mimics the confining stress levels in the subsequent centrifuge tests at the prototype slip plane depths of 0.5, 0.75, 1.0 and 1.25 m. These tests were also used to investigate the effect of confining stress at a given RAR through comparison with series A.

Finally, test series D utilised the straight analogues of series B, with the varied stress levels of series C, again to mimic the confining stress levels in the subsequent centrifuge tests (described in the following section).

Ten acceleration due to gravity centrifuge modelling

Dynamic centrifuge testing, using earthquake ground motion to induce slope failure, was conducted using the 3.5-m-diameter beam centrifuge and servo-hydraulic earthquake simulator at the University of Dundee. Three models, referred to herein as TL03, TL04 and TL07, were flown, representing identical soil slopes with a slope angle of 27° at 1:10 scale in dry HST95 sand at I _D  = 55–60%, and tested at 10g. The model slope TL07 included the 3D model root analogues (Fig. 3), while model TL03 was reinforced with straight ABS rods having the same RAR and spatial distribution at the level of the middle of the 3D root cluster to investigate root morphology, similar to the series B DSA tests. Further detail about the design of model TL03 can be found in Liang et al. (2015). Model TL04 was a fallow reference case of identical slope geometry. A dry soil model was selected to eliminate hydrological effects (specifically any liquefaction) to focus on the pure mechanical reinforcement effect. During pluviation, the soil was instrumented with 13 type ADXL-78 MEMS accelerometers (ACC) (Analog Devices) measuring the horizontal accelerations within the soil specimen and three external linear variable differential transformers (LVDTs) (RDP group) measuring settlement, with one installed at the centre of the crest and the other two placed adjacent to the side walls to determine whether there were any boundary effects (not used herein). Eight successive earthquake motions were applied to each model at 10g, comprising three different records with distinct peak ground acceleration (PGA), duration and frequency content, as shown in Table 4. Full details of the model preparation can be found in Liang (2015).

Table 4 Sequence of input motions in centrifuge tests

Results and discussion

Root mechanical properties

The tensile strength of cylindrical test specimens demonstrated a negative power law relationship when plotted against root diameter, which conforms with typical relationships for plant roots. The tensile strength, T _r (MPa), and Young’s modulus, E (GPa), as a function of the root diameter were found to be

$$ {T}_r=57.886{D}^{-0.523} $$

(6)

$$ E=3.24{D}^{-0.55} $$

(7)

Plotting these values for various root diameters on a material selection chart (Fig. 8) shows that the printed analogues fall within the range of real roots collected from the literature (e.g. Operstein and Frydman 2000; Mickovski et al. 2009; Fan and Su 2008; Mao et al. 2012) and that they are superior to wood and rubber root analogues used previously (results from Duckett (2013) are also shown in Fig. 8).

Fig. 8

Comparison of material properties between tree roots and root anlogues based on material selection chart

Direct shear tests: RAR effects within tap root system

The resulting shear resistance of rooted and fallow soil is plotted against shear displacement in Fig. 9, for the tests from series A. The fallow soil exhibited a strain softening response with peak strength higher than that at larger strain. The reinforced cases had all reached peak strength by this point and remained constant, so the additional shear resistance from the analogues was therefore determined from the differences in peak shear strength. As demonstrated, with different RAR, the additional shearing force was very similar when the shear plane was at the depth of 50, 75 or 100 mm. Given that RAR varied significantly for these three planes (3.99, 2.03 and 2.86%, respectively), the additional shear force should be quite different according to current understanding. The main difference of these three cases to the 125-mm case (RAR 1.03%) was that the tap root crossed the shear plane in the former, but not in the latter. This indicated that the presence of the tap root was very important in mobilising the whole 3D root system to resist shear loading. This would be of great interest for engineering applications, especially as a guide to the planting approach for improving slope stability in terms of the importance of selecting a type of vegetation that has a deep taproot. This has not been reported previously, the reason for which may lie with the difficulty in generating an identical tap root system architecture for comparative testing and the general non-repeatability of shear testing on rooted samples in situ. The ability to produce repeatable model root systems (yet of realistically complex 3D geometry) is therefore a significant benefit of the 3D printing approach described in this paper for future studies on root–soil interaction.

Fig. 9

Results of large DSA tests on 3D root cluster reinforced soil with constant confining stress along different potential slip planes

In general, roots below the tap root broke in quantity (Fig. 10) when the root cluster was subject to moving soil; however, whether the tap root broke or not varied with the test conditions and could be classified into two types of mechanism. The first involved breakage of the root cluster broke at the fork of the tap root (e.g. the 50-mm case). For this test, the lower part disconnected from the upper component and ceased to contribute to the additional shear resistance. Following the disconnection (at ∼17-mm displacement, see Fig. 9, ‘Evidence of root breakage’), the shear resistance decreased rapidly, but there was still a 4.8% increase compared with the fallow case. It is interesting to note that this lower resistance after tap root breakage is approximately the same as in the 125-mm case where the tap root does not cross the shear plane. This indicates that the drop in resistance was due to tap root breakage. The second mechanism involved progressive breakage, in which the increased shear resistance was maintained despite some root breakage of individual deep roots. This is also evident in the resistance–displacement curve (compare the 75 and 50-mm cases in Fig. 9—both initially have the same resistance, but the 75-mm case retains this resistance throughout the test, even with the substantial root breakage shown in Fig. 10). The load displacement data also suggested that deep roots were mobilised and broken in sequence rather than simultaneously (see the numerous sudden small drops on the 75-mm curve, Fig. 9)—this has also been observed previously for straight (vertical) root analogues in large DSA tests (Liang et al. 2015).

Fig. 10

Breakage of root clusters under constant confining stress following careful post-test exhumation

The distinct difference between these two mechanisms can be related to the position of shear plane, as shown in Fig. 11. As the upper part of the root cluster rotated against the moving soil, the tap root acted mechanically as a lever arm. Commonly, ‘lever arm’ structural roots failed at a point where they branch (Tobin et al. 2007), as shown in the tests with tap root breakage. The fulcrum of the lever depended on the depth of the shear plane. As the shear displacement of the top soil was simulated identically between different tests, for cases with a shallow shear plane depth, higher bending moment was transferred to the point of branch.

Fig. 11

Inferred root behaviour as a function of slip plane depth

Direct shear tests: root morphology effect

In the preceding section, it was indicated that the additional root contribution to shear resistance was principally associated with the presence (or not) of the tap root. To better understand the effect of the 3D geometry compared to straight root assumptions typically made in prediction models, the resulting shear resistance of rooted and fallow soil for test series B is plotted against shear displacement in Fig. 12. Following the test, straight analogues were excavated carefully and none were observed to have broken, as shown in Fig. 13. It can also be seen from Fig. 13 that the flexible fine analogues (e.g. of 3, 1.6 and 0.8 mm diameter) showed some irrecoverable bending deformation, while the rigid root behaved within the elastic range (see Liang et al. 2015). Compared with the behaviour observed in Fig. 10, it clearly indicates that the interconnection of roots in the 3D cluster can result in substantial concentration of stress at key points within the root architecture.

Fig. 12

Results of large DSA tests on straight root analogue reinforced soil with constant confining stress along different potential slip planes

Fig. 13

Post-test observation of large DSA tests on straight root analogue reinforced soil

A comparison of the inferred root cohesion between the 3D root cluster (series A) and the straight root tests (series B) is shown in Fig. 14. Root cohesion c _r ′ was determined by dividing the additional shear resistance force above fallow from Figs. 10 and 12 by the area of the CRZ of the 3D root cluster such that the values can be fairly compared. It should be noted here that the additional shear resistance was taken as additional ultimate resistance at the large shear strain (critical state). The variation of c _r ′ with depth for the straight analogues is, unsurprisingly, highly similar in shape to the root CSA distribution (Fig. 6b). This suggests that root cohesion for a bundle of individual root analogues is highly dependent on root density, which has been confirmed by many previous researchers (e.g. Pollen and Simon 2005; Loades et al. 2010). The effect of the interconnection of roots in the 3D cluster, however, can be seen to ‘smooth-out’ the zone of comparatively lower CSA around the 75-mm depth shear plane, due to the ability to structurally transfer load between roots in the 3D case. Figure 14 also shows that c _r ′ of the 3D root cluster is consistently lower than that of straight analogues. This is not unexpected as the individual roots (each root was of identical 150-mm length) of the straight analogue group had higher anchorage length in the lower layer of soil (approximately 445, 800, 128 and 48% higher for 50, 75, 100 and 125-mm cases, respectively) and therefore also a higher overall volume of root material (approximately 60, 124, 190 and 98% higher for 50, 75, 100 and 125-mm cases, respectively).

Fig. 14

Comparison of 3D root cluster and equivalent straight root analogue shear tests, with constant confining stress at slip plane (8 kPa)

It is also consistent with the observation that the straight analogues did not break. The observations made in this section suggest that it may not be suitable to apply theories based on straight vertical rods directly to the tap root systems of trees (as also suggested by Pollen and Simon 2005).

Direct shear tests: confining stress effect

In this section, DSA tests with varied confining stress at the shear plane (test series C) are presented. The results are compared to those of test series A with the same RAR but constant confining stress at all depths (8 kPa) in Fig. 15. There was an initial positive correlation between c _r ′ and the slip plane normal effective confining stress within the soil from 50 to 75 mm. This was consistent with the in situ tests on four common Australian riparian trees (Docker and Hubble 2008) and laboratory tests on a model straight root group (Liang et al. 2015).

Fig. 15

Effect of confining normal stress at slip plane on measured root cohesion, 3D root models

However, as the shear plane moves closer to the bottom of the root cluster and the anchorage length reduces, c _r ′ reduces. This is not surprising and can be interpreted with lateral loaded pile theory. When the anchorage depth is higher than a critical value, a maximum amount of shear resistance will be mobilised and further increase of length beyond this will have no impact on the magnitude of the shear resistance. However, once the root anchorage depth is lower than the critical anchorage depth, the shear resistance of the root will be governed by the anchorage depth. Hence, compared with 75-mm case, 100-mm case had zero main structural tap root anchorage depth and lower values for other analogues, which is why c _r ′dropped so dramatically at this shear plane.

These observations would again suggest that current c _r ′ prediction models (e.g. Wu 1976; Pollen and Simon 2005) which assume that the root cohesion is independent of soil type and stress level may not provide suitable predictions of root contribution to soil strength. Figure 16 shows the breakage of roots during the varied confining stress test series following careful post-test exhumation. Compared with Fig. 10, more roots broke at the higher confining stresses. This was not surprising as the frictional restraint between root and soil will be higher, and hence, it is much more likely that the roots will break, rather than pulling out.

Fig. 16

Breakage of root clusters under varied confining stress following careful post-test exhumation

Validation and implication of scaled root model

In engineering practice, when calculating the mechanical contribution of vegetation on slope stability, anchorage benefit provided by large roots (>20 mm) are always overlooked. The main reason for this is due to the difficulties in measuring or quantifying such benefits considering the large size of whole tree root systems and the limited size of available shear apparatus (Sonnenberg et al. 2011). Using scaled root models in the centrifuge can fill this gap and quantify the potential benefit of the larger structural roots on preventing slope failure. Considering this distinction, for the sizes of analogues considered in this study, 3 mm, 5 mm and tap roots represent architecturally the roots larger than 20 mm, while 1.6 and 0.8-mm-diameter analogues represent the small/fine roots.

Following a detailed data gathering exercise, a series of in situ shear tests on five tree species (linden hibiscus, Japanese mallotus, Chinese tallow tree, ironwood, white popinac) reported by Fan and Chen (2010) is employed here for comparison. Those tests were performed in sandy soil with a shear plane depth of 100 mm. The root area ratio measured was between 0.1 and 5%, which is in the same magnitude range as this study. Considering the similar test conditions (soil, root density and root morphology), the root cohesions should be comparable with the DSA data reported herein. The measured root cohesion c _r ′ values for those five tree species are 11.32 ± 4.42, 20.21 ± 6.51, 31.96 ± 13.08, 46.33 ± 17.04 and 47.87 ± 19.02 kPa for linden hibiscus, Japanese mallotus, Chinese tallow tree, ironwood and white popinac, respectively. For the model root cluster used in this study, the measured root cohesion c _r ′ for the confining stresses of 8 and 16 kPa at the shear plane of 100 mm are 18.95 and 20.52 kPa, respectively (see Table 3), falling towards the lower end of the range for the field tests. This clearly indicates that the printed 3D root cluster does not just provide a representative root mechanical properties for individual roots (see Fig. 8) but also a reliable root–soil interaction response compared with field measurements.

Global performance of rooted slope in 10g centrifuge modelling

In the centrifuge tests, the effect of the 3D root cluster on the shear strength of the rooted soil within the CRZ is most closely represented by the DSA tests from series C, which mimic the higher confining stresses at deeper potential shear planes that may form when the slope becomes unstable. The effect of the straight root analogues under a similar stress regime is represented by series D. A comparison of c _r ′ values between these two cases is shown in Fig. 17.

Fig. 17

Comparison of 3D root cluster and equivalent straight root analogue shear tests, with varied confining stress at slip plane

The key indicator of slope performance considered herein is the crest settlement (proportional to permanent slip) and the acceleration response spectra at the crest (Roy et al. 2016). A comparison of the crest settlement between the three slope models is shown in Fig. 18. The interesting aspects of behaviour were as follows:

• (a) For both root-reinforced slopes and the unreinforced slope, a decreasing trend of settlement was observed when the slope was subjected to successive motions (e.g. aftershocks). This can be associated with slope geometry change (re-grading), as previously reported for fallow and pile-reinforced slopes by Al-Defae et al. (2013) and Al-Defae and Knappett (2014), respectively.

• (b) The presence of root analogues results in a significant reduction (by 61 and 76%, for straight roots and 3D root cluster, respectively) in permanent slope movement compared with the fallow case, especially in the first two motions (EQ1 and EQ2; Fig. 18). This can be interpreted as the rapid mobilisation of root–soil interaction due to the initial soil slip under dynamic loading (see Liang et al. 2015).

• (c) After the first two motions, relatively smaller reductions (in total 14 and 27%, for straight roots and 3D root cluster, respectively) were observed, which indicates that the additional resistance from the roots is largely constant after the initial mobilisation.

• (d) In contrast to the DSA tests, no roots were observed to have broken following careful post-test exhumation, which would infer that the maximum root–soil resistance was mobilised after EQ2 (see previous point) and that this was either associated with (i) yielding of the soil around the root analogues or (ii) the strengthening effect of the roots forcing the slip plane deeper within the slope.

Fig. 18

Comparison of permanent crest settlement of fallow and root-reinforced slopes from centrifuge testing

Figure 19a shows a comparison of acceleration response spectra (ARS) using the motion measured at the crest of the slope (instrument 12) between the slope containing the 3D root clusters and the fallow slope. ARS was normalised by the peak acceleration of the input motion and determined for the case of a structural system having nominal 5% damping. The ARS for the 3D root-reinforced slope at the crest was similar to that of the fallow slope, and less than 3% reduction on average over the entire frequency range due to the root presence was observed in any of the three distinct motions. The ARS of EQ1 for the rooted slope was approximately 13% higher than that of the fallow slope over a limited range of periods. This indicates that the presence of the roots had a limited influence on the general propagation and amplification of earthquake motion from the toe of the slope to the crest of the slope, in contrast to the 76% reduction they had on slip displacements (Fig. 18). Approximately 18 and 10% reductions in ARS magnitude were observed in the near field of the root analogues for EQ2 and EQ5, respectively, as shown in Fig. 19b, which may have reduced settlement at the slope crest. The difference of ARS measured at instrument 11 and instrument 12 suggests that roots can only influence the slope’s response within or close to the CRZ.

Fig. 19

Normalised acceleration response spectra (ARS) of three distinct motions between rooted slope (TL07–3D clusters) and fallow slope. a At the crest of the slope. b At the location of a model root cluster

Root morphology effect observed in centrifuge modelling

Figure 18 shows that the slope reinforced by 3D root clusters provides a similar reduction in slope movement compared to the straight root case. However, the variation of root cohesion c _r ′ with depth of the straight root case was generally much higher, as indicated in Fig. 17 (and as found in the DSA tests for constant confining stress in Fig. 14). It may therefore not be suitable to relate the reduction of permanent settlement (i.e. the reinforcing effect) solely to a change in c _r ′. Given also the observation of no root breakage in either rooted test in the centrifuge, it may therefore be true that in both centrifuge cases, where the position of the shear plane is not constrained as in the DSA, but can form wherever it is energetically most efficient, the root reinforcement (c _r ′) was large enough to force the slip plane into a new deeper position within the slope and that this is both the reason for the similarity in settlements between TL07 and TL03 and the mechanism by which the roots reinforce the slope (i.e. once the root contribution is large enough to force the change in mechanism, there is no additional benefit from further increases in root strength contribution). The ARS at instrument 11 was also checked and showed no considerable difference between 3D and straight root cases in the near field of the root analogues, with the exception of EQ5.

Conclusions

A tree root cluster model that can be used in centrifuge model tests to more realistically model the geometry, spatial distribution and mechanical properties of tree roots was designed and fabricated from ABS plastic using the 3D printing technique and subsequently used to investigate the seismic behaviour of rooted slopes. A series of element tests in a large direct shear apparatus were also conducted to investigate some of the key parameters that may affect the root and soil interaction and support the centrifuge tests. The following conclusions can be drawn from this study:

• 1. When subject to shear loading in soil, the tap root appeared to structurally transfer load within the root system, including to smaller roots and deeper roots which subsequently broke or were pulled out. This suggests that future models must be able to account for the detailed interlinked root morphology of trees. This will make the 3D printing technique developed in this paper increasingly valuable for the validation of such analytical models through further testing of scaled physical models.

• 2. Root reinforcement is not only a function of root mechanical properties but also depends on factors including surrounding confining stress (resulting in depth dependency even for the same RAR), depth of the slip plane and root morphology. The repeatability of the printed analogues was particularly useful to identify these effects through like-for-like comparison, which is impossible in the field due to the natural variability of tree root systems and which will be of significant benefit for future laboratory studies of root–soil interaction.

• 3. The presence of root analogues resulted in a significant reduction (by 61 and 76%, for straight roots and 3D root clusters, respectively) in seismically induced permanent slope movement (crest settlement) compared with the fallow case in the centrifuge tests. The root analogues had a limited influence (less than 3% reduction of spectral acceleration) on the general propagation and amplification of earthquake motion from the toe of the slope to the crest of the slope. However, reductions of spectral acceleration of between 10 and 20% were observed in the near field (CRZ) of the root analogues.

• 4. Comparisons of perpendicular straight root element tests with those for 3D clusters showed significantly higher measured root contribution in the former; however, centrifuge tests of the complete slope stability problem showed similar improvements to performance in both cases. This suggests that the mechanism by which roots reinforce slopes in a first-time slip is not through adding additional root cohesion to the fallow slip plane but involves the stronger zones around the roots forcing the slip plane to a deeper position within the slope (and to a similar position in both straight and 3D cases). Therefore, obtaining the maximal reinforcing effect may require increasing the depth and lateral extent of the rooted zone which has enhanced shear strength, rather than just obtaining biomechanically stronger roots (as conventional models may suggest).

A, area of plane; A _r , total cross-sectional area of roots; a _r , individual cross-sectional area of root; c _r ′, root cohesion; CRZ, critical root zone; CSA, cross-sectional area; d, depth of soil; D, diameter of root; DBH, diameter at breast height; D ₅₀, particle diameter at which 50% is smaller; E, Young’s modulus; e _max , maximum void ratio; e _min , minimum void ratio; F, force in tensile testing; G _s , specific gravity; g, acceleration due to gravity (=9.81 m/s²); I, second moment of area; I _D , relative density; j, weakest root removed at each simulation; k, liner coefficient for ultimate tensile strength; m, exponent coefficient for ultimate tensile strength; n, root number ordered from strongest to weakest; N, total root number; RAR, root area ratio; RAR _n , root area ratio of root size n; R _f , root orientation factor; T _r , ultimate tensile strength; T _rj , tensile strength of weakest remaining root; T _rn , tensile strength of a root of size n; x, shear displacement; Y, cumulative root fraction; ZRT, zone of rapid taper; z, thickness of shear zone; α, exponent coefficient for Young’s modulus; β, depth coefficient; θ, angle of shear distortion of root; ϕ′, effective friction angle of soil