Preliminary analysis on the impacts of the rhizosphere on occurrence of rainfall-induced shallow landslides

Rainfall-induced diffuse shallow landslides are one of the most critical natural hazards as they often evolve into highly destructive flow slides and debris flows. Vegetation is recognised to play a key role in landslide occurrence and is frequently invoked as a potential remedial measure for slope stabilisation at the catchment scale. The beneficial action of vegetation is generally associated with mechanical (root anchoring) and hydrological (suction generated by root water uptake) effects. There is indeed a third effect that has thus far been little explored. The rhizosphere, the portion of soil directly affected by plant roots, is characterised by hydraulic conductivity higher than the underlying soil horizons. This significantly affects hillslope hydrology by promoting lateral diversion of rainwater. This paper presents a case study in Scotland where the rhizosphere is demonstrated to play a major role in controlling shallow landslides. Field investigation and laboratory testing were carried out to characterise the hydraulic conductivity of the rhizosphere and deeper horizons. In turn, this formed the basis for the development of a physically based model for the slope. The model was first validated against its capability to simulate failure of two historical landslides and then exploited to demonstrate the beneficial effect of the rhizosphere. The lesson learned from this study is that shallow landslide hazard can be mitigated by enhancing the capacity of the rhizosphere to act as a natural lateral drainage. This implies that plants with root-system architecture that enhances lateral subsurface flow should be privileged when designing vegetation-based remedial measures.


Introduction
Rainfall-induced landslides represent a natural phenomenon typically observed in mountainous areas after intensive or long rainy periods (Chen and Lee 2003). When heavy and/or prolonged rainfall events occur, water infiltrates into the slope and porewater pressures increase. In turn, these reduce the soil shear strength eventually triggering slope instability (Gonzalez-Ollauri and Mickovski 2017a; Sidle and Bogaard 2016). Their danger is mainly associated to the high velocities and long travel distance developed, causing property and infrastructure damage, injury and death.
It is widely acknowledged that the presence of vegetation can influence landslide occurrence. Plants have been proven to be an effective mitigation measure (Gonzalez-Ollauri and Mickovski 2017a) as they can significantly enhance slope stability thanks to a series of mechanical and hydrological effects (Liu et al. 2016).
A significant amount of literature showed how the roots can reinforce and stabilise the hillslopes via anchoring into deeper soil layers (Kim et al. 2017;Ghestem et al. 2011). Indeed, root biomechanical properties and their distribution in the root zone highly influence the contribution of the vegetation to the mechanical reinforcement on the soil (Ni et al. 2018).
Many studies have been carried out also to quantify the roots' mechanical properties of different plant species (Bordoni et al. 2016;McGuire et al. 2016;Capilleri et al. 2016). Moreover, there is evidence suggesting that the soil shear strength can be improved via the particle bonding associated with root exudates and microbiological activity (Naveed et al. 2017;Ali and Osman 2008;Fan and Su 2008;Zhang et al. 2010).
On the other hand, vegetation enhances the slope stability also through hydrological effects (Gonzalez-Ollauri and Mickovski 2017). It, indeed, promotes soil water extraction via the transpiration process occurring in the soil-plant-atmosphere continuum (soil water is driven from the soil through the plant to the atmosphere). This mechanism contributes to keep the soil in an unsaturated state thus enhancing its shear strength (Liu et al. 2016;Gerten et al. 2004).
There is indeed a third effect that has thus far been little explored. The root system can promote subsurface lateral flow in the rhizosphere by creating networks of self-organised macropores in which preferential flow may be triggered (Shao et al. 2017). Plant roots can, indeed, increase the roughness of the soil and the infiltration capacity of the ground (Wu 2019). It has been tentatively suggested that this hydrological mechanism may help 'preserve' low pore-water pressures in potentially unstable layers thus reducing the susceptibility of slopes to landsliding (Ghestem et al. 2011). However, assessing quantitatively the actual impact of the hydraulic 'diversion' promoted by the rhizosphere on landslide occurrence is still an argument of open discussion.
This paper presents a case study of shallow landslides where the rhizosphere is demonstrated to play a major role in controlling slope instability. Field investigation and laboratory testing were carried out to characterise the hydraulic conductivity of the rhizosphere and the deeper horizons. This characterisation formed the basis for the development of a physically based model for shallow landslides, which was first validated against its capability of simulating two historical landslide events and then exploited to demonstrate the beneficial effects of the rhizosphere.

