Modelling shallow landslide susceptibility by means of a subsurface flow path connectivity index and estimates of soil depth spatial distribution

Topographic index-based hydrological models have gained wide use to describe the hydrological control on the triggering of rainfall-induced shallow landslides at the catchment scale. A common assumption in these models is that a spatially continuous water table occurs simultaneously across the catchment. However, during a rainfall event isolated patches of subsurface saturation form above an impeding layer and their hydrological connectivity is a necessary condition for lateral flow initiation at a point on the hillslope. Here, a new hydrological model is presented, which allows us to account for the concept of hydrological connectivity while keeping the simplicity of the topographic index approach. A dynamic topographic index is used to describe the transient lateral flow that is established at a hillslope element when the rainfall amount exceeds a threshold value allowing for (a) development of a perched water table above an impeding layer, and (b) hydrological connectivity between the hillslope element and its own upslope contributing area. A spatially variable soil depth is the main control of hydrological connectivity in the model. The hydrological model is coupled with the infinite slope stability model and with a scaling model for the rainfall frequency–duration relationship to determine the return period of the critical rainfall needed to cause instability on three catchments located in the Italian Alps, where a survey of soil depth spatial distribution is available. The model is compared with a quasi-dynamic model in which the dynamic nature of the hydrological connectivity is neglected. The results show a better performance of the new model in predicting observed shallow landslides, implying that soil depth spatial variability and connectivity bear a significant control on shallow landsliding.


Introduction
Effective management of the hazard associated with shallow landsliding requires information on both the location of potentially unstable hillslopes and the conditions that cause slope instability.The need for spatial assessment of landslide hazard, along with the widespread use of Geographical Information Systems (GISs), has led to the proliferation of mathematical, GIS-based models (e.g.Montgomery and Dietrich, 1994;Pack et al., 1998;Borga et al., 2002a;Tarolli and Tarboton, 2006;Baum et al., 2008) that can be applied over broad regions to assist forecasting, planning, and risk mitigation.Such models couple a hydrologic model, for the analysis of the pore-water pressure regime, with an infinite slope stability model, for the computation of the factor of safety (i.e. the ratio of retaining to driving forces within the slope) at each point of a landscape.
In particular, with the increasing availability of digital elevation models (DEMs), the topographic wetness index, (TWI - Kirkby, 1975), computed from digital analysis as the ratio between specific upslope contributing area A/b (i.e.upslope contributing area, A, per unit contour length, b) and local slope angle tan β, has been largely used as a general indicator of the influence of topography on soil-water storage dynamics (e.g.Beven and Kirkby, 1979;Lanni et al., 2011) and shallow landslide triggering (e.g.Montgomery and Dietrich, 1994;Wu and Sidle, 1995;Casadei et al., 2003).SHAL-STAB (Montgomery and Dietrich, 1994) and SINMAP (Pack et al., 1998) are among the most popular topography-based slope stability models, where the water table depth is computed based on a steady-state hydrologic balance, and it is expressed as a function of TWI.Borga et al. (2002b) relaxed the hydrological steady-state assumption used in SHAL-STAB by using a modified version of the quasi-dynamic wetness index developed by Barling et al. (1994).This model, called Quasi-Dynamic Shallow Landsliding Model -QD-SLaM, permits us to describe the transient nature of lateral subsurface flow (Grayson et al., 1997).
However, research in the last decade has shown that the establishment of hydrological connectivity (the condition by which disparate regions on the hillslope are linked via subsurface water flow, Stieglitz et al., 2003) is a necessary condition for lateral subsurface flow to occur at a point (e.g.Spence and Woo, 2003;Buttle et al., 2004;Graham et al., 2010;Spence, 2010).Lack of or only intermittent connectivity of subsurface flow systems invalidates the assumptions built into the TWI theory (i.e. the variable -and continuum -contributing area concept originally proposed by Hewlett and Hibbert, 1967).Both field (e.g.Freer et al., 2002;Tromp van Meerveld and McDonnell, 2006) and numerical (e.g.Hopp and McDonnell, 2009;Lanni et al., 2012) studies have shown that subsurface topography (and therefore soil-depth variability) has a strong impact in controlling the connectivity of saturated zones at the soil-bedrock interface, and in determining timing and position of shallow landslide initiation (Lanni et al., 2012).However, despite these evidences, most shallow landslide models do not include a connectivity component for subsurface flow modelling.
Here, we propose a new Connectivity Index-based Shallow LAndslide Model (CI-SLAM) that includes the concept of hydrological connectivity in the description of the subsurface flow processes while keeping the simplicity of the topographic index approach needed to conduct large scale analysis.In our model framework, hydrological connectivity is related to the spatial variability of soil depth across the investigated catchments and the initial soil moisture conditions.Vertical rainwater infiltration into unsaturated soil is simulated by using the concept of drainable porosity (i.e. the volume of stored soil-water removed/added per unit area per unit decline/growth of water table level; Hilberts et al., 2005;Cordano and Rigon, 2008).This allows simulation of pore-water pressure dynamics under the assumption of quasi-steady state hydraulic equilibrium and to estimate the time for development of saturated conditions at the soil/bedrock interface.The model incorporates the computation of a characteristic time for describing the connection of these "patches" of saturation.Specifically, it is assumed that an element (x, y) in a hillslope connects (hydrologically) with its own upslope contributing area A(x, y) when the water table forms a continuous surface throughout A(x, y).Once hydrological connectivity is established, the dynamic topographic index developed by Lanni et al. (2011) is used to describe the transient subsurface flow converging to the element in (x, y).
The hydrological module is then coupled with the infinite slope stability equation to derive CI-SLAM, a shallow landslide model which is able to (a) account for the (positive) effect of the unsaturated zone storage on slope stability, and (b) reproduce pre-storm unsaturated soil conditions.This implicitly helps reducing the fraction of catchment area which is categorized as unconditionally unstable (i.e.failing even under dry soil moisture conditions), improving the confidence in model results (Keijsers et al., 2011).
Model testing is carried out in three study sites located in the central Italian Alps.In this area, shallow landslides are generally triggered by local, convective storms during the summer and early fall seasons.Moreover, accurate field surveys provide a description of hydraulic and geotechnical properties of soils, and a detailed representation of soil depth variation as a function of local slope is reported.An inventory of shallow landslides is also available.Finally, the proposed shallow landslide model is compared with a quasi-dynamic model (QDSLaM) in order to gain insight on the potential improvement brought by the new modelling framework.

