Incorporating the Effects of Complex Soil Layering and Thickness Local Variability into Distributed Landslide Susceptibility Assessments

Abstract

Incorporating the influence of soil layering and local variability into the parameterizations of physics-based numerical models for distributed landslide susceptibility assessments remains a challenge. Typical applications employ substantial simplifications including homogeneous soil units and soil-hydraulic properties assigned based only on average textural classifications; the potential impact of these assumptions is usually disregarded. We present a multi-scale approach for parameterizing the distributed Transient Rainfall Infiltration and Grid-Based Regional Slope-Stability (TRIGRS) model that accounts for site-specific spatial variations in both soil thickness and complex layering properties by defining homogeneous soil properties that vary spatially for each model grid cell. These effective properties allow TRIGRS to accurately simulate the timing and distribution of slope failures without any modification of the model structure. We implemented this approach for the carbonate ridge of Sarno Mountains (southern Italy) whose slopes are mantled by complex layered soils of pyroclastic origin. The urbanized foot slopes enveloping these mountains are among the most landslide-prone areas of Italy and have been subjected to repeated occurrences of damaging and deadly rainfall-induced flow-type shallow landslides. At this scope, a primary local-scale application of TRIGRS was calibrated on physics-based rainfall thresholds, previously determined by a coupled VS2D (version 1.3) hydrological modeling and slope stability analysis. Subsequently, by taking into account the spatial distribution of soil thickness and vertical heterogeneity of soil hydrological and mechanical properties, a distributed assessment of landslide hazard was carried out by means of TRIGRS. The combination of these approaches led to the spatial assessment of landslide hazard under different hypothetical rainfall intensities and antecedent hydrological conditions. This approach to parameterizing TRIGRS can be adapted to other spatially variable soil layering and thickness to improve hazard assessments.

1. Introduction

The prediction of where and when landslide events may occur is greatly needed to reduce fatalities and economic losses. Rainfall-induced shallow landslides commonly occur under conditions of transient infiltration into initially unsaturated soils. Intense rainfall events, especially if preceded by prolonged rainy periods, may cause increasing pore water pressures [1] or loss of apparent cohesion effects [2,3] leading to slope instability. Assessing the timing and potential locations of landslides through mathematical models applied for transient, unsaturated infiltration and slope stability modeling requires understanding unsaturated soil hydrology and mechanics, as well as climate and topography [4].

Distributed modeling of the hydrological response to rainfall in steep layered hillslopes at regional scales represents a challenge for landslide hazard assessments. The timing and spatial distribution of rainfall-triggered landslides depend on both quasi-static and dynamic variables [5]. The quasi-static variables contribute to slope susceptibility and they are related to intrinsic layered soil properties, such as stratigraphic setting, thickness, hydro-mechanical properties, and slope topography affecting drainage and seepage processes throughout soils. In contrast, the dynamic or transitory variables, such as relative soil saturation and strength, represent the landslide-inducing factors. Natural processes, such as climatic or hydrological and human activities, such as cut and fill works or forest clearing, respectively affect dynamic and quasi-static variables, thus conditioning the temporal and spatial patterns of landslides.

Several methods and approaches, qualitative or quantitative, have been proposed and tested for distributed landslide susceptibility assessment [6]. Among qualitative approaches, those heuristic based on geomorphological mapping, analysis of landslide inventories and susceptibility zoning can be cited, while among quantitative methods, physics-based numerical modeling and statistically based classification methods can be recognized [7,8]. Qualitative approaches are subjective and portray susceptibility heuristically and in descriptive terms. Quantitative methods allow for the estimation of probabilities of landslide phenomena numerically in any susceptibility zone [8]. Dynamic variables as well as hydrological behavior of the soils are not considered in empirically based models, which instead are focused on geomorphological evidence [9], statistical analysis of the relationship between slope instability factors and present and past distribution of landslides [10,11,12,13,14], or the analysis of landslide probability considering stability models based on stochastic hydrological modeling [15]. Physics-based models consider dynamic variables showing how landslide triggering is strongly affected by hillslope hydrological and morphological conditions as well as stratigraphic setting of the involved soils [5,16,17,18,19,20]. Although philosophy and usefulness of these approaches have been debated [21], distributed, physics-based modeling of hydrologic response remains a primary method to consider spatial heterogeneity [5] and hydrogeomorphic process interactions [22]. Furthermore, physics-based and statistically based methods are preferred to define landslide susceptibility in quantitative terms [9].

Notwithstanding the scientific advances resulting from distributed hydrologic modeling for slope stability assessment, a realistic representation of hydrologic response on slopes could be further improved by incorporating the influence of soil structure and horizons into hydrologic model parameterizations [23,24]. Often, regional and continental scale hydrologic models treat the entire soil mantle as a single homogeneous, isotropic unit [25,26,27,28], and regardless of scale, hydro-mechanical soil properties are assigned based on textural classification, without evaluating the potential impact of these simplifications. Instead, the analysis of hydrological response at the catchment scale for heterogeneous and layered soil profiles affected by shallow landslides phenomena is more complex and represents a challenging task to be accomplished. In fact, permeability contrasts among soil horizons can impact infiltration processes, thus affecting storage dynamics and pore water pressure distribution [29,30,31,32,33] or causing locally diverted flow [34,35,36,37]. Positive pore water pressure or dissipation of matric suction [38,39,40,41] leads to the decrease of shear strength or increase of driving forces causing slope instability. Subsurface soil horizons influence the geometry of perched water tables within a slope and the corresponding volume of individual slope failures, with important implications for variations in landslide occurrence across landscapes [42]. For this reason, estimating effective anisotropies for unsaturated flow systems [43] and quantifying the influence of subsurface layering on infiltration [44] is relevant for regional scale assessments of landslide potential [37].