Background
Plants can dramatically change the soil environment giving rise to the so-called rhizosphere effect (Jones et al. 2009). The term Brhizosphere^was first coined by the German agronomist Lorentz Hiltner 1904 to describe the area around the plant root inhabited by a unique population of microorganisms and influenced by chemicals released from plant roots. In the years since, the definition of the rhizosphere has been revised and refined. Many authors seem to agree with the idea that the rhizosphere is a region surrounding the roots characterised by a gradient in chemical, biological and physical properties (McNear 2013). Within the Landslides Original Paper rhizosphere, roots actively search for soil-based resources and nutrients flux between organic and inorganic pools (York et al. 2016). For the complexity of the mechanisms involved in it, the rhizosphere has been a crucial point in many studies for years.
The spatial extension of the rhizosphere is influenced by the root system architecture, which may vary, depending on soil and plant species and on their response to changing climatic and biological condition (Schlüter et al. 2018). During their life, roots change in diameter, length, direction, moving through the soil, covering large distances and actively modifying their surroundings (Baluska et al. 2018). Roots growth and decay may lead to changes in both soil hydraulic and mechanical properties (Lehmann and Or 2012;Shao et al. 2017). Thanks to a phenomenon called rhizodeposition (Jones et al. 2004), roots constantly release organic carbon in the surrounding soil, which has been found to produce the most dramatic changes in the physical, biological and chemical nature of the soil (Jones et al. 2009). Craft et al. (Craft and Broom 2002) showed, in example, how vegetated soils show very low value of dry density compared to bare soils (mean value 1.3 kN/m 3 ). Avnimelech et al. 2001 studied the correlation between dry density, porosity and organic matter in vegetated soils. In particular, the dry density was found to decrease with the increasing organic matter while the porosity showed higher values compared to the bare soils.
Root-soil interactions indeed alter the pore structure (Roose et al. 2016;Scholl et al. 2014). In particular, they contribute to the formation of a network of interconnected macropores (Shao et al. 2017;Ghestem et al. 2011;Sidle et al. 2001) with diameters ranging from a few tenths of micrometres (Marshall 1959) to many centimetres (Pierson 1983). Aubertin (1971) studied extensively the nature of the macropores in forested soils and observed their capability to alter subsurface flow, providing pathways or conduits for the rapid movement of free water into and through the soil profile. In particular, the experimental campaign carried out revealed a close correlation between the presence of old root channels and overall conductivity of the soil. Thanks to this hydrological mechanism, roots are capable of influencing water pressures in soils and have therefore a major control on the occurrence of shallow landslides.
There are few but significant examples in the literature about the interplay between root architecture, preferential subsurface flow and landslide occurrence. Gaiser (1952) has pioneered the concept that root channels may serve as large openings for rapid water flow. He mapped the presence of channels formed from decayed roots through the upper horizons of a temperate forest soil and concluded that the effect of root-induced macroporosity in the soil profile is an aspect that should not be neglected in the natural slope stability analysis. Noguchi et al. (1997) also highlighted the similarity between the direction of subsurface flow and the direction of root growth in the upper layers of a semitropical forest. Leung et al. (2015) explored the effect of the change in pore size distribution and consequently, water retention of soil induced by roots. Shao et al. (2017) investigated the impact of plant rootinduced preferential flow on hydromechanical processes of vegetated soils under different planting densities.
Because of their specific properties, root channels need to be considered with reference to their influence on preferential flow. Rhizosphere is therefore an important factor in addressing the complex nature of preferential flow at the hillslope scale, which has been assumed to influence landslide initiation (Sidle et al. 2001).
Case study: rest and be thankful (Scotland)

Site geology
The slopes of the A83 Rest and Be Thankful are located approximately 34 miles North West of Glasgow (Fig. . 1). The study section of the A83 between Ardgartan and The Rest and Be Thankful has a long and well-documented history of slope instability and landslide occurrence. The A83 forms part of the important trunk road network in Scotland and is vital for quick and easy access between the Western region of the country and the Central Belt and Highlands.
The site is underlain by a bedrock layer of steeply sloping metamorphic rocks from the Beinn Bheula Schist Formation (BGS 2012). The formation is composed of pelites, semipelites and a mixture of coarse to fine psammites. The rocks were formed roughly 542 million years ago in the Neoproterozoic where a period of regional metamorphism altered the original sedimentary rocks and caused large-scale folds and faults which can be seen throughout the present-day landscape (BGS 2012).
The superficial deposits on the slopes and the valley floor include mainly glacial till of quaternary age. The valley floor itself is composed of river terrace deposits including gravels silts and clays of quaternary age. The bedrock plays little part in the landslide activity on these slopes. The recent landslides here have largely been associated with slope deposits, including peat and topsoil as well as the underlying layers of colluvium. The colluvium comprises sandy to gravelly silts and clays with varying amounts of cobbles and boulders. The colluvium deposits on this slope represent earlier phases of slope instability (BGS 2012). The slopes of the A83 Rest and Be Thankful are highly populated by local vegetation growing, namely grass, bracken and ferns with almost no tree cover.