The hydrological model
Figure 1 schematizes the hydrological model developed here.During a rainfall event the generation of lateral flow is preceded by the development of a positive pressure head (i.e.perched water table) at the soil-bedrock interface.Several researchers (McNamara et al., 2005;Rahardjo et al., 2005;D'Odorico et al., 2005) have shown that vertical flow in the unsaturated soil zone is reduced when the infiltration front meets a less permeable layer (for example, the bedrock layer).Under this condition, the infiltrating rainwater collects at the less permeable soil layer, inducing rapid increases of pore-water pressure and unsaturated hydraulic conductivity (according to the relationship between matric suction head and unsaturated hydraulic conductivity).As a result, a perched water table will form on the surface of the lowconductive layer, and a subsurface flow will move laterally along the upper surface of this layer (e.g.Weyman, 1973;Weiler et al., 2005).Moreover, in the model it is assumed that a generic hillslope element (x, y) receives flow from the related upslope catchment area A(x, y) only when isolated patches of transient saturation become connected with element (x, y) (Fig. 2).
In the model, unsaturated soil conditions through the whole soil profile (i.e.positive suction head or negative pressure head) are used to initialize our model (step 1 in Fig. 1).For each hillslope element (x, y), the time t wt (x, y) needed to build up a perched zone of positive pore pressure at the soil-bedrock interface is computed by using the following expression (step 2 in Fig. 1): where V 0 [L] is the initial storage of soil moisture through the soil profile before of a rainfall event (Fig. 3); V wt [L] is the e concept of hydrological connectivity.Lateral subsurface flow occurs at point is becomes hydrologically connected with its own upslope contributing area storage of soil moisture needed to produce a perched water table (i.e.zero-pressure head) at the soil-bedrock interface (Fig. 3); and I [LT −1 ] is the rainfall intensity assumed to be uniform in space and time.Computation of V 0 and V wt require the use of a relationship between soil moisture content θ [−] and suction head ψ [L], and a relationship between ψ and the vertical coordinate (positive upward) z [L] (Fig. 3).By using the assumption that the suction head profile ψ(z) changes from one steady-state situation to another over time, where ψ b = ψ(z = 0) is the suction head at the soil-bedrock interface.Bierkens (1998) 1980): (2) and (3) we obtain: Based on Troch (1992), the following relationship between the parameters m and n is used in the model: The storage of soil moisture through the soil profile V is obtained by integrating Eq. ( 4) from the bedrock to the ground surface: where L [L] is the soil depth measured along the vertical.
V wt can be obtained by setting a zero-pressure head at the soil-bedrock interface (ψ b = 0): The suction head value at the soil-bedrock interface at a generic time t < t wt (x, y) (i.e.before development of a perched water table), ψ b t , can be calculated by using the concept of drainable porosity f [−] proposed by Hilberts et al. (2005): By using Eq. ( 8), we can derive an expression for dψ b /dt, useful for estimating the suction head at the soil-bedrock interface at a generic time t, ψ b t (step 5a in Fig. 1): where t is the temporal integration time.For t ≥ t wt (x, y), the generic hillslope element (x, y) exhibits a perched water table at the soil-bedrock interface.
However, this does not guarantee the hydrological connectivity between element (x, y) and its related upslope contributing area A(x, y).In fact, due to the heterogeneity of initial soil moisture and soil depth, isolated patches of saturation which do not necessarily connect with point (x, y) may have developed inside A(x, y).We assume that lateral subsurface flow affects the local soil-water storage of point (x, y) when the water table time t wt indicates continuous saturation through A(x, y).Thus, each point (x, y) has two water table characteristic times: (1) t wt , which indicates the local time for the development of a perched water table; and (2) a connectivity time t up wt -given by the maximum value of t wt in A(x, y) -which indicates the time required by element (x, y) to become hydrologically connected with A(x, y).Therefore, a generic hillslope element (x, y) receives flow from its own upslope contributing area starting from t = t up wt (x, y) (steps 3 and 4b in Fig. 1).Details on the formulation of the connectivity time t up wt are given in Appendix A. The value of the lateral flow rate at element (x, y) is then calculated by using the upslope contributing area A(x, y) as a surrogate for lateral flow (Borga et al., 2002b).In particular, we use the method proposed by Lanni et al. (2011) to describe the variable upslope contributing area which changes linearly with time: where A t (x, y) [L 2 ] and A(x, y) [L 2 ] are the time variable upslope contributing area and the steady state upslope contributing area extending to the divide, respectively; t x, y) is given by the combination of t up wt (x, y) and the time τ c (x, y) required to the lateral flow to reach point (x, y) from the most hydrologically remote location in the corresponding drainage area A(x, y): The value of τ c (x, y) is computed based on the method described in Lanni et al. (2011).Therefore, under the assumptions of constant rainfall intensity I in time and space, the positive pore pressure value at the soil-bedrock interface of point (x, y) for a generic time Hydrol.Earth Syst.Sci., 16, 3959-3971, 2012 www.hydrol-earth-syst-sci.net/16/3959/2012/ t ≥ t up wt (x, y), h b t (x, y), is given by (step 6b2 in Fig. 1): where β is the local ground inclination, K sat (x, y) [L T −1 ] is the saturated hydraulic conductivity and A t /b [L] is the time variable contributing area at time t per unit contour length.