Accordingly, the aim of this research is to advance the assessment of landslide hazards at the regional scale by taking into account both the spatial distribution of soil thickness and layering properties when parameterizing the Transient Rainfall Infiltration and Grid-Based Regional Slope-Stability (TRIGRS) model [19]. TRIGRS is a distributed model that simulates vertical infiltration and slope stability within each homogeneous grid cell in response to temporally variable rainfall intensities. It can account for lateral variations in topographic slope, soil hydro-mechanical properties, and rainfall intensity to provide an advanced framework for assessing the distributed hazard of rainfall-induced shallow landslides. However, it lacks an established approach for incorporating complex stratigraphic settings, since homogeneous failure material is assumed. For this reason, we developed a method for integrating complex soil horizons and variable soil cover thickness in TRIGRS to allow more realistic simulations of hydrological response and related slope stability during major rainstorms.

We implemented the approach to parameterizing TRIGRS for the mountain slopes of the peri-volcanic areas of the Campania region (southern Italy), where the layered ash-fall pyroclastic soil mantle influences hydrological response from seasonal to event time scales [45,46,47,48,49]. In detail, we present a method to account for the vertical heterogeneity of soil layering by combining soil units, with contrasting hydro-mechanical properties, into homogeneous soil columns that vary laterally based on known distribution of thickness and layering properties of the pyroclastic soil cover.

2. Overview of the Testing Area

2.1. Geological and Stratigraphic Settings

The proposed approach was implemented in the Sarno Mountains (Campania region, southern Italy), a NW-SE oriented carbonate ridge culminating at Mt. Pizzo D’Alvano (1133 m above sea level-a.s.l.). The test area is located on the northeastern border of the Campanian Plain, at a mean distance of about 20 km from the Somma-Vesuvius volcano. The entire carbonate ridge, as well as those forming Lattari, Salerno, and Avella Mountains (Figure 1), is characterized by structural slopes formed by the erosional retainment of original fault-line scarps and connected to flat summit palaeo-surfaces. The carbonate mountains surrounding the Somma-Vesuvius volcano are formed by a Mesozoic carbonate platform series which was piled up in the Apennine thrust belt during the Miocene compressive tectonic phases [50] and subsequently faulted during Pliocene-Quaternary extensional one, which determined the current mountainous and steep morphological settings [51]. In the last tectonic phase, carbonate units were sunk along the Tyrrhenian border forming a regional semi-graben, which was filled during the Quaternary by pyroclastic alluvial and marine deposits setting up the actual Campanian Plain. Pyroclastic products were deposited since about 200,000 years until historical times by explosive eruptions of the Phlegraean Fields and Somma-Vesuvius volcanoes. Ash-fall pyroclastic deposits filled the Campanian Plain and covered irregularly the surrounding mountain slopes. In particular, a larger part of the volcaniclastic series mantling the Sarno Mountains was derived by the deposition of ash-fall deposits erupted by the youngest Somma-Vesuvius volcano ([52] and references therein).

Test pits carried out in different sites allowed the reconstruction of the stratigraphic setting of the layered volcaniclastic soils along slopes of mountains surrounding the Campanian Plain, showing complex alternating unweathered pyroclastic deposits and pedogenized soil horizons. These ash-fall pyroclastic soils belong to the andosols class, mostly found in regions where active and recently extinct volcanoes are located [53,54,55,56].

The complex and spatially variable stratigraphic setting of ash-fall pyroclastic soils covering Sarno Mountains and its control on landslide susceptibility were the focus of preceding studies [47,57,58,59,60,61].

Often, the stratigraphic setting of these layered volcaniclastic soils is incomplete along slopes due to denudational processes, such as shallow landsliding, which reduces the total thickness where slope angles are greater than about 28° up to their truncation (bedrock outcropping) at values greater than 50° [57,61,62].

To describe the stratigraphic settings of the volcaniclastic series, also considering the pedogenetic horizons formed during periods between consecutive eruptions, a criterion based on the recognition of the principal pedogenetic horizons was adopted [55,63], combined with the litostratigraphic and both pedological and geotechnical classification of soils by means of the Unified Soil Classification System (USCS).

Complete volcaniclastic series were found in areas of the Avella, Sarno, Lattari, and Salerno mountain ranges where original depositional morphologies are conserved (with slope angles lower than 28°). The representative volcaniclastic series is characterized, from the top, by a pedogenized soil horizon (present-day soil; A + B soil horizons; Pt and SM respectively) covering unweathered lapilli pumiceous horizons (C-Cb soil horizons; GW or GP) and intercalated paleosols (Bb soil horizons; SM). At the bottom, a basal paleosol (Bb_basal soil horizon; SM) is always present on the underlying carbonate bedrock (R). The thinning of the pyroclastic soil mantle, occurring on slope angles greater than 28° up to its disappearance on slope angles greater than 50°, was recognized as having a strong influence also on stratigraphic settings of the volcaniclastic series along the slopes. The variation of ash-fall pyroclastic soil thickness was linked to the local slope angle by an empirical model [49,61,62,64], which was used in this research. The reduction of the total thickness was recognized determining the downslope pinch out of the pyroclastic horizons (both C and Bb). Moreover, a further reduction of the total thickness was observed leading to the direct overlying of the B horizon on the Bb_basal horizon [48,61].

2.2. Flow-Type Landslides Involving Ash-Fall Pyroclastic Coverings

Since the 18th century, many rainfall-induced landslide events have affected the urbanized foot slopes of mountain ridges surrounding the Campania Plain [65,66,67,68,69], such as the Avella, Sarno, Lattari, and Salerno Mountain ranges. Among the deadliest debris flows events are those that occurred along the Salerno coast on 9 October 1910 and 24–25 October and 1954, which caused, respectively, 150 and 318 causalities [70]. Following the most recent widespread landslide event of 5–6 May 1998, which occurred in the Sarno Mountains and caused the loss of 160 lives [71,72], the national and international scientific community focused many studies on better understanding the predisposing factors and landslide triggering mechanisms. Rainfall of the 1998 event was measured by the rain gauge located in Lauro town (at footslope of Sarno Mountains; 193 m a.s.l.) with a total amount of 157.8 mm distributed in 48 h [72]. However, this measure appears sereverely uncertain due to great altitutinal and planimetric distances from the May 1998 landslide source areas (located above 800 m a.s.l.). Although a reliable rainfall characterization of this event is not available for calibrating TRIGRS model in this study an approach based on physics-based rainfall thresholds was adopted for the assessment of distributed slope stability.