Field surveys
A number of field surveys were conducted between 2015 and 2017. Soil profile was examined by inspection of a number of landslide scarps still visible at the site. Figure . 2a shows one of these scarps where three soil horizons can be identified: topsoil, transition soil and subsoil. Overall, the layers overlying the bedrock form a cover that is approximately 1 m.
The soil horizons were also inspected using an open-end sampler as shown in Fig. . 2b. A soil segment was isolated in the tube sampler and trimmed to make it level with the rim of the open-end sampler. The measurement of the length of the segment together with the geometry of the cross-section of the open-end sampler made it possible to estimate the volume of the sub-sample (Online Resource 1). The sub-sample was then cleared away from the sampler, stored in a self-sealed plastic bag and its wet and dry mass measured once back to the laboratory. Following this procedure, in situ soil basic properties (i.e. degree of saturation, dry density, porosity) could be easily determined.
It is interesting to notice that the penetration depth does not coincide with the length of the soil column in the sampler. When forcing the sampler into the ground, the top soil reduced in volume significantly due to its very open structure and, hence, high compressibility.
The top soil is associated with the presence of the rhizosphere, the narrow region of soil that is directly influenced by root secretions and activity of soil microorganisms. Unfortunately, no direct measurement of the rhizosphere could be made at this stage. Detailed measurements of the rhizosphere require more sophisticated equipment (i.e. diffuse reflectance infrared spectroscopy or light microscope capable of detecting the amount of mucilage or exudates defining the rhizosphere). Unfortunately, these were unavailable for the investigations carried out during this work. However, for the purpose of this study, such detail of information was not required as input of the numerical model. The rhizosphere effect on the hydraulic and mechanical properties of the surrounding soil horizon was instead the main core of the numerical modelling. The topsoil extends for approximately 25-30 cm as shown in Online Resource 2 where the depth of the rooting zone can be clearly detected. In terms of grain size distribution, the three horizons are very similar as shown in Fig. . 3, suggesting that the three horizons have the same geological origin (glacial till). The biological activities taking place in the rhizosphere have just promoted soil aggregation giving rise to more open (lower density) texture.
A distinctive feature of the slope morphology is the presence of gullies cutting the planar slope parallel to the dip direction. These gullies are spaced between 30 and 50 m. The gullies are characterised by a relatively thin soil cover and appears to be always saturated. Online Resource 3 shows the saturated nature of the soil cover in the gully as observed after several days without rain.
The landslide case studies Rest and Be Thankful has been a source of problem for mobility in Scotland for years. The A83, which runs along the perimeter of the Rest and Be Thankful slopes, has been closed on a several number of occasions due to rainfall-induced landslides, causing significant disruption to local traffic. Only in the past 4 years, the A83 has been closed five times due to the instability of the slopes and the falling of landslides material on the road.
The two events analysed in the present work are the landslides events of 1 December 2011 and 1 August 2012 ( Fig. . 4). Those events have been chosen because the rainfall data available covered only the period going from January 2011 to January 2013 (Fig. . 5). Both these events happened after a period of heavy rainfall and they have a return period falling between 10 and 15 years.
The December 2011 landslide occurred after 65.8 mm of rain fell in 48 h, a quarter of the December's expected average for the area. The road was subsequently closed in both directions resulting in a 26-mile diversion. The road was kept closed for several days to complete clean-up operation safely.
The landslide is a translational slide that has degraded into a flow as it has passed over a local break of slope, moving approximately 100 t of material. In the lower section of the slope, a small gully was exploited. Figure . 4b shows the debris flow of August 2012 blocking the road. The event was reported to exhibit a translational failure surface, with up to 1000 t of debris falling onto the road. This event happened following a period of heavy rainfall. Again, the road was kept closed for several days to complete clean-up operations safely. Engineering and operational measures are now in place to minimise the impacts of subsequent events. In both cases, the failure surface was found to be placed at the soil-bedrock interface at a depth falling between 0.85 and 1 m.

