Root reinforcement: continuum framework for constitutive modelling

The mechanical contribution of plant roots to soil strength has typically been studied at the ultimate limit state only. Since many geotechnical problems are related to serviceability, such as deformation of infrastructure, a new constitutive modelling framework is introduced. The rooted soil is treated as a composite material with separate constitutive relationships for soil and roots, and a comprehensive stress-strain relationship for the root constituent is presented. The model is compared to direct shear experiments on ﬁeld soil reinforced with gorse, grass and willow roots, as well as an existing root reinforcement model based on Winkler-spring supported beam theory. The results show that both the newly developed model and the beam-type model yield good predictions for the evolution of root-reinforced shear strength as a function of increasing shear displacements. Both successfully capture the large deformations required to reach peak reinforcement, the reduction in reinforcement due to root breakage and the presence of signiﬁcant reinforcement even after very large deformations, associated with root slippage. Since both ﬁbre and beam models only require physically meaningful input parameters, they can be useful tools to study the mobilisation of rooted soil strength and simulate the response of rooted soil in continuum-based numerical simulations.


Introduction
The mechanical reinforcement of soil by plant roots can improve the resistance against landslides or erosion processes (e.g.Stokes et al., 2009).This contribution has typically only been studied at the ultimate limit state, often in the context of slope stability, focussing on the maximum reinforcement the roots may contribute.This peak reinforcement is com-monly expressed in the form of an additional cohesion component c r in the Mohr-Coulomb framework.
It is however difficult to predict the maximum root reinforcement due to the complicated interaction between roots and soil.In addition, root properties such as orientations, architecture, strength or stiffness can vary widely.In practice, the maximum root reinforcement is often linked to the product of root tensile strength (t r,u ), root area ratio (R ra , the fraction of soil cross-sectional area occupied by roots) and a multiplication constant k accounting for the root orientations, as suggested by Waldron (1977) and Wu et al. (1979): where i indicate different root diameter classes.The discrepancy between this model (WWM model) and field measurements is generally accounted for by a empirical multiplication factor k , often attributed to sequential mobilisation of roots and quantified using a fibre bundle model, ranging from k = 0.3-1.0(e.g.Mao et al., 2012;Pollen and Simon, 2005;Bischetti et al., 2009) Although this approach has its practical merits for everyday slope stability calculations, it does not advance our understanding of the underlying dynamics of reinforcement.The focus on the ultimate limit state furthermore overlooks that significant deformations might be required to mobilise root strength.For structures with tight deformation tolerances, for example rooted railway embankments or cuttings (Briggs et al., 2016), these deformations might exceed the serviceability limit state.To investigate the contribution of roots in such cases, greater insight is required into the mobilisation of rooted soil strength.Appropriate constitutive models, linking stresses and strains in rooted soils, are required if numerical simulations are to be conducted for assessing performance at the serviceability limit state.
Several constitutive modelling approaches can be distinguished.In the first approach, rooted soil is modelled using a single constitutive relationship.Świta la and Wu (2018) (validated by Świta la et al., 2018), incorporated the reinforcing effect of roots in a Modified Cam-Clay framework through an adaptation of the hardening law.Two additional empirical parameters are required that must be calibrated by means of element testing.Due to the small number of parameters this model might be useful for practical problems when satisfactorily calibrated, but is less suitable for providing further insight into the underlying reinforcement mobilisation mechanism.
Alternatively, constitutive models for soil reinforced with fibres might provide inspiration.Michalowski and Zhao (1996) proposed a model for fibre-reinforced sand, assuming uniformly distributed, rigid-perfectly plastic fibres that can either slip through the surrounding soil or break.Using a homogenisation technique based on the amount of dissipated energy, the increase in the yield criterion could be expressed in terms of fibre properties.Comparison with experimental tests yielded satisfactory results Michalowski and Cermák (2003).However, while assuming rigid fibres might be a reasonable assumption for the polyamide fibres used, a material with a typical Young's modulus of 5-12 GPa (Jia and Kagan, 2001), roots typically have much lower Young's moduli (in the order of magnitude of 100 MPa, e.g.see Meijer et al., 2018b,a).Therefore the gradual mobilisation of root strength will not be captured well by this model.
Alternative fibre-reinforced soil modelling treated the rooted soil as a composite material, where soil and fibres each behave according to their own constitutive behaviour.The stress in the composite is often calculated using the volumetrically weighted sum of stresses in the components of the composite (upper-bound form of the rule of mixtures, e.g.Sawicki, 1983;Di Prisco and Nova, 1993;Diambra et al., 2013).Most fibre-reinforced soil models assume that fibres only reinforce through tensile action and not through compression, shear or bending (Di Prisco and Nova, 1993;Michalowski and Zhao, 1996;Zornberg, 2002;Diambra et al., 2013), an assumption that is valid for thin, fibrous roots but not for thick, woody roots.
Muir Wood et al. (2016) suggested it might be possible to model root-reinforcement through similar techniques as developed for fibre-reinforced soil and provided a conceptual approach.Such composite modelling framework requires the following components: 1.A constitutive model for the soil constituent.This can be chosen based on specific soil conditions required, e.g.accounting for partial saturation, hardening/softening behaviour as required.
2. A model accounting for the change in pore water pressures resulting from root water uptake.A detailed model was recently presented by Woodman et al. (2020).This model component will not be further addressed in the current study.
3. A constitutive model describing the mobilisation of the mechanical strength of roots.
After introducing this framework in more detail, the current study primarily focusses on developing the third (mechanical root) part.Compared to models for artificial fibres, specific elements that should be addressed for plant roots are: Reorientation of roots during soil deformation: In the model by Diambra et al. (2013), fibres add reinforcement aligned with their non-displaced orientation, and the stress in the fibre is calculated using the macroscopic strain vector in the non-displaced fibre direction.This will be problematic when modelling a simple shear element containing mostly sub-vertical roots.In this case the model will predict hardly any reinforcement due to the unfavourable initial orientation.In reality however, soil strains will change the orientation of roots, therefore affecting their reinforcement contributions.For example, in the case of a simple shear element, roots that originally are loaded in compression (and therefore add no contribution) might eventually mobilise in tension when shear strains are sufficiently large.Since real roots require large tensile strains before peak stress is reached (typically 15-20%), significant reorientation will occur.
Root length: Fibre-reinforcement studies have often used relative short fibres, resulting in relatively straightforward homogenisation of stresses in the fibre to the composite stress.In contrast, roots can be up to several meters long, hence the length scale of soil deformation might be much shorter than the length or the roots, necessitating a different homogenisation procedure.
Incorporation of root breakage.Fibre-reinforcement studies typically have used short (≤ 50 mm), strong fibres, resulting in the tendency to continually slip through the soil rather than break.In contrast, roots are much better anchored in the soil due to their length and branching and have lower strengths, and experimental shear testing often shows many roots breaking.
Natural variation in root biomechanical properties: Roots are a natural biological material with highly variable material characteristics that vary within and between sections of root.
The underlying modelling philosophy is to develop a constitutive model that is based on individually measurable root, root-soil interface and soil characteristics.Such a model will have stronger predictive capabilities compared to existing, more empirical models and will provide a better insight into what is happening within each part of the rooted soil in terms of stresses and strains.
The mechanical root reinforcement part of the model will be compared to direct shear box tests conducted by Liang et al. (2017).In addition, the resulting shear displacementreinforcement traces will be compared to those predicted using the large deformation Winkler beam-spring model recently developed for roots by Meijer et al. (2019a,b).

