Probabilistic Approach to Provide Scenarios of Earthquake-Induced Slope Failures (PARSIFAL) Applied to the Alcoy Basin (South Spain)

: The PARSIFAL (Probabilistic Approach to pRovide Scenarios of earthquake-Induced slope FAiLures) approach was applied in the basin of Alcoy (Alicante, South Spain), to provide a comprehensive scenario of earthquake-induced landslides. The basin of Alcoy is well known for several historical landslides, mainly represented by earth-slides, that involve urban settlement as well as infrastructures (i.e., roads, bridges). The PARSIFAL overcomes several limits existing in other approaches, allowing the concomitant analyses of: (i) ﬁrst-time landslides (due to both rock-slope failures and shallow earth-slides) and reactivations of existing landslides; (ii) slope stability analyses of different failure mechanisms; (iii) comprehensive mapping of earthquake-induced landslide scenarios in terms of exceedance probability of critical threshold values of co-seismic displacements. Geotechnical data were used to constrain the slope stability analysis, while speciﬁc ﬁeld surveys were carried out to measure jointing and strength conditions of rock masses and to inventory already existing landslides. GIS-based susceptibility analyses were performed to assess the proneness to shallow earth-slides as well as to verify kinematic compatibility to planar or wedge rock-slides and to topples. The experienced application of PARSIFAL to the Alcoy basin: (i) conﬁrms the suitability of the approach at a municipality scale, (ii) outputs the main role of saturation in conditioning slope instabilities in this case study, (iii) demonstrates the reliability of the obtained results respect to the historical data.


Introduction
Landslides are one of the most common collateral effect of earthquakes [1]. Seismic action is, indeed, one of the main triggers to be considered for landslide processes [2], especially in the case of their reactivation [3].
Damages caused by earthquake-induced ground effects can exceed the ones directly linked to seismic shaking [4]. Moreover, earthquake-induced landslides represent the result of a "chain effect" [5] which links two natural events (earthquake and landslide) having an independent probability of occurrence (in the case of landslide reactivation) or a conditioned one (in the case of a first-time landslide). As a consequence, the risk can be amplified because of the combined induced damage

•
Step 1: Slope Analysis (SA). This step consists in the identification of slope sectors prone to failure that can be referred to first-time landslides or reactivations of already existing landslide masses, regardless the dynamic input. This step includes a susceptibility analysis for identifying potentially unstable zones (i.e., zones prone to failure for first-time earth slides, rock slides and toppling) and a geomorphological analysis to correctly locate and classify, in terms of failure mechanism, the already existing landslides. The final aim is then to perform a zonation of the study area based on the proneness to a specific failure mechanism. In case of first-time earth slides, areas of outcropping debris cover and colluvial deposits (i.e., the deposits affected by shallow landslides) are subdivided into areas featured by the unique combination of lithology and slope and aspect classes (Unique Condition Units-UCU). The proneness to failure of each UCU is then assessed to select those to be analyzed in detail in the following step. In case of rock slides or toppling, the study area is partitioned in Homogenous Geo-structural zones (HGZ), i.e., zones having similar lithology, morpho-structural setting and geo-structural features (in terms of number and geometrical attitude of discontinuity sets). The kinematic compatibility to wedge, toppling and planar sliding of each elementary unit (raster cell) forming the HGZ is performed to select the cell units to be further analyzed in the following step. In case of already existing landslides (i.e., landslide mass reactivations) they can be identified and inventoried in a GIS geodatabase and analyzed according to their characteristic mechanism (i.e., mainly rotational or translational).

•
Step 2: Slope Stability (SS). This step consists in the assessment of the actual potential of failure/reactivation under combined seismic and hydraulic loading of the areas resulted unstable in the previous phase. Specifically, this step aims at the computation of an exceedance probability of co-seismic displacements based on an assumed threshold for failure (P[D ≥ Dc|a(t), a y ]) and performed through the Newmark's method [22]. In case that the threshold is not exceeded, a safety factor (SF) is computed to provide a stability level of the slope. To perform such an analysis the requested input data consist of: seismic acceleration, response spectra and selected time-histories related to the most severe seismogenic source among the ones that contribute to the seismic hazard of the site. Different conditions of soil/rock-joints saturation can be accounted for.