Laboratory investigation of hydraulic conductivity of soil horizons
The saturated hydraulic conductivity of the three soil horizons was measured in the laboratory on samples collected in the field. Two different procedures were used to measure the hydraulic conductivity, depending on the texture of the soils and therefore the type of the sampler adopted for the specific soil horizon: oedometer cutting ring for the top soil and sampling tube for transition soil and subsoil.
The topsoil is very compressible and any sampler/cutting ring should have a very low height to diameter ratio to minimise the disturbance of the shear stresses developing at the inner surface of the sampler/ring during insertion into the soil. A large oedometer cutting ring cannot be used in the field as the soil sample is difficult to withdraw. Therefore, the sample was cut and trimmed from a block sample transported to the laboratory. The saturated hydraulic conductivity of the topsoil was measured on a specimen 80-mm-diameter and 20mm-high cut in the laboratory from a block taken from the field. The specimen was placed in an oedometer cell and subjected to 3-kPa vertical stress, with the aim to reproduce a grade of compaction close to the one in the field (3 kPa is equivalent to an overburden pressure generated by about 30cm soil column). A head differential of 14 cm was imposed also to reproduce low hydraulic head gradients similar to the ones occurring in the field. For example, this hydraulic head differential may reproduce vertical water flow generated by zero pore-water pressure at the ground surface (ponded water) and a suction of 1.4 kPa in the soil. Water flow was measured by monitoring the change in water mass of the upstream reservoir using a balance. At least three replicates have been investigated; the average value for the density of the samples was 11.57 kN/m 3 .
On the other hand, block samples could not be easily taken from the lower horizons (transition soil and subsoil) due to the cohesionless nature of these soils. Therefore, the saturated hydraulic conductivity of the transition soil and the subsoil were measured on samples still in the sampling tubes used to collect the soil from the field (samples having 3.8-cm diameter and about 10-cm height). The sampling tubes were placed on the pedestal of a triaxial cell and filled with water up to the top. This arrangement reproduces a rigid wall permeameter. The pedestal of the triaxial cell was then connected to a water reservoir to establish a constant hydraulic head differential of 10 cm. Water flow was measured by monitoring the change in water mass of the downstream reservoir using a balance. At least three replicates have been investigated for both the transition soil and subsoil with an average density of 12.7 kN/m 3 and 13.1 kN/m 3 .
The saturated hydraulic conductivity tests were carried out at low pore-water pressures (up to 1-1.5 kPa) to achieve conditions of 'field-saturation', i.e. the maximum level of saturation

Original Paper
Landslides practically achieved in the field due to pressure heads developed by ponding water. The field-saturated hydraulic conductivity measurements are summarised in Table 1.
As previously discussed in BBackground^, the presence of roots together with other biological activity in the rhizosphere (topsoil horizon) modifies the pore structure of the soil. In particular, the continuous growth of roots contributes to the development of macropores within the rhizosphere. The overall porosity of the topsoil is therefore affected by this phenomenon and this results in higher porosity of the top soil compared to the other soil horizons (Table 1). The porosity strongly influences the hydraulic conductivity of a soil; hence, the topsoil was expected to show a higher hydraulic conductivity. This was indeed confirmed by the test performed in the lab. This effect is much less evident in the deeper horizons, hence the hydraulic conductivity appeared to be lower.
Few tests were attempted in order to quantify the difference between the horizontal and vertical permeability. The vertical hydraulic conductivity was found to be half of the horizontal. Therefore, a ratio kv/kh = 0.5 was adopted in the analyses.

Development of a physically based model for the case study
Hydraulic model Geometry Figure . 6 shows a 3D schematic view of a portion of the hillslope. The hillslope is crossed by gullies (schematised as vertical wall channels in the figure) and it can be therefore subdivided into a main hillslope, a side slope and a gully. The soil formation overlies a bedrock formation. The geometrical parameters of the two slopes analysed in this work are summarised in Table 2. The values reported in Table 2 were obtained on the basis of the slope profiles derived from the 5 m × 5 m resolution digital elevation model made available for the site. Thanks to the use of a GIS software, it was easy to investigate the geometries of the slopes both in longitudinal and transversal directions. The thickness of each soil horizon was instead explored during the several surveys and sampling campaign using the open-end sampler. The latter was used to cut mini-bore holes which helped to inspect the stratigraphy at the landslides sites.
where v ! = flow velocity vector, K = hydraulic conductivity, u w = pore water pressure, γ w = unit weight of soil water and z = vertical coordinate increasing upward. The hydraulic conductivity K depends on the degree of saturation, which in turn depends on pore water pressure. The mass balance equation for liquid water can be written as follows: div where θ = volumetric water content (ratio of water volume to total volume) and t = time. By substituting Eq. 1 in Eq. 2, the water flow can be derived. The length of the slope in x′ direction (parallel to the slope) is significantly larger than the dimensions in the cross section orthogonal to the slope (plane y′z′ in Fig. . 6). According to Tarantino and Mongiovi (2003), the water flow equation in the reference system x′y′z′ (Fig. . 6) reduces to a two-dimensional form: The 2D water flow Eq. 3 was solved numerically via the FEM using the module SEEP/W of the software Geostudio. The soil constitutive functions required to solve Eq. 3 include the water retention function θ = θ(u w ) and the hydraulic conductivity function K = K(u w ).

