Probabilistic assessment of seismic-induced slope displacements: an application in Italy

Earthquake-induced slope instability is one of the most important hazards related to ground shaking, causing damages to the environment and, often, casualties. Therefore, it is important to assess the seismic performance of slopes, especially in the near fault regions, evaluating the permanent displacements induced by seismic loading. This paper applies a probabilistic approach to evaluate the seismic performance of slopes using an updated database of ground motions recorded during the earthquakes occurred in Italy. The main advantage of this approach is that of accounting for the aleatory variability of both ground motions and prediction of seismic-induced displacements of slopes. The results are presented in terms of hazard curves, showing the annual rate of exceedance of permanent slope displacement evaluated using ground motion data provided by a standard probabilistic hazard analysis and a series of semi-empirical relationships linking the permanent displacements of slopes to one or more ground motion parameters. The procedure has been implemented on a regional scale to produce seismic landslide hazard maps for the Irpinia district, in Southern Italy, characterised by a severe seismic hazard. Seismic landslide hazard maps represent a useful tool for practitioners and government agencies for a regional planning to identify and monitor zones that are potentially susceptible to earthquake-induced slope instability, thus requiring further detailed, site-specific studies.


Introduction
Earthquakes often produce instability in natural slopes, causing severe human and economic losses. Therefore, several efforts have been devoted in the past decades by the geotechnical community to study the seismic response of slopes aimed at mitigating and preventing such catastrophic events.
The inertial forces induced by seismic loading generate permanent displacements due to the temporary mobilisation of the shear strength during earthquake loading, with accumulation of plastic strains within the sliding surface. More times the shear strength is attained during a severe seismic event, the higher is the permanent sliding experienced at the end of the earthquake. Therefore, the seismic performance of a slope can be described by the permanent displacements developed at the end of the seismic event, which can be quantified through several methods. Stress-deformation modelling (e.g. finite element, discrete element, material point methods) provides the most accurate prediction of the seismic response of slopes as it takes into account some relevant characteristics of soil behaviour, such as the cyclic degradation of shear strength parameters (Bandini et al. 2015) and the development of excess pore-water pressures during earthquake (Di Filippo et al. 2019). However, effort and time involved in the analyses and the necessity of accurate site characterisation and geotechnical investigations reserve these methods to the study of critical projects. Furthermore, the analyses are site-specific and cannot be employed to evaluate slope stability at a regional scale. Conversely, the well-known displacement-based Newmark's method (Newmark 1965) is preferred in common practice as it represents a good compromise between simplicity and accuracy of the results. In the method, the slope is assimilated to a rigid block sliding on a horizontal plane that undergoes permanent displacements only when the acceleration of the input motion is greater than a critical value, the latter depending on slope resistance. Although the Newmark's method ignores the deformation of the sliding mass during seismic shaking, the hypothesis of rigid block is applicable for shallow slides and falls in rock and debris (Jibson 2011).
Many researchers have proposed predictive equations based on Newmark-type computations, using various ground motion databases. Indeed, the integration of a set of acceleration time histories permits to link the permanent displacements to the yield seismic coefficient, representing the slope seismic resistance, and to a series of ground motion parameters (e.g. Del Gaudio et al. 2003;Bray and Travasarou 2007;Jibson 2007;Saygili and Rathje 2008;Fotopoulou and Pitilakis 2015;Song et al. 2017;Tropeano et al. 2017;Bray and Macedo 2019;Cho and Rathje 2020). For instance, Chousianitis et al. (2014) adopted the Arias intensity as descriptor of ground motion as it better characterises the damaging effects induced by the earthquakes. Even  and Cho & Rathje (2020) demonstrate that for predictive equations based on a single ground motion parameter (scalar approach), the smallest standard deviation is obtained using the Arias intensity I A as it encompasses the intensity, the frequency content and the duration of the ground motion. However, the performance of the semi-empirical relationships increases when more than one ground motion parameter are used in the equation (Bazzurro & Cornell 2002) as they can capture additional significant features of the ground motion that affect the slope displacement, typically related to the intensity and frequency content of the input motion. Furthermore, Bray et al. (2018) highlighted the importance of differentiating displacement models depending on the tectonic setting. Gaudio et al. (2020) derived semi-empirical relationships for the Italian territory by performing a parametric integration of an updated version of the Italian seismic database, along the line tracked by Rampello et al. (2010) and Biondi et al. (2011). The database includes 954 records of 208 seismic events with moment magnitude M w ≥ 4, peak ground accelerations PGA ≥ 0.05 g, epicentral distance R ep < 100 km recorded by 297 stations located throughout the national territory referred to a temporal period ranging from 14/06/1972 to 24/04/2017, thus including the seismic records of the recent destructive earthquakes occurred in 2009 L'Aquila (Maugeri et al. 2011), 2012Emilia (Mucciarelli and Liberatore 2014 and 2016 Central Italy seismic sequence (Mollaioli et al. 2019;Luzi et al. 2019). For each station, information about the subsoil class is available according to the classification provided by Eurocode 8, Part 1 (CEN 2003) and by the Italian Building Code (Ministero delle Infrastrutture 2018). Gaudio et al. (2020) also provide information on the tectonic setting of the seismic database: 581 events (61.3% of the total) are characterised by a normal focal mechanism (NF), while the registrations with reverse (TF), oblique reverse (U) and strike-slip (SS) fault mechanisms are 170 (18%), 126 (13.3%) and 65 (6.9%), respectively. Therefore, based on the number of records, the most relevant mechanism in this study is the normal fault. The parametric analyses were performed for different values of the yield seismic coefficient k y using both one (scalar) and two (vector) ground motion parameter empirical models obtained adapting the ones recently proposed by Cho and Rathje (2020) to the Italian seismicity.
In this work, the semi-empirical relationships are further specialised distinguishing the subsoil classes of the recording stations: the registrations are 123 (13%) on rock-like (class A) subsoil, 469 (49.5%) on stiff (class B) subsoil and 294 (31%) on soft (class C) subsoil. Classes D and E are not taken into account in the following due to the small number of records available: 14 (1.5%) and 47 (5%), respectively.
The semi-empirical relationships are useful for a preliminary estimate of the expected seismic-induced displacement of slopes but they cannot account for the aleatory variability of earthquake ground motion and displacement prediction. Therefore, many researchers developed fully probabilistic-based approaches stemming from the above semi-empirical relationships to introduce the concept of hazard associated to the calculated displacement (e.g. Rathje and Saygili 2008, 2011Bradley 2012;Lari et al. 2014;Du and Wang 2016;Rodriguez-Marek and Song 2016;Macedo et al. 2018;Cho and Rathje 2020). The results of the probabilistic approach are typically illustrated in terms of displacement hazard curves, which provide the mean annual rate of exceedance λ d (or the return period T r = 1∕λ d ) for different levels of permanent displacements. In addition to the mentioned works, the probabilistic approach has been applied to the study of real slopes in different seismic areas throughout the world (e.g. Gülerce and Balal 2017;Wang and Rathje 2018), while Du et al. (2018) investigated the influence of the slope properties on the seismic displacement analysis. Within the same theoretical framework, the sliding block procedures and the probabilistic approaches have been also modified to account for the deformability of the sliding mass, to investigate more realistically the seismic performance of deep sliding mechanisms (e.g. Saygili and Rathje 2008;Rathje and Antonakos 2011;Song and Rodriguez-Marek 2015;Bray et al. 2018). Furthermore, in some cases the probabilistic approaches have been combined in a logic tree scheme to account for the epistemic uncertainty of the hazard (e.g. Rathje and Saygili 2009;Wang and Rathje 2018;Macedo et al. 2020).
The probabilistic approach can be profitably implemented on a regional scale using the ground motion hazard information to evaluate landslides hazard maps Wang and Rathje 2015;Chousianitis et al. 2016;Sharifi-Mood et al. 2017). These maps permit to identify the zones that are more susceptible to earthquake-induced slope instability and to evaluate the probability of occurrence of a displacement level in a specific time interval.
The present study aims to contribute towards the evaluation of earthquake-induced landslides hazard for Italy, which is characterised by a high level of seismicity. The displacement hazard curves and maps are developed for the Italian territory to provide a useful tool for practitioners and government agencies for a preliminary evaluation of the seismic performance of slopes, stemming from the updated seismic database and the results of the semiempirical relationships adopted in the study. Both scalar and vector probabilistic approaches are employed: in the first case the peak ground acceleration (PGA) is used as ground motion parameter, while the combination of PGA with the peak ground velocity (PGV) is chosen for the vector approach. Although Cho and Rathje (2020) and Gaudio et al. (2020) demonstrate that the best double-parameters models are obtained for the couples Arias intensity and mean period (I A , T m ) and (I A , PGV), the couple of parameters PGA, PGV is here adopted because: (i) the ground motion prediction equations (sometimes denoted in the literature as attenuation laws) entering the probabilistic method are generally available only for PGA and PGV; and (ii) the standard seismic probabilistic hazard analysis (PSHA) is commonly developed in terms of peak ground acceleration.
The paper is organised as follows: the probabilistic approach proposed by Rathje and coworkers is firstly briefly summarised, highlighting the key ingredients and the procedure of analysis. The analyses results are then shown in terms of displacement hazard curves for different Italian sites; the effect of soil class and yield seismic coefficient of the slope on the results are discussed and the differences between the scalar and vector approaches are critically analysed.
Different ground motion prediction equations are used to obtain the displacement hazard curves for the vector approach and their influence on results is discussed. Finally, a series of hazard maps for the Irpinia district, in the Southern Italy, are illustrated for some threshold values of permanent displacement d y and different yield seismic coefficients k y .