The coupled hydrological slope stability model, CI-SLAM
For hillslopes it is common to define the safety factor as the ratio between maximum retaining forces, F r , and driving forces, The slope is stable for FS > 1, while slope failure occurs when the critical state FS = 1 (such that F r = F d ) is achieved.Lu and Likos (2006) derived a formulation to compute the factor of safety of an infinite slope model that accounts for saturated/unsaturated zones.If the failure surface is located at the soil-bedrock interface, then the Lu and Likos' factor of safety can be written as: with c [FL −2 ] = effective soil cohesion; ϕ [ • ] = effective soil frictional angle; γ w and γ [FL −3 ] = volumetric unit weight of water and soil, respectively; and S e [−] = relative saturation degree.Equations (14a, b) allow the (positive) role played by suction head on the hillslope stability to be taken into account.In this work, locations that are neither unconditionally unstable or unconditionally stable (i.e.locations that are stable when saturated) will be called conditionally unstable as proposed in the pioneering work of Montgomery and Dietrich (1994).
By coupling the hydrological model (Eqs.9 and 12) with the slope stability model (Eqs.14a, b) the factor of safety for conditionally unstable locations (x, y) at a generic time t reads: (15b)

Intensity-duration-frequency relationship for extreme storms
The variability of rainfall intensity with rainfall duration for a specified frequency level is often represented by the intensity-duration-frequency (IDF) relationship proposed by Koutsoyiannis et al. (1998): with I F (d) = rainfall intensity that can be exceeded with a probability of 1−F.ς F and m F are parameters estimated by least squares regression of I F (d) against rainfall duration d.It has been shown (Burlando and Rosso, 1996) that a Gumbel simple scaling model describes well the distribution of annual maximum series of rainfall in the Central Italian Alps.
Based on this model, the rainfall intensity I F (d) can be determined as: with ε = Euler's constant (∼ 0.5772).ς 1 and m can be estimated by linear regression of expectations of rainfall depth against duration, after log transformation, whereas the value of the coefficient of variation (CV) can be obtained as the average of coefficients of variation computed for the different durations, in the range of durations for which the scaling property holds.y T R is given by: where T R [T] is the return period.By combining Eqs. ( 17) and ( 18), T R can be written as a function of rainfall intensity and duration:   Most of the shallow landslides analyzed in this work were triggered during the falls of 2000 and 2002 as a result of relatively short duration storms.The landslide inventory described in this work is part of a more comprehensive archive of shallow landslides which has been described in other papers as well (Borga et al., 1998(Borga et al., , 2002a(Borga et al., ,b, 2004;;Tarolli et al., 2008Tarolli et al., , 2011) ) and executed with a common surveying methodology.The landslides were surveyed during the fall of 2005 and their spatial distribution is reported in Fig. 4, where only the initiation areas are shown.Most of the landslides are found in the lower portion of the basins, where the terrain is steeper and the soil may be deep enough to trigger slope instability.During the survey, the shallow landslides which were evidently induced by forest roads were marked and identified as such and excluded from the model analysis.

