Brief communication “A prototype forecasting chain for rainfall induced shallow landslides”

Although shallow landslides are a very widespread phenomenon, large area (e.g. thousands of square kilometres) early warning systems are commonly based on statistical rainfall thresholds, while physically based models are more commonly applied to smaller areas. This work provides a contribution towards the filling of this gap: a forecasting chain is designed assembling a numerical weather prediction model, a statistical rainfall downscaling tool and a geotechnical model for the distributed calculation of the factor of safety on a pixel-by-pixel basis. The forecasting chain can be used to forecast the triggering of shallow landslides with a 48 h lead time and was tested on a 3200 km 2 wide area.


Introduction
Shallow landslides are a recurrent and dangerous geohazard responsible for casualties and economical damages.To forecast shallow landslide occurrences two different approaches are usually used: empirical approaches based on rainfall thresholds and physically based deterministic models.While the first is widely used for operational applications at the regional scale (Aleotti, 2004;Cannon et al., 2011;Martelloni et al., 2012;Rosi et al., 2012;Lagomarsino et al., 2013), the latter are more commonly applied at the slope or catchment scale (Dietrich and Montgomery, 1998;Baum et al., 2002;Segoni et al., 2009;Apip, 2010).The main constraints thwarting the effective application of such models over larger areas are twofold: (i) poor knowledge of the spatial organisation of the values of the input parameters (Segoni et al., 2012); (ii) need for a relevant computational time which is not compatible with real-time applications (Rossi et al., 2013, this special issue).
In this work, we present a methodology that overcomes these issues and we simulate the prototype of a fully physically based forecasting chain for the timely prediction of the occurrence of rainfall induced shallow landslides over large areas (thousands of km 2 ).

Material and methods
The presented forecasting chain is composed by a series of tools in which the first feeds the subsequent one and so on (Fig. 1a).A numerical weather prediction model provides two kinds of rainfall forecasts, with two different spatial resolutions and two different lead times; a statistical tool downscales these rainfall maps to get a finer spatial resolution compatible with the analysis scale of the subsequent slope stability simulator; the hydrological module of the slope stability simulator uses these rainfalls to model the increase in the water pore pressure; and the geotechnical module of the same slope stability simulator perform a pixel by pixel analysis to identify when and where the increase in water pore pressure is expected to trigger shallow landslides.

Numerical weather prediction model: COSMO LM
The weather prediction model adopted in the simulation chain is COSMO LM (Doms and Schattler, 1998).It is a limited-area non-hydrostatic (i.e.no approximation is used for vertical momentum equation) computational model for weather prediction developed within the framework of the European Consortium for Small-Scale Modelling.This model can be used to solve the physical phenomena in the meso-β-scale (20 to 100 km) and meso-γ -scale (2 to 20 km).Limited area models provide atmospheric forecasts only on a specific region of earth.Therefore, initial conditions and transient boundary conditions are obtained by combining data from global models and observations.In our simulations we used short range weather forecasts made for a time period up to 48 h.Due to the chaotic nature of the atmosphere, in this time range the forecasts are generally more accurate than those obtained with longer time ranges.
In our forecasting chain, the basic horizontal resolution adopted is 7 km.This configuration allows an accurate numerical prediction of near-surface weather conditions, focusing on clouds, fog, frontal precipitation, and orographically and thermally forced local wind systems.Moreover, this is the official resolution currently used in Italy by the National Weather Service (Aeronautica Militare) and by the Functional Centers of the Italian Civil Protection.
However, when the orography of the analyzed area is particularly complex, a 7 km resolution could be considered a too rough approximation.Therefore, a higher spatial resolution (2.8 km) modelling is implemented in the forecasting chain, nested in the previous one.Through this configuration, it is expected to obtain a better simulation of severe weather events triggered by deep moist convection, such as supercell thunderstorms, intense mesoscale convective complexes, prefrontal squall-line storms and heavy snowfall from wintertime mesocyclones.The main drawback of the 2.8 resolution is the computational burden (in terms of time and resources); therefore, in the development of an early warning system, the use of this configuration should be carefully evaluated.