Direct shear testing on rooted soil
Throughout this paper, the solid mechanics sign convention was used for stresses (tension positive) and strains (extension positive).Volumetric strain v is taken positive for expanding soil.Isotropic pressure invariants in the soil (p ), water (p w ) and air (p a ) are positive in compression.

Test details
Liang et al. ( 2017) conducted direct shear tests on 150 mm diameter PVC cores filled with (rooted) soil to quantify mechanical root reinforcements.In all cores, soil was compacted to an initial dry bulk density of ρ d,i = 1.4 Mgm −3 to a height of 500 mm.Four sets of tubes were prepared.One set (2 replicate cores) was planted with Salix viminalis (willow, variety Tora) and tested after 2 months of growth.Another set (3 cores) was planted with Lolium perenne × Festuca pratensis (grass hybrid) and sheared after 2 months of growth.A third set (2 cores) was planted with Ulex europaeus L. (gorse) and tested after three months.A final, unrooted ('fallow') set of 3 cores were prepared to act as a control.All plants were grown under controlled climatic conditions.
Plant shoots were removed prior to shear testing.This removed any effects of plant transpiration and associated hydrological reinforcements on the test results, to focus on mechanical reinforcements only.Cores were subsequently fully saturated and then drained by lowering the water table to 500 mm depth for 48 hours to achieve pore pressure equilibrium.
All cores were sheared at z sh = 100, 200, 300 and 400 mm depths at a rate of 1 mm min −1 to a final shear displacement of u sh = 100 mm.No additional overburden pressure was applied.
After shearing, roots crossing each shear plane were manually counted per 0.1 mm diameter root class (Figure 1).The root area ratio (R ra , fraction of the area of a plane occupied by roots) decreased with depth in all cores.Gorse and willow cores contained thicker roots (root diameter d r > 1 mm) near the surface while for tests at larger depths and all grass tests the roots were mainly thin (d r ≤ 1 mm).
The soil behaviour was fitted using the Modified Cam Clay (MCC) soil model.Because triaxial tests on overconsolidated samples did not show distinct peaks in deviatoric stress, a Hvorslev surface was introduced on the dry side of critical with a gradient M H = M .The full MCC yield surface was used as the plastic potential function, resulting in non-associated flow on the dry side of critical.The best fit for triaxial shearing, isotropic compression and oedometric compression resulted in MCC parameters Γ = 2.08, λ = 0.115, κ = 0.015, M = 1.35 and ν = 0.4 (Figure 2).This M value corresponds with a critical state friction angle of φ cv ≈ 33.4 • .Apart from overestimating the amount of dilation in overconsolidated triaxial tests, the fit was satisfactory.
Samples tested by Liang et al. (2017) were conducted on samples with very low levels of matric suction (1 to 4 kPa, depending on the depth of shearing), and were therefore incorporated into the soil constitutive model using Bishop's effective stress framework.