•
Step 3: Resulting Scenario (RS). This step aims at providing, for each considered scenario of combined seismic and hydraulic conditions, a synthetic and integrated GIS-based graphical representation (mapping) of the slope stability results for both first-time slope failures and landslide reactivations. Such a mapping shows the different territorial units defined as elementary areas able to represent the slope stability results in terms of either exceedance probability or safety factor. In case of rock slopes, these areas are represented by single cells of a raster, while vectorial data are used to represent slope stability results for first time earth-slides and already existing landslides.

Case Study
Alcoy is a Spanish municipality of 70,000 inhabitants, belonging to the Comunidad Valenciana in the province of Alicante, located in the south-eastern part of Spain.
The municipality area is about 130 km 2 extended and lies within a geological basin known as "Hoya de Alcoy". The urban texture shows a close connection with the morphological evolution suffered by this area, where the rivers Polop, Bartxell and Molinar, tributaries of the river Serpis, played a main role, modeling valleys with steep slopes that had a remarkable influence on the urban development over time.
Alcoy is characterized by an internal Mediterranean climate, altered from the altitude which ranges between 450 and 600 m a.s.l. In the city area and between 450 and 1356 m a.s.l. (i.e., at the top of Mt. Menechador) in the whole municipality area. Temperatures range from average 5 • C in winter up to average 26 • C in the summer. Rainfalls are strongly influenced by the proximity to the Gulf of Valencia; the average yearly rainfall amount ranges from 300 mm in the southern arid region up to 500 mm and the maximum rainfalls generally occur in autumn while summer is dry and arid.
The study zone is widely affected by landslides. During the Upper Miocene and Pliocene, this intramontane basin was an endorheic one. During the Quaternary river erosion allowed the connection of this basin with the Mediterranean Sea, a sudden change in the base level of rivers took place, which led to a period of intense incision by erosion of slopes that favors the occurrence of landslides [23]. Currently, shallow landslides mainly occur in relation to intense rainy periods [24,25], although deep and large roto-translational landslides can be also related to different triggers [26]. As it regards earthquake-induced landslides, they are well documented for the Alcoy area [27].
On December 2, 1620, an earthquake (M w 5.5, Imax VII-VIII in EMS-scale) triggered several landslides that affected slopes of the rivers Barchel and Molinar, causing damage (including collapse) to houses located along several streets of the city in that time ( Figure 2). According to historical descriptions, they correspond to moderate sized disrupted landslides (sensu Keefer, [2]) although a large coherent slide was reactivated during this event (instability 1620-III in Figure 2). In historical time, a low magnitude event (M w 4.8) triggered two new, small sized landslides in the urban area of Alcoy [23][24][25][26][27] (Figure 2).
The study area is located in the eastern portion of the Betic Cordillera ridge, south east of the geological Alcoy intramontane basin.
The Betic Cordillera is an Alpine orogen originated by the collision between the African plate and the Eurasian one, which extends 600 km south of the Iberian Peninsula, from the Gulf of Cadiz to the southern province of Valencia and the Balearic Islands [28].
On the basis of the outcropping rocks and their degree of metamorphism, the Betic Cordillera is distinguished in four main geological zones: outer domain; internal domain; complex of Gibraltar and Neogene-quaternary intramontane basins [29]. The Alcoy area is one of these intramontane basins, located within the outer domain of the Betic Cordillera [30].
The geological deposits which outcrop in the area of Alcoy ( Figure 3) include Oligocene bioclastic limestones, interlayered to conglomeratic and marly levels, Miocene marls and clays (locally named "Marne del Tap Formation") widely outcropping in the Alcoy basin, bioclastic and biodetritic limestones that outcrop in the reliefs which surround the basin.
On the basis of the outcropping rocks and their degree of metamorphism, the Betic Cordillera is distinguished in four main geological zones: outer domain; internal domain; complex of Gibraltar and Neogene-quaternary intramontane basins [29]. The Alcoy area is one of these intramontane basins, located within the outer domain of the Betic Cordillera [30].
The geological deposits outcropping in the area of Alcoy ( Figure 3) include Oligocene bioclastic limestones, interlayered to conglomeratic and marly levels, Miocene marls and clays (locally named "Marne del Tap Formation") which widely outcrops in the Alcoy basin, bioclastic and biodetritic limestones that outcrop in the reliefs which surround the basin. At the bottom of the main reliefs alluvial fans of Pleistocene and Holocene ages outcrop, which are composed of reddish gravels and clays. Other quaternary deposits consist of alluvial and travertine terraces and alluvial deposits which fill the central part of the Alcoy basin. Travertine terraces are located above the Miocene marls and clays, with thicknesses varying depending on the suffered erosion while the alluvial terraces are distributed at different levels above the present river channel and are mainly composed of calcareous sediments eroded from the surrounding.

