From empirically to physically based early warning predictions of rainfall-induced landslides in silty volcanic soils: the Lattari Mountains case study

Abstract

The work proposes a procedure to build an early warning predictive tool to assess the occurrence of rainfall-induced landslides in silty volcanic covers. The procedure combines both an empirically and a physically based tool used sequentially: the former is designed to be calibrated using older, highly sized and coarser rainfall data, and the latter to interpret recent and finer weather data. Both approaches need to be informed by a common experimental reference summarising the rainfall history, the rainfall point, defined as the couple made of antecedent 4-month rainfall cumulative value (C_4m) and last-persistent event (C_PLE). The empirical approach aims to identify if, in the (C_4m–C_PLE) plane, the rainfall point falls in a ‘safe’ or ‘potentially unsafe’ zone where the two distinct regions are built by interpreting rainfall data associated or not with landslide events. In the physically based approach, evaporation and runoff are estimated to refine the assessment of ‘effective’ rainfall points. The resulting transformed rainfall point (C′_4 m, C′_PLE) is turned into a prediction of the suction level at the mid-depth assumed as a ‘reference’ for the entire cover. Such value is compared with a suction threshold empirically defined. Suction levels prediction is developed by computing in the C′_4 m–C′_PLE plane the iso-suction lines generated by several rainfall scenarios. The accuracy of the developed procedure is comparable with state-of-the-art literature or operational approaches, properly identifying landslide case events and minimising the number of false alarms. Furthermore, it can inform the preparedness stages more effectively, explicitly accounting for the antecedent slope wetness stage and how it could be far from the incipient slope failure conditions. The developed procedure takes into account the effects of evaporation and antecedent rainfalls that, in dry periods, lead to very dry conditions in the subsoil, making even significant rainfall events inconsequential. Conversely, other procedures already operating in LEWS or highly considered literature background overestimate the effects of rainfalls during dry periods. The developed procedure delivers a simple but robust way to derive landslide thresholds based on the interpretation of past rainfall histories. At the same time, literature methods often require sophisticated approaches to retrieve thresholds.

Introduction

Background

Sloping covers susceptible to rainfall-induced landslides and overstanding populated areas or linear infrastructures represent severe risks for communities and assets (Froude and Petley 2018; CRED 2020; Mazhin et al. 2021). In the last decades, an increasing interest has arisen in strategies based on Landslide Early Warning Systems (LEWSs), which are cheaper, more straightforward to set up, and socially accepted compared to grey structural solutions. From this perspective, well-established international initiatives to share knowledge and scale out good practices are feeding significant progress in understanding the dynamics that lead to landslide events (LANDAWARE, Calvello et al. 2020) .

The success of LEWSs depends mainly on the performance of the predictive models in terms of timing and accuracy (Pagano et al. 2010; Tofani et al. 2017). In general, they can be distinguished in territorial/geographical LEWS (Guzzetti et al. 2020; Piciullo et al. 2020), expecting to be operational at the regional, national and global scale and local LEWS (Lo-LEWS; Pecoraro et al. 2019), dealing with single landslides and limited geographical contexts.

Focusing the attention on the last ones, Pecoraro et al. (2019) framed the state-of-the-art of Lo-LEWS by investigating 29 operational systems (see Table 1). The analysis identifies three main types of alert parameters: precipitation features (P; e.g. intensity, duration, cumulative values), internal-to-slope hydrological variables (H; e.g. pore water pressure, water content, basin stream flows), soil deformation-related variables (S; e.g. displacements, strains, microseismicity, acoustic emissions). Operational alert criteria (functional relationship between the landslide(s) and monitored parameters) seem to be essentially data-driven, exploiting expert-guided thresholds in most cases.

Table 1 Number of operational lo-LEWS according to alert parameters [P, precipitation related; H, internal-to-slope hydrological variables; D, soil deformation] and criteria [heuristic, correlation, probabilistic] classification from Pecoraro et al. (2019)

Being process parameters directly associated with the landslide occurrence, those related to soil deformation are largely the most used. However, they enable a limited lead time (interval between warning issues and forecasted landslide events) and require a ‘perfect’ identification of the potentially unstable areas to be monitored. On the other hand, adopting alert parameters related to the precipitation features or soil wetting state allows longer lead times, which are crucial for proper community and assets protection and a coarser spatial coverage of the investigated area.

Under such premises, most territorial LEWSs (Guzzetti et al. 2020; Piciullo et al. 2018) assume as alert parameters observed and forecasted rainfall parameters such as intensity, duration, cumulations over different periods or a combination thereof while alert criteria are represented by thresholds discriminating conditions potentially leading to landslide occurrence.

Such thresholds may be quantified by exploiting physically based procedures (e.g., Napolitano et al. 2016) or data-driven approaches. These last may consist of visual fits (Posner and Georgakakos 2015), statistical (Althuwaynee et al. 2015), machine learning-driven (Pota et al. 2022) and probabilistic (Piciullo et al. 2017).

Physically based procedures usually rely on numerical approaches to derive threshold values through what-if scenarios or back-analysis of previous events. They can be time-consuming and require adequate knowledge of geomorphological parameters informing the landslide model (Calvello, 2017). On the other hand, data-driven approaches are less time and resource-demanding. However, proper calibration and validation of the approaches need long time series; furthermore, variations in slope features (e.g. geometry or land use) entail re-assessing the alert criteria.

Although the inaccuracy of all the abovementioned approaches may be due to improper identification of rainfall features and quality of rainfall datasets available for calibration and validation (e.g. Gariano et al. 2020), recently, some authors (e.g., Bogaard and Greco 2018) have argued that informing territorial LEWS with internal-to-slope hydrological variables (e.g., soil suction, soil water content) could improve prediction accuracy. LEWS information may be updated by:

• observations collected in situ (Comegna et al. 2016; Mirus et al. 2018; Wicki et al. 2020); they require significant economic and time efforts for data retrieving

• observations through remote sensing techniques (Thomas et al. 2019; Zhuo et al. 2019); they are usually freely available but characterised by coarse spatial resolution and long lag times between two observations

• theoretical computations based on climate reanalysis (Reder and Rianna 2021; Palazzolo et al. 2023), physically based distributed models (Segoni et al. 2018; Zhao et al. 2019; Abraham et al. 2020); they usually exploit authoritative and freely available physically based approaches of opportunity (e.g. made available for other and general purposes), but they can be based on soil parameterisations and spatial/temporal resolutions rarely matching the requirements of LEWS

