Mathematical Problems in Engineering

Mathematical Problems in Engineering / 2021 / Article
Special Issue

Numerical Methods in Geotechnical Engineering

View this Special Issue

Research Article | Open Access

Volume 2021 |Article ID 6623489 |

Hao Bai, Ruidong Li, Wubin Wang, Kang Xie, Xiang Wang, "Investigation on Parameter Calibration Method and Mechanical Properties of Root-Reinforced Soil by DEM", Mathematical Problems in Engineering, vol. 2021, Article ID 6623489, 18 pages, 2021.

Investigation on Parameter Calibration Method and Mechanical Properties of Root-Reinforced Soil by DEM

Academic Editor: Wangcheng Zhang
Received26 Nov 2020
Revised22 Feb 2021
Accepted26 Mar 2021
Published24 Apr 2021


The behavior of root-soil system has raised more and more attention in both ecological and geotechnical fields. In this study, a two-dimensional discrete element method is employed using PFC2D to simulate the root-reinforced soil. The root system is mimicked by chains of bonded discs, while the soil is modeled by granular particles. The tensile strength of the root is modeled by interdiscs’ bonding strength. Three laboratory tests were studied to calibrate the micromechanical parameters of DEM. Finally, direct shear tests on rooted soil are simulated to investigate the influence of different root characteristics on the root reinforcement effect.

1. Introduction

In recent years, the behavior of the root-soil system has raised more and more attention in both ecological [1] and geotechnical fields [24]. On the one hand, from a geotechnical point of view, roots of trees and other plants are able to stabilize the soil against various natural hazards [5] (e.g., landslide, seepage [6], and soil creep) and are considered to perform as soil reinforcement [7] when analyzing the stability of vegetated slopes [8]. More specifically, the roots across the potential slope sliding surface can provide additional shear strength of the soil [9, 10] by transferring the soil shear stress into root tensile forces via the friction surface between roots and soil [11]. On the other hand, from an ecological point of view [12], roots that anchor the upper structure in the earth provide plants with resistance to overturning when they are subjected to wind loads in windstorms. In addition, the external loads (e.g., lateral and rotational forces) that are transferred from the stem into tensile forces in windward roots are subsequently transmitted to soil [13].

Consequently, so as to investigate both the roots reinforcement effect on slope stability and the root anchorage capacity on plants stability, numerous studies have been completed to not only quantify the uprooting behavior and tensile strength of roots embed in soil (namely, the pull-out resistance) by in situ or laboratory pull-out test [14], but also characterize the shear behavior of rooted soil by direct shear test at root bundle [15] or root system scale using either real roots or root analogues [16]. The existing studies demonstrate that the performance of the root-soil system [17] depends on both ground traits (e.g., soil properties, matric suction, and confined pressure) and root traits which can be divided into geometrical characteristics (e.g., branching pattern, bundle type, diameter, length, and orientation) and mechanical properties (e.g., tensile strength and Young’s modulus) [18].

Nevertheless, these experimental methods are difficult to generate replicates and control multiple factors. To fill this gap, numerous numerical models based on FEM were utilized to simulate the uprooting process of plants. These numerical models enable the researchers to study the progressive breakage and slipping out of roots during uprooting. However, modeling such a process using FEM in a detailed manner is technically complex and time-consuming. In addition, it requires a good mesh quality to prevent convergence difficulties when dealing with larger soil strains problem during FEM simulation.

Due to these limitations of FEM in the simulation of the root soil interactions, a two-dimensional discrete element method [19, 20] is employed in this research using PFC to overcome these problems. Some scholars have studied the mechanical properties of root-reinforced soil through DEM method [21, 22]. However, previous investigations mainly analyzed the effect of root on soil strength and did not analyze the effect of root arrangement direction and number on the soil properties. The root system is mimicked by chains of bonded discs, while the soil is modeled by granular particles. Tensile strength of the root is modeled by interdiscs’ bonding strength. The shear strength increment obtained from the direct shear test would be investigated. The primary objective of this research is to define different parameters related to soil, roots, and root-soil interactions correctly in the discrete element model and evaluate the root reinforcement effect on soil through simulation of simple laboratory tests.

2. Calibration of Micromechanical Parameters in DEM

In this section, we first simulate three laboratory tests to examine the microscopic parameters of root and soil models in DEM.

2.1. Tension Behavior of Single Straight Roots
2.1.1. Test Description