Background: probabilistic approach
The displacement models predict the natural log of permanent horizontal displacement d as a function of the natural log of one or more ground motion parameters (GM) or earthquake magnitude.
For the single ground motion parameter model, the expression proposed by Ambraseys and Menu (1988) is adopted: where a 0 , a 1 and a 2 are regression coefficients and d is expressed in cm. The efficiency of the semi-empirical relationships is described by the standard deviation of the natural log of displacement σ ln (Cornell and Luco 2001). The advantage of Eq. (1) with respect to polynomial forms is that of complying with the conditions d → ∞ for k y /PGA = 0 and d = 0 for k y /PGA = 1 expected for the assumption of a rigid sliding block.
The results of a probabilistic analysis are synthesised in terms of seismic displacement hazard curves, providing the mean annual rate of exceedance λ d for different levels of displacements given a yield seismic coefficient and a specific site. According to Cho and Rathje (2020) and Rathje et al. (2014), for the single ground motion PGA displacement model the annual rate of exceedance λ d can be computed as: where P d > x|PGA i is the probability of the displacement exceeding a specific value x, given a peak ground acceleration PGA i , and P PGA i is the annual probability of occurrence of the ground motion level PGA i . The evaluation of the annual rate of exceedance thus requires the integration of the product of these probabilities over all possible levels of PGA for the specific site. The first term is evaluated assuming a cumulative lognormal distribution function for displacement values x, where the mean value is provided by Eq. (1) and σ ln is the standard deviation. The second term is derived from differentiation of the PGA hazard curve obtained through a probabilistic seismic hazard analysis (PSHA). In detail it is (Rathje et al. 2014): where λ is the mean rate of occurrence associated with adjacent PGA values in the hazard curve. Therefore, the probabilistic approach predicts the annual rate of exceedance λ d of different levels of permanent displacements, similarly to the annual rate of exceedance λ for different levels of PGA, obtained by a PSHA. Ambraseys and Menu (1988) suggest to modify Eq. (1) to account for the effects of other ground motion parameters. Therefore, along this track, the following semi-empirical relationship is proposed for the case of two ground motion parameters: Equation (4) represents the key ingredient for the development of the vector approach. For the case of the couple of ground motion parameters (PGA, PGV), Rathje et al. (2014) provided the following expression for the evaluation of λ d : where P d > x|PGA i , PGV j is the probability of the displacement exceeding a specific value x, given the peak ground acceleration PGA i and the peak ground velocity PGV j , and P PGA i , PGV j is the joint annual probability of occurrence of ground motion levels PGA i and PGV j . Again, the first term is calculated assuming a cumulative lognormal distribution function for displacement values x, with the mean value and the standard deviation σ ln provided by Eq. (4). The double summation represents the integration over all levels of PGA and PGV. Although the joint annual probability should be rigorously computed via a vector probabilistic seismic hazard analysis VPSHA (Bazzurro and Cornell 2002), here the simplified approach is adopted to evaluate the former from the output obtained from a conventional PSHA: The term P M k , R l |PGA i represents the probability of occurrence of different earthquake magnitudes M and distances R calculated from the disaggregation of the hazard of PGA while P PGV j |PGA i , M k , R l is computed assuming a lognormal distribution, with the mean value μ ln PGV|PGA, M, R and the standard deviation σ ln PGV|PGA,M,R obtained as: The terms μ ln PGV , σ ln PGV , μ ln PGA and σ ln PGA are the mean values and the standard deviations of the ground motion parameters PGA, PGV and ρ is the correlation coefficient between PGA and PGV. The first quantities are calculated from the most suitable ground motion prediction equation (GMPE) depending on the specific site, whereas for the correlation coefficient the procedure illustrated by Rathje and Saygili (2008) is adopted: where ε PGA , ε PGV , ε PGA and ε PGV are the normalised residuals for the ground motion parameters PGA, PGV and the mean value on the total number of observations, computed as: where ln GM observed is the natural log of the observed ground motion and ln GM predicted is the natural log of the ground motion parameter predicted by the GMPE.

