Temporal evolution of landslide hazard for a road infrastructure in the Municipality of Nocera Inferiore, Italy, under the effect of climate change

In recent years, landslide events have extensively affected pyroclastic covers of the Campania Region in southern 10 Italy, causing victims and conspicuous economic damages. Due to the high criticality of the area, a proper assessment of future variations in landslide occurrences and related risk is crucial for policy-makers, administrators and infrastructure stakeholders. This paper addresses work performed within the FP7 INTACT project, having the goal to provide a risk framework for critical infrastructure while accounting for climate change. The study is a part of the testing and application of the framework in the Campania region, assessing the temporal variation in landslide hazard specifically for a section of the Autostrada A3 "Salerno15 Napoli" motorway, which runs across the toe of the Monte Albino relief in the Municipality of Nocera Inferiore. In the study, hazard is defined as the yearly probability of a spatial location within a study area to be affected by landslide runout given the occurrence of rainfall-related triggering conditions. Hence, hazard depends both on the likelihood of rainfall-induced landslide triggering within the study area and the likelihood that the specific location will be affected following landslide runout. Landslide triggering probability is calculated through the application of Bayesian theory and relying on local historical rainfall 20 data. Temporal variations in triggering probability due to climate change are estimated from present-day to the year 2100 through the characterization of rainfall patterns and related uncertainties using the EURO-CORDEX Ensemble. Reach probability, defining the probability that a given spatial location is affected by debris flows, is calculated spatially through numerical simulation of landslide runout. The temporal evolution of hazard is investigated specifically in the proximity of the motorway, as to provide a quantitative support for landslide risk analysis. 25

As stressed in these works, the issue of climate change issue and its impacts can be considered "a fantastic example of 'very deep' uncertainty". Nevertheless, given the extent of potential impacts on communities (Paris Agreement 2015) including their economic dimension (Stern 2006;Nordhaus 2007;Chancel & Piketty 2015), considerable efforts have been spent in recent 5 years devoted to assessing the variations in frequency and magnitude of weather-induced hazards related to climate changes (Seneviratne et al. 2012). A variety of strategies have been devised and implemented with the aim of detecting the main sources of uncertainty and their extent (Wilby & Dessai 2010;Cooke 2014;Koutsoyiannis & Montanari 2012;Beven 2015).
Among weather-induced hazards, investigations on future trends in the occurrence and consequences of landslides and on the uncertainties in their estimation have received relatively limited interest (Gariano & Guzzetti 2016;Beven et al. 2015). Possible 10 concurrent causes include the mismatch between the usual scale of analysis for landslide case studies and the current horizontal resolutions of climate projections, as well as the extraordinarily relevant role of case-specific geomorphological features, which hinder the generalization of findings to other contexts.

