Tree-root control of shallow landslides

Tree roots have long been recognized to increase slope stability by reinforcing the strength of soils. Slope stability models usually include the effects of roots by adding an apparent cohesion to the soil to simulate root strength. No model includes the combined effects of root distribution heterogeneity, stress-strain behavior of root reinforcement, or root strength in compression. Recent field observations, however, indicate that shallow landslide triggering mechanisms are characterized by differential deformation that indicates localized activation of zones in tension, compression, and shear in the soil. Here we describe a new model for slope stability that specifically considers these effects. The model is a strain-step discrete element model that reproduces the selforganized redistribution of forces on a slope during rainfall-triggered shallow landslides. We use a conceptual sigmoidal-shaped hillslope with a clearing in its center to explore the effects of tree size, spacing, weak zones, maximum root-size diameter, and different root strength configurations. Simulation results indicate that tree roots can stabilize slopes that would otherwise fail without them and, in general, higher root density with higher root reinforcement results in a more stable slope. The variation in root stiffness with diameter can, in some cases, invert this relationship. Root tension provides more resistance to failure than root compression but roots with both tension and compression offer the best resistance to failure. Lateral (slope-parallel) tension can be important in cases when the magnitude of this force is comparable to the slope-perpendicular tensile force. In this case, lateral forces can bring to failure tree-covered areas with high root reinforcement. Slope failure occurs when downslope soil compression reaches the soil maximum strength. When this occurs depends on the amount of root tension upslope in both the slope-perpendicular and slope-parallel directions. Roots in tension can prevent failure by reducing soil compressive forces downslope. When root reinforcement is limited, a crack parallel to the slope forms near the top of the hillslope. Simulations with roots that fail across this crack always resulted in a landslide. Slopes that did not form a crack could either fail or remain stable, depending on root reinforcement. Tree spacing is important for the location of weak zones but tree location on the slope (with respect to where a crack opens) is as important. Finally, for the specific cases tested here, intermediate-sized roots (5 to 20 mm in diameter) appear to contribute most to root reinforcement. Our results show more complex behaviors than can be obtained with the traditional slope-uniform, apparent-cohesion approach. A full understanding of the mechanisms of shallow landslide triggering requires a complete re-evaluation of this traditional approach that cannot predict where and how forces are mobilized and distributed in roots and soils, and how these control shallow landslides shape, size, location, and timing. Published by Copernicus Publications on behalf of the European Geosciences Union. 452 D. Cohen and M. Schwarz: Tree-root control