Similar tension tests can be done by using two separate steel plates to clip the root. As shown in Figure 1(c), the sandpapers and rubber were utilized to prevent the sliding of the root during tension. A motor device locked with the steel plates was employed to pull the root.

2.1.2. Modeling a Single Straight Root in DEM

In this section, single straight roots with different diameters and lengths were modeled by clustering a row of small spheres together (shown in Figure 2). The diameter of these small spheres is taken as the same value as the root diameter in this study. In DEM, the tensile behavior of the roots can be simulated by connecting those small spheres with parallel bonds which represent the stress-strain characteristic of a finite piece of cementitious material between two spheres [23]. The properties of parallel bond corresponding to real roots’ tensile strength and modulus will be firstly back-calculated assuming that the stiffness and strength along the individual root are isotropic and homogeneous.

The total force and moment associated with the parallel bond are denoted by F and M, with the convention that this force and moment represent the action of the bond on sphere B of Figure 2. If the maximum stress exceeds the corresponding bond strength, the parallel bond will break, and this corresponds to the root breakage.

2.1.3. Calibration of the DEM-Simulated Roots Parameters

The micromechanics parameters of the simulated single straight roots in DEM were first assigned according to the theoretical solutions relationship between the micromechanics parameters and the macro tensile properties as determined above. Then, the comparison between the corresponding results obtained from the DEM-simulated tensile test and the theoretical solutions would be done to verify the discrete element model.(1)Physical Model of a Single Straight RootThe real tap root system can be idealized as a single nonbranched straight root with the same circular cross sections along the length. Then, a physical model of the idealized single straight root can be obtained, as shown in Figure 3, which can be considered as a single spring with an equivalent stiffness K. In DEM, this single root can be modeled by joining a row of equal-radius spheres with parallel bonds. Each parallel bond can then be regarded as an individual elastic cylinder to mimic the physical root model.According to the simplified models, the stiffness of each parallel bond can be determined according to the elastic theory. The combination stiffness K of a group of springs is based on whether the collection of springs is acting in series or in parallel and is given byFor the special case in which all of the springs have the same stiffness, k, then the expression for the equivalent stiffness simplifies toWhen the root is under tension, only parallel bonds carry the load. However, when the root is under compression, both the parallel bonds and the contact springs will act together in parallel to carry the load.(2)Theoretical Relations between Parallel Bond and Root Tensile ParametersWhen a single straight root is under tensile force , assuming the root is a perfectly elastic spring, the following equations can be obtained:where is the deformation of the root; L is the initial length of the root; E is the Young’s modulus of the root; A is the cross-sectional area of the root; and K is the equivalent stiffness of the physical root model.If the root is mimicked by only a row of spheres, the following equation can be obtained:where is normal stiffness of the parallel bonds and is the diameter of balls.(3)DEM Simulation ResultsFour groups of root samples are simulated (as shown in Table 1). Tsimu and Ttheor are the tensile strength of numerical simulation and theory, respectively. The first group illustrated that the parallel bond strength is proportional to the tensile strength of the simulated roots. The second group shows that the parallel bond stiffness is proportional to Young’s modulus of the simulated roots. The third group demonstrated that with the same DEM parameters, the length would not influence the macroscopic tensile strength and modulus of the roots. The fourth group presented that, with the same DEM parallel bond stiffness and strength, the simulated Young’s modulus is proportional to the diameter of the root particles which match the theoretical equation.

Group no.Root dimensionsDEM parametersSimulated and theoretical solution
Length (m)Diameter (m)pb_stiff (N/m3)pb_str (MPa)Esimu (MPa)Tsimu (MPa)Etheory (MPa)Ttheory (MPa)





According to the simulated and theoretical results shown in Table 1, the macroscopic properties of the roots under tension load in DEM simulation are almost identical with the theoretical solutions. In addition, when the size of the root element (single sphere) is fixed (2 mm in diameter in this study), the simulated macroscopic properties of the roots (Young’s modulus and tensile strength) are only correlated with the microscopic mechanical parameters (parallel bond stiffness and strength).

Upon the simulation and analysis above, in the following simulation model, the microscopic parameters of the roots in DEM will be assigned according to the theoretical solutions between the bonding stiffness/strength and the root tensile modulus and strength. In addition, the value of the macroscopic root tensile parameters will be determined according to the experimental relationships with the root diameter obtained by Operstein [24].

2.2. Direct Shear Performance of Plain Soil

Before studying the rooted soils, simulations of direct shear test on soil samples without roots are needed to verify the applicability of the discrete element model, and the simulated results would also be valuable for comparison of shear resistance between nonrooted and rooted soils.