Soil depth survey and the relationship with the local slope
Soil depth is an important input parameter in hillslope hydrology (Tromp Van-Merveeld and McDonnell, 2006), but its estimation is usually missing in landslide literature, where often a soil of uniform depth is assumed to stand over an impermeable bedrock.The spatial distribution of soil depth is controlled by complex interactions of many factors (topography, parent material, climate, biological, chemical and physical processes) (e.g.Summerfield, 1997;Pelletier and Rasmussen, 2009;Nicotina et al., 2011).As a result, soil depth is highly variable spatially and its prediction at a point is difficult.Moreover, a soil depth survey is time consuming, and soil depth is difficult to measure even for small basins (Dietrich et al., 1995).Various methods have been explored to allow the estimation of soil depth over landscapes.A process-based approach was suggested by Dietrich et al. (1995) for predicting the spatial distribution of colluvial soil depth.Based on this approach, topographic curvature may be considered a surrogate for soil production.Heimsath et al. (1997Heimsath et al. ( , 1999) ) validated the relationship between curvature and soil production based on observations of cosmogenic concentrations from bedrock in their Tennessee Valley site in California.This approach was incorporated into a landscape evolution model by Saco et al. (2006) to evaluate the dependence of soil production on simulated soil moisture.Roering et al. (1999) supported the idea that soil production follows a non-linear equation and, therefore, modified the dependence of soil depth relationships.However, the various modelling approaches for predicting soil depth over landscapes, described above, showed only partial success (Tesfa et al., 2009).
In contrast to the process-based approaches, a number of studies have applied statistical methods to identify relationships between soil depth and landscape topographic variables (e.g.slope, wetness index, plan curvature, distance from hilltop, or total contributing area) (e.g.Gessler et al., 1995;Tesfa et al., 2009;Catani et al., 2010).Some of these works reported good predictive capabilities for these statistical relationships.For instance, Tesfa et al. (2009) report that their statistical models were able to explain about 50 % of the measured soil depth variability in an out-of-sample test.This is an important result, given the complex local variation of soil depth.
For the purpose of this work, we use a statistical approach for the estimation of the spatial distribution of soil depth over the study catchments.This is based on the availability of a rather dense sample of soil depth measurements.During the fieldwork, conducted in the fall of 2005, 410 direct point measurements of soil depth were made.Survey locations were chosen to represent the range of topographic variation in the areas of model application.Measurements were carried out by driving a 150 cm long, 1.27 cm diameter, sharpened, copper-coated steel rod graduated at 5 cm intervals vertically into the ground until refusal.The advantage of the depth to refusal method is that it is a direct and simple measurement of soil depth.A disadvantage is that the measurement is biased to underestimating the actual depth to bedrock, since there is uncertainty as to what actually causes refusal.Each point measurement is represented by the average of the measures from two or three replicates taken very close each other (at less than 0.5 m distance) (three replicates were taken when the difference between the first two was more than 20 cm).In order to represent soil depth over the extent of grid-size topographic elements, the survey was carried out by taking five point measures 2-3 m apart over a 5 m size grid.A "grid-size" soil depth observation was obtained by taking the average over the five point measures.This permitted us to obtain 82 grid-size observations, which were divided into two subsets: the first (49 grid-size observations) was used to identify and calibrate the statistical model, the second (33 grid-size observations) was used to perform a validation.
The field measurements allowed us to derive the following relationship between soil depth L and local slope tan β: L(x, y) = 1.01 − 0.85 tan β(x, y). (20) Equation ( 20) is limited to local slope less than 45 • (below 2000 m a.s.l.) and 40 • (above 2000 m a.s.l.).In fact, locations with local slope angle larger than these threshold angles are generally characterized by rocky outcrops or very shallow soil thickness and discontinuous soil coverage.The elevation of 2000 m a.s.l.defines a threshold in soil pedological properties (as reported also by Aberegg et al., 2009).The analysis of the soil type distribution showed that Episkeletic Podzols and Dystri-Chromic Cambisols predominantly appear on slopes between 1400 to 1900 m a.s.l.Enti-Umbric Podzols are characteristic for southern exposures at altitudes higher than 2000 m a.s.l.
Other topographic variables, such as plan curvature and specific catchment area, and land cover attributes showed no statistically significant relationship with soil depth.The relationship between soil depth and slope identified for the study watersheds is consistent with findings reported in the literature (Saulnier et al., 1997;Tesfa et al., 2009).
In order to examine the predictive power of the statistical relationship, Eq. ( 20) was applied over the 33 grid-size observations considered for testing purposes.The comparison between modelled and observed soil depth is reported in Fig. 5.We compared the testing data set with the model soil depth values at testing locations using the Nash-Sutcliffe efficiency coefficient (NSE) (Nash and Sutcliffe, 1970).NSE is a normalized model performance measure that compares the mean square error generated by a particular model to the variance of the observations.NSE is equal to 0.4 for our testing case, with a mean error equal to 7 %.This shows that the statistical relationship given by Eq. ( 20) reproduces the variance of the observations considerably better than the simpler model represented by the mean soil depth.While it is clear that a significant uncertainty affects the modelled results, we should note that the reported NSE is in the range of model performance statistics reported in the literature.For instance, Tesfa et al. (2009) reported, for their statistical model applications in Idaho, NSE values ranging from 0.26 to 0.52.