Previous studies of pyroclastic landslides in Campania
Despite the above limitations and, indeed, in the attempt to address them, several recent studies have focused on future 15 variations in the occurrence of landslides affecting pyroclastic covers mantling the carbonate bedrocks in the Campania Region in Southern Italy. These studies considered different test cases; namely: Cervinara (Damiano & Mercogliano 2013;, Nocera Inferiore (Reder et al. 2016;Rianna et al. 2017aRianna et al. , 2017b and Ravello (Ciervo et al. 2016). Several aspects differentiate the case studies and, consequently, the investigations performed in them. For example, depth, stratigraphy and grain size of pyroclastic covers are fundamentally regulated by slope, distance to volcanic centres  Vesuvio), as well as wind direction and magnitude during the eruptions; such differences induce variations in rainfall patterns recognized as effective for slope failure (e.g. intensity, length of antecedent precipitation time window). For these reasons, while daily weather forcing data have been found to result in better assessments for the Cervinara and Nocera Inferiore test cases, sub-daily data have been found to improve the quality of assessments for the Ravello test case. Consequently, daily observations modified according to projected anomalies (Damiano & Mercogliano 2013) or daily data provided by climate 25 simulations subjected to statistical bias correction are used in the former cases, while a stochastic approach is adopted with bias-corrected data to provide assessments at hourly scale for the latter. Moreover, in some studies (Reder et al. 2016;Ciervo et al. 2016;Rianna et al. 2017aRianna et al. , 2017b, slope stability conditions are assessed through expeditious statistical approaches referring to rainfall thresholds, while physically based approaches are preferred in other cases. Finally, climate projections at 8km in the optimized configuration over Italy (Bucchignani et al. 2015) andthe Zollo et al. (2014) configuration of 30 COSMO_CLM model (the highest resolution currently available for Italy up to 2100) are used as inputs in the aforementioned case studies, while climate projections from the Euro-CORDEX multimodel ensemble (Giorgi et al., 2016) are adopted in Rianna et al. (2017b). Nat. Hazards Earth Syst. Sci. Discuss., https://doi.org/10.5194/nhess-2017-416 Manuscript under review for journal Nat. Hazards Earth Syst. Sci. Discussion started: 5 December 2017 c Author(s) 2017. CC BY 4.0 License.

Object of the study
The present study focuses again on the Nocera Inferiore site, and also makes use, as will be discussed, of rainfall data from the sites of Gragnano and Castellammare di Stabia. The geographical collocation of the three towns in Italy and in the Campania region is illustrated in Fig. 1.
The study presents significant elements of novelty. For instance, through a Bayesian approach, it characterizes precipitation 5 values cumulated on two time windows as proxies for the triggering of landslides affecting pyroclastic covers in the Monti Lattari mountain chain. The resulting quantitative model returns temporal variations in triggering probability, thus accounting for the effect of climate change on rainfall trends. Uncertainties in variations in rainfall patterns are taken into account recurring to the EURO-CORDEX ensemble. Projections provided by climate simulations are bias-adjusted, allowing the comparison with available physically-based rainfall thresholds while adding further assumptions and uncertainties in simulation chains. 10 Landslide runout is also investigated probabilistically through a frequentist estimate of reach probability performed in a GIS environment, thus allowing the seamless mapping of landslide hazard under current and future climate change scenarios.

Geographic and geomorphological description
Most of the territory of the Nocera Inferiore municipality belongs geomorphologically to the Sarno river valley. The most 15 urbanized area of the town is located at the toe of the northern slopes of the Mount Albino relief, pertaining to the Monti Lattari chain (Fig. 2, sector A); other more sparsely populated areas are located at the foot of the Torricchio hills (Fig. 2, sector B). These reliefs are constituted by carbonate rocks covered by air-fall pyroclastic deposits originated from volcanic eruptions (Somma-Vesuvio complex) during the last 10,000 years (Pagano et al. 2010). Such covers in loose pyroclastic soils have been historically affected by multiple types of flow-like rainfall-induced landslides; among the most relevant events: Gragnano, 20 1997;Sarno & Quindici, 1998;Nocera, 2005;Ischia, 2006 (moreover, see Table 2 for a complete list of events affecting the area investigated in the work during 1960-2015 time span), including: (a) hyper-concentrated flows, which are generally triggered by washing away and/or progressive erosive processes along rills and inter-rill areas; (b) channelized debris flows, which can be generated by slope failure in ZOB areas (Dietrich et al. 1986;Cascini et al. 2008); and (c) un-channelized debris flows, which are locally triggered on open-slopes areas propagating as debris avalanches. The latter type characterized the 25 most recent event which affected the city in March 2005, causing three fatalities and extensive damage to buildings and infrastructures (Pagano et al. 2010;Rianna et al. 2014). This study focuses specifically on a section of the Autostrada A3 "Salerno-Napoli" motorway, which runs across the toe of the Monte Albino relief as shown in Fig. 3. Nat. Hazards Earth Syst. Sci. Discuss., https://doi.org/10.5194/nhess-2017-416 Manuscript under review for journal Nat. Hazards Earth Syst. Sci. Discussion started: 5 December 2017 c Author(s) 2017. CC BY 4.0 License.

Method of analysis
The study is conducted by coupling mathematical software with GIS in order to obtain spatially referenced estimates and allow mapping of hazard. The study area was modelled into the GIS software through a digital terrain model (DTM) having a resolution of 15x15 m. The original resolution adopted in the Regione Campania ORCA project (2004) was 5x5 m. A variety of DTM resolutions were tested for the case study. The adopted resolution proved to be sufficient to adequately represent the 5 surface morphology and landslide runout as detailed in Section 6. Hazard is estimated quantitatively for each cell of the GISgenerated grid through the following model: in which is the probability of landslide triggering, and is the reach probability for the cell. It should be noted that the probability of triggering is related to rainfall parameters and, thus, is assumed to be uniform for the entire area, while reach probability depends on geomorphological factors, and is thus cell-specific and spatially variable within the area. 10 The study is conducted according to a modular approach initially involving the disjoint estimation of triggering probability (including its temporal variation) as described in Sec. 5, and reach probability, as detailed in Sec. 6. Subsequently, hazard is calculated using the model described above. The operational flowchart of the study is shown in Fig. 4. Source datasets are presented in Sec. 4.

Observed precipitation data
Observed datasets are used to identify time windows used as proxies for landslide triggering, to calibrate the Bayesian approach (see Sec. 5) and finally data from the Nocera Inferiore station are used to bias adjustment of climate projections (see Sec. 6).
Although the study is focuses on Nocera Inferiore landslide events, data from the neighbouring towns of Gragnano and Castellammare di Stabia are considered in order to increase the size of the event database, thus increasing the statistical 20 significance of the approach. At both sites, landslide events affecting pyroclastic covers were observed to be very similar to those of the Nocera Inferiore slopes (De Vita & Piscopo 2002) as described in Sec. 4.2.
The dataset related to daily precipitation spans across the time window January 01, 1960-December 31, 2015. Unfortunately, for the three towns, no weather stations were in operation throughout the entire period. Consequently, the dataset is reconstructed by merging data provided by different weather stations. Prior to 1999, the network of monitoring stations was 25 managed by Servizio Idrografico e Mareografico Nazionale (SIMN, Hydrographic and Tidal National Service) network at national level. In that period, the selected reference weather station is that located within the town and identified with the town's name as can be found in the SIMN yearbooks. Subsequently, the management was delegated to regional level, with the Regional Civil Protection managing the dataset for the Campania region. Since 1999, the reference weather stations are selected among those adopted for the towns in Regional Early Warning Systems against geological and hydrological hazards (Sistema 30 Nat. Hazards Earth Syst. Sci. Discuss., https://doi.org /10.5194/nhess-2017-416 Manuscript under review for journal Nat. Hazards Earth Syst. Sci. Discussion started: 5 December 2017 c Author(s) 2017. CC BY 4.0 License.
di Allertamento Regionale per il rischio idrogeologico e idraulico ai fini di protezione civile, 2005). Checks for the homogeneity of time series and for the unwarranted presence of breakpoints between the two periods were carried out for this study through the Pettitt (1979) and CUSUM (CUmulative SUM) (Smadi & Zghoul 2006) tests. Source weather stations, location, installation time and main (i.e., at least four months in a year) out-of-use periods are reported in Table 1.

Landslide inventory 5
The inventory of landslide events was compiled using three main references: Vallario (2000), De Vita & Piscopo (2002) and, for the more recent events, the "Event Reports" drafted by the Regional Civil Protection. It is worth recalling that only events affecting pyroclastic covers have been considered and included in the dataset. Sixteen events were observed in the period 1960-2015; these are reported in Table 2.

Climate projections 10
Climate projections for Nocera Inferiore were conducted as a preliminary step to the quantitative characterization of the temporal evolution of triggering probability as detailed in Sec. 5. The adopted simulation chain includes several elements. (compared to preindustrial era) of about 2.6, 4.5, 6 and 8.5 W/m 2 . Such scenarios force Global Climate Models (GCM). These are recognized to reliably represent the main features of the global atmospheric circulation, but fail to reproduce weather conditions at temporal and spatial scales of relevance for assessing impacts at regional/local scale. In order to bridge such gap, GCMs are usually downscaled through Regional Climate Models (RCMs). These are climate models nested on GCMs, from which they retrieve initial and boundary conditions, but which work at higher resolution (including a non-hydrostatic 20 formulation) on limited area. The dynamic downscaling from GCMs to RCMs allows a better representation of surface features (orography, land cover, etc.) and of associated atmospheric dynamics (e.g., convective processes). Nevertheless, persisting biases can hinder the quantitative assessment of local impacts.
In order to cope with such shortcomings, a number of strategies can be adopted. For instance, to characterize uncertainty associated to future projections, climate multi-models ensemble can be utilized where different combinations of GCM and 25 RCM run on fixed grid and domain. Furthermore, statistical approaches (e.g., Maraun 2013;Villani et al. 2015;Lafon et al. 2013) can be pursued to reduce biases assumed as systematic in simulations. More specifically, quantile mapping approaches have been applied with satisfactory results in recent years for different impact studies. In these applications, the correction is performed as to ensure that "a quantile of the present day simulated distribution is replaced by the same quantile of the presentday observed distribution" (Maraun 2013). However, limitations and assumptions associated to these approaches should be 30 clear to practitioners (Ehret 2012;Maraun & Widmann, 2015). In the present study, climate simulations included in EURO-CORDEX multi-model ensemble at 0.11 ' (approximately 12 km) are considered under the RCP4.5 and RCP8.5 scenarios as described in Table 3. Climate simulations are bias-adjusted through an empirical quantile mapping approach (Gudmundson et al. 2012) using data from Nocera Inferiore weather stations from the period 1981-2010.

Calculation of triggering probability 5
Landslide triggering probability was estimated quantitatively as a function of two cumulative rainfall thresholds; namely, the 1-day rainfall 01 and the 59-day rainfall 59 . Such cumulative parameters were calculated using a moving window procedure for and associated with each day from January 01, 1960 to December 31, 2015 from the observed precipitation data described in Sec. 4.1. The number of landslide events observed for each day at the Nocera Inferiore, Gragnano and Castellammare di Stabia as reported in the landslide inventory was associated with the rainfall data. The conditional probability of triggering of a landslide given the simultaneous occurrence of 01 ( ) and 59 ( ) is obtained from the Bayesian approach is estimated using a Bayesian approach as suggested by Berti et al. (2012). The procedure refers to 15 Bayes' theorem, formulated as follows: The above Bayesian probabilities can be computed in terms of relative frequencies as follows: in which   Fig. 7a shows the temporal variation of triggering 10 probability for scenario RCP4.5 for each of the 10 CORDEX models, along with the sample mean, minimum and maximum. Figure 7b shows the corresponding data for scenario RCP8.5. Figure 7c plots the outputs of both scenarios, while Fig. 7d plots the temporal variation of the sample coefficient of variation for both-scenarios. Such statistic is given by the ratio of the sample standard deviation and the sample mean, and describes quantitatively the scatter among CORDEX model outputs.
For the RCP4.5 scenario, considering the running 30-year averages, all available projections predict a moderate increase, with 15 few differences. In this regard, a higher spread among the models is recognizable at the middle of the XXI century. Such increased spread is mainly due to the outputs of two models constantly representing, respectively, the upper and bottom boundaries of the ensemble throughout the entire investigated period. Nevertheless, the increase in mean terms worth about Nat. Hazards Earth Syst. Sci. Discuss., https://doi.org/10.5194/nhess-2017-416 Manuscript under review for journal Nat. Hazards Earth Syst. Sci. Discussion started: 5 December 2017 c Author(s) 2017. CC BY 4.0 License.
35-40% at the end of the century. It is interesting to note how, in average terms, minor differences are recognizable between the two scenarios with the largest variations in 2041-2080 time span where the triggering probability is slightly higher under the more pessimistic scenario. Also in terms of coefficient of variations, in the first part of the century, slight differences are recognizable under the two scenarios; however, also in terms of concentrations of gases (green house of chemically active), higher differences arise in the second part of the century. 5 In this perspective, it is interesting to note how in the end of the century differences nullify probably due to joint effect of contrasting trends in proxies estimated, in average terms, under RCP8.5 (increase in values of β01 and reduction in β59). To this aim, it is worth to note that the effects of evapotranspiration, reducing the soil wetting and that could have a significant role because of increased warming, are neglected in the developed framework. On the other hand, in this case, the projections returned by a model included in the ensemble significantly differ from the other ones with PT values attaining about 3x10 -3 10 against a mean value of about 1.5x10 -3 . In this way, the coefficient of variation is significantly higher than this estimated under RCP4.5. However, several assumptions and constraints affect retrieved findings; among the other ones, land use/land cover features are assumed not significantly affecting proxy values during the calibration and future period; nevertheless, local conditions could substantially modify the landslide susceptibility (e.g. fires destroying vegetation) altering, at the same time, exposure (also assumed constant during the entire period). 15

Estimation of reach probability
Investigation of the spatial variability of landslide hazard entails the modelling of its downslope propagation (runout). Reach probability is the probability (from 0: certainty of no reach to 1: certainty of reach) of each point in the spatial domain being affected by the landslide during the runout process. Several morphological, empirical and physically-based approaches are available for quantitative runout analysis (Hürlimann et al. 2008); each of these may present advantages or weaknesses in 20 relation to site-and/or phenomenon-specific attributes, data availability and scale of the analysis. Consistently with the methods previously used to define triggering-rainfall scenarios, the approach used to define downslope runout scenarios is based on an algorithm involving stochastic modelling.

Reach probability calculation method
Landslide reach probability was computed spatially using Flow-R, a DTM-based distributed empirical model developed in the 25 Matlab ® environment (Horton et al. 2013). The flow-slide spreading is controlled by a flow direction algorithm that reproduces flow paths (Holmgren 1994) and by a persistence function to consider inertia and abruptness in change of the flow direction (Gamma 2000). The flow direction algorithm proposed by Holmgren (1994), in the setting used in this study (x=1, see Eq. (3) in Horton et al. 2013) is similar to the multiple D8 of Quinn et al. (1991Quinn et al. ( , 1995. The multiple D8 distributes the flow to all neighbouring downslope cells weighted according to slope. The algorithm tends to produce more realistic looking spatial 30 patterns than the simple D8 algorithm by avoiding concentration to distinct lines (Seibert & McGlynn 2007). The maximum Nat. Hazards Earth Syst. Sci. Discuss., https://doi.org/10.5194/nhess-2017-416 Manuscript under review for journal Nat. Hazards Earth Syst. Sci. One-run propagation simulation provides possible flow-paths generated from previously identified triggering/source areas. In this work, source areas were identified by means of the official geo-morphological map of the "Campania Centrale" River Basin Authority (PSAI 2015) and coincide with the union of the "zero order basin" (ZOB) and current "niche/failure" areas. 5 This hypothesis is in accordance with the requirement of consistency with accounts of historical events and with the aim to consider the most pessimistic possible triggering scenarios (i.e., those with maximum mass potential energy).
The reach probability for any given cell is calculated by the following equation: where and are the flow directions; is the probability value in the u-th direction; is the flow proportion according to the flow direction algorithm; is the flow proportion according to the persistence function; and 0 is the probability 10 determined in the previous cell along the generic computed path. The values are subsequently normalized. Runout routing is stopped when: (1) the angle of the line connecting the source area to the most distant point reached by the flow-slide along the generic computed path is smaller than a predefined angle of reach (Corominas 1996); and (2) the velocity exceeds a user-fixed maximum value, or is below the value corresponding to the maximum energy lost due to friction along the path. The values which do not fit the above-mentioned requirements are redistributed among the active cells to ensure conservation of the total 15 probability value.

Reach probability outputs
The propagation routine was applied to the DTM described in Sec. 3. An angle of reach of 4° was calibrated based on the geomorphological information (i.e., the extension of the slope fan deposition) and the official hazard maps of the Landslide Risk Management Plan of the River Basin Authority (PSAI, 2015). Consistently with the mean values reported by the scientific 20 literature (Faella & Nigro 2001;Revellino et al. 2004) for the same phenomena and in the same region, the maximum runout velocity was set at 10 m/s. Figure 8 illustrates the spatial distribution of reach probability at hillslope scale. Source areas are also indicated.
In this area, the highway runs mostly on a soil embankment. The road level is generally elevated with respect to the paths of the downslope flows. The propagation impacts the embankment and stops in front of -or laterally continues according to -the 25 topographic information and the model setting. Differently, in some points, the highway runs approximately at the same level of the fans, thereby allowing the propagating flow to invade the road. In both cases, damage or disruptions may be caused to the infrastructure. In order to overcome this distinction and to cover both scenarios, only flow propagation to the upstream boundary of the infrastructure are considered in the study. An illustrative example is shown in the magnified focus area in Fig.   9. Due to the reasons mentioned above, the road surface is only partially affected by the flow-slides. Confining the study just 30 Nat. Hazards Earth Syst. Sci. Discuss., https://doi.org/10.5194/nhess-2017-416 Manuscript under review for journal Nat. Hazards Earth Syst. Sci. Discussion started: 5 December 2017 c Author(s) 2017. CC BY 4.0 License. to a part of the infrastructure (e.g., from A to B), the runout values to be considered in the risk assessment should be taken along the section A-B (Fig. 10). Hazard can be calculated directly for a given year and RCP scenario by applying Eq. (1).
The results shown in Fig. 10 attest for the marked spatial variability of reach probability (and, therefore, of hazard) along the investigated section of the A3 motorway infrastructure.

Concluding remarks 5
This paper has illustrated a methodology for the quantitative estimation of rainfall-induced landslide hazard. An example application of the proposed method was conducted for a short section of a motorway. Despite the limited extension of the study area, the results displayed a marked temporal and spatial variability of hazard. The temporal variability of hazard is a consequence of climate change as parameterized through quantitative projections for concentration scenarios RCP4.5 and RCP8.5. Significant temporal variability was predicted for both concentration scenarios. The considerable spatial variability 10 resulting from the case study stems from the spatial variability of reach probability as modelled in the runout analysis.
The hazard outputs obtained by the method can be used directly in the quantitative estimation of landslide risk. The latter also requires the quantitative estimation of the vulnerability of human-valued assets (i.e., vehicles, persons, etc.) and the exposure (i.e., the number and/or degree of presence) of the assets themselves in the study area in a reference time period.
The quantitative estimates of hazard as obtained in this paper are pervaded by significant uncertainty. Among the main sources 15 of uncertainty are the climate change projections, the runout model and the Bayesian model developed to quantify triggering probability. These uncertainties are epistemic in nature, as they stem from the inherent difficulty in compiling climate change projections, the inevitable degree of approximation and imperfection in runout modelling capabilities, the limited rainfall and landslide occurrence data used to develop triggering probability curves. As such, increased modelling capability and improved databases could reduce the magnitude of uncertainty associated with hazard estimation. 20 Notwithstanding the above uncertainties, the quantitative estimation and assessment of the spatial and temporal variability of hazard provide an important decision support tool in the disaster risk management cycle; specifically, in the planning and prioritization of hazard mitigation and risk mitigation measures. The availability of quantitative methods allows a more rational decision-making process in which the costs and effectiveness of risk mitigation can be compared and assessed in terms of convenience. 25  1964,1965,1967,1981,1982 Tramonti (