2.2.1. Test Description

The direct shear test has been widely used in the industry for testing sands and continues to attract research interest. The direct shear box is split horizontally into two halves. Each half contains a certain amount of soil during the test. Since the upper half and the lower half are separated from each other, they can move relatively if one-half of the box is under one horizontal force T. And the horizontal displacement will be caused, denoted by . On both the upper and lower surfaces, there is a porous stone pad. The usage of these two porous stone is to allow the drainage of water during the test. However, in DEM simulation, the pore water pressure is not considered. Consequently, in the laboratory test, the test sample is actually dry sandy soil. But in order to follow the convention of the test, we still apply the upper and lower porous stone to the dry sample in the test.

As it is the strain-controlled test, a shear rate of 1 mm per minute is applied to the upper shear box. During the test, the horizontal displacement can be observed by the horizontal dial gauge. And the resisting shear force is indicated by the horizontal proving ring. The height evolution of the sample (indicated as volume change in 2D case) is recorded by dial gauge which is used to measure the vertical displacement of the upper loading plate. Compared to the stress-controlled test, which has a uniform loading increment, the strain-controlled test has its advantage. As for dense sand, both the peak value of the shear resistance and the residual shear resistance (after failure) can be obtained.

2.2.2. DEM Model of Direct Shear Box

The simulated sample preparation scheme is the same as that in laboratory experiments. To begin with, an assemblage of spheres with very small diameters without overlapping was randomly distributed in the shear box which consisted of 8 rigid walls. Then, the spheres are expanded to their final size according to the desired void ratio.(1)Wall Model in DEMThe model of the simulated direct shear box consisted of an upper box and a lower box with the same height of 30 mm along the y-axis and the same width of 100 mm along the x-axis (see Figure 4). The top wall of the upper box is free to move vertically, while the side walls of the upper box could be translated horizontally. The lower box is fixed and all the walls do not interact with other walls but with contacted balls. These walls are rigid and do not deform during shearing.(2)Soil Particles in DEMThe initial state of the soil sample is generated by Radius Expansion Method. Firstly, a specified number of particles with artificially smaller radii are placed at random coordinates within the shear box. Then, the particles are expanded to achieve the desired porosity of the soil sample.(3)Initial Stress State by Servo ControlTo simulate the real case of a laboratory direct shear test in which the dead load (σn) was applied at the top of the upper box during the shear process, the vertical displacement of the top wall of the upper box is servo-controlled to maintain constant vertical stress. Firstly, the applied overburden pressure is computed according to the stress states of the soil particles that are contacted with the top wall. The function iterate is called to obtain the requested stresses within a specified tolerance of 0.5%.(4)Stress Monitoring during ShearingThe ratio of the mean static unbalanced force to the mean contact force within the sample can be utilized to illustrate whether the shearing speed of the test is low enough to maintain the quasistatic state. Additionally, the stress acting on the top wall of the upper box can also be an indicator. The shear stress at the shear band between the upper box and the lower one is determined by computing the mobilized horizontal forces acting on the wall.

2.2.3. Simulation Results Interpretation

The simulation results are presented in this part. The actual particle friction will be determined after loading the overburden pressure and before the shear process. A series of numerical simulations will be conducted on dense soil samples in which real interparticle friction coefficient (µ) and normal stress (σn) were varied to examine their effects on the shear banding behavior. In this investigation, assemblies are prepared by assigning a relatively low initial porosity of 10%.(1)Influence of Overburden PressureTypical stress-displacement relations from simulations of dense particle assemblies under different overburden pressure are shown in Figure 5(a). The simulation results show that higher overburden pressure gives higher shear stress during the shearing process.(2)Influence of Interparticle Friction CoefficientsFigure 5(b) illustrates the stress-displacement curve of which different interparticle friction coefficients (µ) were assigned to the soil particles. Microscopic interparticle friction coefficient (µ) strongly influences the stress-strain relations. With larger interparticle friction values, the shear stress increase with a stiff slope until a peak value is reached followed by a softening to a limit residual shear stress.

The peak values of the shear stress are recorded as the ultimate shear strength of the soil sample. As shown in Figure 6(a), the peak shear stresses are plotted versus the overburden pressures, and the linear regressions are conducted for each sample with different friction coefficients. Then, the apparent friction angles of each soil sample are calculated according to the slope of the straight line. As illustrated in Figure 6(b), the apparent friction angles for each soil sample increase with the interparticle friction coefficients which describe the microscopic friction angle between two individual soil grains.

