The impact of soil-vegetation-atmosphere interaction on a landslide initiation under torrential storms: a case study

. Vegetation has long been used in geotechnical engineering to protect slopes from surface instability phenomena. However, its effects remain difficult to quantify, given the diversity of plants in situ, the variability of their morpho-mechanical characteristics and their impact on soil hydraulic states. For this purpose, the case of a landslide that occurred following torrential rains in Costa Rica was analysed. After the geotechnical characterisation of the soils affected by the movement, the plant species and their main morpho-mechanical characteristics were identified. These characteristics were used to evaluate the changes in the hydromechanical parameters of the soil due to the presence of the roots, and the transpiration rates generated by the plants themselves. In the FE model, a “continuum material” approach was adopted to simulate the vegetated soil numerically, jointly with a failure criterion for partially saturated soils. Using available meteorological data, the evolution of the hydraulic state of the slope in the weeks before and during the storm that caused the landslide was simulated, comparing the cases of vegetated and hypothetically not vegetated slopes. Results validated the observations carried out on the site and confirmed the hydromechanical effects of plants in delaying in time the triggering of the landslide.


Introduction
Global warming and the need for sustainable environmental development represent new challenges for engineers. In the field of geotechnical engineering, very intense isolated rainfall events occurring after long periods of drought, are causing more landslides and erosion events [1][2]. In Costa Rica, heavy rainfalls, and the presence of steep slopes in its mountain system result in many landslides. Among the events that affected the country most significantly in the last decade are Hurricane Tomas in 2010, Hurricane Otto in 2016, Storm Nate in 2017 and Hurricanes Eta and Iota in 2020 [3]. Using vegetation against soil instabilities is a sustainable alternative and may be a possible solution to mitigate these phenomena. Despite being an attractive solution, the impact of vegetation on slope stability has not been studied much in Costa Rica. The research that has been generated so far has focused more on erosion phenomena and from the perspective of agronomy [4,5]. Porras and Laporte (2002) [4] indicate that vegetation has multiple hydraulic effects on the underlying soil. Plant leaves behave as protection in the interception of drops, cushioning the impact on the soil. In addition, transpiration from vegetation increases the level of matric suction in the soil, which is reflected in an indirect increase in the shear strength of the soil, particularly after a dry season. * Corresponding author: alessandro.fraccica@isprambiente.it The mechanical soil reinforcement (increase in shear strength) provided by roots has been extensively studied in geotechnical laboratory investigations, generally as an increase in cohesion [6-9 among others]. Regarding soil hydraulic behaviour, an increase in saturated permeability has been observed in soil with roots compared to soil without vegetation [10][11][12][13]. This increase is associated with the generation of fissures and due to root growth, which create macropores, root-soil interfaces and preferential pathways for water. During the soil drying process, plants produce mucilage, which acts as a surfactant and generates some pore occlusion. In addition, root-soil interfaces are expected to enlarge on drying due to the differentiated shrinkage of roots and soil. All these phenomena induce changes in the water retention and permeability of the soil [13][14]. Compared to the numerous laboratory studies, quantification of vegetation effects at the site scale is more complex and scarcer [15][16][17], given the costs and the number of variables involved.
In view of what was said above, the main objective of the research is to understand how the roots and leaves cover of vegetation affected the stability of a slope for a landslide case in Costa Rica following a storm. This analysis will be based on evaluating the main effects of vegetation and matric suction on the hydromechanical behaviour of the soils within the slope through constitutive laws described in the literature. Finally, the stability of a slope, with and without vegetation, will be simulated using numerical analyses considering the hydraulic boundary conditions, the initial state of the ground, atmospheric factors and plant parameters observed at the study site.
2 Study site characteristics 2.1 Location, rainfall, and atmospheric variables The study site is located in Costa Rica, in the province of San José, in the Aserrí district, on national route 209 at kilometre 22, in the San José -San Ignacio de Acosta direction. The daily precipitation that occurred during storm Nate (October 2017) was recorded from the meteorological stations "Aserrí 1" and the station "Repetidora Abejonal" of the Costa Rican Electricity Institute (Instituto Costarricense de Electricidad) and ranged between 136 mm/d and 217 mm/d. Since the storm was concentrated over 4 hours, an hourly rainfall intensity between 34 mm/h and 54 mm/h occurred. Information from the National Meteorological Institute (Instituto Meteorológico Nacional) of Costa Rica has been used to reconstruct the typical meteorological conditions of the study area [18]. This institute provided the climatic variables necessary for the numerical modelling of this study. Mean values of these variables during September are: wind speed: 2.22 m/s, atmospheric temperature: 17.4 °C, relative humidity: 88.9%, and maximum hourly solar radiation: 1.42 kW/m 2 .