Water retention functions
A 2-parameter Van Genuchten function (Van Genuchten 1980) was selected to simulate water retention behaviour of the three soil horizons: where α and n are empirical soil parameters. Weather conditions in Scotland made it difficult to explore a wider suction range experimentally and fitting the parameters α and n against a relatively large data set (the soil is most of the time close to the condition of zero suction). A hybrid approach was therefore adopted with one parameter determined via a pedotransfer function (n) and one parameter (α) fitted against a single experimental data point. This method falls in the class of one-point methods widely used in literature (Rajkai et al. 2004;Chin et al. 2010;Han et al. 2017;Yin and Vanapalli 2018). The water

Original Paper
Landslides retention data point was measured in field after a period of 3 days of no rain, which is quite uncommon at Rest and Be Thankful, to measure the highest suction and, hence, the widest suction range as possible. It is worth mentioning that, in the first instance, the direct measurement of a wider suction range was attempted in laboratory via tensiometers. However, the whole procedure resulted to be extremely time demanding, as the instrument took a long time to stabilise at the suction value, most probably due to the internal macroporosity of the specimens. The parameter n essentially controls the slope of the tangent to the water retention curve at the inflexion point. To derive this parameter, the water retention data points for the three horizons were first derived from the Arya and Paris (1981) pedo-transfer approach. Each set of data points were fitted using a Van Genuchten function and the average value of the parameter n was used for all water retention curves of the three soil horizons (n = 2.3).
The parameter α was then estimated by constraining the water retention function to pass through a single water retention data points derived from field measurement. At given depth, suction was measured with a high-capacity tensiometer (Tarantino and Mongioví 2003) inserted at the bottom of a 22-mm-diameter borehole. Degree of saturation was determined on sub-samples withdrawn using the open-end sampler following the procedure described in the BField Surveys^section. The three water retention curves passing through the single experimentally determined water retention data point are shown in Fig. . 7a. As consequence of the hysteresis phenomenon, the soil's hydraulic behaviour in field is far distant from the main drying-wetting curves. It follows instead several scanning paths, which, over time, are responsible of the changing of spatial connectivity within the pores and air entrapment which reduces the water content of newly wetted soils (Izady et al. 2009). Depending on the soil type, the degree of saturation at zero matric suction could vary from 15% (sandy soils) up to 35% (fine soils) between the main SWRC Lu and Khorshidi 2015; Likos et al. 2014). Seeing the nature of the soils object of this study, the degree of saturation at zero suction was then set equal to 85%. The choice of the degree of saturation at zero suction, although reasonable, is inevitably arbitrary. However, this parameter is less critical than it seems as the water flow is controlled by the derivative of the volumetric water content with respect to the pore-water pressure (see Eq. (3)) rather than the volumetric water content itself.

Hydraulic conductivity functions
The hydraulic conductivity functions for the three soil horizons were determined by combining the field-saturated hydraulic conductivity k sat measured from laboratory testing (Table 1) and the relative hydraulic conductivity k r estimated from the water retention parameters. The unsaturated hydraulic conductivity was modelled via a Mualem-Van Genuchten function (Mualem 1976).
where k s is the saturated hydraulic conductivity, S r is the degree of saturation and n is the water retention parameter appearing in Eq. (4). & Water outflow concentrated at 10 cm below the ground surface (to simulate evapotranspiration from the root system).

Hydraulic boundary conditions
& Impermeable bottom boundary (to simulate the bedrock).
The first two boundary conditions are discussed in more detail hereafter. It is worth noticing that the runoff was not directly considered in the present analysis. However, the top boundary was given a potential seepage face boundary condition. At the end of each iteration, the conditions along the boundary are 'reviewed' to assure that the pore water pressure at the top boundary never exceeds zero. In particular, when the given boundary condition (rainfall) leads to the condition u w = 0, then the software automatically switch from a flux-type boundary condition to a head type (u w = 0), hence reducing the infiltration rate. In these circumstances, the difference between the reduced infiltration rate and the original boundary condition could potentially be seen as surface runoff.
Rainfall Rainfall data were made available by Transport Scotland, who kindly shared the data collected by their instruments. The rainfall data from 1 January 2011 to 1 January 2013 were previously shown in Fig. .