Introduction
Shallow landslides are hillslope processes that play a key role in shaping landscapes in forested catchments (Istanbulluoglu and Bras, 2005;Sidle and Ochiai, 2006). Many studies have highlighted the importance of roots and their mechanical properties for the stabilization of hillslopes (e.g. Schwarz et al., 2015), but usually, only basal root reinforcement is considered. When considering how roots reinforce soil, however, three different mechanisms of root reinforcement must be recognized: (1) Basal 5 root reinforcement acting on the basal shear surface of the landslide. This is the most efficient mechanism, if present. In many cases, however, this mechanism is absent because the position of the failure surface is deeper than the rooting zone. (2) Lateral root reinforcement that acts on lateral surfaces of the landslide. The magnitude of the contribution of this mechanism depends on the type of deformation of the landslide mass. If the landslide behaves as a rigid mass, lateral reinforcement may act almost simultaneously along all the edges of the sliding mass (in tension, shear, and compression). In cases where there is differential 10 deformation of the soil mass, this leads to the progressive activation of lateral reinforcement, first in tension at the top of the landslide, and then in compression at the toe at the end of the triggering. The magnitude of lateral root reinforcement depends on the spatial distribution of the root network. (3) Roots stiffen the soil mass. The presence of roots in the soil increases the macroscopic stiffness of the rooted soil mass leading to a larger redistribution of forces at the scale of the hillslope through small deformations. This mechanism increases the effects of the previous two (basal and lateral root reinforcements). 15 On top of these considerations on root reinforcement mechanisms acting on a single landslide, it is important to highlight that those mechanisms assume different meaning when considering the more global context of landslide processes at the catchment scale. Specifically, the effects of root reinforcement on landslide processes are considered limited by: (i) the magnitude of root reinforcement (a function of forest structure and tree species composition). Root reinforcement needs to reach values of the order of a few kPa in order to be significant . (ii) The heterogeneity of root distribution (tree species, 20 topography, local soil condition, etc.). Root reinforcement must be active in specific places and at specific times to have any effect on slope stability: mean values of apparent cohesion across the entire hillslope are not representative and not sufficient for considering the specifics of actual root reinforcement effects. (iii) The depth of the landslide shear surface (effects of basal root reinforcement). The deeper the shear surface is, the less important the effect of basal root reinforcement is. (iv) The length and volume of the landslide (lateral root reinforcement and buttressing/arching mechanisms and stiffening effects). The larger 25 the landslide is, the lower are the effects of lateral root reinforcement. In order to characterize the efficacy of roots for slope stabilization, a spatio-temporal quantification of root reinforcement is needed.
In view of the importance of root reinforcement and of shallow landslides to landscape evolution and to human societies, mechanistic models that include the processes linked to the triggering of shallow landslide and the influence of root reinforcement on it are needed. In the large majority of cases, slope stability models add apparent cohesion to the soil to simulate root 30 reinforcement (e.g. Milledge et al., 2014;Bellugi et al., 2015;Hwang et al., 2015). Few models includes the effects of root distribution heterogeneity (Stokes et al., 2014), and none consider the stress-strain behavior of root reinforcement and the strength of roots in compression. Recent field observations show that shallow landslide triggering mechanisms are characterized by differential deformation that indicates localized loading of soils in tension, compression, and shear (Schwarz et al., 2012a).
These observations contradict common assumptions used in models until now, yet the direct coupling of these different root reinforcement mechanisms, and their mobilization during the triggering of shallow landslides, has not yet been made.
Here we present a new model for shallow slope stability calculations that specifically considers these important effects. To fully understand the mechanisms of shallow landslide triggering, a complete re-evaluation of the traditional apparent cohesion approach is required. To do so, it is important to consider the forces held by roots in a way that is entirely different than 5 done thus far. Moreover, measurements and models indicate that the assumptions of constant elasticity and homogeneous root properties, as applied in typical finite element geotechnical model, cannot reproduce the mechanisms leading to the triggering of forested slope failures .
The SOSlope (for Self-Organized Slope) model presented here fills this gap by developing a mechanistic model for predicting shallow landslide sizes across landscapes, considering the effects of root reinforcement in a detailed quantitative manner 10 (spatio-temporal heterogeneity of root reinforcement). The SOSlope model allows to explore the activation of root reinforcement during the triggering process and helps to shed light on the contribution of roots to the slope stability. The SOSlope model is used in this work to test the following main hypotheses: -Both tensional and compressional forces resulting from mobilization of forces in the roots and the soil are efficient in stabilizing slopes but have higher effectiveness when occurring simultaneously. 15 -Weak zones in the root network (Schwarz et al., 2010b(Schwarz et al., , 2012a determine the effectiveness of root reinforcement at the slope scale if no basal reinforcement is present.
-Coarse roots dominate reinforcement and its efficacy, when present.
In what follows we first present general background on the importance of vegetation for geomorphic processes in the context of hillslopes and landslides (Section 2). We then describe the SOSlope model in details (Section 3), present the data set for 20 roots and soil used in simulations (Section 4), show and discuss results (Section 5), and synthesize a typical force redistribution process during landslide triggering (Section 6). Conclusions are given in Section 7.

Background and motivation
Understanding the role of shallow landslides to the geomorphic evolution of landscapes is of prime importance and motivates the present work. In some regions, shallow landslides are the dominant regulating mechanisms by which soil is delivered from 25 the hillslope to steep channels or fluvial systems (Jakob et al., 2005). The magnitude and intensity of these phenomena also has important societal impacts both in the long (landscape evolution and soil resource availability (Istanbulluoglu and Bras, 2005;Montgomery, 2007)) and short terms (risks due to landslides, debris flows and sediment transport, water quality, soil productivity (Wehrli et al., 2007;Hamilton, 2008)).
On long time scales, shallow landslides are important geomorphic processes shaping landscapes of both vegetated and 30 non-vegetated basins. For vegetated basins, the spatio-temporal distribution of root reinforcement has a major impact on the dynamic of sediment transport at the catchment scale (Sidle and Ochiai, 2006) and on the availability of productive soil, a key resource for human needs. At the hillslope scale, the presence of vegetation generally increases soil thickness, lowering the frequency of landsliding events but increasing their magnitudes (Amundson et al., 2015). At the catchment scale, vegetation causes slopes to steepen and sediment mobilization is then often dominated by deep landslides driven by fluvial incision (Larsen and Montgomery, 2012). The influence of shallow landslides on shaping the landscape on long time scales is, in part, masked by continuously-changing factors influenced by human activities, climate change, and other disturbances such as storms, fires, 5 etc. Under these constant disturbances soils never reach an equilibrium state that would otherwise require between 10 and 1000 years (Blume et al., 2010;Bebi et al., 2017). Nevertheless, the presence of soils on steep slopes are a necessary condition for preserving important functions of mountain environments, such as water supply, nutrient production, biodiversity, landscape aesthetics, and cultural heritage.
While soil as a resource is gaining increasing attention in the context of global sustainable development (Nature Editorial, 10 2015), risks related to shallow landslides and to processes linked to them (debris flows, bedload transport, large wood transport during floods) as well as the availability of quality water are issues that impact human societies in the short term (Miura et al., 2015), particularly in mountainous regions. Water quality is linked to shallow landslides because sediments mobilized by landslides are transported as suspended sediments in streams.
While sustainable resource management in forestry and in agriculture targets to keep the frequency of shallow landslide 15 events to pseudo-equilibrium conditions at the catchment scale and to reduce the overall erosion rate (Li et al., 2016), disturbances such as those due to human activities may lead to a rapid and dramatic increase of shallow landslide frequency and magnitude. For instance, deforestation and intensive agriculture may lead to an increase of the overall erosion rate by one order of magnitude. Marden (2012) reports that in the 17 km 2 catchment of Waipaoa (New Zealand), erosion rate increased from 2.7 to 15 Mt year −1 after deforestation and conversion of slopes to pasture land. In this new environment, shallow landslides 20 contribute ∼60% of the sediment yield of the Waipaoa river during floods and 10 to 20% of total erosion. Similar conditions occurred in the European Alps until the first half of the 20th century that led to a considerable increase in erosion rates (Mariotta, 2004). Meusburger and Alewell (2008) reported that, in a catchment in the central Alps, the increase of landslide area by 92% within 45 years was likely due to dynamic factors like climate and land-use changes and had a decisive influence on landslide patterns observed today. 25 Risks due to shallow landslides are associated with different types of phenomenon ranging from hillslope debris flows (example of process causing a direct risk to infrastructures and individuals) to various channel processes such as large sediment transport during floods, wood debris transport, channelized debris flows, etc. (examples of processes causing an indirect risk to infrastructures and individuals). It is estimated that landslides triggered by heavy rainfall cause damages upwards of several millions each year and more than 600 fatalities per year (Sidle and Ochiai, 2006). 30 Next to the constellation of factors well-known to influence the triggering of shallow landslides, vegetation has been recognized to play an important role (Sidle and Ochiai, 2006;Schwarz et al., 2010c;McGuire et al., 2016) and its function is considered an important component of ecosystem services provided in mountain regions. The importance of the effects of vegetation is, in some cases, recognized at a political level. For instance, the global forest area managed for protection of soil and water is 25% of all global forested areas (Miura et al., 2015). In Switzerland, protection forests occupy more than 50% of all 35 forested areas (Wehrli et al., 2007). Moreover, bio-engineering measures are often considered an important part of integrated risk management and disaster risk reduction strategies. The management of such protection forests and bio-engineering measures needs quantitative tools to optimize the effectiveness of such important ecosystem services for society. The formulation of such tools needs to be based on quantitative methods applicable to a large range of situations. Moreover, these methods need to consider different time and spatial scales at which vegetation influences processes. To put the motivation for the present 5 work in the appropriate context, we shortly summarize the effects of vegetation on long and short term geomorphic processes.
In the long term, the presence of vegetation (i) increases soil production rates through mechanical and chemical processes (Wilkinson et al., 2005;Phillips et al., 2008) (100-1000 years); (ii) increases soil residence time on hillslopes due to root reinforcement and protects against runoff erosion (Istanbulluoglu and Bras, 2005) (10-100 years). Note that in the case of natural or human driven disturbances, the response time of the system (i.e. root decay) is of the order of a few years (Vergani 10 et al., 2016); (iii) enhances soil diffusion rates on hillslopes due to tree wind-throw (Pawlik, 2013;Roering et al., 2010), root mounds (Hoffman and Anderson, 2014), and biological activity (Gabet and Mudd, 2010) (100-1000 years).
In the short term, vegetation mainly influences root reinforcement and regulates water fluxes. At the hillslope scale, the hydrological effects of vegetation are assumed to play a small role on slope stability compared to the contribution of root reinforcement (Sidle and Bogaard, 2016;Sidle and Ziegler, 2017). At the catchment scale, however, the regulation of water 15 fluxes may have important implications for the stability of those slopes that drain large areas, particularly for short and intense rainfall events.
Root are considered the hidden half of plants due to the difficulties in characterizing and quantifying their distribution and mechanical properties. In slope stability, the process of root reinforcement remains hidden because direct observations have not yet been made on steep hillslopes. Field and laboratory experiments (e.g. Zhou et al., 1998;Ekanayake and Phillips, 20 1999;Roering et al., 2003;Docker and Hubble, 2008) generally explore only a small part of the complex root reinforcement mechanisms.
Methods for the quantification of different types of root reinforcement mechanisms went through a succession of models in the last few decades, starting with the assumption of the simultaneous breakage of all roots (Wu et al., 1979;Waldron and Dakessian, 1981) to the application of fiber bundle models that consider the progressive failures of roots (Pollen and Simon, 25 2005;Schwarz et al., 2010a;Cohen et al., 2011). Fiber bundle models may be differentiated on the basis of the type of loading, whether it is by stress (Pollen and Simon, 2005) which does not allow for the calculation of displacement, or by strain Cohen et al., 2011), which does. We enumerate below some aspects of root reinforcement models important for slope stability.
1. Breakage versus slip out. Field observations show that in tree-root bundles, the dominant failure mechanism of roots is 30 by breakage (Schwarz et al., 2012a). Slippage is limited to small roots that usually contribute only a small fraction of the total root reinforcement. For this reason, numerical models usually assume that all roots fail by breaking Cohen et al., 2011). 2. The contribution of root reinforcement must be differentiated between different types of stress conditions: tension, compression, and shearing. While most of the literature has focused on the shear behavior of rooted soils (e.g. Docker and Hubble, 2008), some works have investigated the contribution of root reinforcement under tension (Zhou et al., 1998;Schwarz et al., 2010aSchwarz et al., , 2011 and compression . In general the contribution of maximum root reinforcement under tension and shearing is of the same order of magnitude, whereas under compression the contribution 5 of roots is about one order of magnitude smaller. However, roots contribute significantly to increase the stiffness of soil under compression. This may play an important role in the re-distribution of forces during the triggering of a shallow landslide . 3. The mechanical interactions of neighboring roots in a bundle are usually neglected. Giadrossich et al. (2013) showed with laboratory experiments that the failure mechanisms of single roots is influenced by neighboring roots only at high 10 root density that are usually reached only near tree stems (0-0.5 m).
4. The mechanical and geometrical variability of roots was recently considered using survival functions ) that represent the complexity of several factors contributing to the variable stress-strain behavior of roots. Specifically, these factors are root tortuosity (Schwarz et al., 2010a), root-soil mechanical interactions , and position of root breakage along the root. Pulled roots break at different distances from the point of force application 15 because of branching, root geometry, changes in root diameter due to soil properties, presence of stones, etc.
5. The spatial and temporal heterogeneity of root reinforcement is related to several factors such as topography, soil water content, soil disturbances, resistance and resilience of forest cover to disturbances, animal browsing, etc. (Schwarz et al., 2012a;Vergani et al., 2016) 3 The SOSlope model 20

