Physically based modeling of rainfall-triggered landslides: a case study in the Luquillo Forest, Puerto Rico

This paper presents the development of a rainfall-triggered landslide module within a physically based spatially distributed ecohydrologic model. The model, Triangulated Irregular Networks Real-time Integrated Basin Simulator and VEGetation Generator for Interactive Evolution (tRIBS-VEGGIE), is capable of a sophisticated description of 5 many hydrological processes; in particular, the soil moisture dynamics is resolved at a temporal and spatial resolution required to examine the triggering mechanisms of rainfall-induced landslides. The validity of the tRIBS-VEGGIE model to a tropical environment is shown with an evaluation of its performance against direct observations made within the Luquillo Forest (the study area). 10 The newly developed landslide module builds upon the previous version of the tRIBS landslide component. This new module utilizes a numerical solution to the Richards equation to better represent the time evolution of soil moisture transport through the soil column. Moreover, the new landslide module utilizes an extended formulation of the Factor of Safety (FS) to correctly quantify the role of matric suction in slope stability 15 and to account for unsaturated conditions in the evaluation of FS. The new modeling framework couples the capabilities of the detailed hydrologic model to describe soil moisture dynamics with the Inﬁnite Slope model creating a powerful tool for the assessment of landslide risk.


Introduction
Landslides represent one of the most important landform shaping process.All 50 states in the United States have recorded some landslide activity with minor or major socio-economic consequences (Schuster and Highland, 2001;Swanston and Schuster, 1989).Human activity has accelerated the occurrence of landslides by modifying vegetation distribution (i.e.root-soil cohesion), altering the hydrologic regime and disrupting the natural equilibrium.Recent studies of the impacts from landslide activity in the US alone predict 25-50 deaths per year and an annual average total cost of up to $ 3.2 billion (including direct costs and indirect costs) (Sidle and Ochiai, 2006).
Landslides are categorized by both the mode of movement as well as the material being transported during an event.Some movements are connected to specific climatic triggering factors, such as rainfall.These movements are generally shallow and rapid when triggered by large precipitation events on slopes with wet antecedent conditions; deeper and slower when the accumulated precipitation volume gradually builds over an extended time period eventually leading to stability failure.This study focuses primarily on very fast and shallow earth displacements triggered by intense and short rainfall events, such as debris slides (Sidle and Dhakal, 2002).For these movements, the assumption that the failure plane is oriented approximately parallel to the soil surface can be accepted, allowing for the use of the infinite slope hypothesis, which introduces greater simplification to the modeling process.
The initiation of rainfall-triggered landslides is strongly coupled to various hydrological processes, which impact the mechanisms that ultimately cause slope failures.Hydrologic triggering of landslides is a result of both a decrease in the shear strength due to an increase of pore water pressure and an increase of the soil weight that increases the destabilizing forces on a slope.For these reasons, the knowledge of spatio-temporal dynamics of soil water content, groundwater, infiltration and vegetation processes are of considerable importance for the understanding and prediction of landslides (De Vita and Reichenbach, 1998).
The most common practice followed when developing landslide forecasting models is to couple the hydrological model with a stability model to assess the factor of safety (FS).The Infinite Slope model is commonly used in the literature (Arnone et al., 2011;Burton and Bathurst, 1998;Capparelli and Versace, 2010;Montgomery and Dietrich, 1994;Pack et al., 1998;Rosso et al., 2006;Simoni et al., 2007;Wu and Sidle, 1995).It assumes that the plane of failure is parallel to the slope resulting in a static well defined equation for the calculation of FS (Taylor, 1948).Most models can be grouped according to two main features: (1) the hydrological framework used to describe the soil moisture dynamics, i.e. based on a modified approach of wetness index (Arnone et al., 2011;Borga et al., 2002;Burton and Bathurst, 1998;Claessens et al., 2007;Montgomery and Dietrich, 1994;Pack et al., 1998) or based on an approximation of Richards' equation (Baum et al., 2002;Capparelli and Versace, 2010;Iverson, 2000;Simoni et al., 2007;Tsai and Yang, 2006); and (2) the failure criterion used to derive the FS equation, which sometime accounts for unsaturated conditions (Meisina and Scarabelli, 2007;Montrasio and Valentino, 2008;Yeh et al., 2006) or not (Arnone et al., 2011;Tarolli et al., 2011;Simoni et al., 2007;Baum et al., 2002;Iverson, 2000;Burton and Bathurst, 1998;Pack et al., 1998;Montgomery and Dietrich, 1994).
A detailed description of literature and the infinite slope model equations are provided in Arnone (2011) and Arnone et al. (2011).The approach used in these two works is based on a modified formulation of the relative wetness index (Montgomery and Dietrich, 1994) as a function of soil moisture in the unsaturated zone and water table dynamics.This methodology neglects the cohesive effect of negative pore pressure, i.e. matric suction, of unsaturated soils, which should be accounted for in the failure criterion.In fact some failure events can be caused solely by the decrease of matric suction occurring after a rainfall event (Lu and Godt, 2008).The approach proposed by Iverson (2000), based on the Richards' equation, tried to overcome this limitation introducing matric suction into the FS equation.However, since the failure criterion is still based on a either dry or fully saturated scheme, this approach can lead to errors in the calculation of the FS.
These formulations use the soil-water characteristic curve to predict the shear strength function for an unsaturated soil, in addition to saturated shear strength, This work builds on the previous version of the tRIBS (Triangulated Irregular Networkbased Real-time Integrated Basin) landslide component (Arnone et al., 2011).The main differences with the previous version are both in the soil moisture dynamics model and in the failure model.First the tRIBS-VEGGIE (VEGetation Generator for Interactive Evolution), compared to tRIBS, utilizes a numerical solution to the Richards equation to better represent the time evolution of soil moisture transport through the soil column resulting in a more accurate representation of the soil moisture profile.Particular attention is given to both vertical and lateral redistribution, and most of those hydrologic factors important in the assessment of landslides are resolved.Secondly, the landslide module is extended using the Bishop formulation to correctly quantify the role of matric suction in slope stability and to account for unsaturated conditions in the evaluation of FS.
Being a process driven model, it can be applied to any type of basin where weather and basic soil information exist.This model can be utilized to test different rainfall scenarios at the catchment scale to identify high-risk storm characteristics, and possibly be integrated into an early warning system.Because of its vegetation component, the model can be also used to test the effects of land use changes in terms of hydrologic regime's variations and its effects on slope stability.
The newly developed model is here used in the Rio Mameyes basin, Puerto Rico.This basin displays particular features that make the territory prone to landslides (Lepore et al., 2012).
The structure of the manuscript is the following: in Sect. 2 we introduce the two main components of the model and in Sect. 3 the study area.Section 4 describes the data used in both the validation and in the modeling exercise.Section 5 shows the results of a model validation based on direct observation of soil moisture.Section 6 reports an in depth analysis and discussion of the response of the model to an intense rainfall event both at point and catchment scale.Section 7 concludes the manuscript.