5.
Potential and water-limited evapotranspiration Figure . 8 shows the evapotranspiration fluxes in the energy-limited regime (potential evapotranspiration). They were calculated using the Penmann-Monteith equation (Monteith 1965) as detailed in the Appendix. Potential evapotranspiration only occurs if the soil-plant system can deliver the water flow demanded by the atmosphere. For the case of high potential evapotranspiration rate and/or low soil moisture content, this condition cannot be met and the actual water outflow is dictated by soil-plant system rather than the meteorological conditions (water-limited regime). The reduction of water outflow in the water limited regime was modelled via a reduction function that relates the ratio between actual and potential evapotranspiration to the suction at the water extraction as shown in Fig. . 9 (Feddes et al. 1978). The approach developed by Balzano et al. (2018) was adopted here to identify the suction value s 0 . An outward flux (Fig. . 8) equal to the potential evapotranspiration was imposed to a soil column reproducing the soil profile at the Rest and Be Thankful landslide site (Fig. 10). Suction in the upper portion of the soil profile (rhizosphere) increases at a rate that becomes abruptly very high and tends to infinite after a certain time. The parameter s 0 was associated with the value of the suction where the suction rate becomes very large. Mathematically, this is defined by the vertical asymptote of the suction rate versus time as discussed in Balzano et al. (2018). The value s 0 is independent of the initial condition assumed in the water-flow analysis as discussed in Appendix.
The parameter s 1 was then fixed at 1500 kPa, which is the value most authors would suggest as associated to the permanent wilting point (Feddes et al. 1978;Ghorbani et al. 2017;Mohanty et al. 2015) The values derived for s 0 and s 1 are reported in Table 3.

Initial condition and validation of the hydraulic model
The numerical analysis requires an assumption about the initial condition in terms of pore-water pressure profile at the start of the analysis. This initial condition is unknown and in principle cannot be assumed a priori due to its significant influence on the numerical results. However, if the simulation is run for an antecedent period sufficiently long, the memory of the arbitrary initial condition is erased. Figure . 11 shows the evolution of suction at a control point for three different arbitrary initial conditions at January 2011. By February 2011, the influence of the arbitrary initial condition was already lost. The suction regime at the time of the landslides events (December 2011 and August 2012) was therefore unaffected.

Shear strength
The saturated shear strength for the top and the transition soils was determined experimentally in the laboratory on samples collected in the field using a square cutter 60-mm side. Samples were transferred to the shear box frame using a wooden piston and saturated overnight by flooding the shear box container. Specimens were then loaded in steps to 100 and 200 kPa and sheared at the shearing displacement rate of 0.033 mm/min.  The topsoil showed it to have friction angle ϕ′ = 12°and an effective cohesion c′ = 50 kPa. The high cohesion c′ is generated by the fine root system in the specimen. The transition soil exhibited a failure envelope in the Mohr-Coulomb plane (normal stress σ′ versus tangential stress τ) passing though the origin. Shear strength parameters were therefore characterised by an effective cohesion c′ = 0 kPa and a friction angle ϕ′ = 36°. The friction angle of the transition soil could not be directly measured. Seeing the very similar characteristics (i.e., grain size distribution, dry density) between the subsoil and the transition soil, a value of 36°was then considered for the subsoil as well. All the other mechanical properties were directly measured.
The shear strength τ in the unsaturated range was modelled according to Tarantino and El Mountassir (2013): where σ is the total normal stress, u w is the pore-water pressure, S r is the degree of saturation, c′ is the cohesion and ϕ′ is the (saturated) friction angle. The mechanical properties of the three soil horizons are reported in Table 4.

Landslides
It is worth noticing that the soil horizons exhibited very low values of dry density. However, this was not surprising, as permanent vegetation is expected to highly modify soil porosity and dry density due to the addition of organic matter (Udawatta et al. Udawatta and Anderson 2008;Avnimelech et al. 2001).

