1 Introduction

Plant roots are natural soil anchors that resist external stress from wind loading or soil movement (e.g. landsliding), and also act as reinforcement of soil. They could play a prominent role in risk mitigation strategies (e.g. as a nature-based solution, or NBS) through modifying the groundwater regime [42] and/or through direct mechanical soil strength improvement [19, 36, 43]. With the increased frequency of windstorms [35] responsible for recent damage in European forests [14, 26], root–soil interaction is of interest in Forestry, where it is important to understand how forest management plans can be adapted to increased natural hazards [12, 30]. It is also of interest in Civil Engineering, where trees could provide protection against earthquake-induced landslides [22] or be employed as cost-effective protective barriers [36] to protect sensitive structures from wind loading or debris flow [6].

As it is time-consuming and expensive to conduct large field trials and impossible to find identical root systems for systematic studies, an increasing body of research has used scaled (1:N) root models, which can be tested in the laboratory, either under 1g [16, 28, 46] or centrifuge conditions [21, 47], to investigate root–soil interaction. To avoid unwanted boundary effects from the model container, it is desirable to select a high value of N. However, this in turn may cause scale effects related to particle size [20] where the surrounding soil could no longer be considered as a continuum in coarse granular soils. Even where root–soil interaction tests can be conducted at 1:1 (full) scale, small diameter roots may similarly be subject to particle size effects in coarse-grained soils.

Previous experimental studies have led to a good understanding of particle size effects on foundations, from centrifuge testing. The threshold where particle size effects can be neglected depends on the stress level [4], shear band formation [3, 17, 18] and structure-soil interface roughness [41]. Garnier et al. [15] proposed that there was some deviation from continuum behaviour when the ratio of shallow footing diameter to grain size was less than approximately 35. However, shaft resistance of a pile subjected to loading or stress relief in dilatant soils was affected by its diameter even at large scale according to cavity expansion theory [32, 48]. Compared with foundations, individual roots could result in much lower structure-grain size ratios, which has been rarely investigated [28].

Considering the tree-uplift problem shown schematically in Fig. 1, under external lateral loading (e.g. from a windstorm), individual roots (excluding the taproot) are idealised to be either pulled-out or pushed into the soil. The aim of this study is to numerically model the vertical sinker roots of such a root system, which are loaded axially, validating performance against laboratory testing, and then using the model to investigate the influence, if any, of potential particle size effects that may arise in 1:N scale physical model tests of sinker roots of a particular prototype. For a higher N, chosen to avoid boundary effects in designing such tests, the model roots become smaller while the soil particle size remains the same; this changes the ratio of root diameter (\({d}_{\rm r}\)) to mean particle size (\({D}_{50}\)) under the same effective stress conditions. To this end, 3D DEM simulations were employed to replicate the behaviour of sinker roots of different diameters interacting with a granular bed. The results were then compared with equivalent FEM simulations, where the soil was considered as a continuum and hence did not incorporate particle size effects, to identify any deviations from continuum behaviour. The numerical simulations were then used to inform a simplified practical approach for estimating axial root resistance, incorporating particle size effects.

Fig. 1
figure 1

Simplified tree root system under lateral loading

Full size image

2 Summary of physical model tests of sinker roots for validation

Laboratory testing of root analogue pull-out reported by Liang et al. [23] provided useful data of sinker root–soil interaction and was used to validate the numerical model in this study. A brief description of the experimental set-up and properties are given below, to support the numerical modelling decisions that were made.

2.1 Soil

A well-graded medium to fine sand (HST95 Congleton silica sand) was used, which was air-pluviated around vertically orientated root analogues (mechanical simulants of living roots), within cylindrical containers (test chambers). HST95 sand is a specific fraction of the sand extracted at Bent farm, Congleton, Cheshire, and has been widely used at the University of Dundee in physical testing e.g. [1, 20]. The particle size distribution (PSD) of the HST95 sand is shown in Fig. 2, with some other index properties summarised in Table 1. The coefficients of uniformity and curvature of the soil were 1.5 and 1, respectively, and its maximum and minimum density were 1.8 and 1.5 g/cm3, respectively. The critical state friction angle of the sand was 32° based on direct shear tests across a range of relative densities (9–93%) and effective confining stresses (5–200 kPa), as reported by Al-Defae et al. [1].

Fig. 2
figure 2

Particle size distribution of the soil

Full size image
Table 1 HST95 sand properties and key numerical parameters for DEM (after Sharif et al. [39])
Full size table

2.2 Physical testing of axial sinker root–soil interaction