General framework
SOSlope is a hydro-mechanical model of slope stability that computes the factor of safety on a hillslope discretized into a twodimensional array of blocks connected by bonds. Bonds between adjacent blocks represent mechanical forces acting across the blocks due to roots and soil (Cohen et al., 2009). These forces can either be tensile or compressive depending on the relative displacements of the blocks. A digital elevation model (DEM) is used to divide the hillslope into squares in plan view where 25 the centers of the squares are points of the DEM (Fig. 1). Three dimensional blocks are created by extruding the squares to the bottom of the soil layer along the vertical. The center of mass of a block is connected to the four lateral blocks by four force bonds (Fig. 1). Initially, bond forces between blocks are set to zero. Rainfall onto the slope will increase the mass and decrease the soil shear strength of the blocks. At each time step, the factor of safety is calculated for each block using a force balance (resistive force over active force, see equations below). If the factor of safety of one or more blocks is less than one, 30 those blocks are moved in the direction of the local active force (defined below) by a predefined amount (usually 0.1 mm) and the factor of safety is recalculated for all blocks. Because of the relative motion between blocks that have moved and blocks that remain stationary, mechanical bond forces between blocks are no longer zero and the force balance changes. This relative motion triggers instantaneous force redistributions across the entire hillslope similar to a self-organized critical (SOC) system of which the spring-block model (Bak et al., 1988;Hergarten and Neugebauer, 1998;Cohen et al., 2009) is a subset. Looping over blocks and moving those that are unstable is repeated until all blocks are either stable (factor of safety greater or equal to 5 1) and the system reaches a new equilibrium, or some blocks have failed (their displacements are greater than some set value, usually a few meters), triggering a landslide.

Factor of safety
The factor of safety for each block is calculated as the ratio of resistive to active forces. Resistive forces include the soil basal shear strength and the strength of roots that cross the basal slip surface, assumed to be located at the bottom of the soil layer. 10 The active forces include the gravitational driving force due to the soil mass and the push or pull forces between blocks that include the effects of soil and root tension and compression. These later forces are the bond forces between the blocks described above. Including all these forces in a force balance yields the factor of safety where F s is the soil basal resistive force that includes soil cohesion and friction, F r is the basal root resistance, F d is the driving force vector due to gravity, and F j , j = 1, · · · , 4, are the 4 bond vector forces that quantify soil and root tension or compression between the block and its four neighbors. The vertical bars in the denominator denote the norm of a vector. This factor of safety 5 is calculated for each block but an index for the block number is not included so as not to clutter the equations.
Soil basal resistance is where A is the surface area of the block along the failure surface and τ b is the basal shear stress (described below). In the present model, we set F r = 0, focusing on lateral root reinforcement. This is justified in many cases where the depth of the slip surface 10 is 1 meter or greater and very few roots are present (e.g. Bischetti et al., 2005;Tron et al., 2014). Basal root reinforcement can easily be added using a formulation similar to lateral root reinforcement (discussed below) with values of root reinforcement a function of the shear displacement and the density of roots crossing the slip surface.
The driving force is where γ is the specific weight of the wet soil, D is the depth to the shearing surface, perpendicular to slope, andt is the unit tangent to the slope in the direction of the maximum slope. The specific weight of the wet soil is calculated based on water content and solid fraction, i.e., where ρ s and ρ w are the solid (grain) and water densities, respectively, φ s is the solid volumetric fraction, θ the volumetric 20 water content, and g is gravity.
Bond forces are given by where F soil j and F roots j are the soil and root components of the four bond forces, respectively, andb j are unit vectors along the bond axes pointing outward of the block. These quantities are detailed below.

Bond forces due to roots
The force in bond j between a block and its neighbor due to roots (F root j ) depends on four factors: the root density and the root-diameter distribution at the bond center, the strength of roots which depends on root diameter, and the change in length (elongation) of the bond with respect to its initial length. Changes of root density with depth (e.g. Bischetti et al., 2005) are not taken into account. This force is computed using the Root Bundle Model (RBM) of Schwarz et al. (2013) with Weibull statistics, called RBMw. For the sake of completeness, the full details of the model are given below.

Root density and root-diameter distribution
Roots are binned according to their diameters in 1-mm size bins from 0.5 mm to an upper limit given by data. A bin is usually 5 referred to as a root-diameter class, with φ i denoting the mean root diameter of class i, i = 1, · · · , i max . At each of the four faces of a block, the total number of roots for each root-diameter class i that crosses a face j is the sum of the number of roots for that root-diameter class from each surrounding tree in the stand. Summing roots from each tree implies no competition for resources. Following the empirical model of Schwarz et al. (2010a) in its version described by Giadrossich et al. (2016), the number of roots depends on the distance of the face center to the tree trunks, the tree trunks diameters, and the tree specie. 10 For simplicity all trees in the stand are assumed to belong to the same specie. The model assumes a linear allometric relation between trunk size and root density, a power-law decay of root density with distance from the tree trunk, and a logarithmic decrease of root density with root-diameter size. The number of roots of class diameter φ i crossing face j is where A j is the surface area of face j, T is the number of trees in the stand (more specifically the number of trees whose roots 15 reach face j of the cell), and ρ j k is the density of fine roots of tree k for face j. This later quantity is given by where N k , the total number of fine roots of tree k, is d max k , the maximum rooting distance for tree k, is and φ max k , the maximum root diameter class of tree k, is In these equations, φ o = 1 mm is the size of the smallest root diameter class, d j k is the distance between face j and tree k, and φ k is the tree diameter (usually diameter at breast height or simply DBH). This model contains four fitting parameters (µ, η, 25 ψ, and γ) that must be determined from data Schwarz et al., 2016).