2.3. Interfacial Friction Characteristic of Rooted Soil

In order to determine the appropriate micromechanical parameters in DEM to simulate the interfacial interaction between root and surrounded soil particles, a laboratory lateral pull-out test has been accomplished.

2.3.1. Test Description

A typical laboratory pull-out test of single straight roots was conducted in a relatively small box with dimensions of 444 mm wide × 355 mm long × 400 mm deep. Due to the lack of information, in this study, the pull-out test done by Schwarz [25] was considered as a reference. A uniformly distributed vertical pressure was activated on the top of the root-soil sample. The pull-out load would be applied by using a hydraulic jack which will react against the rigid structure. The pull-out displacements can be measured by a Linear Variable Differential Transformer (LVDT).

2.3.2. Simulation Scheme and Process

To simplify the model, only one single root is simulated within a relatively small pull-out box with a dimension of 60 mm height and 100 mm width. The process of sample preparation is also the same as the experiments conducted in the laboratory .Step 1. Preparation of the sample in the lower half-box.As shown in Figure 7, at the first step, half of the spheres (2 mm on average) were randomly generated within the whole box without overlapping. A relatively large density (1e6 kg/m3) and large gravity (98.1 m/s2) are assigned to all balls in order to speed up the simulation process. Then, all balls were directly deposited in the pull-out box and cycled to equilibrium under a changing density and gravitational acceleration which was reduced gradually to 2.7e3 and 9.81 m/s2, respectively. After that, the balls outside the lower half box are deleted in order to place the root.Step 2. Generation and placement of the root.At the second step, a row of spheres with the same diameter as the simulated root is then generated. The effective length of the simulated root is the distance between the center of the second ball to the center of the ended ball (the first ball is outside the pullout box and is considered as part of the pulling system).Step 3. Preparation of the sample in the upper half-box.At the final step, the residual soil was randomly generated by almost the same scheme in the first step. At this time, another half-spheres are randomly generated above the root without overlapping.

In order to analyze the influence of different parameters in DEM on the pull-out behavior of the simulated root, parametric studies are conducted following the scheme shown in Table 2.

Simulate groupRoot dimensionsRoot physical para.Soil physical para.
Length (m)Diameter (m)Eroot (MPa)Troot (MPa)fricroot (no unit)Esoil (MPa)fricsoil (no unit)σvertical (kPa)

(1) Soil friction0.080.0041e111e70.02.5e80.325

(2) Root friction0.080.0041e111e70.12.5e80.625

(3) Vertical pressure0.080.0041e111e70.02.5e80.350

(4) Root modulus0.080.0045e111e70.02.5e80.325

2.3.3. DEM-Simulated Result

The maximum friction mobilized along the interface of the root and soil can be obtained by dividing the normal force with shear force. The total normal force acting on the root balls can be obtained according to the following equation:where is the equivalent total normal force acting on the root balls; is the vertical stress; is denoted as the effective length of the root embedded in the soil, which can be obtained by tracing the displacement of the last ball at the tail of the whole root. Thus, the maximum normal force can then be obtained by substituting the length of the root at the largest strain into the above equation. Then, the friction coefficient can be obtained:

Figure 8(a) demonstrates the influence of overburden pressure on the pull-out force. The resultant friction under each overburden pressure is shown in Figure 8(b).

As can be seen from Figure 8(b) that, when the friction coefficient of the root ball is 0 and the friction coefficient of soil particles is 0.3, the friction angle of the root soil interface is around 22.7 degrees.

As shown in the figure above, with the increase of the root friction coefficient (from 0.1 to 0.6), the pull-out force increases a little bit. Figure 9(b) above demonstrates the increment of pull-out force when the soil friction coefficient increases from 0.3 to 0.9. Comparing the two figures above, it is indicated that the soil particles friction coefficient plays a more important role in root-soil interaction and governs the pull-out behavior in DEM model. Figure 9(c) is obtained by varying the stiffness of the parallel bond. It is obvious that with a larger parallel bond stiffness, the peak pull-out force would be mobilized earlier at relatively smaller strain. Figure 9(d) is applied to the soil. This indicates that a higher pull-out force will be mobilized in the soil with a larger cohesion.

3. Study on Shear Behavior of Root-Reinforced Soil

In this section, direct shear tests on rooted soil are simulated in order to investigate the influence of root number and root inclination on root reinforcement effect.

3.1. Description of Laboratory Direct Shear Tests on Rooted Soil

