ELSEVIER

Contents lists available at ScienceDirect

Engineering Geology

journal homepage: www.elsevier.com/locate/enggeo

A probabilistic approach for landslide hazard analysis

S. Lari * P. Frattini, G.B. Crosta

Dept. Earth and Environmental Sciences, Université degli Studi Milano Bicocca, Piazza della Scienza 4,20126 Milano, Italy

ARTICLE INFO ABSTRACT

We present a general framework for probabilistic landslide hazard analysis. With respect to other quantitative hazard assessment approaches, this probabilistic landslide hazard analysis has the advantage to provide hazard curves and maps, and to be applicable to all typologies of landslides, if necessary accounting for both their onset and transit probability.

The method quantifies, for a given slope location, the exceedance probability of being affected by a landslide with a specific local intensity within a reference time interval, i.e. the hazard curve, under the common assumption that landslides behave as a Poisson process. Hazard maps are calculated, reducing the hazard curve to single values by choosing a fixed probability of exceedance following standards or regulation requirements. The method is based on the assessment of a landslide onset frequency, a runout frequency for long-runout landslides, and the local definition of landslide intensity, which can be expressed through different parameters, according to landslide typology. For long runout landslides, the runout and spatially-varying intensity and uncertainty are considered.

Hazard curves and maps play a fundamental role in the design and dimensioning of mitigation structures, in land planning and in the definition of risk and hazard management policies. Starting from the general framework, we apply the methodology for rockfall hazard analysis, and we test it in an area affected by the Christchurch 2011 earthquake, New Zealand, which triggered a large number of rockfalls, killing five people.

© 2014 The Authors. Published by Elsevier B.V. This is an open access article under the CC BY-NC-SA license