3D printed straight analogues fabricated from Acrylonitrile Butadiene Styrene (ABS) plastic with 150 mm anchorage length were used to simulate individual sinker roots under vertical pull-out conditions. The diameters of these plastic rods were 12 mm, 3 mm and 1.6 mm. A steel hook was attached to the top of each root analogue so that an Instron 5985L7706 loading frame could apply a vertical pull-out force to the model root. The test chamber was 150 mm in diameter and 500 mm in depth, which could minimise unwanted boundary effects during uprooting for even the largest rods considered [34]. HST95 sand was air-pluviated at a relative density of 55–60% (corresponding target dry density: 1.65 g/cm3) to a depth of 140 mm beneath a suspended single root analogue. The analogue was then pushed vertically into the soil by 20 mm. The process of pluviation was then continued until the analogue was completely buried. Three types of slotted circular surcharge weights (Fig. 3) were placed at the soil surface, in order to simulate four different levels of confining effective stress at the mid-depth of the root analogue, namely 1.2 kPa (no surcharge weight), 4 kPa, 8 kPa, and 12 kPa. The root analogue in each case was pulled out at a velocity of 10 mm/min using the aforementioned load frame (Instron 5985L7706, Instron Inc., USA). The capacity of the load cell used was 30 kN with an accuracy of 1 mN. Further details relating to the physical modelling can be found in Liang et al. [23].

Fig. 3
figure 3

Diagram of physical test set-up for sinker root uplift test (all dimensions in mm)

Full size image

3 Numerical modelling and validation

3.1 Overview of numerical modelling

DEM, which uses discrete particles to simulate large deformation problems in granular soils, has proven to be an adequate tool to simulate centrifuge conditions and properly capture soil interaction behaviour of varied complexity [8]. The great advantages of DEM are the high level of control of the initial conditions which allow rigorous comparative parametric studies and the automatic capturing of soil behaviour through simple particle contact laws [9]. This helps interpret the results of physical model tests and can also be used to conduct virtual laboratory tests which require specialist facilities (e.g. centrifuge) that are not always readily accessible [2]. In addition to DEM simulations, comparative FEM simulations were also performed, in which case the soil was modelled as a continuum that was free from particle size effects.

3.2 Numerical models and parameters

To model the root–soil interaction in DEM, Particle Flow Code (PFC) 3D 5.0.35 (Itasca Consulting Group, 2016) was employed. The mechanical interface interactions among soil particles and between the soil particles and root were modelled by a simplified Hertz-Mindlin contact model [29]. The contact model was calibrated to replicate Congleton HST95 sand behaviour. The particle size distribution (PSD) adopted was representative of the HST95 sand tested in the laboratory (see Fig. 2). The adopted contact parameters, summarised in Table 1, for particle–particle interaction were calibrated against existing drained triaxial element test data under an effective confining pressure of 60 kPa (see Fig. 4; after Sharif et al. [37]). Following Arroyo et al. [2], particle rotation was inhibited to simplify model calibration whilst obtaining a realistic macroscopic behaviour, as has been previously demonstrated both at element scale [7] and for boundary value problem simulations such as cone penetration tests [10, 40] and centrifuge modelling of screw piles [39].

Fig. 4
figure 4

Validation of drained triaxial tests of HST95 sand under a confining stress of 60 kPa: a Deviatoric stress—axial strain; b volumetric strain—axial strain (DEM and laboratory data from Sharif et al. [37]; FEM from this study)

Full size image

To calibrate the root–soil contact parameters, a series of interface constant normal load direct shear tests were conducted in the present study to determine the interface shear strength properties between HST95 sand at a pluviated relative density of 55–60% and printed ABS plastic specimens with surface texture and print layering matching the model root analogues along the shaft. The interface friction angle at critical state under effective normal stress levels from 8 to 200 kPa was 25°, suggesting a coefficient of particle-structure interface friction of 0.466 in DEM (Table 1).

To understand the soil behaviour when treated as a continuum, which could be compared with the DEM results to quantify any particle size effects, root–soil interaction was also simulated by two-dimensional (2D) axisymmetric FEM simulations using [33], using a nonlinear elasto-plastic constitutive model, ‘hardening soil with small strain stiffness’ (‘HS small’), where the mobilised stress conditions associated with typical limiting displacements, along with the influence of soil compressibility on these distributions, can be investigated. The key soil parameters summarised in Table 2 (parent soil) were previously calibrated for the HST95 sand by Al-Defae et al. [1] based on direct shear and oedometer testing. Independent triaxial tests (mentioned previously) were simulated numerically as shown in Fig. 4. When axial strain was lower than 5%, the shear behaviour in the laboratory test was generally well captured by both DEM and FEM, while as expected, beyond this point soil softening could not be captured by the HS small model (though this was captured in DEM).

Table 2 HST95 sand key numerical parameters for FEM (corresponding to 60% relative density; parent soil after Al-Defae et al. [1] and reference parameters defined at p′ = 100 kPa)
Full size table

3.3 Modelling procedure

