The role of plants in the prevention of soil-slip: the G-SLIP model and its application on territorial scale through G-XSLIP platform

This paper discusses the role of plants in the prevention of shallow landslides induced by rain (soil slips); these phenomena, related to “hydrogeological instability,” are among the most feared because their evolutionary processes can cause huge damages and losses of human lives when interacting with anthropized areas and infrastructures. The paper first highlights how the plants interact with the soil; then introduces the G-SLIP (Green – Shallow Landslides Instability Prediction) model, i.e., the simplified physically-based SLIP model, modified to predict soil slips at punctual and large scale taking into account the vegetation effects. The G-SLIP model is thus applied to a case study of the Parma Apennines (Northern Italy) by using the G-XSLIP platform. In this area, during the intense events of rain between the 4th and 5th of April 2013, numerous landslides occurred, provoking huge damages to structures and infrastructures, and consequent economic losses. The stability analyses carried out with G-XSLIP demonstrate that the presence of vegetation in the study area led to a significant reduction in the triggering of shallow landslides. Finally, an attempt at soil slip mitigation through naturalistic techniques (planting of specific vegetation) is presented.


Introduction
The relationship between plants and soil is widely studied in natural sciences and agronomy; much less in soil mechanics, even if it is well known that the plants improve the strength characteristics of the soil and reduce rainfall infiltration, preventing landslides (O' Loughlin and Pearce 1976;Megahan et al. 1978;Sidle and Wu 1999;Dhakal and Sidle 2003;Imaizumi et al. 2008;Imaizumi and Sidle 2012). The interest in assessing the plant contribution to the soil has derived since the 1960s, from the evidence of a considerable increase in landslides following the deforestation of some areas of Alaska (Bishop and Stevens 1964;Burroughs and Thomas 1977;Bordoloi and Ng 2020;Montrasio and Gatto 2020). Moreover, several agronomic studies have focused on the role of plants in rainfall interception, i.e. the "canopy effect" (Carlyle-Moses and Price 1999; Hemmati et al. 2012;Mużylo et al. 2012;Corti et al. 2019;Lin et al. 2020;Dong et al. 2020), and the root water uptake (Gardner 1964;Feddes et al. 1976;Prasad 1988;Yu et al. 2007). A detailed discussion on this topic has been recently presented by Montrasio and Gatto (2020). In any case, as the roots rarely reach excessive depths, their contribution to the soil strength and to the slope stability is meaningful for phenomena involving modest soil thicknesses, such as erosion or shallow landslides; the latter will be the focus of this paper.
Shallow landslides (or soil slips) involve thin debris layers (topsoil), less than two meters thick, triggered by intense and/or prolonged rainfalls that penetrate the topsoil and cause instability; once triggered, they can evolve into flows of debris or mud, characterised by sustained speeds (up to 9 m/s) and a highly destructive power. Soil slips are widespread not only in Italy (Montrasio and Valentino 2007;Montrasio et al. 2009;Martelloni et al. 2012), but also worldwide (Ko Ko et al. 2004;Staley et al. 2012;Chung et al. 2017;Cogan et al. 2018;Dang et al. 2019;Li et al. 2020;Lee et al. 2021;Kim et al. 2021;Shen et al. 2021). The rapid triggering and their spatial-temporal unpredictability make soil slips dangerous and high-risk phenomena, especially when they occur near anthropized areas (Crozier 2010;Kirschbaum et al. 2012Kirschbaum et al. , 2015Dowling and Santi 2014;Froude and Petley 2018). Therefore, the prediction of soil slip triggering and the mitigation of the associated risks are of particular interest.
Forecast models of rainfall-induced shallow landslides can be classified into two broad categories: "black box" and "physicallybased" models. The "black box" models relate the triggering to the precipitation expected or measured, through relationships that are mainly empirical-probabilistic (Caine 1980;Cannon and Ellen 1983;Govi et al. 1985;Wieczorek 1987;Crosta 1990Crosta , 1994Crosta , 1998Guzzetti et al. 1992Guzzetti et al. , 2007Guzzetti et al. , 2008Ceriani et al. 1994;Crosta and Frattini 2001;Marchi et al. 2002;Giannecchini 2006;Segoni et al. 2018); thanks to their simplicity, these models are widely applied but their predictive capability is often unsatisfactory. From the analysis of the phenomena that really occurred, it has been well demonstrated that the triggering of soil slips depends not only on the amount of precipitation but also on the geometry, morphology, and mechanical and hydraulic characteristics of the soil. Physically-based models introduce these dependencies and, using the classical mechanics of soil approaches, are able to model the real phenomena taking into account geometry and mechanical and hydraulic parameters of soil, in addition to rainfall data. Among these, there is the SLIP model. SLIP (Shallow Landslides Instability Prediction) is a physicallybased model that defines the evolution with time of the safety factor of slopes potentially at risk of soil slip, using the equilibrium limit method, associated with simplified hypotheses on both the rainfall infiltration phenomenon and the mechanical behaviour of the unsaturated soil (Montrasio 2000). The SLIP simplicity makes it suitable for application on a territorial scale. The model, applied to different territories, showed excellent predictive skills (Montrasio andValentino 2007, 2008;Montrasio et al. 2009Montrasio et al. , 2011Montrasio et al. , 2012Montrasio et al. , 2014Montrasio et al. , 2015Schilirò et al. 2016;Montrasio and Schilirò 2018); specifically, its use on a territorial scale was possible thanks to its implementation into the multi-risk platform owned by the Italian Civil Protection. Recently, Gatto and Montrasio (2023) developed a proper platform based on MATLAB, named X-SLIP, where they implemented SLIP. Their results underlined the necessity to properly model the vegetation effects and consider them in SLIP formulation and X-SLIP algorithm. This gap is overcome by this paper, aimed at investigating how vegetation reduces the triggering of shallow landslides.
After an overview of SLIP and X-SLIP (the "Overview of SLIP model and X-SLIP platform" section), plant contributions are examined and introduced into the G-SLIP model and G-XSLIP platform ("G" standing for green), respectively SLIP and X-SLIP update (the "Vegetation contributions: the G-SLIP analyses" section). Finally, the application of the modified G-SLIP model to a study area located in the Parma Apennines is shown (the "Application of the G-XSLIP to the Parma Apennines study area" section). In this area extensive shallow landslides occurred from the 4 th of April to the 5 th of April 2013 after intense rainfalls; the analyses are performed by considering the cases of "no vegetation" and "with vegetation" (Vegetation Map provided by ISPRA, i.e., the Italian national network for the environmental protection). G-SLIP can be the starting point for the selection of specific vegetation for the active protection of areas that are vulnerable to rainfall-induced shallow landslides.

Overview of SLIP model and X-SLIP platform
SLIP model SLIP is a simplified physically-based model aimed at analysing the triggering of rainfall-induced shallow landslides. It was introduced by Montrasio (2000) and has been validated and applied at any scale, from the lab scale (Montrasio and Valentino 2007;Montrasio et al. 2015) to the in situ scale: single slope (Montrasio et al. 2009; and territorial areas (regional and national) (Montrasio et al. 2011(Montrasio et al. , 2015Schilirò et al. 2016;Montrasio and Schilirò 2018). The SLIP model could be defined as a mechanical-hydraulic-dynamic simplified model that couples the loss of equilibrium with the variation of the soil's degree of saturation due to the rainfall infiltration, taking into account the infiltrated rain varying with time.
The model is based on the "infinite slope" scheme applied to the topsoil, whose thickness H ranges around 1-1.5 m (Godt et al. 2008;Montrasio et al. 2011) and represents the depth at which the soil permeability radically changes: above, due to the widespread fracturing, a great heterogeneity is present and facilitates the water infiltration; below, the soil is no longer affected by fracturing and the permeability is much lower. The safety factor F S is defined for a soil characterised by the Mohr-Coulomb criterion T = � tan � +c � tot as: is the slope angle, φ is the soil friction angle, γ is the soil's unit weight, and c � tot is the cohesion of saturated/unsaturated soil. For the unsaturated soil, c � tot is the sum of the effective cohesion c' and the apparent cohesion c (Montrasio 2000): The evaluation of c follows the schematisation of the physical phenomena represented in Fig. 1.
At time t j , the soil (degree of saturation S 0j ) involved in the generic rainfall event (cumulative rainfall h j ) is subjected to the rainfall infiltration that can be represented by a single saturated layer of m j H thickness ( Fig. 1b): The temporal evolution of the saturated area related to a single rain episode is expressed as: It depends on the velocity of desaturation k t , which takes into account permeability and fracturing of topsoil, substrate permeability, and external conditions (affecting the evapotranspiration). At generic time t, the total m(t) is the result of the previous rainfall, and it is evaluated from the superimposition of all generic dampened m j as: With respect to the shear strength model, the soil results in being made up of two different layers: a saturated one (thickness m(t)) and a partially saturated one (thickness 1-m(t)); by assuming no residual volumetric water content and a sort of "hydraulic elasticity" (Montrasio (2000)), the shear strength of the latter is defined by c * : Fig. 1 Schematisation of the effects of rainfall infiltration accounted for by the SLIP model. a Creation of "saturation bubbles" associated with the generic jth rainfall event; b simplification of bubbles with a saturated layer where A is an experimentally derived parameter, dependent on soil particle size (increasing with the decrease of the soil's particle size) and λ is a modelling parameter. The total cohesive strength c of the slice represented in Fig. 1b is derived by homogenising the strength characteristics: The validity of Eq. (7) has been confirmed by experimental analyses carried out on stratified soils (Montrasio and Schilirò 2018).
The evolution of safety factor F S with time is obtained by combining Eq. (1) and Eqs. (5-7) and depends on geometric factors (slope angle β, thickness of topsoil H), soil and water state parameters (specific gravity G S , porosity n, saturation degree S r , water unit weight γ w ), shear strength parameters of the soil (friction angle φ'; effective cohesion c'), the modelling parameters A, λ, α, and the "slope drainage capacity" k t , and rainfall (pluviometric diagram h j ); a detailed description of the SLIP model is presented by Montrasio (2000), Montrasio and Valentino (2008), Montrasio et al. (2016). Montrasio et al. (2014) showed the implementation of the SLIP model in the platform DEWETRA, owned by the Italian Civil (6) c * = AS r (1−S r ) (7) c =c * (1 − m) Protection, capable of performing stability analyses on a territorial scale and its application at many places located in different Italian regions. Recently, the SLIP model has been implemented in a free and easy-to-use MATLAB-based platform, named X-SLIP, which allows spatial-temporal stability analyses to be performed by elaborating territorial data through the MATLAB Mapping Toolbox and its built-in functions, together with specific implemented scripts (Gatto and Montrasio 2023). The stability is assessed by evaluating Eq. (1) for different events (temporal analysis) in a study area discretised with a regularly spaced grid (spatial analysis). The X-SLIP platform runs three computational phases: (i) preprocessing, (ii) processing, and (iii) post-processing; in each phase, the spatial distribution of the parameters, as well as the results, are managed through a matrix approach. A schematic representation of the X-SLIP functioning is illustrated in Fig. 2.

X-SLIP platform
The pre-processing phase is finalised to build the parameter matrices by elaboration on some raw data: the Digital Terrain Model (DTM) and the Lithological Map for the study area, and the rainfall recordings in gauging stations adjacent to the study area. This phase is divided into three subphases: "morphology", "soil" and "rainfall". During the "morphology subphase", the matrix of slope angles β is built by applying a gradient function to the altimetry matrix extracted from the DTM. The matrices of soil parameters (φ', c', n, k t , A, λ, α) are obtained in the "soil subphase", through an appropriate implemented script that performs a "single-point parameter assignment": a point-in-polygon problem is solved to define which lithology polygon contain each point of the reference grid; a soil parameter set (deriving from experimental characterisation or back-analysis) is then associated with each lithology polygon and parameters are distributed to the grid accordingly. Finally, the rain matrix h (one for each rainfall event affecting the stability) is constructed in the "rainfall subphase", through the scatteredInterpolation function, which interpolates recorded rain data with a "natural neighbour" approach and associates a rainfall value to each grid point (Sibson 1981;Amidror 2002).
In the processing phase, the rainfall and soil parameters are combined to evaluate first the thickness of the saturated layer m, through Eq. (5), and then the safety factor F S . The post-processing phase finally allows the visualisation of the susceptibility distribution in the study area (susceptibility maps), the temporal variability of results in target points, and the evaluation of the prediction quality if past soil slips have been detected and their position is known. The validation of the platform was obtained by comparing the results of X-SLIP to the ones of DEWETRA applied to several areas of Emilia Romagna (Gatto and Montrasio 2023).

Role of plants in slope stability
The vegetation contributes to the stability of a slope from both a mechanical and a hydraulic point of view, through roots and foliage respectively. Both contributions are plant-specific; i.e., they depend on the vegetation species.
The contribution of roots to soil strength is related to the root's intrinsic characteristics (size, tensile strength, stiffness) and it is taken into account by adding the extra term c r to the Mohr-Coulomb criterion: c r is defined as root cohesion and it has been linked to the soil failure mechanism since Wu (1976). A simplified expression of c r can be represented by Eq. (9): where k' takes into account the geometry of the root (inclination) and the soil's friction angle, while t r is associated with the rootbreaking mechanisms (tensile, pulling out, or elongation) and the dimensions in relation to the soil's volume. Wu (1976) originally identified the resistance t r as "tensile", while Waldron (1977) discussed the possibility of three root-breaking mechanisms, tensile, pulling out, or elongation; in all these cases, the resistance also depends on the diameter of the root. c r can be measured directly, by means of specific tests (Endo and Tsuruta 1969;Waldron 1977;Waldron and Dakessian 1981;Ziemer 1981;Waldron et al. 1983;Abe and Iawamoto 1986), or derived by back analysis or empirical correlations (Endo and Tsuruta 1969;Swanston 1970;O' Loughlin 1974;Gray and Megahan 1981;Sidle and Swanston 1982;O' Loughlin and Ziemer 1982;Buchanan and Savigny 1990;Bischetti et al. 2002).
Over the years, other theoretical models of increasing complexity have also been proposed for c r evaluation (Waldron 1977;Wu et al. 1979;Waldron and Dakessian 1981;Ekanayake and Philipps 1999;Operstein and Frydman 2000;Pollen and Simon 2005). The (8) τ = � tg � +c � +c r (9) c r =k � t r concept of c r has been extended to complex root systems, as well as to different types of plants (Waldron 1977;Gray and Ohashi 1983;Watson and O' Loughlin 1990;Sainju and Good 1993;Bischetti et al. 2005;Dupuy et al. 2005;Fan and Chen 2010;Cecconi et al. 2012); the Wu and Waldron model for root systems, in addition to Eq. (9), introduces the density of roots (varying with depth). One of the most common expressions for a system of multiple roots is: t r,i is the resistance of the i-th root, while the RAR is the root area ratio and takes into account the density of roots, varying with depth; both t r and RAR depend on the species. RAR is evaluated by means of several techniques, including the root-wall technique based on image analysis (Burke and Raynal 1994;Vinceti et al. 1998;Schmid and Kadza 2001;Mattia et al. 2005).
From the hydraulic point of view, the vegetation provides the rainfall interception by the "vegetation canopy" and the "forest floor"; both of them reduce the amount of water infiltrating into the topsoil. Due to the canopy effect, only the throughfall (a rate of precipitation) infiltrates the soil; this effect depends on the size of the foliage and the type of leaves and varies seasonally, being prevalent in spring/summer. In autumn/winter the leaves of the trees form the "forest floor," i.e., a bed of leaves that also contributes to intercepting and reducing the infiltrating precipitation.

G-SLIP model to take into account the vegetation effects
The SLIP model is modified in the G-SLIP model (Green-SLIP model), updated to consider the double contribution provided by vegetation to slope stability; a first one represented by the increase in soil strength, due to the root systems, and a second due to foliage, which reduces the rainfall infiltration and ameliorates indirectly the shear strength, affecting the partial saturation of the soil.
The first contribution is represented by the root cohesion c r , which is added to the soil cohesion ( c ′ and c of Eq. (2)). The total cohesion c � tot is therefore expressed as: Depending on the root system, c r varies with depth (see Eq. 10); however, a single value is used for Eq. (11) and it is assumed as the average evaluated in the soil layer thickness H.
In terms of infiltration, a reduction coefficient 1-β* (less than unity) is introduced in the m(t) expression, which becomes: The vegetation contributions introduced in the G-SLIP model are summarised in Fig. 3.

G-XSLIP platform: implementation of the "vegetation subphase" in X-SLIP
Following the modification of the SLIP model described in the previous subsection, the X-SLIP platform is updated and renamed as G-XSLIP (Green-XSLIP). G-XSLIP allows to include not only information on morphology, lithology, and rainfall (as in X-SLIP) but also about vegetation through specific maps provided by public entities that contain information on plants, their type, and extent. Vegetation maps are converted into matrices of parameters c r and β* through a simplified methodology based on the Ray Cast (RC) algorithm, already described in Gatto and Montrasio (2023) to consider the spatial variability of soil parameters. In this case, the procedure is vegetation-based: vegetation parameters are assigned to each pixel based on the vegetation type in which the pixel is located. For each vegetation type, a set of correlated parameters is considered, so that they are not spatially independent from each other. The input of the vegetation map and of the user-defined parameters associated with the plant species are performed in a specific "vegetation subphase". Other than the map-based assignment, it is also possible to introduce uniform vegetation parameters. The features of G-XSLIP are summarised in Fig. 4.

Application of the G-XSLIP to the Parma Apennines study area
This section describes the application of the G-XSLIP platform to a study area located in the Parma Apennines (Emilia Romagna Region, Northern Italy), which includes the municipalities of Neviano degli Arduini, Corniglio, Tizzano Val Parma, and Palanzano (Fig. 5a); the total extension of the area is 516 km 2 . In the following, we analyse how the vegetation present in the area affected the triggering of rainfall-induced shallow landslides, which occurred in large numbers between 4 and 5 April 2013, following a period of intense rainfall. The locations of some of the soil slips detected by Montrasio et al. (2015) are shown in Fig. 5b.
The SLIP model has already been applied to the area during a cooperative activity between the University of Parma, the Emilia-Romagna Region, and the Civil Protection Department, the latter allowing the use of the model in their territorial platform DEWETRA (Terrone 2014;Montrasio et al. 2015). The results of the elaborations are shown in Fig. 6; they refer to the analyses conducted taking into account the vegetation contribution only for the runoff and the canopy effects, in a first parametric approach consisting of two cases: i) uniform 1-β* (equal to 70%), assigned to the whole study area (Fig. 6a); ii) specific 1-β* values for forested areas, grass, and bare soil, respectively 60%, 70%, and 80% (Fig. 6b). The results demonstrated the benefits of vegetative coverage, even if related only to infiltration effects. The current study considers both the mechanical and hydraulic contributions of plants, as described in the ""G-XSLIP platform: implementation of the "Vegetation subphase" in X-SLIP" section, with the plant-specific approach implemented in G-XSLIP ("G-XSLIP platform: implementation of the "vegetation subphase" in X-SLIP").

G-XSLIP evaluation of the vegetation-independent parameters for the study area
G-XSLIP analysis is performed in the time window 4-5 April 2013 (109 analysed stability events, one for each hour); the reference grid for the spatial analysis comes from the intersection of the Digital Terrain Model (DTM) and the study area polygon, shown in Fig. 5a. The raw data required by G-XSLIP in the pre-processing phase are the DTM, the lithology map and the rainfall recordings at specific rain gauging stations.
The DTM is provided by the Emilia Romagna region with a 5-m resolution; 19 DTMs, each extension of 6280 m × 7320 m, are required to cover the study area. Since there are no clear criteria for choosing the most suitable grid scale to evaluate susceptibility maps (Iwahashi et al. 2012;, the DTM resolution is changed to 20 m during the G-XSLIP morphology subphase to increase the computational efficiency with an extension compatible with the soil slips occurred in the area. The elevation map obtained from the 20-m resolution is shown in Fig. 7a, while Fig. 7b illustrates the 3D representation of the study area. As for the DTM, the lithology map is also provided by the Emilia Romagna region database; it is essential to assign the soil parameters to each grid point. Figure 7c shows the spatial distribution of the lithology polygons in the study area, belonging to nine lithological classes: massive stone rocks (A), layered stone rocks (As), Rocks made up of alternation with stone dominant Fig. 5 The study area. a Geographical framework; b Detail of the municipalities included with the soil slips detected on April 5, 2013, at 12:00 after intense rainfall (Terrone 2014) Fig. 6 Results of previous elaborations obtained through the SLIP model, on the study area after Terrone (2014) and Montrasio et al. (2015). Analyses were performed with a uniform vegetation parameters; b mean plant-specific parameters levels (Bl), Rocks made up of alternation of stone and pelitic levels (Blp), Rocks made up of alternation with pelitic dominant levels (Bp), Conglomerates and clast-supported breccias (Cc), Marlstone (Dm), Claystone breccias (Dol) and scaly clays (Dsc). Montrasio and Valentino (2008) presented a simplified methodology for soil parameter assignment based on the clustering of lithologies in soil types, each of which is characterised by averaged physical, mechanical, and SLIP modelling parameters. The lithology-soil types association, with the set of soil parameters, is reported in Table 1; these values are at the basis of   the single-point assignment procedure described in the "X-SLIP platform" section followed by X-SLIP to build the matrices of soil parameters. The rainfall data related to the period March-April 2013 were supplied by Dexter Arpae; specifically, the recordings of 12 gauging stations included in the study area, whose positions are shown in Fig. 7d, are considered.

A Massive stone rocks Excluded
After having imported each raw data, G-XSLIP performs the procedures summarised in Fig. 2, to obtain the matrices of parameters required by the SLIP model, involved in the F S evaluation through Eq. (1). The spatial distributions of some parameters are shown in Fig. 8.

Assessment of vegetation parameters
The new "vegetation subphase" implemented in G-XSLIP allows us to assign the vegetation parameters to each grid point of a study area; in particular, the single-point assignment procedure discussed in "G-XSLIP platform: implementation of the "vegetation subphase" in X-SLIP" is used, where the raw data is a vegetation map.
For this study area, the vegetation map provided by ISPRA reveals the presence of forests and herbaceous species, with seven main vegetation types, whose distribution is shown in Fig. 9. Forests are made up of Fagus sylvatica (19.79%), Quercus cerris (14.28%), Quercus pubescens (13.90%), Ostrya carpinifolia (5.82%), Castanea sativa (2.28%), and allochthonous conifers with prevalence of Pinus nigra (2.26%), while the herbaceous species consist of extensive crops (18.31%), scything meadows (8.78%) and intensive crops (3.08%), which are mainly made up of Barley, Wheat, Pea, and Alfalfa in the Emilia Romagna region; the percentages reported in brackets are related to the whole study area. According to the vegetation type, specific root cohesion c r and throughfall 1-β* are assigned, with values obtained from previous studies. Note that the root cohesion is introduced only to the forested area because the rooting apparatus of the herbaceous species is less extensive and their contribution to the soil reinforcement should be net of the harrowing effects (Ziemer 1981;Comino et al. 2010).
For the c r of forested areas, the reference studies are Burylo et al. (2011), for the Quercus pubescens and the Pinus nigra; Cazzuffi et al.  Fig. 10. In all these studies, a c r variation with depth was derived by applying the Wu and Waldron model, taking into account the root system architecture. Since the SLIP stability analysis is performed by considering a soil slice of thickness H (the "G-XSLIP platform: implementation of the "vegetation subphase" in X-SLIP" section), equal to 1.2 m in our case based on experimental observations, an average c r value is introduced. When the root length is smaller than 1.2 m (e.g., for Quercus pubescens), a thickness weighted average is evaluated, considering the mean of non-null values and zero; in the other case, the arithmetic mean is computed. The c r -mean values are reported in Fig. 10.
The quantification of β * is a common subject of agronomic studies, having the main interest in irrigation; its values change depending on the foliar apparatus, according to size, age, and species, as well as the vegetative state of the plants (Aston 1979;Fleischbein et al. 2005). For the species of the study area, the β * reference values used in the analysis come from the studies of Haynes (1940) and Butler and Huband (1985), for the herbaceous species; Staelens et al. (2008) Lin et al. (2020) for Castanea sativa, and Dong et al. (2020), for Pinus nigra. Table 2 summarises the vegetation parameters adopted for each vegetation type, in the analyses performed in the presence of vegetation. Figure 11a and b illustrates the distribution of the vegetation parameters in the analysis with no vegetation, while Fig. 11c and d shows the parameter distribution based on the vegetation map. A runoff coefficient computed as the cosine of β (slope angle) is assigned as 1-β * to the points external from vegetation polygons, as suggested by Lei et al. (2020). Figure 12 illustrates the maps of unstable points (F S < 1) obtained by applying the SLIP model over the period 4-5 April 2013, in cases "no vegetation" and "with vegetation". The results refer to three reference events: (i) 4th April at 13:00, corresponding to when the rain cumulate began to grow after a period of stall ( Fig. 12a and b); (ii) 5th April at 00:00, at the rainfall peak occurrence (Fig. 12c and d); (iii) 5th April at 12:00 when the rain cumulate reached the maximum value ( Fig. 12e and f).

Susceptibility in the absence and presence of vegetation
The susceptibility of the area is observed to be reduced thanks to the plant contributions, which increased the soil strength, directly through c r (according to the map shown in Fig. 11d), and indirectly thanks to 1-β * , which decreases the thickness of the saturated layer m (see Eq. 12), determining an increase in c (see Eq. 7); as an example, in Fig. 13 the spatial distribution of m is shown, relative to the event of 5 April at 12:00 (Fig. 13a, in the absence of vegetation; Fig. 13b, with vegetation). Figure 14 reports the number of unstable points obtained for the three reference events and in the two analysis cases; the hydraulic and/ or mechanical vegetation contributions lead to a maximum reduction of the unstable points equal to 45% for the event of 5th April at 12:00.
To assess the performance and prediction quality of the SLIP model in the two cases aforementioned (with and without vegetation), the receiver operating characteristic curves (ROC) (Metz   not catch landslides occurred. These values allow to define the true positive rate (TPR) and the true negative rate (TNR) as: ROC curves are defined on the TNR-TPR plane, varying the critical safety factor F S,critical that discriminates stable from unstable points; F S,critical values between 1 and 10 are considered, with steps of 0.5.
As reported by Montrasio et al. (2011), a polygon surrounding each of the surveyed points has to be considered, before comparing with the outputs of the SLIP model; the area of each polygon is 40 × 40 m 2 , consistent with the typical extension of the soil slip. When at least one point of the model with F S ≤ F S,critical is contained within the surveyed-landslide polygon, a TP is counted; if all points inside the polygon have F S > F S,critical , the result of the analysis is considered FN. FPS and TNS correspond instead to the number of points outside the polygons of the landslides surveyed, respectively having F S ≤ F S,critical and F S > F S,critical. The forecasting performance of the model is then defined on the basis of the value of area under curves (AUC ), calculated integrating the area under the ROC curve; according to Hosmer and Lemeshow (2000) AUC = 70-80% means acceptable, while AUC = 80-90% excellent performance. Figure 15 shows the ROC curves for the two cases (with and without vegetation) created using the model results for April 5 at 12:00; a better accuracy is shown when vegetation is considered. This meant that the case "with vegetation" is more representative of reality and confirms that during the analysed event, the vegetation helped to mitigate the triggering of some soil slips.

Stability improvement of susceptible areas through naturalistic techniques
Besides the vegetation that is already present on slopes, planting specific types may be a solution for active protection and risk mitigation of susceptible areas (Wang et al. 2019). Since planting itself requires agronomic considerations, such as conditions for root growth, time for establishing root reinforcement, etc., the latter is neglected, and the stability improvement provided by post-planted vegetation in steady-state conditions is evaluated. Our analysis is focused on areas with F S < 1 in Fig. 12f, i.e., the ones related to 5 April 13:00 after introducing the vegetation included in the ISPRA map ( Fig. 9), characterised by 23,045 unstable points. The vegetation polygons containing such points are reported in Fig. 16: more than 60% of these points are occupied by Quercus pubescens, while in the remaining points there are herbaceous species (Alyssoides utriculata and Saxifraga paniculata) or no coverage at all.
The main observation is that the reinforcement considered in our analysis for Quercus pubescens (derived by Burylo et al. 2011) may be not representative of the real reinforcement, being also very low compared with the one provided by Quercus cerris. Further measurements should be done to support the adopted value. Moreover, removing a tree and replacing it with species providing greater improvement is not a rational choice; we propose therefore to intervene in the remaining types. Herbaceous species contribute to slope stability mainly through rainfall interception; the use of vegetation widespread in the area (native species) giving also root reinforcement could be a successful choice.
Even if not present in the municipalities of our study area, neighboring areas show the presence of vineyards and olive groves: these vegetative species are thus considered in our analyses. Their scientific names are Vitis vinifera and Olea europeae; Moreover, Gómez et al. (2001) quantified the rainfall interception by olive trees equal to 20% for high-density plants. For olive trees, we perform three analyses, one with the contributions together (root and foliage) and two considering them separately. In such a way, it is possible to discuss the different effects. Figure 17 summarises the number of unstable points in all the new cases, compared with the original situation.
In all cases, vegetation is observed to improve slope stability decreasing the number of unstable points: from about 23,000 to 8216 (64% reduction) in the worst situation, (i.e. olive trees modelled through the provided rainfall interception), and to almost 1800 in the remaining cases (92% reduction). It is interesting to remark that when introducing vegetation through only root reinforcement, stability is clearly better for greater c r (1736 points with V. vitifera and 1712 points with O. europeae). By modelling both the contributions provided by olives (c r and β * ), the number of unstable points is larger (1831). This may be explained by the adopted interception rate which gives throughfall equal to 80%: such a value corresponds to a slope of about 36° following the relationship introduced by Lei et al. (2020), adopted when vegetation is neglected or missing (see the "Assessment of vegetation parameters" section). For greater slopes angles, the assumed throughfall without vegetation is smaller, causing a minor reduction in the soil's apparent cohesion due to rainfall infiltration (Eq. 12). To focus on this aspect, Fig. 18 shows the F S values recorded from 4 April 12:00 to 5 April 12:00 in a point which is among the unstable ones of Fig. 12f and where the slope angle is almost 50° (1-β * = cos(50°) = 0.64).
It can be observed that stability can be guaranteed with the root reinforcement derived by planting V. vinifera and O. europeae. For O. europeae with no root reinforcement (orange line), the safety factor is lower than the Initial configuration (black line) because of the assumption of a larger rainfall amount infiltrating into the topsoil. For the same reason, even the case with both vegetation contributions is worse than the one with only root reinforcement provided by O. europeae (in the latter case, 1-β * is equal to 0.64, i.e., the value without vegetation). It is worth pointing out that further consideration must be done regarding the possibility to plant the mentioned species in such steep areas and how the vegetation throughfall should be modified to consider the slope angle.

Conclusion
This study has analysed how vegetation improves soil strength, with a specific focus on the role of plants in the mitigation of rainfall-induced shallow landslides. The main contributions, identified in the root reinforcement and rainfall interception due to foliage, have been introduced into the G-SLIP model and G-XSLIP platform. The application of G-XSLIP to a study area affected by soil slips during intense rainfall that occurred between the 4th and 5th of April 2013 has demonstrated that the number of phenomena was reduced by the vegetation covering the area. Moreover, it has been discussed how the planting of specific species may improve the stability of susceptible areas. The goodness of the platform makes it suitable for use in deriving the vegetation parameters that would guarantee the stability of susceptible areas during rainfall. Future studies will consider the use of the updated G-XSLIP to plan mitigation measures in high-risk areas by a selection of proper vegetation to be planted in non-vegetated areas or providing a better reinforcement (if necessary) in areas already vegetated.

Funding
Open access funding provided by Università degli Studi di Brescia within the CRUI-CARE Agreement.

Data availability
Data available on request from the authors.

Conflict of interest The authors declare no competing interests.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/. Fig. 18 Trend of safety factor F S in a point that is unstable in the initial configuration, considering the vegetation effects derived from planting olive trees (O. europeae) or grapevines (V. vinifera)