(http://creativecommons.org/licenses/by-nc-sa/4.0/).

CrossMark

Article history: Accepted 22 July 2014 Available online 1 August 2014

Keywords:

Probabilistic landslide hazard analysis

Hazard maps

Hazard curves

Rockfalls

New Zealand

1. Introduction

Landslide hazard expresses the probability that a landslide with a certain intensity can occur in a certain location within a given period of time [ÍSSMGE Glossary of Risk Assessment Terms, http://140.112.12. 21/issmge/2004Glossary_Draft1.pdf]. This definition underlines that hazard is a function of intensity. This function is usually known as hazard curve, in the literature generally related to seismic risk (Frankel et al., 1996), windstorms and floods (Grünthal et al., 2006) and tsunamis (PTHA, González et al., 2009; Annaka et al., 2007; Liu et al., 2007; Geist and Parsons, 2006). While the concepts of intensity and magnitude are well defined for these threats, the terms are not always really and easily formalised for landslides. A formalization was proposed by Hungr (1997), and some clarifying advices have been expressed in some recommendations for landslide risk assessment (e.g. OFAT-OFFE-OFEP, 1997; AGS, 2007; Fell et al., 2008; Corominas et al., 2014). In most cases, however, intensity is used as a general term, which can include different concepts, such as size, volume, velocity, energy. Magnitude is frequently used to describe the size of a landslide in terms of volume (e.g. Hungr et al. 1999; Marchi and D'Agostino, 2004; Jakob and Friele, 2010; Santana et al, 2012) or area

* Corresponding author. E-mail address: serena.lari1@unimib.it (S. Lari).

(e.g.: Hovius et al, 1997; Stark and Hovius, 2001; Dussauge et al, 2003; Malamud et al., 2004; Guthrie and Evans, 2004).

A low consensus on the use of terms derives from the fact that landslides include different phenomena, which can be described by different parameters. Due to the objective difficulty to generate hazard curves for landslides, a reason why they are extremely rare in the landslide risk literature, the selection of an intensity parameter is still an important issue. Intensity should correspond to a measure of "damage potential". Hence, it should not express the size of a landslide, but its destructive power. On the other hand, the frequency of landslides is often related to their size and not necessarily to their destructive power, expressing the magnitude of the events, more than the intensity (e.g., Hungr et al, 1999; Dussauge et al, 2003).

The concepts of magnitude and intensity applied to landslides can be clarified referring, in analogy, to earthquake engineering. For earthquakes, the magnitude expresses the energy released by the single event, and can be considered a description of earthquake "size". Magnitude-Frequency relationships (also known as Gutenberg-Richter's law) are used to characterize the frequency of occurrence of earthquakes with different magnitude (Gutenberg and Richter, 1942). However, ground motion parameters or functions (e.g., peak ground acceleration, peak ground displacement, spectral acceleration), expressing the local intensity of the earthquake, are needed to assess damages by using fragility curves (i.e., the probability of exceeding a given

http://dx.doi.org/m1016/j.enggeo.2014.07.015

0013-7952/© 2014 The Authors. Published by Elsevier B.V. This is an open access article under the CC BY-NC-SA license (http://creativecommons.org/licenses/by-nc-sa/4.0/).

damage state as a function of a ground motion parameter) (Syner-G, 2011). To calculate the local ground motion, attenuation relationships as a function of distance from the earthquake epicenter and magnitude, are used, also accounting for uncertainty.

Techniques to derive hazard for each location along a slope can be different as a function of the typology of the landslide and the scale of the analysis. For local scale analysis of single landslides it is possible to simulate various scenarios considering different volumes (but also displacement or velocity, especially for already existing landslides with known volume) and associated probabilities (i.e., M-F relationships) through numerical models in order to determine the spatial distribution of intensity during landslide movement (Archetti and Lamberti, 2003; Friele et al, 2008). Hence, for each location along the slope it is possible to build the hazard curve by adopting the frequency values provided by M-F relationships and the intensities calculated by the model. This approach is also adopted for snow avalanches (Keylock et al, 1999; Keylock and Barbolini, 2001). In these methodologies, however, the uncertainties involved in modeling the landslide dynamics are not taken into account. If uncertainty is considered, the intensity at each location along the slope cannot be expressed as a single value for each magnitude scenario, but as a frequency distribution of values. To characterize this distribution, a simple statistic is normally used, such as the arithmetic average (Agliardi et al 2009), the maximum value (Gentile et al, 2008; Calvo and Savi, 2009), or a specific percentile (95th in Spadari et al., 2013; 90th in Lambert et al., 2012). In these cases, the hazard curves are obtained by associating these unique values of intensity to corresponding scenario frequencies derived from M-F relationships. However, these approaches introduce a strong assumption about the distribution of intensity, because the arithmetic mean is representative only for normally distributed intensity, the maximum value can reflect outliers of the distribution, and the percentiles may strongly overestimate the actual hazard.

The aim of the paper is to propose a probabilistic methodology for the assessment of hazard connected to all typologies of landslide, quantifying the probability of exceeding various intensities at a site (or a map of sites) given all possible events. In the second section we present the general framework for landslide probabilistic hazard analysis. In the third section, we discuss its applicability to all landslide typologies. In the fourth section, we decline the methodology to rockfall hazard analysis and we investigate the nature of the intensity distribution for rockfalls by means of parametric numerical modeling. In the fifth section, we present an application of the methodology to the area of Richmond Hill, Port Hills, Christchurch, New Zealand.

2. A general framework for landslide hazard analysis

The landslide hazard assessment methodology here proposed is conceptually derived from the numerical/analytical approach formalized by Cornell (1968) for probabilistic seismic hazard analysis (PSHA), which integrates over all earthquake scenarios, allowing to estimate the likelihood of exceeding selected ground motion parameters (generally peak ground acceleration, PGA) at a given site, within a reference time interval.

For each position along the slope, z, the probability of exceeding a certain value of landslide intensity, i, is

P(I>i) = f p(I)dI (1)

where p(I) is the probability density function of landslide intensity at the position z. This function reflects the stochastic nature of intensity, whose values can vary for each position along the slope, due to the uncertainty about the models used to simulate the intensity, and the temporal and spatial variability of the landslide behavior. The shape of the probability density function can be different (e.g., normal, log-

normal), based on both the nature of the physical processes and the types of uncertainty.

Multiplying the exceedance probability by the annual frequency of occurrence f, we obtain the annual rate at which i is exceeded, F(I > i) as:

F(I>i)= f ■ P(I>i) (2)

The annual frequency of occurrence,f can be calculated for landslides by direct or indirect approaches (Picarelli et al., 2005; Corominas and Moya, 2008). Direct approaches are based on the analysis of available historical data of past landslides, which can also be related to geology, geomorphology, and other factors (Moon et al., 1992; Cruden, 1997; Jaiswal et al, 2011; Geist et al, 2013). Indirect approaches derive the landslide frequency from triggering factors, such as rainfall intensity and duration (Sidle et al., 1985; Crozier, 1997; Dai and Lee, 2001; Schuster and Wieczorek, 2002; D'Odorico et al., 2005; Rosso et al., 2006; Salciarini et al., 2008; Frattini et al., 2009), or earthquake (Del Gaudio et al., 2003; Rathje and Saygili, 2008).

In case landslides scenarios with different magnitude potentially occur in a certain position along the slope, the total annual rate at which i is exceeded, Ftot (I > i), derives from the sum of all scenarios, s.

Ftot (I>i) = E fsPs(I>i) (3)

By assuming a homogeneous, stationary Poisson process for the occurrence of the events (Crovelli, 2000), the probability of exceeding each intensity i in the next Tyears from this annual rate,

Pp0iss, is:

Ppoiss(I>i, T) = 1 —e-F'°'T (4)

This represents the hazard curve for each position along the slope.

In order to represent hazard through a hazard map, it is necessary to reduce the hazard curve to a single value for each position. This is typically done by choosing the intensity value having a 10% (or 2%) chance of exceedance in 50 years (as done for earthquakes, Frankel etal, 1996). As a consequence, a map of these values for the region of interest can then be generated.

In the literature, the probability of landslides is frequently expressed in terms of annual frequency (e.g. Hungr et al., 1999) or return time, implicitly assuming a binomial occurrence probability model, for which, in fact, the exceedance probability equals the annual frequency. While this assumption holds in case of rare events (e.g. large rock avalanches), it can be violated for frequent events (e.g. rockfalls, debris flows, landslide reactivations) (Crovelli, 2000).

The use of a stationary Poisson process for the occurrence of the events implies the assumptions that the rate of occurrence of landslides is constant in time, and that the probability of more than one event in a small time interval is order of magnitudes lower than the probability of just one event (Straub and Schubert, 2008). The recurrence time of landslides deviates from that expected for a stationary Poisson process at short recurrence times, due to a temporal clustering effect, resulting from climatic conditions and/or seismic triggers (Geist and Parsons, 2008; Tatard et al, 2010; Witt et al, 2010). In this case, the occurrence rate may increase slightly immediately after an event, but is that of a stationary Poisson process on a larger timescale.

A time-varying occurrence rate however may be introduced, with an increase of the complexity of the analysis, as done for other hazards (Ogata, 1999; Parsons, 2008).

3. Application of the methodology to different typologies of landslide

According to landslide typology, the treatment of intensity can be different as a function of landslide dynamics. In particular, we can distinguish two main categories: slow moving landslides and fast long-runout landslides (rock avalanches, debris flows, rockfalls).

For slow-moving landslides, intensity is generally expressed in terms of displacement rate (Hungr, 1981; Mansour et al., 2011 Frattini et al, 2013), total displacement (Rathje and Saygili, 2008 Picarelli, 2011), or differential ground deformation (Bird et al, 2005 Negulescu and Foerster, 2010). Displacement rate, however, represents a velocity, and has to be referred to a specific time interval to express the effective intensity of a landslide. It could rather be useful in assessing the potential lifetime of facilities built on slowly moving sloping grounds, or in defining the time interval spanning consecutive maintenance operations (Picarelli, 2011). For each location, intensity can be described as a distribution of values (Eq. (1)) to account for the uncertainty related with the simulation of the behaviour of the landslide. In case of active landslides showing a continuous and slow downslope movement, in some cases due to periodic reactivations of existing instabilities (e.g., earth flows, earth slides, and DSGSDs), the intensity has necessarily to be referred to a defined time interval. This has to be consistent with the time T adopted for the Poisson probability (Eq. (4)). For hazard analysis of potential unknown landslides, multiple scenarios with different volumes and/or sliding surfaces should be integrated (Eq. (3)). For each scenario, a frequency of occurrence has to be calculated (Eq. (2)). On the other hand, for hazard analysis of single known landslides, the integration of different scenarios (Eq. (3)) is not required.

Fast long-runout landslides (e.g., debris flows, rockfalls, rock- and debris- avalanches) are characterized by onset and propagation processes, controlling the probability and the intensity of a landslide at a specific point in space (Pierson et al., 1990; Cancelli and Crosta, 1993; Jaboyedoff et al., 2005; Agliardi et al, 2009).

Intensity can be expressed with different parameters according to the landslide typology: kinetic energy for rockfalls (Corominas et al, 2005; Jaboyedoff et al, 2005; Agliardi et al, 2009); peak discharge (Jakob, 2005), velocity (Hungr, 1997; Bovolin and Taglialatela, 2002; Calvo and Savi, 2009) or depth (Shih et al, 1997; Borter, 1999; Fuchs et al, 2007) for debris flows, rock avalanches or debris avalanches. Landslide intensity significantly varies along the path of movement, or trajectory, and from one path to another for the same landslide size and source area. Moreover, landslides with

different intensity can converge from different sources. Hence, for a certain position, the intensity needs to be described as a distribution of values (Eq. (1)).

Onset frequency, fo, is a function of the volume of landslides (i.e., MF relationship) and it depends on geological, morphological, geomechanical and hydrological conditions controlling the stability of the slope. Propagation or transit frequency depends on the path of movement, or trajectory, and on where moving material will come to rest. These are controlled by several factors such as slope geometry, slope materials, slide volume, size and shape of blocks, and changes in material properties (e.g. water content, entrainment). Due to a large epistemic uncertainty about these factors, the simulation of landslide propagation is affected by a certain degree of uncertainty, which should be considered for the assessment of transit frequency. This can be done by using stochastic models that, given a landslide detachment, could simulate a number of possible paths. Moreover, the convergence of different paths to a certain point can be simulated by 3D models (Chen and Lee, 2003, 2007; Crosta et al., 2004; Hungr et al., 2007; Hungr, 2009). The modeled transit frequency, ft, should be calculated as the number of potential paths passing through a position, t, normalized by the total number of simulated paths from single or multiple sources, ttot:

f - (to-,

Fig. 1. Flow chart of hazard assessment methodology applied to rockfalls.

As a consequence, the frequency at a certain position, f, for long-runout landslides in (Eq. (2)) is:

f = fo ■ ft (6)

4. Methodology for rockfall hazard assessment

An interesting and complex example of the probabilistic landslide hazard analysis can be the case of rockfall process. Rockfalls are characterized by a high mobility, and the intensity calculated at each point along the slope is affected by a large uncertainty about the elasto-plastic properties of material at impact, the roughness of the slope, and the effect of vegetation. Therefore, they represent a class of phenomena, which are affected by most of the discussed aspects and can be a good example for testing the proposed landslide hazard assessment method. Starting form Frattini et al. (2012), we present a methodology that can incorporate both uncertainty and multiple sources. The methodology consists of three steps: (1) computation of onset annual frequency, (2) modeling of rockfall runout and spatially varying kinetic energy, and testing its distribution at each location along the slope and (3) performing a probabilistic analysis, which integrates over all rockfall volume scenarios to produce rockfall hazard curves including different probabilities for different block sizes and different source areas. These give an estimate of the probability of exceeding values of kinetic energy, at each location, within a reference time interval. Eventually hazard maps can be derived from hazard curves, choosing a fixed annual exceedance probability value.

The methodology for rockfall hazard analysis is summarized in Fig. 1.

4.1. Onset frequency, fo

Onset frequency of a rockfall, fo, is a function of magnitude, and can be expressed in terms of magnitude frequency relationships, MCF, where magnitude is expressed in terms of rockfall volume, v (Hungr et al. 1999; Dussauge et al. 2003). These relationships can be developed from surveyed historical events, or derived from the existing literature (e.g. Hungr et al. 1999; Chau et al., 2003; Guzzetti et al., 2003; Guthrie and Evans, 2004).

The occurrence rate for a given volume range can be derived by using:

f. _ NCV >v). ^ VQ -b

where N(V > v) is the cumulative number of events with volume larger than v for a catalogue of duration T; N is the total number of the observed events with volume larger than the lower bound V0. N0 depends on both the area extent and the overall susceptibility of the cliff, whereas the power law exponent, b, mainly depends on lithology and geological structure (Hungr et al., 1999).

For a rigorous analysis, it should be necessary to define an MCF for each potential source area, because this can vary according to local geologic and morphologic characteristics. Practically, this local analysis is generally impossible because the historical catalogues used to obtain MCF curves are defined for large areas with different geological and geomechanical settings as well as type and size of external perturbations in order to have a statistically significant number of events.

4.2. Transit frequency and hazard curve

For the assessment of transit frequency it is necessary to model rockfall dynamics either by using empirical methods (e.g. the "shadow cone" method, Evans and Hungr, 1993) or by simulating free-fall, impact, bouncing and rolling motions by the use of kinematic (e.g., Stevens, 1998), hybrid (Pfeiffer & Bowen, 1989; Jones et al. 2000; Crosta et al. 2004) and dynamic mathematical models (Descouedres and Zimmermann, 1987; Azzoni et al. 1995) in a 2D or 3D space.

The simplest intensity parameter to consider is total kinetic energy, Ek. Its distribution at each position on the slope (i.e. each cell) should be evaluated, in order to apply Eq. (1) in a correct form.

Table 1

Values of the model parameters for rockfall modeling on the artificial slope. en = normal restitution coefficient; et = tangential restitution coefficient; Tan ®r = rolling friction coefficient.

Model parameters

Block shape Radius et en

sphere 4.5 m 0.8 0.4 0.3

4.3. Distribution of intensity

The kinetic energy of a rockfall significantly varies along the trajectory, and from one path to another. Hence, for a given initial rockfall magnitude (i.e., volume), a continuous change of intensity (i.e., kinetic energy or dynamic pressure) along the path occurs. At the same time, due to the possible convergence of many trajectories in a specific place, different values of local intensity can be observed at a certain location. The uncertainty associated to runout is accounted in the model by using a stochastic approach for the damping parameters used in the simulations. At each impact and during rolling motion, the rockfall loses energy according to the values of restitution and friction coefficients, which are selected randomly from pre-defined stochastic distributions. The propagation of each block along the slope can be seen as a series of impacts and rolling, for which a stochastic variation is assigned to the coefficients of restitution and friction, using a normal distribution. Hence, the whole propagation is a sequence of products among the initial kinetic energy of the block and the restitution coefficients randomly sampled from a normal distribution. For the Central Limit Theorem, the product of random variables that are normally distributed approaches a log-normal distribution. Hence, from a theoretical point of view and for the simplest cases, we could expect the kinetic energy to be log-normally distributed. However, since these distributions in the proposed methodology have a strong importance, conditioning the correct use of Eq. (1) in the correct form, the hypothesis of log-normality of kinetic energies of blocks at a specific transit spot must be verified before the implementation of the probabilistic analysis.

• • * * Kinetic Energy (KJ) < 100,000 • 100,000 - 150,000 • 180,000 - 200,000 * > 200,000

* • # • •

• • • • • • * • • •

• • • • • •

¿3 •M •Î..I Hf.1 » • « irk • 4 • >

2 m 3 m

Horizontal width

Fig. 3. Kinetic energies and location of block impacts along the sampling plane # 5, two point source, transversal disposition.

Fig. 4. Percentages of sampling planes showing a log-normal distribution of kinetic energy of transiting blocks for different variations of restitution coefficients values and different source type: 1P = one point source, 2PH = two point sources horizontally distributed, 2PV = two point sources vertically distributed, L = linear source, A = areal source.

4.4. Parametric analysis

To verify the distribution of kinetic energy of blocks falling along a slope under controlled conditions, a set of rockfall simulation experiments was performed by means of the numerical model HY STONE (Crosta et al., 2004). The code provides kinetic energy and velocity of the blocks at each sampling point. It allows to simulate the motion of non-interacting rocky blocks in a three-dimensional framework. It is based on a hybrid (mixed kinematic-dynamic) algorithm with different damping relationships available to simulate energy loss at impact or by rolling. Slope topography is described by a raster DEM, and all the relevant parameters are spatially distributed. The stochastic nature of rockfall processes and the variability of parameters are introduced by random sampling most parameters from different probability density distributions (e.g. uniform, normal, lognormal, exponential). Block fragmentation and visco-plastic impact algorithms (Di Prisco and

0.01 -

1E-4 0.01

' ■ ....... ' 1 ....... 1 limit l . ; • T 500 years

f =0.1316 v ' °7 o • Tr 2500 years _ R2=0.993 ;

à f =0.263 V 0

Volume (m

Fig. 5.325 mapped blocks and rockfall sources in the area of Richmond Hill, Christchurch, N.Z. From Christchurch City Council.

Fig. 6. Magnitude cumulative frequency distributions (MCF) of block volumes fallen during 2011 Christchurch earthquake (from Aurecon, courtesy from CERA, Canterbury Earthquake Recovery Authority), with respect to 500 and 2500 years return time earthquakes. Probably, the smallest blocks were undersampled, such as the largest ones, due to block fragmentation. Blocks within the 0.5-10.5 m3 volume range are supposed to be correctly sampled. 950 m3 outlier block is not represented.

Vecchiotti, 2006) are also included in the code. Effects of impact on structures (i.e. buildings), nets (countermeasures) and trees can be evaluated. The rocky block is described as a solid geometric shape with a certain volume and a certain mass.

To test the distribution of kinetic energy of blocks, a simple artificial slope was considered. This 3 m x 3 m gridded slope is 620 m long, with a slope gradient of 37°, followed by a flat terrain sector (Figure 2). The use of such a simple slope allows to exclude the effects of: complex morphology (i.e. curvatures, convergences, slope changes), roughness variation, presence of vegetation in order to explore the effects of stochastic variation of restitution and friction coefficients, and the effect of convergence of trajectories of blocks from different sources. Different simulations were performed considering different source cells located at the top of the slope. 1000 blocks were launched from a single point source area, 500 blocks from two transversal or longitudinal point source areas, 100 blocks from each source cell of a linear and a polygonal source area (Figure 2, Table 1 ).

The distribution of kinetic energy values is analyzed at different positions, along 10m long transversal planes equally distributed along the slope every 60 m (Figures 2 and 3). Kolmogorov-Smirnov test was used to explore the log-normality of the results.

To explore the effect of the stochastic nature of the rockfall process and the uncertainty of the relevant parameters, different values of standard deviation (0.05, 0.2 and 0.5) were introduced in the experiments by randomly sampling from a normal distribution the restitution and rolling friction coefficient values.

To explore the effect of convergence of multiple trajectories, coming from different sources, the position and typology of the source area were varied. Several experiments have been performed for different source area configurations: a single point-source (Figure 2b), two transversally distributed point sources (Figure 2c), two longitudinally distributed point sources (Figure 2d), a linear source (Figure 2e), and an areal source (Figure 2f).

In most experiments, the distributions of kinetic energy values at each sampling plane resulted to be log-normal, as hypothesized before (Figure 4). The higher the variation of the restitution and friction coefficients, the higher the number of lognormal distributions. This confirms that a log-normal distribution derives from a stochastic process, as hypothesized before from a theoretical point of view.

When the source areas are linear or are represented by two points, the percentage of cells presenting a log-normal distribution of kinetic

Table 2

Rockfall volume scenarios and related frequencies of occurrence, fo, referring to earthquakes with 500 and 2500 years return time (Tr).

Scenarios 1 2 3 4 5 6

Range of 0.001-0.01 0.01-0.1 0.1-1 1-10 10-100 100-1000

volume (m3)

fo (Tr 500) 195.28 16.243 1.351 0.11237 0.009347 0.00077747

fo (Tr 2500) 39.42 3.279 0.273 0.02268 0.001887 0.000157

energies is higher. For a single-point source area, the behavior of trajectories is very similar, especially with a small variation of the parameters. Hence, the kinetic energy values tend to be almost constant. On the other side, areal sources introduce an effect of superimposition of numerous trajectories characterized by very different paths. This effect tends to produce a large variety of distribution shapes (e.g. Weibull, logistic, log-logistic).

In conclusion, the distribution of rockfall intensity falling along an ideal slope can be considered, with a good approximation, as lognormal, especially considering linear source areas.

In real slopes, further morphometric parameters and slope variability (e.g. heterogeneous materials, different land uses) could influence this result, which will have to be tested in the specific case.

In this case study, under the evidence that intensity approximately has a lognormal distribution, the equation for the probability density, Ps (lnEk) can be expressed as:

Ps( lnEk) =

( ln Ek-Ms)2

a sV2n

and the probability of exceeding each lnl, Ps(> ln I) as:

Ps(> ln Ek) =

a sV22n

( lnEk-ßs)2

e 2as d ln Ek

The total annual rate of exceeding each i, Ftot (I > i) is obtained summing over all N scenarios, according to Eq. (3).

5. Case study

On February 22, 2011, a 6.2 Mw earthquake took place in New Zealand, having its epicenter in a zone with high seismicity, approximately 6-10 km south-east of Christchurch. Hundreds of rockfalls were triggered in Port Hills (65 km2), around Christchurch, where the population density is almost 260/km2. At least five people were killed by falling rocks, and several hundred homes were evacuated. The specific seismic history, the slope morphology, the lithology, and the significant population density make of this area a potential rockfall high-risk hot spot. This claims for quantitative hazard assessment, which is necessary to evaluate risk for people and structures. The availability of quantitative data related to earthquakes, fallen block positions and volumes, allows a straightforward application of the proposed methodology for quantitative probabilistic hazard analysis.

Table 3

Parameters used in rockfall modeling for the Richmond Hill area, Christchurch, N.Z.

Model parameters

Block shape

Radius

Density

Sphere

6 scenarios (Table 2) 2700 kg/m3

Spatially distributed (Supplemental Figure 1s) Spatially distributed (Supplemental Figure 2s) Spatially distributed (Supplemental Figure 3s)

In this study, we test the methodology in the southern slope of Richmond Hill, which is located 6 km far from the epicenter, and covers an area of 0.6 km2. After the earthquake, 325 fallen blocks have been mapped along the slope (Christchurch City Council) (Figure 5) with volumes ranging from 0.027 to 15.8 m3, with an outlier block of 950 m3. By using this dataset, we derived the MCF curve (Figure 6) by fitting through Least Square Regression the cumulative frequency of blocks ranging in volume between 0.6 and 14 m3, obtaining a power-law exponent b = —1.07 (R2 = 0.993), which is similar to others found in the literature (Malamud et al, 2004). We fitted the frequency within this range because we observed a deviation from power-law for larger and smaller blocks (Massey et al., 2012).

The onset frequency was obtained by combining the size distribution of surveyed blocks with the Christchurch earthquake recurrence responsible for rockfall triggering. We assumed that the probability of rockfalls in the study area is controlled by earthquake occurrence probability, neglecting other triggering causes.

To assign a return time to the Christchurch earthquake, the response spectra measured at three recording stations located close to the study area and derived from EQC-GNS Geonet, were compared with the NZS1170 (2004) elastic design spectra (Kam and Pampanin, 2011). The Christchurch earthquake acceleration spectra lie between 500 and 2500 years return time. To account for the uncertainty related to the definition of the return time, these return times were used to calculate upper and lower rockfall hazard (Table 2 and Figure 6). For both of them, fo was calculated for six scenarios of block volume applying (Eq. (7)), with N° equal to 213 and 43 blocks for 500 and 2500 years as T, respectively. These upper and lower frequency bounds were calculated by dividing the total number of blocks obtained using the MCF curve for a volume ranging between 0.001 and 1000 m3 by the two return times. This volume range reflects the size of blocks mapped after the event in the entire Port Hills area (Massey et al., 2012). In the frequency calculation, we assume that only earthquakes with the magnitude of the event (Mw = 6.2) can trigger rockfalls. This could induce an underestimation of the upper frequency bound (N>=213 blocks), since rockfalls can also be triggered by different magnitude earthquakes, either smaller or bigger (Keefer, 1984; Massey et al., 2012). However, the upper frequency bound is calculated with a conservative return time of 500 years, thus partially compensating this underestimation. In fact, the expected number of blocks here calculated is larger than in Massey et al. (2012).

The frequency, fs, and the spatially variable distribution of kinetic energy were obtained through 3D runout simulations by using HY-STONE code (Crosta et al., 2004; Agliardi et al, 2009) with a 1-m gridded LiDAR DEM. As source areas, the rocky cliffs characterized by outcropping basaltic lavas have been considered. A simulation was performed for each of the 6 volume scenarios, launching 20 blocks from each source cell (1,529 cells). Modeling parameters have been back-calibrated from the 2011 event (Table 3, Figures 1s, 2s and 3s).

For each scenario, s, and for each grid cell of the model, HY-STONE provides raster and vector outputs. Among the raster outputs, the number of block transits per cell was used to calculate the transit frequency, ft. The vector outputs consist in points with values of rotational, transla-tional and total kinetic energy, velocity, and height, sampled at a distance of 5 m.

In order to obtain for each position in the slope a probability density of kinetic energy to be used in Eq. (7), we subdivided the study area in 5 x 5 m squares, and we sampled, within each square, the computed values of total kinetic energy. To further investigate the distribution of the kinetic energy and to support the previously discussed assumption of log-normality, we tested this hypothesis for several squares randomly selected along the slope by using Kolmogorov-Smirnov's test (Figure 7 and Supplemental Figure 4s). In more than 67% of the cases, the best fit of the data was obtained by a lognormal distribution. Exceptions are mainly located in areas which are of secondary importance for hazard assessment, (e.g. very close to cliffs, or in distal portions with

Fig. 7. a) Kinetic energy values sampled along the trajectories for the simulation scenario with blockvolume ranging from 100 to 1000 m3; b) examples of ln(Eknn distribution (squares 1 to 4). Probability density is shown in y axis. Further distributions are provided in Suppl. Fig. 4s.

fewer trajectories). For these areas, no dominant distribution can be found for block kinetic energies. This supports the adoption of lognormal distribution thus allowing the use of Eq. (1) in its lognormal form (Eqs. (8) and (9)). Hence, for each square, the arithmetic mean and the standard deviation of the logarithm of kinetic energy have been calculated (Figure 8).

By using Eqs. (8) and (9), we calculated, for each position along the slope, the hazard curve, where the exceedance probability in 50 years is expressed as a function of rockfall kinetic energy, Ek. By assuming an earthquake return time of 500 years, we obtain slightly higher hazard curves with respect to a return time of 2500 years (Figure 9).

By selecting specific values of probability of exceedance (e.g. 10% in 50 years) we sampled the corresponding value of kinetic energy from each of the hazard curves. This represents the kinetic energy which is exceeded with that probability, and it is referred to each specific 5 m x 5 m square, with a specific location on the slope. Mapping kinetic energy values fully probabilistic hazard maps are obtained, to be used for land planning and risk management (Figures 10 and 11, Figures 5s and 6s). As expected, hazard maps derived for a 500 years return time earthquake show, for a given exceedance probability, higher values of kinetic energy with respect to a return time of 2500 years. The modal values of kinetic energy with an

a e„<j) b ejj)

ln<Ek) ln(Ek)