Typically, these shallow landslides are characterized by the involvement of very thin failures (1–2 m) and source areas of small extent (some tens to hundreds of square meters). Three fundamental consecutive evolutionary stages were recognized in characterizing these complex landslides [73,74,75,76,77]: (1) initial debris slides (soil slips), involving a few tens or hundreds of cubic meters of material; (2) debris avalanches, involving progressively greater volumes of pyroclastic deposits along open slopes entrained by a dynamic liquefaction mechanism [78]; and (3) debris flows, characterized by the channeling of flow-like debris masses into the hydrographic network [78]. Due to the common initial stage, this type of landslide can also be defined as “landslide triggered debris flow” [79]. Given the complex sequence of these evolutionary stages [80], the initial slide stage is consistent. However, depending on slope morphological conditions and availability of ash-fall pyroclastic deposits along the downslope path, the initial slide may also (a) arrest along the slope, not evolving into other phases; (b) evolve directly into a debris flow, when falling into an established channel and entraining the channel-fill deposits; and/or (c) evolve to a debris avalanche only, in the case of open slopes with no channels or with downslope diminishing pyroclastic soil thickness and/or slope angle.

Over recent years, some authors identified the basic factors controlling the susceptibility to initial landsliding, which were recognized in small-scale specific geomorphological and anthropogenic features, such as knickpoints, morphological discontinuities related to outcropping carbonate rocky cliffs and artificial cuts into the pyroclastic mantle, excavated for the construction of mountain roads [67,75,76,81,82,83,84,85,86]. Given the undoubted cause/effect relationship between landslide triggering and rainfall patterns [87,88,89,90,91], the comprehension of temporal and spatial hydrological dynamics is a key factor to understand the landslide triggering processes into these layered soils, which can store a large amount of water over the pressure-head interval ranging between the field capacity and the permanent wilting point [33]. Contrasting conceptual models describing the increase of pore water pressure into the ash-fall pyroclastic soil mantle during heavy and/or long-prolonged rainfall events were proposed. Initially, they were based on considering an outflow from the fractured carbonate [70,92,93] and then on a throughflow process occurring within the ash-fall pyroclastic soil mantle [57,74,94,95]. According to the complex layering of these soils and spatially varying thicknesses, an increase of pore water pressure, critical for slope stability, may occur in specific sectors of a slope leading to a local decrease of shear strength, as a consequence of the reduction of both apparent cohesion and effective stress, and an increase of driving forces due to the increase of the unit weight. In such a framework, antecedent-hydrological conditions [96] are considered as fundamental in predisposing or preventing a slope from failing under a heavy rainfall event [77,97].

3. Data and Methods

3.1. Overview of TRIGRS

The TRIGRS model [19] represents an important advance in modeling slope hydrological response and stability for the spatial and temporal prediction of landslide susceptibility. By the assumption of a single-layered, homogeneous soil cover with spatially variable thickness and initial moisture conditions, the governing equations of TRIGRS are based on a linearized solution of Richards [98] equation proposed by Iverson [99]. The vertical distribution of pressure head through time, for infinite depth is thus computed with:

$$\mathsf{\Psi }\left(\mathrm{Z},\mathrm{t}\right)=\left(\mathrm{Z}-\mathrm{d}\right)\mathsf{\beta }+2\frac{{\mathrm{I}}_{\mathrm{nZ}}}{{\mathrm{K}}_{\mathrm{s}}}\sqrt{{\mathrm{D}}_1}\mathrm{tierfc}\left(\frac{\mathrm{Z}}{2\sqrt{{\mathrm{D}}_1\mathrm{t}}}\right)$$

(1)

where Ψ is the groundwater pressure head; t is time; Z = z/cos δ is the vertical depth; z is the slope-normal coordinate direction; δ is the slope angle; d is the steady-state depth of the water table; β = cos²δ − (I_ZLT/K_S) where: K_S is the saturated hydraulic conductivity in the Z direction; I_ZLT is the steady (initial) surface flux that can usually be approximated by the average precipitation in recent weeks or months that is needed to maintain the initial conditions; I_nZ is the surface flux of a given intensity for the nth interval (i.e., model timestep); and D₁ = D₀/cos²δ where D₀ is the saturated hydraulic diffusivity (D₀ = K_S/S_S; S_S is the specific storage); N is the total number of time intervals for the defined duration of the simulation; H(t–t_n) is the Heaviside step function; and t_n is the time at the nth time interval in the rainfall sequence. The function ierfc is of the form ierfc (η) = $1/\surd \mathsf{\pi }$exp (−η²) − ηerfc(η), where erfc(η) is the complementary error function.

The analytical solution for infiltration in TRIGRS partitions the problem into an upper unsaturated zone, with a capillary fringe developing above a lower saturated zone. Infiltration is characterized by absorbing part of the infiltrated water within the unsaturated zone and percolation of the remaining water at the base, leading to the rise of the shallow water table. Under this condition, TRIGRS uses four parameters including residual water content (θ_r), saturated water content (θ_s), inverse of capillary fringe (α), and hydraulic conductivity (K_sat) to approximate the Soil Water Retention Curve-SWRC [100] and thus the one-dimensional infiltration flux [101], with no lateral flow/throughflow. In this case, the vertical water pressure head changes in the unsaturated zone are thus computed:

$$\mathsf{\Psi }\left(\mathrm{Z},\mathrm{t}\right)=\frac{\cos \mathsf{\delta }}{{\mathsf{\alpha }}_1}\ln \left(\frac{\mathrm{K}\left(\mathrm{Z},\mathrm{t}\right)}{{\mathrm{K}}_{\mathrm{s}}}\right)+{\mathsf{\Psi }}_0$$

(2)

where α₁ = α cos2δ; Ψ₀ is the pressure head at the water table (Ψ₀ = 0) or at the top of the capillary fringe (Ψ₀ = −1/α); and K(Z,t) is the hydraulic conductivity as a function of time and depth in the unsaturated zone [90,101,102].

Below the saturated zone, TRIGRS computes the pressure head increasing for finite depth basal boundary using a formula based on a Fourier series solution:

$$\mathsf{\Psi }\left({\mathrm{Z}}_{\mathrm{w}},\mathrm{t}\right)={\mathsf{\Psi }}_{\mathrm{hn}}\left\{1-\frac{4}{\mathsf{\pi }}{\displaystyle \sum }_{\mathrm{m}=1}^{\infty }{\left(-1\right)}^{\mathrm{m}-1}\frac{1}{2\mathrm{m}-1}\exp \left(-\frac{{\left(2\mathrm{m}-1\right)}^2{\mathsf{\pi }}^2{\mathrm{D}}_1\mathrm{t}}{4{\mathrm{d}}_{\mathrm{LZw}}^2}\right)\cos \left(\frac{\mathsf{\pi }}{2}\left(2\mathrm{m}-1\right)\left(\frac{{\mathrm{Z}}_{\mathrm{w}}}{{\mathrm{d}}_{\mathrm{LZw}}}-1\right)\right)\right\}$$

(3)

where Z_W = Z − d is the vertical depth below the initial water table; Ψ_hn = β_hn is the pressure head applied after the accumulation of water above the initial water table; and d_LZw is the vertical height of the saturated layer (d_LZw = d_LZ − d).

Finally, a one-dimensional infinite-slope-stability analysis [103] is used by TRIGRS to determine the Factor of Safety (FoS), following Equation (4):

$$\mathrm{FoS}\left(\mathrm{Z},\mathrm{t}\right)=\frac{\tan {\mathsf{\varphi }}^{\prime }}{\tan \mathsf{\delta }}+\frac{{\mathrm{c}}^{\prime }-\mathsf{\Psi }\left(\mathrm{Z},\mathrm{t}\right){\mathsf{\gamma }}_{\mathrm{w}}\tan {\mathsf{\varphi }}^{\prime }}{{\mathsf{\gamma }}_{\mathrm{s}}\mathrm{Z}\sin \mathsf{\delta }\cos \mathsf{\delta }}$$

(4)

where ϕ′ and c′ are respectively the soil friction angle and cohesion for effective stress; γ_w is unit weight of water; and γ_s is unit weight of soil. Ψ (Z, t) is the transient pressure head at depth Z and time t obtained from either Equations (1), (2) or (3) depending on the particular conditions modeled. Regarding the unsaturated zone, TRIGRS computes the FoS above the water table multiplying the matric suction Ψ(Z, t)γ_w by χ = (θ − θ_r)/(θ_s − θ_r) [4].

3.2. Parameterizing TRIGRS

Considering results of preceding studies at site-specific scale within the framework of the Sarno Mountains, a multi-scale approach was adopted to assess the temporal and spatial hazard to onset of initial landslides (debris slides) with the TRIGRS model.

The first scale is site specific, related to three observed cases of initial landslides (L1, L2, and L3; Figure 2) that occurred in the Sarno Mountains in May 1998. These three landslides were chosen due to their representativeness of typical geomorphological conditions of initial landslides and studied for the assessment of Intensity-Duration (I-D) rainfall thresholds by coupled hydrological and slope stability modeling [76,77]. Moreover, at the same scale, results of soil hydrological monitoring carried out at an instrumented station installed upslope of the source area of a fourth initial debris slide (L4; Figure 2) [47,58] were also considered. The second scale is distributed across the steep slopes of the Sarno Mountains to assess the region impacted by the debris flow event of May 1998, thus encompassing areas of the site-specific scale (Figure 2). The application of TRIGRS to the site-specific scale was used to calibrate its extrapolation for distributed landslide susceptibility assessment across the Sarno Mountains region.

Considering the small extent of initial source areas (<100 m²), a 5-m resolution Digital Elevation Model (DEM) of the area obtained by LIDAR data made available from the Italian Environment Minister was used for TRIGRS modeling. This high-resolution DEM allowed the reconstruction of a new map of the distribution of ash-fall pyroclastic soil thickness according to the empirical model linking ash-fall pyroclastic soil thickness and local slope angle, which replicates observations that slope is inversely proportional to deposit thickness [59,61,62,64]. Furthermore, the same elevation data were used to calculate inputs required for TRIGRS simulations, including flow direction and slope angle maps.

Since TRIGRS represents the soil as a single-layered column for each grid cell, a key aspect of applying the model to the Sarno Mountains slopes was to translate the properties of the spatially variable ash-fall pyroclastic soil stratigraphy into uniform parameters to assign to each grid cell. To apply this approach, two main aspects of the volcaniclastic deposits were taken into account: (1) their variable spatial thickness and (2) their representative unsaturated/saturated hydraulic and geotechnical properties. At this scope, starting from typical stratigraphic settings carried out at the site-specific scale [47,58,59,60,76] as well as laboratory and field characterizations of unsaturated/saturated hydraulic and mechanical properties, the complex stratigraphic settings were simplified in a single soil horizon with spatially variable thickness.

Data on hydraulic unsaturated/saturated and geotechnical properties of individual soil horizons (Figure 3), and their thicknesses throughout the Sarno Mountains were measured by laboratory and in situ tests carried out on L1, L2, L3 and L4 sites by De Vita et al. [76] and Napolitano et al. [77]. For each of these four test sites, three vertical profiles at the upslope, midslope, and downslope location were selected to account for observed variability in soil horizon thickness.