Operstein and Frydman [24] carried out the direct shear tests on both samples of chalky soil and chalky soil reinforced by roots. Single shear tests were carried out using the round pipe samples. Two 300 mm square aluminium blocks with a 200 mm dia. Ball bearing tracks minimized friction between the blocks during shearing. The soil-filled pipe was inserted vertically through the central hole such that the precut pipe section on which shearing was to be performed was located at the interface of the blocks. These shear planes were located at depths of between 50 and 750 mm below the soil surface.

The vertical force normal to the shear plane was provided by the weight of the overlying soil, and in some tests, additional dead loads were applied to the surface of the sample. Prior to shearing at a particular level, the clamp band was removed, the section of pipe above the cut was slightly raised and locked to the top of the upper aluminium block (to eliminate friction between the upper and lower surfaces of the cut), and the inner nylon wrapper was sliced.

Shearing was applied by the hydraulic jack (item 7), acting through a load cell (item 6), to the upper aluminium block, while the lower block was held in place. The horizontal displacement of the upper block was measured by two dial gauges (item 9).

The direct shear tests were executed to a maximum displacement of 40 mm (the limit of the apparatus), at a constant shear displacement ratio of about 0 : 1 mm = min, allowing a test to be completed in a 6 ± 7 h working day. Following shearing, the plant roots were exposed and their dimensions were measured.

3.2. Simulation Scheme and Process

Firstly, the procedure of sample preparation is demonstrated in the following steps. Both the soil and roots are generated layer by layer .Step 1. Preparation of the sample for the first layer.According to the length (l_root) and diameter (d_root) of the root, the root soil system is divided into n layers (n=l_root/d_root). As demonstrated in Figure 10, the first layer of the root balls is generated at first. Then, the soil particles are randomly distributed and then deposited under the gravity reduced from 98.1 m/s2 to 9.81 m/s2. After the whole system was cycled to equilibrium, the soil particles beneath the roots would all be deleted.Step 2. Preparation of the sample for the second layer.The second step follows the same process as the first step to generate the second layer of the root soil system. It should be emphasized that, during the generation process, the root balls are all fixed to prevent unexpected displacement during the deposition of the soil particles.Step i. Preparation of the sample for the rest layers.Loops Step 2 to generate the rest layers of the root soil system, as shown in Figure 11.Step n. Generation of the top wall.Finally, the whole system is cycled to equilibrium, and the top wall of the upper shear box is generated to seal the whole sample.

3.3. DEM-Simulated Result
3.3.1. Progressive Mobilization of Root Tensile Force

Before shearing, due to the process of activating overburden pressure, part of the vertical stress is sustained by the roots; thus, compression stress is induced through the root body (shown in Figure 12). Thus, when shear strain is 0, the induced normal stress acting on the shear plane is reduced due to the presence of roots, as shown in Figure 13(a).

During the initial shearing stage (shear strain smaller than about 0.01), since the initial normal stress acting on the shear plane of the rooted soil is smaller than that of nonrooted soil (shown in Figure 13), the initial shear stress mobilized in rooted soil is much smaller than that nonrooted one (shown in Figure 13). Then, with furthermore shear displacement, as shown in Figure 13, normal stress acting at the shear plane of the soil sample increases gradually due to the decrease of compression loading balanced by roots (during this period, the roots are still under compression loading, shown in Figure 12). Once the shear displacement went beyond a certain level, the tensile force began to be mobilized in roots bodies (see Figure 12(c)). Since almost all the compression loadings are sustained by soil, normal stress balanced by soil in the rooted case became larger than the nonrooted case at this stage (shown in Figure 13).

Moreover, with further shearing, since larger tensile forces are induced due to extraction of root (see Figure 12(d)), normal stress acting at the middle of the shear plane continues to be larger (see Figure 13). In this stage, shear stress in rooted soil became much larger than in nonrooted case, as demonstrated in Figure 13. In addition, the dilation behavior of the soil is also limited by the presence of root reinforcement. It is shown in Figure 14 that much smaller volumetric change is measured in rooted soil than in nonrooted one.