Model application
CI-SLAM has been applied on the study catchments by using the hydraulic and mechanical soil parameters reported in Table 1, identified based on a field survey carried out during the summer season 2005.The soil properties are assumed to be the same for all the three catchments.Based on observations on the erosion crowns and along the forest roads, the survey revealed that the trees in this area are characterized by shallow root systems that spread laterally with small vertical sinker roots that penetrate deeper into the soil.Owing to these observations, we decide not to consider the root strength contribution into the shallow landslide stability analysis.Since the factor of safety calculated by the infinite slope stability equation is fairly insensitive to the values of tree surcharge (Borga et al., 2002a), we omitted considering this factor too.
The soil moisture initial conditions were assumed to represent average climatic conditions based on estimated evapotranspiration fluxes and inter-storm duration statistics which are typical of the seasons where shallow landslides were recorded (summer season and first half of the fall season).These unsaturated soil moisture conditions correspond to considerable cohesion which is due to capillarity, as conceptualized in the generalized principle of effective stress (Lu and Godt, 2008;Godt et al., 2009).
We used the procedure reported by Borga et al. (2005) to estimate the following scaling parameters of the IDF relationship (Eq.18): CV = 0.42, m = 0.48, ς 1 = 13.7 mm h 1−0.48 .These parameters are kept spatially uniform across the three catchments.
The model was not implemented on areas likely characterized by bedrock outcrops, i.e. on topographic elements characterized by slope exceeding either 45 • (for elevation < 2000 m a.s.l.) or 40 • (for elevation ≥ 2000 m a.s.l.).These areas amount to low percentages on Cortina and Fraviano (with 0.7 % and 1 % of catchment area, respectively), whereas they are much more considerable on the steeper Pizzano catchment, where they amount to 11.7 %.
Two general procedures may be considered for shallow landslide model application: diagnostic and predictive (Rosso et al., 2006).With the first procedure, terrain stability is simulated for a given temporal pattern of rainfall intensity and for given initial soil moisture conditions.This allows exploration of the pattern of instability generated by specific storms and could be used to make real-time forecast of shallow landslides.In the predictive mode, the lowest return period of the critical rainfall is computed for each conditionally unstable cell in the landscape.This procedure is well suited for generating maps of shallow landsliding susceptibility, and can be specifically adapted to assess CI-SLAM's capability of predicting shallow landslides which are generated by multiple storms with different storm depth and intensities.The predictive procedure is adopted in this work based on the following steps.First, the critical duration d c of rainfall which generates instability (i.e.FS = 1) is computed for a range of constant rainfall intensities I (ranging from 5 to 60 mm h −1 at a 5 mm h −1 step) which are kept uniform both in time and in space.Then, the return period T r is computed for each considered (I, d c ) pair by using Eq. ( 19).Finally, the lowest return period for each conditionally stable location is selected.The map of the return period of the critical rainfall will provide a representation of the susceptibility to shallow landsliding across the landscape.