Considering the variability of unsaturated/saturated and geotechnical properties of each soil horizon, values representative for the 12 vertical profiles were assigned by harmonic means, weighted by thickness of single soil horizons. The harmonic mean of the soil hydraulic properties was considered as representative for simulating vertical infiltration through the horizontally layered soil profile, since the hydraulic resistance is measured perpendicularly to the flow direction. Moreover, the weighted harmonic mean was also considered for representing values of geotechnical soil properties since it will result in a value lower than the arithmetic mean and therefore provides a more conservative approach for slope stability analyses. Finally, median values of weighted harmonic means, obtained for each of the 12 vertical profiles, were then assigned to the ash-fall pyroclastic soil cover model. In particular, the median operator was considered the most representative statistical parameter for incorporating the variability of soil profiles in a single value.

To set the TRIGRS model, the hydrological modeling at the site-specific scale area was carried out using the rainfall intensities (I) adopted for deterministic I-D rainfall thresholds (2.5, 5.0, 10.0, 20.0 and 40.0 mm/h) determined at site-specific scale [77]. Thus, soil water pressure head (h) and FoS time series, simulated under a rainfall of given constant intensity, were analyzed for three cells composing the L4 test site (Figure 2). In particular, these cells (5 × 5 m resolution) correspond to a section that intersects the initial slide area (cell 1) and part of the monitoring station area (cells 2 and 3). Rainfall duration (D) with FoS values equal to 1 of sample cells corresponding to the monitoring station area (L4) were considered. In detail, this approach relies on two fundamental assumptions, derived from in situ observations and VS2D (Variably Saturated Two Dimensional; version 1.3) [104] modeling results (Figure 4): (1) slope areas upslope of morphological discontinuities (rocky cliffs, knickpoints and road cuts), leading to the abrupt reduction of ash-fall pyroclastic soil thicknesses, represent the source areas of initial debris slide phenomena; and (2) in these zones, saturated or near-saturated conditions occur during heavy and/or long prolonged rainfall events due to the concentration of unsaturated flows along the slope.

This initial modeling step allowed the calibration of the TRIGRS model to reach a good match between estimations of FoS derived by the coupled hydrological modeling, performed by VS2D code, slope stability analysis [77], and the TRIGRS results.

According to the TRIGRS structure, the vertical distribution of initial pressure head values under unsaturated conditions are expressed in terms of an equilibrium profile that is defined by assigning the depths of the hypothetical groundwater table below the soil for all individual cells in the model domain. In such a sense, pressure head time series obtained at different depths for each assigned constant rainfall intensity (2.5, 5.0, 10, 20, and 40 mm/h) show how the infiltration process controls the dynamics of the virtual water table as it rises closer to the overlying soil, following sustained rainfall.

It is worth noting that in the study area the regional groundwater table that persists throughout the year is found only within the carbonate bedrock at a depth of hundreds of meters [105,106].

The results of this modeling were compared with the inventory map of debris flows that occurred in May 1998. During this calibration phase, typical winter pressure head values of the ash-fall pyroclastic soil cover, equivalent to those recorded by monitoring activities [47,58], were set as initial conditions. Based on TRIGRS setting, an initial vertical distribution of soil water pressure head within the ash-fall pyroclastic soil cover was defined by an iterative approach varying water table depth from the bottom of the modeled cover. In the second step, by implementing representative soil properties (median value of weighted harmonic means) and initial winter soil hydrological conditions [33], distributed hydrological and slope stability modeling for the region were carried out. Susceptibility maps based on given antecedent hydrological conditions defined by the hypothetical initial water table depths and rainfall I-D thresholds were calculated by considering spatially variable thickness of ash-fall pyroclastic soil cover. These maps were set to show the distribution of unstable areas, or pixels with FoS ≤ 1 values, and likely unstable areas, or pixels with FoS values ranging between 1.01 and 1.05, according to different rainfall intensities and durations. The choice of the FoS range to as-sume as indicative of slope instability were derived by the consideration that these values can be approximated to the limit equilibrium condition (FoS = 1) by truncation. The mapping results were compared with source areas of debris flows occurred in May 1998.

4. Results

4.1. TRIGRS Model Calibration