The root analogues from the physical tests were first simulated in DEM, where the virtual soil bed, with density consistent with that of the laboratory tests, was created using the periodic cell replication method (PCRM, [11]). Following previous studies [38, 39], the particle size was scaled up in the radial direction with increasing radial distance from the root to minimise the computational time (Fig. 5a). The scale-up factors for particle size within different zones are listed in Table 3. Due to the significant increase in number of particles required for 12-mm-diameter root cases, with the computational time taken into consideration, only cases of 3 mm and 1.6 mm root under 4 different confining stresses (8 simulations in total) were utilised in the DEM validation. The periodic cylindrical cell (5.7 mm-high) used for vertical replication was designed to have a central non-scaled zone suitable for both of these two sizes of roots. The test chamber was designed to be wider than \(18{d}_{\rm r}\) (i.e. chamber radius \({B}_{\rm s}>{9d}_{\rm r}\)), with soil beneath the bottom (\({D}_{\rm s}\)) of the root deeper than \(2{d}_{\rm r}\), to minimise boundary effects. A rigid cylinder with interface friction coefficient 0.466 was used to model the sinker root. To model quasi-static pull-out in DEM, it was displaced upward vertically, at a constant speed of 4, 7, 10 and 12 mm/s in the 1.2, 4, 8, 12 kPa confining stress cases, respectively, to set inertial numbers in each case lower than 10–3 [8]. To reduce runtime, the length of the root analogues in DEM was scaled down by a factor of ten (i.e. scaled down from 150 to 15 mm) to limit the number of particles required along the length of the root, and g-level scaled up by a factor of ten (i.e. scaled up from 1 to 10g) to guarantee the same stress level (without contribution from the surcharge) at homologous points along the root length with the physical tests. It was systematically checked that the reduction of root length (i.e. \({l}_{\rm r}/{D}_{50}\) reduced from 1064 to 106.4) did not cause any particle size effects related to length. However, root diameters were kept identical to the physical tests (i.e. they were not scaled by 1:10) to avoid further grain size effects by the reduction of root diameter. An extra periodic cylindrical cell with zero friction and the required specific gravity (30.7, 19.3, 8 and 0.01 for 12, 8, 4 and 1.2 kPa, respectively) was placed on the top surface to simulate the surcharge, which aimed to achieve the same stress distribution (with contribution from the surcharge) along the root length with the physical testing. To avoid the interaction between the root and the surcharge a rigid frictionless cylinder was assembled on the top of the sinker root throughout the surcharge zone (Fig. 5a) and displaced at the same velocity as the root. The velocity and displacement of particles at boundaries (side and bottom boundaries) were fixed to be zero, except for the top boundary which was treated as a free surface.

Fig. 5
figure 5

Diagram of simulation models for validation of pull-out tests in: a DEM (shading indicates the different particle size scale-up factors applied); b FEM

Full size image
Table 3 Particle scale-up factor in different zones within the DEM test chamber of sinker roots
Full size table

In FEM, as the soil was considered as a continuum, only the 12 mm root cases, where \({d}_{\rm r}/{D}_{50}\) (= 85) was the largest amongst the testing conducted and where grain size effects could be assumed to be negligible, were simulated in FEM, where the domain was discretised with 8715 15-node elements. Following the geometry of the experiment, an axisymmetric model was developed (Fig. 5b). The radius of the test chamber (\({B}_{\rm s}\)) was 75 mm, while \({D}_{\rm s}\) was greater than \(2{d}_{\rm r}\). Consistent with the method adopted in the DEM, the length of the root was scaled down at 1:10 with a corresponding increase of g-level. A surcharge was applied on the soil surface as a uniform distributed pressure, having the same values as those applied in the experiments. The root was simulated by three weightless rigidly connected plate elements, to which the upward displacement was applied. These elements enclosed a rigid material having the same unit weight as the ABS plastic used in the simulations.

As a common, simple interface element applying a strength reduction factor (\({R}_{\rm inter}\)) to the adjacent soil properties cannot consider soil dilatancy, a zone of soil adjacent to the root was used to model the (non-negligible) thickness of the shear band, using the ‘HS small’ model with adjusted parameters (described below). According to Mascarucci et al. [25], the thickness of the shear band was modelled to be 3 mm (approximately \({20D}_{50}\)). The input friction angle for this interface zone (\({\phi }_{i}\)) was 25° according to the laboratory interface shear tests, while the dilation angle (\({\psi }_{i}\)) was determined to be 6° based on the assumption below:

$$\frac{\tan \psi_i }{{\tan \psi }} = \frac{\tan \phi_i }{{\tan \phi^{\prime} }}$$
(1)

where \(\phi^{\prime}\) and \(\psi\) are the soil peak friction and dilation angles, respectively. In the FEM modelling, the compression modulus of the zero-thickness interface element is proportional to the interface shear modulus \({G}_{i}\), which is scaled according to:

$$G_i = \left( {\frac{\tan \phi_i }{{\tan \phi^{\prime} }}} \right)^2 G$$
(2)

where \(G\) is the operative soil shear modulus [33]. Based on these definitions, for the case of a 3 mm thick shear band, the input interface reference oedometer modulus \({E}_{{\rm oed},i}^{\rm ref}\) of the soil was therefore scaled in the same manner as:

$$E_{{\rm oed},i}^{\rm ref} = \left( {\frac{\tan \phi_i }{{\tan \phi^{\prime} }}} \right)^2 E_{\rm oed}^{\rm ref}$$
(3)

