Hydrological control of soil thickness spatial variability on the initiation of rainfall-induced shallow landslides using a three-dimensional model

Thickness and stratigraphic settings of soils covering slopes potentially control susceptibility to initiation of rainfall-induced shallow landslides due to their local effect on slope hydrological response. Notwithstanding the relevance of the assessment of hazard to shallow landsliding at a distributed scale by approaches based on a coupled modelling of slope hydrological response and slope stability, the spatial variability of soil thickness and stratigraphic settings are factors poorly considered in the literature. Under these premises, this paper advances the well-known case study of rainfall-induced shallow landslides involving ash-fall pyroclastic soils covering the peri-Vesuvian mountains (Campania, southern Italy). In such a unique geomorphological setting, the soil covering is formed by alternating loose ash-fall pyroclastic deposits and paleosols, with high contrasts in hydraulic conductivity and total thickness decreasing as the slope angle increases, thus leading to the establishment of lateral flow and an increase of pore water pressure in localised sectors of the slope where soil horizon thickness is less. In particular, we investigate the effects, on hillslope hydrological regime and slope stability, of irregular bedrock topography, spatial variability of soil thickness and vertical hydraulic heterogeneity of soil horizons, by using a coupled three-dimensional hydrological and a probabilistic infinite slope stability model. The modelling is applied on a sample mountain catchment, located on Sarno Mountains (Campania, southern Italy), and calibrated using physics-based rainfall thresholds derived from the literature. The results obtained under five simulated constant rainfall intensities (2.5, 5, 10, 20 and 40 mm h−1) show an increase of soil pressure head and major failure probability corresponding to stratigraphic and morphological discontinuities, where a soil thickness reduction occurs. The outcomes obtained from modelling match the hypothesis of the formation of lateral throughflow due to the effect of intense rainfall, which leads to the increase of soil water pressure head and water content, up to values of near-saturation, in narrow zones of the slope, such as those of downslope reduction of total soil thickness and pinching out of soil horizons. The approach proposed can be conceived as a further advance in the comprehension of slope hydrological processes at a detailed scale and their effects on slope stability under given rainfall and antecedent soil hydrological conditions, therefore in predicting the most susceptible areas to initiation of rainfall-induced shallow landslides and the related I-D rainfall thresholds. Results obtained demonstrate the occurrence of a slope hydrological response depending on the spatial variability of soil thickness and leading to focus slope instability in specific slope sectors. The approach proposed is conceived to be potentially exportable to other slope environments for which a spatial modelling of soil thickness would be possible.