Results and comparison with QDSLaM
By using the predictive procedure discussed in the previous paragraph, we derived the shallow landslide susceptibility map of Fig. 6, where the surveyed landslides are also reported.The criterion of shallow landslide susceptibility is based on the return period of the critical rainfall; higher return period values represent medium (T R = 30-100 yr) and low (T R > 100 yr) shallow landslide propensity, and lower return period values represent high (T R = 10-30 yr) and very high (T R < 10 yr) shallow landslide propensity.A "very low" level of shallow landslide susceptibility is assigned to unconditionally stable points (i.e.locations that are stable when completely saturated).
The landslide susceptibility map indicates that the catchments can be subdivided into two geomorphological units.Topographic elements with very high shallow landsliding susceptibility (T R < 10 yr) are reported only in the lower portion of the three catchments, where also surveyed shallow landslides are reported.Conversely, areas characterized as "unconditionally stable" or with very low landsliding susceptibility (T R > 100 yr) are found mostly in the upper portion of the catchments, where no landsliding activity was observed.
The assessment of the predictive power of CI-SLAM is carried out by superimposing the map of the observed landslides onto the map of return period of critical rainfall causing slope instability.The analysis is performed by comparing the proportion of catchment area placed in the various   critical rainfall ranges with the corresponding fraction of the landslide area (Table 2).Good model performances are expected when a large percentage of observed landslides and a small percentage of catchment area occur for low values of return time.For example, this is the case for the Pizzano basin where the percentage of catchment area with T R in the range of 0-10 yr is equal to 1.8 % (34.1 % in the range of 10-30 yr), while the corresponding fraction of observed landslide area is equal to 51.6 % (41.4 % in the range of 10-30 yr).On the other hand, the percentage of landslide area with T R > 100 yr is only 1.1 % versus 43.9 % of the catchment area (including the locations classified as unconditionally stable).Therefore, CI-SLAM is able to correctly classify with a high or very high levels of shallow landslide susceptibility most of the observed landslide areas.This is confirmed by the results for the Cortina and the Fraviano catchments, with this last one showing the best model predictions (63.8 % of landslide area falling in the 2.5 % of catchment area with T R ≤ 10 yr).
Our results suggest that model predictions capture a high percentage of observed landslides at the expense of some overprediction of slope instability.However, one should note that overprediction of landsliding susceptibility has been observed in other applications of topographic index-based shallow landsliding models (e.g.Dietrich et al., 2001) and may be due to a number of reasons, including: (i) inaccurate topographic data, (ii) legacy effects of previous landslides, and (iii) limitation of the landslide surveys.