• slope or basin-specific hydrological models (Pagano et al. 2010; Reder et al. 2018); time and resource requirements could prevent their adoption as operational landslide models in LEWS

For all mentioned cases, thresholds can be established according to physically based, data-driven or hybrid approaches (Reder et al. 2018).

An interesting virtual laboratory for LEWSs development and testing has been gradually established in the Campania Region (Italy), where, over centuries, rainfall-induced landslides have repeatedly affected widespread silty volcanic slopes (Picarelli et al. 2020a, b). Specifically, several phenomena occurred at the turn of the new millennium (e.g. San Pantaleone in 1997, Budetta 2002; Sarno in 1998, Del Prete et al. 1998; Cervinara in 1999, Fiorillo et al. 2001; Nocera Inferiore in 2005, Revellino et al. 2013; Ischia in 2006, Iovino and Perriello Zampelli 2007) significantly stimulated scientific investigations based on field monitoring, laboratory tests, experiments by physical models and numerical interpretations. According to such studies, several LEWS predictive tools have been developed based on rainfall interpretation. For example, some authors focused on features of triggering rainfall events, setting empirical intensity-duration thresholds (e.g., Guadagno 1991; Calcaterra et al. 2000). In the same vein, the approach adopted by the official Regional Alert System since 2005 and applied to the whole region defined alarm thresholds set on rainfalls cumulated over the last 24, 48 and 72 h, based on return periods of two (alert), five (pre-alarm) or ten (alarm) years. Even more recently, Piciullo et al. (2017) used an inventory of 96 shallow landslides and data from 58 rain gauges to set rainfall threshold equations identifying a potential trigger event by adopting a frequentist approach. However, such approaches disregard the significant role antecedent rainfall plays on soil conditions.

After the San Pantaleone landslide (1997), Fenelli (1998 not published) proposed an alert system for the Napoli-Salerno highway by interpreting a fifty-year dataset of rainfalls and landslides. Thresholds set on rainfalls cumulated over the last 2 months and 24-h drive alarm. The system relies on two indicators, i.e., rainfalls cumulated during the previous 2 months and 24 h, to be overtaken for issuing the alarm. Similarly, De Vita and Piscopo (2002) defined thresholds based on daily and cumulated rainfalls over the previous 59 days. The effects of the seasonality were taken implicitly into account by Napolitano et al. (2016). The authors suggested two different thresholds for the cold-wet and hot-dry seasons. Fusco et al. (2019) derived triggering thresholds through a physically based hydrologic model assuming an infinite slope stability approach, which also implicitly accounted for the effects induced by antecedent rainfalls.

The present work illustrates a slope hydrological model assuming, as the internal-to-slope variable, the suction levels predicted considering the effects of antecedent rainfalls and evaporation and requiring hourly rainfall datasets at least for calibration and validation. This procedure allows us to develop a threshold based on merely rainfall parameters, adequately cleaned up from runoff and evaporation according to physically based criteria. Its accuracy is assessed against two historical events that affected the area of interest and several rainfall histories that did not induce main landslide events. The paper also combines this model with empirical processing of rainfall-landslide or -no-landslide links observed over past decades, referring to poor time-resolute rainfall datasets. A deepening of the rationale of paper objectives and followed procedures is provided in the next paragraphs.

Rationale

The role of antecedent and triggering rainfalls was experimentally investigated by Rianna et al. (2014), interpreting monitoring results from a physical model of a lysimeter. Such a physical model consists of a tank filled with silty volcanic soils, standing on a capillary barrier (a geotextile, simulating a typical site condition of fractured bedrock) and exposed to weather elements over the years. The monitoring activity indicated that, at the beginning of the hydrological year, rains with different intensities infiltrated the cover almost entirely due to the topsoil's high potential infiltration. This last then reduced, making wetting-effective only in rainfall events characterised by long duration rather than significant intensity. Antecedent rainfalls progressively reduce over months suction levels in the cover, leading to a linear suction distribution, with zero at the base and positive values above. If an event decreases suction below the threshold of this linear trend, the capillary barrier breaks, generating drainage. A potentially triggering rainfall identifies with an event able to reduce suction up to approach a zero value throughout the cover. As the cover experiences drainage and low potential infiltration, hourly wetting increases are quantitatively limited but long-lasting. The persistency plays a significant role because interruptions longer than 1–2 h would determine a drainage-driven return to the initial linear distribution of suction induced by antecedent rainfalls. The same work also highlighted the effects of evaporation, which effectively limits the wetting increases of antecedent rainfalls.

By interpreting 10 years of rainfall records, including those triggering the Nocera Inferiore 2005 landslide in the Lattari Mountains of the Sorrento Peninsula (Campania Region, Southern Italy), Pagano et al. (2010) compared the performances of empirically against physically based approaches. The former was a two-threshold-based approach built on rainfall cumulated over 2 months and 24 h. The latter was a Richards equations-based approach, turning the rainfall history into suction evolution at different depths of a continuum domain. Adopting the landslide time as the reference time for calibrating the alarm thresholds of both approaches, the physically based one performed much better in terms of number of false alarm. The work also indicated that rainfalls of the last 4 months affected the cover hydrological state at the landslide time, with infiltrating rainfall amounts computing around fifty per cent of total rainfall. At the triggering time, suction was estimated to vanish throughout the cover under a persistent event lasting 14 h following an exceptionally wet antecedent period.

Coppola et al. (2020) postulated that suction nullification throughout the cover leads to landslide susceptibility as it entails the failure of the different root strength mechanisms at any depth, simultaneously making the cover susceptible to liquefaction, according to Olivares and Picarelli (2003). Reder et al. (2018) re-analysed the Nocera Inferiore 2005 case study by coupled thermal–hydraulic analyses accounting for evaporation. Comparison with results obtained running the Richards equation highlighted how predicting suction by simulating outgoing fluxes (evaporation) other than incoming ones (infiltrating rainfalls) nullified the number of virtual false alarms.

Aim

This work capitalises on data and knowledge from previous scientific investigations to develop a physically based LEWS procedure. This procedure attempts to cope with the following questions:

• Can it be possible to benefit from the coupled thermal–hydraulic approach accuracy without resorting to its online numerical solutions?