From these observations from DEM simulation results, a conclusion can be made that the effect of root becomes significant only when the shear strain is larger than the critical value, denoted by εcrit. Additional shear resistance provided by root is mobilized gradually due to the progressive development of root tensile force. Normal stress activated at the shear plane for rooted soil was smaller than the nonrooted case at relatively small strain, which results in a decrease in shear resistance when shear strain is smaller than εcrit. Then, the normal stress increases gradually and even goes beyond the value in the nonrooted case since the vertical component of the root tensile force gives additional vertical loading to the soil nearby. For frictional granular soils, this critical value ε might be dependent on root-soil interface property and the modulus of soil and root. In this simulation, with soil stiffness of 1e6 Pa, root tensile modulus of 1e8 Pa, root friction coefficient 0.5, and soil friction coefficient 0.3, the critical strain value is about 0.025.

3.3.2. Influence of Roots Number

For frictional granular soil simulated in this study, the shear resistance of the rooted soil is significantly influenced by the number of roots. According to the shear displacement-shear stress curve for different numbers of roots, the influence of roots number on reinforcement effect is quite different at the different shearing stages.

As can be seen from Figure 15, when shear strain is smaller than about 2.3 mm, the soil reinforced with larger root numbers tends to have lower shear resistance. The reason can be obtained according to Figure 16. At relatively low shearing strain, normal stress sustained by soil at shear plane decreases with the increment of roots since a larger number of roots would balance much more normal loadings.

However, when the shear displacement is larger than εcrit (around 0.023 in this case), all rooted soil exhibit higher shear resistance since tensile forces began to be mobilized in the root body. However, in this stage, since only part of roots are under tension force, a larger number of roots do not exhibit higher shear resistance until more roots are under tension. Then, with further shear strain, more roots start to be mobilized progressively. When shear displacement is larger than around 3.2 mm and 4.5 mm, respectively, larger shear resistance continues to be mobilized in soil sample with larger root numbers. Moreover, the dilation behavior of rooted soil is limited by roots number. As demonstrated in Figure 17, when the roots’ number is getting larger, the rooted soil tends to dilate later and exhibits a much smaller dilation angle.

As can be seen from Figure 17, in general, the volumetric strain first decreases, then increases when it reaches the minimum, and finally tends to be stable. With the increase of the number of roots, the shear strain corresponding to phase transition state increases and the dilatancy angle decreases. When there are 8 roots, the sample will only show shear shrinkage. The results show that the root reinforcement makes the stress in the specimen more uniform, thus restraining the dilatancy.

3.3.3. Influence of Different Roots’ Inclined Angle

In order to investigate the influence of root inclined angle, roots with different inclinations (45, 60, 75, 90, 105, 120, and 135 degrees) are simulated in a shear box. As can be seen from Figure 18, generally, roots with an inclined angle larger than 90 degrees tend to have smaller additional shear resistance than those with an inclined angle smaller than 90 degrees.

To reveal the reinforcement mechanism, normal stress acted at the shear plane of the soil sample is measured by the measuring circle in PFC2d. As shown in Figure 19, at the beginning of shearing, rooted soils with different inclined angles all have smaller initial normal stress acting on the shear plane.

Then, with furthermore shearing strain, the normal stress measured in soil shear plane for those samples with roots’ inclination angles smaller than 90 degrees tends to increase rapidly. The reason can be inferred that the vertical component of the tensile loading is transferred from the root body to the nearby soil as an additional normal reinforcement. Thus, larger vertical stress is observed at the shear plane.

For roots with inclined angles larger than 90 degrees, normal stress measured in the shear plane tends to decrease because more normal stress is balanced by root body; when the root orientation is opposed to the shearing direction, more compression loadings are sustained by root body.

As can be seen from Figure 20, for roots with inclined angles larger than 90 degrees, the roots are almost parallel to the contact force chain of the soil sample, in this case, according to limit equilibrium analysis, normal loadings acting at the shear plane of the soil are reduced due to the increment of the normal component of compression force mobilized in root body. Thus, although the horizontal component of root compression force tends to provide additional shear resistance, the reinforcing effect is not that obvious as compared to those with inclined angles smaller than 90 degrees in which the vertical component of the root tensile force tends to provide additional normal loading to the nearby soil.

4. Conclusions