Physically distributed hydrology model: tRIBS-VEGGIE
The tRIBS-VEGGIE is an eco-hydrological model that consists of a spatially distributed physically based hydrological model coupled to a model of plant physiology and spatial dynamics (Ivanov et al., 2008a,b).
The tRIBS model has the capability to explicitly consider the spatial variability in precipitation fields, land-surface descriptors and the corresponding moisture dynamics, stressing the role of topography in lateral soil moisture redistribution and accounting for the effects of heterogeneous and anisotropic soil.Basin hydrological response can be simulated at very fine temporal (sub-hourly) and on an irregular spatial mesh (10-1000 m).The irregular mesh allows for the use of smaller computational elements to represent regions within the basin of interest while using larger elements for less critical regions.This framework offers a flexible and reliable computational structure that reduces the number of computational elements without significant loss of information (Donati and Turrini, 2002).Some of the hydrological components presented in the previous tRIBS model were modified, the most critical to this study is in relation to the infiltration scheme.The original tRIBS is based on the kinematic wave approximation of Cabral et al. (1992) and Garrote and Bras (1995), which models the evolution of the wetting front and top front and leads to four saturation states (unsaturated, perched, surface, and completely saturated); this is the framework used for the development of the previous version of landslide module (Arnone et al., 2011) and as well as of an erosion module (Francipane et al., 2012).The tRIBS + VEGGIE model instead is based on a numerical approximation of the one dimensional Richards equation (Hillel, 1980) that governs the fluid flow into the unsaturated soil.
The dynamics of each computational element are simulated separately, but spatial dependencies are introduced by considering the surface and subsurface moisture transfers among the elements, which affect local dynamics via the coupled energy-water interactions.Consequently, when applied to a catchment, the model offers a quasi-three-dimensional framework by which lateral moisture transfers and difference in topographic characteristics may lead to the spatio-temporal variability of states (Ivanov et al., 2008b).
VEGGIE models the spatial distribution of vegetation types (e.g.trees vs. grasses), and pertinent properties of the vegetation state (e.g.vegetation fraction, leaf area index, biomass etc.) (Ivanov et al., 2008a,b) and enriches the hydrology model by accounting for vegetation characteristics that affect both water and energy balance dynamics.Above ground, the canopy has the ability to intercept and retain precipitation that affects the soil infiltration process, and the energy budget.Below ground, vegetation influences the vertical distribution of soil moisture via plant water uptake, and to the cohesion of the soil-bedrock-root matrix, depending on the distribution of roots.VEG-GIE utilizes the Plant Functional Type (PFT) vegetation classification scheme (Smith et al., 1997), and allows for a single computational element to simulate multiple PFTs that may differ in life form (e.g.tree, shrub, grass), vegetation physiology (e.g.leaf optical and photosynthetic properties), and structural attributes (e.g.height, leaf dimension, or root), each interplaying differently with the hydrological response of the basin (Bonan et al., 2002).

Slope stability module
The stability model involves the assessment of FS by applying limit equilibrium analysis.The Infinite Slope model is commonly used in literature and it provides a well determined static equation for the calculation of FS (Taylor, 1948), as described in details in Arnone et al. (2011).The approach used in Arnone et al. (2011) for the landslide module accounts for the soil moisture in the unsaturated zone and water table dynamics, but it neglects the effect of matric suction when computing FS and the corresponding shear strengths in unsaturated soils.
In the last decades, the literature has proposed new procedures to predict the shear strength of unsaturated soil (Vanapalli et al., 1996;Fredlund et al., 1996;Khallili and 1339 Discussion  , 1998;Oberg and Sallfors, 1997;Bao et al., 1998) based on an extended failure criterion of Mohr-Coulomb, (Bishop, 1955), which utilizes the soil-water characteristic curve in addition to the saturated shear strength parameters (cohesion and friction angle).According to this formulation, the shear strength is expressed as:

Khabbaz
where τ f is the shear strength of unsaturated soil; c is the effective cohesion; σ n is the normal stress state on the failure plane; u a is the air pressure and thus (σ n − u a ) is the net normal stress state on the failure plane; ϕ is the friction angle; u w is the water pressure and thus (u a − u w ) is the matric suction on the failure plane; χ is a parameter which depends on the degree of saturation of the soil.The value of χ is assumed to vary from 1 to 0, which represents the variation from a fully saturated condition to a total dry condition.A comparative assessment of several formulations to evaluate the χ term has been carried out by Vanapalli and Fredlund (2000).They considered three soil types and different ranges of matric suction to compare the measured and predicted values of unsaturated shear strength for each of the formulations, and investigate the performance for different soil types and matric suction ranges.It followed that performance of each formulation is more sensitive to the matric suction ranges than the type of soils.As soils approach the saturation, the performances are almost the same.Vanapalli et al. (1996) parameterization is utilized for this study because of its simplicity.For a given soil the χ term is a function of the volumetric water content θ and the residual and saturated soil moisture contents, θ r and θ sat , which are both parameters of the soil texture: The FS equation can be derived by substituting Eq. ( 2) in Eq. ( 1) and using equilibrium of forces according to the infinite slope scheme.Then, assuming the pore-air pressure pressure head (ψ), the FS reduces to: with ψ equal to matric suction in unsaturated soils and pressure head in saturated soils.Equation ( 3) is equal to the equation proposed by Iverson (2000) but for the term in parenthesis ((θ −θ r )/(θ sat −θ r )), not present in his work.In Eq. ( 3), this term modulates the overall contribution of the suction to FS, suggesting an overestimation of FS by the Iverson formulation.Under saturated conditions (θ = θ sat ), the two approaches reduce to the same expression.
As a result of the multi-layer representation of soil moisture within tRIBS + VEGGIE, the FS can be computed locally with no need to make an a priori hypothesis on the depth of the failure surface.The final product of the module is thus a spatially distributed vertical FS profile that takes into account the local moisture and soil conditions within the computational element: the first depth at which the condition FS ≤ 1 is obtained, is designated as the failure surface depth and consequently represents the depth used to estimate the volume of detached material.During a rainfall event, this modeling framework can monitor the time evolution of FS, and by setting different warning thresholds, can dynamically define a spatial distribution of instability levels of the basin.