To reduce the number of independent parameters, \({E}_{50,i}^{\rm ref}=1.25{E}_{{\rm oed},i}^{\rm ref}\) and \({E}_{ur,i}^{\rm ref}=3{E}_{{\rm oed},i}^{\rm ref}\) were assumed, as in the parent soil [1]. All other parameters were unchanged from the parent soil. According to a preliminary parametric study, \({\psi }_{i}\), \({E}_{{\rm oed},i}^{\rm ref}\), \({E}_{50,i}^{\rm ref}\) and \({E}_{ur,i}^{\rm ref}\), determined based on the discussion above, were strongly correlated with the stiffness of root pull-out behaviour, while the pull-out capacity was dominated by the interface friction angle, which was physically measured in the present study. These key parameters are summarised in Table 2. Simulations were also conducted with the reference small strain shear modulus and the reference strain scaled based on the method above, yet the results showed no significant difference.

3.4 Validation of sinker root–soil interaction

In order to compare with results from the laboratory tests, the uplift force and displacement post-peak obtained from the DEM and FEM numerical simulations were both scaled up by a factor of 10 due to the reduction of root length from 150 to 15 mm. However, the displacement pre-peak was scaled 1:1 as the uplift displacement where the peak force appeared does not change even with the reduction of the root length. This was because the shear strain (i.e. the ratio of vertical displacement to the shear band thickness) at the peak force was unchanged and the shear band size was proportional to \({D}_{50}\) [24, 25], which was the same between physical testing and (scaled) numerical simulation. By normalising the upward displacement by the root diameter \({d}_{\rm r}\) in Fig. 6a, it was found that the response of the 3 mm root under 12 kPa confining stress, which represented the largest root under the highest confining stress considered in DEM, was consistent with physical modelling especially during the process pre-peak. The peak force (pull-out capacity), and forces post-peak, differed from the measured data by no more than 10%. Expanding the scale over the first 4 mm of displacement (at prototype scale) can be seen in Fig. 6b. The light grey lines in Fig. 6 represented the original laboratory testing data, where there was an initial low stiffness due to load string compliance in tensile testing using the Instron with wedge-type tensile grips. To eliminate this effect, the curves were leftward shifted (black solid lines) for comparison with DEM results. As the peak force describes the pull-out resistance, to save runtime, subsequent simulations were stopped shortly after peak force had reached. The pull-out curve of the 1.6 mm root under 1.2 kPa confining stress, representing the smallest root under the lowest stress level, is shown in Fig. 6c. While the secant stiffness was higher than that in the laboratory test at larger displacement pre-peak, the two curves matched in capacity. As depicted in Fig. 6d, all eight DEM validation tests showed similar peak forces to the physical modelling with a slight underestimation of the uplift force for 3 mm root cases under 1.2 kPa confinement.

Fig. 6
figure 6

Validation of sinker roots modelling in DEM expressed at prototype scale: a complete force–displacement curve of 3 mm root under 12 kPa; b force—displacement curve of 3 mm root under 12 kPa over the first 4 mm of displacement; c force—displacement curve of 1.6 mm root under 12 kPa over the first 4 mm of displacement; d summary of peak forces from all simulations

Full size image

Results of FEM validation against 12 mm roots under the four different confining stresses are shown in Fig. 7, with physical test curves again shifted to eliminate the effect of load string compliance. Simulations were stopped after the peak force had been reached and the mesh was becoming more severely distorted. The difference of peak resistance between laboratory tests and numerical simulations was less than 10%, indicating that using a zone of soil to model the shear band for larger roots could generally capture the root uplift capacity. However, the displacements at which peak force were reached in FEM was generally smaller than those in the laboratory testing.

Fig. 7
figure 7

Validation of 12 mm sinker roots in FEM expressed at prototype scale

Full size image

4 Parametric study of particle size effects

Particle size effects on sinker roots during both upward and downward movements were investigated using the validated numerical models. The principal of these simulations was to scale down the length, diameter and cover depth (denoted as c, the distance between top of the sinker root and the soil surface) of the roots by 1:N and scale up the gravitational acceleration by N to maintain the same prototype virtual model (in both length and diameter) and effective confining stress level but without scaling particle size. By doing so, the number of soil particles surrounding the roots was intentionally decreased. This physically represents choosing to model a particular sinker root at different scale factors on the centrifuge, therefore having relevance to the design of physical model tests, but also is a convenient way of changing only the \({d}_{\rm r}/{D}_{50}\) ratio for a given root while using the same underlying DEM particle assembly. Following conventional centrifuge scaling laws [31], force at prototype scale could be obtained by multiplying the force of the numerically simulated scaled root by N2. The uplift capacity at prototype scale should be the same in all cases if no particle size (scale) effects are present. The prototype diameter (Ndr) and cover depth were selected to be 60 mm and 130 mm, respectively (\({c/d}_{\rm r}\)= 2.2), according to a Pinus pinaster tree root system, while diameters of scaled roots (\({d}_{\rm r}\)) investigated were based on the size of individual roots in the 1:20 centrifuge model of this root system (\({d}_{\rm r}\)= 3, 2.4, 2 and 1.6 mm) reported in Zhang et al. [45], with additional smaller \({d}_{\rm r}/{D}_{50}\) ratios (< 10) also considered (\({d}_{\rm r}\)= 1.2 and 0.8 mm), as shown in Table 4. Drained pull-out was considered under fully saturated soil conditions in the prototype; to avoid explicitly modelling water in DEM, a lower g-level, \({N}_{\rm d}\) (i.e. smaller than the scale factor N) was applied to the dry particle arrangement so that the effective stress was representative of that at prototype scale in the saturated soil, following Eq. (4):