Initial stress conditions
The exact value of vertical stress acting on each shear plane in tests by Liang et al. (2017) was difficult to predict due to 1) the soil being partially saturated, and 2) unknown side friction along the core walls.Using this critical state friction angle measured previously, the effective stress level acting on each plane in each core was back-calculated from the measured shear resistance in fallow core shear tests, resulting in an effective stress profile with depth z: Figure 2: Element test results (solid lines) and predictions using the Modified Cam Clay with the added Hvorslev surface (dashed lines).p 0 is the effective isotropic pressure at the start of the shearing phase.
with σ s,zz the vertical effective stress in the soil.The increase in effective stress with depth due to the self-weight of the soil was counteracted by a reduction in effective stress due to a reduction in suction with depth and friction along the sides of the tube.Assuming the compaction procedure by Liang et al. (2017) could be modelled as oedometric compression, the MCC (with Hvorslev surface) soil model was used to estimate the initial horizontal effective stress in the soil σ s,xx in the soil and therefore the coefficient of lateral earth pressure at rest (K 0 ), resulting in K 0 = 0.44.

Root biomechanical properties
The biomechanical behaviour of the three root species was measured in uniaxial tension tests.Peak tensile strength and Young's moduli were previously reported by Liang et al. (2017).However, because the current study is concerned with the mobilisation of root strength, the full tensile stress-strain behaviour needs to be characterised.These are described in terms of engineering stresses (t r ) and strains ( r ), as is common practice in root biomechanics research.
Root tensile strengths (t r,u ) may vary as a function of root diameter.This is commonly fitted using a power-law relation (e.g.Mao et al., 2012): where d r is the root diameter and d r,0 = 1 mm (a reference diameter).For the tensile strength of gorse and willow roots, no clear diameter trends were detected (Figure 3), so β t = 0 was assumed.Similar curves were used to fit the root tensile strain to peak r,u : Because no diameter trends were observed, for all species β = 0 (Figure 3) was assumed.
Both the tensile strength and strain to peak showed considerable scatter.This was captured by fitting probability distribution functions to normalised tensile strengths ( tr,u = t r,u /t r,u,f it ) and normalised strains to peak (ˆ r,u = r,u / r,u,f it ).Uniform distributions (characterised by a minimum and maximum value) could be fitted satisfactorily, see Figure 3. Therefore, the probability of root occurrence p r , given normalised strength and strain to peak: Despite the large variation in measured strength and strain to peak, tensile stress-strain curves all had similar shapes (Figure 4).Roots were modelled as linear elastic to simplify the development of the constitutive model.Tensile stiffness E r was defined as the secant stiffness at peak strength, rather than the Young's modulus, as it was considered more important to approximate the root tensile stress accurately at relative large deformations:    3 and 4).The blue shaded area indicates the domain of the best uniform distribution fits of normalised strength and strain to peak.The number of successful test is 39, 147 and 52 for gorse, grass and willow, respectively.

Gorse
Grass Willow 0.00 0.25 0.50 0.75 0.00 0.25 0.50 0.75 0.00 0.25 0.50 0.75 0.00 0.25 0.50 0.75 1.00 In summary, the mechanical behaviour of a bundle of roots can be characterised using tensile strength-diameter parameters (t r,u,0 and β t ), tensile strain to peak-diameter fit parameters ( r,u,0 and β ) and probability densities of the normalised strength/strain.All of these parameters have a clear physical meaning and can be simply obtained from a series of uniaxial tensile tests.

Constitutive model development
Rooted soil might require large deformations to reach peak capacity because of the low stiffness in the root.Deformations of the composite material are therefore described using the deformation gradient tensor F , which enables integration into finite-deformation framework.J indicates the Jacobian determinant (J = det F ), a measure for the volumetric expansion/contraction of the material (J ≈ 1 + v )

Phase relationships
The rooted soil was treated as a four-phase material, consisting of soil grains (volume fraction φ g ), water (φ w ), air (φ a ) and roots (φ r ), see Figure 5.The volume balance satisfies: For the purpose of defining phase relationships, both solid phases (soil grains and roots) are considered incompressible, therefore: where φ r,i and φ g,i are the root and soil grain volume fractions in the initial state, respectively.Assuming the mass of roots is negligible compared to the mass of soil particles: where ρ d,i is the initial dry bulk density of the soil material and ρ g the density of the soil grains.
The four-phase material was modelled as a composite consisting of two constituents: (a) 'soil' (subscript 's'), i.e. the mixture of air, water and soil grains, and (b) roots (subscript 'r').The presence of roots makes some void space inaccessible to the volumetric behaviour of the soil through two mechanisms: (a) roots take up void space, and (b) roots may prevent nearby void space from being used by the soil ('stolen voids', as hypothesised by Muir Wood et al., 2016).The volume fraction of voids belonging to the root phase φ vr ('stolen voids') is expressed as a fraction ξ r of the root volume fraction: Thus it follows that only the voids that are 'available' to the soil (φ vs ) are considered compressible.The void ratio of the soil (e s ) is therefore defined as the fraction of the 'available' pore space over the volume of soil grains (Figure 5): The presence of roots reduces the specific volume v s (v s = 1 + e s ) of the soil, and therefore makes the soil behave 'more densely'.
'Stolen' voids: Total stress in soil constituent Soil void ratio: Total stress in composite:

Soil-root composite stresses
Each of the two constituents ('soil' and 'root') in the composite material followed its own constitutive law (Figure 5).The rule of mixtures (e.g.Sawicki, 1983) was used to calculate the total (Cauchy) stress σ in the composite material: Where σ s and σ r are the (total) stresses in the soil and root constituent respectively, and φ s = 1 − φ r the volume fraction of the soil constituent.
The major advantage of the composite modelling approach is its flexibility.For example, the soil constitutive model can be interchanged to suit specific soil conditions, without having to alter or recalibrate the model used to calculate the root reinforcement effect.

Soil stresses
Any description for the total stress in the soil can be used, for example Bishop's effective stress formulation: where σ s is the effective stress in the soil, I indicates the second order unit tensor, p w and p a are the pressures in the water and air phases respectively (positive in compression), and χ is the Bishop parameter, often assumed to be related to the saturation ratio S w .(Unsaturated) pore fluid flow models can be used to describe the evolution of air and water pressures p a and p w , see Woodman et al. (2020).Subsequently in this paper, the effective stress in the soil is calculated using the Modified Cam-Clay type model described in Section 2.2, which was calibrated appropriately to measured soil behaviour for the test conditions considered.The MCC model inherently assumes parity between increments in specific volume and increments in volumetric strain multiplied by the initial specific volume ( v = v i ˙ v ), which becomes invalid when using the updated definition for specific volume accounting for root space (Equation 11).Small changes are therefore required in the definition of (a) soil elastic bulk stiffness K, linking changes in elastic volumetric strain e v and isotropic effective stress p , and (b) the volumetric hardening behaviour, linking plastic volumetric strains ( p v ) and pre-consolidation pressure (p c ): where bracketed terms are additions to the original model formulations.This demonstrates that the presence of root volume and stolen voids increase the soil bulk stiffness.Conveniently, this effect can be incorporated using an unmodified version of the MCC model that uses the standard (unrooted) definition of specific volume (v = 1/φ g ) using modified stiffness parameters κ * and λ * : This also necessitates a change in the critical state line intercept in v-ln p space: where v i is the initial specific volume of the soil according to standard (unrooted) definition.

Root stresses
Key differences between previous modelling of fibre-reinforced soil and modelling of rooted soil that were addressed: 1. Roots are typically much more flexible than polypropylene fibres modelled previously (E r ≈ 10 2 MPa instead of 10 4 MPa), necessitating modelling of the mobilisation of their strength.Because large deformations might be required to mobilise the full potential of the roots, finite deformations and root reorientation need to be considered; 2. Roots are typically much longer than polypropylene fibres, calling for a different homogenisation procedure of root strength into continuum behaviour at a single integration point; 3. Roots may break as well as slip through the soil when not sufficiently anchored; 4. Root biomechanical properties typically show large variability.
It was assumed that roots behave like long fibres, and therefore only reinforce in tension (i.e.zero bending stiffness).Roots were assumed to be cylindrical, straight and unconnected to other roots, with diameter d r , cross-sectional area A r and length L r .Roots were modelled as linear elastic (E r = t r,u / r,u ) with tensile strength t r,u .

Single root: elastic
Consider a single root whose initial orientation is described by unit vector n r and its initial volume fraction by φ r,i .The stretch in the composite material in the direction of the (deformed) root (Λ s ): where s is the (engineering) strain in the soil in the direction of the displaced root.Since roots are only assumed to reinforce in tension, their reinforcing effect is only taken into account when s > 0. The unit vector m r describing the orientation of the deformed root: Previous stress homogenisation techniques for fibres relied on the assumption that the deformations in the soil around the fibre are homogeneous along the entire fibre length.This is unlikely to be the case for the roots, which can be up to several meters long.Here, it is therefore assumed that soil deformations are only present along part of the root length, henceforth called the soil-loaded root length L s = ζ r L r .The coefficient ζ r (0 < ζ r ≤ 1) accounts for the fact that the typical length scale of the deformation in the composite (e.g. a shear band within a slope) might be (much) smaller than the length of the roots.Near the middle of the root, soil strain in the direction of the root is assumed to be s , and outside Length scale of soil displacement L s this region it is assumed zero.In the case of very short roots or fibres, or when modelling an element test, it follows that ζ r ≈ 1.
The matter at hand is now to find the root stresses in the middle of the root.Roots mobilise tensile stresses by mobilising soil-root interface shear resistance τ i .The interaction between soil and root in the rhizosphere is complicated and may vary depending on soil and root conditions.τ i is difficult to quantify, for example by an axial pulling test, as roots may be tortuous and branched.It was therefore incorporated in the model by means of a generic (and therefore flexible) cohesive-frictional approach.The interface friction τ i (defined per unit area of undeformed root) is assumed to fully mobilise as soon as there is differential soil-root displacement (rigid-perfectly plastic behaviour): where a i is the interface adhesion, δ i the interface friction angle and σ s,n the the average normal effective soil stress acting on the deformed root: Mobilisation of root-soil interface friction changes the tensile stress and strains along root axis s: Solving this differential equation using the assumed interface behaviour and zero strain boundary conditions at both root ends, the profiles of root tensile stress t r , root tensile strain r and relative soil-root displacement (u s − u r ) along the root can be determined (Figure 6).Depending on soil conditions and root properties, a root will respond in one of three ways: 'Anchored': The root is sufficiently well-anchored in the soil so that it closely follows soil deformations.No relative soil-root displacement occurs at either root end or near the middle of the root (Figure 6a); 'Partially slipping': The root slips within the soil near the middle of the roots (slip is defined as finite relative soil-root displacement), but not near either root end (Figure 6b); 'Fully slipping': Root-soil interface friction is mobilised along the entire root length, and tensile stresses can therefore not increase any further.Slipping occurs at both root ends (Figure 6c).
The maximum tensile stress in the root (t r,m ), occurring in the middle of the root, and the criteria to determine the type of response are (Figure 7b): Where t r,s is the tensile stress that would develop in the middle of the root in case soil-root interface friction is mobilised along the entire root length: Soil deformations will also change the number of roots present on a cross-section normal to the deformed root orientation (Figure 8).Taking this into account, the Cauchy root stress tensor for a single root, defined in the deformed state, can be written as: And after substituting Equation 8:

Single root: Root breakage and variation in root biomechanical properties
Roots will break when the tensile stress t r,m exceeds the tensile strength t r,u .Alternatively expressed, 'anchored' roots break when s ≥ r,u , 'partially slipping' roots break when ζ r s t r,s ≥ r,u t r,u , and 'fully slipping' roots break when t r,s ≥ t r,u ,.Once roots are broken, they are assumed to lose all reinforcement potential.Since roots cannot 'unbreak' during unloading of the soil, it will be necessary when modelling stress and/or strain cycles to keep track of the largest soil strain ( s,y = max s (t)) and root slippage stress (t r,s,y = max t r,s (t)) experienced by the root thus far.
A root breakage parameter f break is defined (Figure 7c) such that:   The average tensile stress in the root t r,m can by found by integrating the product of (a) probability of root occurrence p r , (b) elastic tensile stress in the root t r,m , and (c) the root breakage parameter f break , over the domain of root tensile strength and root tensile strain to failure (Figure 7): Although the double integral in Equation 29 appears complex, due to the simple (piecewise) nature of the expressions of p r , f break and t r,m , this integral and its derivatives can be expressed analytically, enabling rapid computation.An example of how t r,m might develop with increased soil strain s is presented in Figure 7e.At zero strain, all roots will behave as 'anchored', but with increasing levels of strain roots will start to break and/or slip.Even at large deformations, a 'residual' root reinforcement might be present due to a fraction of roots continuing to slip through the soil.
The (averaged) root stress tensor for a single root thus equals:

Multiple roots
The rule of mixtures (Equation 12) can be easily expanded to include multiple roots: Conceptually, it would be elegant to express the total root stress tensor as an integral over all orientation, volume fractions, root properties etc.This approach was for example followed by Diambra et al. (2010) in triaxial conditions, aided by the fact that in such conditions the orientations of principal stresses and strains are known.Here, due the the inclusion of large deformations and defining the root stress tensor in a generalised, three-dimensional stress and strain state, such an integral cannot be expressed analytically, and therefore has to be approximated by summation (Equation 31).Lebedev quadrature can be used to approximate the (three-dimensional) distribution of root orientations, and a finite number of root classes (each with their own length, diameter, biomechanical properties etc.) can be used to approximate the natural variation in root properties.
4 Numerical simulation of direct shear -Input

Fibre-based model (this paper)
The developed constitutive model was tested against direct shear tests performed by Liang et al. (2017).Several assumptions were required.All roots were considered vertically oriented initially.This is justified by the fact they were grown in the cores, promoting vertical growth.Every root diameter class (see Figure 1) and their corresponding root volume were included in the model.Lacking exact measurements of root length for each root, every root was assumed to have grown all the way to the bottom of the core (root length L r = 500 mm).Although the model allows for inclusion of 'stolen' voids, these were not included in simulations (ξ r = 0) Figure 9: Schematising the behaviour in direct shear using a drained, simple shear element.
as suitable values for ξ r have not yet been determined by means of systematic experimental studies.The root-soil interface was assumed perfectly rough (δ i = φ cv ) with no adhesion (a i = 0).The statistical variation in both root peak tensile strength (t r,u ) and strain to peak strength ( r,u ) was taken into account (data from Figure 3).The behaviour in the direct shear tests was approximated by modelling a drained, simple shear element at the shear plane, similar to Muir Wood et al. ( 2016), see Figure 9.The two displacement directions free to move are the vertical deformation (F zz ) and the simple shear displacement F zx (x-axis points in the direction of shearing, z axis points in the normal direction).Normal total stress σ zz and pore pressures were kept constant during shearing.
A key parameter required is the height of the simple shear element, which was assumed as h sh = 30 mm, based on the thickness of the shear zone as observed by Bull et al. (2019) for the same soil during direct shear tests conducted in a X-Ray CT scanner.Direct shear displacements (u sh ) were directly linked to simple shear deformations F zx (Figure 9): The soil-loaded root length L s (i.e. the length of root along which soil strains are assumed non-zero) was initially assumed equal to the shear zone thickness (L s = h sh ).The influence of this parameter is investigated later.The constitutive behaviour of the rooted soil was solved incrementally, solving increments in total stress in the composite σ as a function of increments in the deformation gradient tensor Ḟ .The rule of mixtures in incremental form: A small strain formulation of the soil model was used.To account for the significant rotations that might occur when shear strains were large, the Jaumann objective stress rate was used to calculate an increment in Cauchy soil effective stress σ s : Where D ep s is the incremental MCC elasto-plastic soil stiffness matrix, spin matrix W = L − L T /2, rate of deformation tensor D = L + L T /2, and velocity tensor L = Ḟ • F −1 .
The amount of root stress mobilised is a function of both the deformations in the composite material (F ) as well as the amount of confinement pressure the soil applies on the root (since the root-soil interface friction depends on the effective stresses in the soil).Therefore, the incremental root stress tensor σr : where D F r and D σ s r are the incremental 'stiffness' matrices linking increments in the deformation gradient tensor ( Ḟ ) and effective stress in the soil ( σ s ) receptively to increments in the root stress tensor.