• Because the accuracy of the coupled thermal–hydraulic approach derives essentially from cleaning up rainfalls from evaporation and runoff, can this model be used to compute apart these amounts to quantify effective rainfalls more straightforwardly linked to slope safety conditions?

• Being slope safety conditions strictly related to slope suction levels, can effective rainfalls be linked to suction through simplified relationships?

This LEWS predictive tool assumes cumulations over the previous 4 months (C_4m) and the last (persistent) event (C_PLE) as variables linked to suction. Parametric studies running the coupled hydro-thermal model provide simplified criteria to clean up rainfall parameters (C_4m and C_PLE) from runoff and evaporation. The transformed variables C′_4 m- C′_PLE are linked to suction levels to build iso-suction lines by processing simplified rainfall scenarios with the coupled hydro-thermal model.

Rainfall records with reasonable time resolution, hourly at least, have to supply the procedure to discriminate between infiltrating and runoff amounts, individuate rainfall interruptions, and, in case of landslides, establish the time of the phenomenon. The empirical procedure synergically used with the physically based one interprets total amounts C_4m-C_PLE and bases on rainfall-landslide links recorded over the last decades, being rainfalls often provided with poor time resolution. The two procedures adjust typical variable conditions of rainfall data that may potentially inform LEWS for calibration in the Campania Region. A large amount of older rainfall data at the daily resolution, collected for several decades since the Fifties up to around 2000, are provided by rain gauges not always close to the susceptible slopes; in these circumstances, poor representative of slope rainfalls. These highly sized datasets can be exploited for empirical interpretations but are supposed to not result in accurate calibrations of forecasting tools. After 2000, automatic rain gauges enabled the collection of hourly (or finer) data while, at the same time, the spatial resolution of the observation network significantly increased. Such data allows better detection of rainfall duration and discrimination of runoff from infiltrating rainfall, making these data suitable for quantifying boundary conditions in a physical interpretation of the rainfall history. However, they are poorly sized and not suited for empirical interpretation unless merged with less recent daily data.

This work presents a procedure to estimate slope safety conditions trying to address and overcome the limitations identified in already proposed approaches for the Lattari Mountains area; in the specific:

- the procedure is developed to work on a hybrid scale between that represented by the single slope scale (lo-LEWS) and those at the regional scale (territorial LEWS); it can represent an added value for the investigated test cases where several slopes (§2) are characterised by very similar geomorphological features (exposure, bedrock, stratigraphy).

• the approach exploits the coupled thermo-hydraulic approach, calibrated and validated for the back-analysis of previous landslide events in the area (Reder et al. 2018), to retrieve simple relationships suitable for the time and resource constraints of LEWS by recurring scenario analysis

• the alert criteria adopt a hybrid approach according to which physically based approaches are used to derive iso-suction locus of points, but those acting as the threshold are identified by the back-analysis of previous events

• the same theoretical approach is used for defining alert parameters over the previous 4 months (C_4m) and the last (persistent) event (C_PLE) but also for determining the locus of points assumed as the threshold, overcoming the physical inconsistencies affecting different literature approaches.

Then, the ambition is that of using both an empirically and a physically based approach. Older data calibrate the former; recent data supply the latter to estimate suction from rainfall points. In this way, the LEWS prediction benefits from the accuracy provided by a physically based interpretation of reliable data without neglecting the information content of less recent data. The general attempt is to overcome the traditional dichotomy between empirically and physically based approaches, considering them synergic rather than an alternative.

The reference test-site and case histories

The geological and geomorphological context

The reference test site is the Lattari Mountains of the Sorrento Peninsula (Campania Region, Southern Italy), in the area where the Penisula connects to the mainland (Fig. 1). Over millenniums, materials erupted by the volcanic area of Campi Flegrei, 50 km away from the site, and, more recently, by the Vesuvius volcano, 20 km away, mantled repeatedly the calcareous slopes. These materials form cover thicknesses that tend to decrease with increasing slope gradients, ranging from a few meters on gentle slopes to a few decimetres on steep slopes. Soil grain size distribution is primarily silty sand and sandy silt, with the finest component being non-plastic. Pumices layers are sometimes interbedded between silty layers tending to thin, up to disappear, when the slope gradient increases. The porosity of silty covers is often very high, exceeding 70%, with peaks up to 80%. The hydraulic conductivity at saturation exceeds 10^–7 m/s, often resulting in around 10^–6 m/s. It implies a draining behaviour upon saturation, which prevents the formation of aquifers. On slopes tilted between 30 and 45 degrees, rainfalls triggered several landslides over centuries. The kinematics of these phenomena are typical of debris avalanches, according to classification by Hungr et al. (2014). They are essentially due to liquefaction phenomena (Olivares and Picarelli 2001), sometimes followed by fluidification.

Fig. 1

Location of landslides occurred throughout Lattari Mountains (labels link the position to the rainfall points of Fig. 2, Fig. 3, Fig. 4) and of weather stations (S) providing rainfall data for this work (S2 falls outside the map, about 3.5 km east of S1 and S3)

Insights from past events

Figure 1 indicates the location of the main landslides censed since the Fifties. The most impacting phenomena are N. 17 (2 January 1971), N. 3 (10 January 1997) and N. 10 (4 March 2005). In the following, a brief description of the geomorphological context and features of triggering rainfalls is reported. Four rain gauges provided rainfalls in the area (Fig. 1). The first three operated in the Nocera Inferiore district: Station-1, placed in the Nocera Inferiore town 3 km away from susceptible slopes, provided daily records over 1950–1996; the second one, very close to the first one, giving hourly records since 1997 up to now; Station-3, placed at the slope toes, provided hourly records over 1998–2010. The fourth one (Station-4 in Gragnano town) provided daily rainfall from 1950 to 1971.

The 1971 phenomenon occurred on the Monte Pendolo slopes near Gragnano town (Fig. 2). Unfortunately, the unavailability of hourly rainfalls at Station-4 does not permit drawing up details on the evolution of the last event. Thus it is only identified with the cumulative value over the last three days, C_3days, resulting in 190 mm. The cumulative value over the last 4 months resulted in C_4m = 780 mm. The landslide dynamic entailed a downward propagation that progressively increased the wideness of the landslide area, drawing the observed triangular shape. The slope attains 37° at the triggering apical zone, and the cover thickness slightly exceeds two meters. The landslide mobilised around 25,000 m³ of material along a slope development of about 300 m. The soil mass impacted a hotel and some houses causing six victims.

