Mathematical and computational modelling of vegetated soil incorporating hydraulically-driven finite strain deformation.

Abstract In this paper a new model for the hydro-mechanical behaviour of rooted soils is developed. It is a physically-based model that couples finite strain soil deformation with unsaturated water and air flow, while improving on existing cohesion-based approaches to mechanical root reinforcement and empirical soil water-uptake approaches typically used to deal with rooted slopes. The model is used to show that the dynamics of soil-water pressure and soil deformation depend strongly on the physics of the root-water uptake and the elasto-plastic soil mechanics. Root water uptake can cause suctions and corresponding soil shrinkage sufficiently large to necessitate a finite-strain approach. Although this deformation can change the intrinsic permeability, hydraulic conductivity remains dominated by the water content. The model incorporates simultaneous air-flow, but this is shown to be unimportant for soil-water dynamics under the conditions assumed in example simulations. The mechanical action of roots is incorporated via a root stress tensor and a simulation is used to show how root tension is mobilised within a swelling soil. The developed model may be used to simulate both laboratory experiments and full-scale vegetated slopes.


Symbols
seasonal shrinkage and swelling of soils) and increase the stability slopes through both the generation of pore water suctions and mechanical reinforcing action of the roots.
Field and laboratory experiments have demonstrated the mechanical reinforcement effects of plant roots on slopes [11]; [12], with roots able to increase significantly the apparent shear resistance of the host soil [13]; [14]; [15].Relatively thick woody roots are able to withstand bending moments and can act as 'soil nails' [9], while thinner more flexible roots can provide tensional reinforcement.In taking up water from the soil, plant roots can develop pore water suctions up to 1.5 MPa, increasing the effective direct soil stresses substantially [11].Pore water pressure changes transmit beyond the rooting zone, effectively extending the influence of the roots to greater depths and thereby influencing more deepseated failure mechanisms [16]; [4].In [17] it was demonstrated that the trees growing on earthworks can maintain pore water suctions beneficial to slope stability, particularly during the winter months in a temperate climate, when transpiration is low and rainfall high.In [18] and [19] researchers measured vertical ground surface displacements of several centimetres in clay railway embankments, as a result of seasonal changes in plant transpiration and following tree felling.
Rigorous simulation of the soil-root composite requires the mechanical behaviour of the soil and roots, and the fluid (air and water) flow to be integrated and accounted for.There is considerable demand for models that can replicate triaxial, shear-box and centrifuge tests of vegetated soil [20], and account for 1 We want to point out that different communities use different terminology to refer to soil that has both air and water in the pore space.Some use the term nonsaturated soil and others use a more accurate term of partly saturated soil.We will use the term partly saturated soil throughout this paper.
the influence of vegetation on typical engineered slopes such as railway and highway embankments [21].
A number of numerical models of vegetation effects on slopes have been developed or implemented in commercial software, but these often have limitations.For example, some approaches do not consider mechanical deformation and/or failure [22]; [17]; [23], or do not have mechanical-hydraulic coupling [24]; while other models do not account for partial saturation of the soil [25], or simplify hydraulic boundary conditions for example by ignoring rooting depth [23].
A useful model of rooted soil may need to: (i) incorporate mechanical ground movements and failure, and be able to deal with relatively large (finite) strains; (ii) couple this to partly saturated water and air flow; and (iii) represent key mechanical and hydrological aspects of the plant roots and overall plant physiology relevant for plant transpiration.This paper describes the development of a coupled physically-based soil deformation model that accounts for the influence of plant roots and is parameterised by independently measurable quantities.The aim of this paper is to show the development of this model, explaining how the root mechanical stress-strain behaviour can be incorporated, and examining the key assumptions required.Particular consideration is given to examining whether a finite-strain approach is necessary (as opposed to a typical infinitesimal engineering strain approximation); the choice of constitutive model for estimating soil deformation under root-water uptake; whether the fluids need to be treated as a two-phase, i.e. air and water, flow system; and whether deformation could be sufficient to modify the intrinsic permeability hence affect soil water flow.The individual constitutive relationships (including the soil behavioural model) can be readily changed to suit a particular scenario.Standard and relatively simple components have been adopted, so the focus is on the dynamics of coupling rather than a critique of the components, each of which is already widely accepted individually.
The basis for the modelling assumptions is set out below, with reference to the literature, and is then appraised though the investigation of simple soil-water interaction scenarios.

