The Role of Cover Thickness in the Rainfall-Induced Landslides of Nocera Inferiore 2005

: In the context of rainfall-induced landslides involving pyroclastic soils, the present work analyzes the inﬂuence of cover thickness on slope stability conditions. To this aim, the slope failure that occurred in Nocera Inferiore (4th March 2005) is selected as a reference test case, providing the actual weather forcing history that preceded the event, the hydraulic characterization of the soil involved, and the lowermost boundary condition (variously fractured calcareous bedrock underlying the cover). By maintaining unchanged soil hydraulic properties, the relationship between domain thickness, initial soil suction distribution, and slope instability induced by critical rainfall is investigated by numerical analyses. These refer to a rigid unsaturated domain subject to one dimensional ﬂow conditions under the e ﬀ ects of incoming (precipitation) and outcoming (evaporation) ﬂuxes applied at the uppermost boundary. The main outcomes indicate that critical event duration increases signiﬁcantly with increasing the domain thickness. This relationship is strongly inﬂuenced by initial suction distribution. A linear relationship results for soil suction that is assumed to be constant at the beginning of the critical event. However, this relationship is quadratic if, by simulating the actual antecedent meteorological conditions, suction at the beginning of the critical event is the main function of the domain thickness. Additional numerical analyses were carried out to characterize the inﬂuence of a di ﬀ erent lowermost boundary condition. Outcomes indicate that, for the same thickness, critical duration is substantially longer if the cover contact is with the same material as that of the cover.