Fig. 2

Monte Pendolo landslides occurred in 1971 (phenomenon N. 17 in Fig. 1)

In 1997, a landslide occurred along the northern slope of St. Pantaleone hill, a small limestone hill rising 200 m in closeness to Nocera Inferiore town (Fig. 3). The landslide was triggered at the hilltop, where the slope gradient is about 37°, and the silty volcanic cover is 2 m thick. Rainfalls recorded at Station-2 indicated C_3days = 163 mm and C_4m = 1008 mm. The hourly resolution set at that time allowed detecting that the last event occurred persistently over 39 h, with rainfall intensity below 5 mm/h over 20 h and overtaking 10 mm/h over 3 h, with a peak of 17 mm/h. The road climbing on the involved slope was affected by significant deformation 9 h before the landslide occurrence. The plan view of the landslide area also indicates a tendency to enlarge downward, according to two rectangularly shaped zones of different wideness. The landslide mobilised around 4500 m³ of material along a slope development of about 200 m. The soil mass invaded the highway Napoli-Salerno (Fig. 3), causing one victim.

Fig. 3

San Pantaleone landslides occurred in 1960 (N. 1 in Fig. 1), 1972 (N. 2 in Fig. 1) and 1997 (N. 3 in Fig. 1) (The picture refers to event N. 3)

On 4 March 2005 at 4:00 p.m., rainfalls (recorded at Station-3 with hourly resolution) triggered a main phenomenon along the Monte Albino hill (Fig. 4) under C_3days = 149 mm and C_4m = 1038 mm. Rainfall intensity resulted between 10 mm/h and 16 mm/h over 6 h and below 10 mm/h over 10 h. C_3days approached 200 mm at the end of the landslide day, generating other minor landslides all around the zone of the main phenomenon.

Fig. 4

Nocera Inferiore landslide occurred in 2005 (N. 10 in Fig. 1); a landslide area; b picture of the landslide taken few days after the event; c detachment front of the landslide at the apical zone 3 months after the event; d grain size distribution of the involved material

The landslide area (Fig. 4a, 4b) was 24,600 m², and the mobilised soil mass was 33,000 m³. The propagation dynamic was similar to that introduced for the first case study (N. 17). At the triggering apical zone, the slope gradient (39°) significantly exceeds the average slope value (36°) (de Riso et al. 2007). In addition, the bedrock was 2 m deep, while, in the landslide area, the depths ranged between 1 and 2 m. The cover arises as homogeneous at the triggering zone (Fig. 4c) with the grain size distribution of silty sand (Fig. 4d).

Empirical link between rainfall point and slope safety conditions from available datasets

For the three cases mentioned above, as in most cases of Fig. 1, the cover thicknesses at the triggering zone were about two meters. Smaller thicknesses should be even more susceptible (Coppola et al. 2020) and likely to have been involved in numerous non-fatal (and, for this reason, hardly inventoried) events of reduced magnitude.

Stratified deposits originating from several past eruptions typically characterise susceptible slopes. They arise relatively uniformly from a hydraulic and mechanical point of view (Forte et al. 2019). The geological inhomogeneity turns into a simplified geotechnical homogeneity. For all cases, landslides triggered under the impact of a prolonged wet period synergic with a significant last event. The cases documented by hourly records indicate how the triggering event occurred persistently and was not always intense.

Figure 5 plots rainfall points recorded at the four mentioned stations since the Fifties, including only those arising with C_3days exceeding 85 mm. The fill large-sized points indicate landslide (their location, displayed in Fig. 1, includes the three phenomena already presented), while empty small-sized points indicate no-landslide. Due to the daily resolution of most of the recordings, it is impossible to quantify the persistent rainfall amounts C_PLE, so C_PLE is again identified with roughly C_3days as done in the three case studies. Bounded by threshold lines at C_4m* and C_3days*, it is possible to identify a safe region where rainfall points indicate a ‘no-landslide zone’ and an unsafe region where points associated with landslides get mixed with several no-landslide points. If a LEWS adopted this criterion as a predictive tool, major accuracy would meet in assessing stable conditions if the point fell in the safe region. In contrast, scarce accuracy would pursue in issuing an alarm if the point fell in the unsafe region. In this latter case, the prediction would be excessively conservative, as many no-landslide points would correspond to false alarms. Therefore, this predictive tool may be preliminary used to establish if a landslide danger exists (point falling in the unsafe region) based on the broad experience provided by the long-term dataset. However, a more accurate prediction is needed to establish if a point located in the unsafe region is associated with a wetting state predisposing to a landslide.

Fig. 5

Representation of significant rainfall histories in terms of antecedent rainfalls C_4m (4-month cumulation) against last events C_PLE (3-day cumulation). (Symbols relate to rain gauge locations; greater symbols relate to landslide events). A list of main rainfall events is in the Appendix 1

A physically based approach to assess the slope-wetting state

General description of the procedure

When weather observations at proper time resolution are available, the slope wetting state may be assessed following a physically based procedure, turning real-time weather evolution into internal-to-the-domain hydrological variables.

Nevertheless, according to an ‘off-line’ procedure (Greco and Pagano 2017), records, rather than being processed in real-time as boundary conditions, are matched with a pre-analysed weather scenario to borrow the suction it predicts. This is because real-time processing records could result in time-expensive computations unsuitable for LEWS requirements of prediction promptness. The equivalence between a recorded weather history and an analysed scenario has been established in terms of transformed rainfall point (C′_4 m, C′_PLE), expressing the algebraic sum of infiltrating rainfall and actual evaporation. According to such simplified procedures, rainfall records, cleaned up of evaporation (§4.2) and runoff (§4.3), have been turned into cumulative values (C′_4 m, C′_PLE) that identify the pre-analysed weather scenario yielding suction levels.

Suction levels have been identified with the computed mid-depth suction, s_h/2 (= − u_w-h/2). The choice of s_h/2 as the proxy variable synthesising the domain hydrological state results from shapes that suction profiles assume upon persistent wetting, starting from initial conditions of hydraulic equilibrium. Profiles approach a zero-distribution according to bell-shaped distribution (Coppola et al. 2020), in which suction nullifies for the latest around the domain middle part. Mid-depth suction is then used as the proxy for the entire suction distribution.