Root mechanical forces
Roots are assumed elastic in both tension  and compression . The linear elastic force in a root is expressed using a spring constant (i.e., Hooke's law) that depends on the root diameter class. For a root in diameter class i on bond j, that elastic force is 5 where the superscript E indicates either tension (E = T) or compression (E = C) and x j is the elongation of the bond from its initial length (positive for tension, negative for compression). Based on data (e.g. Schwarz et al., 2013Schwarz et al., , 2015 we assume the spring constant depends linearly on root diameter, i.e., with k E 0 and k E 1 constants to be determined from data. Other formulations based on a power-law relation can also be used 10 Giadrossich et al. (2016).
The variability of root bio-mechanical properties (e.g. maximum tensile or compressive strength, elastic moduli in tension or compression) due to the presence of biological or geometrical weak spots is handled probabilistically. The probability of failure of a root in tension (or in compression) is captured by multiplying the elastic force by a Weibull survival function (S) that depends on a dimensionless bond elongation. Then, the total root-bond force is obtained by summing over all roots of each 15 diameter class, i.e., where N j φi is given by Eq. 6, F E i,j by Eq. 11, and 20 where λ E and ω E (E = T or C) are two scale and two shape parameters to be determined from field or laboratory experiments (see Schwarz et al., 2013Schwarz et al., , 2015. F E i,max is the maximum force held in a root at breakage (in tension) or at the critical buckling condition (in compression, see Schwarz et al., 2015) for a root of diameter φ i and is given by the commonly used power-law equation with α E the power-law exponent and F E o a pre-exponential factor for tension or compression (E = T or C). The scaling of the displacement with the maximum strength of a root eliminates the effect of root diameter on maximum displacement. Similarly, the parameter λ E scales the root strength variability to the root diameter. Equation 13 has a maximum (F roots j,max ) called the maximum root reinforcement and occurs at a bond elongation x j,max .
3.4 Bond forces due to soil 5 The soil bond force (F soil j , Eq 5) depends on whether the soil is in tension or in compression. For tension, we assume that resistance scales with soil apparent cohesion (including the effects of suction stress for unsaturated soils) as a function of displacement using a logarithmic function (Win, 2006) where c a is the apparent cohesion, ε T max is a strain threshold above which soil loses any tensional resistance, and L j is the 10 length of bond j. In compression, following the work of Schwarz et al. (2015) we assume that the soil compressional resistance is mobilized across the shear plane that forms during the failure of a downslope wedge, similar to the earth pressure force in the geotechnical engineering literature that develops during passive state when a retaining wall moves downslope toward the adjacent backfill (e.g. Milledge et al., 2014). According to Schwarz et al. (2015), the mobilized force on the downslope wedge scales with the maximum passive earth pressure force F p and with the displacement, i.e., where and K pγ and K pc are the passive earth pressure coefficients due to soil weight and to cohesion, respectively, obtained from a fitting of equations given in Soubra and Macuh (2002), c is effective soil cohesion, and P w1 and S w2 are the Weibull cumulative 20 density and the Weibull survival functions, respectively, given by and with µ 1 , κ 1 , µ 2 , and κ 2 four parameters determined from compression experiments. The first Weibull function, P w1 , serves 25 to scale the maximum passive earth pressure force with displacement during initial block motion while the second one, S w2 , reduces that same force as the wedge is overridden by the block and the failure surface area of the slip plane decreases (see Schwarz et al., 2015, for details). We neglect the active earth pressure force on upstream faces of cells because the magnitude of the active force is small in comparison to other forces.

Hydrological triggering
Rainfall-triggered shallow landslides can fail under saturated conditions during increases of pore-water pressure and/or loss of suction under unsaturated conditions (Lu and Godt, 2013). Our objective here is not to reproduce the detailed physical mechanisms by which changes in subsurface hydrology trigger a landslide but to develop a simple empirical model that realistically mimics observed changes in pore-water pressure and water content during rainfall infiltration. Although diverse hydrologic 5 triggers have been observed and described (e.g. Reid et al., 1997;Iverson, 2000), here we use, as a representative example for the hydrological conditions triggering a shallow landslide in our model, pore-pressure measurements during the artificial triggering of the Ruedlingen shallow landslide experiment in Switzerland (Askarinejad et al., 2012;Lehmann et al., 2013).
Data from Lehmann et al. (2013) indicate that high pore-water pressures were attained relatively quickly and remained steady across the slope long before failure occurred, and that the decrease in the standard deviation of the water saturation prior to 10 failure indicated an increase in the connectivity of water saturated regions that reduced soil shear strength across the full length of the slip surface leading to failure. Other data in different localities (e.g. Matsushi et al., 2006;Bordoni et al., 2015) have also shown high, steady pore-water pressure prior to failure. Because our model focuses on the effects of roots and soil strength on slope stability rather than on the details of hydrologic triggering, we choose a simplified, empirical, dual-porosity model for our slope hydrology. Our objective is only to reproduce reasonable pore-water pressure distribution and water content evolution in 15 both the matrix and the preferential flow domains, but not to model the physics of evolving subsurface hydrology. The model embodies the rapid increase in positive pore-pressure in a preferential flow domain (representing macropores) and the slow decrease of suction in the soil matrix caused by slow water transfer from the macropores to the matrix. This decrease of suction is the equivalent of the increasing connections of water saturated regions represented by the decrease in the standard deviation of water saturation observed by Lehmann et al. (2013) that eventually caused slope failure in the Ruedlingen experiment. 20 We assume that water flow in soils during a rainfall event is a combination of slow matrix flow (also called immobile water with capillary number lower than 1) and fast preferential flow (mobile water, capillary number higher than 1) (Sidle and Ochiai, 2006;Beven and Germann, 2013). While slow matrix flow influences the change in suction stress, the fast preferential flow directly influences pore-water pressure in the macropores. Our formulation of this concept is empirical and is a simplification of the more common dual-porosity models that employ two flow equations (e.g. Richards' equation) that exchange moisture 25 between the two domains, and mixture equations for water content, hydraulic conductivity, rainfall partitioning based on the volumetric ratio of the fast and slow flow domains (e.g. Gerke and van Genuchten, 1993;Shao et al., 2015). In accord with continuum mixture theory for effective stress (e.g. Borja and Koliji, 2009), we write the mean pore-water pressure of the soil (matrix + macropores),p, as 30 where ψ 1 and ψ 2 are the pore fractions along the potential failure surface of the landslide of the matrix and the macropores, respectively (volume of pore in matrix or macropores over total pore volume, with indices 1 for matrix and 2 for macropores) with ψ 1 + ψ 2 = 1, and wherep i , i = 1, 2, are the matrix and macropores intrinsic mean pore pressures. Pore fractions ψ 1 and ψ 2 s are related to the partial porosities of the matrix and the macropores, φ 1 and φ 2 , respectively, by where φ 1 + φ 2 = n, n the total porosity of the soil. The solid volume fraction of the matrix (macropores have only pore space) The superscripts and subscripts in these equations and in equations below refer to partial and intrinsic quantities, respectively. Partial and intrinsic water content of the matrix and macropores are related as follows: where θ 1 and θ 2 are the partial water contents of phase 1 and 2 (volumetric water content of phase 1 or 2 over total soil volume) and θ 1 and θ 2 are the intrinsic water contents of each phases (volumetric water content of phase i over volume of 10 phase i, i = 1, 2). At saturation θ 2 = φ 2 since the macropore phase contains only void space and thus θ 2 = 1. The total water content of the soil is and is used in Eq. 4 to compute the soil specific weight. Equations similar to Eq. 26 can be written for saturated and residual water contents. 15 We assume that the time evolution of the intrinsic pore-water pressure in the macropores,p 2 , and of the partial water content in both the macropore (θ 2 ) and the matrix phases (θ 1 ) can be modeled using cumulative distribution functions. For the macropore phase, we writē where p max is a constant here but ultimately depends on rainfall infiltration rate and upstream contributing area (e.g. Montgomery and Dietrich, 1994), t is a dimensionless time, F is the normal cumulative distribution function with mean µ p and standard deviation σ p , and θ r 2 is the intrinsic residual water content for the macropores (we have used the fact that the intrinsic saturated water content θ s 2 = 1 since macropores have no solid fraction). For the water content in the matrix we assume that where θ o and θ s are the soil initial and saturated water contents, respectively, and F fold is the folded normal cumulative distribution with mean µ θ and standard deviation σ θ . The pore-water pressure in the matrix is given by (Borja and Koliji, 2009) where p 1 is the intrinsic pore-water pressure in the matrix and S e is the equivalent degree of saturation (also called effective saturation) in the matrix. Following Lu et al. (2010), we have used the equivalent degree of saturation (S e ) in Eq. 30 instead of the more commonly used degree of saturation. Under unsaturated condition,p 1 is a matrix suction stress (Lu et al., 2010). The 5 equivalent degree of saturation in the matrix is defined as where θ r 1 and θ s 1 are the intrinsic residual and saturated water content of the matrix phase with θ s 1 = φ 1 . Using Eqs. 24 and 25, and equations for the residual and saturated water content equivalent in for to Eq. 26, Eq. 31 can be rewritten as 10 where θ 1 is given by Eq. 29. Using van Genuchten formulation (Van Genuchten, 1980), we can write the suction stress as (Lu et al., 2010) p where and α vg and n vg are the soil parameters.
Pore-water pressure in the macropores (Eq. 27), matrix water content (Eq. 29), matrix suction (Eq. 33), and mean pore-water  Table 1. The standard deviations are chosen so that macropore water pressure reaches its maximum before matrix water content, to mimic, but not reproduce, the behavior observed by Lehmann et al. (2013).

Basal shear stress
Basal shear resistance along the slip surface is calculated using the Mohr-Coulomb failure criterion including contributions from both the suction stress and the pore-water pressure using the mean pore-water pressurep of Eq. 22, i.e., where σ n is the normal stress and φ is the soil friction angle.

Soil
Mechanical soil parameters from Schwarz et al. (2013Schwarz et al. ( , 2015 and other parameters used in simulations are listed in Table 2.

Roots
Model parameters for roots (Table 3) are taken from field and laboratory data of Schwarz et al. (2010aSchwarz et al. ( , 2012bSchwarz et al. ( , 2013Schwarz et al. ( , 2015 for Picea abies (Norway spruce). Figure 4 shows root reinforcement as a function of bond elongation (both in tension and compression) for four values of tree diameter (DBH, diameter at breast height) and for three distances (d) from the tree trunk  lateral root reinforcement upwards of 10's of kPa is typical (Schwarz et al., 2012b). In tension, root reinforcement becomes negligible once the bond has stretched over 0.1 m, regardless of the distance from the tree trunk. In tension, the bond elongation 5 over which reinforcement is active depends on the distance from tree and range from 0.15 m close to the tree trunk to about 0.05 m at 2.5 m distance from the tree trunk.

Results and discussion
To  Table 2.
Negative values of displacement indicate compression.      During loading, cells in the clearing move downhill more than cells in the stand (Fig. 6a-d). A discontinuity in displacement appears near the top of the clearing. This gap, 12 m long and slope parallel, occurs where the surface slope is about 0.62 (ca. 32 • ). This gap represents the formation of a vertical tension crack at the upper edge of a soil slip that has yet to fail completely.

Displacement and force redistribution
With increasing loading, displacement across the crack grows to exceed 1 m prior to failure (Fig. 6c). Although this crack is in the clearing in a zone devoid of trees, a few small roots from trees above the crack are present and extend across this vertical 5 tension crack (see Fig. 5b). Cells above the crack show barely perceptible displacements (< 0.1 m). The situation is different in the forested area where, up to failure, displacement is significantly smaller (about 10 times smaller), uniform (no discontinuity), highest in the steepest portion of the slope (not visible in Fig. 6), with no evidence of a crack forming in the upper part of the slope. The slope in the stand remains stable after the clearing fails for the remaining of the simulation (2200 minutes). In the forested area, cells that have undergone displacement extend further uphill than in the clearing. We attribute this effect to the 10 connected root system of trees that activates tensional forces uphill and pull rooted cells downhill. These tensional forces are absent in the clearing due to lack of roots and negligible soil tensional strength. loading as the slope slowly slips downhill (Fig. 6i-k). Root tension perpendicular to slope is highest on both edges of the vertical crack. This is where the largest displacements are observed generating the highest tensional forces in the roots. In that zone, tension in roots reaches almost 20 kN just before the clearing fails (Fig. 6k). Simultaneously, some roots of trees in the lower part of the slope are in compression, relieving some of the compression in the soil. Across-slope (also referred to as lateral or slope-parallel) root forces are shown in Fig. 6m-p. Downward motion of soil in the clearing causes a lateral tension in roots that span the transition zone from clearing to forested area. This zone is about 6-10 meters wide. It is across this boundary that displacement gradients are high and across-slope root forces highest. The lateral tension increases up to about 6 kN with increasing downhill motion of the clearing and stays high after failure because the relative downslope displacement of cells across the slope remains. clearing (pink symbols, y = 5 m). With increasing load, both tension and compression in the slope increase. Tension is highest where root density is highest (x = −12) and slightly lower at the edge of the clearing (x = −9 m). At t = 1200 min, the edge of the clearing has formed a crack and forces downhill of that crack are now in compression. Roots that cross the crack at the edge of the clearing are now broken and no longer provide any tensional resistance (pink symbols in Fig. 7f).
Bonds that were in tension in the upper part of the slope at t = 900 min are now in compression owing to the failure of 5 roots across the widening crack near the edges of the clearing. The clearing is now entirely held by compressive forces and by lateral (across-slope) tensile forces shown in Fig. 7i-l. These lateral forces are due to root tensile strength and are highest near the transition from forest to clearing (pink symbols, x = −9 m) where the relative downslope displacement between adjacent cells is highest. Along the first row of trees (red symbols), cells that host a tree have larger values of across-slope tensional forces than cells that do not giving rise to a saw-tooth pattern of tensional force. In the clearing (black symbols), positive 10 lateral tensional forces are entirely due to soil apparent cohesion which reaches values of almost 1 kN. With increasing load and decreasing soil shear strength due to increasing mean pore-water pressure, the clearing eventually fails at t = 1359 min but the forested area remains stable for the remaining of the simulation (up to 2200 min).
Results from this simulation demonstrate that maximum tensional and compressive forces in rooted slopes do not contribute simultaneously and equally to the stability of the slope during the initiation of a shallow landslide. Roots provide reinforcement in tension. This tensional root force can disappears once displacement across a vertical crack becomes sufficiently large. In our example, this occurs when the crack grows to about 0.1 m (see Fig. 4). Compression is higher in the clearing (no roots) than in the vegetated area. Where present, when slope-perpendicular root tensional reinforcement is eliminated, soil stability 5 is entirely accommodated by soil compressive resistance and by lateral tension held by roots. Lateral root forces provide additional stability to the clearing by redistributing slope-perpendicular forces laterally across the slope. The clearing fails when soil strength at the base can no longer be held by the combination of the lateral root bond forces and downslope soil compression, and compression in the soil exceeds the maximum strength.
We can summarize the redistribution of forces during the loading of a rooted hillslope into three distinct phases: 1. Increasing load and weakening of soil strength along the basal failure plane (not shown) without any soil motion (factor of safety above 1).

Initiation of downward motion after some cells reach critical condition (factor of safety equal to 1). Force redistributions
(compression in soil, tension and compression in roots) prevent the slope from failing. These forces increase with increasing load and increasing mean pore-water pressure (e.g. Fig. 7, t = 900 min). The culmination of slope-perpendicular 15 tensional forces across the crack (t = 900 min, Fig. 7e) occurs with: (1) less-than-maximum compressive forces in the lower-half of the slope, and (2) lateral tensional forces activated at the edge of the forested area (Fig. 7i).
3. Culmination of compressive forces leading to failure when exceeded (t = 1359 min, Fig. 7, last column). This occurs after tensile, slope-perpendicular forces due to roots are lost across the vertical crack and when lateral root tensile forces reach their maximum values. 20 The timing and duration of these three phases will vary with soil mechanical properties, slope inclination, slope morphology, root distribution, and hydrology, resulting in an increase or decrease in the stability of the slope. These three phases of force redistribution are used as criteria to define the triggering of a landslide. In civil engineering, calculations using infinite slope analysis, for example, must yield a factor of safety greater than 1 for the slope to be deemed stable. Any values below one imply an unstable slope with the possibility of a landslide, even if slope motion subsequently stops with no occurrence of a 25 runout. This definition of a landslide corresponds to the second phase of force redistribution where motion has initiated but complete failure has not yet occurred. Many such occurrences of a failed landslide (at least temporarily) exist; one is shown in Fig. 8. In risk analyses, or when studying geomorphological processes, a landslide occurs by definition only when the soil mass fails completely and is followed by a runout, corresponding to the third phase of our force redistribution process. In that case, the transition from phase 2 to phase 3 and the accompanying redistribution of forces, is the critical process. 30 Changes in the values of the factor of safety (FS) over time help understand the processes of landslide triggering and illustrates the three phases of landslide initiation and force redistribution. Figure 9 shows the evolution of displacement, the factor of safety, and the mean pore-water pressure with time at the center of the clearing. Initially, FS is larger than 1 and decreases with increasing mean pore-water pressure up until about 400 minutes. This corresponds to the phase 1 described above. Beyond 400 minutes, the value of FS oscillates rapidly just above the value of 1. These oscillations, which last until failure, correspond to the critical state of a self-organized system before global failure (e.g. Bak et al., 1988). Here, these oscillations correspond to our phases 2 and 3. During this critical period, the number of cell moves (not shown) increases dramatically as a results of force redistribution between bonds of connected cells and as the number of cells with factor of 5 safety less than 1 increases. The increase of the number of redistribution with loading is similar to the process of avalanching in load-controlled self-organized systems like fiber bundle models (Cohen et al., 2009;Lehmann and Or, 2012). The increase in force redistribution across the slope corresponds to the progressive slope failure stage of coalescence of local failure surfaces that eventually leads to global failure (Petley et al., 2005;Cohen et al., 2009). This is equivalent to our phase 3.
The decrease in the factor of safety is linked to the increase in mean pore-water pressure in the soil (Fig. 9). A detailed 10 analysis of how hydrology impacts slope stability is beyond the aim of this paper. Here we wish to point out that our simple dual-porosity model, with the coexistence of pore-water pressure in the macropores and suction stress in the matrix (see Fig. 2), is realistic and can model a wide range of hydrological situations that can lead to shallow landslide triggering. In the simulation shown in Fig. 9, there is an imperceptible increase in the factor of safety during the first phase of the simulation until about 100 minutes. This increase is due to the increase (in absolute value) of the suction stress that increases the soil apparent cohesion 15 (see Fig. 2). The increase in pore pressure after about 200 minutes causes the soil to weaken with an associated decrease in the factor of safety eventually leading to the critical state (FS close to 1). Decrease of matrix suction linked to flow of water from the macropores to the matrix increases the mean pore-water pressure (Fig. 2) and eventually causes soil to weaken sufficiently for a landslide to occur. Depending on the application of the model and on the local hydrological properties, choices of different values of hydrological parameters than those used in this example could lead to different hydrological triggering. For example, triggering could be due to the rapid increase in macropore water pressure and the saturation of the soil from top to bottom with little time for changes in matrix pore pressure to occur. In our example, preferential flow paths lead to local increases of pore-water pressure that, in combination with a loss of suction stress in the soil matrix, resulting in a critical drop of soil 5 shear strength typical of forested soils on compacted bedrock (Lehmann et al., 2013). Yet in another situation, high pore-water pressure can originate from ephemeral springs or water exfiltration from fractured bedrock (Montgomery and Dietrich, 1994).

Effects of root tensile and compressive strength
Our results show that force mobilization and redistribution in the soil and in the root system during the triggering of a shallow landslide is a complex process. Our model can be used to investigate the effects of the various components of the bond 10 force system (roots and soil) on the dominant reinforcement mechanisms (tension or compression, lateral or downslope) and how these forces control the stability of the slope. Understanding which of these forces control slope stability under certain conditions is important for making appropriate simplifications when the full level of details is not needed or not known. Figure 10 shows the displacement for 9 hillslope simulations where trees, spaced 3 meters apart, cover the entire slope. In The slope behaves differently depending on the tree size and the type of root reinforcement. Root reinforcement for the 50-cm diameter trees is sufficiently large that the slope does not fail regardless on the type of root reinforcement (tensile, compressive, or both). For the 40-cm diameter trees, there is a threshold: the tensile strength of roots is needed to keep the slope stable. Without root tensile strength (compression only), the slope fails (Fig. 10b). Finally, for the 30-cm diameter trees, all root reinforcement configurations lead to slope failure, but at different times, with compression-only roots failing first and roots with both compression and tension last.
Results indicate that roots with only tensile strength limit downward slope slip under loading and delay slope failure more than roots that have only compressive strength. Roots that have both tensile and compressive strength offer the best protection against slope motion and slope failure. Neglecting root compression in the simulations results in only a couple of centimeters 5 difference in slope displacement or less than one hour in the timing of the landslide. Neglecting root tension, however, can result in predicting a false slope failure. Also, neglecting tension misses the jump in displacement during the early initiation of the landslide when roots across the tension gap in the upper part of the slope fail under tension (see Fig. 10c). Note also that when roots across the vertical crack fail in tension (as is the case for the 30-cm diameter trees), the slope eventually fails. This appears to be the case for all simulations we tested. However, simulations with roots that do not break across the widening 10 crack do not necessarily remain stable over the duration of the simulation (2200 min). Figure 11 illustrates the conditions under which the slope fails for the different tree-size diameters and root-strength configurations. Each graph in Fig. 11 shows a bond force along the downslope section at the center of the slope (x = 0). Downslope root bond force along that section (Fig. 11a, d,g) indicates that when roots have no tensile strength (C only), roots in the lower section of the slope bear a higher compressive load. Similarly, roots that have no compressive strength (T only) bear higher 15 tensile loads in the upper part of the slope, but only slightly higher than roots that have both tensile and compressive strength (T + C). As expected, roots of larger trees can bear higher tensile and compressive forces owing to higher root densities and more roots of larger diameters. Downslope soil bond forces (Fig. 11b,e,h) indicate that soils in slopes covered by smaller trees must take more of the compressive force caused by the slope downhill motion. For the 30-cm diameter trees (Fig. 11h), the soil compressive force eventually reaches its maximum value and the slope fails, regardless of the root configuration because 20 roots hold only a small fraction of the tensile or compressive resistance that helps maintain the slope stable: roots are too few and too small for this tree size. This is also the case for the 40-cm diameter trees when roots have only compressive strength ( Fig. 11d,e). Because roots do not hold any tension in the upper part of the slope, and root compression is insufficient to support much load, soil bond in compression eventually reaches a maximum and the slope fails at t = 1960 minutes. Figure 11c,f,i shows the lateral root force across the slope. This force is one order of magnitude smaller than the downslope force and has 25 only a limited role in the slope stability for the cases shown here. These simulations indicate that downslope root and soil forces control slope stability which is regulated by the maximum soil compression. Roots can reduce soil compression by taking up some of the force in tension in the upper part of the slope, preventing or delaying failure. Root compression alone is insufficient to offset soil compression in the lower part of the slope.

Effects of weak zones 30
The structure of the stand (dimension, density, and relative position of trees) plays an important role on root reinforcement and slope stability. Moos et al. (2016) found that susceptibility to landslide was higher in plots with longer downslope gaps in the tree stand and in locations where the distance to nearby trees was higher. Inversely, Moos et al. (2016) also found that susceptibility to landslide was smaller where root reinforcement, based on tree diameter and distance from tree, was high. simulations that did not fail, or at the time step just prior to failure for those that did (see Fig. 10).
Weak zones, zones with low values of root reinforcement, can serve as initiation points for slope movement and control the location and size of a landslide (Schwarz et al., 2010b(Schwarz et al., , 2012a). An example of a weak zone where a soil slip initiated is shown in Fig. 12. Roots around a tree provide sufficient stiffness to make the soil around the tree behave as a rigid body. The zone in between the tree and its neighbors does not provide sufficient root reinforcement and a gap opens as a result of loading (here rainfall). Here, we explore independently how tree size (diameter) and tree spacing can affect landslide initiation and hillslope 5 stability.

Tree diameter
Our base scenario is the simulation presented earlier with trees 50 cm in diameter spaced 3 meters apart on a sigmoid hillslope with a 20 × 50 m 2 clearing in the center. Five other simulations were run with the clearing area planted with trees of diameter 10,20,30,40, and 50 cm, all spaced 3 meters apart as in the forested area surrounding the clearing. These 6 simulations are 10 referred to as 50/0, 50 /10, 50/20, 50/30, 50/40, and 50/50, where the first and second number indicate the stand tree diameters and the tree clearing diameters, respectively. Figure 13 shows the computed factor of safety and displacement at the center of the slope for these 6 simulations. The 50/10 simulation fails earlier (t = 1266 min) than the 50/0 simulation (t = 1359 min). The time evolution of the factor of safety depends on the tree size inside the clearing. Simulations 50/40 and 50/50 have values of factor of safety that remain significantly higher than the remaining simulations, although their values sometimes oscillate very close to 1. Although these two configurations have undergone some downhill motion, it is limited to a few 5 centimeters, significantly less than the other cases. These two slopes with large trees are in critical condition because their factors of safety is nearly equal to 1 (< 1.1). Slope motion is limited to a small area near the center of the clearing and to very few cell moves owing to the large tensional resistance of roots that limit downslope movement. Figure 14 shows the distribution of displacement, root and soil bond forces across the slope just before failure (or at the last time step of the run for simulations that did not fail), and displacement at failure for all 6 simulations (one simulation 10 per row). For cases where a landslide occurred (0, 10, 20, and 30-cm), only the clearing area fails except for the 50/30 case where the entire slope fails (see Fig. 14, last column). The clearing with the 30-cm trees pulls down the slope with the stand of 50-cm diameter trees. Lateral root forces in the clearing and across the stand/clearing transition, and downslope tensile forces in the stand are significantly higher for the 30-cm simulation than for any other simulations (see Fig. 14r,s). Despite smaller displacement before failure (Fig. 14p), 30-cm diameter tree roots mobilized more force than the simulation with smaller trees 15 owing to higher root density and sizes, and larger root stiffness. This caused high downslope root forces at the upper edges of the clearing. Also, lateral force in the clearing are higher and extend across the full width of the clearing (Fig. 14s). As a result, unlike simulations with smaller trees inside the clearing, tensile root failure does not occur inside the clearing but outside in the stand, resulting in the collapse of the stand, pulled down by lateral forces originating from the 30-cm diameter trees in the clearing. This is a case where lateral tensile forces plays a crucial role: by extending spatially across a larger area, lateral forces of the 30-cm trees eventually pull roots of 50-cm trees when the clearing fails. This numerical result helps explain the field observations of Rickli and Graf (2009) who found that the mean landslide area was greater for forested slopes than for nonforested slope in the same catchment. This behavior is also illustrated in Fig. 15 which shows downslope soil and root forces as well as across-slope root force along two downslope sections (x = 0 and x = −9 m) for the 6 simulations at the time step just prior failure (for simulations that resulted in a landslide, i.e., 10, 20 and 30-cm diameter trees in clearing) or at the end of the 5 run (40 and 50-cm diameter trees). Results in that figure clearly show that soil bond forces along the slope center (Fig. 15a) for small trees (0 to 30-cm in diameter) reach significantly higher values than for large trees (factor of 5 to 6), eventually reaching the soil maximum compressive strength just before failure. Downslope root-bond force is smaller for these smaller trees owing to smaller densities and smaller root sizes (Fig. 15b). For larger trees (40 and 50 cm), roots take up some of the load on the soil, reducing compressive forces in the soil downslope. The situation is nearly similar along the clearing-stand transition ( Fig. 15d-f) except that the 30-cm diameter trees have the highest downslope root bond force (Fig. 15e, ), and across-slope root force (F x root ), at the time step just before failure and displacement at failure (dfail) for the 6 simulations with different tree size diameters inside the clearing shown in the smaller displacements of soil for the 30-cm trees than for smaller trees (see Fig. 14a,f,k,p), yet significantly larger than for the 40 and 50-cm trees (Fig. 14p,u,z, note the different scale), downslope root force is maximized for the 30-cm trees owing to the combination of displacement and root-diameter sizes that are mobilized. This is also observable on the across-slope root bond force which is highest for that tree size (Fig. 15f). The across-slope root force is nearly zero for the large trees, all the load being handled via the downslope bond forces (Fig. 15c,f). For the smaller trees, the across-slope root force is significant 5 at the clearing-stand transition (Fig. 15f) but small or close to zero at the center of the clearing (Fig. 15c). The 30-cm diameter configuration stands out from the others in having the largest across-slope root bond forces which eventually fail outside the clearing area, entraining the large trees in the stand during the collapse. Figure 15 also helps understand why the 50/10 simulation fails before the 50/0 simulation. This is counterintuitive but is the result of the force balance and the effects of higher root stiffness with increasing root diameter. For the 50/10 simulation, slope-10 parallel (across-slope) root reinforcement is slightly smaller (order 100 N) than for the 50/0 simulation and limited to areas around the large, 50-cm, tree trunks at the edge of the clearing ( Fig. 15e and also compare Fig. 13d to Fig. 13i). Reinforcement values are smaller because displacement is smaller (compare Fig. 13a to Fig. 13f) and root stiffness small. As a result, slightly more load is taken into downslope soil compression in the 50/10 case than in the 50/0 case (again order 100 N, see Fig. 15a) and the 50/10 case reaches maximum soil compression before the 50/0 case, failing first despite a higher root density and a 15 higher potential for root reinforcement. For the larger trees (simulations 50/20 and 50/30), although displacement is less prior to failure (e.g., Fig. 13k,p), the larger root stiffness associated with these larger tree roots produces larger root forces at smaller displacements resulting in less downslope soil compression for the same time and thus delaying the time to failure. Inset shows details at early stage of displacement and the failure of roots across the tension crack at 580 minutes for the 5-m tree spacing.

Tree spacing
Effects of tree spacing on slope stability also yielded some unexpected results. Trees where spaced evenly on the slope using the center of the slope (x = y = 0) as the reference point for a tree. All other trees are located at equal intervals along the x and y axes from this central tree. Figure 16 shows displacement as a function of time at the slope center for five simulations with tree spacing of 3, 5, 7 (2 simulations), and 10 meters. Intuitively, one would expect that increasing tree spacing would 5 decrease root reinforcement away from trees and increase the likelihood of a weak zone to fail. Results, however, show a different behavior. Slopes with trees spaced 3 or 7 meters (no offset, simply called 7-m spacing in Fig. 16) apart were stable but the slope with tree spacing of 5 meters was not. Despite having higher tree density than the 7-m spacing simulation, and thus having higher root density and root reinforcement values, the slope with the 5-m tree spacing failed at t = 1781 min while the 7-m spacing did not fail. Figure 17 shows the slope displacement and tree density (a-e), downslope root force (f-j), and 10 downslope soil force (k-o). Because trees are spaced at regular intervals around y = 0, tree positions for the 5-m spacing are at y = 0, 5, 10, and 15 meters. For the 7-m spacing, trees are positioned at y = 0, 7, and 14 meters. The vertical crack that forms upslope occurs at the smallest root reinforcement location in between two rows of trees. This is at about y = 13 m for the 5 and 10 meter spacing, and at about y = 11 m for the 7-m spacing (see Fig. 17f-j). The 3-m spacing had sufficiently high root reinforcement that a crack did not form. The vertical crack is 2 meters higher up the slope for the 5-m tree spacing than for the 15 7-m spacing. Because the crack is higher upslope, the number of cells along the y axis that move downhill due to loading is larger for the 5-m spacing than for the 7-m spacing. As a result, near the bottom of the hill, compression is significantly higher for the 5-m spacing (Fig. 17l). With increasing load, the 5-m spacing slope reaches its ultimate value of compression and fails while the 7-m spacing never reaches that point. (n) 7 m with offset Figure 17. Effect of tree spacing on hillslope behavior. (a-e) Tree location (black circles) on the hillslope over slope displacement at the last stable time step or last time step for simulations where no landslide occurs (see Fig. 16). To explore the effect of crack location on slope stability, trees in the 7 meter spacing slope were offset 2 meters uphill, so that a vertical crack would form at a higher elevation than without the offset. This simulation is shown with a dashed curve in Fig. 16. Figure 17d,i,n show the hillslope for this simulation and indicate that a crack forms at y = 13, like in the 5 meter simulation, thus resulting in high soil compression forces downslope. In this configuration, the slope eventually fails. High values of soil compression forces that lead to failure (5 m, 7 m with offset, and 10 m spacing) are clearly visible in Fig 17l-o, 5 and contrast with lower soil compression forces in simulations that did not fail (3 m, 7 m without offset, Fig. 17k,m).

Effects of maximum root diameter
SOSlope was used to test the influence of the range of root diameter classes on the stability of a slope. Figure 18 shows the displacement at the center of the slope (x = y = 0) as a function of time for six simulations with different maximum root diameter: 5,7,8,10,20, and 100 mm. The simulations with 20 and 100 mm maximum root size diameter have no landslide and 10 are practically indistinguishable. This is because the number of roots larger than 20 mm is insignificant and contributes little additional strength to the root bundle. The 8 and 10-mm simulations also do not fail and have only slightly larger displacements (6 to 7 cm instead of 5 cm for the 20 and 100-mm simulations). The two simulations with a maximum root diameter class of 5 and 7 mm, however, fail at 1400 and 1500 minutes, respectively. The threshold for stability is thus obtained by including root size up to 8 mm in diameter. Root reinforcement that includes only smaller roots is significantly smaller than if the entire 15 bundle is included. Not including large roots can yield incorrect predictions of slope behavior. Figure 19 shows the downslope and across-slope forces at the center line for several simulations with different maximum root-size diameter. The 5-mm simulation has the smallest amount of downslope root force but the highest across-slope root force and downslope soil compression, explaining why this simulation fails while others (10,20, and 100-mm maximum root diameter) do not. Insufficient root density and lack of large roots compromises the stability of the slope by offering little 20 resistance to loading and declining shear strength of soil. Lateral root forces are small for all cases and has negligible impact here on slope stability (Fig. 19d). 6 Synthesis of force redistributions during triggering of shallow landslides Figure 20 summarizes the typical evolution of forces during landslide initiation of a forested slopes for the case 50/30 described in Section 4.3 (see Figs. 14-16). In this simulation, a clearing is planted with trees 30-cm in diameter while the rest of the slope has trees 50-cm in diameter.
The largest force that contributes to slope stability is soil compression in the area above the landslide toe. There, soil 5 compression increases initially rapidly until it plateaus at about 700 minutes. During this increase, root tension across a growing crack increases and also plateaus. Root compression downslope similarly increases and stops increasing but is significantly smaller than either root tension upslope or soil compression downslope. This time period is defined as Phase 2 of our landslide initiation process which starts when many areas of the slope have a factor of safety that has decreased to 1 (Fig. 20b). Phase 1 of the initiation was the decrease of the factor of safety due to loading and soil weakening without any slope motion. 10 At t = 720 minutes, the roots across the tension crack fail and that tensional resisting force goes to zero. Instantaneously, the slope moves downhill and the force lost by tree roots is taken up by both soil and root compression downslope with the soil taking up most of the increase. This is the beginning of Phase 3. With continued loading, soil compression increases but root compression slowly decreases. Lateral root forces at the edge of the clearing begin to take some of the load to resist downslope movement. Eventually the soil maximum compressive strength is reached and the clearing fails just before 1800 minutes. The time span of the three phases vary with tree size, tree spacing, maximum root diameter, and of course soil and hydrological properties (here fixed for all simulations). Looking back at Figs. 13, 16, and 18, phase 2 can last from several hours to less than one. Sometimes, no crack forms, there is no crack-root failure, and phase 2 and 3 overlap. When the slope has no clearing (as in simulations shown in Figs. 17 and 18), these same three phases exist but lateral forces play no role. Force redistribution and force balance is dominated by soil compression, adjusted by root tension in the upslope area and to a lesser 5 extent root compression downslope. Root forces modify the force balance significantly but soil compression, due to its magnitude, dominates and controls the slope stability and its time to failure. Simulations with smaller soil depth will change this balance: smaller depth will decrease the absolute values of soil compression (see Fig. 3) and tree roots will then support tensile and compressive forces equal or greater to soil compression. In such a situation, roots may be the main factor controlling slope stability.

Conclusions
There are growing evidences that the effects of root reinforcement on slope stability are the results of complex interactions of different factors in which individual contributions are difficult to isolate using classical methods (e.g. infinite slope calcula-tions). The model SOSlope presented here is the final element of a series of related studies aiming to quantitatively upscale the stress-strain behavior of rooted soils under tension, compression, and shearing. In this framework, SOSlope represents the final module where previously investigated aspects of root reinforcements are combined to quantify the macroscopic influences of root reinforcement on slope stability considering spatial heterogeneities of root distribution. The model can produce a systematic analysis of the factors influencing the contribution of root reinforcement on slope stability, yielding a quantitative basis for 5 discussion of root reinforcement mechanisms for slope stabilization and support for the assumptions or simplifications needed to implement such effects in simpler approaches for slope stability calculations (Dorren and Schwarz, 2016). Specifically, simulation results obtained with SOSlope highlight the potential of the model to investigate fundamental questions such the role of forest structure (e.g. tree size, tree spacing), root distribution, and root mechanical properties on the triggering mechanisms of shallow landslides. Based on the results presented here the following general statements can be made:

10
-Maximum root reinforcement under tension and compression does not take place simultaneously; -Root tensile strength is more effective than root compressive strength in preventing or delaying a landslide; -The stabilization effect of roots depends on their spatial distribution: the presence of a"weak zone" leads to behavior similar to bare soils. With little or no root reinforcement, slope failure is more likely and occurs earlier.
-Root reinforcement at the macroscopic scale is dominated by intermediate to coarse roots when present. For the specie 15 considered here and based on available data, roots between 5 and 20 mm contribute the most to root reinforcement.
-Tree positions in the tension zone of a potential landslide influence the stability of the slope. In general, the effect of lateral root reinforcement in tension contributes most to stability along the transition between stable and unstable zones of the hillslope where a crack can form.
These observations indicate that the standard, slope-uniform, constant apparent cohesion approach for rooted soil is often in-20 appropriate, especially for forested slopes where roots contribute significantly to the balance of forces. For example, our model shows that the specific locations of trees on a slope (Fig. 17) is important for predicting slope failure, a conclusion that cannot be reached with the apparent cohesion model. Also, root force distribution on the slope may result a larger landslide for trees with higher root densities (Fig. 14), a result impossible to predict with the apparent cohesion model. Also, root stiffness can modify the time to failure of rooted soils by either increasing or decreasing forces mobilized in roots at different displacement 25 (Fig. 13). Finally, our simulations quantify the importance of considering the heterogeneous distribution of tensional as well as compressional root and soil forces, an element that is entirely missing from traditional infinite slope stability models.
To our knowledge, SOSlope is the first model to implement a new approach that characterizes the force-displacement behavior of rooted soils under both tension and compression. Including this fundamental behavior is key for understanding and modeling shallow landslide triggering. Further work is needed to extend the applicability of standard geotechnical methods 30 (e.g. Schwarz et al., 2015) for the quantification of those soil and root forces.
The SOSlope model can be applied at the hillslope scale to investigate the effect of single factors such as root distribution and root mechanical properties (species specific) on slope stability, and quantification of bio-engineering measures and pro-tective effects of forests. An important application at the hillslope scale is the testing of hypotheses that would support the simplification of calculations in problem-specific applications, e.g. for slope stability model at a regional scale.
The use of the SOSlope model at the catchment scale will be useful for studying the effects of vegetation on slope stability processes in the short and long terms. In the long terms, root strength can vary orders of magnitude (Vergani et al., 2016), and estimation of slope stability and landslide initiation is necessary for an integrated management of mountain catchments for risk 5 reduction and control of sediment balance. In the short term, estimations of safety factors for rooted slopes provide important data for risks assessment in forested mountain catchments. Future work will focus on both these short and long time scales.