Fig. 9. Examples of hazard curves for 10% exceedance probabilities in 50 years, corresponding to earthquake with return time 500 (a) and 2500 years (b). Kinetic energy values exceeded with 10% exceedance probabilities in 50 years are extracted in correspondence of the black dotted line. Each colored line represents a 5 m x 5 m square hazard curve.

exceedance probability of 10% in 50 years are about 60 kJ and 8 kJ for 500 and 2500 years return time, respectively (Figure 11).

6. Discussion

With respect to other quantitative hazard assessment approaches (e.g. Hovius et al, 1997; Hungr et al. 1999; Stark and Hovius, 2001; Dussauge et al, 2003; Guthrie and Evans, 2004; Malamud et al., 2004; Marchi and D'Agostino, 2004; Jakob and Friele, 2010; Santana et al,

2012), the probabilistic landslide hazard analysis here proposed explicitly expresses hazard as a function of landslide destructive power, considering landslide intensity rather than magnitude. This supports the use of vulnerability and fragility functions, and the assessment of risk. With a few exceptions (Straub and Schubert, 2008; Spadari et al.,

2013), the above mentioned approaches do not account explicitly for

uncertainty. For intensity, simple statistics of the expected value are used (e.g. arithmetic average, maximum value, or a specific percentile), neglecting the dispersion and introducing strong assumptions about its distribution. In the proposed methodology, we explicitly account for the uncertainty by using the whole intensity distribution in the hazard analysis. This also allows to obtain hazard curves, which explicitly represent the probability of exceeding a certain level of landslide intensity within a defined time period, accounting for uncertainty and integrating different magnitude scenarios.