Introduction
In recent decades, rainfall-induced landslides have caused increasing damage and causalities globally [1]. Practitioners and scientists have consequently attempted to characterize the main mechanisms involved and to develop mitigation strategies. As communities and decision makers have increased their interest in environmentally compliant solutions, an increasing number of investigations have been focused on improving the efficiency of nature-based approaches [2] and non-structural approaches (early warning systems) [3][4][5].
Among these phenomena, debris avalanches and debris flows occurring in silty volcanic cover [6][7][8] have been greatly considered, as they have increased in frequency and magnitude globally during the past three decades. The Campania Region (Italy) is a much-studied area, due to the extensive presence of volcanic cover that originated from Vesuvius and Campi Flegrei eruptions over tens of thousands of years [9]. A typical susceptible stratigraphy comprises cover of a few meters placed over flysch deposit, tuff, or calcareous bedrock, with a slope of 30-45 degrees, and characterized by high porosity (70% or above). Sliding rates of tens of kilometers per second have sometimes characterized the kinematics of hundreds of thousands of cubic meters of debris, resulting in several catastrophic events.
For these phenomena, several recently published works have focused on identifying the factors affecting the triggering stage and the post-failure kinematics of the sliding mass. Many studies have attempted to detect susceptible slopes in terms of geomorphological indicators, particularly stratigraphy, slope tilting, and the thickness of the volcanic cover overlying tuff, calcareous, and flysch bedrocks [10][11][12][13]. Several works illustrate the intrinsic and state properties of the involved materials [14]. Other works characterize mechanical soil behavior in unsaturated or saturated conditions, in terms of strength parameters [6,[15][16][17], compressibility [14], collapsibility upon wetting [18,19], state conditions inducing liquefaction susceptibility [20,21], and soil hydraulic properties [22][23][24]. Still others illustrate and interpret results by field or physical model monitoring [25,26], with the aim to assess the impact of suction and water content evolution on pyroclastic covers under different weather scenarios and, in addition, how these evolutions affect slope stability conditions [25,[27][28][29][30][31]. In the same field, some works attempt to gain insight into the complex topics of the soil-atmosphere interaction and boundary conditions generated at the uppermost surface by weather conditions in terms of rainfall infiltration and evapotranspiration [32][33][34]. Some studies have even focused on the effect of the lowermost boundary conditions on different contact types standing at the cover bottom [32,35].
Despite its importance, the role that cover thickness plays in stability conditions has received little attention by researchers in the field of the stability of volcanic covers. Some theoretical investigations have been made about covers made of coarser and more pervious soils than those in hand, in which it was found that cover thickness does not significantly affect stability conditions if compared with soil hydraulic and mechanical properties [36] or that a thicker cover worsens stability conditions [37,38]. In several geological works concerning silty pyroclastic covers, the cover thickness is adopted as a key factor to assess susceptibility to debris avalanches [39]. This susceptibility is related to the existence of a sloping mass that is sufficiently greater than a critical amount to cause damage in the case of sliding. Thus, a thicker cover in these a works seems implicitly related to a worsening of stability conditions. This assumption is certainly true if the strength contribution provided by suction, s u , is, for instance, assumed constant throughout the cover and not affected by the cover thickness.
Indeed, considering, for example, a cohesionless sloping cover of thickness s z , relatively small compared to the slope length, lying on an infinitely resistant bedrock, inclined by ε, the expression of safety factor (FS) under partially saturated conditions delivered by limit analysis of an infinite slope [40] can be written as: FS = tanϕ tanε + s u χtanϕ γ s z senε cosε (1) where ϕ represents the angle of internal friction, γ the soil unit weight, and χ a parameter, ranging from 0 to 1, representing soil wetness conditions. Under the assumption of constant suction within the cover, a virtual increase in cover thickness would determine an increase in driving forces (denominator of second term) with resistance forces held constant. On the other hand, monitoring and theoretical findings indicate that covers experience a continuous change in outcoming and incoming water fluxes that systematically lead to transient conditions and suction changes, with the latter strongly affected by the cover thickness. In this way, stability conditions may even benefit from a thicker cover, as it facilitates the spreading of the same rainfall volume within a larger soil volume, limiting the average rainfall-induced suction drops. In the second term of Equation (1), the numerator might increase by more than the denominator, resulting in an increase in FS.
The aim of the current research is to provide quantitative insights about the effects of cover thickness on stability conditions for slopes subjected to rainfall and other meteorological conditions. The study exploits a previous theoretical investigation of the hydrological evolution that resulted in the Nocera Inferiore 2005 landslide [8], which occurred on the slopes of Mount Albino, situated in the Lattari Mountains (Figure 1).
Such a case study is used as a reference regarding material properties, evolution of antecedent meteorological conditions, and lowermost boundary conditions, to carry out the investigation within a realistic context. The effects of a hypothetical increase or decrease in cover thickness from the actual value (2 m) measured at the triggering area are investigated by analyses carried out by means of numerical simulations. These are aimed at interpreting the hydrological response of domains of different thicknesses to weather forcing histories able to potentially induce slope failure events. The numerical analyses are carried out by Finite Element Method (FEM) resolution of Richards' equation. The aim is to characterize the relationship between cover thickness and a virtual rainfall event inducing instability (critical event). The influence of antecedent meteorological conditions on this relationship is also characterized.
The paper initially describes the relevant geomorphological and climatological contexts. It then describes soil properties of covers involved in the test case area. Subsequently, it provides details about theoretical approaches adopted to study the influence of the cover thickness. Finally, it illustrates and discusses results yielded by the analyses.

Geomorphological Context
The study area is the district of Nocera Inferiore, located in Southern Italy, 15 km north-west of the town of Salerno. The southern part of the town is bordered by the Lattari Mountains, where sloping pyroclastic covers have often been affected by debris avalanches.
Several works have attempted estimations of the thickness distribution of pyroclastic covers over fractured carbonaceous slopes of the Lattari Mountains. Matano et al. [10] derived pyroclastic cover thicknesses throughout the Mount Albino slopes using in situ investigations, a geostatistical interpolation technique, and a heuristic approach. In steeper and craggier slopes, cover thicknesses were found to vary greatly between 1.5 and 4.5 m, while in the foothill zones they were found to be larger. Lucà et al. [12] adopted a stochastic approach to build a map of cover thickness for the Mount Pendolo slopes (Gragnano), located around 10 km from Mount Albino slopes. Their find were in a range fully consistent with that obtained for Mount Albino slopes. In several investigations, cover thicknesses are estimated at a scale (regional) larger than that of a single slope [11,13]. Del Soldato et al. [13] compared cover thicknesses evaluated using different approaches for areas surrounding the volcanic centers of Somma-Vesuvius and Campi Flegrei, and including the slopes mentioned above. For the Lattari Mountains, the work confirms the thickness range noted earlier. On gentle slopes, the stratigraphy includes, moving from bedrock to ground surface, a basal paleosol of silty sand or sandy silt about 1 m thick, resulting from ancient eruptions; a layer of pumice (coarse sand and gravels) with a low degree of weathering up to 1 m thick; and a weathered and pedogenized volcanic ash layer, up to 1 m thick, made of silty sands or sandy silt [41][42][43]. For the same area, a similar layering was also indicated by Di Crescenzo & Santo [44] for those slopes known to be previously affected by debris avalanches.
For steeper slopes, characterized by slope angles higher than 30 • -35 • , the layer of pumice tends to vanish. In this way, the stratigraphy of the pyroclastic cover, still made of two layers from a geological point of view, becomes almost homogeneous in hydraulic properties. Such homogeneity is observed for the landslide area of the Nocera Inferiore landslide at the apical zone, where the phenomenon was thought to have been triggered [45]. The same homogeneity was also observed in numerous other areas of the Lattari Mountains affected by debris avalanches in the past, such as the San Pantaleone hill (phenomena occurred in 1961, 1972, and 1997) and Monte Pendolo (phenomenon occurred in 1971).
For these silty covers, the finest component of ash soils (silt with some clay) is generally non-plastic due to its volcanic origin. Soil porosity can be as high as 70% with peaks attaining 80%.

Climatological Context and Impact on Landslides
Within the reference area, the mean annual rainfall was about 1140 mm for the period 1981-2010, with mean maximum daily rainfall of 92 mm. Over the same time span, the mean temperature was 18 °C, with maximum mean value of about 22.8 °C and minimum of 13.4 °C. According to the Köppen-Geiger climate classification [46], the investigated site falls into the Csa group-that is, a humid temperate climate with dry summers and the warmest month characterized by temperatures higher than 22 °C. Two different climatological regimes can be recognized each year: the first, from October to April, features significant rainfall cumulative values (about 70% of the total annual value) and lower temperatures, and the latter, from May to September, features warmer temperatures and low cumulative rainfall values.
Landslides affect covers around two meters thick. The debris avalanches that occurred in the past 60 years were induced by rainfall amounts exceeding 500 mm over the antecedent period considered in the same hydrological year, and a critical rainfall event that was exceptional for its persistency (as long as 15 hours) rather than its intensity or accumulated amount. Rianna et al. [47] attempted to provide a physical interpretation of the rainfall patterns generating slope instability in these soils.
The last and most significant landslide occurred on 4th March 2005. A weather station located very close to the landslide area provided the evolution of meteorological conditions that resulted in the phenomenon. Hourly values of precipitation, air relative humidity, and air temperature are available over the timespan 1st January 1998-1st November 2008. Pagano et al. [8] provided a first numerical interpretation of the case based on an account of rainfall evolution. Reder et al. [34] and Pagano et al. [48] re-analyzed the case following different approaches that also accounted for evaporation and evapotranspiration evolutions. This case is quite representative of many other debris avalanches that have affected the Lattari Mountains in terms of involved soils, antecedent meteorological evolution (800 mm since 1st September 2004), and duration of the critical event (16 hours). The rainfall amount of the critical event was 143 mm. After the landslide, a triangularly

Climatological Context and Impact on Landslides
Within the reference area, the mean annual rainfall was about 1140 mm for the period 1981-2010, with mean maximum daily rainfall of 92 mm. Over the same time span, the mean temperature was 18 • C, with maximum mean value of about 22.8 • C and minimum of 13.4 • C. According to the Köppen-Geiger climate classification [46], the investigated site falls into the Csa group-that is, a humid temperate climate with dry summers and the warmest month characterized by temperatures higher than 22 • C. Two different climatological regimes can be recognized each year: the first, from October to April, features significant rainfall cumulative values (about 70% of the total annual value) and lower temperatures, and the latter, from May to September, features warmer temperatures and low cumulative rainfall values.
Landslides affect covers around two meters thick. The debris avalanches that occurred in the past 60 years were induced by rainfall amounts exceeding 500 mm over the antecedent period considered in the same hydrological year, and a critical rainfall event that was exceptional for its persistency (as long as 15 h) rather than its intensity or accumulated amount. Rianna et al. [47] attempted to provide a physical interpretation of the rainfall patterns generating slope instability in these soils.
The last and most significant landslide occurred on 4th March 2005. A weather station located very close to the landslide area provided the evolution of meteorological conditions that resulted in the phenomenon. Hourly values of precipitation, air relative humidity, and air temperature are available over the timespan 1st January 1998-1st November 2008. Pagano et al. [8] provided a first numerical interpretation of the case based on an account of rainfall evolution. Reder et al. [34] and Pagano et al. [48] re-analyzed the case following different approaches that also accounted for evaporation and evapotranspiration evolutions. This case is quite representative of many other debris avalanches that have affected the Lattari Mountains in terms of involved soils, antecedent meteorological evolution (800 mm since 1st September 2004), and duration of the critical event (16 h). The rainfall amount of the critical event was 143 mm. After the landslide, a triangularly shaped area was clearly discernible. It spread over a 36 • inclined open slope from the apical zone downward, mobilizing a soil mass of 33,000 m 3 . The kinematic was that of a debris avalanche. At the apical zone, the slope angle approached the soil friction angle of 39 • [49] and the cover was 2 m thick. The bedrock consisted of highly fractured limestone located at a depth ranging from 1 to 2 m, approaching the maximum values at the apical zone. Figure 2 shows the evolution of the antecedent rainfall and reference evaporation from bare soil. The latter was obtained by interpreting observations of air temperature and relative humidity over the time span between the start of the hydrological year and the landslide time, according to the approach proposed by Allen et al. [50].
Geosciences 2020, 10, x FOR PEER REVIEW 5 of 21 approaching the maximum values at the apical zone. Figure 2 shows the evolution of the antecedent rainfall and reference evaporation from bare soil. The latter was obtained by interpreting observations of air temperature and relative humidity over the time span between the start of the hydrological year and the landslide time, according to the approach proposed by Allen et al. [50].

Materials
The material involved in the Nocera Inferiore landslide of 2005 was a non-plastic sandy-silt volcanic soil ( Figure 3) having a porosity as high as 70%. At such a porosity, it was tested in a physical model of lysimeter to retrieve the soil hydraulic properties. An extensive description of the physical model and experimental data may be found in Rianna et al. [33].
The Soil Water Characteristic Curve (SWCC) was assessed ( Figure 4a) directly from measurements taken within the lysimeter as a best fit function of water content-suction points recorded at all depths by Time Domain Reflectometry probes (TDRs)and jet-fill tensiometers. The SWCC adopted interpolates available experimental data up to suction values of around 90 kPa and extrapolates them to much higher suction levels, adopting as an asymptote the lowest volumetric water content measured locally by TDRs.
The hydraulic conductivity function ( Figure 4b) resulted from the back analysis of the monitored evolution of suction and volumetric water content [34]. The hydraulic conductivity of the material drops by three orders of magnitude (from 1 × 10 −6 to 1 × 10 −9 m/s) for suction increasing from 0 to 100 kPa.

Materials
The material involved in the Nocera Inferiore landslide of 2005 was a non-plastic sandy-silt volcanic soil ( Figure 3) having a porosity as high as 70%. At such a porosity, it was tested in a physical model of lysimeter to retrieve the soil hydraulic properties. An extensive description of the physical model and experimental data may be found in Rianna et al. [33].  The Soil Water Characteristic Curve (SWCC) was assessed ( Figure 4a) directly from measurements taken within the lysimeter as a best fit function of water content-suction points recorded at all depths by Time Domain Reflectometry probes (TDRs)and jet-fill tensiometers. The SWCC adopted interpolates available experimental data up to suction values of around 90 kPa and extrapolates them to much higher suction levels, adopting as an asymptote the lowest volumetric water content measured locally by TDRs.
The hydraulic conductivity function ( Figure 4b) resulted from the back analysis of the monitored evolution of suction and volumetric water content [34]. The hydraulic conductivity of the material drops by three orders of magnitude (from 1 × 10 −6 to 1 × 10 −9 m/s) for suction increasing from 0 to 100 kPa.

Stability Conditions
Characterizing the critical event requires that state conditions determining susceptibility to debris avalanches are preliminary defined. The concomitant presence of unconventional strength mechanisms acting within a pyroclastic cover makes accomplishing this task difficult.
One such mechanism is the strength contribution provided by abundant vegetation standing on these slopes [43,45]. Dense root systems of different types could cross the cover down to penetrate the fractured bedrock, thus tying the cover to bedrock. Soil suction supplies the action normal to the root surface needed to activate soil-root friction (weight-induced stresses are instead very low), allowing roots to work as reinforcements within the cover. The root system may also act as a traverse reaction and resistance against upstream-downstream oriented movements. The strength contribution is, in this case, based on soil passive strength, which is still related to soil suction levels throughout the cover.
Another factor affecting soil strength, in this case reducing it, is the possibility for spontaneous liquefaction upon wetting. In these soils, liquefiable conditions exist when full saturation is approached, and material porosity is above 70%. These state conditions may generate triggering and propagation of liquefaction phenomena for slope angles higher than 30°, hence for values even much lower than material friction angle at the critical state (usually, about equal to 39°) [7,20].

Stability Conditions
Characterizing the critical event requires that state conditions determining susceptibility to debris avalanches are preliminary defined. The concomitant presence of unconventional strength mechanisms acting within a pyroclastic cover makes accomplishing this task difficult.
One such mechanism is the strength contribution provided by abundant vegetation standing on these slopes [43,45]. Dense root systems of different types could cross the cover down to penetrate the fractured bedrock, thus tying the cover to bedrock. Soil suction supplies the action normal to the root surface needed to activate soil-root friction (weight-induced stresses are instead very low), allowing roots to work as reinforcements within the cover. The root system may also act as a traverse reaction and resistance against upstream-downstream oriented movements. The strength contribution is, in this case, based on soil passive strength, which is still related to soil suction levels throughout the cover.
Another factor affecting soil strength, in this case reducing it, is the possibility for spontaneous liquefaction upon wetting. In these soils, liquefiable conditions exist when full saturation is approached, and material porosity is above 70%. These state conditions may generate triggering and propagation of liquefaction phenomena for slope angles higher than 30 • , hence for values even much lower than material friction angle at the critical state (usually, about equal to 39 • ) [7,20].
Such strength factors are difficult to quantify and consequently complicate the task of identifying which hydrological conditions induce an instability mechanism to become a rapid post-failure kinematics.
The numerical interpretation of the Nocera Inferiore landslide allowed the identification of a simple hydrological criterion associated with the slope failure. The analyses were carried out under one-dimensional flow conditions, at first by solving Richards' equation to simulate rainfall infiltration [8], and then solving Richards' equation to account for evaporation as a boundary phenomenon [34]. Finally, solving the equation according to the hydrothermal coupled approach developed by Wilson et al. [51], to model evaporation as an internal phenomenon [48]. Soil hydraulic behavior is characterized by SWCC and the conductivity function displayed in Section 2.3. Numerical results yielded by the different types of analyses of the actual meteorological evolution recorded during the landslide hydrological year, all identified the same hydrological singularity at the triggering time, which was that of suction almost vanishing throughout the cover. This should represent a hydrological state consistent with a cover prone to generate a debris avalanche due to two reasons: 1.
An extremely rapid kinematic was observed, which indicates unequivocally the occurrence of extensive liquefaction phenomena, in turn necessarily related to a state approaching full saturation and suction vanishing at any depth.

2.
A dense root system is present [43,45], which implies that an unavoidable requirement for sliding is the "shutdown" of all root strength mechanisms mentioned previously, and this occurs if suction vanishes throughout the cover.
Suction vanishing throughout the cover will be identified in the next section as the hydrological state relating to the occurrence of a debris avalanche.

Mathematical-Numerical Model
With the objective to translate a given meteorological evolution ( Figure 2) into a suction evolution within a continuum domain schematizing the cover, a simple model was used. Its main features are: One-dimensional water flux is assumed; this hypothesis should be realistic for the case of both homogeneous and heterogeneous covers, as pointed out by previous theoretical and experimental works [8,53,54].

•
Hourly rainfall is translated into hydraulic boundary conditions; a water flux condition equal to rainfall intensity is maintained at the uppermost surface if surface pore water pressure is less than zero; otherwise, a null pore water pressure condition is assumed. The latter is then maintained until entering water fluxes are less than the rainfall intensity, otherwise entering flux is again switched at the rainfall intensity. • Air temperature and air relative humidity are used to estimate actual evaporation flux, AE. For the sake of simplicity, all terms of internal evaporation are neglected, including that related to transpiration, although the presence of vegetation was implicitly taken into account in the criterion identifying cover instability. AE is then considered only acting at the boundary and is computed through several steps: evaluation of reference evaporation, ET 0 , quantified from weather data only (air temperature and relative humidity); computation of potential evaporation, PE, from bare soil, by applying the crop coefficient to ET 0 (PE = k crop × ET 0 ; k crop = 1.15, as found by Rianna et al., [33]); and characterization of a falling law for PE, consistent with available pore water at the uppermost surface and with a suction threshold, which may not exceed the value associated with thermodynamic equilibrium [51]. • Surface seepage is applied at the bottom boundary. For the case in hand, this condition corresponds to the effects induced by the intensely fractured bedrock standing at the bottom of the silty volcanic layer. This same condition simulates, in general, the presence of a capillary barrier at the bottom of a silty layer and may hence be used to simulate the lowermost boundary condition of a silty layer overlaying a pumice layer. This is the case of the gentle slopes of the Lattari Mountains (see Section 2.1) and of many other geomorphological contexts surrounding volcanic areas [35].
The numerical resolution of the equation was performed by exploiting the FEM code HYDRUS-1D [55]. This adopts a finite element solution of the partial differential flow equation for unsaturated soils.
Two main set of analyses were conducted to characterize the relationship between cover thickness and critical event. The first set, set A, characterizes this relationship in the simplest way, neglecting the antecedent meteorological evolution and applying directly on the domain the critical event started by a fixed initial suction (Figure 5a). The steps followed in this case are: The s uA0 trend was assumed to be constant at 20 kPa throughout the domain, corresponding to the average suction level typically observed after the rainfalls of the autumnal season. Four different suB0 were adopted. The first, suB0 = 70 kPa, considers the occurrence of preceding dry years. The second, suB0 = 20 kPa, considers the occurrence of a preceding very wet year. Between these extremes, intermediate initial conditions of suB0 = 35 kPa and suB0 = 50 kPa were also analyzed.
Two further set of analysis were additionally carried out to investigate the sensitivity of results to variations in saturated hydraulic conductivity (set C) and to a change in the antecedent rainfall window, obtained by backdating the potential triggering event (set D) The antecedent meteorological evolution selected for the analyses was that recorded over the hydrological year of the Nocera Inferiore 2005 landslide. This was to study the effects of a realistic meteorological evolution and, at the same time, one that is consistent with the occurrence of a landslide. It was applied since the beginning of the hydrological year (1st September 2004) until the day preceding that of the landslide (3rd March 2005 for set A, B, and C, and 1st January 2005 for set The critical event was defined of constant intensity, at 10 mm/h. This intensity is higher than the potential infiltration of the domain at full saturation (3.6 mm/h). In this way it generates a constant incoming flux when the domain approaches the near saturated state related to the occurrence of instability conditions, and should thus represent a realistic boundary condition for cases potentially leading to instability conditions. The critical event is hence characterized by a single parameter, i.e., the rainfall duration that induced the vanishing of the suction, consistent with findings by Rianna et al. [33].
Under the effects of a same antecedent meteorological evolution, different cover thicknesses are assumed to give rise to different suction distributions, in turn affecting the duration of the critical event.
To investigate these dynamics, a second set of analyses, set B, were carried out. For these analyses, the critical event was again applied with a constant intensity at 10 mm/h, after the simulation of the antecedent meteorological evolution (Figure 5b). In this case, the steps followed are: The s uB0 trend was assumed to be constant throughout the domain. This synthetizes the suction distribution resulting from the hydrological years preceding the year that was analyzed.
Four different s uB0 were adopted. The first, s uB0 = 70 kPa, considers the occurrence of preceding dry years. The second, s uB0 = 20 kPa, considers the occurrence of a preceding very wet year. Between these extremes, intermediate initial conditions of s uB0 = 35 kPa and s uB0 = 50 kPa were also analyzed.
Two further set of analysis were additionally carried out to investigate the sensitivity of results to variations in saturated hydraulic conductivity (set C) and to a change in the antecedent rainfall window, obtained by backdating the potential triggering event (set D) The antecedent meteorological evolution selected for the analyses was that recorded over the hydrological year of the Nocera Inferiore 2005 landslide. This was to study the effects of a realistic meteorological evolution and, at the same time, one that is consistent with the occurrence of a landslide. It was applied since the beginning of the hydrological year (1st September 2004) until the day preceding that of the landslide (3rd March 2005 for set A, B, and C, and 1st January 2005 for set D) (Figure 2). Figure 6 shows the final part (last three months) of the antecedent meteorological evolution (in blue bars) preceding the potential critical rainfall event (in red), applied on 4th March 2005.  Figure 7a plots the suction evolution induced by critical events. Suctions are averaged throughout domains having four different thicknesses si (s0.5 = 0.5 m, s1.5 = 1.5 m, s2.5 = 2.5 m, s4 = 4 m). All curves move from the same initial suction level (suA0 = 20 kPa). They then evolve by falling to reach null suction. As expected, it can be seen that the rate of suction reduction decreases, and critical event duration (di) increases, with increasing domain thickness. The relationship between critical duration and domain thickness (di-si) is linear (Figure 7b), ranging between d0.5 = 10 and d4 = 80 h.   All curves move from the same initial suction level (suA0 = 20 kPa). They then evolve by falling to reach null suction. As expected, it can be seen that the rate of suction reduction decreases, and critical event duration (di) increases, with increasing domain thickness. The relationship between critical duration and domain thickness (di-si) is linear (Figure 7b), ranging between d0.5 = 10 and d4 = 80 h.

Analysis of Set B
For the two domains s1.5 and s4, Figure 8 plots suction profiles predicted after the simulation of the antecedent meteorological evolution. The four profiles refer to the four initial suction values adopted: suB0 = 20 kPa; suB0 = 35 kPa; suB0 = 50 kPa and suB0 = 70 kPa.
All s1.5 profiles show a linear hydrostatic independent from suB0, with zero suction at the domain base and increasing suction upward. This is because the small size of the domain means that antecedent rainfall amounts are sufficient to increase water storage close to the maximum value, independently from the initial imbibition state. At the end of last precipitation of the antecedent period, suctions are reduced to very low levels, placed below the air entry value su* (about 10 kPa). When unsaturated conditions exist at the lowermost surface, the seepage surface acts as an

Analysis of Set B
For the two domains s 1.5 and s 4 , Figure 8 plots suction profiles predicted after the simulation of the antecedent meteorological evolution. The four profiles refer to the four initial suction values adopted: s uB0 = 20 kPa; s uB0 = 35 kPa; s uB0 = 50 kPa and s uB0 = 70 kPa.
All s 1.5 profiles show a linear hydrostatic independent from s uB0 , with zero suction at the domain base and increasing suction upward. This is because the small size of the domain means that antecedent rainfall amounts are sufficient to increase water storage close to the maximum value, independently from the initial imbibition state. At the end of last precipitation of the antecedent period, suctions are reduced to very low levels, placed below the air entry value s u* (about 10 kPa). When unsaturated conditions exist at the lowermost surface, the seepage surface acts as an impervious barrier, with accumulating precipitation water flowing downward and promoting the disappearance of suction. Once saturation has been reached, the surface seepage sets pore water pressure to zero and activates drainage. During the last three days of the antecedent period, no precipitation occurs (see Figure 6). During this period, suction tends to increase at any depth to re-establish hydraulic hydrostatic equilibrium, with water content falling due to base drainage. This process is rapid as suction changes, occurring below the air entry value (for s u < s u* ), require very small changes in volumetric water content.
At the end of the antecedent period, the thicker domain s 4 shows the suction profiles are markedly affected by initial suctions. Suction distributions obtained for the two initial applied limit conditions, s uB0 = 20 kPa and s uB0 = 70 kPa, differ from each other and both depart from the hydrostatic trend retuned by the thinner domain (s 1.5 ).
For s uB0 = 20, the analyses yield a profile less than that of the hydrostatic value, with zero suction at the base, while for s uB0 = 70, analyses yield an irregular profile, fluctuating around 30 kPa.
For s uB0 = 20kPa, a sequence like that described for the s 1.5 domain takes place. Because of the initial relevant imbibition state, the amount of antecedent rainfall pushes suction below the air entry value at the end of the last precipitation of the antecedent period. Subsequently, suction travels at different rates towards a hydraulic hydrostatic equilibrium over the last three dry days. At elevations showing s u < s u*, suction changes are rapid (as explained above), while at elevations showing s u > s u* , changes are slow, requiring significant water content changes. Suction hence easily reaches the hydrostatic value only there where s u < s u* . In contrast, the re-equilibrium process remains uncompleted where hydrostatic equilibrium requires s u > s u* . In this latter case, suction remains close to the air entry value, departing from the hydrostatic value of an amount decreasing with depth.
For s uB0 = 70kPa, suction cannot approach values close to s u* under the effects of antecedent rainfall, because of the initial relevant dry state of the domain. The domain remains, therefore, in an unsaturated state, which is characterized by low unsaturated hydraulic conductivity and suction changes related to significant water content changes, according to the active branch of the SWCC.
Both of these features promote re-equilibrium delay. Slow transient processes result in an irregular suction distribution far from that that of hydrostatic equilibrium.  Figure 9 plots suction evolution induced by critical events. Suction is averaged throughout each domain si, for all four initial conditions assumed. For domains thinner than s1.5, the trends are unaffected by initial conditions, consistent with the previous insights about the effects of antecedent rainfall, which induce a same (hydrostatic) suction distribution at the beginning of the critical event ( Figure 8). Domains thicker than s2 instead show an evolution clearly affected by initial conditions. At the end of the antecedent period, s2.5 yields suction values of su = 13 and 17 kPa, for suB0 = 20 and 70 kPa, respectively. For these same initial conditions, s4 yields su = 16 and 29 kPa.   Figure 10 plots the relationships between domain thickness, si, and duration of critical event, di, for the four different initial conditions taken into consideration. It is possible to notice that taking into account the antecedent meteorological evolution implies that the di-si relationships follows a quadratic trend in place of the linear trend previously obtained (see Figure 7) without considering the effects of the antecedent meteorological evolution. For the actual cover thickness of the case study (si = s2), the plots obtained for initial conditions of 20 and 35 kPa deliver a critical event duration (16 h) consistent with that recorded.   Figure 10 plots the relationships between domain thickness, s i , and duration of critical event, d i , for the four different initial conditions taken into consideration. It is possible to notice that taking into account the antecedent meteorological evolution implies that the d i -s i relationships follows a quadratic trend in place of the linear trend previously obtained (see Figure 7) without considering the effects of the antecedent meteorological evolution. For the actual cover thickness of the case study (s i = s 2 ), the plots obtained for initial conditions of 20 and 35 kPa deliver a critical event duration (16 h) consistent with that recorded.  Figure 10 plots the relationships between domain thickness, si, and duration of critical event, di, for the four different initial conditions taken into consideration. It is possible to notice that taking into account the antecedent meteorological evolution implies that the di-si relationships follows a quadratic trend in place of the linear trend previously obtained (see Figure 7) without considering the effects of the antecedent meteorological evolution. For the actual cover thickness of the case study (si = s2), the plots obtained for initial conditions of 20 and 35 kPa deliver a critical event duration (16 h) consistent with that recorded.

Analysis of Set C
To investigate the sensitivity of results to the main parameter characterizing the hydraulic response of the system, the hydraulic conductivity at saturation, its value (1 × 10 −6 m/s) was altered within a range of reasonable uncertainty. To this aim, the B analyses (Section 4.2) at s uB0 = 20 kPa were repeated twice by lowering the parameter to 0.8 × 10 −7 m/s and increasing it to 1.2 × 10 −6 m/s. These two alterations result ( Figure 11) in critical durations unaffected at all of the smallest thicknesses, and slight effects for largest thicknesses.

Analysis of Set C
To investigate the sensitivity of results to the main parameter characterizing the hydraulic response of the system, the hydraulic conductivity at saturation, its value (1 × 10 −6 m/s) was altered within a range of reasonable uncertainty. To this aim, the B analyses (Section 4.2) at suB0 = 20 kPa were repeated twice by lowering the parameter to 0.8 × 10 −7 m/s and increasing it to 1.2 × 10 −6 m/s. These two alterations result ( Figure 11) in critical durations unaffected at all of the smallest thicknesses, and slight effects for largest thicknesses.

Analysis of Set D
In order to better investigate the effects of antecedent rainfall, several analyses were carried out using the setting adopted in B analyses but backdating the starting time of the critical events to 1st January 2005, thus reducing the antecedent period of 62 days and the antecedent rainfall amount by about 500 mm (from about 1100 mm to 600 mm). The analyses were carried out for initial suctions (suB0) of 20 and 70 kPa. For suB0 = 20 kPa, antecedent rainfall amount is still enough to reduce suction levels ( Figure 12) to the same values previously computed on 4th March (Figure 8). Obviously, this entails an unvaried relationship between critical duration and domain thickness ( Figure 13). In contrast, for suB0 = 70 kPa, suction distributions ( Figure 12) differ significantly from those previously obtained (Figure 8), and this results in a di-si relationship characterized by critical durations that are considerably longer at any thickness value.

Analysis of Set D
In order to better investigate the effects of antecedent rainfall, several analyses were carried out using the setting adopted in B analyses but backdating the starting time of the critical events to 1st January 2005, thus reducing the antecedent period of 62 days and the antecedent rainfall amount by about 500 mm (from about 1100 mm to 600 mm). The analyses were carried out for initial suctions (s uB0 ) of 20 and 70 kPa. For s uB0 = 20 kPa, antecedent rainfall amount is still enough to reduce suction levels ( Figure 12) to the same values previously computed on 4th March ( Figure 8). Obviously, this entails an unvaried relationship between critical duration and domain thickness ( Figure 13). In contrast, for s uB0 = 70 kPa, suction distributions ( Figure 12) differ significantly from those previously obtained (Figure 8), and this results in a d i -s i relationship characterized by critical durations that are considerably longer at any thickness value.

Instability of Subdomains
A further elaboration permits accounting for the possibility that, for each domain si, instability could take place partially due to suction vanishing throughout a shallow subdomain having thickness si k < si. This is assumed to occur under the effects of a critical sub-event lasting di k < di. Figure 14 plots the results yielded by several analyses carried out to characterize this feature. Each point di-si previously obtained for suB0 = 70 kPa is now connected to a line built with points (di k -si k ), representing the subdomains of si and their critical sub-events. Regarding the domain s4, for which the critical event had been previously quantified in d4 = 124 h, Figure 14 shows that the subdomain

Instability of Subdomains
A further elaboration permits accounting for the possibility that, for each domain si, instability could take place partially due to suction vanishing throughout a shallow subdomain having thickness si k < si. This is assumed to occur under the effects of a critical sub-event lasting di k < di. Figure 14 plots the results yielded by several analyses carried out to characterize this feature. Each point di-si previously obtained for suB0 = 70 kPa is now connected to a line built with points (di k -si k ), representing the subdomains of si and their critical sub-events. Regarding the domain s4, for which the critical event had been previously quantified in d4 = 124 h, Figure 14 shows that the subdomain

Instability of Subdomains
A further elaboration permits accounting for the possibility that, for each domain s i , instability could take place partially due to suction vanishing throughout a shallow subdomain having thickness s i k < s i . This is assumed to occur under the effects of a critical sub-event lasting d i k < d i . Figure 14 plots the results yielded by several analyses carried out to characterize this feature. Each point d i -s i previously obtained for s uB0 = 70 kPa is now connected to a line built with points (d i k -s i k ), representing the subdomains of s i and their critical sub-events. Regarding the domain s 4 , for which the critical event had been previously quantified in d 4 = 124 h, Figure 14 shows that the subdomain s 4 0. 5

Discussion
Rainfall-induced debris avalanches in densely rooted pyroclastic covers are associated with extensive liquefaction phenomena, relating to suction nullifying throughout the cover. Under this assumption, the results shown above clearly indicate how the duration of the critical rainfall significantly increases with increasing cover thickness. This outcome overturns previous prevailing literature findings about the topic (see Section 1, [36][37][38]), based on a traditional assumption about the instability criterion and often obtained under a full saturation hypothesis for the analyzed domain. According to the instability criterion proposed in this work, numerical results clearly show that, starting from the same initial suction distribution, a thicker cover better mitigates the drop in suction induced by rainfall, distributing the same rainfall volume through a greater depth. This results in a longer critical duration required for a debris avalanche. In this case, the thickness-critical duration relationship is linear. However, this enhanced capability to spread rainfall volumes also reflects the way that antecedent meteorological conditions affect suction distribution at the beginning of the critical event. Mean suction levels are in this case higher for thicker domains, further enhancing the capability of a thicker cover to mitigate a critical event. This explains why the previously mentioned relationship becomes quadratic, amplifying the beneficial effect on stability exerted by a larger thickness.
The influence of thickness on the critical event is substantially modified if the hydraulic condition applied to the lowermost boundary is varied. For example, if the shallowest portion of the cover is viewed as a subdomain that could be subject to instability, the subdomain's lower boundary does not act as a discontinuity in the hydrological behavior and the capillary barrier results shift downward with respect to the subdomain base. This contact of the subdomain base with the same material entails a significant increase in the duration of this sub-critical event ( Figure 14). Such duration is essentially affected by the distance between depths of the subdomain base and the capillary barrier, because the latter generates water accumulation, promoting increments in water content above it that reduce moving upward.
For a given capillary barrier depth, the relationships between subdomain thickness and duration of subcritical events ( Figure 14) are linear. The relationships are also flattened, with subcritical events comparable with the critical event of the entire domain. This result may be explained through the analysis of the shapes of suction profiles and their evolution ( Figure 15). These profiles approach a zero-distribution forming a sort of bell-shaped distribution, with suction values that tend to nullify through a partial branch just before it nullifies throughout the entire domain.  The relationships referred to subdomains ( Figure 14) are shifted towards longer duration for deeper capillary barriers. This implies that subdomains having the same thicknesses, but different capillary barrier depths, are characterized by completely different subcritical events. For instance (Figure 14), the critical duration d 4 0.5 (97 h), is much longer than that of d 0.5 (10 h).
In Figure 10, the departing point of the curves marks the thickness (1.5-2 m) over which the antecedent period, which influences suction at the beginning of the critical event, becomes progressively longer than that analyzed (6 months). The analysis of stability conditions of covers thicker than 3 m should then go so far as to consider the influence of an antecedent meteorological window extended to some years. For these thicknesses, a landslide may be supposed to trigger under a critical event of realistic duration only if consecutive wet years have earlier dropped mean suction levels to very low values (about lower than 20 kPa).
A thicker domain is hence less hazardous, as it requires the occurrence of several consecutive very wet hydrological years to be affected by a landslide. On the other hand, it should also be noted that this lower likelihood of occurrence corresponds with a greater magnitude of the induced event.
It has been remarked that, for the steep slopes of the Lattari Mountains, pyroclastic cover frequently includes two layers that, although differing from a geological point of view, can be assumed to be very similar from a geotechnical point of view. Hence, the hypothesis of homogeneous medium is reliable, allowing for the results of the study to be extended to different contexts. The study may also be a reference for inhomogeneous stratigraphy if the shallowest silty layer overlays a pumice layer that acts as a seepage surface (capillary barrier) instead of the carbonaceous fractured bedrock. In such a frequent case, the results of the study may be referred to for the thickness of the shallowest silty layers.

Conclusions
The current study investigated the link between the domain thickness and critical event by parametrizing the thickness of a previously analyzed case study of a rainfall-induced landslide that occurred in Italy and involved a silty pyroclastic cover. Findings can be summarized as follows: 1.
if the antecedent meteorological evolution is neglected, starting from initial conditions in terms of domain imbibition typical of late autumnal periods, the relationship can be assumed to be linear; 2.
if the effects of an antecedent meteorological evolution, severe in terms of rainfall amount, is taken into account, the impact of the domain thickness is significantly more relevant and the relationship follows a quadratic trend; 3.
the influence of antecedent meteorological evolution implies that for relevant domain thicknesses the date of the critical event plays a crucial role; 4.
slight changes in hydraulic conductivity at saturation do not significantly affect results; 5.
the correct quantification of the critical event for domains thicker than 2 m requires that the antecedent meteorological evolution longer than one year is considered; 6.
under the same thickness, the duration of the critical event is strongly affected by the hydraulic boundary condition assumed at the lowermost boundary.
Funding: This research received no external funding.