Results of the probabilistic approach
The results of the probabilistic approach can be synthesised in terms of displacement hazard curves, plotting the mean annual rate of exceedance λ d for different levels of permanent displacements and in terms of hazard maps, showing the distribution of the return period (6) T r = 1∕λ d for a given permanent displacement. Hazard curves and maps are developed for different zones of the Italian territory and the main factors affecting the results are critically analysed. The results have been obtained using the ground motion database described by Gaudio et al. (2020), that provide details on the seismic events considered in the analyses. According to Gaudio et al. (2020), the permanent displacements are computed with the rigid sliding-block model (Newmark 1965) for the simple scheme of an infinite slope assuming a constant shear strength, that is a constant k y , during earthquake loading. Computation have been carried out for different values of the yield seismic coefficient, using an updated version of the Italian strong motion database. Four values of k y = 0.08, 0.1, 0.12, 0.15 are taken into account as they represent the average values of the yield seismic coefficient for the scheme of an infinite slope characterised by angle of slope inclination in the range of 5-20° and angles of shearing resistance φ′ = 22°-28°. Table 1 reports the regression coefficients of the semi-empirical relationships evaluated for computed displacements greater than 0.0001 cm. Despite very small values of displacements are meaningless from an engineering perspective, they have been considered to develop semi-empirical relationships capable to fit satisfactory the data as the PGA approaches k y . The coefficients of the semi-empirical relationships are further specialised to account for the influence of the subsoil class on the probabilistic approach.
The permanent displacements computed through the Newmark's method using the whole database are plotted in Fig. 1 versus PGA, for four values of the seismic coefficient k y .
The grey lines represent the predictions of the scalar semi-empirical relationship of Eq. (1) specialised for the different values of the yield seismic coefficient. It is seen that the adopted model not only predicts correctly the conditions at the extrema but also nicely reproduces the non-linear variation of the permanent displacements with PGA.
The computed permanent displacements and the semi-empirical relationships are plotted in Fig. 2 against the ratio k y /PGA for the whole database and the three subsoil classes (A, B and C). The vector predictive equation is plotted for two values of PGV within the range associated to the recorded ground motions: PGV = 10 and 30 cm/s. As reported in Table 1, the use of the two parameters semi-empirical relationship reduces significantly the standard deviation associated to the Newmark's displacements as compared to the scalar approach: an expected result as the couple of ground motion parameters PGA, PGV are more representative of the strong motion database that the PGA only. Furthermore, to assess the quality of the predicted equations, the residuals of  Figures 3 (a-c) show the residuals against k y /PGA while Fig. 3 (b-d) show the mean of residuals for different k y /PGA bins and the standard deviation of the mean values. It is seen that the vector approach predicts about halved values of the residuals with respect to the scalar one, with a trend characterised by negligible biases for increasing values of k y /PGA and no relevant variation of the standard deviations with k y /PGA. To better compare the efficiency of the scalar and vector models, Fig. 4 shows the variation of the standard deviation with k y /PGA: σ ln d is significantly reduced when the PGA, PGV model is used, its values being nearly constant with k y /PGA.
The residuals of the displacements plotted against other relevant ground motion parameters show negligible bias, especially for the vector approach, hence demonstrating that the regressions adopted in this study are appropriate for the Italian ground motion database. For the sake of conciseness, the results are reported in the Appendix (Fig. 21).