In this research, a series of discrete element modeling of experimental tests on root soil systems are carried out by using PFC2D software. Firstly, three laboratory tests are simulated to calibrate the microscopic parameters of root and soil in DEM model including tension test on a single straight root, direct shear test on plain soil, and lateral pull-out test of single straight root embedded in soil matrix. Then, direct shear tests on rooted soil are simulated considering different root inclination angles and root numbers. The major conclusions can be made as follows:(1)When the size of the root element (single sphere) is fixed (2 mm in diameter in this study), the simulated macroscopic properties of the roots (Young’s modulus and tensile strength) are only correlated with the microscopic mechanical parameters (parallel bond stiffness and strength).(2)The friction angle of the root soil interface is determined by both the friction coefficient of the root body, but also determined by the friction coefficient of grain particles. In this study, with root friction of 0 and soil particles friction coefficient of 0.3, the friction angle of the root soil interface is around 22.7 degrees. The roots mobilize the tensile strength progressively from the pulled end to the anchored end, with the stretching of the root occurring before the peak pull-out force and contraction of the root occurring after the peak.(3)The overall root pullout behavior is dependent on three factors: the root behavior, soil behavior, and the interface properties.(4)As for the direct shear test on rooted soil, the root reinforcement effect is mobilized gradually with the increment of shear strain. Comparing to plain soil, the root reinforcement effect is only significant when shear strain is larger than the critical value. In this study, for soil particle stiffness of 1e6 and root tensile stiffness of 1e8, the critical strain is around 0.23. In addition, the dilation behavior of the soil sample is also limited due to the presence of roots.(5)Both root inclination and root number have a significant influence on shear stress-strain curves. For roots with different roots’ number, at the initial shearing stage, the rooted soil with larger roots’ numbers tends to have smaller shear resistance since more vertical loadings are sustained by root body and result in the decrease of shear resistance when shear strain is relatively small. Then, with larger shear strain, more roots are under tensile loading which provides the soil sample with larger shear resistance.(6)As for root inclination angles, for root’s inclined angle larger than 90 degrees (the root orientation opposed to the shearing direction, roots are under compression loading during the test), although the horizontal component of root compression force tends to provide additional shear resistance, the vertical component of root compression force tends to decrease the vertical stress acting on the shear plane of the soil and results in a significant decrease in shearing resistance. For roots’ inclined angle smaller than 90 degrees, vertical normal stress acting on the shear plane increases with shear train due to the vertical component of the tensile loading and results in larger shear resistance measured from the shearing test.

Data Availability

The data used to support the findings of this study are included within the article.

Conflicts of Interest