Stratigraphy and geotechnical data
SPT tests were carried out in 2 verticals down to 5 m from the ground [19]. Soil physical properties (grain size distribution curves, natural water contents and Atterberg's limits) and residual strength parameters (torsional ring shear tests) were evaluated on the remoulded material extracted from the SPT sampler [19]. A variable head permeability test in a surface well was carried out in the upper silty soil layer [3]. A seismic prospection was carried out to explore the geometry and the disturbance of the deeper rock layers [20]. The geotechnical model with the different soil layers is presented in Fig. 1, while the geotechnical parameters obtained upon the interpretation of the geophysical information and the geotechnical in-situ/laboratory test results, are summarised in Table 1. As observed in the figure, the water table is located within layer 3 and emerges on the slope surface in correspondence the interface with layer 4. The water retention curves attributed to the soil layers were taken from [13] as they were consistent with those inferred, following the approach proposed by [21], from the grain size distribution curves of the site's materials. The water permeability in partially saturated conditions was evaluated, from the retention curves, following the approach proposed by [22,23]: kw = kw,sat Sr 3 1

Vegetation characteristics and their impacts on soil parameters
The current vegetation species on the slope are predominantly coffee crops combined with Avocado and Sambucus trees. Few/no vegetation was observed at the toe of the slope due to an erosive phenomenon caused by a past debris flow event propagating close to the lower part of the studied slope. Some values of the Leaf Area Index (LAI) were directly inferred on the trees laying on the slope by image analysis through ImageJ [25]. LAI is defined as the sum of the areas of the leaves present in a section of the tree canopy, divided by the area of the same section. Plants of the same species presenting similar values of LAI in the literature were chosen to adopt the root parameters needed for the numerical analyses. These parameters: Root Area Ratio (RAR -ratio indicating the sum of the cross-sectional areas of all the roots laying in a given soil section parallel to the ground surface), Root Length Density (RLD -ratio indicating the sum of the single root lengths within a unit volume of soil) and the average Root Tensile Strength (σt,r -evaluated by tensile tests on single roots with diameters between 1.9 and 2.5 mm). Values of the morphological and mechanical characteristics of the plants considered in this study are summarised in Table 2.   The water saturated permeability variation induced on soil by roots was calculated based on correlations between RLD, kw,sat and macropore void ratio over total void ratio proposed by [13]). In particular, knowing the void ratio e, and the saturated permeability of the bare soil, the macropores void ratio eM with higher connectivity in the bare soil was inferred using equation 3 [13]. Macropores correspond to the void spaces between aggregates. An increase in macropores void ratio was observed by [13] on vegetated soils and depends on the generation of macropores during root growth and a partial closure of the remaining micropores due to complex chemo-hydro-mechanical interactions. Once the new eM was inferred for the vegetated soil (assuming that vegetated and bare soils have the same void ratio at saturation), its permeability was estimated using equation 3.

3
The water saturated permeability estimated in the vegetated soil layers 1 and 3 in Table 1  bare soil case, the vegetated soil water retention curves were taken from [13], consistently with the RLD values observed for the vegetation of this study case. These curves exhibited lower air entry values and lower retention capacity than those evaluated on the same soil in bare conditions due to the increase in macropores observed after root growth in the matrix. The bare and vegetated soil water retention curves were fitted using the model of [31], based on a modified form of the van Genuchten expression [32]: Where s is matric suction, Sr the degree of saturation, C(s) a function used to adjust the higher suction range of the curve. The air entry value is approximately linked to 1/α, the slope of the curve to n and the residual degree of saturation to m. The maximum suction corresponding to the residual degree of saturation is indicated by parameter a. The calibrated values of these parameters are provided in Table 3.