Displacement hazard curves
The displacement hazard curves represent a very attractive tool in the common practice to obtain a preliminary indication of the seismic performance of a slope. In fact, they allow to evaluate the frequency of occurrence of a given displacement for a slope in a specific  1) and (4), which provide the mean value and the standard deviation of the displacement d.
In addition, Eqs. (2) and (5) require specific site information of the seismic hazard. Therefore, the results of the probabilistic approach depend on the slope characteristics, via the yield seismic coefficient, and the collected strong motion database, as well as on the seismic hazard of the site at hand. Figure 5 shows the PGA hazard curve and the displacement hazard curve obtained through the scalar probabilistic approach for the Italian site of Amatrice, afflicted by the 2016 Central Italy seismic sequence, assuming k y = 0.1 . The PGA hazard curves have been extracted from the INGV interactive seismic hazard maps (http:// esse1. mi. ingv. it/ d2. html) that permits to query probabilistic seismic hazard maps of the Italian territory on a regular grid spaced by 0.05°. The PGA hazard curves are provided for rock-like (class A) subsoil with reference to the 16th, 50th and 84th percentile only, representing multiple fractiles of the hazard curve to account for the epistemic uncertainty in the ground motion hazard.
The computed displacement hazard curves show a decreasing annual rate of exceedance with increasing permanent displacements, as expected, with the dashed lines showing the effect of the epistemic uncertainty on the curve. In the following, the mean hazard curves are considered in computations. With reference to the mean hazard curve, for threshold displacement d y = 2 cm, related to a rock-like subsoil, d y = 5 cm, related to a brittle soil behaviour, and d y = 15 cm, related to a free-field ductile soil behaviour (Idriss 1985;Wilson and Keefer 1985), return periods T r = 1∕λ d of about 850, 2250 and 12,500 years are obtained. Figure 6a-c show the displacement hazard curves computed using the scalar approach with k y = 0.1 for three Italian areas: Central Italy, Emilia (Northern Italy) and Irpinia (Southern Italy).
Despite the significant number of sites investigated, a relatively narrow variability of the displacement hazard curves is observed, so that homogeneous macro-areas can be recognised in terms of seismic hazard, with lower and upper bounds of the displacement  Fig. 6d. For a given displacement level d, the annual rate of exceedance is greater for Irpinia and Central Italy than Emilia as the latter is characterised by a lower seismic hazard. In other words, the probability of occurrence of a permanent sliding displacement in two similar slopes located in two different sites solely differs for the site-specific seismic hazard.
According to the Newmark's sliding block analysis, the seismic resistance of a slope, described by the yield seismic coefficient, is a crucial ingredient for the evaluation of the permanent displacements of a slope subjected to an earthquake. As a consequence, Eqs. (1) and (4) explicitly depend on the values of k y adopted for the integration of the acceleration time histories, thus affecting the results of the probabilistic analysis in terms of displacement hazard curves.
This aspect is illustrated in Fig. 7 for the site of Amatrice, stemming from the scalar probabilistic approach: the curves are associated to four different values of the yield seismic coefficient by virtue of the dependence of Eq. (1) on k y . The more stable is the slope (i.e. greater is k y ), the lower is the annual rate of exceedance associated to a given displacement (e.g. Rathje and Saygili 2009;Macedo et al. 2018).
The results discussed till now refer to the signals recorded at all the stations located in the Italian territory, irrespective of the subsoil class associated to each station. Figure 8 illustrates the displacement hazard curves computed with the scalar approach for the Amatrice site, distinguishing the three subsoil classes and using k y = 0.1 . It is seen that higher annual rate of exceedance λ d is computed for soft (class C) subsoil as compared to that of stiff (class B) subsoil or, in other terms, for class C sites the model predicts the largest displacements at fixed rate of exceedance. This is consistent with the fact that the ground motions recorded on soft soils are richer in low frequencies than those recorded on stiff soils, hence leading to larger displacements. It is less clear the trend observed for the case of rock-like (class A) subsoil, for which the smallest displacements are expected, while the obtained hazard curve lies between those computed for classes B and C. This could be attributed to the fact that for rocklike subsoils about 50% of the records with PGA > 0.2 g are characterised by high values of the mean period T m > 0.4 s, that could explain the unexpectedly high displacements computed for subsoil class A, leading the hazard displacement curve to the right of that obtained for subsoil class B (see Figs. 17,18,19,20 in the Appendix). The displacement hazard curve obtained for the whole database is also plotted in Fig. 8. The displacement hazard curves can be also computed through the vector probabilistic approach using the couple of ground motion parameters PGA, PGV. The PGA hazard curve and the disaggregation information for probability of exceedance equal to 2, 5, 10,22,30,39,50,63, 81% for the Italian territory, used in Eqs. (3) and (6), are available at the INGV interactive website (http:// esse1. mi. ingv. it/ d2. html). The disaggregation is provided for each point of the grid for epicentral distance R in the range of 0 to 200 km and magnitude M = 3.5-9. The ground motion prediction equation of Lanzano et al. (2019) is adopted in this work as it improves the predictive capability for the Italian seismicity with respect to the GMPE proposed by Sabetta and Pugliese (1996) and Bindi et al. (2010). It provides the median values and the standard deviation used in Eq. (7) and leads to a correlation coefficient ρ = 0.843 calculated through Eq. (8). Figure 9 compares the displacement hazard curves computed with the scalar and the vector probabilistic approaches for the site of Amatrice using a yield seismic coefficient k y = 0.1 and distinguishing the three subsoil classes. The vector approach always predicts smaller values of λ d in the whole range of permanent displacements considered in the analysis, thus reducing the hazard estimate. For instance, with reference to case of all subsoil classes together, for a 10% probability of exceedance in 50 years (return period T r = 475 years) the displacements are equal to about 1.0 cm and 0.65 cm using the scalar and the vector approach, respectively, whereas for a 2% probability of exceedance in 50 years (return period T r = 2475 years) displacements of about 5.2 cm and 4.0 cm are obtained. This result is consistent with the outcomes of Rathje and Saygili (2008) and Cho and Rathje (2020). In fact, the semi-empirical relationship accounting for PGA and PGV leads to smaller standard deviation and median displacements as compared with the PGA prediction displacement equation, as reported in Table 1.
Similarly to the scalar approach, Fig. 10 shows the displacement hazard curves computed with the vector approach for the site of Amatrice using a yield seismic coefficient k y = 0.1 and distinguishing the three subsoil classes. Again, the curve obtained for the class C subsoil leads to larger displacements at fixed annual rate of exceedance as compared to that of stiff soils (class B). A less clear trend is still observed for the class A subsoil, as  mentioned above. Nevertheless, the effect of the subsoil classes on the hazard curves provided by the vector approach is less pronounced than that observed for the scalar approach.
As for the scalar approach, the displacement hazard curves obtained from a vector probabilistic analysis depend on the seismic resistance of the slope via the critical seismic coefficient k y and on the site seismicity. However, the GMPE represents an additional feature with respect to the single parameter approach and should be selected on the base of the seismicity of the investigated site. This is why the law proposed by Lanzano et al. (2019) is adopted in this work, that accounts for the relevant earthquakes occurred in Italy during the last decade. Nevertheless, to investigate the role of the ground motion prediction equation  on the displacement hazard curves, the probabilistic analyses have been also performed using other GMPEs: the one proposed by Sabetta and Pugliese (1996), also derived for Italian seismicity and, for comparison, the well-known predictive equations of Boore and Atkinson (2008) and Campbell and Bozorgnia (2008). Table 2 summarises the values of the correlation coefficient calculated through Eq. (8).
The correlation coefficients assume positive values, proving that the ground motion parameters PGA and PGV are positively correlated, in agreement to what obtained by Rathje and Saygili (2008), and show small changes among the different prediction equations, the higher value of ρ being provided by the Lanzano et al. (2019) GMPE.
The GMPE also affect the values of μ ln PGV|PGA, M, R and σ ln PGV|PGA, M, R provided by Eq. (7). Figure 11 shows the displacement hazard curve computed for the site of Amatrice and k y = 0.1 using the same GMPEs listed in Table 2, demonstrating that the influence of different ground motion prediction equations is negligible.

Displacement hazard maps
Hazard maps are more effective than the hazard curves to identify the hazard of earthquake-induced landslides for a specific area. The probabilistic approach is used in the following to develop the hazard maps for the district of Irpinia (Campania), in the South Italy, about 50 km East of Naples. The area has an extension of about 40 × 40 km 2 , is largely mountainous and centred on the section of the Apennines, a zone of Italy characterised by a severe seismic hazard and affected by many severe earthquakes that induced instability phenomena in several natural slopes (Porfido et al. 2002;Del Gaudio and Wasowski 2004).

3
The seismic hazard of the area under study is plotted in Fig. 12 in terms of distribution of PGA for the return period of 475 years (10% of probability of exceedance in 50 years). Information pertaining to the seismicity of the zone have been extracted from the INGV interactive seismic hazard maps (http:// esse1. mi. ingv. it/ d2. html) and refer to the rock-like (class A) subsoil.
The displacement model used for the hazard maps is based on a generic site class, considering the whole set of records of the area irrespective of the subsoil class associated to each station. In fact, due to the uncertainty related to the variability of soils at the regional scale, it is more convenient to consider the displacement model accounting for all the subsoil classes.
The hazard maps depict the contours of the return periods T r associated to different levels of seismic-induced displacement and yield seismic coefficient. It is worth mentioning that the maps do not account for the morphology of the region and the real distribution of the landslides. However, they are provided for different values of k y , so that they can be used for a preliminary evaluation of the seismic hazard of the region for different slope scenarios, each one characterised by a different value of k y . The hazard maps are developed stemming from both the scalar and the vector probabilistic approaches. Figure 13 shows the contours of the return period obtained from the two approaches for a threshold displacement d y = 15 cm and k y = 0.1 . The Cartesian coordinates of the map are East and North according to the reference coordinate system WGS84. The return periods are smaller for the scalar than the vector approach, again demonstrating that the vector probabilistic approach reduces the displacement hazard levels because of the more complete information provided by the combination of multiple ground motion parameters. Figure 14 shows the contour maps obtained from the vector approach for a threshold displacement d y = 2 cm and different values of the yield seismic coefficient. As expected, the return periods associated to a given value of d y increases with k y . The variability of T r is Fig. 12 Ground motion hazard maps for the Irpinia district: distribution of PGA for T r = 475 years clearly associated to the PGA hazard curves and disaggregation information of the region, that is more severe for the zone of the Apennines extending from North-Western to South-Eastern corners of the map. Therefore, the contour lines are nearly parallel to the NW-SE diagonal of the map, with larger return periods in the South-Western area, where hills take place of the mountains.
Similar results are shown in Figs. 15 and 16 for permanent threshold displacements equal to 5 and 15 cm. In agreement to what discussed for the hazard curves in the previous section, the larger is the admissible level of displacement, the larger is the related return period.

Summary and conclusions
The increasing adoption of probabilistic methods in engineering practice experienced in last decades has been the incentive to revisit the problem of seismic slope stability and highlight the advantages of this methodology to assess the seismic performance of natural slopes.
In this paper, the probabilistic approach formalised by Rathje and co-workers is adopted to evaluate the permanent slope displacement induced by earthquake loading for the Italian territory.
The probabilistic analyses are performed stemming from the earthquake-induced displacements evaluated through a Newmark-type computation under the simple hypothesis of a rigid sliding block, employing the Italian strong motion database updated by Gaudio et al. (2020). The method benefits from the large amount of recently available earthquake data and thus better characterises the variability in the strong ground motion. The results of the probabilistic approach are first summarised in terms of displacement hazard curves, showing the annual frequency of exceedance associated to the seismic-induced slope displacements.
Two probabilistic methods have been used: a scalar approach, characterised by the single ground motion parameter PGA and a vector approach, in which the ground motion is described by the couple of parameters PGA, PGV. The semi-empirical relationships (a) (b) Fig. 13 Hazard maps for the Irpinia district for d y = 15 cm and k y = 0.1: a scalar approach; b vector approach adopted in this study, linking the slope displacements to the ground motion parameters mentioned above, permit to evaluate the earthquake-induced displacements for different values of the yield seismic coefficient k y and different subsoil classes. The displacements provided by the PGA, PGV model are characterised by a lower standard deviation and lower residuals than the scalar one, also proving that the biases of the residuals with some significant ground motion parameters are negligible. In both cases, seismic hazard information pertaining to the specific site is a key ingredient requiring a standard probabilistic hazard analysis, whose results are available throughout the Italian territory and synthesised in terms of PGA hazard curves and disaggregation information. The influence of the site-specific seismic hazard information on the results of the probabilistic analysis have been first investigated with reference to the scalar approach. The displacement hazard curves, computed for three different sites in Italy, demonstrate that it is possible to identify homogeneous macro-areas in terms of displacement hazard curves, highlighting, as expected, that a more severe seismic hazard increases the hazard associated to the permanent slope displacements. In addition, with reference to a specific site, the hazard reduces with the yield seismic coefficient that represents the seismic (c)

(a) (b)
Fig. 14 Vector approach: hazard maps for the Irpinia district for d y = 2 cm: a k y = 0.1, b k y = 0.12, c k y = 0.15 resistance of a slope. In principle, the probabilistic approach does not account for the subsoil class. In this paper, the influence of site conditions on the hazard curves has been investigated referring the semi-empirical relationships for seismic-induced slope displacements to the subsoil class of the recording stations of the Italian strong motion database. It has been shown that the hazard increases for the soft (class C) subsoil due to the lower frequency content of the ground motions, while it decreases for signals recorded on stiff soils. The results of the scalar and the vector approaches have been also compared, showing that the latter provides a lower estimate of the seismic hazard, as a result of the more reliable description of input motion through two ground motion parameters. Indeed, the scalar approach, despite its simplicity, cannot account for the complexity of the input motion in that the ground motion parameter PGA alone is not sufficient to describe the main features of a seismic event. Use of the couple PGA, PGV parameters also showed to be less sensible to site conditions, in that the displacement hazard curves computed for subsoil of class A, B and C plotted in a narrower band than that observed for the scalar approach. Therefore, for what mentioned above, the vector approach should be preferred to the scalar one to compute the seismic hazard maps at a regional scale. The proposed probabilistic approach has been applied to the study area of Irpinia district, located in the Southern Italy and characterised by a high seismic hazard, computing seismic landslide hazard maps that provide the return periods associated with different values of threshold earthquake-induced displacements and seismic yield coefficient. Consistently with the results obtained for the hazard curves, the seismic hazard maps show increasing return periods as the threshold displacements and the seismic yield coefficient increase.
The results presented in this paper show that the probabilistic approach represents a useful tool for a preliminary evaluation of the seismic performance of slopes. The displacement hazard curves provide site-specific information, whereas the hazard maps are suitable for the analysis and the evaluation of the seismic risk of natural slopes at a regional scale. The probabilistic nature of these maps also enables them to be combined with other hazards for a more complete hazard analysis, useful for disaster management to minimise damages caused by earthquake-induced landslides.
The hypothesis of rigid sliding block should be intended as appropriate for the study of shallow translational mechanisms. Nevertheless, the procedure can be extended to account for the deformability of the sliding mass. Other possible aspects that could be considered to improve the accuracy of seismic landslide maps are the spatial variability of soil properties and the topographic amplification of ground shaking, thus improving further the predictive capability of seismic landslide hazard maps.

Appendix
The seismic database were divided into five groups depending on the subsoil of the recording stations, according to the classification provided by Eurocode 8, Part I (CEN 2003): rock-like subsoil (class A), dense coarse-gained and stiff fine-grained subsoils (class B), medium dense coarse-grained and medium stiff fine-grained subsoil (class C), loose coarse-grained and soft fine-grained subsoil (class D) and subsoils with stiffness of class C or D underlain by stiffer material (class E). Figure 17 shows the relative frequency distribution of M w , R ep , PGA, PGV, T m and D 5-95 for the subsoil classes A, B and C. The records on classes D and E were neglected as they represent only 1.5% and 5% of the total strong motions, respectively.
The datasets are similarly distributed over the considered ground motion parameters with a lower number of records for subsoil class A, especially for values of PGA and PGV greater than 0.2 g and 10 cm/s, respectively. The distribution shows that values of PGA and PGV are strongly correlated, with most frequent peak ground acceleration and velocity falling between 0.05 and 0.1 g and between 1 and 5 cm/s, respectively.
The relationships between different ground motion parameters are distinguished for the three subsoil classes A, B and C. The moment magnitude M w is plotted in Fig. 18 against the epicentral distance, while Figs. 19, 20 show the peak ground velocity PGV and Arias intensity I A , as well as the mean period T m and spectral acceleration S a (1.5T s ) against the peak ground acceleration PGA. As expected, PGV, I A and S a (1.5T s ) increase with increasing PGA, while no remarkable influence is observed for T m . However, for subsoil class A it is seen that the 50% of the records with values of PGA > 0.2 g are characterised by high values of the mean period T m > 0.4 s that could explain why the semi-empirical relationships provide unexpectedly high displacements for subsoil class A, leading the hazard displacement curve to the right of that computed for subsoil class B. Figure 21 shows the variation of the residuals obtained from the one-parameter (PGA) and the two-parameters (PGA, PGV) models with M w , PGV, S a (1.5T s ) and R ep ; the grey lines in the figures represent the mean of residuals. For all the considered ground motion parameters, the vector approach provides much lower residuals with nearly constant mean values, proving that the semi-empirical relationship of Eq. (4) is suitable for the application to the Italian territory.