The proposed probabilistic landslide hazard analysis is a flexible approach, applicable to all typologies of landslides, if necessary accounting for both their onset and transit probability. However, its mathematical formulation has to be declined according to the distribution of landslide intensity values, which has to be verified in each specific case.

Fig. 10. Rockfall hazard maps with 10% exceedance probabilities in 50 years, for an earthquake with 500 years (a), and 2500 years (b) return time.

7 50 400 3 10a 2.2 10' 1.6 10s 1.2 10s

tr 500

in (Ek)

Fig. 11. Frequency of hazard values of hazard maps with 10% exceedance probabilities in 50 years, for earthquakes with 500 and 2500 years return time.

In general, the requirements for the successful application of the probabilistic landslide hazard analysis are:

• the frequency of landslide scenarios, usually as a function of landslide magnitude;

• an appropriate measure of landslide intensity, which depends on the type of phenomenon (e.g. kinetic energy, displacement, velocity, depth);

• an appropriate probability density distribution of intensity; the general methodological framework can be applied with different intensity probability density distribution, modifying Eq. (1). The probability density distribution could vary according to the landslide typology and the geological and morphological context. For rockfall, a lognormal probability density function seems suitable. For slow moving landslides, where the intensity distribution is controlled only by uncertainty, a normal probability density function should be suitable.