The final predictive tool consists of an abacus displaying iso-suction lines of s_h/2 (ISLs) in the C′_4 m- C′_PLE plane, computed theoretically by analysing a large number of weather scenarios (see §4.4). During a rainfall event, the (C′_4 m, C′_PLE) evolution extracted from records crosses sequential ISLs according to a progressive drop in domain suction levels. Alarm conditions are associated with approaching the point to a well-defined iso-suction line assumed as the ‘alarm line’ (AL). Criteria for selecting and validating AL are provided in §4.5.

The computational schemes adopted to achieve the different mentioned steps are reported in Fig. 6. According to Pagano et al. (2010) and Pirone et al. (2015), a one-dimensional flow problem has been assumed through a domain two meters thick. In addition, a seepage surface has been applied at the domain bottom (Reder et al. 2017), simulating the presence of fractured bedrock.

Fig. 6

1D domains analysed to determine a the K_c coefficient, b the potential infiltration, c suction levels

The hydraulic and thermal properties of the models have been calibrated, assuming the parameters of the material involved in the Nocera Inferiore landslide (Reder et al. 2018, Table 2).

Table 2 Hydraulic and thermal soil parameters (from Reder et al. 2018)

Elaboration of a simplified procedure to estimate evaporation

Evaporation over the antecedent period is calculated according to the approaches suggested by the FAO guidelines (Allen et al. 1998) suitably adapted to the case of pyroclastic soils. As the topsoil water content decreases, the actual evaporation, E, is estimated by reducing the so-called reference evaporation, Er.

E_r represents the maximum outgoing flux driven by meteorological forcing for a reference vegetation cover, assuming the topsoil wetting state can fully satisfy the atmosphere evaporative demand. FAO guidelines compute it as a function of atmospheric variables (net radiation, air temperature, air relative humidity) (Rianna et al. 2018). Actual evaporation is derived by exploiting as:

$$E={K}_{c} {K}_{r} {E}_{r}$$

(4)

where K_c and K_r account for the cover type and the actual topsoil wetting state, respectively.

Over a reference time span Δt_a equal to 1 month, K_r is assumed primary depending on:

• the monthly mean evaporative demand K_cE_r

• the mean time interval, Δt_r, between wet days, d_j:

Δt_r=Δt_a/ sum(d_j)

• the mean intensity, i_r, of rainfall events:

i_r=sum(i_j ×d_j)/sum(d_j)where K_c increases under low Δt_r, higher i_r and lower E_r.

To be considered a wet day, the cumulative precipitation over the day has to exceed the mean daily potential evaporation K_cE_r.

Predictive plots yielding K_r as a function of Δt_r, i_r and E_r have been arranged according to the representation proposed by Allen et al. (1998) to quantify evaporation. These variables have been parametrised in a set of numerical analyses carried out by integrating the Wilson equations (Wilson 1990; Wilson et al. 1994) over the domain of Fig. 6a. These equations mathematically describe the process of liquid and vapour water flow through an unsaturated medium under non-isothermal conditions (see Appendix 2 for details). Thermal and geotechnical parameters have been calibrated experimentally according to values presented previously (see Table 1). The analyses assumed three main rainfall intensities i_r at the top for the event, which are 10, 30 and 60 mm/day. Such a choice aims to refine the plots reported in FAO 56 guidelines, where only 10 and 40 mm/event are considered. For different E_r and Δt_r values, the outgoing fluxes E have been computed by the analyses over Δt_a and K_r (ratio E/ K_cE_r) has been quantified (Fig. 7). Operationally, the specific plot is selected based on i_r. It then yields K_r depending on Δt_r for the selected iso-K_cE_r curve. Finally, E_r is estimated by exploiting the relation suggested in FAO guidelines.

Fig. 7

Numerically computed Kr coefficient for silty volcanic soils against mean time interval between rainfall events and different potential evaporation intensities (values in-between 1 and 9 mm/day) for three daily rainfall intensities: a i_r = 10 mm/day; b i_r = 30 mm/day; c i_r = 60 mm/day

If the estimated i_r is in-between two available plots, i_r1 and i_r2 (i_r1 < i_r2), K_r may be quantified interpolating linearly by the corresponding K_r1 and K_r2:

$${K}_{r}=\frac{{\mathrm{i}}_{r}-{\mathrm{i}}_{r1}}{{\mathrm{i}}_{r2}-{\mathrm{i}}_{r1}}\left({K}_{r2}-{K}_{r1}\right)+{K}_{r1}$$

(5)

Plots of Fig. 7 allow for estimating the magnitude order of K_r expected during the different seasons. For example, the summer season is associated with long dry-hot periods, implying large K_cE_r and long Δt_r. All these conditions contribute to scarce water availability for evaporation, resulting in K_r lowering at around 0.3–0.4, that is, actual fluxes occurring to 30–40% of the potential one. On the other hand, the winter season is associated with long cold-wet periods, implying low K_cE_r and short Δt_r. As a result, the significant availability of water for evaporation rises K_c at around 0.8–1, with actual outgoing fluxes occurring near or at the potential levels.

Elaboration of a simplified procedure to estimate rainfall infiltration

The first step towards an estimation of rainfall infiltration is to redefine the criterium to quantify C_PLE based on knowledge of the hourly rainfall history. According to Rianna et al. (2014) and the technology usually adopted for rain gauges, the minimum hourly rainfall value is assumed to ensure an event’s continuity is 0.2 mm. This parameter is about an order of magnitude lower than the saturated permeability. As a general rule, it should be verified when calibrating and validating the approach. For Station-3, Fig. 8 shows how rainfall points are modified, replacing the C_PLE = C_3days criterion. Furthermore, the same records are modified (Fig. 9) by ruling out evaporation losses from C_4m by discriminating according to the different periods of the year of interest for slope stability in the Campania Region.

Fig. 8

Representation of significant rainfall histories in terms of antecedent rainfalls C_4m (4-month cumulation) against last events, C_PLE, represented as both 3-day cumulation and cumulation over time without interruptions

Fig. 9

Rainfall points against rainfall points net of evaporation; the points are clustered according to different periods of the year

The leftward shifting of rainfall points strictly correlates to the season (see Fig. 9). It is significant for the Autumn points (Sep-Oct and Nov-Dec), as the 4 months are typically featured by large cumulative atmospheric evaporative demand. However, in the remaining months, the magnitude of reduction varies and is driven by interannual variability.