Geometry
The stability problem is 3-D in principle. However, failures in the transverse direction (plane y-z) were not considered due to the lower slope inclination compared to the slope dip direction. Instabilities were therefore assumed to occur only in the direction of the maximum slope (hillslope) as corroborated by experimental observations of historical landslides. Considering that the thickness of the soli cover (less than 1 m) is much   Original Paper lower that the length of the slope (a few hundredths of meters), an infinite slope model seemed appropriate to interpret hillslope movements. As a result, failure was modelled as a translational movement in the x-z plane (Fig. . 12). Factor of safety The factor of safety at any depth can be derived via the limit equilibrium method. By considering the shear strength criterion given by Eq. 8, the following equation can be derived where H is the depth of the failure surface, β is the inclination of the slope and γ is the average unit weight given by: where γ s and γ w are the unit weight of solids and water, respectively, and n * is the porosity. The bedrock interface was assumed to be rough, hence failure occurs in a soil layer just above the interface. As a result, the friction angle at the soil-bed rock interface equals the friction angle of the soil. Field survey confirmed the validity of this assumption as the bedrock was never visible at the landslides sites.
The role of the rhizosphere on the occurrence of shallow landslides Validation of the physically based model for the shallow landslide The physically based numerical model was first challenged to reproduce failure of the two historical landslides analysed in this work. The evolution of the profile of factor of safety (FoS) over time at the centre of the hillslope (axis of symmetry in Fig. . 6) is shown in Fig. . 13a, c for the landslides 'December 2011' and 'August 2012', respectively. It can be observed that the minimum factor of safety always occurs at the base of the hillslope at the interface with the bedrock. This is in accordance with the field surveys carried out after the landslide's occurrence.
The evolution of the FoS at the base of the slope over time is shown in Fig. .   The two failures are associated with heavy rainfall occurring over a period of a few days preceding the landslide event. It is interesting to observe that a heavy rainfall event also occurred in February 2011 without triggering slope failures (in both the numerical simulation and the real world). It is also interesting to observe that the same rainfall event on December 2011 triggered failure at one location only (landslide December 2011).
It appears evident that the amount of rainfall cumulated over a few days is not the only triggering factor and other aspects control the landslide mechanism. It is interesting to interrogate the physically based numerical model to explore the hydro-mechanical processes leading to failure initiation. Fig. . 14 shows that pore water pressure contours and the flux vectors for landslide December 2011 at two distinct times, February 2011 (heavy rainfall not causing failure) and December 2011 (heavy rainfall causing failure), respectively.
It can be observed that pore-water pressures on February 2011 (not causing failure) are lower (more negative) than the porewater pressures on December 2011, where pore-water pressure becomes positive at the base of the slope as shown by the formation of a phreatic surface (Fig. . 14). At the same time, it can be observed that flux vectors on February 2011 highlight a predominance of lateral flow in the rhizosphere whereas flux vectors on December 2011 shows a significant component of downward flow, which led eventually to the formation of a phreatic surface at the base of the slope.
This different hydrological responses of the slope is controlled by the antecedent rainfall, which is significantly lower on February 2011. The lower antecedent rainfall maintains low the degree of saturation and, hence, the hydraulic conductivity of the transition soil and subsoil. This promotes the diversion of rainwater ' channelled' by the topsoil (rhizosphere) on the main rainfall event and preserves the pore-water pressure in the bottom layers.
It is also interesting to compare the pore water pressure contours and the flux vectors on December 2011 for the two landslides, landslide December 2011 (heavy rainfall causing failure) and landslide August 2012 (heavy rainfall not causing failure), respectively (Fig. . 15).
It can be observed that the pore water pressure distribution and the flow direction pattern is very similar. However, a phreatic surface only forms for landslide December 2011. This is associated with the slightly lower 'drainage' capacity of the side-slope for landslide December 2011, in turn associated with the lower inclination of the side-slope (25°for landslide December 2011 against 30°for landslide August 2012).

Effect of the rhizosphere
To highlight the role of the rhizosphere in creating a subsurface lateral drainage, an exercise was developed for comparison by considering the case where rhizosphere is not present. The slope was assumed to be formed by a homogenous soil profile with hydraulic properties equal to the ones of the subsoil.
The comparison of the pore water pressure contours and the flux vectors for the two cases is shown in Fig. . 16 (by considering antecedent and incident rainfall on February 2011). It can be clearly observed that the rhizosphere diverts the incident rainfall thus preserving low water pressures in the deeper horizons. On the other hand, flux vectors are nearly vertical in the rhizosphere-less soil. This accelerates the formation of a phreatic surface at the bottom of the slope triggering premature failure of the slope.
If one compares the evolution of the FoS versus time for the rhizosphere profile and the rhizosphereless profile, it is observed that rhizosphereless profile fails on the first significant rainfall event. In other words, if the rhizosphere at Rest and Be Thankful had not acted as a lateral subsurface drainage, the soil cover would have been already swept away and one would observe the bedrock entirely exposed.