Landslide hazard curves and maps allow a quantitative comparison with the hazard connected with other threats existing in the study area.

As an example, we quantitatively compare rockfall and seismic hazards in the Christchurch area by assessing the expected level of damage with a 10% exceedance probability in 50 years. For rockfalls, the expected level of damage induced by a 1 m radius block potentially impacting a two storey reinforced-concrete building is calculated by using fragility curves proposed by Mavrouli and Corominas, (2010) for slight, moderate, extensive and complete damage levels. For the same type of building, the expected damage due to seismic hazard is calculated by using NZS1170 (2004) earthquake actions for Christchurch, considering a building natural period ranging from 0.2 to 0.33 s (Filiatrault et al.,

2009), and New Zealand fragility curves (Uma et al., 2008). The expected levels of damage for the two hazards are similar, and amount to about 10% of the value of the building (Table 4). In particular, among the 30 buildings potentially impacted by rockfalls, the majority shows higher levels of expected damage due to earthquakes, and only 8 due to rockfalls (Figure 12). This would suggest that seismic hazard is more relevant than rockfall hazard, except for few buildings lying below the most dangerous sector of the rocky cliff. However, the analysis was performed considering reinforced concrete buildings, which are not common in the study area. This is the only building type for which fragility curves for rockfalls are available. For more common woodframe buildings, the fragility to rockfall would be larger with respect to earthquakes, leading to different expected damages.

A validation of the hazard maps obtained by the applications of the methodology here presented is quite hard to be performed directly. It should require the availability of kinetic energy data referred to fallen blocks, and should be supported by the presence of quantitative data related to damages and specific fragility curves. Generally, both these types of data are extremely difficult to be collected. However, the probabilistic landslide hazard analysis can be validated in its different steps: onset and propagation. In the proposed case study, a validation of the runout modeling has been performed through the localization of fallen blocks and the verified compatibility with the trajectories simulated by the runout model, and through the number and position of effectively damaged buildings.

7. Conclusions

The hazard curves are calculated for each position along the slope, with a spatially distributed approach, leading to the possibility to realize hazard maps, which represent, for a defined exceedance probability, the spatial distribution of intensity along the slope. This opens great possibilities in the use of these hazard maps for land planning and management, since they allow both the recognition of areas impacted by landslide with different level of intensity, given a certain probability, and the presence of hazard hot spots.

Since the intensity in hazard curves is defined as a measure of damage potential, it can be used directly in vulnerability or fragility curves for a probabilistic risk analysis, in which the choice of an appropriate time period over which the probability of exceedance is calculated (e.g., 1 year, 50 years) may be guided by the lifetime of exposed elements for which the analysis is performed. Probabilistic risk analysis allows a cost-benefit analysis based prioritization of mitigation strategies aimed at minimizing damages and danger for people, buildings and infrastructures. Moreover, risk analysis could be extremely useful for insurance purposes.

Beyond the hazard assessment purposes, a probabilistic analysis of landslide processes can be used for landscape evolution models. For rockfalls, a probabilistic analysis can help in modeling the accumulation of blocks over time, thus simulating lowering of rocky cliffs and the formation of talus or debris fans.

Table 4

Application of vulnerability curves to the earthquake and to a 1 m radius rockfall, for a two storey reinforced-concrete building, considering a building natural period ranging from 0.2 to 0.33. P(D) is the probability of exceeding each damage level DL, and D is the contribution of each damage level to the total expected damage, which results to be 0.0906 and 0.0968 for the earthquake and for the rockfall, respectively.

Damage level DL Earthquake Rockfall

% loss P(D) D P(D) D

DL1 slight 0.02 0.4 0.0046 1 0.0158

DL2 moderate 0.1 0.17 0.006 0.21 0.006

DL3 extensive 0.5 0.11 0.03 0.15 0.075

DL4 complete 1 0.05 0.05 0 0

0.0906 0.0968

Acknowledgements