The potential hourly infiltration rate, h*, depends on the surface hydraulic gradients, increasing with the topsoil suction levels, and on topsoil hydraulic conductivity, conversely decreasing with topsoil suction levels. The primary factor affecting h* is hence topsoil suction levels. As these levels fluctuate over time, h* also does. The simplified procedure considers only two constant h* values, respectively h*_PLE, associated with the last event from its beginning, t*, to its end, t_e, and h*_4 m, associated with the antecedent rainfalls, from their beginning, t₀, to t*.

In quantifying h*_PLE, the wet soil state over t*– t_e is assumed to be landslide-prone, hence approaching saturation. Accordingly, the entering gradient is supposed to match the gravity one, with h*_PLE thus identifiable with the saturated hydraulic conductivity (1 × 10^–6 m/s = 3.6 mm/h).

In quantifying h*_4 m, h* has been averaged over the time t₀– t*. The value h₀* at t₀ has been computed by analysing the scheme of Fig. 6b by Richards equation (Richards 1931) using SEEP/W code (Geo-Slope 2008). A value of 40 kPa has been assumed as the initial condition throughout the domain, according to the monitoring data provided by various authors (e.g., Rianna et al. 2014; Reder et al. 2018) that found this value representative of mean conditions acting at the end of the first wet period of the hydrological year (conventionally starting in September). Zero pore water pressure has been applied at the top domain over 1 h, computing the infiltrating hourly amount h₀*. Averaging between h₀* and h*_PLE, h*_4 m resulted equal to 5.8 mm/h.

For rainfall points provided by Station-3, Fig. 10 compares rainfall points C_4m-C_PLE with points of infiltrating rainfalls. As cumulation reductions affect C_4m and C_PLE, the points shift leftward and downward. It may be observed that, on average, around 50% of C_PLE infiltrates the domain, as ascertainable by comparing the values over the y-axis.

Fig. 10

Rainfall points against rainfall points net of runoff amount

Drawing iso-suction lines

To set up the predictive tool, the general assumption is that the rainfall point (C′_4 m, C′_PLE) of a given weather history relates to mid-depth suction (s_h/2(t_e) =  − u_w-h/2(t_e)). This relationship has been built numerically by solving the Richards equation for the scheme of Fig. 6c, subject to several weather scenarios. Several ‘wetting rules’ or water flux relationships (WFRs) are adopted for C′_PLE and C′_4 m-PLE, being the latter defined as C′_4 m net of C′_PLE (C′_4 m-PLE = C′_4 m – C′_PLE).

A constant-continuous time function, WFR₁ is applied over the time span t_e–t* (Fig. 11). As the applied infiltrating rate has been set invariably at the domain potential infiltration (h*_PLE), the duration of the event, t_e–t*, varies with C′_PLE according to the linear relationship t_e–t* = C′_PLE / h*_PLE. Three different linear functions, WFR_2i, have been adopted to spread C′_4 m-PLE over the antecedent time span t*– t₀ (Fig. 11). These functions envisage that the same C′_4 m-PLE could result in different wetting evolutions if (1) regularly spread over time, according to a constant relationship (WFR₂₁); (2) acting more severely farther to C′_PLE, according to a linear decreasing relationship (WFR₂₂); and (3) acting more severely closer to C′_PLE, according to a linearly increasing relationship (WFR₂₃).

Fig. 11

Water flux relationships, WFR, applied at the top of the domain (Fig. 6c) simulating antecedent rainfalls and last event

For each couple (C′_4 m, C′_PLE), WFR_2i has been varied according to these three relationships, and the maximum computed u_w-h/2(t_e) has been assumed as the final effect of (C′_4 m, C′_PLE):

$$\begin{array}{ccc}{U}_{\mathrm{w}-\mathrm{h}/2}({t}_{\mathrm{e}})\hspace{0.17em}=\hspace{0.17em}\mathrm{max}[{g}_{1}({{C}^{\mathrm{^{\prime}}}}_{\mathrm{PLE}}& {{C}^{\mathrm{^{\prime}}}}_{4\mathrm{ m}}({\mathrm{WFR}}_{21}))& \begin{array}{ccc}{g}_{2}({{C}^{\mathrm{^{\prime}}}}_{\mathrm{PLE}}& {{C}^{\mathrm{^{\prime}}}}_{4\mathrm{ m}}({\mathrm{WFR}}_{22}))& \begin{array}{cc}{g}_{3}({{C}^{\mathrm{^{\prime}}}}_{\mathrm{PLE}}& {{C}^{\mathrm{^{\prime}}}}_{4\mathrm{ m}}({\mathrm{WFR}}_{23}))]\end{array}\end{array}\end{array}$$

For the analysis, an initial constant suction throughout the domain of 40 kPa has been assumed, consistently with the procedure previously described to compute rainfall infiltration.

Figure 12 plots u_w-h/2(t_e) against C′_4 m for three different C′_PLE. Each curve shows that u_w-h/2(t_e) first significantly increases with C′_4 m and then tends asymptotically to a constant value, identifying the maximum effect induced by the antecedent weather period. Worth noting is that for low C′_4 m values, u_w-h/2(t_e) decreases with increasing C′_PLE, while the reverse occurs for high C′_4 m (above 450 mm). In the former case, increasing C′_PLE at the same small amount C′_4 m implies a significant percentage decrease in C′_4 m-PLE and a delay of the wetting process, with less effective imbibition at mid and high depths.

Fig. 12

Numerically computed u_w(h/2-t_e) against C′_4 m for different C′_PLE

The same numerical results are rearranged by plotting the iso-u_w(h/2)(t_e) lines in the plane C′_4 m- C′_PLE (Fig. 13). Again, lines arise downward concavely and provide each point of the plane with the physical meaning of the suction level. This way, a rainfall history, post-processed as a sequence of (C′_4 m, C′_PLE) points, displays a series of suction levels.

Fig. 13

Iso-u_w(h/2-te) curves numerically computed in the C′_4 m – C′_PLE plane

Represented lines at 0, − 1, − 2 and − 3 kPa are linked to possible landslide-prone states for slope degrees approaching the soil friction angle (Pagano et al. 2010; Balzano et al. 2019), according to the case studies presented. The line at 0 kPa lower bounds the region where pore water pressure is invariably null under infiltrating intensity matching the draining one.