Discussion
The lesson learned from this study is that the rhizosphere 'protects' the slopes by hampering downward rainwater infiltration through the promotion of subsurface lateral flow towards the gullies. As a result, shallow landslide hazard can be potentially mitigated by enhancing the capacity the rhizosphere to act as a natural lateral drainage.
Diffuse landslides require 'diffuse' mitigation measures and, for this reason, vegetation is often seen as one of the few practicable solutions. This is also the case of R&BT where a few preliminary studies have been carried out to examine ecological mitigation measures, i.e., tree planting and other forms of re-vegetation (Winter and Corby 2012). In particular, the authors suggested a mix of native broadleaf tree and shrub species. This would give a mix of root spread and depth, including potentially to bedrock, maximising the root reinforcement effect.
The final design of the vegetation-based engineered slope would likely require further studies (including the development of a field trial) to clearly identify the hydro-mechanical processes to be 'engineered' by vegetation. In turn, this should stem from the understanding of the hillslope hydrology and its interplay with failure mechanisms.
The study presented in this paper potentially suggests that one of the criteria to be considered in the selection of the vegetation to be re-planted should stem from the role played by the rhizosphere in diverting rainwater from the deeper soil horizons. This would imply that plants with root-system architecture that enhances lateral subsurface flow should be privileged. Ghestem et al. (2014aGhestem et al. ( , 2014b have discussed the influence of plant root systems on subsurface flow and their implication on slope stability. In particular, they proposed a classification of plant root architecture with respect of the mode of subsurface flow they generate. Figure . 17 shows the different scenarios of the effects of root architecture on preferential flow according to Ghestem et al. (2011).
For example, roots that are oriented perpendicular to the slope gradient are able to capture downward water flow (Fig. .  17c) and would be ideal for plantation on the hillslope. On the other hand, tuft-root systems (Fig. . 17b) allow water infiltrating into upper soil layers and would be ideal for plantation on the side-slope.
On the other hand, taproot systems conveying water to deeper soil layers (Fig. . 17a) should be avoided as these will short-circuit water flow towards the interface between the soil cover and the bedrock cancelling any beneficial effect of the rhizosphere. It is therefore possible that remedial measures based on tree planting might have a detrimental effect on slope stability. Fig. 16 Comparison of the pore water pressure contours and the hydraulic velocity vectors on February 2011 for the case with rhizosphere and in absence of rhizosphere Fig. 17 Illustrations of different scenarios of the effects of root architecture on preferential flow. a Taproot systems convey water to deeper soil, where it can drain into cracks in the bedrock if cracks are present. However, if the soil-bedrock limit is impermeable, zones of high water pressures can be created. b Tuft-root systems allow water to infiltrate into upper soil layers. c Roots oriented perpendicular to the slope gradient capture downward water flow

Conclusion
Vegetation is often seen as one of the few practicable solutions to mitigate landslide hazard at the catchment scale. Vegetation may reinforce slopes via the anchoring of root system to deeper stable layers (mechanical effect) or by promoting the increase in suction via the transpiration process (hydrological effect). This paper has investigated a third effect thus far little explored, i.e. the role played by the high hydraulic conductivity of the rhizosphere on hillslope hydrological processes and, hence, on the occurrence of shallow landslides.
The case study of Rest and Be Thankful in Scotland has been the object of this study. Soil profile of shallow landslides has been reconstructed by inspecting scarps at the upper edge of recent landslides and hydraulic and mechanical properties of the rhizosphere and the underlying horizons have been characterised via a combination of laboratory testing and field measurements. A physically based hydro-mechanical model has been derived thereof and validated against its capability of reproducing time of occurrence of two historical landslides.
The physically based model has been then exploited to highlight the role of the rhizosphere on slope hydrology. For comparison, the case of a soil profile where rhizosphere is not present has also been considered and water flow was simulated for the case of bare and vegetated soil. Numerical simulations clearly show that the rhizosphere hamper vertical infiltration of rainwater. This allows the pore-water pressures to remain generally negative in the deeper layers at the interface with the bedrock. When prolonged rainfalls exceed the drainage capacity of the rhizosphere in diverting water flow laterally towards the gullies, downward infiltration of rainwater is promoted potentially triggering slope instability.
The lesson learned from this study is that the rhizosphere has beneficial effects on slope stability by preventing downward rainwater infiltration. As a result, shallow landslide hazard can be potentially mitigated by enhancing the capacity the rhizosphere to act as a natural lateral drainage. In turn, this implies that plants with rootsystem architecture that enhances lateral subsurface flow should be privileged when designing vegetation-based remedial measures.