The authors would like to thank Hannes Salzmann (Free Fall Geotechnical Engineering) and CERA, Canterbury Earthquake Recovery Authority for providing Christchurch rockfall data, Federico Agliardi for stimulating discussion about the methodology, Elena Valbuzzi and Stefania Zenoni for the preparation of database and rockfall inventory. The research has been partially funded by the PRIN-MIUR 2010-2011 - 2010E89BPY_007 and SafeLand — Living with landslide risk in Europe: Assessment, effects of global change, and risk management strategies, 7th Framework Programme.

Fig. 12. Location of buildings with a higher level of expected damage due to rockfalls and to earthquakes.

Appendix A. Supplementary data

Supplementary data to this article can be found online at http://dx. doi.org/10.1016/j.enggeo.2014.07.015.

References

Agliardi, F.,Crosta, G.B.,Frattini, P., 2009. Integrating rockfall risk assessment and counter-measure design by 3D modelling techniques. Nat. Hazards Earth Syst. Sci. 9, 1059-1073. http://dx.doi.org/10.5194/nhess-9-1059-2009.

AGS, 2007. Guidelines for landslide susceptibility, hazard and risk zoning for land use management. Australian geomechanics society landslide taskforce landslidingzoning working group. Aust Geomech. 42 (1), 13-36.

Annaka, T., Satake, K., Sakakiyama, T., Yanagisawa, K.,Shuto, N., 2007. Logic-tree approach for probabilistic tsunami hazard analysis and its applications to the Japanese coasts. In: Basel, Birkhäuser (Ed.), Tsunami and Its Hazards in the Indian and Pacific Oceans, pp. 577-592 http://dx.doi.org/10.1007/s00024-006-0174-3.

Archetti, R, Lamberti, A., 2003. Assessment of risk due to debris flow events. Nat. Hazards Rev. 4 (3), 115-125. http://dx.doi.org/10.1016/0148-9062(95)00018-C.

Azzoni, A., La Barbera, G., Zaninetti, A., 1995. Analysis and prediction of rock falls using a mathematical model. Int. J. Rock Mech. Min. Sci. Geomech. Abstr. 32, 709-724. http://dx.doi.org/10.1016/0148-9062(95)00018-C.

Bird, J.F.,Crowley, H.,Pinho, R.,Bommer, J.J., 2005. Assessment of building response to liquefaction induced differential ground deformation. Bull. N. Z. Soc. Earthq. Eng. 38 (4), 215-234.

Borter, P., 1999. Risikoanalyse bei gravitativen Naturgefahren. Bezugsquelle: Bundesamt für Umwelt, Wald und Landschaft, Dokumentation.

Bovolin, V.,Taglialatela, L., 2002. Proceedings Scritti in onore di Lucio Taglialatela, CNR-GNDCI Pubblicazioni no. 2811, Naples, Italypp. 429-437.

Calvo, B.,Savi, F., 2009. A real-world application of Monte Carlo procedure for debris flow risk assessment. Comput. Geosci. 35 (5), 967-977. http://dx.doi.org/10.1016/jxageo. 2008.04.002.

Cancelli, A.,Crosta, G., 1993. Hazard and risk assessment in rockfall prone areas. In: Skipp, B.O. (Ed.), Risk and Reliability in Ground Engineering. Thomas Telford, Springfield, pp. 177-190.

Chau, K.T.,Wong, R.H.C.,Liu, J., Lee, C.F., 2003. Rockfall hazard analysis for Hong Kong based on rockfall inventory. Rock Mech. Rock. Eng. 36 (5), 383-408. http://dx.doi. org/10.1007/s00603-002-0035-z.

Chen, H.,Lee, C.F., 2003. A dynamic mode for rainfall-induced landslides on natural slopes. Geomorphology 51, 269-288. http://dx.doi.org/10.1016/S0169-555X(02)00224-6.

Chen, H., Lee, C.F., 2007. Landslide mobility analysis using MADFLOW. Proc The 2007 International Forum on Landslide Disaster Management. , 2, pp. 857-874.

Cornell, C.A., 1968. Engineering seismic risk analysis. Bull. Seismol. Soc. Am. 58 (5), 1583-1606.

Corominas,J.,Copons, R.,Moya, J.,Vilaplana,J.M.,Altimir, J.,Amigo, J., 2005. Quantitative assessment of the residual risk in a rockfall protected area. Landslides 2 (4), 343-357. http://dx.doi.org/10.1007/s10346-005-0022-z.

Corominas, J., Moya, J., 2008. A review of assessing landslide frequency for hazard zoning purposes. Eng. Geol. 102 (3), 193-213. http://dx.doi.org/10.1016/j. enggeo.2008.03.018.

Corominas, J., van Westen, C.,Frattini, P.,Cascini, L.,Malet, J.-P.,Fotopoulou, S.,Catani, F.,Van Den Eeckhaut, M.,Mavrouli, O., Agliardi, F.,Pitilakis, K., Winter, M.G.,Pastor, M.,Ferlisi, S.,Tofani, V.,Hervas, J.,Smith, J.T., 2014. Recommendations for the quantitative analysis of landslide risk. Bull. Eng. Geol. Environ. 73 (2), 209-263.

Crosta, G.B., Agliardi, F., Frattini, P.,Imposimato, S., 2004. A three dimensional hybrid numerical model for rockfall simulation. Geophys. Res. Abstr. 6,04502.

Crovelli, R.A., 2000. Probability models for estimation of number and costs of landslides. US Geological Survey.

Crozier, M.J., 1997. The climate landslide couple: a Southern Hemisphere perspective. Paleoclim. Res. 19, 333-354.

Cruden, D.M., 1997. Estimating the risk from landslide historical data. In: Cruden, D.M., Fell, R. (Eds.), Landslide Risk Assessment. Proc. Inter. Workshop on Landslide Risk Assessment, Honolulu, 19-21 February 1997. Balkema, Rotterdam, pp. 177-184.

Dai, F.C.,Lee, C.F., 2001. Frequency-volume relation and prediction of rainfall-induced landslides. Eng. Geol. 59 (3), 253-266. http://dx.doi.org/10.1016/ S0013-7952(00)00077-6.

Del Gaudio, V.,Pierri, P.,Wasowski, J., 2003. An approach to time-probabilistic evaluation of seismically induced landslide hazard. Bull. Seismol. Soc. Am. 93 (2), 557-569. http ://dx.doi.org/10.1785/0120020016.

Descouedres, F., Zimmermann, T., 1987. Three-dimensional dynamic calculation of rockfalls. Proceedings of the Sixth International Congress of Rock Mechanics, Montreal, Canada, pp. 337-342.

Di Prisco, C., Vecchiotti, M., 2006. A rheological model for the description of boulder impacts on granular strata. Geotechnique 56 (7), 469-482.

D'Odorico, P., Fagherazzi, S., Rigon, R., 2005. Potential for landsliding: dependence on hyetograph characteristics. J. Geophys. Res. Earth Surf. 110 (F1). http://dx.doi.org/ 10.1029/2004JF000127 (2003-2012).

Dussauge, C., Grasso, J.R., Helmstetter, A., 2003. Statistical analysis of rockfall volume distributions: implications for rockfall dynamics. J. Geophys. Res. 108 (B6), 2286. http://dx.doi.org/10.1029/2001JB000650.

Evans, S.G.,Hungr, O., 1993. The assessment of rockfall hazard at the base of talus slopes. Can. Geotech. J. 30,620-636. http://dx.doi.org/10.1139/t93-054.

Fell, R.,Corominas, J.,Bonnard, Ch.,Casrini, L.,Leroi, E., Savage, W.Z., on behalf of theJTC-1 Joint Technical Committee on Landslides and Engineered slopes, 2008. Guidelines for landslide susceptibility, hazard and risk zoning for land-use planning, Commentary. Eng. Geol. 102, 99-111.