Rock Slopes
To analyze the susceptibility to failure of rock slopes ( Figure 4), a GIS procedure on a cell-by-cell basis was implemented to verify the compatibility of different failure mechanisms (wedge, planar sliding and toppling), based on data collected through geomechanical surveys according to the ISRM standards [31]. A zonation of the outcropping rock masses in terms of orientation and spacing of main joint sets is required for this analysis. The attribution of geo-structural parameters over the whole territory (continuous information) starting from punctual survey sites (discrete information) is a challenging issue since too many variables influence the distribution of such parameters. Furthermore, it has to be considered that the clustered location of survey sites, in its turn related to the availability and accessibility of significant outcrops, could greatly limit the reliability and accuracy of a rigorous geostatistical approach. For these reasons an approach based on expert judgment, i.e., the capacity of a geologist to assess to what extent a punctual datum can be attributed to an area, was adopted.
The so derived synthethic stereplots were analyzed to catch the relationships with the morpho-structural setting; the related synthethic information were then attributed to zones considered homogeneous in terms of morphological, lithological and structural features, thus partitioning the study area in the HGZ that are shown in Figure 5.
Based on geomechanical data assigned to each HGZ, tests to check the compatibility of different failure mechanisms were carried out according to Markland tests [32] through a GIS procedure. More specifically the kinematic compatibility to wedge and planar sliding as well as to toppling was verified by applying a conditional analysis of the discontinuity data with respect to slope angle and slope aspect of each single cell pertaining to a given HGZ. This procedure consists in a cell-by-cell analysis on a 5 m × 5 m square mesh grid, where data referred to dip and dip-direction of discontinuity sets are extrapolated from the corresponding HGZ.
The Markland tests [32] performed only assess the kinematic compatibility of the main discontinuity set (for the planar sliding mechanism) or intersection of sets (for wedge sliding mechanism) or release surface (for toppling mechanism). A further filter is then required to verify if other discontinuity sets acting as lateral and/or rear release surfaces are present: if all these necessary conditions in the HGZ are present, the mapping units that resulted kinematically compatible by the Markland tests are finally considered as potentially unstable.
By considering the average spacing and the orientation of discontinuity sets which isolate potentially unstable rock blocks, representative volumes and shapes can be attributed to each mapping unit to be analyzed in the subsequent step (Slope Stability).
In the Alcoy case study, 6 typical blocks for planar sliding, wedge sliding and toppling have been identified at all, nevertheless, among these blocks the ones having volumes lower than 0.03 m 3 were excluded for slope stability analyses ( Figure 6).