Representative values of unsaturated/saturated hydraulic and geotechnical soil properties (

Table 1. (a) Unsaturated and saturated hydraulic and geotechnical soil properties determined for principal ash-fall pyroclastic soil horizons [76,77]; (b) Unsaturated/saturated hydraulic and geotechnical soil properties used for setting TRIGRS model and derived as median value of weighted harmonic means (WHM), estimated for 12 representative soil columns of four test sites (L1, L2, L3, L4) [76,77]. Keys to symbols: (K_sat) saturated hydraulic conductivity; (θ_s) saturated volumetric water content; (θ_r) residual volumetric water content; (α and n) van Genuchten’s fitting parameters of soil water retention curve; ($\mathsf{\varphi }$’) = effective friction angle; (c’) = effective cohesion.

(a)			Hydro–Mechanical Properties
			Ksat (m/s)	θs (ad.)	θr (ad.)	α (cm−1)	n (ad.)	$\mathsf{\varphi }$’ (°)	c′ (kPa)
Soil horizons (USDA)		B	4.82 × 10−5	0.505	0.083	0.884	1.307	32.0	4.500
		C	2.82 × 10−3	0.500	0.001	20.39	1.081	37.0	0.000
		Bb	6.00 × 10−6	0.663	0.001	0.884	1.307	34.0	1.800
		Bbbasal	2.48 × 10−7	0.505	0.083	0.884	1.307	35.0	8.100
(b) Test site	Soil column and thickness (m)		Weighted Harmonic Mean
			Ksat (m/s)	θs (ad.)	θr (ad.)	α (cm−1)	n (ad.)	$\mathsf{\varphi }$’(°)	c′ (kPa)
L1	1a	4.85	8.53 × 10−6	0.590	0.132	3.059	1.294	34.0	3.820
	1b	4.34	4.78 × 10−6	0.581	0.111	3.851	1.316	33.8	4.350
	1c	1.62	1.68 × 10−6	0.584	0.118	3.590	1.307	33.9	4.186
L2	2a	3.23	5.68 × 10−6	0.589	0.143	2.689	1.299	34.0	3.493
	2b	3.13	5.30 × 10−6	0.565	0.122	3.618	1.372	33.5	3.884
	2c	1.75	2.06 × 10−6	0.567	0.099	4.384	1.356	33.5	4.585
L3	3a	3.99	2.65 × 10−6	0.580	0.073	5.174	1.302	33.8	5.496
	3b	3.35	3.34 × 10−6	0.573	0.098	4.347	1.335	33.7	4.659
	3c	2.43	5.20 × 10−6	0.557	0.122	3.689	1.399	33.3	3.806
L4	4a	3.21	2.88 × 10−6	0.574	0.089	4.673	1.328	33.7	4.959
	4b	3.09	3.78 × 10−6	0.556	0.080	5.109	1.385	33.3	5.027
	4c	1.47	7.16 × 10−6	0.558	0.097	4.512	1.386	33.3	4.539
TRIGRS model			4.28 × 10−6	0.574	0.105	4.099	1.332	33.7	4.445

), which were estimated as median of weighted harmonic means respectively of previous characterizations [76], allowed us to set the simplified regional-scale parameterization of the ash-fall pyroclastic soil cover used for TRIGRS modeling (Figure 5).

The values were estimated considering the entire soil cover, conbining soil properties of all horizons and their thickness variability from 12 representative soil columns of four test sites (L1, L2, L3, L4; Figure 2).

The first conceptual step was to understand how the TRIGRS model can simulate at the site-specific scale the hydrological response of ash-fall pyroclastic soil-mantled slopes during a rainfall event. At this scope, a local-scale application of TRIGRS for L4 site (Figure 2) was calibrated based on previous results from hydrological and slope stability modeling. In particular, previous research based on the VS2DTI model [104] illustrated that vertical infiltration and lateral throughflow may lead to near-saturated or saturated zone formation upslope of morphological discontinuities, where a reduction or truncation of the ash-fall pyroclastic soil cover occurs [76,77]. This process, which is strongly dependent on stratigraphic setting, local morphological conditions, and unsaturated/saturated hydrological properties, controls the regime of soil water pressure head in ash-fall pyroclastic soils, leading to slope instability [76]. Therefore, considering that only vertical flows (infiltration) are modeled in TRIGRS without lateral hydraulic connection between adjacent cells (or pixels), lateral flows (throughflow) were accounted by assuming a vertical rise or lowering of a hypothetical water table, below the bottom of the pyroclastic cover, rather than to explicitly simulate slope-parallel flows. Thus, the influence of lateral flow is accounted for through the initial conditions, though it is not possible to include the event-based lateral-flow response within the current TRIGRS framework.

During the simulated rainfall, a gradual increase of pore water pressure head (h) and a decrease of FoS values are simulated in cell 1 and to a lesser extent in cells 2 and 3 (Figure 6). For constant rainfall intensities of 2.5, 5.0 and 10.0 mm/h, rising of the water table leads to saturated and/or near-saturated conditions with different durations, while saturation is reached at roughly the same duration for 20 and 40 mm/h.

These conditions strongly affected cell 1 where unstable conditions were simulated due to the higher slope angle. In addition, FoS > 1 characterized cells 2 and 3, despite the increasing of soil water pressure head up to 0 or with slightly positive values at depth, presumably due to higher cover thicknesses.

The second step applied to calibrate the TRIGRS model was to consider rainfall intensity (I) and duration (D) values leading to slope instability conditions in cell 1, which allows definition of a deterministic rainfall threshold for the L4 test site. The latter was compared with those obtained by using the VS2D code and slope stability of an infinite slope model for the other three test sites (L1, L2 and L3) [77]. The results show a good correlation between rainfall I-D values calculated by means of TRIGRS and those belonging to the other three thresholds, especially for lower rainfall intensities (Figure 7), allowing the extrapolation of hydrological and slope stability modeling from the site-specific to the scale of the Sarno Mountains region.

4.2. Slope Stability Maps for Initial Landslides

Based on the calibration of TRIGRS, distributed slope stability maps of the Sarno Mountains slopes were developed considering the same constant rainfall intensities applied at the local scale (2.5, 5.0, 10.0, 20.0 and 40.0 mm/h) and variable durations. The resulting maps show the spatial distribution of stable and unstable pixels, computed for the same time step of the slope failure obtained previously for the site-specific scale (Figure 8A, Figure 9A, Figure 10A and Figure 11A). Pixels with computed FoS values <1.0, and ranging between 1.0 and <1.05, were considered unstable or likely unstable, respectively. Finally, slope stability maps were compared with the location of initial landslides that occurred in May 1998 as well as with pre-May 1998 landslide events, including the four source areas used for the local-scale analysis (Figure 8C, Figure 9C, Figure 10C, and Figure 11C). However, since the rainfall time-series for May 1998 was not considered explicitly due to the lack od reliable recordings, these comparisons are only to demonstrate that the potentially unstable zones identified with TRIGRS correspond to observed locations of previous slope failures.