Filiatrault, A., Christovasilis, I.P., Wanitkorkul, A., van de Lindt, J.W., 2009. Experimental seismic response of a full-scale light-frame wood building. J. Struct. Eng. 136 (3), 246-254. http://dx.doi.org/10.1061 /(ASCE)ST.1943-541X.0000112.

Frankel, A., Mueller, C.,Barnhard, T., Perkins, D., Leyendecker, E.,Dickman, N., Hanson, S., Hopper, M., 1996. National seismic hazard maps: documentation June 1996. USGS Open-File Report 96-532 (70 pp.).

Frattini, P.,Crosta, G.B.,Allievi, J., 2013. Damage to buildings in large slope rock instabilities monitored with the PSInSAR™ technique. Remote Sens. 5 (10), 4753-4773. http://dx. doi.org/10.3390/rs5104753.

Frattini, P., Crosta, G.B., Lari, S., Agliardi, F., 2012. Probabilistic rockfall hazard analysis (PRHA). In: Eberhardt, E., Froese, C., Turner, A.K., Leroueil, S. (Eds.), Landslides and Engineered Slopes: Protecting Society Through Improved Understanding. Taylor & Francis, London, pp. 1145-1151.

Frattini, P.,Crosta, G.B.,Sosio, R., 2009. Approaches for defining thresholds and return periods for rainfall-triggered shallow landslides. Hydrol. Process. 23 (10), 1444-1460. http://dx.doi.org/10.1002/hyp.7269.

Friele, P., Jakob, M., Clague, J., 2008. Hazard and risk from large landslides from Mount Meager volcano, British Columbia, Canada. Georisk 2 (1), 48-64. http://dx.doi.org/ 10.1080/17499510801958711.

Fuchs, S.,Heiss, K.,Hübl, J., 2007. Towards an empirical vulnerability function for use in debris flow risk assessment. Nat. Hazards Earth Syst. Sci. 7 (5), 495-506. http://dx.doi. org/10.5194/nhess-7-495-2007.

Geist, E.L., Parsons, T., 2006. Probabilistic Analysis of Tsunami Hazards. Nat. Hazards 37 (3), 277-314. http://dx.doi.org/10.1007/s11069-005-4646-z.

Geist, E.L., Parsons, T., 2008. Distribution of tsunami interevent times. Geophys. Res. Lett. 35 (2), L02612. http://dx.doi.org/10.1029/2007GL032690.

Geist, E.L.,Chaytor, J.D., Parsons, T.,ten Brink U., 2013. Estimation of submarine mass failure probability from a sequence of deposits with age dates. Geosphere 9 (2), 287-298. http://dx.doi.org/10.1130/GES00829.1.

Gentile, F.,Bisantino, T.,Liuzzi, G.T., 2008. Debris-flow risk analysis in south Gargano watersheds (Southern-Italy). Nat. Hazards 44 (1), 1-17. http://dx.doi.org/10.1007/ s11069-007-9139-9.

González, F.I., Geist, E.L.,Jaffe, B.JKanoglu, U.,Mofjeld, H.,Synolakis, C.E.,Titov, V.V., Arcas, D., Bellomo, D., Carlton, D., Horning, T., Johnson, J., Newman, J., Parsons, T., Peters, R., Peterson, C., Priest, G.,Venturato, A., Weber, J., Wong, F.,Yalciner, A., 2009. Probabilistic tsunami hazard assessment at Seaside, Oregon, for near and far field seismic sources. J. Geophys. Res. Oceans 114 (C11). http://dx.doi.org/10.1029/2008JC005132 (1978-2012).

Grünthal, G.,Thieken, A.H.,Schwarz, J.,Radtke, K.S.,Smolka, A.,Merz, B., 2006. Comparative risk assessments for the city of Cologne—storms, floods, earthquakes. Nat. Hazards 38 (1-2), 21-44. http://dx.doi.org/10.1007/s11069-005-8598-0.

Gutenberg, B.,Richter, C.F., 1942. Earthquake magnitude, intensity, energy and acceleration. Bull. Seismol. Soc. Am. 32,163-191.

Guthrie, R.H., Evans, S.G., 2004. Analysis of landslide frequencies and characteristics in a natural system, coastal British Columbia. Earth Surf. Process. Landf. 29 (11), 1321-1339. http://dx.doi.org/10.1002/esp.1095.

Guzzetti, F.,Reichenbach, P.,Wieczorek, G.F., 2003. Rockfall hazard and risk assessment in the Yosemite Valley, California, USA. Nat. Hazards Earth Syst. Sci. 3 (6), 491-503.

Hovius, N., Stark, C.P., Allen, P.A., 1997. Sediment flux from a mountain belt derived by landslide mapping. Geology 25, 231-234. http://dx.doi.org/10.1130/0091-7613.

Hungr, O., 1981. Dynamics of rock avalanches and other types of mass movements(Ph. D. thesis) Univ. of Alberta, Edmonton, Canada (500 pp.).

Hungr, O., 1997. Some methods of landslide hazard intensity mapping. Landslide risk assessment. Balkema, Rotterdam, pp. 215-226.

Hungr, O., 2009. Numerical modelling of the motion of rapid, flow-like landslides for hazard assessment. KSCE J. Civ. Eng. 13 (4), 281-287. http://dx.doi.org/10.1007/s12205-009-0281-7.

Hungr, O., Evans, S.G., Hazzard, J., 1999. Magnitude and frequency of rock falls and rockslides along the main transportation corridors of southwestern British Columbia. Can. Geotech. J. 36, 224-238. http://dx.doi.org/10.1139/t98-106.

Hungr, O.,KcKinnon, M.,McDougall, S., 2007. Two models for analysis of landslide motion: application to the 2007 Hong Kong benchmarking exercises. Proc. The 2007 International Forum on landslide disaster management. vol. 2, pp. 919-932.

Jaboyedoff, M.,Dudt, J.P.,Labiouse, V., 2005. An attempt to refine rockfall hazard zoning based on the kinetic energy, frequency and fragmentation degree. Nat. Hazards Earth Syst. Sci. 5 (5), 621-632. http://dx.doi.org/10.5194/nhess-5-621-2005.

Jaiswal, P., van Westen, C.J.,Jetten, V., 2011. Quantitative assessment of landslide hazard along transportation lines using historical records. Landslides 8 (3), 279-291. http://dx.doi.org/10.1007/s10346-011 -0252-1.

Jakob, M., 2005. Asize classification for debris flows. Eng. Geol. 79 (3), 151-161. http://dx. doi.org/10.1016/j.enggeo.2005.01.006.

Jakob, M., Friele, P., 2010. Frequency and magnitude of debris flows on Cheekye River, British Columbia. Geomorphology 114 (3), 382-395. http://dx.doi.org/10.1016/j. geomorph.2009.08.013.

Jones, C.L., Higgins, J.D., Andrew, R.D., 2000. Colorado Rockfall Simulation Program Version 4.0. Colorado Department of Transportation, Colorado Geological Survey (March, 127 pp.).

Kam, W.Y., Pampanin, S., 2011. General Performance of Buildings in Christchurch CDB After the 22 Feb 2011 Earthquake: a Contextual Reportprepared for the Department of Building and Housing, University of Canterbury (62 pp.).

Keefer, D.K., 1984. Landslides caused by earthquakes. Geol. Soc. Am. Bull. 95,406-421.

Keylock, C.J., McClung, D.M., Mar Magnusson, M., 1999. Avalanche risk mapping by simula-tion.J. Glaciol. 45 (150), 303-314. http://dx.doi.org/10.3189/002214399793377103.

Keylock C.J., Barbolini, M., 2001. Snow avalanche impact pressure-vulnerability relations for use in risk assessment. Can. Geotech. J. 38 (2), 227-238. http://dx.doi.org/10. 1139/cgj-38-2-227.