Study area
The island of Puerto Rico is located in the northeastern Caribbean and is the smallest island of the Greater Antilles.With its roughly rectangular shape, the island is characterized topographically by flat coastal areas and two mountain ranges, the Cordillera Central, which spans east to west with a peak of 1338 m, and the Sierra de Luquillo which will be the focus of this study, with Pico del Este, or East Peak, at 1075 m.The Sierra de Luquillo was chosen for this study as it contains the Luquillo Experimental forest (LEF) which is part of both the Long-Term Ecological Research (LTER) and of the Critical Zone Observatory (CZO) networks.This study will focus on the Rio Mameyes basin (hereto referred to as "Mameyes Basin"), in the northeast of the Sierra de Luquillo and within the LEF boundaries (see Fig. 1).The LEF has been a focal point for studies in landslide impacts on ecology, geomorphology, biology, disturbance and recovery of vegetation (Myster et al., 1997;Scatena and Lugo, 1995;Shiels et al., 2008;Walker and Shiels, 2008;Walker et al., 1996).The Mameyes Basin was selected for the modeling effort based on the availability of specific data and parameters needed to implement tRIBS + VEGGIE.
The Mameyes basin has an area of 16.7 km 2 , and is characterized by a rapid change in elevation from 104.2 m to 1046 m across a horizontal distance of 3 km.An analysis of the slope distribution derived from a DEM (digital elevation model) (see Fig. 2a, d) reveals 10 % of the basin area being characterized by slopes greater than 30 • and 30 % of the basin area with slopes greater than 25 • (see Table 1).The climate of the island is controlled by the easterly trade winds from the Atlantic ocean as well as the pronounced topographic characteristics during synoptically calm days (Garcia-Martino et al., 1996).The basin falls in the windward climatic region of the Island (see Table 1, Aspect) making it one of the wettest basins on the island.Because of the strong gradient in elevation, rainfall, cloudiness and temperature vary consistently throughout the basin.The mean annual precipitation (MAP) varies from approximately 3000 mm, measured at an elevation of 352 m (Bisley Tower), to 5000 mm at higher elevations (Lepore et al., 2011).Table 1 reports some of the climatic characteristics of the basin.
In terms of vegetation, the Luquillo forest has been described as a mix of lower montane wet tropical, wet subtropical and rain forest (Ewel and Whitmore, 1973).The changes in the forest structure, composition and productivity, denote three major forest types: the Tabonuco forest (Dacryodes excelsa) in the wet subtropical and subtropical rain forest life zones, typically within the 150-600 m elevation range; the Colorado forest (Cyrilla racemiflora) in the lower montane wet and rain forest life zones, within 600 and 900 m, and the Dwarf (cloud) forest, above 900 m (Waide et al., 1998) (see Fig. 2c).Another class present throughout the Tabonuco and the Colorado forests, is the Palm  (2002-2007 and 2009-2010); the available daily rainfall data was disaggregated ad-hoc into storms in order to be consistent with the measured soil moisture data (see Sect. 5).Table 1 reports the main characteristics of the station for the complete available time series (2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009).
Surface data and parameters: The basin was delineated from a 30 m DEM of the island (US Geological Survey (USGS) National Elevation Dataset, NED).The hydrographic TIN method was used to derive the TIN mesh from the raster DEM (Vivoni et al., 2004).This method reduces the total number of elevation nodes while preserving the distribution of topographic attributes as slope and curvatures, as well as hydrographic features.
The TIN is characterized by 2213 points and 1831 voronoi cells, whose areas range from 212 to 63 000 m 2 (with average computational element area of 9000 m 2 ) (Arnone,

2011
). Figure 2a shows the TIN together with the elevation distribution for the basin; the blue lines illustrate the drainage network of the basin, whereas the red line delineates a sub-basin used in the model confirmation (see Sect. 5).The soil data for the LEF was obtained from the USDA Forest Service's International Institute of Tropical Forestry of San Juan, the descriptions of the soil types are presented in USDA/NRCS (2002).The USDA defines 12 soil complexes within this basin, however for this study these were further simplified to 4 soil types (Fig. 2b).The generalized soil properties database compiled by Rawls et al. (1982) is used here as a reference for most of the hydraulic properties of the soils; however, it underestimates the saturated hydraulic conductivity when compared to other studies and measurements taken within the region of interest (NCRS, 2002;Harden and Delmas Scruggs, 2003;Silver et al., 1994).This difference in saturated hydraulic conductivity, which are generally much higher than those listed in Rawls et al. (1982), has been attributed to the abundance of organic matter which increases soil drainage and infiltration volumes in wet montane and tropical environments (Harden and Delmas Scruggs, 2003;Knapp, 1978).Other specific characteristics of the soils present in the Mameyes basin have been documented in previous studies: Harden and Delmas Scruggs (2003) observed a very high lateral subsurface flow parallel to the slope, Silver et al. (1994) measured a decrease of hydraulic conductivity with depth, coupled with an increase in bulk density values, and Simon et al. (1990) documented a clay-rich horizon at 50 cm in dioritic soils.The soil bulk density was also adjusted to match observations made within this basin indicating a range from 600 kg m −3 for the first 30 cm and 1300 kg m −3 for deeper layers (Silver et al., 1994).
Based on this literature, in our modeling we (1) allow for a variable saturated hydraulic conductivity, which decreases with depth, and (2) enhance the lateral redistribution of soil moisture using the saturated and unsaturated anisotropy ratio a r [−] (Ivanov et al., 2008a,b;Sivandran, 2012), a parameter in the model.
tRIBS defines the anisotropy ratio as the ratio of the saturated hydraulic conductivities in the directions parallel to the slope K sp and normal to the slope K sn and thus it is partially responsible of the lateral subsurface flux transfer.Therefore the lateral flux transfer is a topographically controlled process; for example an isotropic soil ratio of 1 (i.e. a r = 1) implies same hydraulic conductivity both vertically and laterally, but it does not imply that the same volume of water will flow in these two directions.Further details on the moisture fluxes model are given in Ivanov (2006).Simon et al. (1990) and Lohnes and Demirel (1973) reported values for cohesive strength and friction angle for some of the geological units, and illustrated the expected high variability of these two quantities.The uncertainty associated with these parameters is large when we consider their spatial distribution across four soil types within the domain.For simplicity and to ensure the occurrence of a useful number of failures, this study will assume spatial homogenous values of cohesive strength (3 kPa) and friction angle (25 • ) over the entire basin.
Ecological data and parameters: In this study vegetation is assumed to be static and VEGGIE does not simulate the temporal variability in Leaf Area Index (LAI), Net Primary Production (NPP), carbon pools, etc. of the forest.In preliminary analyses that allowed for dynamic vegetation, the resulting variability of these quantities over a short time period (less than 2 yr) was negligible; this was also supported by an analysis of satellite observations (MODIS Leaf Area Index) where there was little observed season dynamics in forest structure.However, in future applications, the dynamic vegetation component can be utilized to investigate the recovery of the forest after landslide events, and how the seasonal root dynamics may play a stabilizing role in shallow landslides.Figure 2c reports the actual vegetation map as in Helmer et al. (2002); however, our work assumes Tabonuco forest only, because it is the predominant vegetation type of the basin and it is present where both the meteorological and the soil moisture measurements were taken.The Tabonuco forest is modeled as Broadleaf Evergreen been modeled with a rectangular and homogenous density through a depth of 40 cm, their cohesive effect is included into the soil cohesion term.The remaining parameters used in the vegetation modeling have been extracted from literature (Ivanov et al., 2008a,b;Schellekens, 2000;Wang et al., 2003;Weaver and Murphy, 1990).