Regardless of the criteria followed in selecting the alarm iso-suction line (they will be discussed in the next paragraph), it could be stated that for C′_4 m lower than a threshold C′_4 m*, hardly a C′_PLE could induce a landslide-prone state. Conversely, for C′_4 m exceeding a second threshold, C′_4 m**, maximum effects caused by antecedent weather conditions are attained. Consequently, the triggering event is minimised and no longer affected by further C′_4 m increments. In short, while in the former case, neither C′_PLE should be feared to trigger a landslide, in the latter, a not exceptional C′_PLE is expected to induce a landslide-prone state.

If C′_4 m is in-between C′_4 m* and C′_4 m**, the C′_PLE threshold strongly relates to the antecedent weather conditions.

Selection and validation of the alarm iso-suction line (AL)

The LEWS prediction requires turning real-time recorded rainfalls in rainfall point (C′_4 m – C′_PLE) and evaluating the point distance from an hypothetical Alarm iso-suction line (AL).

To minimise missing and false alarms, AL may be selected as the line closest (upper-bounding) to the safety region displayed by the cloud of all recorded no-landslide points. The cloud of points provided by Station-3 results satisfactorily bounded by the line at 0 kPa (Fig. 14), hence identified as trial AL. It is worth noting that if runoff and evaporation were neglected, the cloud of rainfall points would not cluster in a clear pattern. As a result, developing an alarm criterion would become challenging (several points would fall in the zero suction zone). The seasonal dynamics introduced by cleaning up rainfall points of rainfall infiltration and evaporation permit a clear identification of safe and unsafe zones. In particular, all autumnal rainfall points are almost invariably far from the zero suction zone, flattened in the safe area at zero C′_4 m due to evaporation fully tackling rainfall infiltration. Three rainfall points falling in winter approach the zero suction zone, while only one point enters the zero suction zone.

Fig. 14

Rainfall points, rainfall points net of runoff amount and rainfall points net of runoff and evaporation amount, compared with AC curve (iso-u_w(h/2-t_e) curves)

AL may be validated considering rainfall points associated with weather histories resulting in landslides. Fig. 15 shows the evolution of hourly rainfall points of events triggering the Nocera Inferiore 2005 (Fig. 15a) and S. Pantaleone 1997 (Fig. 15b) landslides, respectively. Rainfall points are expressed as total rainfalls (C_4m, C_PLE), rainfalls net of runoff and rainfalls net of runoff and evaporation (C′_4 m, C′_PLE).

Fig. 15

Iso-u_w(h/2-te) curves in the plane C′_4 m – LE′ and processed rainfall points of a Nocera Inferiore 2005 landslide and b San Pantaleone 1997 landslide

The Nocera Inferiore 2005 (Fig. 15a) triggering event started on 4 March 2005 at 0.00. The first plotted point refers to 3.00 a.m.. After that, the C′_4 m, C′_PLE points moved progressively towards AL, approaching it just before the triggering time, which occurred at 4.00 p.m.. Therefore, this case study fully validates the selected AL.

The S. Pantaleone 1997 (Fig. 15b) triggering event started on 9 January 1997 at 6.00 a.m. The first point represented refers to 8 a.m. of 9 January. The C′_4 m, C′_PLE points moved progressively towards AL, overtaking it at 8.00 a.m. of the following day, 10 January, when damages to the road climbing the S. Pantaleone hill were notified. The point continued to travel within the zero-suction zone until the main instability was observed at 8.00 p.m., with C′_4 m, C′_PLE quite away from AL. This case study validates the selected AL against the time of the first visible effects induced by the triggering event.

Comparison with operational and state-of-the-art predictive models proposed for pyroclastic covers of Campania Region

Three consolidated approaches are selected among those listed in §1 to compare their predictions with that provided by the approach hereabove presented:

• the operational alert regional system relying on rainfall thresholds associated with a fixed return period (here set at a 10-year return period)

• the rainfall threshold proposed by De Vita and Piscopo (2002), processing daily data and accounting for antecedent precipitations cumulated over 59 days;

• the probabilistic intensity-duration threshold suggested by Piciullo et al. (2018) with the exceedance probability level here fixed at 50%.

Precipitations provided by Station-3 over the 10 years, including the Nocera Inferiore landslide in 2005, are used to evaluate the predictive performances of approaches, assessed in their ability to yield threshold exceedance (alarm) at the landslide time and below-threshold predictions (no false alarms) over the other time. Results are plotted in Fig. 16 for the selected literature approaches and in the previously illustrated Fig. 14 for the proposed approach.

Fig. 16

Application of reference EWS thresholds for the area of study: a operational alert regional system; b De Vita and Piscopo (2002) threshold; c Piciullo et al. (2018) threshold. These thresholds have been used considering rainfall histories in Fig. 14 and precipitation data from Station-3 to detect the event of March, 4 (2005). Points are classified according to specific months of the year

The first and third approaches (Fig. 16a and Fig. 16c) show that the threshold was exceeded only at the time of the Nocera Inferiore landslide (4 March 2005). The predictions yielded by the second approach (Fig. 16b) are instead of threshold exceedance at the landslide time, and at several different times, probably due to the coarse resolution of precipitation data it processes (daily scale). The first and third approaches seemingly show slightly better performances than that here-above proposed, which is affected by a false alarm (Fig. 14). Nevertheless, it is important to point out how they indicate several close-to-landslide conditions for rainfall events of September–October (Fig. 16a and Fig. 16c, green symbols), falling in a period that, because following the dry season, results in the past unaffected at all by landslides. Predictions are so close to threshold lines to require, in any case, the activation of alert procedures preceding alarm if the models were implemented in a LEWS. This erroneous assessment stems from neglecting antecedent rainfall and evaporation phenomena that, in these months, lead to very dry conditions in the subsoil, making even significant rainfall events inconsequential. Instead, this feature is well considered by the proposed approach that delivers rainfall points for the same events (Fig. 14, green symbols) significantly far from the threshold line. Furthermore, the two literature methods require the adoption of sophisticated (hardly explainable to the practitioners/civil protection operators) approaches to retrieve thresholds. At the same time, the developed procedure tries exploiting, at any stage, simple physically based but robust laws to derive the resulting rainfall histories.

Conclusions

The proposed framework allows the evaluation of the occurrence of weather-induced landslides using an offline strategy. According to it, monitored weather histories permit depicting several scenarios by exploiting two proxies: the soil water budget in the antecedent period and the cumulative infiltrated rainfall in correspondence with the potential triggering event.

The alternative of performing an online solution of the Wilson model, updating time-by-time boundary conditions consistently with the recorded weather variables, is not manageable under the time constraints of a LEWS.