Consideration of mechanical and fluid behaviour of rooted soil
The behaviour of the soil matrix is simulated using a Modified Cam Clay model [26]; [27] with non-linear elasticity, hardening/softening, state dependency and a critical state.This well-established model is chosen because it is computationally efficient and requires only a limited number of input parameters, compared with other models simulating similar aspects of behaviour.
There is an increasing body of evidence for the nature of the mechanical behaviour of roots in soil (e.g. [9]; [28]; [15]).A common method for representing the action of the roots is to adjust the "apparent cohesion" parameter on a Mohr-Coulomb failure envelope [29]; [14]; [1]; [11]; [30]; [2]; [4].This increase in "apparent cohesion" is generally estimated using the volume fraction (root area ratio) of roots and their tensile strength [31].Fibre-bundle models are often used because in an element of soil containing roots of different diameters, roots may fail progressively rather than simultaneously [32]; [33].To account for non-negligible bending stiffness (e.g.large diameter tree roots), [34] used a beam-on-a-nonlinear-Winkler-foundation model to estimate an appropriate equivalent apparent cohesion.A different approach was followed by [3],[ [35]], who simulated the hardening effect of roots empirically by means of strain-driven expansion of the pre-consolidation pressure within a Modified Cam Clay model.In essence, these approaches try to capture the behaviour of a composite (soil and root) into a single material model.The approach adopted here is different, in that the soil and roots are treated as separate phases, each with their own constitutive behaviour and both contributing to the stress state within the composite material.Thus, the parameters and history of each material may be physically characterised individually.This approach has been adopted in modelling soil reinforced by short fibres [36]; [37].In such models, the development of root tension during strain results in an increase in the effective stress within soil around the roots, thereby increasing the shear resistance of the soil.
The root contribution is determined from knowledge of the geometry, orientation and biomechanical properties of in-situ root bundles via a fibre-based root-soil interaction approach [38].This paper describes how this phase is incorporated within the model, and illustrates it by means of a simple example.
It is assumed that air and water movement in the soil can be described by two-phase Darcy-Richards' equation.Root water uptake (RWU) may be calculated on the basis of the geometry of the root system (e.g.[39]; [40]; [41]; [21]; [42]; [43]; [44]; [45]; [46], which can increasingly be identified using 3D X-ray Computed Tomography (CT) scanners [47].However, for slope-scale application it is likely that root water uptake will need to be considered macroscopically at the level of the whole planted soil [48,49] as considering each individual root on this scale is computationally prohibitively expensive and it is unlikely each and individual root matters; rather it is their collective hydro-mechanical behaviour that matters.
RWU is frequently used as a continuum sink term with Richards' equation [50].It is commonly taken to be proportional to the suction produced by the roots and distributed over the depth profile according to the density of roots (as parameterised by the length of roots per unit volume of soil [2]; [51]).These semi-empirical approaches typically adopt an inhibition term with increasing water suction in the soil (e.g.[52]; [3]; [1]; [53]; [54]).Despite its practical application, this approach has been criticised as unrepresentative of the actual process and it can be difficult to measure the parameters directly [55].In addition, most of these models do not account for the overall plant transpiration stream from soil to roots to stems and leaves.
To account for the local pore pressure distribution around each root, researchers [49] used semianalytical solutions (a method subsequently adopted by [56]).In [57] researchers used a localised Richards equation for the near-root flow coupled to a macroscopic water flow model.This is appropriate where roots are sparse, but [58] demonstrated that the hydraulic equilibration times are likely to be considerably shorter than other key timescales of interest such as typical growth time-scales, and environmental changes over weeks or longer.Thus, it is considered that the pressure between root subbranches will relatively quickly reach a quasi-equilibrium value provided the root density is not too low.
Radial water flow to the root has previously been examined with reference to the detailed anatomy of the root (e.g.[59]; [60]; [61]; [62]).Water must pass through the cellular structure of the root cortex to reach the xylem (water-carrying capillaries).The cortex is a complex structure, but flow has been categorised into (i) apoplasmic, pressure-driven flow in the extracellular space, and (ii) symplasmic, an osmotic flow passing though cells and cell walls.In this study the bulk macroscopic water movement is taken to be dominated by the apoplasmic pathway, neglecting osmotically-driven flow [57].However, it is relatively easy to incorporate this in our model, but in addition to the water uptake term in the Richards' equation needing to be modified, we would also need to introduce an extra equation to describe the chemical changes in the soil and within the plant.We consider this beyond the scope of this paper and refer an interested reader to the literature that deals with this [63].The focus of this paper is to integrate the mechanical and hydrological effects plant roots have on soil.
Vertical (axial) flow inside the root xylem can be described by Poiseuille's law [64] if the distribution of the radii of functioning vessels can be estimated.However for thin roots, the axial conductivity may be considerably greater than the radial conductivity, so some studies have assumed a single constant pressure in the root system [21]; [41]; [48].[58] showed that for maize, this may be valid for zeroth order roots (main axes), but not for first order roots (lateral branches), which could be simulated by means of a distribution coefficient on the zeroth order root, resulting in neither the radial nor the axial conductance dominating the pressure distribution in the root.Therefore, in this study, the xylem pressure is allowed to vary axially along the root.
Both the mechanical and hydraulic approaches to root behaviour will relate self-consistently to the geometric parameters of a root bundle (i.e. the distribution of diameters and lengths).