Model evaluation
tRIBS + VEGGIE has been previously validated for semi-arid areas showing the capabilities of the model to reproduce most of the hydrologically relevant fluxes and their very intricate feedback relationships with vegetation dynamic (Bisht, 2010;Bisht et al., 2010;Ivanov et al., 2008a,b;Sivandran et al., 2008).The Mameyes basin is very different from the previous validation sites, in term of climate (tropical vs. semi arid), morphology (very steep vs. moderate steep), and vegetation type (BET vs. C4 or C3 grasses/shrubs).
For these reasons, we present here a confirmation of the hydrological model, based on soil moisture data.The available data are nine soil moisture time series (Fig. 3a), measured within 500 m of the Bisley tower.These measurement were part of a broader research study carried out at the LEF focused on the climate change effects on humid tropical forests (T.Wood, personal communication, 2010).The soil moisture measurements are available from May to November of 2008.Three locations each with three time-domain reflectometry (TDRs) Campbell Scientific Model CS616 were used to measure the volumetric content hourly at a 30 cm depth.
For the validation exercise we considered a smaller basin within the Mameyes basin, for faster simulation time and a finer mesh.The sub-basin, with 339 computational elements with an average area of 5000 m 2 , is shown in Fig. 1a  Among the soil parameters presented in Table 2, the value of θ sat corresponds to the upper value suggested by Rawls et al. (1982) for the corresponding soil type or, where available, the highest soil moisture content measured.Several simulations were run varying saturated hydraulic conductivities within the upper and lower limit (see Table 2) and anisotropy values, which varied between 100 and 300.These two parameters showed to be the most important in order to reproduce the observed data.Figure 3a shows the nine soil measurements taken near the Bisley tower.The observed variability is very high in terms of both soil moisture average value and dynamics, even within measurements taken from the same plot.
When comparing model output of two computational elements in the vicinity of the Bisley tower, each with an area of 6000 m 2 , the difference between the point scale measurements and output area scale must be kept in mind.Despite the lack of measured soil hydraulic properties at the observation locations, tRIBS + VEGGIE does a very good job in reproducing the measured data.At some points it is not possible to distinguish between the observed and the modeled series; in particular the maximum value of soil moisture and the wetting dynamic is well reproduced.However, at times, the model cannot correctly reproduce the soil moisture decay, but the overall dynamic is very well captured.This type of analysis cannot be extended to other soil types present in the basin because of the lack of soil moisture measurements.However, given that the hydraulic parameters based on Rawls et al. (1982) used in the validation yielded good results, the rest of this work will utilize the same generalized parameters data set to represent and simulate other soil types in the rest of the basin.

Slope stability analysis
The previous section showed that both anisotropy and hydraulic conductivity were the most important parameters and required careful selection in order to reproduce the observed data.The model captured most of the observed soil moisture dynamics, but the high variability of soil measurements suggested that is not, at present, possible to identify a unique value of anisotropy to use in other locations.This section will focus on the sensitivity of the timing and depth of failures obtained from the slope stability model to the parameterization of anisotropy coefficient.
In order to explore the sensitivity to anisotropy ratio of the slope stability model, three anisotropy coefficients, a r1 , a r2 and a r3 , equal to 1 (isotropic soil), 300 and 500 respectively, were simulated.The same meteorological, ecological and hydrological data used for the validation was used.The initial conditions of each experiment were obtained through a spin-up time of one year, with the meteorological forcing of 2007.
Results are analyzed both at element scale, to evaluate the behavior of slope stability over the whole column of soil, and spatially, to evaluate the distribution of unstable elements.