Root beam model
The predictions of the fibre model were compared to the 'beam model' recently developed by Meijer et al. (2019a), This model assumes that roots act as beam elements with a finite bending and axial stiffness, that can move independently from the soil, and allows for large beam deformations and rotations.The displacements of the soil are an input parameter in the model, and the same simple shear mechanism as used for the fibre model was assumed.
The same assumptions for root quantity, root diameters, root-soil interface friction and root biomechanical behaviour as used in the fibre model were assumed.One exception is that the average root tensile strength t r,u,f it and tensile strain to peak r,u,f it were used for every root, rather than statistical distributions, as the root beam model does require root stiffness to be given as a scalar value rather than a statistical distribution.
The model used for calculating transverse and axial soil-root interface resistances was the same as used by Meijer et al. (2019a).Equation 2was used as the required profile of vertical effective stress input, and the soil angle of internal friction required for calculation of both resistances was taken as φ cv .
Roots were assumed to have broken once their tensile capacity was exceeded.Once broken, their contribution to root reinforcement was set to zero.The reinforcement response of all individual roots was summed to acquire the total root reinforcement.The calculated reinforcements were added to the measured average fallow soil resistance to obtain predictions for the root-reinforced strength.
For more details about this model and its assumptions, see Meijer et al. (2019a).

Results
The fallow soil behaviour in direct shear could be fitted well with the MCC+Hvorslev surface soil model previously calibrated against triaxial and oedometer tests.Similar to element tests results, no clear peaks in shear strength were observed although there are many fine-scale sharp-peaks corresponding to small scale variation (Figure 10).The direct shear tests on rooted soil showed that significant shear displacements, often of the order of 50 mm, were required to reach peak root-reinforced shear strength (Figure 11).Some tests showed subsequent softening behaviour, associated with root breakage.
Both the fibre model and the beam model yielded generally good predictions for the stress-strain behaviour of the rooted soil (Figure 11).In most cases, both models predict (a) the peak strength of the rooted soil, (b) the gradual mobilisation of root reinforcement and (c) the presence of significant reinforcement even at the end of the shear test (shear displacement u sh > 90 mm) reasonably well.Note that the predicted traces are the same for each depth because of the assumed initial stress conditions.
Varying the length of root along which soil strains exceed zero (soil-loaded root length L s ) in the fibre model shows that reducing this size has little effect on the magnitude of peak root reinforcement, but 'delays' mobilisation of root reinforcement, resulting in larger displacements to peak and larger predicted root reinforcements at the end of the shear tests (Figure 12).This highlights the importance of assumptions made during the homogenisation of (long) roots into continuum behaviour.When choosing a length scale similar to the experimentally observed shear zone thickness (L s ≈ h sh ) good predictions were achieved, supporting the hypothesised link between soil-loaded root length and the typical length scale of the soil displacement mechanism.
Comparing the observed peak strength to model predictions shows that both the beam and fibre model yield good predictions (Figure 13).Reduction factors k required to match the WWM predictions (Equation 1) to the experimentally measured results varied between species and was 0.3 ≤ k ≤ 0.6.The largest deviations in predicted peaks strengths occurred for gorse roots.It was hypothesised that this is due to the presence of a few thick roots, for which little biomechanical test data was available.For example, in the gorse tests at 100 mm depth, two individual roots with diameters exceeding 2 mm account for more than half of the root volume (Figure 1).There is a realistic chance these roots may have been weaker than the assumed average, given the large observed variation in root tensile strength (coefficient of variation: 50%, Figure 3).This highlights the importance of sufficient replication, both in root biomechanical testing as well as shear testing.
The predicted stress-strain behaviour in a single test was investigated further to highlight the complicated interaction between the soil and root components, see Figure 14.Initially, strains are low and root reinforcement therefore negligible.The overconsolidated soil therefore dilates.With increased shear strains, roots start to mobilise their tensile strength.The shear resistance of the composite now increases due to two effects: (a) resistance in the roots, and (b) soil hardening, as the roots pull the top and bottom of the simple shear element closer   together, explaining the predicted soil contraction ( → ).At some point, roots begin to gradually break, resulting in a reduction of reinforcement.The soil is now free to dilate as it naturally tends to (softening: → ×).Because some roots will continue to slip, the root-reinforced shear resistance remains higher than the unrooted resistance even at large shear strains.This example showcases the complex interaction between soil and roots, but also the power of the composite modelling approach.