These maps show the distribution and FoS of unstable cells across the entire study area. In particular, they are located both close to the main drainage network and along open slope areas, where morphological discontinuities exist causing a local increase in slope angle and a thinning or downslope truncation of the ash-fall pyroclastic soil cover. The areal extension of initial slope instabilities (debris slides) are inversely related to rainfall intensity (Figure 8B, Figure 9B, Figure 10B, and Figure 11B) as the expansion of unstable areas increases with the reduction of the rainfall intensity, albeit with a corresponding increase in event durations. In detail, unsaturated throughflow can converge in small zones, where thickness and drainage section of the ash-fall pyroclastic soil cover reduce due to morphological factors, leading to the localized formation of near saturated/saturated conditions (Figure 4; [76,77]). In these areas, the FoS values became lower than 1.05 and then slope instability conditions were simulated, due to the combined effects of slope angle values higher than 35° and the reduction of the cover thickness from 5 m to 1 m or less, and the applied antecedent conditions in these convergent zones with thinner soils. Moreover, the map showed that stable zones are located in areas of low slope angles and high ash-fall pyroclastic soil thicknesses.

Firstly, unstable/likely unstable cells were observed starting from 89 h for 2.5 mm/h, 46 h for 5.0 mm/h, 20 h for 10 mm/h, and 13 h for 20 mm/h and 40 mm/h, respectively, with cumulative rainfall of 222.5 mm, 230.0 mm, 200.0 mm, 260.0 mm, and 520.0 mm. The increasing rainfall intensities, which eventually exceed the infiltration capacity of the soil, would result in some runoff that is not considered explicitly by TRIGRS. Thus, the higher intensities may require increased durations before landslides are triggered within this modeling framework, which has been identified as a possibility in other settings [107]. The spatial variability of the unstable/likely unstable cells, obtained for each rainfall intensity and duration modeled, are more evident in the enlarged area of the maps (Figure 8C, Figure 9C, Figure 10C and Figure 11C). The comparison of the landslides events that occurred on 5th and 6th May 1998, in particular with those affecting the four test sites, showed that the model simulated many of the cells for the triggering areas (initial debris slide), but performed more accurately for lower rainfall intensities than for higher ones. In addition, in some cases, the distribution of the unstable cells along slopes falls within, or are coincident with, the border of the debris avalanche and flow zones. Other potential landslide-prone areas, not affected by slope instability phenomena during May 1998 and preceding debris flows events, were correctly identified by the TRIGRS simulations as remaining stable during the rainfall conditions considered. Moreover, other unstable cells not coinciding with the landslide inventory were identified.

A second analysis was carried out considering failure times of initial slides characterizing the 5th and 6th May 1998 landslide event, starting from obtained distributed slope stability maps for each rainfall intensity. In detail, initial slides areas (geometries) were combined with those coinciding with only FoS ≤ 1 cells. The resulting times at failure were identified as corresponding to the condition of one or more cells with FoS ≤ 1 contained within geometry of each source area. Therefore, frequency of initial slides, depending on rainfall intensity and duration, were analyzed and compared with deterministic rainfall thresholds estimated at the local scale by Napolitano et al. [77] (Figure 12). The comparison showed that simulated unstable cells coinciding with triggering areas (initial debris slides) did not occur before the duration of the I-D thresholds (

Table 2. Duration of rainfall with constant intensity leading to slope failure [77] and durations to slope failure for given frequencies of May 1998 source areas, obtained by TRIGRS modelling.

	Rainfall Intensity (mm/h)	2.5	5.0	10	20	40
Winter Threshold [77]	Duration (hours)	77	47	22	11	6
Frequency of May 1998 source areas, unstable under TRIGRS modelling	1%	90	47	22	14	14
	5%	92	49	23	15	15
	50%	109	55	33	27	27
	95%	170	120	77	49	34

) but in the range of 90 h for 2.5 mm/h, 47 h for 5.0 mm/h, 22 h for 10 mm/h, 15 h for 20 and 40 mm/h, respectively, with cumulative rainfall values of 225.0 mm, 235.0 mm, 220.0 mm, 300.0 mm, and 600.0 mm.

Considering both analyses, the results observed for susceptibility maps obtained for 20 and 40 mm/h appear dependent on TRIGRS modeling capability that does not explicitly consider runoff generation potentially occurring for highest intensity rainfall conditions.

5. Discussion

Mathematical models of rainfall infiltration consider dynamic variables showing how landslide triggering is strongly affected by hillslope hydrological and morphological conditions as well as stratigraphic setting of the involved soils [5,16,17,18,19,33]. These models can be applied to predict the temporal and spatial variation of debris flow-type landslide susceptibility, as already done in the Oregon Coast Range area, USA [4], in Umbro-Marchean Appennine area, central Italy [108,109], or in peri-Vesuvian areas, southern Italy [81,85]. However, all of them reveal how difficult it is to incorporate the influence of soil structure and horizons into parameterization of distributed physically based models.

Multi-layered soil coverings are usually discretized as a single homogeneous and isotropic soil unit with uniform hydro-mechanical properties whose area is assigned by empirical relationships based on textural classification. According to such an issue, the aim of this research was to advance the assessment of distributed landslide hazard at a regional scale by taking into account both spatial variations in soil thickness and composite stratigraphic settings and properties of soil horizons to define uniform soil properties for each TRIGRS grid cell. Considering the results of previous hydrological monitoring and modeling activities [33,58,77], a coupled distributed hydrological and slope stability model was tested for ash-fall pyroclastic soils covering the Sarno Mountains area, extending the approach from the site specific to the entire mountain range. Other attempts to model hydrological and slope stability behavior of Sarno Mountain slopes have considered layered soils, but they have not taken into account the combined influence of the spatially variable soil thickness and antecedent soil hydrological conditions, as well as the influence of soil layering on the model parameterization for distributed applications [81,85,94].