Comparison with QDSLaM
To assess the role of the variable soil depth and infiltration in generating landslides, we found it useful to compare the predictive capability (i.e. the ability to discriminate areas of  high slope failure hazard) of CI-SLAM with that of the quasidynamic model QDSLaM (Borga et al., 2002b).QDSLaM is based on coupling a hydrological model to a limit-equilibrium slope stability model to calculate the critical rainfall necessary to trigger slope instability at any point in the landscape.The hydrological model assumes that flow infiltrates to a lower conductivity layer and follows topographically-determined flow paths to map the spatial pattern of soil saturation based on analysis of a "quasidynamic" wetness index.With respect to the model proposed in this paper, QDSLaM does not consider the following aspects: (i) vertical rainwater infiltration into unsaturated soil; (ii) analysis of the connectivity to compute the quasidynamic wetness index; and (iii) soil depth variability.All remaining aspects of the modelling framework are described in a consistent way by the two models.Both models have been applied to the three catchments by using the same parameters set.A map of shallow landsliding susceptibility obtained by using QDSLaM is reported in Fig. 7, whereas Table 3 reports the corresponding percentages of slope-stability categories in terms of catchment area and observed landslide area in each range of critical rainfall frequency (i.e. return period T R ).Examination of Fig. 7 and of Table 3 shows the considerable impact of the areas considered unconditionally unstable in QDSLaM.The percentage of topographic elements considered unconditionally unstable ranges from 7.7 % (Fraviano) to 9.9 % (Pizzano).This overrepresentation is reported both in the lower and in the upper portions of the catchments, where high local slope values are present.CI-SLAM does not predict unconditionally unstable locations.In fact, the contribution of negative pressure head (Eq.15a) ensures the stability of steeper topographic elements (i.e.locations with tan β ≥ tan ϕ for cohesionless soils) that would be otherwise classified as unconditionally unstable by QDSLaM (as well as by other traditional landslide models, e.g.Montgomery and Dietrich, 1994;Wu and Sidle, 1995;Pack et al., 1998, among others), which does not account for the role of negative pressure head on soil shear strength.
We also carried out a comparison of the areas which are considered as conditionally unstable by both models.The procedure consists of comparing the proportions of catchment area which fall beneath various return time levels to the corresponding fraction of the landslide area.To compare the two models, the T R values are set so that the same percentage of terrain elements falls beneath the values.The two sets of unstable regions, resulting from the application of the models, are partially overlapping but are not the same.Then the percentage of observed landslide area within each unstable region is computed.The model with the higher percentage provides a better prediction of landslide hazard.This assessment may be repeated for various T R values, obtaining two empirical distribution functions F B (T R ) and F L (T R ),defined for the terrain elements and for the observed landslide cells, respectively.F B (T R ) and F L (T R ) represent the fraction of catchment area and of landslide area, respectively, characterized by return time less than T R .
A function is defined by reporting F L (T R ) versus F B (T R ) (Fig. 8).In the figure, better model performance would be reflected as a steeper curve, which indicates larger differences between fractions of catchment and observed landslide area corresponding to a given T R value.The "naive" model, which predicts that the distribution of slope instability occurs in direct proportion to the terrain area mapped for each threshold, is represented in Fig. 8 by the 1 : 1 line.The figure shows that, for a given percentage of basin area with value of T R less than the threshold, many more observed landslides fall in the area mapped by CI-SLAM.For instance, 60 % of the observed landslides fall in the 10 % less stable fraction of the basin defined by using CI-SLAM, whereas only 26 % of the landslides fall in the corresponding fraction of the basin defined by using QDSLaM.This shows that CI-SLAM provides a better representation of the susceptibility to shallow landsliding with respect to QDSLaM.At the same time, this quantifies the impact that the combination of various hydrological processes (transport in the unsaturated zone, connectivity and soil depth) that are not taken into account by QDSLaM has on shallow landslide triggering.