Table 4 Sizes of sinker roots for parametric study (upward movement = SU; downward movement = SD)
Full size table
$${N}_{\rm d}=\left(1 - \frac{1}{{G}_{\rm s}}\right)N$$
(4)

A series of preliminary uplift tests of sinker roots was conducted. A 3-mm-diameter root with 12 mm length was uplifted at 12g, while a 3 mm root with 3.2 mm length was uplifted at 45g. Note that the effective stress distributions along these two roots were identical and the only difference was the root length (Fig. S1). It was found that the peak force of the 12-mm-long root was approximately 3.75 (= 12/3.2) times higher than that of the 3.2-mm-long root, suggesting that there were no further grain-size effects associated with the reduction of the root length from 12 to 3.2 mm. For the parametric study, the length of 3-mm-diameter sinker roots was determined to be 12 mm in length (at model scale, representing a field scale sinker root 60 mm in diameter and 240 mm in length), which guaranteed that no further particle size effects would result from reduced length on the smallest sinker root investigated in this study (i.e. the smallest 0.8-mm-diameter root became 3.2 mm in length).

During uplift only the shaft resistance was active. To simulate the shaft resistance only, a frictionless cylinder was assembled on the top of the sinker root throughout the cover depth (similar to Fig. 5a, with the zone of surcharge weight changed to the soil with cover depth c) in uplift cases.

The uplift velocity was set to be 4 mm/s for all simulations to ensure quasi-static conditions in DEM. The methods to generate the sample were the same as those described in the validation study. The domains were set to be large enough to avoid boundary effects, being \({{D}_{\rm s}/d}_{\rm r}\) = 2, \({{B}_{\rm s}/d}_{\rm r}\) = 9 for pull-out (informed by the validation study) and \({{D}_{\rm s}/d}_{\rm r}\) = 8, \({{B}_{\rm s}/d}_{\rm r}\) = 9 for push-in (\({{D}_{\rm s}/d}_{\rm r}\) increased to minimise boundary effects on base resistance at the tip of the roots). Corresponding FEM simulations were also conducted using these same dimensions for comparison, and the methodologies were also similar to those introduced in the validation study.

4.1 Pull-out

As shown in Fig. 8a, the prototype uplift capacities of sinker roots with different g-levels applied were all approximately 0.15 kN in FEM. In DEM, there existed a negative linear relationship between the pull-out force and \({d}_{\rm r}/{D}_{50}\), indicating the significance of particle size effects. The uplift force of a 0.8 mm sinker root (\({d}_{\rm r}/{D}_{50}\) = 6) was approximately twice that of a 3 mm root (\({d}_{\rm r}/{D}_{50}\) = 21), indicating that the decrease in \({d}_{\rm r}/{D}_{50}\) from 21 to 6, when considering the same prototype, would overestimate the vertical root pull-out capacity by 100%. However, when \({d}_{\rm r}/{D}_{50}\) reached 21, the results from DEM converged to those from FEM.

Fig. 8
figure 8

Uplift force of sinker roots with a the same prototype dimension but different scale factors N; b the same length but different diameters

Full size image

Since no particle size effects resulted from the reduced root length, the pull-out capacity of roots with the same length (12 mm) and stress level but different diameters could also be obtained by multiplying by N only to the measured force at model scale (Fig. 8b). It can be seen that in FEM, the uplift force, mainly from the root shaft resistance (\({Q}_{\rm s}\)), increased linearly with the increment of \({d}_{\rm r}/{D}_{50}\), which was reasonable according to Eq. (5) as the root surface area (\({A}_{\rm s}\)) was increased linearly in proportion to \({d}_{\rm r}\):

$$Q_{\rm s} = \pi d_{\rm r} \int_0^{l_{\rm r} } {\tan \delta \sigma_n^{\prime} \left( z \right){\text{d}}z}$$
(5)

where \({\sigma }_{n}^{{\prime}}\left(z\right)\) is the normal effective stress applied on the shaft of the root, and \(\delta\) is the root–soil interface friction angle (equal to 25° in this study). The pull-out capacity obtained in DEM was significantly higher, especially when \({d}_{\rm r}/{D}_{50}\) was lower than 15. Note that the simulation of larger roots (i.e. \({d}_{\rm r}/{D}_{50}\) > 21) was not conducted due to the limited size of the central non-scaled zone considered in the DEM.