The common approaches used to define susceptibility maps of rainfall-induced flow-type landslides in this area consider geomorphological and stratigraphic factors including slope angle, cover thickness and presence of natural or artificial discontinuities [84,86]. However, the highly variable and contrasting hydraulic properties of ash-fall pyroclastic soil horizons can impact the modeling of infiltration processes and thus storage dynamics and pore water pressure distribution [33,110,111,112]. To advance this critical aspect, the approach used in this research represents a further step toward a more consistent hazard assessment of the temporal and spatial occurrence of initial landslide triggering debris flows [78] considering effective layering of involved materials. Accordingly, the TRIGRS code was parameterized using an innovative method which considers the variability of ash-fall pyroclastic soil thickness and hydraulic properties, as well as the role of slope position on antecedent conditions. Therefore, the essential point of the research was the definition of a representative physical model of ash-fall pyroclastic soil cover. At this scope, for the parameterization of TRIGRS, unsaturated/saturated hydrological and geotechnical properties were defined as the median of harmonic means of values determined by field and laboratory testing [76,77]. A first calibration at the site-specific scale was obtained taking into account deterministic rainfall I-D threshold previously estimated by a slope hydrological modeling, carried out by VS2D code, coupled with hydrological and slope stability modeling for winter antecedent hydrological conditions [77]. The match between rainfall threshold and failure conditions modeled by TRIGRS at the local scale was used as a calibration, which justified a distributed extrapolation across the slopes of the entire Sarno Mountains range. The results obtained from our approach can be considered an advance of those previously obtained [81,85], since they are based on consistent physical modeling of the ash-fall pyroclastic soil cover, considering both the spatial variability of soil cover thickness as well as reliable estimates of layered pyroclastic soil hydraulic and geotechnical properties.

The obtained slope stability maps, which can be applied to evaluate both spatial and temporal landslide hazards, show zones affected by slope instability, which are consistent with the whole landslide inventory of the Sarno Mountains (Figure 8B, Figure 9B, Figure 10B, and Figure 11B). The results obtained are expressed by potentially unstable areas whose spatial resolution is finer than that identified in literature [81]. Furthermore, a different distribution of FoS values leading to slope instability was observed for each critical rainfall intensity modeled. Areas close to cover discontinuity (rocky cliffs, knickpoints or road cuts) or characterized by high slope angles (more than 35°) resulted in potentially more unstable areas, which is consistent with the prior literature [76,77,81,85]. Cover thickness, as well as rainfall intensity and duration strongly affect the infiltration processes within the cover and then the occurrence of saturated or near-saturated conditions in those “critical” areas [33,113]. These conditions were clearly observed by the comparison of slope stability maps obtained for each rainfall intensity value. In fact, clustered and more localized unstable areas for higher rainfall intensity were observed.

By the comparison of times to slope failure in source areas of the 5th–6th May 1998 landslide event, obtained from distributed slope stability modeling, with the winter rainfall thresholds of Napolitano et al. [77], a good matching was found. In fact, instabilities were recorded after overcoming the threshold time values (Figure 7), for all the rainfall intensities modeled for the point locations (Figure 2). This result validates the proposed approach to extrapolating site-specific findings across broader scales based on the definition of effective hydrological and geotechnical properties of the spatially variable columns of a single-layered soil cover, derived from weighted averaging of corresponding multi-layered pyroclastic horizons. Given that TRIGRS applies the infinite slope approach, which often overestimates unstable areas in comparison to that resulting from field obeservations, the proposed approach can be conceived as conservative.

The results of the distributed slope stability modeling in the Sarno Mountains test area also emphasized the importance of the spatial variability of a rainfall event (cells), whose distribution is strongly affected by morphological conditions. Considering the May 1998 Sarno landslide event, the maps showed many potentially unstable areas located in the western part of the mountains ridge not affected by these instabilities, which instead were confirmed by previous landslide events (for example, the Palma Campania 1986 event) [114]. This means that the triggering rainfall event in May 1998 was probably localized and influenced by local morphological setting and not characterized by a homogeneous distribution across the mountain ridge [115]. Details of spatially and temporally variable rainfall were not captured locally, as described previously. Other unstable areas identified by the modeling, not coinciding with source areas of May 1998 landslides can be related to site specific conditions (such as stratigraphic setting and land cover) which were not modeled explicitly.

Finally, the results obtained both at the local and distributed scales confirmed that slope failures are caused by saturated or near-saturated conditions occurring in areas close to downslope discontinuities of the ash-fall pyroclastic soil mantle controlling reduction of thickness or abrupt interruption, as previously demonstrated by results of hydrological slope modeling at the local scale only [76,77,113].

6. Conclusions

The distributed modeling of the hydrological response at large spatial scales of slopes characterized by layered and thickness-varying soils represents a challenge for the assessment of hazard to rainfall-induced shallow landslides. Due to the difficulty in characterizing and discretizing complex stratigraphic settings and spatially variable thickness of soils into numerical models, assumptions of a single homogeneous, isotropic unit are commonly considered. According to such an aim, a new approach for assessing landslide hazards from site-specific to distributed scales was presented by taking into account both the spatial distribution of soil thickness and layering properties.

The aim of this research was to advance the existing approaches when parameterizing physically based models of hydrological response for applications on steep, landslide-prone hillslopes. The utility of the method was implemented for volcaniclastic soils mantling the steep mountain slopes of the Sarno Mountains where complex layered soils influence seasonal and event scale hydrology.

The proposed physics-based approach highlights the importance of using both laboratory and field data including soil hydrological field measurements to define and improve input data for physically based modeling. In this way, results obtained provide guidelines for estimating effective soil hydrological and geotechnical properties in complex stratigraphic settings and for advancing slope hydrologic and stability models under varying initial and boundary conditions. Moreover, the results emphasize how topography and the type of precipitation phenomena can strongly control the spatial variability of a high-intensity rainfall events, determining its precise and narrow spatial position.

Finally, the definition of the spatial distribution of the landslide occurrence probability, related to a specific value of rainfall threshold, is essential to obtain landslide hazard maps at a distributed scale over a mountainous district. In this way, results of this research could facilitate the development of a consistent early-warning system, taking into account seasonally variable landslide hazards.