Discussion
Both experimental and model results showed both hardening behaviour (associated with root reinforcement mobilisation), softening behaviour (associated with root failure), as well as significant root reinforcement at large shear displacements near the end of the test, indicating that a significant large fraction of roots slipped through the soil rather than broke.This demonstrates the importance of incorporating both breakage and slippage mechanisms, in contrast to previous fibre-reinforced soil models that did not include breakage (e.g.Świta la and Wu, 2018;Diambra et al., 2013).
The increase in shear deformations required to reach peak shear stress in root-reinforcemed soil are in line with previous direct shear studies, e.g.75-84 mm in soil reinforced by willow and 32-35 in fallow soil (Mickovski et al., 2009), or field tests performed by Ekanayake et al. (1997) on Kanuka trees (22-52 / 6-20 mm) or Monterey pine (22-52 / 6-20 mm), and Docker and Hubble (2008) on rooted riverbanks (50-100 mm in rooted soil).This shows the importance of including large rotations/deformations in the constitutive model, and emphasises the importance of accounting for the serviceability limit state calculations when for example designing root-reinforced embankments with tight tolerances.Incorporating root reinforcement into finite element analyses by adapting the soil yield criterion (e.g. through an increase in soil cohesion) is unable to capture this 'delayed reinforcement' effect, and is also unable to account the observed soil softening due to root breakage.However, such approaches may be sufficient if only the ultimate limit state is of interest.The results also show the importance of assumptions made during homogenisation of (long) roots into a continuum response.The hypothesis that the assumed length of root loaded by soil (L s ) is related to the characteristic length scale of the deformations in the soil was confirmed by the direct shear tests.However, this length scale is not known in a typical boundary value problem.Previous work on direct shear testing on fibrous and rooted soil have indicated that fibres/roots may increase the thickness of the shear zone, but this effect has not been studied quantitatively (Jewell and Wroth, 1987;Shewbridge and Sitar, 1989;Abernethy and Rutherfurd, 2001;Fan and Su, 2008).There is a need for further studies investigating the failure mechanism in rooted soil without applying such rigid displacement boundary conditions as in the direct shear test.Field investigations of failures in rooted slopes, or physical model testing of rooted slope stability, for example using a centrifuge (Askarinejad and Springman, 2015;Liang et al., 2017), will be able to provide such data.
The fibre-based model predicts similar reinforcements compared to the beam model by Meijer et al. (2019a) for the shear experiments considered.While the latter accounts for potential bending and shear effects in roots, these were negligible for most roots in the shear experiments because of their relatively small diameters (Figure 1).
Comparison of the peak root reinforcement results to popular Wu/Waldron type models (WWM) confirms the necessity of an additional reduction factor k in these models, as previously found by many authors.Required values for the tests conducted were different for each species, and varied between k = 0.3 and k = 0.6, in line with previous studies (e.g.Mao et al., 2012;Pollen and Simon, 2005;Bischetti et al., 2009).
The fibre model shows that roots and soil interact in a complex fashion.Not only will soil displacement affect root reinforcement, but root reinforcement will in turn also affect the behaviour of the soil (Figure 14).This root-soil-root interaction may in part explain the current persisting difficulties in predicting root reinforcement when using a WWM-type framework, incorporating root properties only in reinforcement predictions (e.g.Wu et al., 1979;Pollen and Simon, 2005).
The developed root-reinforced constitutive modelling framework allows incorporating root reinforcement into finite element analyses.Its generic formulation allows for straightforward implementation of the coupling between mechanical and hydrological effects in the soil, for example additional suction pressures introduced by roots (Woodman et al., 2020).A major advantage of the adopted composite framework is that it can be used with any existing effective stress-based soil model.This is important for practical applications as soils may vary widely, and varying climatic conditions and ongoing plant transpiration may, in specific applications, result in larger matric suctions compared to the experiments by Liang et al. (2017).
The constitutive model for the root constituent requires a sizeable amount of input data (root quantities, root geometry, root biomechanical properties, interface properties), contrary to the approach developed by Świta la and Wu (2018).However, all parameter have a clear, physical meaning, and can be measured in the lab or predicted using existing data.
The proposed framework opens up a clear pathway to untangle some of the complexity of the mechanical behaviour of rooted soil.Further experiments to verify the model are required.
Triaxial tests on large samples (their size carefully chosen to allow for sufficient root lengths) would be extremely valuable as they would allow simultaneous investigation of stress and volumetric behaviour.Such testing will help to verify key assumptions (for example those surrounding 'stolen voids') and help to identify where model simplifications and assumptions can be made safely without losing accuracy, systematically working towards a model which is both accurate and more practical (requiring fewer input parameters) for rooted soils.

Conclusions
A constitutive framework was developed for root-reinforced soil, adopting a four-phase composite modelling approach using separate constitutive models for the root and soil phases; A formulation for the stress in the root phase was developed following concepts previously developed for fibre-reinforcement, accounting for the gradual mobilisation of root reinforcement, large rotations and deformations, and both root breakage and slippage effects; The model showed good comparison with experimentally measured direct shear displacementshear stress data measured for rooted soils.It matched both the magnitude and the gradual mobilisation of root reinforcement, as well as the significant 'residual' root reinforcement at large deformations associated with slipping roots, making it suitable for both serviceability and ultimate limit state analyses within finite element calculations; A previously developed beam model, modelling roots as spring-supported beams with both axial and bending stiffness, matched both experimental results and fibre model results well; Exploration of the model revealed a complex two-way interaction between soil and root behaviour.The model provides a suitable framework for future investigations into the key mechanisms governing the hydro-mechanical behaviour of rooted soil.