Lambert, C., Thoeni, K., Giacomini, A., Casagrande, D., Sloan, S., 2012. Rockfall hazard analysis from discrete fracture network modelling with finite persistence discontinuities. Rock Mech. Rock. Eng. 45 (5), 871-884. http://dx.doi.org/10.1007/ s00603-012-0250-1.

Liu, Y.,Santos, A., Wang, S.M.,Shi, Y.,Liu, H.,Yuen, DA, 2007. Tsunami hazards along Chinese coast from potential earthquakes in South China Sea. Phys. Earth Planet. Inter. 163 (1), 233-244. http://dx.doi.org/10.1016/j.pepi.2007.02.012.

Malamud, B.D.,Turcotte, D.L.,Guzzetti, F.,Reichenbach, P., 2004. Landslide inventories and their statistical properties. Earth Surf. Process. Landf. 29 (6), 687-711. http://dx.doi. org/10.1002/esp.1064.

Mansour, M.F.,Morgenstern, N.R.,Martin, C.D., 2011. Expected damage from displacement of slow-moving slides. Landslides 8 (1), 117-131. http://dx.doi.org/10.1007/s10346-010-0227-7.

Marchi, L.,D'Agostino, V., 2004. Estimation of debris-flow magnitude in the eastern Italian Alps. Earth Surf. Process. Landf. 29, 207-220. http://dx.doi.org/10.1002/esp.1027.

Massey, C.I.,McSaveney, M.J.,Heron, D.,Lukovic, B., 2012. Canterbury earthquakes 2010/11 Port-Hills slope stability: pilot study for assessing life-safety risk from rockfalls (boulder rolls). GNS Science Consultancy Report, 311, p. 115.

Mavrouli, O.,Corominas, J., 2010. Rockfall vulnerability assessment for reinforced concrete buildings. Nat. Hazards Earth Syst. Sci. 10 (10), 2055-2066. http://dx.doi.org/10. 5194/nhess-10-2055-2010.

Moon, A.T., Olds, R.J., Wilson, R.A., Burman, B.C., 1992. Debris flow zoning at Montrose, Victoria. Proceedings of 6th International Symposium on Landslides, vol. 2, pp. 1015-1022.

Negulescu, C.,Foerster, E., 2010. Parametric studies and quantitative assessment of the vulnerability of a RC frame building exposed to differential settlements. Nat. Hazards Earth Syst. Sci. 10,1781-1792. http://dx.doi.org/10.5194/nhess-10-1781-2010.

NZS 1170.5, 2004. Structural design actions. Part 5: Earthquake Actions New ZealandStandards New Zealand, Wellington, NZ.

OFAT, OFEE, OFEFP, 1997. Recommendations: prise en compte des dangers dux aux meuvements du terrain dans le cadre des activités de l'aménagement du territoire. p. 42 (Berne, Switzerland).

Ogata, Y., 1999. Estimating the hazard of rupture using uncertain occurrence times of paleoearthquakes. J. Geophys. Res. 104,17,995-18,014. http://dx.doi.org/10.1029/ 1999JB900115.

Parsons, T., 2008. Monte Carlo method for determining earthquake recurrence parameters from short paleoseismic catalogs: example calculations for California. J. Geophys. Res. 113 (B3), B03302. http://dx.doi.org/10.1029/2007JB004998.

Pfeiffer, T.J.,Bowen, T., 1989. Computer simulation of rockfalls. Bull. Assoc. Eng. Geol. 26, 135-146.

Picarelli, L., 2011. Discussion to the paper "Expected damage from displacement of slow-moving slides" by MF Mansour, NR Morgenstern and CD Martin. Landslides 8 (4), 553-555. http://dx.doi.org/10.1007/s10346-010-0227-7.

Picarelli, L., Oboni, F., Evans, S.G., Mostyn, G., Fell, R., 2005. Hazard characterization and quantification. In: Hungr, O., Fell, R., Couture, R., Eberhardt, E. (Eds.), Landslide Risk Management. Balkema Publishers, London, UK, pp. 27-62.

Pierson, L.A.,Davis, S.A., Van Vickle, R., 1990. Rockfall Hazard Rating System: Implementation Manual (No. FHWA-OR-EG-90-01).

Rathje, E.M.,Saygili, G., 2008. Probabilistic seismic hazard analysis for the sliding displacement of slopes: scalar and vector approaches. J. Geotech. Geoenviron. 134 (6), 804-814. http://dx.doi.org/10.1061/ASCE1090-02412008134:6804.

Rosso, R.,Rulli, M.C.,Vannucchi, G., 2006. A physically based model for the hydrologic control on shallow landsliding. Water Resour. Res. 42 (6). http://dx.doi.org/10.1029/ 2005WR004369.

Salciarini, D., Godt, J.W., Savage, W.Z., Baum, R.L., Conversini, P., 2008. Modeling landslide recurrence in Seattle, Washington, USA. Eng. Geol. 102 (3), 227-237. http://dx.doi. org/10.1016/j.enggeo.2008.03.013.

Santana, D., Corominas, J.,Mavrouli, O.,Garcia-Sellés, D., 2012. Magnitude-frequency relation for rockfall scars using a Terrestrial Laser Scanner. Eng. Geol. 145, 50. http://dx. doi.org/10.1016/j.enggeo.2012.07.001. Schuster, R.L., Wieczorek, G.F., 2002. Landslide triggers and types. Landslides: Proceedings of the First European Conference on Landslides. Taylor & Francis, Prague, pp. 59-78.

Shih, B.J.,Shieh, C.L.,Chen, L.J., 1997. The grading of risk for hazardous debris-flow zones. Debris-Flow Hazards Mitigations Mechanics, Prediction, and Assessment. ASCE, pp. 219-228. Sidle, R.C., Pearce, A.J., O'Loughin, C.L., 1985. Hillslope Stability and Land Use. Water Resources Monograph no. 11American Geophysical Union, Washington (DC). Spadari, M.,Kardani, M.,De Carteret, R.,Giacomini, A.,Buzzi, O.,Fityus, S., Sloan, S.W., 2013. Statistical evaluation of rockfall energy ranges for different geological settings of New South Wales, Australia. Eng. Geol. 158, 57-65. http://dx.doi.org/10.1016/j.enggeo. 2013.03.007.

Stark, C.P.,Hovius, N., 2001. The characterization of landslide size distributions. Geophys.

Res. Lett. 28 (6), 1091-1094. http://dx.doi.org/10.1029/2000GL008527. Straub, D., Schubert, M., 2008. Modeling and managing uncertainties in rock-fall hazards. Georisk2 (No. 1).

Stevens, W., 1998. RocFall: a tool for probabilistic analysis, design of remedial measures and prediction of rockfalls(M.A.Sc. Thesis) Department of Civil Engineering, University of Toronto, Ontario, Canada (105 pp.).

SYNER-G project, 2011. Systemic Seismic Vulnerability and Risk Analysis for Buildings, Lifeline Networks and Infrastructures Safety Gain". Deliverable 3.1: fragility functions for common RC building types in Europe (March, 209 pp.).

Tatard, L.,Grasso, J.R.,Helmstetter, A., Garambois, S., 2010. Characterization and comparison of landslide triggering in different tectonic and climatic settings. J. Geophys. Res. 115, F04040. http://dx.doi.org/10.1029/2009JF001624.

Uma, S.R.,Bothara, J.,Jury, R.,King, A., 2008. Performance assessment of existing buildings in New Zealand. Proceedings of the New Zealand Society for Earthquake Engineering Conference, Wairakei, New Zealand, 11-13 April, p. 45.

Witt, A.,Malamud, B.D., Rossi, M.,Guzzetti, F., Peruccacci, S., 2010. Temporal correlations and clustering of landslides. Earth Surf. Process. Landf. 35 (10), 1138-1156. http:// dx.doi.org/10.1002/esp.1998.