4.2 Push-in

Unlike the uplift cases, there was no observable peak resistance as the soil reaction continued increasing with the root downward displacement (i.e. increased burial and effective confining stress). In this case, the prototype force at 10 mm prototype displacement (16.67%\({d}_{\rm r}\)) was extracted from curves of roots with different diameters for comparison (Fig. 9). It can be seen that the decrease of root diameter from 3 mm (\({d}_{\rm r}/{D}_{50}\) = 21) to 0.8 mm (\({d}_{\rm r}/{D}_{50}\) = 6) resulted in 40% higher push-in resistance, suggesting non-negligible particle size effects. Given that the force in these cases was contributed to by both the shaft and base resistances and the shaft resistance was largely affected by particle size effects (Fig. 8), the base bearing resistance was extracted separately for further investigation. To eliminate the contribution of shaft resistance and consider only the base resistance in FEM, an interface element with strength reduction factor \({R}_{\rm inter}\) = 0 was applied in this case rather than using soil elements to model the shear band. As shown in Fig. 9, while the change in base resistance was almost negligible, the changes of the shaft resistance were substantial due to the previously described particle size effects. Toyosawa et al. [44] reported particle size effects on the bearing capacity of model footings with \(c/{d}_{\rm r}\) equal to 0 and 0.5 when \({d}_{\rm r}/{D}_{50}\) was lower than 50 and 33, respectively. In this study, the limiting \({d}_{\rm r}/{D}_{50}\) = 15 for \(c/{d}_{\rm r}\) = 2.2, which might suggest that particle size effects are more pronounced at lower cover depths.

Fig. 9
figure 9

Push-in force of sinker roots with different scale factors

Full size image

4.3 Deformation mechanisms

Comparison of the influenced zone of the root (in terms of displacement) during an uplift simulation can be seen in Fig. 10. The sinker root is shown in its final position at the end of the uplift phase (i.e. upward displacement was 0.7 mm at model scale) when the soil around the root reached maximum mobilisation of strength. Particles with vertical displacement smaller than 1% of the root displacement are not shown. It is noted that although the 4 mm root was 5 times wider than the 0.8 mm root, the sizes of the influenced zones were quite close in lateral extent, indicating the localised interface mechanism and appropriate application of the fixed-width ‘shear band’ zone in the FEM simulations. This also implies a thinner root can obtain a higher uplift force per diameter than a thicker root as a higher proportion of surrounding soil is mobilised.

Fig. 10
figure 10

Soil vertical displacement field of sinker root at the end of uplift phase: a 0.8 mm root; b 4 mm root (all double-headed arrows have the same length)

Full size image

The deformation mechanism from sinker root push-in simulations is illustrated in Fig. 11. As the slenderness ratio (\({l}_{\rm r}/{d}_{\rm r}\)) was the same in these simulations, the roots are shown at the same prototype size. Figure 11 shows that the smaller root affected a proportionally larger amount of soil around the shaft compared to its diameter, explaining the trend of shaft resistance contribution in Fig. 9. Particles above the root also moved downward with the root; however, their contribution to the resistance was negligible. In terms of particles beneath the root, the displacement field had a similar size compared with the root diameter, consistent with the observation in Fig. 9 that a change in scale factor did not show a large influence on the root end-bearing capacity.

Fig. 11
figure 11

Deformation mechanism for sinker root push-in cases: a 0.8 mm root; b 2.4 mm root

Full size image

5 A simple model for uplift resistance incorporating particle size effects

The results depicted in Fig. 8b suggested that for sinker roots with the same length and confining stress, surprisingly, the uplift force from DEM did not increase linearly with the increase of the diameter. To seek an explanation, a series of complementary further FEM and DEM simulations were performed (Table 5). In these simulations, the length of roots, cover depth (\(c\)) and g-level were kept constant at 12 mm, 0 mm and 120g, respectively, with the diameter (\({d}_{\rm r}\)) changed only. To obtain a wider range of \({d}_{\rm r}/{D}_{50}\) (particularly at the higher end), a new DEM test chamber with a 9-mm-diameter central non-scaled zone was generated according to the methodology described previously. The reason for using zero cover depth and 120g was to reduce runtime to compensate for the larger central zone and avoid extremely low stress levels.

Table 5 Sizes of sinker roots for further uplift simulations with fixed root length
Full size table

As shown in Fig. 12, the trend of model force in the range of \({d}_{\rm r}/{D}_{50}\) smaller than 20 was similar to that in Fig. 8b, while at larger \({d}_{\rm r}/{D}_{50}\) the DEM tracked the linear relationship suggested by the continuum model (FEM). Peak average shear stresses (τpav), derived as the peak applied tension force divided by the shaft surface area, were determined and are plotted against \({d}_{\rm r}/{D}_{50}\) in Fig. 13. It is noted that using the following equation would largely underestimate τpav:

$${\tau }_{\mathrm{pav}}= {\sigma }_{n0}^{{\prime}}\mathrm{tan}\delta$$
(6)

where \({\sigma }_{n0}^{{\prime}}\) is the average initial horizontal effective stress (i.e. value at mid-depth, approximately 5 kPa as obtained from DEM). This was because the dilation of the interface soil resulted in a significant increase in lateral effective stress on the sinker root, which in turn influenced the mobilised shear stress along the root–soil interface. The amplification of root-normal effective stress (\({{A}_{n}= \sigma }_{n}^{{\prime}}/{\sigma }_{n0}^{{\prime}}\), where \({\sigma }_{n}^{{\prime}}\) is the current lateral stress) along the root for some typical diameters at 0.4 mm upward displacement (approximately where the peak resistance occurred) are presented in Fig. 14, where the root was uniformly divided into 10 segments, with \({A}_{n}\) marked at the mid-depth of each segment. \({A}_{n}\) generally reached a peak value around the middle of the root and became significantly increased at low \({d}_{\rm r}/{D}_{50}\) values.

Fig. 12
figure 12

Dependence of peak uplift force on ratio of root diameter to mean grain size

Full size image
Fig. 13
figure 13

Dependence of average shear stress on ratio of root diameter to mean grain size

Full size image
Fig. 14
figure 14

Amplification of lateral stress from in-situ conditions along the displacing root

Full size image

The mobilised coefficient of lateral earth pressure K depicted in Fig. 15a was calculated using:

$$K= {\tau }_{\mathrm{pav}}/({\sigma }_{v0}^{{\prime}}\mathrm{tan}\delta )$$
(7)

where \({\sigma }_{v0}^{{\prime}}\) is the initial vertical stress at mid-depth of the root. It can be seen that the lateral earth pressure coefficient was always larger than \({K}_{0}\) (at rest condition; \(K\) > 3.5 \({K}_{0}\)), when using the continuum peak friction angle of 41° (Table 1). Previous studies [23, 27] have used coefficient of passive earth pressure \({K}_{\rm p}\) (with the same continuum peak friction angle) to quantify the maximum soil dilation during the root pull-out; however, this was seen to only be a reasonable approximation at low \({d}_{\rm r}/{D}_{50}\), e.g. for small diameter roots (Fig. 15a). As shown in Fig. 15b, the normalised coefficients of lateral earth pressure, i.e. \(K/{K}_{0}\), from DEM are consistent with those back-calculated from physical dry pull-out tests conducted by Mickovski et al. [28], where the average normal effective stress adjacent to the model roots was 0.8 kPa. Note that although a single interface friction angle of 21° was used in Mickovski et al. [28], a higher value of 41° was adopted for interpretation in this study based on interface shear tests from Mickovski et al. [28] and the known normal effective stress in the pull-out tests (Fig. 15c).

Fig. 15
figure 15

a Dependence of coefficient of lateral earth pressure on ratio of root diameter to mean grain size; b comparison of normalised coefficients of lateral earth pressure with Mickovski et al. [28]; c interface friction angle versus normal effective stress in Mickovski et al. [28]

Full size image

The significant increase of lateral earth pressure with the reduction of root diameter could be explained using a cavity expansion analogy [5]. Given that radial expansion \(\Delta r\) occurred near the root during the uplift as a result of dilation, Boulon and Foray [5] deduced the following approximate expression for the increase in normal effective stress (\(\Delta {\sigma }_{n}^{{\prime}}\)):

$$\Delta {\sigma }_{n}^{{\prime}}= \frac{4G}{{d}_{\rm r}}\Delta r = {k}_{n}\Delta r$$
(8)

where \(G\) is the equivalent linear shear modulus of the soil mass constraining the dilation and \({k}_{n}\) is the equivalent normal stiffness. In the DEM simulations, \(\Delta \sigma_n^{\prime}\) was obtained from the summation of the normal component of the particle-root contact forces and \(\Delta r\) was measured at 0.4 mm upward displacement via the average radial displacement of particles in contact with the shaft of the root. Knowing \(\Delta \sigma_n^{\prime}\) and \(\Delta r\), \({k}_{n}\) at this instant in each simulation was calculated according to Eq. (8). Peng et al. [32] fitted the normalised \(\Delta {\sigma }_{n}^{{\prime}}-{k}_{n}\) relationship from numerically simulated interface shear tests using the following equation:

$$\left( {\frac{{\Delta \sigma_n^{\prime} }}{{\Delta \sigma_n^{{\prime} *} }}} \right)^2 + \left( {1 - \frac{k_n }{{k_n^* }}} \right)^2 = 1$$
(9)