Conceptual model
Based on the considerations outlined in Section 2, an isothermal 4-phase porous system is conceptualised (Figure 1), comprising three overlapping continua: water (w), a homogeneous granular/porous matrix (s) and air (a).The porous component is assumed to have a single set of homogeneous mechanical properties.The water phase has overlapping continua for flow in the soil and flow within the xylem in the root, interacting via root water uptake (  ).The soil volume fraction comprises the solid particles, air and water phases (i.e.  +  +  ), but not the roots.Mechanical coupling is via the effective stress (′), root stress (  ) and pressure terms (  ,   for water and air respectively).A capillary pressure difference (  ) exists between the water and air phases.
The roots are envisaged as a further overlapping continuum, affecting both hydraulic and mechanical behaviours.Hydraulic interaction between the roots and the soil is via the source/sink term in the soil water equation.Mechanical interaction is via a root stress tensor, which modifies the soil effective stresses.The solid phase is coupled to the hydraulics via the effective stress.The volume taken up by plant roots is here considered as a homogenised part of the solid (i.e. for simplicity a soil-root composite is assumed).It would be straightforward in future studies to further partition the volume to include the root fraction explicitly.
The assumption of a bulk mechanical response relies on the averaging scale significantly exceeding the scale of inhomogeneity produced by each root.Thus the approach would not be valid for individual roots having a significant individual effect, for example a single long, thin root transmitting tensile force over a distance or a small number of large-diameter roots capable of carrying significant bending stress.The homogeneity assumption also neglects macropores, formed for example by cracking or biological activity.Infrastructure slopes in the UK are often fairly uniformly vegetated; however, caution should be exercised, particularly for individual large plants or trees, which may possibly need to be modelled explicitly (e.g.[65]), and in cases where root systems are very sparse.
The compressibility of water and of the solid particles is assumed to be negligible.Owing to the high diffusivity of air, gas compressibility is only important under special circumstances (for example compression of individual entrapped air bubbles beneath a saturated wetting front); hence, following [66], it is neglected in this paper.
The difference in pressure between the two fluid phases is described by a single capillary pressure curve.
For simplicity in introducing the model, wetting and drying hysteresis effects, and coupling between void ratio and moisture retention [67], are neglected.As data and necessity dictate, such constitutive relationships could be easily be incorporated if needed.The possibility of mixed phases or isolated (static, non-percolating) phases is also neglected.
A finite (i.e.large or non-infinitesimal) strain approach is adopted, since it does not require the assumption of infinitesimal engineering strain.Thus the model is not restricted to small strain applications, allowing it to be used in conditions where significant deformation occurs e.g.shear box tests and extensively plastically-deformed slopes.Furthermore, roots typically require large strains to mobilise their full capacity (15-20% typically).Therefore, to investigate the behaviour of rooted soil at peak root capacity, soil behaviour may have to be cast in a large-deformation framework to remain accurate.
The concept of combined modelling presented in this paper is novel in that it integrates a large-strain mechanical model of soil coupled with partly saturated air and water flow and physically-based root water uptake equations.This moves away from more empirical approaches for root water uptake (e.g. [52]), for which the parameters are not as easily obtainable by direct measurement.

Model development
The kinematic behaviour of a solid continuum at large strains requires rigorous treatment [68]; [69]; [70]).Within the Eulerian (spatial) reference frame, a point in space is described by vector  and within the Lagrangian (material) reference frame, a point in the material by .Displacement of the material in relation to the reference frame is  =  − .The displacement vectors can be decomposed into = ̂+ ̂+  � and = ̂+ ̂+  � , where ̂ and ̂ are unit vectors on the horizontal plane and  � is the upward vertical unit normal.
The governing equations are cast in Eulerian co-ordinates, in which physical parameters are readily interpreted.However, the domain itself changes size and shape and the boundaries move, which needs careful numerical treatment.For numerical implementation, the model is therefore converted to Lagrangian co-ordinates, where nominal parameters referenced to the initial reference state and boundaries are fixed (as detailed in Supplementary Material A).