First Time Earth-Slides
First-time earth-slides were accounted for by analyzing proneness to soil sliding involving slope covers. At this aim, an engineering-geological survey was carried out, supported by a GIS dataset. The mapping unit selected for assessing the susceptibility to shallow earth-slides is represented by the Unique Condition Units (UCU), which represents areas that are homogeneous in terms of slope, aspect, lithology and thickness of the outcropping soil covers. At this aim the study area was firstly divided in small basins by performing GIS-based hydrological analyses consisting in the automated extraction of watersheds from the local DEM. This zoning was then overlaid and intersected with other relevant thematic layers and zoning such as slope and aspect classes, polygons representing outcrops of colluvia and debris in their turn differentiated on the basis of the soil thickness. UCU are then the result of such an intersection of geo-thematic data. For the Alcoy case study 35 UCU were identified (Figure 7), having a total extension of 0.18 km 2 ; the 60% of these UCU correspond to a slope dip ranging from 20 • to 40 • , the 14% from 40 • to 50 • , while the remnant percentage is referred to slope dip lower than 20 • . The susceptibility analysis was performed according to the Frequency Ratio (FR) method [33] which measures the reliability of the causative relation between given predisposing factors and landslide processes. A FR value is computed for each class of each considered factor as the ratio of proportion of landslide in a factor class to proportion of area of the same class. Such a relation is significant if the FR value is higher than 1. By summing on a cell-by-cell basis all the Frequency Ratios attributed to each class of all the considered factors, a Landslide Susceptibility Index (LSI) can be defined for each mapping unit.
As it results from this analysis, "lithology" and "aspect" factors do not represent discriminating factors and they were, therefore, neglected for following analysis. More in particular, while the "lithology" factor is too homogeneous the aspect cannot be regarded as a discriminant one since all the derived classes show quite similar FR values.
As a consequence, the "slope" factor resulted the only relevant one, so a LSI threshold of 1 has been chosen to select the potentially unstable UCU.

Already Exsisting Landslides
Already existing landslides were inventoried in the Alcoy municipality area from literature data, aerial photo analysis, and online catalogs as well by field surveys [23,24,34]. The surveyed landslides were mapped (Figure 8), with different failure mechanisms for a total area of 0.92 km 2 .  Table 1. Table 1. Geomechanical parameters. γ n is the total weight for unit volume, c is the cohesion, Φ is the friction angle, n is the normal stress acting on joint, JRC is the Joint Roughness Coefficient and the JCS the Joint Compressive Strength parameter [31]. The HGZ are referred to Figure 5.

PARSIFAL Step 2: Slope Stability (SS)
Slope stability analyses were carried out to derive critical pseudostatic acceleration (a y ), suitable for computing co-seismic displacements following the Newmark's method [22]. At this aim, sensitivity analyses were performed by varying both saturation conditions and pseudostatic actions. Different failure models were assumed for rock slides and toppling, first-time earth-slide and already existing landslides.
The slope stability of the identified rock blocks was carried out by a pseudostatic approach performed according to conventional computational methods based on Global Limit Equilibrium (GLE) approach [35]. The yielding acceleration coefficients, a y , were derived by a sensitivity analysis performed for each rock block and by assuming joint saturation conditions (expressed by the ratio between wet and total height, h w /h) that range from dry to completely wet. First-time earth-slides were analyzed according to a shallow sliding model represented by infinite slope conditions and assuming a pseudostatic seismic action for different saturation conditions of the earth mass to compute a y .
For the already existing landslides, engineering-geological cross-sections were obtained to derive geometrical features of the landslide mass ( Figure 9) and attributing geotechnical properties to perform slope stability analyses. Figure 9. Example of simplified engineering-geology cross-section obtained for the Molinar landslide (see Figure 8 for location). The slices used for the GLE slope stability are also shown.
The inventoried landslides were classified by considering morphological and geometrical features, according to Dewitte and Demoulin [36]: Janbu method [37] was applied for slides characterized by 00.1 ≤ H/L ≤ 0.1, i.e., with a mainly translational mechanism; Bishop method [38] was applied for a roto-translational mechanism (where H/L > 0.1). Also in this case, a pseudostatic seismic action was assumed at different saturation conditions expressed by the Bishop coefficient for water pressure distribution (r u ) to compute the a y coefficient.
The degree of mobility is expressed by the exceedance probability of a critical displacement (P[D ≥ Dc|a(t), a y ]) responsible of the slope collapse, herein assumed to be 10 cm for earth slides and 5 cm for rockslides [4,20]. The specific values considered for the slope failure threshold could be of slightly different [12] and this could obviously condition the results, but only in terms of computed (P[D ≥ Dc|a(t), a y ]) values, i.e., not changing the severity ranking of the expected slope failures for each considered scenario. At this regard, it is worth stressing that the PARSIFAL approach mainly aims at restituting a synthetic mapping, suitable to define a priority ranking for planning more detailed studies and/or engineering designs; therefore, the computed P[D ≥ Dc|a(t), a y ] values assume a secondary role. The distribution of co-seismic displacements from which probabilities are computed, derives from the application of a set of real accelerograms matching, for each site or zone, the proper seismic hazard scenario according to the selection criteria described in the following sections.