Introduction
Landslide occurrence is a continuous and widespread geomorphological process that occurs under specific environmental conditions (Glade 2003) as a result of different triggering mechanisms (such as earthquake, precipitation, volcanic activity, etc.) (Sidle and Ochiai 2006). Shallow landslides (debris or earth flows) triggered by rainfall events represent one of the most catastrophic natural hazards because they occur suddenly and can potentially travel a long distance and at high velocity (Iverson 2000;Formetta and Capparelli 2019). Because of their disastrous effects in terms of fatalities, casualties and damage to infrastructure, they represent a relevant individual and societal risk that must be evaluated by a modern approach to natural risk management (e.g. Guzzetti 2000;Dai et al. 2002;Bell et al. 2005;Salvati et al. 2010). Unfortunately, the assessment of landslide hazard and risk is a complex and uncertain task which requires the combination of different scientific disciplines, techniques and methodologies. In recent years, two approaches have been commonly used for predicting landslide events and implementing landslide early warning systems at local and regional scale: (1) empirical, lumped-statistical approaches, which relate rainfall information to the observed occurrence (e.g. Caine 1980;Crozier and Eyles 1980;Glade et al. 2000;Godt et al. 2006;Guzzetti et al. 2008;Berti et al. 2012;Mirus et al. 2018;Tufano et al. 2019); (2) spatially distributed physics-based deterministic modelling (e.g. Montgomery and Dietrich 1994;Baum et al. 2008;Godt and McKenna 2008;Frattini et al. 2009;Thomas et al. 2018). The first approach allows to investigate wider areas (Galanti et al. 2018), but it is typically affected by significant uncertainties depending on the limited landslide inventories, the unrepresentativeness of rainfall measurements, and the scarce knowledge about antecedent soil hydrological conditions (Nikolopoulos et al. 2014). Moreover, this approach does not consider the dynamic predisposing factors such as the wetness state of the potentially unstable slopes (Bogaard and Greco 2018).
Otherwise, physics-based deterministic modelling, based on hydro-geomorphological and geomechanical models, allows one to estimate the hydrologic response of soil cover under simulated rainfall as well as the related slope stability, but it is usually applied only on small areas due to detailed spatial data required. In particular, results obtained strictly depend on geological frameworks and mechanisms, by which rainfall infiltration conditions slope stability, spatial scale and the related quality and quantity of hydrologic, hydraulic and geotechnical data (Formetta et al. 2014a(Formetta et al. , 2014b. In addition, slope angle, slope shape and soil (regolith) Landslides 18 & (2021) 3367 thickness represent important controlling factors of shallow landsliding (Iida 1999).
Many studies have used the physics-based deterministic modelling to advance the understanding of hydrological response of soil cover involved in shallow landsliding (e.g. Reid et al. 1988;Greco et al. 2013;De Vita et al. 2018), through the implementation of one-(e.g. Srivastava and Yeh 1991;Terlien 1998;Greco et al. 2013), two-(e.g. Crosta and Dal Negro 2003;Capparelli and Versace 2014;Napolitano et al. 2016;Tufano et al. 2016) or threedimensional (e.g. Ng et al. 2001;Frattini et al. 2004;McDougall and Hungr 2004; Lizárraga and Buscarnera 2019) hydrological slope model. Commonly, these approaches use the analytical solution of Richards' equation, firstly implemented by Iverson (2000) to describe near-surface soil water pressure head values, and subsequently applied to simulate the transient infiltration process for unsaturated soils (Baum et al. 2010).
However, hydrological response of natural slopes is very complex because it is strongly influenced by infiltration into unsaturated soil, surface runoff, slope-parallel unsaturated/saturated throughflow, subsurface flows from upstream areas, effect of vegetation and flows through fractured bedrock. All these issues affect the predictive ability of the models and make difficult the interpretation of the results obtained (Capparelli and Versace 2014).
The GEOtop distributed hydrological model Rigon et al. 2006;Endrizzi et al. 2014) solves the 3D Richards' equation and extends it to the case of saturated soils. Moreover, GEOtop is able to model the transient infiltration processes in layered soils as characterised by different thicknesses. Also, GEOtop was found useful for investigating hydrological response and probability of occurrence of rainfall-induced shallow landslides in complex topographic conditions (e.g. Formetta et al. 2016b), when stratigraphy of soils is known (e.g. Formetta et al. 2014a;Formetta and Capparelli 2019) or/and for well-defined bedrock positions (e.g. Simoni et al. 2008;Formetta et al. 2014a).
In this paper, we have tested GEOtop model capabilities to investigate the hydrological slope response of a test area of Sarno Mountains (Campania, southern Italy) for which a validated empirical model, describing the spatial distribution of ash-fall pyroclastic soil thickness, is available (De Vita et al. 2006aFusco et al. 2017). This area was selected given the recurrence of rainfall-induced, very-to-extremely rapid debris flows, involving ash-fall pyroclastic soils covering steep slopes. The most important landslide event occurred on 5-6 May 1998 and caused the loss of 160 lives. After this event, many studies were carried out to recognise the initiation mechanisms of these events to be used for landslide susceptibility assessment. In particular, several authors (e.g. Celico and Guadagno 1998;Guadagno et al. 2005; recognized factors that directly control landslide susceptibility and location of the triggering areas, the role of stratigraphic and morphological settings of natural slope morphological features as well as artificial modifications, including man-made cuts and fills. Moreover, unsaturated properties and hydrological regime of ash-fall pyroclastic soils were also considered as a key factor for understanding the hydrogeological behaviour of the soil mantle (Picarelli et al. 2004;Sorbino 2005;Fusco et al. 2017). Finally, soil antecedent-hydrological conditions were recognised as strongly controlling the effects of the rainfall events on the slope stability (Napolitano et al. 2016).
Distributed models that analyse both vertical and lateral fluxes responsible for landslide triggering at the Sarno Mountains were previously developed (e.g. Frattini et al. 2004;Lizárraga and Buscarnera 2019;Formetta and Capparelli 2019). These models are mainly based on vertical discretisation of domain and on the effect of alternation of pyroclastic soil horizons characterised by different hydraulic conductivity, which are responsible for a local increase of soil pore water pressures (Reid et al. 1988;Johnson and Sitar 1990). However, the effects of spatial variability of soil thickness and stratigraphic setting on hydrological regime and slope stability have not been tested yet by three-dimensional numerical modelling.
Accordingly, to advance the knowledge on how the abovementioned factors condition the local hydrological slope response and the probability of slope failure, in this work we combine (1) spatial variability of thickness and stratigraphic settings of the ash-fall pyroclastic soil mantle (De Vita et al. 2006a, 2006bDe Vita and Nappi 2013;Fusco et al. 2017), including bedrock outcrops, with (2) a three-dimensional hydrological model (GEOtop 2.0, Bertoldi et al. 2006;Rigon et al. 2006;Endrizzi et al. 2014). At this scope, the proposed GEOtop modelling is conceived to be useful to complement large-scale hydrological studies in shallow landslide alert systems.

Description of the study area
The study area ( Fig. 1) is located in the Sarno Mountains which, along with the Avella, Salerno and Lattari mountain ranges, surround the Campanian Plain. The latter represents the most extended alluvial plain of the region, where the Somma-Vesuvius and Phlegrean Fields volcanoes are located. The bedrock of Sarno Mountains comprises a Mesozoic carbonate platform series which during the Miocene was thrusted by Apennines compressive tectonic events (D'Argenio et al. 1973;Mostardini and Merlini 1986;Patacca and Scandone 2007). During the Quaternary, these tectonic units were affected by normal faulting that caused their tectonic lowering along the Tyrrhenian border and the formation of a semigraben where a back-arc volcanic activity began, first forming the Phlegrean Fields (39 k-yrs B.P.) and then the Somma-Vesuvius volcano (25 k-yrs B.P.). These two volcanic centres were the sources of pyroclastic flow and ash-fall deposits that filled the structural depression of the Campanian Plain and covered the nearby mountain ranges. A complete stratigraphic sequence of pyroclastic series, considered as a reference, was identified at the western foot of the Sarno Mountains by Rolandi et al. (2000) where two main pyroclastic complexes were distinguished: the older one, named 'Ancient Pyroclastic Complex' (APC), which is mainly constituted of ash-flow deposits of 'Campanian Ignimbrite' (39 k-yrs B.P.) and deposits of other important eruptions of Phlegrean Fields and Ischia Island; the younger one, named 'Recent Pyroclastic Complex' (RPC), which is formed by deposits of the Somma-Vesuvius volcano eruptions: Codola 25 k-yrs B.P. (Rolandi et al. 2000), Sarno 17 k-yrs B.P. (Rolandi et al. 2000), Ottaviano 8 k-yrs B.P. (Rolandi et al. 1993a) and Avellino 3760 yrs B.P. (Rolandi et al. 1993b). In the RPC deposit of historical and modern eruptions of the Somma-Vesuvius volcano are also included A.D. 79 (Lirer et al. 1973), A.D. 472 (Rolandi et al. 1998), A.D. 1631 (Rosi et al. 1993) as well as some minor eruptions, including the last one that occurred in 1944. The Somma-Vesuvius eruptions that constitute the RPC had dispersion axes

Original Paper
Landslides 18 & (2021) generally oriented eastward, thus mainly involving Sarno Mountains with their ash-fall deposits. Only the A.D. 79 eruption had a southward-oriented dispersion axis, involving mainly the Lattari Mountains. However, along Sarno Mountain slopes, only ash-fall deposits of the RPC are usually found, typically characterised by intercalated unweathered pumiceous lapilli and pedogenised soil horizons (paleosols), formed during the stages between consecutive eruptions. Considering both pedological classification (USDA 2014 and the USCS system), the succession of soil horizons can be summarised as (De Vita et al. 2006a): (1) A horizon, consisting of prevailing humus (Pt); (2) B horizon, mainly characterised by highly pedogenised pumiceous pyroclasts (SM); (3) C horizon, formed by pumiceous pyroclasts, slightly weathered (GW or GP); (4) Bb horizon, corresponding to a B horizon buried by a successive depositional event and thus considered a paleosol (SM); (5) Cb horizon, representative of a buried C horizon (GW or GP); (6) Bb basal horizon, corresponding to a residual pyroclastic deposit, highly weathered by pedogenesis (SM); and (7) R horizon, fractured carbonate bedrock with open joints, infilled by the overlying paleosol for the first few meters. Moreover, the slopes are covered with woods, mainly represented by deciduous chestnuts (Castanea sativa) with sparse deciduous beeches (Fagus sylvatica).
Ash-fall pyroclastic soil coverings are involved in complex landslides (sensu Cruden and Varnes 1996). Specifically, landslide movements are characterised by up to three fundamental evolutionary phases: (1) initial soil slips (Campbell 1975) or debris slides (Cruden and Varnes 1996) involving few tens or hundreds of cubic metres of ash-fall pyroclastic soils (debris, with more than 20% of gravel); (2) debris avalanches (Hungr et al. 2001), involving progressively greater volumes of pyroclastic deposits along open slopes by a dynamic liquefaction mechanism; (3) debris flows (Hungr et al. 2001), characterised by the channelling into the hydrographic network of rapid to very rapid flow-like debris masses. Depending on morphology of the slope and on continuity of pyroclastic cover, the evolutionary stages series can be different: the first phase (debris slide-soil slip) can directly develop into a channelled flow (debris flow) and thus defined as 'landslide-triggered debris flow' (Jakob and Hungr 2005). In some cases, debris slides-soil slips can develop into an avalanche only (debris avalanche) or do not evolve into debris flows, thus stopping their movement along the slope or at the foot slope.
These initial landslides are caused by the rapid growth of soil water pressure head (Campbell 1975;Wilson 1989) and the loss of apparent cohesion (Fredlund 1987).

Distribution of ash-fall pyroclastic deposits
The implementation of a three-dimensional geomorphological model in GEotop requires an accurate knowledge about the spatial distribution of ash-fall pyroclastic soil thickness. To achieve this goal, data and results of preceding studies (De Vita et al. 2006b;De Vita and Nappi 2013) were considered in this research. They were focused on developing a regional distribution model of ash-fall pyroclastic deposits along slopes and, in particular, on assessing total theoretical thickness of ash-fall pyroclastic soils by the algebraic sum of isopach maps of principal Plinian eruptions of the Phlegrean Fields and Somma-Vesuvius volcanoes (Fig. 1).
The latter were known from previous scientific literature (De Vita et al. 2006b;De Vita and Nappi 2013). This modelling approach was supported by test pits along slopes less than 30°in the Sarno and Lattari Mountains. These data were used to compare real thickness values against the theoretical distribution provided by the regional model (De Vita and Nappi 2013). According to fall depositional mechanism, ash-fall pyroclastic soil deposits tend to mantle the slopes generating a bedding parallel to the slope (Fisher and Schmincke 1984), whose theoretical total thickness can be estimated as where z (m) is the stratigraphic thickness, z 0 (m) is the total theoretical thickness fallen and α (°) the local slope angle. In contrast, for slope angle values higher than 30°, real stratigraphic thicknesses were observed as being lower than the theoretical ones in the same area, with a progressive decrease as the slope angle increases, up to negligible values for slope angle values greater than 50°. Such a distribution model was considered chiefly controlled by landsliding, occurring mostly for slope angle values >28° (de Riso et al. 1999). In addition, the increase of slope angle causes the reduction of the number of pumiceous lapilli horizons: for lower values of slope angle (below 30°), up to two pumiceous lapilli horizons can be usually found along slopes of the Sarno Mountains and only one for the Lattari Mountains. Instead, for greater values of slope angle, C and Bb soil horizons pinch out downstream, thus determining the direct contact between the B and the Bb basal horizons. This condition generates a vertical hydrological discontinuity resulting in a localised building up of soil water pressure head (Reid et al. 1988) during heavy rainfall.
Specifically, the empirical relationship between the expected maximum stratigraphic thickness of pyroclastic soil mantle and slope angle (Fig. 2) (De Vita et al. 2006b;De Vita and Nappi 2013) was used here to implement the vertical soil variability of ash-fall soil pyroclastic cover.

Methods and data
GEOtop-FS 2.0 framework The hydrological and slope stability analysis is based on the modelling framework of GEOtop-FS 2.0, which couples (1) the open source three-dimensional fully distributed hydrological model GEOtop 2.0 (Endrizzi et al. 2014) and (2) a space-and time-varying hillslope stability component based on the infinite slope model. The framework is presented in detail in Formetta et al. (2016b) and Formetta and Capparelli (2019). GEOtop models surface and subsurface flows of variably saturated soil, snow cover dynamics and soil freezing, using a three-dimensional finite volume approach that couples heat and water flow equations following the time-lagged approach proposed in Panday and Huyakorn (2004). Unsaturated and saturated flows are computed by solving the three-dimensional Richards equations (Richards 1931) and the van Genuchten (1980) Soil Water Retention Curves (SWRCs). Overland flow is based on the extension of Darcy's law to surface flow as proposed in Gottardi and Venutelli (1993), and channel routing is modelled by the shallow water equation neglecting inertia (Endrizzi et al. 2014). The hillslope-stability model is based on the soil moisture (soil suction) provided for each cell of the computational domain and for different vertical layers by the GEOtop model. It involves the computation of the suction stresses and probability of failure using the infinite slope model for each cell of the computational domain (and assuming in turn the depth of each layer as potential failure surface). The slope stability model implements the general effective stress framework developed by Lu and Likos (2006) and Lu et al. (2010) to compute the soil effective stress: where σ (kPa) is the total stress, u a (kPa) is the pore air pressure and σ s (kPa) is the suction stress (Lu and Likos 2004). The latter represents all the inter-particle stresses such as capillary stress, the electric double-layer force, the van der Waals attractive force and the matric suction of soil, and is defined as in Lu and Likos (2004): where u w (kPa) is the soil pressure head, θ is the volumetric water content, θ r is the residual volumetric water content and θ s is the saturated volumetric water content. The factor of safety FS (ratio of stabilising to destabilising forces) is computed by combining the generalised effective stress and strength failure criteria, for a uniform, infinite slope as expressed in Lu and Godt (2008) and Formetta and Capparelli (2019): where ϕ′ (°) is the effective internal friction angle, β (°) is the slope angle, c′ (kPa) is the cohesion at zero normal stress due to the intergranular bonding stress, γ (kNm −3 ) is the unit weight of soil and Z (m) is the soil thickness. Z in this paper is considered not constant but varying pixel by pixel as specified in the previous section.   through OMS ranges from modelling river flows and snow-melting evolution (Formetta et al. 2011;Formetta et al. 2014bFormetta et al. , 2014cAbera et al. 2017aAbera et al. , 2017b, earth energy balance (Formetta et al. 2013(Formetta et al. , 2016a, to soil moisture and soil temperature modelling (Formetta et al. 2016b(Formetta et al. , 2016c. The modelling framework GEOtop-FS 2.0 receives as inputs meteorological data (rainfall, air temperature, solar radiation and air humidity); raster maps such as digital elevation model, river network, soil types, land cover and soil depth; and physical parameters of soils such as SWRCs, hydraulic conductivity functions, soil cohesion and friction angle for each soil horizon. The model outputs are time-varying maps of soil moisture, soil pressure head, suction stress and failure probability for each soil layer of the computational domain.

Hydrological modelling and slope stability analysis
The setting up of a GEOtop model needs the input of spatially distributed data and the definition of hydrological and geotechnical parameters. At this scope, available LIDAR data, resampled for a resolution of 5 × 5 m, have provided the Digital Elevation Model (DEM) of the computational domain, corresponding to a mountain catchment with an extension of 22,274 m 2 and a mean slope angle of 37° (Fig. 3), which was implemented in the code. The calculation domain and grid size coincide with those of DEM. Moreover, slope angle and aspect maps, needed for domain characterisation, were calculated from DEM by terrain analysis tools of QGIS (release 2.18.11) (Fig. 3). Instead, the land use map, called 'land cover', has been set as a unitary value map without distinction among cultivation types. Finally, the soil cover thickness map (De Vita and Nappi 2013) was implemented in the model.
As already specified, the thickness of ash-fall pyroclastic soil covers is spatially varying (De Vita et al. 2006b), with a depth to the bedrock ranging between 0 (bare soil) and 5.6 m. Higher thickness occurs mostly on the higher and middle elevations of the computational domain (Fig. 4). In addition, the stratigraphic series of ash-fall pyroclastic soil deposits was discretised into layers 100 mm thick, for five horizons, from the top to bottom. The total thickness of ash-fall pyroclastic soil cover as well as that of each soil horizon was set on the basis of empirical relationships, which were found existing with the slope angle (Tufano et al. 2016;De Vita et al. 2018).
For the maximum thickness and the most complete stratigraphic series, five soil horizons were adopted: (1) B horizon (1100 mm); (2) C horizon (1100 mm); (3) Bb horizon (1700 mm); (4) C horizon (1100 mm); (5) Bb basal horizon (600 mm). Instead, for thinner total soil thickness, soil horizons were set with a reduced thickness or, alternatively, as disappearing after their pinching out (Fig. 4). This stratigraphic setting allows one to consider the vertical heterogeneity in terms of texture and bedrock depth and to investigate vertical distribution and spatial variation of soil water pressure head under specific rainfall events. In fact, several studies have highlighted how, under these conditions, soil water pressure head increases downstream along the more permeable pumice soil horizons, without the formation of an extensively continuous perched water table (e.g. Frattini et al. 2004) and undergoes an abrupt increase where these horizons pinch out. Soil thickness maps and vertical stratigraphic discretisation were implemented in the model to emphasise the role of stratigraphy (e.g. downstream pinching out) and effects of morphological discontinuities (e.g. rocky cliffs) on thickness, stratigraphic setting and hydrological response of ash-fall pyroclastic soil cover.
Three observation points were identified at upslope positions, going from a greater to smaller thickness (Fig. 4), which were used as control points for visualising plot scale evolutions of modelled variables, such as soil water content and soil water pressure head. Soil hydraulic conductivity functions and SWRCs, characterised by field and laboratory tests in previous research for each pyroclastic soil horizon , were implemented in the GEOtop in terms of horizontal K h (mm s −1 ) and vertical K z (mm/s) hydraulic conductivity, van Genuchten parameters n (-) and α (m −1 ), and residual and saturated water content θ r and θ s ( Table 1). The time step of simulation was set equal to 30 s, which was considered as a good compromise between computational time and numerical stability.
Rainfall input was derived from physics-based rainfall thresholds, specifically developed by Napolitano et al. (2016) for the wet season (winter). In particular, we considered different constant rainfall intensities, corresponding to 2.5, 5.0, 10.0, 20.0 and 40.0 mm h −1 for a duration, respectively, of 78, 41, 22, 12 and 6 h. Finally, the initial hydrological conditions were set according to values measured in the year 2005 on the analysed hillslope (Fusco et al. 2013(Fusco et al. , 2017. The soil strength parameters used in the slope stability analysis were derived from direct shear tests determined for each soil horizon   (Table 2).
Specifically, effective cohesion (c′) and friction angle (ϕ′) were respectively set to the 5th and 50th (median value) percentiles of the dataset to cope with the significant variability of results of laboratory tests. In particular, friction angle strictly depends on the existence of coarser pyroclasts (pumiceous lapilli) dispersed in a sandy-silty matrix (coarse and fine ash) of the soil samples, whereas the effective cohesion (c′) depends on the effects due to crushing of coarser pyroclasts (cohesion intercept). To eliminate the effect of the latter, a lower percentile was chosen for c′.
To complete the assessment of the model results, soil water pressure head and water content maps were used to implement the slope stability analysis in term of probability of failure maps.

Results
Model calibration and analysis of soil water pressure head distribution GEOtop-FS 2.0 provided maps and time series of soil moisture and soil water pressure head at different depths of the computational domain (Formetta and Capparelli 2019). To calibrate the model, different initial pressure head values were set on the basis of soil water pressure head maps resulting from the first time-step   Napolitano et al. (2016) determined an extremely rapid saturation of pyroclastic soil covering. Therefore, to obtain the optimal initial conditions, we set up initial soil water pressure head conditions equal to −0.50 m for B horizon, 0.00 m (saturation conditions) for C and Bb horizons, and −1.0 m for Cb and Bb basal horizons. Next, we set a period of 9.0 h of free lateral and bottom drainage, preceding the rainfall input. Such conditions determined mean values of soil water pressure head, at the beginning of the simulated rainfall event, equal to −1.10 m for B horizon, −2.30 m for C horizon, −1.60 m for Bb horizon, −0.70 m for Cb and −0.30 for Bb basal soil horizons (Table 3). These values were found satisfactorily matching those resulting from field hydrological monitoring (Fusco et al. 2017) and in particular, those recorded in the rainy periods (typically spanning from November to March).
To analyse the role of the stratigraphic setting on the hydrological regime, time series of soil water pressure head for the three observation points were plotted considering a rainfall intensity of 40 mm h −1 (Fig. 5).
Values assumed by soil water pressure head and soil water content were found to progressively increase as the soil thickness decreases. For the three observation points, characterised respectively by two for the first and five horizons for second and third, the deepest soil horizons reach the highest soil water pressure head values (−0.30 to −0.10 m), whereas the surficial soil horizons suffer mostly from the rainfall effects, especially in observation point 1. In this case, a nearsaturation condition of the entire ash-fall pyroclastic soil cover is reached and a lateral downstream throughflow (Kirkby 1978), with a prevailing component parallel to the slope, occurs. To visualise such a phenomenon across the whole computational mountain catchment, Fig. 6 shows the soil water pressure head variation at a constant depth of 1.1 m and for the five different rainfall conditions (2.5 mm h −1 for 72 h; 5 mm h −1 for 41 h; 10 mm h −1 for 22 h; 20 mm h −1 for 12 h; 40 mm h −1 for 6 h). Specifically, for rainfall intensities of 10, 20 and 40 mm h −1 , soil water pore pressure head increases particularly in correspondence to zones characterised by a downstream pinching out of the ash-fall pyroclastic soil covers, close to the outcropping bedrock. In these zones, Table 3 Mean soil water pressure head (and SD) for each soil horizon obtained after free drainage in absence of simulated rainfall  soil water pressure head increases until reaching near-saturation conditions and making these zones the most susceptible to landslide triggering. In particular, where the ash-fall pyroclastic soil cover thickness is thicker, soil water pore pressure head shows a damped and delayed response to the rainfall input. Instead, as total soil thickness reduces, the wetting front tends to involve the whole thickness with an extension strictly connected with rainfall conditions. Another important result achieved from simulations is the recognition of the extension of the area involved in the increase of soil water pressure head, which is limited to few metres upstream of the pinching out of ash-fall pyroclastic soil horizons.

Probability of failure
GEOtop-FS 2.0 also provided, in addition to soil moisture and soil water pressure head, maps of the probability of failure at different depths of the computational domain, which were calculated by maps of soil pressure head derived by the previous computational phase and the soil shear strength parameters for each soil horizon (Table 2). In this phase, stratigraphic settings were found to control the vertical location at which slope instability occurs, providing additional information in the rainfall-induced debris-slide hazard assessment. In this case, we analysed the failure probability maps at various depths, paying particular attention to those located at the interface between two horizons of ash-fall pyroclastic soil cover. In particular, the maps considered, as the most interesting from a slope stability point of view, are those computed for the depths between 1.10 ( Fig. 7) and 2.20 m (Fig. 8) of the ash-fall pyroclastic soil cover, corresponding to the contacts between B-C and C-Bb soil horizons, respectively.
To carry out a sensitivity analysis, we considered as unstable pixels those characterised by a failure probability value between the interval of 0.75 and 1.00, which is corresponding to 75-100%.
Depending on the depth, the extension of ash-fall pyroclastic soil cover involved in slope instability changes and specifically is reduced in Fig. 8, where the bedrock is more widespread in the mountain catchment, respect to Fig. 7. At 2.2 m, the unstable pixels are confined in areas of ash-fall pyroclastic cover. However, in both cases, the unstable pixels increase as rainfall intensity increases. Table 4 shows the number of unstable pixels and their ratio in comparison to the total number of pixels of the whole computational catchment. In particular, for rainfall intensity of 2.5 mm h −1 , unstable pixels are 49 at 1.1 m and 0 at 2.2 m; for 5 mm h −1 , we obtained respectively 19 and 0 pixels; for 10 mm h −1 , 76 and 0 pixels, respectively; for 20 mm h −1 , 300 and 547 pixels; and, finally, for 40 mm h −1 , 2182 pixels at 1.1 m and 123 pixels at 2.2 m. The total number of pixels making up computational catchment (18,028 at 1.1 m and 14,428 at 2.2 m) allows one to determine the unstable percentage in the entire computational catchment and to estimate the most critical depth to formation of rupture surfaces.
The instability percentage is generally equivalent for the two depths and for the lowest simulated rainfall intensities of 2.5 and 5 mm h −1 , but for intensity values of 10 and 20 mm h −1 , a greater percentage was recorded at 2.2 m of depth. The surface of rupture was found not always at the contact between the same horizons, but at variable depths depending on local factors and rainfall intensity as well. However, the small number of unstable pixels in comparison to that of the whole computational catchment proves very narrow triggering areas, which is spatially localised in close proximity to morphological drops in total thickness of ash-fall pyroclastic soil cover and pinching out of pyroclastic soil horizons, corresponding to local increases of slope angle.

Discussion
The modelling of transient infiltration and seepage in ash-fall pyroclastic soil coverings by taking into account their spatially variable thickness and stratigraphic setting is demonstrated to provide more detailed information about slope hydrological response under simulated rainfall events and antecedent soil hydrological conditions. The first results obtained by the calibration phase are represented by soil water pressure head distribution maps for different rainfall inputs and at different depths. Considering that at the depth of 1.1 m, where a high contrast in hydraulic conductivity between current soil (A + B soil horizons) and pumiceous lapilli (C horizon) exists (Fig. 6), soil water pressure head varies from −4.0 m to the near-saturation condition. In particular, the lowest value was recorded for rainfall input of 2.5 mm h −1 in the area characterised by the maximum thickness of ash-fall pyroclastic soil cover due to a prevalence of a vertical water flow. Instead, highest soil water pressure head conditions result as being localised in narrow zones (5-10 m wide) where the formation of lateral throughflow leads to a local increase of the soil water pressure head and water content up to near saturation conditions, as already demonstrated by the 2D hydrological models developed by De  and by distributed modelling of Frattini et al. (2004).
The ash-fall pyroclastic soil cover thickness (a) and soil water pressure head (b-c) profiles developed in Fig. 9 for a rainfall input of 40 mm h −1 show clearly how in correspondence to abrupt reductions in thickness of the ash-fall pyroclastic soil mantle, closely to the outcropping bedrock, the soil water pressure head tends to reach higher values. This phenomenon is clearly recognisable both at initial conditions of simulation (b) derived Fig. 9 Ash-fall pyroclastic soil cover thickness (a) and pressure head maps at start (b) and end (c) of hydrological modelling at 1.1 m of depth. On the right, profiles of stratigraphic thickness of ash-fall pyroclastic soil cover and of soil water pressure head. The increase of soil water pressure head can be observed as well correlated to variation of ash-fall soil cover thickness by redistribution of water content after a free drainage period and at the end of rainfall period (c) after 6 h. This trend is in agreement with results of some studies (Celico and Guadagno 1998;Crosta and Dal Negro 2003;Di Crescenzo and Santo 2005;Guadagno et al. 2005;Cascini et al. 2008;Napolitano et al. 2016), thus affirming the influence of topographical factors on landslide susceptibility for ash-fall pyroclastic deposits on steep slopes. Moreover, due to the control of slope angle on total thickness of soil coverings and the related stratigraphic setting, the obtained modelling results can be conceived as representing a hydro-geomorphological model (Sidle and Onda 2004;Fusco et al. 2017) which indicates the local increases of slope angle, leading to knickpoints or outcrops of rocky scarps, as the slope sectors most susceptible to landslide initiation.
Outcomes of modelling deliver very important information concerning hydrological and slope stability mechanisms leading to landslide initiation. The first is that the stratigraphic setting not only controls the spatial distribution of soil water pressure head, but also the vertical location at which instability can occur, thus providing additional information that can be used to evaluate the model performance. For such a purpose, the slope stability analysis was conducted in GEOtop-FS 2.0 and the results, derived by probability of failure maps, were evaluated at two depths: 1.1 and 2.2 m corresponding respectively to the contact between (A + B)-C and C-Bb soil horizons. Results shown in Figs. 7 and 8 indicate a variable behaviour. For 2.5 mm h −1 , there are 49 unstable pixels at 1.1 m but no pixels at 2.2 m. For 5 mm h −1 , the slope instability percentage is the same, but for 10 and 20 mm h −1 , the number of unstable pixels is greater for 2.2 m, which is contrasting to what happens for 40 mm h −1 . Also interesting is the instability distribution; at first layer interface, the unstable pixels are located in correspondence to pinching out of horizons and to bedrock outcrop. Such behaviour is promoted by the high contrasts of permeability and by a temporary impermeable barrier established in the paleosols (Mancarella et al. 2012). If failure does not occur during this soil water pressure transition, ash saturates and starts to slowly allow drainage, thus enabling the propagation of the infiltration front to deeper layers (Lizárraga and Buscarnera 2019). This is visible going from 1.1 to 2.2 m, where on the lower part of the maps the ash-fall pyroclastic soil cover is very discontinuous.

Conclusions
For the assessment of hazard to initiation of shallow landslides, we developed a physics-based model aimed at investigating the effect of spatial variability of soil thickness on slope hydrological response and next on the related probability of slope instability. This approach was tested in a sample area of the Sarno Mountains (Campania, southern Italy) whose slopes are mantled by ash-fall pyroclastic soils erupted mainly by the Somma-Vesuvius volcano with varying thickness and stratigraphic settings depending on the slope angle. The physics-based model was implemented in GEOtop-FS 2.0 and based on the following datasets: (1) thickness map of ash-fall pyroclastic soils, obtained by an empirical model of the spatial distribution of ash-fall pyroclastic soil thickness (De Vita et al. 2006aFusco et al. 2017 GEOtop-FS 2.0 was found able to solve the 3D Richards' equation and to consider a variable soil thickness as well as a complex stratigraphic setting to investigate their effect on hydrological regime of ash-fall pyroclastic soil cover. The computational catchment was calibrated in a small area falling within the portion of the territory affected by the event of May 1998. In particular, the soil water pressure head maps obtained highlight that the most susceptible areas to triggering of initial debris slides occur closely to morphological discontinuities, where the ash-fall pyroclastic soil thickness reduces until the outcrop of bedrock. In addition, the failure probability maps gave further information about the depth to which the slope instability can occur.
The results achieved in this research can be conceived as relevant advances in the topic of the assessment of hazard to initiation of rainfall-induced shallow landslides because of the identification of the relevant role played by spatial variability of soil thickness and stratigraphic settings.
In the context of early warning systems for rainfall-induced landslides, the application of the proposed approach to larger areas presents some limitations linked to the need of (1) very detailed input data (such as in situ measured soil-hydraulic parameters, soil thickness maps, etc.) and (2) high computational power. However, these limitations are somehow intrinsic to all physics-based approaches. In contrast, empirical, lumped statistical approaches continue to have a strong relevance and a practical value for landslide hazard assessment (e.g. early warning system), being less demanding both in terms of input data and computational power, even if limited by significant uncertainties due to incompleteness of landslide inventory considered, unrepresentativeness of rainfall measurements and unknown antecedent soil hydrological conditions. Finally, the most important result reached in this study is a more comprehensive understanding of the slope hydrological response controlled by spatially variable soil thickness and stratigraphic setting, indicating a spatial link among slope sectors most susceptible to the landsliding onset and rainfall conditions leading to landslide triggering. Therefore, the proposed approaches can be used as a tool to understand and reduce the uncertainty that potentially affects rainfall thresholds of shallow landslides, estimated by empirical or physics-based methods, and conceived to be potentially exportable to other slope environments for which a spatial modelling of soil thickness would be possible.

Acknowledgements
Authors acknowledge the funding of the PhD Program of Dipartimento di Scienze della Terra, dell'Ambiente e delle Risorse (DiSTAR), Università di Napoli Federico II (32 nd cycle).

Funding
Open access funding provided by Università degli Studi di Napoli Federico II within the CRUI-CARE Agreement.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's

Original Paper
Landslides 18 & (2021) Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/ 4.0/.