Landslide analysis at element scale
Model results are shown for three selected elements falling within three different soil types: a clay element (C) with a slope of 28 • (0.49 rad), clay loam (CL) with a slope of 52 • (0.90 rad), and sandy loam (SL) with a slope of 47 • (0.82 rad) (small red polygons in Fig. 2a).The CL and C elements are located upslope, while the SL element is closer to the river network; these elements have been selected in order to have the most representative results in term of slope condition, other elements falling throughout the basin have been analyzed and gave similar results.simulation period.This period corresponds to a rainfall event recorded between the 27 and 28 April 2008 with a peak rainfall intensity of 100 mm h −1 .
Figure 4 shows results relative to the three anisotropy values used for the simulation for the CL element.The corresponding infiltration dynamics and soil moisture pattern are reported on the left panels.The initial wedge along the profile depicts homogeneous soil water content before the rainfall event.Rainfall peaks at t = t a and at t = t b produce an increase of water content in surface layers (red wedge) up to saturated values which then propagates vertically through the profile.Fully saturated conditions are achieved only for a r1 in the first 300 mm at time t c .In fact with a r1 moisture is retained in the soil column for the longest time, when compared to simulations with other anisotropy ratios.Increasing the anisotropy alters the redistribution of soil moisture, and consequently water leave the column faster, depth of infiltration is shallower; the two rainfall events impact only the superficial layers with a depth of about 1250 mm for a r2 and 1000 mm for a r3.
Right panels of Fig. 4 illustrate instantaneous FS vertical profiles at selected time periods (t a , t b , t c , t d ) for the three anisotropy cases.Abrupt changes in the shape of the curves denote a sudden change in soil moisture.When soil moisture is constant, the FS curve simply decreases with depth according to Eq. ( 3), whereas with variable moisture both depth and moisture content play an important role.In particular, Fig. 4a shows that the selected element fails at time t c (red line) and depth of about 400 mm, when the soil is fully saturated.At time t c , critical FS values (< 1.3) are also reached at a depth of about 1000 mm.At time t d (purple line), the soil column is not fully saturated, however the combination of soil weight and moisture level makes this condition unstable deeper in the soil column (FS < 1).
In contrast, for a r2 or a r3 the element never reaches a critical value and FS values are higher for higher anisotropy.Increasing the anisotropy involves higher rate of lateral fluxes, which result in more lateral redistribution of water and in this case lower values of soil moisture, which impacts FS.However, it is still possible to see the influence of the wetting front in the column in the superficial layers, where lower FS values are recorded after the rainfall peak, at t c and t d .Figure 5 shows the profile results for the SL element as a function of the same three anisotropy values, a r1 , a r2 and a r3 .As in the previous case, fully saturated conditions are achieved for a r1 at time t c (Fig. 5a, left panel) which at a depth of 300 mm shows a drop in FS to values only slightly greater than 1 (Fig. 5a right panel) and a second decrease of FS values from 600 mm and deeper.The feature at 300 mm of depth is mainly due to the instantaneous peak of rainfall at t c , whereas the deeper FS feature is due to the previous rainfall amounts propagated to greater depths.In fact, due to the higher hydraulic conductivity and gentler slope of the element, the wetting front infiltrates faster than in the CL case, resulting in a narrow wedge, and values of soil moisture close to saturation down to a depth of 850 mm, where the element fails two hours after the first rainfall peak (t c ). Moving from a r1 to a r2 and a r3 , propagation of the infiltration wedge is limited to the more superficial layers (Fig. 5b, c left panel) and FS is always greater than 1, and thus the element is always stable.
Figure 6 illustrates the results for the selected clay element (C).Water content dynamics for a r1 are significantly different from the previous cases.Infiltration rate is lower because of the lower value of hydraulic conductivity and the great volumes of precipitation generating surface runoff, however clay retains more water and thus the initial condition is wetter than the other two cases.The increase of water content due to rainfall contribution at time t a and t b up to quasi saturated condition is limited to the surface layers down to a depth of about 250 mm.In term of stability (right panel), the most severe value of FS is achieved for a rl at a depth of about 900 mm two hours after the event (t d ), but critical FS values (less than 1.3) are present down to a depth of 2000 mm.Increasing the anisotropy ratio (second and third row) results in a slower percolation rate in the normal direction and limits the vertical propagation of the wetting front (left panels).With regard to slope stability (right panels), FS values fluctuate in response to the wetting front, but never fall below 1.5 and no failure occurs within the 2 m of soil.

Landslide analysis at basin scale
In this section the links between soil moisture patterns over the whole basin and the FS values are evaluated, while changing the anisotropy coefficient.Model outputs return spatially distributed results at selected time period and selected layer depth, providing us the distribution of instantaneous values of a selected variable of interest (in this case FS and soil moisture).The value of the FS at each element is taken either as the minimum value within the soil column or as the first value, from top to bottom, resulting in a failure.It is important to mention that post events dynamics are not simulated here.Whenever an element fails, it is re-processed at the next time step based on the updated soil moisture.
Figures 7 and 8 show the results obtained for a r1 and a r2 .Similar analyses, not reported here, have been carried out for a r3 .They depict the spatial distribution of soil moisture at depths of about 485 mm (top row) and 950 mm (middle row), and the FS spatial distribution with the corresponding histogram of absolute frequency of failures at different depths (bottom rows), for four time steps (from left to righ t a to t d ).
Figure 7a (a r1 ) shows that soil moisture pattern is strongly dependent on soil type distribution over the basin and values of soil moisture are quasi homogenously distributed within each type (t a ).The entire soil clay area (bottom corner of the basin map), because of its properties, shows higher values of soil moisture at both times t a , t b ; the variation through time is small and it rarely reaches full saturation due to the slow dynamics of infiltration.The deeper we move and the longer is the lag in the response of the basin to rainfall forcing (Fig. 7b).The resulting FS spatial distributions are shown in Fig. 7c: black elements denotes areas where FS ≤ 1; from red to white, FS values gradually increases up to values greater than 5. Since the geotechnical parameters (c and ϕ) are assumed homogeneous over the whole basin, the clay zone is expected to be more stable during our simulations because: (a) in this area soil rarely achieves the fully saturated conditions and, as a consequence, (b) the matric suction effect plays its stabilizing role; this effect is particularly strong for a clay due to its fine texture.Together with clay elements, flat areas also result in larger FS values.In contrast, the steepest elements commonly fail.The number of failing voronois cells for each depth (and up to a depth of 2 m) is shown in the absolute frequency histogram reported on top of each map (Fig. 7c).At t = t a and t = t b around 121 voronois cells fail, all at depths larger than 1000 mm; at the time of the highest rainfall peak (t = t c ), the number of failing elements rises to 278 mainly at depths between 385 mm, with 79 failing elements, and 485 mm, with 81 failing elements, due to the high water content present in shallow layers and mostly distributed in clay loam and sandy loam.As we move from t = t c to t = t d , the number of failing elements drops down to around 145 but they now occur at larger depths, between 700 and 1100 mm.
Results for a r2 are shown in Fig. 8 for the same four time steps.The soil moisture pattern at time t a (Fig. 8a, left panel) is similar to Fig. 7.At the next time steps (Fig. 8a, center and right panel) the effect of soil type distribution on the soil moisture pattern is less evident due to the higher anisotropy that is redistributing water laterally.Moving from depth 450 mm to depth 950 mm (Fig. 8b), the soil moisture pattern is more influenced by the topography (slope distribution), which determines the lateral flow directions; in fact, the drainage of the steepest areas is larger than in the previous case, and thus the soil moisture is locally lower.Slope instability occurs only in a few elements at depths deeper than 1500 mm, with a frequency distribution almost invariant in time (Fig. 8c), suggesting that the failing elements are always unstable for these set of parameters.