Previous investigations have already used the two proxies to account for antecedent rainfall and triggering event. The work rearranges these ingredients as a common information basis for a twofold predictive tool referring sequentially to an empirically and a physically based approach. Specifically, a traditional ‘weather-based’ approach, empirical and disregarding the specific slope features, is synergised with a ‘slope-tailored’ approach, physically based and accounting for geometry and hydraulic properties of the slope.

Furthermore, in the physically based approach, slope safety conditions are evaluated through an empirical approach. It is due to several constraints: first, there is an unbeatable difficulty in characterising soil strength under the presence of root systems and liquefiable conditions; second, the different simplifications, affecting the procedure, lead the predicted suction to represent a virtual variable working as a comparative indicator of different weather histories rather than an accurate estimation of the actual suction distribution. As such, the suction threshold does not affect the assessment of the slope safety factor. The implicit assumption is that possible inaccuracy in predicting suction may also influence threshold estimation to a similar extent without affecting the accuracy of the LEWS prediction.

Overtaking the AL and entering the prone-to-landslide region does not necessarily imply landslide occurrence. Indeed, the investigated landslides take place when (1) the cover attains a wetting state predisposing to landslide propagation and (2) local factors sparking the sliding mehanism occur (Pagano et al. 2010). Local factors determining the landslide to happen in one place rather than another (Greco and Pagano 2017) may consist of internal-to-slope singularities (e.g., weaker layers, inhomogeneities in hydraulic properties, fractures, concentrated erosion and other peculiarities concentrating fluxes in the triggering zone). They may also include external factors, such as rock elements impacting the slope, blasts, animal passage, etc. In this work, the LEWS attempts to predict the predisposing wetting state, as it may somehow be rationalised. Instead, it avoids guessing what local factor could trigger the landslide once the predisposing state has been attained. This is because local factors are challenging to be detected and strongly influenced by aleatory. This explains why, for instance, it should be not surprising that for the San Pantaleone case rainfall point enters well inside the prone-to-landslide region ahead of landslide initiation or that several neighbouring slopes, although similar for stratigraphy and geometrical features, remain stable under the impact of the same weather history. The perspective to alarm, as soon as a predisposing state to landslide is recognised, is very similar to that successfully pursued to predict the occurrence of snow avalanches (in this case, the proxy for the predisposing state is the snow temperature).

The procedure targets a well-defined slope typology for soil properties and geometry. Reference has then been made implicitly to a slope scale. As silty pyroclastic covers mantling slopes around a volcanic area usually feature limited variability for hydrological and strength properties, extending the procedure to a regional scale should not require significant effort. Furthermore, geometrical features in terms of thicknesses, stratigraphy and slope degrees reduce the need to analyse a few representative schemes of the type presented in this work to shift from slope to regional scale.

Being the procedure based on simplified assumptions, each one is susceptible to criticisms indicating possible sources of inaccuracy. For most steps, the accuracy may be improved at the expense of a higher degree of complexity and heavier development efforts. This increasing complexity should be well worth the effort. For instance, a needed improvement could be further improving time resolution by attaining sub-hourly values. It could prevent great mistakes in computing rainfall infiltration: intense-short lasting events of few minutes, for most flowing as runoff, may lead to overestimating infiltration amounts if coarse time resolution data are available.

Appendices

Appendix 1

Table 3

Table 3 Date, antecedent rainfall over 4 months (C_4m), rainfall over 3 days (C_3days), reference rain gauge, and if a landslide event occurred or not (from Fig. 5)

Appendix 2

The formulation adopted for the water balance equation is:

$$\frac{1}{{\gamma }_{w}}\frac{\partial \left({u}_{a}-{u}_{w}\right)}{\partial t}=\frac{1}{{\gamma }_{w}{m}_{2}^{w}}\left[\frac{\partial }{\partial z}\left({k}_{w}+\frac{{k}_{w}}{{\gamma }_{w}}\frac{\partial \left({u}_{a}-{u}_{w}\right)}{\partial z}\right)+\left(\frac{{P}_{a}+{u}_{v}}{{P}_{a}{\rho }_{w}}\right)\frac{\partial }{\partial z}\left({D}_{v}\frac{\partial {u}_{v}}{\partial z}\right)\right]$$

(A1)

where ${u}_{w}$ = liquid pore water pressure [kPa]; ${u}_{a}$ = pore air pressure [kPa]; ${u}_{v}$ = partial pressure of pore water vapour [kPa]; ${m}_{2}^{w}$ = the gradient of the water retention curve [1/kPa]; ${P}_{a}$ = total atmospheric pressure [kPa]; ${k}_{w}$ = the hydraulic conductivity function [m/s], ${D}_{v}$ = the soil vapour diffusivity function [(m kg)/(kN s)]; ${\rho }_{w}$ = liquid water density; ${\gamma }_{w}$ = liquid water specific weight [kN/m³]; $t$ = the transpiration flux [m/s].

The formulation adopted for the heat balance equation is:

$${C}_{h}\frac{\partial {T}_{s}}{\partial t}=\frac{\partial }{\partial z}\left(\lambda \frac{\partial {T}_{s}}{\partial z}\right)-{L}_{v}\left(\frac{{P}_{a}+{u}_{v}}{{P}_{a}}\right)\frac{\partial }{\partial z}\left({D}_{v}\frac{\partial {u}_{v}}{\partial z}\right)$$

(A2)

where ${T}_{s}$ = soil temperature [°C]; ${C}_{h}$ = the volumetric specific heat function [J/(m³°C)]; $\lambda$ = the thermal conductivity function [W/(m°C)]; ${L}_{v}$ = latent heat of water vaporisation [J/kg].

The formulation adopted for the thermodynamic equilibrium, which forces the energy balance between vapour and liquid water linking partial vapour pressure, matric suction and temperature, is:

$${u}_{v}={u}_{v0} exp\left(\frac{\left({u}_{a}-{u}_{w}\right){M}_{w}g}{R{T}_{s}{\gamma }_{w}}\right)$$

(A3)

where ${u}_{v0}$ = saturated partial pressure of pore vapour [kPa]; ${M}_{w}$ = water molecular weight [= 0.018 kg/mol]; $R$ = the ideal gas constant [= 8.314 J/(mol °C)]; $g$ = gravitational acceleration [= 9.81 m/s²].