Numerical analyses
The modelling of slope stability under the rainfall event was carried out using GeoStudio 2019 R2 [33]. A comparison of stability analyses between the case of vegetated and non-vegetated slopes was carried out. In the simulations of the vegetated soil, two strata have been considered with thicknesses of 0.30 m and 0.70 m at the top of the slope, affecting layers 1 and 3, as can be seen in Fig. 2. Geotechnical parameters (c'r and kw) of the improved soil with roots, presented in chapter 2.3, were assigned to these layers. The road asphalt lying on the slope's upper part was modelled as a very low-permeability boundary (i.e. with the same geotechnical properties of layer 1 and a saturated permeability of 10 -15 m/s). Below layer 4, an impervious boundary condition was attributed, assuming a lower disturbance of the rock with depth. The boundary conditions controlling soilvegetation-atmosphere interactions are implemented in the hydraulic suite of GeoStudio (SEEP/W) and are based on the Penman-Monteith (FAO) equation [34,35]. These flow boundary conditions are represented by green (on vegetated soil layers) and blue (on nonvegetated layers) arrows in Fig. 2: these conditions only differ in the values of LAI used. Indeed, the atmospheric variables presented in chapter 2.1, a LAI = 4 (average value of the estimates in Table 2) for the vegetated layers and a LAI = 0.5 for the bare soil layers, were used in the equation. The initial groundwater level was positioned, in SEEP/W, as presented in Fig. 1 and by imposing a condition of water pore pressure = 0 on the groundwater surface. The 14 days before the storm (with no rainfalls observed) and the storm itself were modelled in the Finite Element hydraulic suite (SEEP/W) under transient conditions to reproduce the evapotranspiration at the ground surface, the initial hydraulic state of the slope before the storm and the rainfall hourly intensity that triggered the landslide. The stability analysis was carried out at different time intervals to follow the evolution of the factor of safety throughout the rainfall event. The factor of safety was calculated using the limit equilibrium method (SLOPE/W) and considering the approach of [36]. The failure criterion adopted is the Mohr-Coulomb criterion for partially saturated soils based on the following constitutive stress (σ', see e.g. [37]): σ' = (σ-ua)+(ua-uw)•(θ-θr)/(θs-θr) 6 τf = c' + c'r+ [(σn-ua)+(ua-uw)•(θ-θr)/(θsθr)]tanϕ' 7 Where σ and σn are the total stress and the total stress orthogonal to the plane on which the maximum shear stress (τf) lays, ua and uw the air and pore water pressures, θ, θr and θs are the actual, residual and saturation volumetric water contents, estimated on the soil water retention curves presented in Table 3 (for bare and vegetated soil, according to the case analysed). The evolution in volumetric water contents during the rainfall simulations was calculated by the software using Richard's equation [38], assuming the water permeability values defined by equation 1. Finally, c' and c'r are the soil cohesion and the additional cohesion due to root reinforcement, respectively. The rainfall intensity used in the model was 44 mm/h. The duration of the event was prolonged until failure was observed (safety factor FS =1) in the model. Fig. 3 shows the saturation states in the modelled slopes after the 14 preceding days before the storm and without rainfall. Due to the higher LAI value in the vegetated slope, a higher evapo-transpiration rate was computed by the Penman-Monteith equation. This made the shallow vegetated soil reach a drier condition than the bare soil in the same position along the slope. These hydraulic states were the initial conditions upon which the rainfall boundary condition (44 mm/h) was applied. The evolution of the minimum safety factor recorded at different rainfall times is presented in  It is interesting to note that in both slopes analysed, the failure surface resulting with FS = 1 was very similar to that shown in Fig. 5a, and also observed by c' and ϕ' reduction through a finite element approach in Plaxis 2D [39] (Fig. 5b). Moreover, the critical failure surface appears deeper than the vegetated layer and almost not intersecting it, apart from a short part at the slope crest. Hence, mechanical reinforcement exerted by roots on soil remained almost un-exploited, while the only factor that delayed the failure time was the drier hydraulic state, improving the shear strength of soil via the product of suction and effective volumetric water content (or, alternatively, degree of saturation) formalised in equations 6 and 7.

Results and discussion
The modelled failure surface is in good agreement with the evidence of the landslide observed on-site and presented in Fig. 6. Some settlements were observed on the berm above the crest of the slope (point A in the figure), and tensile cracks / shear deformations along the road in the upper part of the slope (point B).

Conclusions
The results of the simulations showed that, due to evapotranspiration at the beginning of the storm, the slope with vegetation had a lower degree of saturation and higher matric suction than the slope without vegetation, which translated into a greater soil shear strength, in part increased by the tensile strength of the roots themselves. Due to this difference in the initial hydraulic state, the slope with vegetation reached instability a few hours later than the slope without vegetation, despite an increase in soil water permeability due to root penetration/growth in the medium. The simulated time evolution and the kinematic of the landslide were in good agreement with local observations, despite the availability of few geotechnical data, little information on the root parameters of the plants and the use of simple constitutive equations to infer the changes in soil shear strength and permeability induced by vegetation. This study shed light on the importance of considering both hydraulic and mechanical effects of roots in the soil, as the matric suction induced by plant transpiration plays a key role (in this case, the main role), together with the tensile strength of roots, in delaying slope instability during heavy rainfall. In this way, more attention should be paid to monitoring such kind of variables within slopes to get an accurate assessment of the actual factor of safety and the probabilities of landslide activation. Fig. 6. Evidence of the landslide at the crest of the slope (point A in Fig. 5) and on the road surface (point B in Fig. 5).