Analysis of failing elements
Further analyses between failing elements and topographical properties are presented in this section where we look at the relationship between failure depth, number of failing elements, slope and soil moisture (SM).and a r2 respectively.The analysis is repeated for time t b , t c , and t d and presented in the three rows.
The two red lines designate the limits of the region of all possible instability cases; particularly, for Fig. 9, column 2, the upper red line describes the worst conditions obtained from Eq. (3) at a fully saturated state of the most permeable soil (sandy loam, which is also the dominant soil type in the domain) and the lower line is obtained from Eq. (3) at unsaturated conditions for the minimum value of soil moisture at which a failure occurs (0.29).At time t b for column 2, the distribution of the circles within the defined region hints at a curve that can be described by Eq. ( 3) at constant soil moisture value.It follows that over the basin, failures occur at a specific soil moisture value, i.e. for same values of the term (θ − θ r )/(θ sat − θ r ), and, according to Eq. ( 3), shallower failures occur at steeper slopes.Failure at 1000 mm occurs wherever slope is greater than a value of about 46 • (∼ 0.80 rad), except for a few steeper elements (> 52 • ,∼ 0.90 rad) that fail at 1800 and 2000 mm.Gentler slopes fail at depths between 1200 and 2000 mm.At the time corresponding to the rainfall peak (t c ), most of the failures occur at shallower depths and follow the upper red line, meaning that the soil is fully saturated; other elements also fail at deeper depths and unsaturated conditions.
For the last time step, t d , depths increase again, as the wetting front moves deeper.
Column 3 describes the relation between depths and soil moisture.The region of possible unstable states was defined from Eq. ( 3) at fixed slope values, by considering the steepest slope ( 63• ,∼ 1.1 rad, upper red line), and two gentler slopes (35  3, failures occur now for a wider and lower range of soil moisture values, corresponding to unsaturated conditions; consequently, circles are not distributed along a curve.In fact, despite the failing elements mainly fall within the same soil type (see Fig. 8c), the soil moisture values at which failures occur are different and thus each circle belong to different theoretically curves at fixed (θ − θ r )/(θ sat − θ r ).Moreover, the frequency distribution highlights that failures occur at higher depths than those found in Fig. 9, since the weight of soil is more important in triggering failure than the soil moisture state.

Discussion
The initiation of shallow landslides is highly controlled by hydrological processes that govern the amount of water contained in the soils; however the role of hydrology is not always evident in existing coupled hydrological-stability models.The results of this study provide an in-depth investigation of how soil moisture affects the slope stability, in relation to topography (slope distribution) and soil type characteristics.The impact of the most important hydrological parameters (soil properties, hydraulic conductivity, and anisotropy coefficient) is also elucidated.Particularly, the analysis highlighted the following findings: in a permeable soil, initiation of instability mostly occurs at the more superficial layers, where a full saturation condition is most often achieved during precipitation events; in a impermeable soil, the effects of the rainfall event on slope stability are delayed; at equal degree of saturation, i.e. relative unsaturated soil moisture content, the stabilizing effect of matric suction on slope stability is greater in finer soils, where ψ values are higher; This study presents the development and analysis of a rainfall-triggered landslide modeling tool that includes both saturated and unsaturated states and can monitor the FS dynamics at high spatial and temporal resolutions.The modeling of landslide susceptibility is highly dependent on the parameterization of hydrologic characteristics, such as soil moisture dynamics, and on the knowledge of the geotechnical properties of a specific area.In particular, this model is capable of capturing a dynamics of soil moisture, and it has been shown to perform well in reproducing observations in both arid and tropical climates.However, a detailed knowledge of some key geotechnical parameters is still very difficult.This study investigated the response of the landslide model to different ranges of a key hydrological parameter such as the anisotropy, which showed the importance of ( 1 The complex yet flexible framework of this model makes its applications to other areas simple.In fact, the detailed soil moisture description, based on the numerical solution to the Richards equation, allows for computing the FS locally with no need to make a priori assumptions on the failure surface depth.Moreover the extended failure criterion formulation, makes it suitable to applications where failures can occur both in saturated and unsaturated conditions.
Many more open questions could be investigated within this framework: analyze the effect of a small localized disruption on a larger area; the effects on slope stability of changes in the rainfall pattern and intensities due to climate change; the effects of deforestation or reforestation on slope stability.Future extensions of the work will focus on including vegetation dynamics, to show the importance of vegetation in slope stability both through its role in moisture fluxes partitioning and through the stabilizing effects of the root biomass.
The detailed descriptive capability of the model, with the FS evaluated at each time step and at various soil depths, together with the flexibility of the hydrology model in describing soil moisture dynamics in a variety of settings, make the presented work a promising tool for landslides assessment.Environ. Geol., 50, 525-534, 2006. Vanapalli, S. K., Fredlund, D. G., Pufahl, D. E., and Clifton, A. W.: Model for the prediction of shear strength with respect to soil suction, Can. Geotech. J., 33, 379-392, 1996 Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper |

1335
Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper |cohesion and friction angle.The physical properties of the soil, such as texture and pore size distribution, become very important in the factor of safety dynamics.

1337
Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | 2 Methods Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Paper | Discussion Paper | Discussion Paper | Discussion Paper | and the pore-water pressure u w as the product of water density (γ w ) and the Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper |

1341
Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper |

1343
Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper |

Discussion
Paper | Discussion Paper | Discussion Paper | Discussion Paper | Tropical (BET) class, with vegetation height of 20 m and a LAI of 6 m 2 m −2 (Wang et al., 2003; Weaver and Murphy, 1990).The roots of the Tabonuco forests are usually shallow, reaching about 40 cm in depth, and with the particularity of being connected with intraspecific root grafts (Basnet et al., 1993).In our application the root component has 1345 Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | delineated by the red line.The sub-basin covers two soil types of the Mameyes basin, Clay-Loam and Silty-Clay, has a slope distribution limited within 0 • and 30 • and aspects predominately N and NE facing.Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Figure 3b-d illustrates the comparison of three of the 9 soil measurements, black circles, with model results, red solid lines.The two grey vertical lines indicate the time period over which hourly data was not available and was reconstructed (see Sect. 4 Meteorological data for details).The different simulated series are obtained by varying the anisotropy ratio values only, all the other parameters are kept constant.The three different simulations shown in the figure well reproduce the three chosen observations, especially for the second part of the series.

1347
Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Temporal dynamics of the soil moisture profile and the FS are analyzed by focusing the discussion on a time window of 48 h encompassing the most severe event of the Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper |

1349
Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper |

1351
Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Figures 9 and 10 show the distribution in frequency of failure depths (column 1) and scatterplots of slope vs. depth (column 2), and SM vs. depth (column 3) of the failing elements, for the anisotropy ratio values a r1 Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Fig. 10.As clearly shown by column 2 and 1353 Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | ) the correct modeling of the 1355 Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | lateral moisture redistribution and (2) the inclusion of the unsaturated processes in the definition of the instability of a hillslope.

Fig. 1 .Fig. 2 .Fig. 3 .Fig. 4 .Fig. 5 .Fig. 6 .Fig. 7 .
Fig. 1.Puerto Rico is an island of the Caribbean with a roughly rectangular shape.The Mameyes basin is located in the North East and on the most elevated areas of the island.
Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | forest (Pretoea montana), usually on steep and poorly drained sites or where climatic conditions (i.e.high winds) are undesirable for the other two species.We present here the parameter sets required by the tRIBS-VEGGIE model.Meteorological data: Hourly meteorological data is available at few locations within the LEF (http://www.sas.upenn.edu/lczo/).This study used the Bisley tower meteorological data (lat.18.31, long.65.74, 352 m) (see Fig. 2a and Table 1), which has been relatively continuous since 2002 with the exception of brief interruptions.The Bisley tower measures all the inputs needed by tRIBS-VEGGIE at an hourly resolution (wind speed and direction, air temperature, cloud cover, relative humidity, rainfall, incoming shortwave radiation) except for the atmospheric pressure, which was obtained from the NCDC weather station at San Juan airport (http://weather.noaa.gov/weather/current/TJSJ.html).During the period from 7 July to 4 August 2008 only daily data was available at the Bisley tower.To fill this data gap, all the meteorological forcings, with the exception of precipitation and cloud cover, were replaced with average values observed for those days during other available years • , ∼ 0.6 rad, lower solid line and 23 • , ∼ 0.4 rad, dashed line).Failure at steep slopes, can occur for any value of moisture, the wetter is the condition the shallower is the failure.A slope equal to 35• represents the minimum slope at which most of the elements fail, as seen in column 2. At time t b , most of failures occur at a soil moisture of ∼ 0.30, while at time t c and t d , when the wetting front has reached deeper layers, failures occurs at saturated and quasi-saturated conditions (∼ 0.55 or ∼ 0.56 depending on the soil type).With anisotropy ratio a r2 , the soil moisture values of failing elements are more spread out and lower in magnitude, as presented in

Table 1 .
. Waide, R. B., Zimmerman, J. K., and Scatena, F. N.: Controls of primary productivity: lessons from the Luquillo Mountains in Puerto Rico,Ecology, 79, 31-37, 1998.Walker, L. R. and Shiels, A. B.: Post-disturbance erosion impacts carbon fluxes and plant succession on recent tropical landslides, PlantSoil, 313, 205-216, 2008.Physiographic and climatic characteristics of the Mameyes Basin and of the mete-orological stations, Bisley (352 m a.s.l.) tower: Climate -mean diurnal, refers to the average variation of the variables throughout the day (numbers in parenthesis refer to the minimum and maximum of the series); Climate -mean annual, refers to the yearly average of the variables for the two stations.

Table 2 .
Hydraulic soil properties for the four soil types present in the Mameyes basin.