Probabilistic Seismic Hazard
The seismic inputs to be used in the Newmark's analysis were selected after: (1) defining the local seismic hazard through the identification of seismogenic sources responsible of most severe earthquakes that struck the study area and the relative seismological parameters; (2) evaluating the possible level of local ground shaking due to each seismogenetic source by the calculation of response spectra attenuated to the site by the Sabetta and Pugliese [39] attenuation law, their ranking and selection of the most severe as target spectrum (TS); (3) selecting Time Histories (TH) from accelerometric databases queried both in term of seismological parameters and with a spectrum-similarity respect to the TS; scaling the selected TH to the probabilistic PGA calculated for the Alcoy area for a Return Time (RT) of 475 years [34].
The selected seismogenic sources ( Figure 10  All the seismogenic sources in the neighborhood of Alcoy are located within a distance of about 10 km. Ranking the response spectra attenuated at the site by the Sabetta and Pugliese [39] attenuation law, although the closest seismogenic source is Alcoy-Cocentaina (E623), the most sever effects in Alcoy are expected to be due to by the seismogenic source of Benasau (E624).
The expected ground motion at each landslide site was obtained from natural time-histories (TH) stored in the strong motion database ITACA (Italian ACcelerometric Archive-version 2.0).
At this aim, TH were chosen through the following criteria. regarding the seismogenic source properties, the focal mechanism was considered "normal" or "normal-oblique", the M w interval ± 0.3 of the maximum earthquake inferred from the seismic hazard analysis (i.e., 6.4 to 7.0), the hypocentral depth ≤ 30 km. As regards the recording-site conditions, free-field" or "ground" installation were considered as well as the Eurocode 8 ground type "A" (hard rock/rock or Vs_30 ≥ 800 m/s) or "B" (soft rock/very dense soil or 360 ≤ Vs_30 < 800 m/s) and flat-slope conditions (i.e., slope angle < 15 • ).
As it regards the ground-motion parameters, the intensity of Arias of the time-history (AI TH ), was select in the range: average log 10 (AI) − 2σ log 10 ≤ (AI TH ) ≤ average log 10 (AI) + 2σ (1) Furthermore, a spectrum-similarity analysis was obtained comparing the shape of the TH to the shape of response spectra of each TH so far selected ( Figure 11). Finally, the resulting TH was scaled to the PGA calculated by the Istituto Geografico National, Seismic Hazard of Spain [34], for a return period of 475 years (i.e., 10% probability of exceedance in 50 years). Figure 11. Response spectra of the selected input (grey) in comparison with the average response spectrum (blue) and target response spectrum (red).

Results and Mapping
The computed P[D ≥ Dc|a(t), a y ] for each category of landslide (i.e., rock slides, topples, first-time earth-slides, existing landslides) were mapped for the considered saturation scenarios.
To restitute a comprehensive mapping, a GIS solution was adopted providing a raster zonation for rocky landslides, by contouring the grid cell susceptible for different failure mechanisms (i.e., planar sliding, wedge sliding, toppling). A weighted probability was initially computed for each failure mechanism respect to the estimated volume of typical blocks and, subsequently, a probability combined was obtained considering different mechanisms. It was possible to map a unique probability value for grid cell representative for the most critical slope condition. Critical threshold for computing P[D ≥ Dc|a(t), a y ] was assumed at 10 cm for earth-slides and 5 cm for rock-slides, according to Romeo [20].
To provide a unique value of P[D ≥ Dc|a(t), a y ] into the same cell unit, the following criteria were assumed:

•
For the same mechanism, the P[D ≥ Dc|a(t), a y ] computed for different block volumes are solved in a unique value by performing a "weight averaging"; • "combined probability" is computed for each different mechanism (i.e., toppling, wedge sliding, planar sliding) using the probability values computed for each block volume; • In case that only SF values are computed for each block volume and landslide mechanism the "lowest SF value" is attributed to the cell unit.
In a first scenario dry conditions and seismic action corresponding to a return time of 475 years (i.e., corresponding to a PGA of 0.16 g) were considered.
The results obtained for this scenario are reported in Figure 12 and Tables 2-4. In a second scenario, middle saturation conditions and the same seismic action of the dry scenario were considered; the obtained results are reported in Tables 5-7.     In the third scenario, full saturation conditions and the same seismic action considered for the previous scenarios of dry conditions and middle saturation were assumed. The obtained results are reported in Tables 8-10 and in Figure 13. Table 8. Results obtained in terms of exceedance probability for the full saturation scenario (RT = 475 years).  Table 9. Safety factor values computed in the full saturations scenario for different landslide typologies.

Discussion
The analysis of the resulting scenarios by PARSIFAL at Alcoy highlights that the area involved in landslide in the Alcoy municipality increases with increasing saturation. In the most severe scenario (i.e., total saturation) unstable conditions are reached by among the 85% of the areas susceptible to failure.
These results demonstrate the main role of rainfalls which represent the predominant and most hazardous triggering action in the basin of Alcoy since they occur seasonally. The seasonal rainfalls are generally concentrated in Autumn, after dry and arid Summer seasons.
All the earth-slides involve the Marne del Tap, marly formation that suffer shrinkage during the arid summer generating tension cracks; during the Autumn rainfalls these cracks fill of water so favoring infiltration, soil saturation and earth-slides. The computed P[D ≥ D c |a(t), a y ] are consistent with specific situations that can be monitored over time since landslide activity and stabilization works are well documented; this is the case of the Molinar and La Beniata landslides within the Alcoy urban area. The results obtained for rock landslides indicate that very hazardous conditions interest an extended zone of the Alcoy urban area located NW (namely Sincra Carrer Barranc), close to the mountain relief where the carbonatic succession widely outcrop, so producing high risk scenarios.
The PARSIFAL applied to the Alcoy's urban area demonstrates its suitability for managing landslide hazard over extended areas leading to higher awareness in environmental planning. In addition to the safety factor, commonly used for rating landslide masses which are prone to failure [40] as well as to the co-seismic displacement computed according to conventional approaches [41], the probability of exceedance used by PARSIFAL allows a thematic mapping which is more suitable in the framework of seismic microzonation studies. Since the earthquake scenarios are linked to the local seismic hazard, the computed probability of exceedance multiplied for the seismic hazard can provide a combined and unconditional probability of slope failure occurrence. In Italy the PARSIFAL is going to be proposed for deriving the mapping of unstable areas due to landslide following the guidelines proposed by the National Civil Protection [21,42]. More in particular, "attention zones" for landslide proneness can be mapped through PARSIFAL and "areas of respect" can be selected among them after further and more detailed investigations.

Conclusions
The PARSIFAL was tested in the area of Alcoy municipality, demonstrating its suitability in restituting scenarios of earthquake-induced landslides for a specific seismic hazard and at different saturation conditions. A comprehensive mapping, in terms of exceedance probability of a critical displacement allowed to output high risk areas located close to the carbonatic relief located NW of the Alcoy urban areas, due to the susceptibility to rock slides and toppling as well as along the Molinar river valley due to the widely outcropping marls of the Marne del Tap formation responsible for diffused earth-slides. As it regards the Marne del Tap formation, a relevant effect due to seasonal rainfalls, which generally occur in Autumn after the arid Summer season, was demonstrated.
Future applications of PARSIFAL could take into account the intrinsic variability of geomechanical parameters attributed to rock joints or soil to compute safety factors as well as exceedance probability of co-seismic displacements by performing stochastic parametric analyses. An automatic procedure could be improved to speed up computation and processing for the comprehensive thematic mapping. Moreover, since roads re generally not evident in DEM, because of too low resolution, further applications could be focused on taking into account the role of man-made cuts along roads in earthquake-induced landslide scenarios.