Eulerian equations
As shown in Figure 1, the total volume is divided into four distinct phases of air (a), root (r), soil (s) and water (w).The volume fractions, , of each phase sum to unity, i.e.: The solid volume fraction   (, ) is the sum of the soil and root components, i.e.   =   +   and is related to the initial solid fraction,  ,0 =   (, 0) via the Jacobian determinant  = det(), i.e.: where the deformation-gradient tensor, , defined in relation to the Eulerian frame [70], is: The solid velocity can also be derived in relation to : The mass (volume) balance of the water phase is given by: where   (, ) is a water-sink term representing root water uptake (RWU).The mass conservation equation for the air phase is: where gaseous exchange between the roots and the soil is neglected.Equation [5] enables Darcy's law to be implemented for water flow relative to the soil, i.e.   (  −   ) = −(  ⁄ )�∇  +  � �, within the water balance: where   is soil pore water pressure,  is water density,  is gravitational acceleration,  is the Darcy permeability of the soil to water and  is the water (dynamic) viscosity.Similarly, manipulating Equation [6] enables implementation of Darcy's law for air flow relative to the solid, i.e.   (  −   ) = −(  ⁄ )�∇  +    � �, within the air balance: where   is the air pressure,   the Darcy permeability of the soil to air,   the air density and   the air (dynamic) viscosity.Additionally, a van Genuchten [71] water-retention relationship is assumed, in which the water saturation, defined as , is a function of the air and water pressures: where  is an exponent, ,  is a soil-water retention parameter.More generally, the moisture retention relationship observed in experiments has a separate wetting / drying path (i.e.exhibits hysteresis).It is also altered by changes to the distribution of pore sizes and shapes, which may be brought about by both roots occupying pore space, organic matter content changes and deformations to the solid skeleton [67].In principle, equation [9] can be modified to reflect these effects, where sufficient data are available for the soil; however this is neglected here to retain the simplest possible treatment.
In partly saturated conditions, a proportion of the pore space becomes unavailable for water flow.This is addressed by introducing an effective permeability for water flow, k w , which, following [71], depends on the relative saturation ratio   as: , for   >   , else  =   , [10] where  is the intrinsic permeability, and  and  are constants.Similarly, the effective permeability for air, k a , is given by [72]: The model includes a further level of coupling, whereby the initial intrinsic permeability  0 may vary with the void space geometry [73] according to a Kozeny-Carman type relationship [70]; where k 0 is the intrinsic permeability at a reference stress of 4 kPa, the total (fluid, i.e. air and water filled) porosity is   =   +   and the subscripts 0 denote initial conditions.Refuting earlier assertions that this relationship does not apply to very fine-grained materials, [74] show that, once artefacts and experimental error are eliminated, the equation predicts the saturated hydraulic conductivity of most soils fairly well.
A simple upscaled approach for water pressure in an equivalent single root is adopted.This is a firstorder RWU relationship for steady state radial pressure-driven flow [75]: where   is the root diameter,   is the root radial conductivity (m/Pa/s) (a bulk measure of the crosssectional geometry and effective diffusivity over the ensemble of roots) and   is the root internal xylem pressure.This is a quasi-steady state simplification of unsteady pressure transmission through the anatomy of the root [59].[58] showed that, since the cortex dimension scales with the size of the root,     may be considered to be independent of root diameter and therefore potentially applicable to a range of root sizes.In general   is a function of the distance down the root, but it is held constant for the purposes of this paper, aggregating all depth dependence in the root transfer term into   , which is the length of roots in each unit volume of soil.Treating any changes to internal water storage as negligible, the root pressure along a non-branching vertical root is given by: where   is the root axial conductivity (m 4 /Pa/s).Assuming   to be constant with depth gives: Since xylem capillary flow may be reasonably estimated using Poiseuille flow, the bulk rate is given as the summation of all individual xylem vessel contributions [58]: where   is the radius of the xylem vessel,   is the number of open xylem vessels with cross-sectional radius   and  is the index of each radius category.  may therefore be either measured [64] or estimated based on the xylem radius distribution using equation [16].This single cylindrical equivalent root model can be used to aggregate the RWU from higher order lateral roots by assuming the flow rate from higher order roots to be smoothly distributed over the zeroth order roots [58].The upper bound for RWU occurs when the soil water is at atmospheric pressure and the root base suction pressure is a maximum.The plant root to stem to leaves transpiration stream is captured as a boundary condition to equation [15] as it links leaf water potential to the atmospheric pressure and ultimately appears as a boundary condition at the base of the root system for equation [15].
The final component of the model is the mechanical behaviour of the solid skeleton.Neglecting inertial effects, the equation for force equilibrium in Eulerian co-ordinates is: where compressional stress is defined as positive on the inward normal, the (Cauchy) stress tensor is denoted by  and the gravity body force per unit volume,   = −(    +     +     +     ) � .
Expressing the force balance in terms of the effective stress,  ′ , gives: where  is the identity matrix.Coupling of the solid mechanics with the fluid pressures thus arises via the effective stress.Here, he effective soil stress is represented by the so-called Bishop stress,   ′ = ( −   ) + (  −   ), with the parameter  =    [76] where  is an exponent that varies as a function of plasticity index, after [77].It is known that   ′ does not represent a true effective stress for soils below a certain critical initial degree of saturation (which may be as high as 90 % for clays), as it does not account for the volumetric compression or collapse that can occur on wetting.However, it is shown by [78] that the Bishop stress expression with χ replaced by the saturation ratio S w does have physical significance.In any case, this simplification can be revised if data show  to be a non-linear function of   at high and low water contents [79]; [80]).  is the (Cauchy) stress tensor within the root volume fraction.Roots can develop tension as they are extended, and thicker roots may hold shear and bending stress.A key facet of the model is the assumption that the constitutive behaviour of the soil is unaffected by the presence of the roots, other than through modification of the effective stress tensor.In essence, the soil is treated as a fallow material, with the roots influencing its behaviour through modification of the effective stress state.
The constitutive relationship between effective stress and strain for the soil can be computed in terms of elastic behaviour as  ′ = :   , where  is the (fourth-order) elastic tensor.An appropriate elastic strain measure is the Euler-Almansi strain, defined by: where the elastic deformation gradient (Eq. 3) is given by   =   − and   is the plastic deformation gradient tensor.Rearranging Equation [18] gives: For simplicity, residual saturations of air and water are assumed to be zero.
The Modified Cam-Clay (MCC) model [26]; [27] is used to populate the elasticity tensor , with stressdependent stiffnesses.The elastic response (within the yield surface) is expressed succinctly using a decomposition into volumetric and deviatoric components [27], as: where    is the elastic volumetric strain,    is the elastic deviatoric strain,  ′ = 1 3 ( ′ ) is the mean effective stress,  is the deviator (von Mises) stress and  is the slope of the elastic (unload-reload) line in the :ln′ plane.There is a choice as to whether to fix the shear modulus  or the Poisson's ratio  [27];  is fixed here, hence the shear modulus is calculated as: . [22] The elliptical yield surface is given by: where  is the slope of the critical state line in the : ′ plane.It is assumed that the soil obeys the normality condition, i.e. that the vector for plastic strain increments is in the direction of the outward normal to the yield envelope.The plastic response is given as: where  =  ′ ⁄ and  and  are the slopes (in the :ln′ plane) of the normal compression line and unload-reload lines, respectively.
The constitutive relationship,   , between the root stress and root strain can be formulated based on upscaled root data (including the distribution of root orientation, thickness and breaking strain, consideration of soil-root friction and inclusion of root failure mechanisms).A detailed methodology for this process is given by Meijer et al. (2020), in relation to new mechanical test data.
For the purposes here of demonstrating the mechanical coupling under 1D loading conditions, a very simple root model is assumed as follows: there is a rigid root-soil interface (i.e. the root strains are equal to the soil strains,   = ); the roots provide only a vertical tensional force; root stress occurs only in tension; the roots act after time   , otherwise hold no tension; the cross-sectional area of the roots is unchanged by stretching (i.e.Poisson's ratio for the roots is zero); there is a linear relationship between tensional stress and strain, cast in an intermediate reference frame that defines the initial state of the roots; see Supplementary Information A; and the roots do not reach their breaking strain.
The complete model thus comprises 6 independent equations for 6 independent variables (,   ,   ,  ,   ,   ), which are solved together with initial and boundary conditions on the external surface, .The equations are summarised in Table S1 in Supplementary Material A. Initial and boundary conditions for a particular scenario are given in Section 5.1.The proportions of each phase are calculated as The model is implemented as a system of partial differential equations (PDEs) using COMSOL 5.3a (2016) [81], via the inbuilt 'general form PDE'. The COMSOL Structural Mechanics and Geomechanics modules were used to solve the equations of the mechanics of the deforming solid and for the COMSOL implementation of Modified Cam Clay.
Simulations to verify the approach are given in Supplementary Material B.

5.1
Model set-up

Boundary and Initial Conditions
The model is used to examine finite elasto-plastic deformation of a vegetated soil during a dryingwetting cycle.The soil is first consolidated under self-weight (for 10 days), dried for 120 days by root water uptake (RWU) and downward drainage, and then re-wetted by continuous rainfall.Three variants of this scenario are simulated.In Simulation 1 the air pressure is kept constant at atmospheric pressure and the intrinsic permeability is assumed to be constant.Simulation 2 is identical to Simulation 1, except that the intrinsic permeability of the soil is varied to account for void contraction.Simulation 3 is also identical to Simulation 1, except that the air pressure is computed.In Simulations 1 and 2, the air phase pressure -which should be governed by equation [8] -is held constant at   =   ; a typical simplification applied to partly saturated soils.In Simulation 4, the root stress provides a tensional stress proportional to the extension experienced during the re-wetting phase.
The soil is bounded by vertical boundaries allowing frictionless vertical movement, but preventing lateral movement (analogous to an idealised oedometer or a one-dimensional (1D) soil column).Thus the soil solid undergoes 1D consolidation and subsequent swelling, driven by its weight and changes in effective stress due to changing pore water pressures.The consolidation under self-weight is a modelling necessity and is not intended to represent a real geological process, although the ramping of the self-weight load has similarity with "spin-up" in a geotechnical centrifuge.
The soil profile is initially 5 m deep with the upper surface flat.The upper boundary is defined as the free surface,  ∈ (, ), and moves as the solid deforms, i.e. as  =   + , where the initial surface   is a horizontal plane at elevation  =  above the base.On the upper boundary, the column has a water flux boundary condition for the infiltration of rainfall, , which is set to 0 for the first 130 days and to 5 mm/d thereafter.For air, the upper boundary is given by the fixed atmospheric pressure condition,   =   .The lower boundary is restrained from moving in all degrees of freedom.At the base ( = 0) the water pressure is fixed to   =   , ensuring that the lower boundary is always saturated with water, making it impermeable to air and allowing drainage from the base.The root base pressure (i.e., the pressure within the root at its base, at the soil surface)  = -1 MPa during root water uptake (from 10 to 130 days) and is otherwise 0 MPa.Root base pressure  integrates all above ground plant transpiration effects into one single parameter accounting for the atmospheric temperature, pressure, humidity, windspeed etc. outside of the plant and the plant above ground stem and leaf stomatal conductance all of which can and do vary with time owing to the interaction between vegetation and the weather [82], [83].But for illustration purposes we assume it to be constant during the day and zero during the night.
There is a no-flow condition at the tip of the root (initially located at a distance   below the base of the root), as owing to root physiology roots do not take up or bleed water through their tips.The initial condition for the whole solid is a displacement of zero,  = 0.The initial condition for the pore water pressure is for the soil to be just less than fully saturated, i.e.,  =   − 1 kPa.The initial xylem pressure   is in equilibrium with the initial local pore water pressure, i.e.,   = .Prior to root water uptake, the soil is consolidated under self-weight, with the body force increased linearly over 10 days starting from zero.

Parameters
The parameters used in the simulations, based broadly on a clayey silty sand known as Bullionfield soil [12], are given in Table 1.As in conventional soil mechanics practice, it is assumed that the soil grains are incompressible.A wide variety of values for the root uptake parameters   and   have previously been identified, depending on species [85]; [64]; [21].  = 1×10 -13 m 4 /s/Pa is adopted here as central to the range reported in the literature.  is set to 1×10 -14 m/s/Pa, which gives a maximum root water uptake of 4.8 mm/d when   = 1 MPa and   = 0 MPa (this is based on observations at Newbury, UK, by [86] and gives a July monthly average potential evapotranspiration of between 2 and 3 mm/d).For the period of zero RWU,   is set to 0 m/s/Pa for t>130 days.
In Simulations 1-3, all terms in the root stress tensor   and the volume fraction allocated to the roots   are zero.However, in Simulation 4, the root stress tensor is computed based on the deformation that takes place after 130 days (i.e. the start of the wetting period).The root is assumed to be anchored from this point in time, with zero slip between the root and the soil.This is primarily a modelling necessity, although it is likely that by this time the root has been able to fix itself to the soil through the growth of short lateral branches.

Model results
Simulation 1: fixed pore air pressure   and intrinsic permeability   Figure 2b shows relatively modest reductions in pore water pressures after ten days of consolidation and downward drainage.RWU then generates stronger negative pore water pressures (plotted every 30 days after the initial consolidation period).By 130 days from the start of the simulation, the pore pressure in the soil adjacent to the roots (i.e. in the top 1 m of soil) has approached the xylem pressure of -1 MPa.Thereafter, pore water pressures recover as a wetting front moves down from the surface in response to rainfall.
The degree of saturation (Figure 2c), mapped from the pore water pressure via the soil water characteristic curve or water retention relationship (Equation [9]), follows a similar pattern.At 250 days the soil is drier than the initial state as a result of the vertical downward drainage through the base of the column.The pressures are tending towards a steady-state value of ~4.7 kPa, whereat the hydraulic conductivity is equal to the 5 mm/d rainfall rate.The void ratio reduces at a gradient (with increase in the natural logarithm of the mean effective stress) of =0.015,steepening to =0.1 along the isotropic normal consolidation line during plastic deformation.
On re-wetting the soil unloads elastically.At 164 days, the unloading stress path crosses the x-axis, giving zero deviatoric stress.Owing to the continued increase in pore pressure due to rainfall, the vertical effective stress then rises above the lateral effective stress (thereby increasing the deviator stress back above zero).Figure 3c shows that by the end of the simulation, the soil has contracted relative to the initial condition (i.e., the re-wetting has not returned the soil to the starting point, but to a void ratio of 0.915 compared with 1.0 initially).
The right-hand column in Figure 3 shows the behaviour deeper in the soil, which is dominated by selfweight.The soil at 0.5 m elevation encounters the yield surface after 7 days and the consolidation pressure ′  increases to 71 kPa.Thereafter, the soil undergoes a small elastic swelling.This expansion is due to the reduction in total stress at this depth, owing to the removal of water (hence reduction in soil unit weight) from the soil above by RWU. ′ decreases as the roots remove water, then recovers after 130 days in response to the continuous rainfall.Thus RWU causes the soil to follow opposite loading directions at different depths, the effective stress being dominated by changes to the pore pressure at shallow depth (resulting in a large increase) and by changes in the total stress associated with a reduced unit weight of the overlying soil at greater depth (resulting in a small decrease).
For comparison, the responses given by a linear-elastic (LE) soil model are plotted in the top row of The largest (Euler-Almansi) volumetric strain encountered at the half-rooting depth (0.5 m below the surface) is -4.2%.The engineering strain (based on an assumption of infinitesimal strain) is -4.0%(i.e. there is a 5% relative error between the strains).Given that the root system base pressure is set to -1 MPa, yet plant roots are capable of inducing soil pressures down to -1.5 MPa [9], there is potential for even greater strains generated purely by root water uptake than observed in this simulation.Hence the errors introduced by assuming small strains may be significant for some applications.

Simulation 2: variable intrinsic saturated permeability k S
Given the potentially significant changes in effective stress due to RWU it may be possible that the deformation of the solid skeleton is sufficient to affect the permeability, hence the coupled fluid flow and mechanical response.This has been shown to affect consolidation in saturated media at large strains [87]; [70].However, Simulation 2 showed no significant differences in pore pressures compared with those assuming a constant intrinsic permeability.Adjacent to the root, the soil pressures approached -1 MPa after 120 days and the volumetric solid fraction (  ) increased from 0.5 to 0.52.From Equation [12], the corresponding intrinsic permeability is reduced to 80% of the uncompressed value.The reduction in  due to desaturation calculated by Equation [10] is eight orders of magnitudes greater (to k = 8×10 -21 m 2 ); it is this that dominates water flow, hence pore pressures and displacements.Thus given the assumptions made and for the conditions simulated, the degree of saturation (i.e.Equation [10]) dominates the hydraulic conductivity -any changes arising from changes in porosity are very much a second order effect.It is possible, however, that if a sufficient surcharge is applied to a soil close to saturation, changes in  due to changes in porosity may become significant.This might be important, for example during spin-up in a geotechnical centrifuge (e.g.Liang et al., 2017).
Simulation 3: dynamic air phase (variable pore air pressure   ) The air-phase pressure is not kept constant in this simulation, but is calculated using Equation 8.The lower boundary is impermeable to air (i.e., a zero-flux condition applies) and the upper boundary is set to atmospheric pressure.The initial air pressure condition is a positive aerostatic distribution, i.e.   =   ( − ), giving a positive air pressure of 0.06 kPa at the base of the model.Downward drainage of water by gravity followed by RWU starting at day 10 causes air to move into the soil, under a relative suction, to replace some of the volume of water drained.Conversely, when rainfall is wetting the soil, the internal air pressure rises above the aerostatic line and air is discharged from the surface of the soil.The soil conductivity for air is higher than for water, and the pressure gradient for a given volumetric flow rate is proportionately lower.The air pressure at the base of the model reaches a maximum at the end of the simulation of 0.21 kPa.This is negligible in comparison with the changes in pore water pressure, which were up to 1 MPa.For simulations under such circumstances, it may be expedient (to save on computational time) to fix the air pressure, as is frequently done in standard unsaturated zone water flow simulations.

Conclusions
Given that water-soil-root interactions are potentially very significant, detailed simulation of vegetated slopes will require a coupled approach that includes physically realistic root water uptake dynamics.To achieve this, a new physically-based macroscopic model for vegetated soils has been introduced, characterised by well-established parameters that are relatively straightforward to measure.On the basis of the simulations presented, the following conclusions may be drawn.
1. Root water uptake may lower the pore water pressure sufficiently for the associated change in effective stress to cause finite strain in the soil.This implies that a linear-elastic model or an assumption of infinitesimal strain should be adopted with caution when modelling the mechanical behaviour of the rooted zone.
2. For the clayey silty sand modelled, permeability is dominated by the degree of water saturation, not by deformation (reduction in void ratio).
3. For the conditions simulated, air pressure changes make little difference to the estimated pore water pressure changes and soil deformations caused by root water uptake and rewetting by rainfall.It may be expedient (in terms of computational efficiency) to assume that the air phase remains at constant pressure.
4. For situations where there is soil expansion or shear, the stress carried by the roots may be included via a root stress tensor (which will be a function of stress and strain) to provide a rigorous basis for making improved predictions of the behaviour of vegetated slopes.
Amit Patil from COMSOL is thanked for his help in remedying discrepancies we found in the Modified Cam-Clay code (the corrections will appear in version 5.5).

Figure 1 :
Figure 1: Schematic of multi-continuum rooted soil system.The volume is divided into four phases (i.e.  ,   ,   ,   ).The water phase has overlapping continua for flow in the soil and flow within the xylem in the root, interacting via root water uptake (  ).The soil volume fraction comprises the solid particles, air and water phases (i.e.  +  +  ), but not the roots.Mechanical coupling is via the effective stress (′), root stress (  ) and pressure terms (  ,   for water and air respectively).A capillary pressure difference (  ) exists between the water and air phases.

Figure 2a shows the
Figure2ashows the period of soil drying by RWU, followed by rain wetting.The rate of drainage from

Figure 2 -
Figure 2 -simulation 1. 10 days of settlement, 120 days of drying due to RWU (solid lines), followed by 120 days of rewetting due to 5 mm rainfall (dashed lines) (a) Flows out of the column by downwards drainage, RWU and Rainfall into the column (negative outflow) (b) pore-pressure profiles with depth, (c) the corresponding degree of saturation and (d) vertical displacement.Legend gives times (in days) from the start of simulation (shade darkens with time, dashed line for wetting period, solid lines for drying period and dot-dash lines for first 10 days).

Figure
Figure2dshows how ~100 mm of vertical surface displacement accumulates progressively over the 5 m

Figure 3 .
Figure 3.The LE bulk modulus is matched at the initial condition, i.e.  = ′   ⁄ .The MCC model

Elevation 4 .Figure 3 :
Figure3: Simulation 1. 10 days of settlement due to self-weight, followed by 120 days of drying due to RWU, then 120 days of rewetting due to 5 mm rainfall.Left hand column is for 0.5 m below the surface (i.e. at 4.5 m elevation, midway down the root zone): (a) Stress path for the soil (p' is effective stress and  is deviator stress; defined in the nomenclature) (c) void ratio, where solid line is the Modified Cam-Clay model (MCC) state path and (e) effective stress p' and pore water pressure p w against time.Right hand column is for 0.5 m elevation (i.e.4.5 m from the surface): (b) stress path for the soil, (d) void ratio and (e) effective stress p' and pore water pressure p w against time.LE denotes Linear-Elastic, CSL is the Critical State Line.

Figure 4 :Simulation 4 :
Figure 4: Pore air pressure for Simulation 3. 10 days of settlement, 120 days (d) of drying due to RWU (solid lines), followed by 120 days of rewetting due to 5 mm rainfall (dashed lines).Simulation 4: mobilised root reinforcement Figure5contrasts the vertical displacement of the soil surface where the mechanical reinforcement effect of the roots is neglected (simulation 1), with simulation 4, in which it is included.The root tension resists re-expansion of the wetting soil, giving rise to a smaller upward movement of the soil surface over the final 120 days of the simulation.

Figure 5 :
Figure 5: Vertical displacement of the soil surface for Simulation 4 where root strength is mobilised.The unmobilised simulation has otherwise identical parameters, but with zero root stiffness (i.e.  =0).