Figure 1 :
Figure 1: Cumulative root area ratio (R ra ) for each direct shear test.

Figure 3 :
Figure3: Root tensile strength (t r,u ) and tensile strain to peak ( r,u ) as a function of root diameter.Blue, dashed lines indicate the best power-law fit (Equations 3 and 4).The blue shaded area indicates the domain of the best uniform distribution fits of normalised strength and strain to peak.The number of successful test is 39, 147 and 52 for gorse, grass and willow, respectively.

Figure 4 :
Figure 4: Normalised root tensile stress-strain curves.Dashed lines indicate the assumed root linear-elastic behaviour.

Figure 5 :
Figure 5: Phase relationships and definition of soil-root composite behaviour

Figure 6 :
Figure 6: Profiles of tensile strain, tensile stress and relative soil-root displacement along a stretched root.

Figure 7 :
Figure 7: Calculating the average tensile stress in a root as function of soil strain in the direction of the root s , accounting for (a) statistical variation in root biomechanical behaviour, (b) different root response types and (c) root breakage.This results in (d) tensile stresses in roots at different strains, from which the (e) average tensile stress in roots can be expressed as a function of soil strain.

Figure 8 :
Figure 8: Soil element containing root (left) undergoes deformation F (right), resulting in modified root area ratio.
r,u,min p r f break t r,m dt r,u d r,u

Figure 10 :
Figure10: Experimentally measured fallow (unrooted) direct shear resistance (average of two replicates at each depth) and numerical soil model predictions for every shear plane depth.Note that the predicted traces are the same for each depth because of the assumed initial stress conditions.

Figure 11 :
Figure 11: Experimentally measured and predicted root-reinforced shear resistances.Solid, grey zones indicate the average shear resistance in non-rooted tests.

FallowFigure 12 :
Figure12: Effect of the length scale of soil loading along the root (L s ) on the shear resistance predictions by the constitutive model.Dark grey shaded area indicates the experimentally measured average fallow resistance, and the light grey shaded area the experimentally measured root-reinforces resistance (All experimental data used taken from z sh = 200 mm depth; Gorse replicate 2; grass replicate 3; willow replicate 1).

Figure 13 :Figure 14 :
Figure13: Comparison between experimentally measured root reinforcement and predictions using various models.The solid black line indicates parity.k is the gradient of the fit for each model.This is the multiplication factor required to match model results with experimentally measured results.
r,s Partially slipping roots: ζ r t r,s < E r s < Greek symbols:β : Root diameter-tensile strain to failure power law coefficient β t : Root diameter-tensile strength power law coefficient Γ: Critical state line intercept in v-ln p space in Modified Cam Clay model δ i : Interface friction angle r : Root tensile (engineering) strain r,u : Root tensile (engineering) strain to failure r,u,0 : Tensile strain to failure for root with diameter d r = d r,0 ˆ r : Root tensile (engineering) strain, normalised by tensile strain to failure s : Soil tensile (engineering) strain v : Volumetric strain ζ r : Ratio of length scale of soil deformation L s and root length L r κ: Modified Cam Clay model parameter Λ: Stretch parameter λ: Modified Cam Clay model parameter M : Modified Cam Clay model parameter M : Oedometer primary compression parameter C s : Oedometer unload-reload parameter c r : Root reinforcement d r : Root diameter d r,0 : Root reference diameter E r : Root tensile stiffness e: Void ratio F : Deformation gradient tensor f break : Fraction of root currently unbroken I: Unit tensor J: Jacobian determinant K: Bulk stiffness parameter K 0 : Coefficient of lateral earth pressure at rest k : Wu/Waldron multiplication factor k : Sequential root mobilisation reduction factor L r : Root length L s : Length scale of soil deformation m r : Root orientation unit vector in displaced state n r : Root orientation unit vector in initial state p : Isotropic pressure p c : Isotropic preconsolidation pressure p r : Probability of root occurance q: Deviatoric stress R ra : Root area ratio t r : Toot tensile stress tr : Root tensile stress, normalised by root tensile strength t r,s : Root tensile stress in slipping roots t r,m : Root tensile stress in the middle of a root t r,u : Root tensile strength t r,u,0 : Tensile strength of root with diameter d r = d r,0 u H : Hvorslev yield surface parameter ν: Poisson's ratio ξ r : Fraction of 'stolen' voids with respect to root volume fraction ρ d : Dry density σ: Cauchy stress tensor σ : Cauchy effective stress tensor τ i : Interface shear stress φ: Volume fraction φ cv : Soil critical state friction angle Latin symbols: a i : Interface adhesion C c sh : Direct shear displacement v: Specific volume z sh : Shear plane depth Commonly used subscripts: a: Air g: Soil grains i: Initial r: Root s: 'Soil', i.e. the mixture of soil grains, water and air v: Voids / volumetric y: Yield