Summary and conclusions
The shallow landslide model developed in this study, CI-SLAM, takes into account both the vertical infiltration in unsaturated soil and the lateral flow in the saturated zone in modelling of the local pore-water pressure, and introduces in its framework the concept of hydrological connectivity.
CI-SLAM's hydrological module is based on the idea that lateral flow occurs when a perched water table develops over the whole upslope area, and identifies a connectivity time for such a condition to be achieved.This connectivity time represents therefore the time lag (from the onset of rainfall) required for a point in the basin to become hydrologically connected with its own upslope contributing area.For time less than the connectivity time, vertical infiltration is simulated by using the concept of drainable porosity under the assumption of quasi-steady state hydraulic equilibrium.For times greater than the connectivity time, a dynamic topographic index allows us to describe the transient lateral flow dynamics.Moreover, unlike the traditional, lateral flow-dominated, topographic index-based models, CI-SLAM is able to account for partial soil saturation, which, in turn, affects the soil shear strength used for modelling slope stability.A spatially variable soil depth is the main control of hydrological connectivity in the model.
Model performance was evaluated over three catchments located in the central Italian Alps, where reliable intensityduration-frequency relationships of extreme rainfall events, detailed inventories of shallow landslides, and estimates of soil depths are available.We found that in the case studies CI-SLAM provides a reasonably correct estimation for failure initiation probability.CI-SLAM has also been compared to the simpler QDSLaM conceptual model, showing better performances in predicting shallow landslide activities, as quantified by the statistical analysis of the results.
of the water table time t wt encountered along this flow path.This highest value is assigned to each new cell encountered downslope until a higher value is encountered.This can be done because of the use of the D8 flow algorithm which assumes that each cell has a unique downslope flow direction.Therefore, when a flow path P2 converges in a pre-processed path P1, P2 is terminated if it contains a water table time lower than the encountered water table time in P1.On the other hand, P2 continues downslope to modify P1 with the highest upslope water table time.
Thus, each grid cell in the basin has both a t wt value, which indicates the local time for the development of a perched water table, and a connectivity time t up wt , which defines when a cell is hydrologically connected with its own upslope contributing area.

Fig. 1 .
Fig. 1.A flow chart depicting the coupled saturated/unsaturated hydrological model developed in this study.

Fig. 2 .
Fig. 2. The concept of hydrological connectivity.Lateral subsurface flow occurs at point (x, y) when this becomes hydrologically connected with its own upslope contributing area A(x, y).

Figure 3 .
Figure 3. i(z) and i(z) are, respectively, the initial water content and the initial suction5

Fig. 3 .
Fig. 3. θ i (z) and ψ i (z) are, respectively, the initial water content and the initial suction head vertical profiles.θ wt (z) and ψ wt (z) represents the linear water content and suction head vertical profiles associated with zero-suction head at the soil-bedrock interface.

Figure 4 .
Figure 4. Study catchments.The map shows the location of the three catchments, and the3

Fig. 4 .
Fig. 4. Study catchments.The map shows the location of the three catchments, and the landslide distribution.

Figure 6 . 8 Fig. 6 .
Figure 6.Patterns of Return period T R (years) of the critical rainfalls for shallow landslide

Figure 7 .Fig. 7 .
Figure 7. Patterns of Return period T R (years) of the critical rainfalls for shallow landslide

Fig. 8 .
Fig. 8.Comparison of CI-SLAM and QDSLaM: relationship between cumulative frequencies F L (T R ) and F B (T R ) for the study area.

Table 1 .
Hydraulic and mechanical soil parameters relative to the three investigated catchments.

Table 2 .
Percentages of catchment area (C) and observed landslide area (L) in each range of critical rainfall frequency (i.e. return period T R ) for CI-SLAM.
a C = catchment area; b L = landslide area

Table 3 .
Percentages of catchment area (C) and observed landslide area (L) in each range of critical rainfall frequency (i.e. return period T R ) for QDSLaM.
a C = catchment area; b L = landslide area