The authors declare that they have no conflicts of interest regarding the publication of this paper.


  1. C. O’Loughlin, “The effect of timber removal on the stability of forest soils,” Journal of Hydrology (New Zealand), vol. 13, no. 2, pp. 121–134, 1974. View at: Google Scholar
  2. K. Abe and R. R. Ziemer, “Effect of tree roots on a shear zone: modeling reinforced shear stress,” Canadian Journal of Forest Research, vol. 21, no. 7, pp. 1012–1019, 1991. View at: Publisher Site | Google Scholar
  3. D. M. Bishop and M. E. Stevens, “Landslides on logged areas in southeast Alaska,” US Forest Service Research Paper NOR, vol. 1, 1964. View at: Google Scholar
  4. T. Endd and T. Tsuruta, “The effect of the tree's roots upon the shear strength of soil,” Forestry and Forest Products Research Institute, vol. 18, pp. 167–182, 1986. View at: Google Scholar
  5. D. H. Gray, Effects of Forest Clear Cutting on the Stability of Natural Slopes, University of Minnesota, Minneapolis, Minnesota, 1969.
  6. T. Manbeian, The Influence of Soil Moisture Section, Cyclic Wetting and Drying, and Plant Roots on the Shear Strength of a Cohesive Soil and Stability of Slopes, University of California, Berkeley, CA, USA, 1973.
  7. R. A. Jewell, Some Factors Which Influence the Shear Strength of Reinforced Sand, University of Cambridge, Department of Engineering, Cambridge, UK, 1980.
  8. T. H. Wu and A. Watson, “In situ shear tests of soil blocks with roots,” Canadian Geotechnical Journal, vol. 35, no. 4, pp. 579–590, 1998. View at: Publisher Site | Google Scholar
  9. D.-G. Lin, B.-S. Huang, and S.-H. Lin, “3-D numerical investigations into the shear strength of the soil-root system of Makino bamboo and its effect on slope stability,” Ecological Engineering, vol. 36, no. 8, pp. 992–1006, 2010. View at: Publisher Site | Google Scholar
  10. D. G. Lin, W. T. Liu, and S. H. Lin, “Estimating the effect of shear strength increment due to root on the stability of Makino bamboo forest slopeland,” Journal of GeoEngineering, vol. 6, no. 2, pp. 73–88, 2011. View at: Google Scholar
  11. S. B. Mickovski, A. Stokes, R. Van Beek, M. Ghestem, and T. Fourcaud, “Simulation of direct shear tests on rooted and non-rooted soil using finite element analysis,” Ecological Engineering, vol. 37, no. 10, pp. 1523–1532, 2011. View at: Publisher Site | Google Scholar
  12. S. B. Mickovski, P. D. Hallett, M. F. Bransby, M. C. R. Davies, R. Sonnenberg, and A. G. Bengough, “Mechanical reinforcement of soil by willow roots: impacts of root properties and root failure mechanism,” Soil Science Society of America Journal, vol. 73, no. 4, pp. 1276–1285, 2009. View at: Publisher Site | Google Scholar
  13. L. Peng, K. Shengzu, Z. Wenwei et al., “Numerical simulation on single perpendicular root reinforcement mechanism in slope stability,” in Proceedings of the 2011 International Conference on Electric Technology and Civil Engineering (ICETCE), pp. 508–511, IEEE, Lushan, China, April 2011. View at: Google Scholar
  14. T. H. Wu, W. P. McKinnell III III, and D. N. Swanston, “Strength of tree roots and landslides on prince of wales Island, Alaska,” Canadian Geotechnical Journal, vol. 16, no. 1, pp. 19–33, 1979. View at: Publisher Site | Google Scholar
  15. S. B. Mickovski, A. G. Bengough, M. F. Bransby, M. C. R. Davies, P. D. Hallett, and R. Sonnenberg, “Material stiffness, branching pattern and soil matric potential affect the pullout resistance of model root systems,” European Journal of Soil Science, vol. 58, no. 6, pp. 1471–1481, 2007. View at: Publisher Site | Google Scholar
  16. M. Yang, P. Défossez, F. Danjon, and T. Fourcaud, “Tree stability under wind: simulating uprooting with root breakage using a finite element method,” Annals of Botany, vol. 114, no. 4, pp. 695–709, 2014. View at: Publisher Site | Google Scholar
  17. L. J. Waldron, “The shear resistance of root-permeated homogeneous and stratified soil,” Soil Science Society of America Journal, vol. 41, no. 5, pp. 843–849, 1977. View at: Publisher Site | Google Scholar
  18. N. S. Nilaweera and P. Nutalaya, “Role of tree roots in slope stabilisation,” Bulletin of Engineering Geology and the Environment, vol. 57, no. 4, pp. 337–342, 1999. View at: Publisher Site | Google Scholar
  19. F. Bourrier, F. Kneib, B. Chareyre, and T. Fourcaud, “Discrete modeling of granular soils reinforcement by plant roots,” Ecological Engineering, vol. 61, pp. 646–657, 2013. View at: Publisher Site | Google Scholar
  20. S. Frydman and V. Operstein, “Numerical simulation of direct shear of rootreinforced soil,” Proceedings of the Institution of Civil Engineers - Ground Improvement, vol. 5, no. 1, pp. 41–48, 2001. View at: Publisher Site | Google Scholar
  21. M. J. Jiang, Y. G. Zhu, and B. L. Xi, “Investigation into the soil-root composites using distinct element method,” in Proceedings of the International Conference on Discrete Element Methods, Santa Fe, NM, USA, September 2017. View at: Google Scholar
  22. Z. Mao, M. Yang, F. Bourrier et al., “Evaluation of root reinforcement models using numerical modelling approaches,” Plant & Soil, vol. 381, no. 1-2, pp. 249–270, 2014. View at: Publisher Site | Google Scholar
  23. L. Xu and Q. W. Ren, “Optimization inversion of the micro-properties of bonded particle model in PFC based on the asynchronous particle swarm algorithm,” Advanced Materials Research, vol. 1082, pp. 240–245, 2014. View at: Publisher Site | Google Scholar
  24. V. Operstein and S. Frydman, “The influence of vegetation on soil strength,” Proceedings of the Institution of Civil Engineers - Ground Improvement, vol. 4, no. 2, pp. 81–89, 2000. View at: Publisher Site | Google Scholar
  25. M. Schwarz, D. Cohen, and D. Or, “Root-soil mechanical interactions during pullout and failure of root bundles,” Journal of Geophysical Research: Earth Surface, vol. 115, no. F4, 2010. View at: Publisher Site | Google Scholar

Copyright © 2021 Hao Bai et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Related articles

No related content is available yet for this article.
 PDF Download Citation Citation
 Download other formatsMore
 Order printed copiesOrder

Related articles

No related content is available yet for this article.

Article of the Year Award: Outstanding research contributions of 2021, as selected by our Chief Editors. Read the winning articles.