Statistical downscaling tool
Even the 2.8 km resolution can sometimes be inappropriate to accurately model the complex spatial patterns of precipitation, which are markedly influenced by the morphology and the geographic location of the area.Moreover, slope stability analyses are performed at very finer spatial resolutions (usually from 5 to 20 m, 10 m in this work), bringing the necessity to feed stability models with rainfall maps at comparable resolutions.
Currently, no physical model is able to provide reliable estimations at this resolution; therefore, the simulation chain makes use of statistical downscaling techniques, following a mixed approach that combines radial basis function (RBF) and multivariate regression method (MRI) techniques.The first is a deterministic interpolation creating surfaces from measured points based on the degree of smoothing with distance; the latter performs a disaggregation of precipitation on the basis of topographical variables (Baillargeon, 1989;Vigier, 1981;Hay et al., 1998).
As a first step, at a given time, the precipitation value R predicted by the numerical weather prediction model (NWP) in a grid point is disaggregated according to the equation: where X 1 , X 2 ...X n represent a set of topographic variables (to which the spatial variation of the precipitation is assumed to be' related); a 0 through a n are the calculated regression coefficients and residual represents the rainfall rate that does not depend on the chosen variables.Based on Antofie (2009), we used a set of topographic variables constituted by slope, aspect and elevation.Regression coefficients, assumed as invariants over the domain, are retrieved through least square method.Further details on the approach are reported in Burrough and Mc Donnell (1998).
As second step, the residual values in each point of the coarse grid are interpolated by using the RBF statistical method proposed by Ninyerola et al. (2000) and Agnew and Palutikof (2000).The rainfall values are estimated at each grid point (i,j) of the higher resolution domain (required by the slope stability model) through the equation: R i,j = a0 +a 1 Elevation i,j + a 2 Aspect i,j +a 3 Slope i,j + residual (RBF) (2) In this way, part of precipitation is physically interpolated by using topographic variables (MRI) and part by a statistical approach (RBF).

Landslides prediction model (HIRESSS)
The landslide prediction model used in the forecasting chain is the HIgh REsolution Slope Stability Simulator (HIRESSS) (Rossi et al., 2013).
The HIRESSS code is a physically based distributed slope stability simulator for analyzing shallow landslide triggering in real time, on large areas, using parallel computational techniques.The physical model proposed is composed of two parts: hydrological and geotechnical.The hydrological one receives the rainfall data from the downscaled COSMO-LM model as dynamical input and computes the pressure head as perturbation to the geotechnical stability model, which provides results in factor of safety (FS) terms.
The hydrological model is based on an analytical solution of an approximated form of Richards equation under the wet condition hypothesis and it is introduced as a modelled form of hydraulic diffusivity to improve the hydrological response.The geotechnical stability model is based on an infinite slope model and it takes into account the increase in strength and cohesion due to matric suction in unsaturated soils, where the pressure head is negative.The soil mass variation on partially saturated soil caused by water infiltration is also modelled.
HIRESSS computes the factor of safety at each selected time step (and not only at the end of the rainfall event) and at different depths within the soil layer.In addition to rainfall, the model input data are constituted by slope gradient, geotechnical and hydrological parameters and soil thickness (Rossi et al., 2013).Furthermore, in order to manage the problems related to the uncertainties in the main hydrological and mechanical parameters, a Monte Carlo simulation has been implemented.HIRESSS can theoretically work at any spatial resolution; in this work we adopted a 10 m resolution because it corresponds to the grid size of the DEM of the area and because it represents a fair compromise between spatial accuracy and computational resources needed.
HIRESSS code was developed to run over multiprocessor systems and was tested for performances with an increasing number of processing units to design an optimal cost/benefit approach covering the entire prediction chain, from rainfall data acquisition to the factor of safety computation.Further details on HIRESS can be found in Rossi et al. (2013; in the same special issue).

Test case description
The test case selected is located in Tuscany (Italy), it includes part of the northern Apennines and has an extent of 3103 km 2 (Figs. 2 and 3a).The study area is mainly mountainous and shows two different geological settings.On the west sector carbonaceous rocks prevail and the slope gradients are frequently higher than 60 • , the east sector is characterised by flysch formation rock-type (Macigno) and the slope gradients vary from 0 • in the alluvial plains to 55 • .In the midand upper sections of the valleys, where most landslides usu-ally occur, the stratigraphy consists of a 1.5 to 5 m thick layer of colluvial soil overlying the bedrock (Casagli et al., 2006).
The forecasting chain was tested by simulating the rainfall event of 30 October 2003 -1 November 2003, during which an intense rainstorm affected part of the study area triggered 50 shallow landslides (Fig. 3).
Since it is well established that soil thickness spatial organisation has a marked influence on slope instability (Segoni et al., 2012), the GIST model (Catani et al., 2010) was applied to the study area to get a spatially distributed soil thickness map.Geotechnical and hydraulic parameters were differentiated on a lithological basis: 40 survey points (Fig. 3a) were analysed with a series of in situ and laboratory tests, including grain size analysis, measurement of Atterberg limits, phase relationship analysis, constant head permeameter test using an Amoozemeter (Amoozegar, 1989), borehole shear test (Lutenegger and Halberg, 1981;Tofani et al., 2006) and matric suction measurements with a tensiometer.

Results
The output of the forecasting chain is a series of instability maps (one for each time step) where each pixel is associated with a probability of instability (FS < 1).In this simulation, we obtained 58-hourly time steps, covering a timeframe from 30 October 2003 00:00 to 1 November 2003 10:00 UTC.
Figure 3b-e show the temporal evolution of the factor of safety in the study area.For a better visualisation of the results, pixel values were aggregated for each 1st order watershed, according to the methodology proposed by Rossi et al. (2013); moreover, the single time steps are aggregated in time frames with 12 and 3 h widths.
The various time frames allow for following the forecasted evolution of the instability.The simulation results are shown in Fig. 3.During 30 October and part of the subsequent day (Fig. 3b), a generalised situation of stability can be observed: this is in accordance with the experimental data since no landslides were triggered during this period.Official reports state that all landslides were triggered in a time window ranging from the evening of 31 October to the early hours of 1st November.The simulation matches quite well with the ground truth observed during this period: instability conditions are approached within the first hours of the 31 October, they consistently and sharply increase at the end of the day (Fig. 3c) and reach the peak during the first hours of 1 November (Fig. 3d); afterwards, during the morning, the simulated stability conditions progressively tend to normalise (Fig. 3e).The night between the 31 October and 1 November is therefore correctly identified as the most critical period; moreover, the simulation indicates the western part of the region as the most probable for landslide triggering, followed by the central northern area.The spatial distribution of landslides is partially in accordance with this outcome.The eastern part of the test area and the spot in the northern-central sector appear affected by an overestimation of the slope failure probability if it is compared with the number and the position of landslides contained in the official reports.A possible partial explanation is a lack of reporting because these areas are scarcely populated mountainous regions, while most of reported landslides involved infrastructure or water streams.
At this stage of the research, a quantitative validation of the results cannot be performed, because it would require the definition of calibration procedures to translate the probabilistic outputs into warning levels.However, the results of the simulation demonstrated the technical feasibility of the implementation of the forecasting chain into an early warning system.

Implementation of the simulation chain within early-warning systems
The organisation of the warning system is unavoidably influenced by the computational time of the single tools.The basic rainfall forecasts provided by the Cosmo-LM model are two: a rough forecast (7 km spatial resolution), with 48 h lead time, and a detailed forecast (2.8 km spatial resolution) with 24 h lead time.The different time range depends on the additional computational time required to obtain the weather forecasts with a finer spatial resolution.The downscaling algorithm needed a few minutes to accomplish its task, while HIRESSS returned a series of 58 factor of safety maps (58hourly time steps to simulate the whole period) in 32 min, with a spatial resolution of 10 m and for a total of over 50 million pixels.The computational time of the whole forecasting chain is compatible with an operative use for civil protection purposes.
Figure 1b shows the organisation of a prototype early warning system where the computations of the forecasting tools are organised based on the runtime needed to obtain the output to be used by the subsequent tools of the chain.A first level of warning is available with a 48 h lead time directly using the rough (7 km) precipitation forecasts.A second level of warning can be obtained with a 24 h lead time using the detailed (2.8 km) precipitation maps.
We set up a prototype forecasting chain for the prediction of rainfall induced shallow landslides over a large areas (e.g.thousands of km 2 ).The chain is composed by a physically based weather prediction model, a statistical tool for the downscaling of the rainfall maps and an advanced slope stability simulator composed by a hydrological and a geotechnical tool.The simulations performed in this work demonstrated the feasibility of the chain and allowed a preliminary estimation of the degree of reliability of the forecasts.
Nevertheless, before integrating the forecasting chain into a fully operative warning system, further improvements are needed to overcome the main weak points identified in this case of study: -Simulation of other well-documented case studies to calibrate a procedure for the interpretation of the probabilistic results (e.g.definition of probability values corresponding to alert thresholds).
-Deeper comprehension of the spatial organisation of the input parameters (e.g.rainfall, soil thickness and geotechnical parameters) to further improve the quality of the results.

Fig. 1 .
Fig. 1.Flow chart of the simulation chain (a) and time windows of the prototype warning system (b).

Fig. 3 .
Fig. 3. Study area (a) and simulations of slope instability forecasts according to different time frames (b to f).