where \({k}_{n}^{*}\) of 10,000 kPa/mm was a reference normal stiffness and \(\Delta {\sigma }_{n}^{{{\prime}}*}\) was a reference increased normal effective stress to represent constant volume (CV) conditions. Due to the lack of interface CV tests at 5 kPa in this study, \(\Delta \sigma_n^{{\prime} *}\) used was estimated to be 20 kPa to fit the DEM simulations with diameter larger than 3 mm (corresponding \({k}_{n}\) < 1300 kPa/mm) well, though discrepancies were observed on smaller roots (corresponding \({k}_{n}\) \(\ge\) 1300 kPa/mm). This suggests that besides scale effects associated from cavity expansion, there were additional effects related to particle size.

Figure 16 expresses the number of particle contacts on each root at 0.4 mm displacement (~ peak force), and a significant linear relationship can be identified. The gradient of the linear fitted line was 1365.6 mm−1, which physically means the number of contacts per mm diameter. The fitted line has a vertical intercept of 634.9, implying an equivalent root diameter of 0.47 mm (calculated from 634.9/1365.6 ≈ \(3.3{D}_{50}\)) that always contributes to the shaft resistance. This suggests an increase of effective root diameter by \(3.3{D}_{50}\), such that by combining Eqs. 5, 8 and 9, the peak uplift force (\({F}_{\mathrm{up}}\)) arising from shaft friction can be modelled by:

$$F_{{\text{up}}} = \pi (d_{\rm r} + 3.3D_{50} )l_{\rm r} \left( {\sigma_{h0}^{\prime} + \Delta \sigma_n^{{\prime} *} \sqrt {1 - \left( {1 - \frac{4G}{{d_{\rm r} k_n^* }}} \right)^2 } } \right)\tan \delta$$
(10)
Fig. 16
figure 16

Number of contacts on the sinker root at 0.4 mm (model scale) upward displacement (~ peak uplift force)

Full size image

This equation is consistent with findings from previous studies, where Foray et al. [13] found that an increase of \({D}_{50}\) would result in higher pile uplift resistance. The additional term ‘3.3D50’ in the effective diameter \({(d}_{\rm r}+3.3{D}_{50})\) is significant when \({d}_{\rm r}\) is small but becomes negligible for a thicker root. As shown in Fig. 12, the consideration of particle size effects through the increased effective diameter and enhanced normal stress in Eq. (10) captures (1) the insensitivity of \({F}_{up}\) to particle size effects for \({d}_{\rm r}/{D}_{50}\) > 20; (2) the apparent plateau for 10 \(\le\) \({d}_{\rm r}/{D}_{50}\) \(\le\) 20 as particle size effects become more apparent; and (3) the significantly enhanced \({F}_{up}\) compared to a continuum model for \({d}_{\rm r}/{D}_{50}\) < 10. Note that this equation is valid for vertical root analogues embedded in HST95 sand, with \({d}_{\rm r}/{D}_{50}\) ranged between 3 and 57. Caution is advised when applying it to other granular soils with different particle size distributions (e.g. well-graded or gap-graded materials) when interacting with roots of different sizes for which further validation is warranted. Also, estimating realistic interface behaviour between real roots and soil might be more complicated in field applications due to the presence of root hairs and exudates, affecting the soil local to the root (rhizosphere) and controlling equivalent interface properties. Using ABS plastic analogues may to some extent lead to the underestimation of pull-out strength. Future work is required to fully understand how these factors affect root–soil interface properties and provide improved input parameters for this model.

6 Conclusions

This study has presented an investigation into potential particle size effects that may arise in the axial pull-out or push-in of vertically orientated plant roots, such as the sinker roots of trees. 3D DEM numerical simulations were validated against laboratory pull-out tests and were then employed to work as a ‘virtual centrifuge’ to replicate the behaviour of sinker roots of different diameters interacting with a granular bed and allow the ratio of root to particle diameter to be changed with the same prototype geometry and confining stress conditions. Results were compared with validated FEM simulations, where the soil was considered as a continuum.

It was found that a change in \({d}_{\rm r}/{D}_{50}\) did not significantly affect the end-bearing resistance. This implies that estimates of this element of resistance could be based on conventional continuum analytical approaches (e.g. end-bearing of closed-ended piles). In terms of the shaft resistance, significant particle size effects were observed, partly arising from the effects of cavity expansion in the surrounding soil enhancing normal effective confining stresses and partly from an increase in the effective diameter of the root by approximately \(3.3{D}_{50}\). A simple analytical model was developed for estimating shaft/pull-out resistance incorporating these effects which was able to capture three key regimes of behaviour observed in DEM simulations, defined by strong particle size effects for \({d}_{\rm r}/{D}_{50}\) < 10 and a transition to continuum response for \({d}_{\rm r}/{D}_{50}\) > 20.

Further research using different particle shapes, particle size distributions, and thus varying \({d}_{\rm r}/{D}_{50}\) through different values of \({D}_{50}\), along with corresponding physical tests, are required to generalise these findings. Nevertheless, the results provide useful guidance for estimating axial soil–root interaction properties for use in rooted soil models underpinning the use of vegetation in nature-based solutions such as slope stabilisation or in crop lodging analyses. They can also be used to inform the design of physical model (e.g. centrifuge) tests of root–soil systems in terms of the selection of appropriate scaling factors.