Hazard-Consistent Earthquake Scenario Selection for Seismic Slope Stability Assessment

: Design ground shaking intensity, based on probabilistic seismic hazard analysis (PSHA) maps, is most commonly used as a triggering condition to analyze slope stability under seismic loading. Uncertainties that are associated with expected ground motion levels are often ignored. This study considers an improved, fully probabilistic approach for earthquake scenario selection. The given method suggests the determination of the occurrence probability of various ground motion levels and the probability of landsliding for these ground motion parameters, giving the total probability of slope failure under seismic loading in a certain time interval. The occurrence hazard deaggregation technique is proposed for the selection of the ground shaking level, as well as the magnitude and source-to-site distance of a design earthquake, as these factors most probably trigger slope failure within the time interval of interest. An example application of the approach is provided for a slope near the highway in the south of Sakhalin Island (Russia). The total probability of earthquake-induced slope failure in the next 50 years was computed to be in the order of 16%. The scenario peak ground acceleration value estimated from the disaggregated earthquake-induced landslide hazard is 0.15 g , while the 475-year seismic hazard curve predicts 0.3 g . The case study highlights the signiﬁcant di ﬀ erence between ground shaking scenario levels in terms of the 475-year seismic hazard map and the considered fully probabilistic approach.


Introduction
Several approaches are used to evaluate the seismic stability of a slope: deterministic, statistical, and, last but not least, probabilistic.
Physically-based models, as a type of deterministic approach, are applied for site-specific earthquake-induced landslide hazard assessment. The numerical simulation of internal stresses and strains gives the most realistic view of the behavior of slopes during earthquakes [1]. Stress-deformation analysis requires well-identified and measured soil properties; at the same time, ground motion selection remains a subjective matter.
Most of the statistical and probabilistic approaches are designed to produce hazard maps on regional or global scales. Because material parameters are difficult to identify in detail for large areas, slope parameters in many cases are estimated using Global Information System (GIS) analysis tools [2]. These maps help evaluate the general hazard level and specify sites, in which detailed geotechnical investigations are needed. Regional, seismically-induced landslide hazard maps can be produced in terms of susceptibility or probability [3]. Susceptibility maps [4] highlight the areas susceptible to landsliding, based on the physical Figure 1. The seismicity map and seismic source model of the studied area. The studied area is marked by the filled red circle. The association of the area and line sources with their corresponding seismicity parameters is given in Table 2.

Fully Probabilistic Approach
The fully probabilistic assessment of a seismically-induced landslide hazard requires two crucial stages: (1) Evaluation of the probability of occurrence of a certain value of ground motion for the studied area within the time interval of interest; (2) evaluation of the conditional probability that given ground motion parameters trigger a landslide. Several geometrical and mechanical parameters of the slope should be taken into account in the second stage. We use peak ground acceleration (PGA) as a ground shaking intensity.
After the probabilities defined in stages (1) and (2) are known, the total probability of slope failure in the next years according to the total probability law is calculated as slope failure PGA • slope failure| , model , where PGA is the probability of occurrence of PGA in the next years, slope failure| , model is the conditional probability of slope failure under seismic loading within the given slope model and • slope failure| , model is the probability of slope failure in the next years by the scenario of ground motion intensity = within the given slope model .

Fully Probabilistic Approach
The fully probabilistic assessment of a seismically-induced landslide hazard requires two crucial stages: (1) Evaluation of the probability of occurrence of a certain value of ground motion for the studied area within the time interval of interest; (2) evaluation of the conditional probability that given ground motion parameters trigger a landslide. Several geometrical and mechanical parameters of the slope should be taken into account in the second stage. We use peak ground acceleration (PGA) as a ground shaking intensity.
After the probabilities defined in stages (1) and (2) are known, the total probability P T of slope failure in the next T years according to the total probability law is calculated as where P T (PGA = a i ) is the probability of occurrence of PGA = a i in the next T years, P slope failure a i , model j is the conditional probability of slope failure under seismic loading a i within the given slope model j and p ij = P T (a i )·P slope failure a i , model j is the probability of slope failure in the next T years by the scenario of ground motion intensity = a i within the given slope model j. Equation (1) should be computed for all possible discrete values of PGA that may trigger landslides and for all available geomechanical slope models ranked by weights w j , where j w j = 1. ( Del Gaudio et al. [18] introduced a method to calculate the full occurrence probability of a slope failure at a given site during a time interval of interest. The authors called this the time-probabilistic approach, which is very similar to Equation (1).
Olsen et al. [19] used the modified form of Equation (1) for regional, seismically-induced landslide hazard analysis and mapping in Oregon (USA). A novelty of the given method is the use of the deaggregation of the total probability of slope failure as a target in the selection of design source parameters (source-to-site distance, magnitude) and ground shaking level (see Section 2.1.4). The same approach is applied in structural analysis, aiming at the determination of the response of facilities to seismic loading. The parameters of the design earthquake may further be used by geotechnical engineers to specify the seismic loading in stress-deformation analysis.

Conditional Probability of Slope Failure Under Seismic Loading
The dynamic model of slope stability, based on the Newmark approach [24], considers landslide mass as a solid block sliding along the planar surface. The static factor of safety can define the capacity of the slope to resist the ground shaking. In our work, we use a rather simplified limit-equilibrium model of an infinite slope [25].
The assumptions of the model are as follows: (1) The sliding mass is assumed to be a rigid solid; (2) the static factor of safety is stress-independent, and therefore, is constant; (3) the coefficients of static and dynamic friction are constant and equal; (4) the effects of pore pressure of geological medium are negligible; and (5) the soil characteristics are homogenous.
According to the limit equilibrium theory, the static factor of safety F S , defined as the relationship between the force keeping the sliding mass on the slope to the force moving the sliding mass down the slope, might be defined as [9] where c is the effective cohesion, z is the slope-normal thickness of the potential sliding mass, α is the dip angle of the potential sliding surface, γ is the material unit weight, γ w is the unit weight of groundwater, ϕ is the effective friction angle and m is the proportion of the sliding mass thickness that is saturated. The sliding mass is stable if F S > 1, and unstable when F S < 1. According to Newmark's approach, a landslide will not start until the sliding mass accumulates some internal deformations. The accumulation of internal deformations takes place when the seismic acceleration exceeds some critical value.
Newmark [24] has shown that critical acceleration of slope failure a c is the simple function of the static factor of safety and the slope geometry: where g is the factor of gravity. The time interval of the accumulation of internal deformations depends on the level and duration of ground shaking. The cumulative Newmark displacement D N was proposed by Newmark [24], and later applied by Jibson [9,26] to produce the seismic landslide hazard index.
Several regression curves between the Newmark displacement and ground motion parameters (PGA, Arias intensity and others) are known [11,25,26]. To assess the slope stability under seismic loading, measured as the peak ground acceleration, the following equation was used [27]: where D N is in cm, the last term is the standard deviation and a is the peak ground acceleration. The Newmark displacement D N does not necessarily correspond to the coseismic deformations measured directly. The predicted sliding displacement in (5) describes the probability of slope failure given the triggering conditions. The accurate calibration of such probabilities would require detailed information regarding the regional features of the seismically-induced landslide processes and their control factors. The recorded data were collected after the 1994 Northridge earthquake (CA, USA) [5]. Jibson et al. [9] have compared the spatial distribution of seismically-induced landslides with the predicted Newmark displacement and developed a probabilistic model of slope failure under seismic loading, which corresponds well to the Weibull distribution. Unfortunately, the authors have rare data for calibrating the probability model for Sakhalin Island; thus, Jibson's probabilistic model for California [9] was imported into Equation (1): The displacement D N in (6) is defined by parameters of the slope model and ground motion intensity a.

Probability of Occurrence of Ground Shaking Intensity
The main goal of probabilistic seismic hazard analysis (PSHA) is the determination of the probability of exceeding a certain PGA level within the time interval of interest [28]. The result of such analysis is a seismic hazard curve: The relationship between the mean annual frequency λ of exceeding acceleration a: where n sources is the number of seismic sources, n M is the number of possible magnitudes, n R is the number of possible distances, ϑ(M i > m min ) is the annual frequency of earthquakes with magnitude greater than m min in the source i, P PGA > a m j , r k is the probability of exceeding a certain acceleration a for an earthquake with magnitude m j at the source-to-site distance r k , P M i = m j and P R i = r j is the probabilities of occurrence of discrete magnitudes and distances in source i, respectively. The probability of exceeding the given PGA value for an earthquake with magnitude m at a source-to-site distance r corresponds to the standard normal cumulative distribution function: where g(m, r) predicts the mean natural logarithm of PGA for a given m and r and σ is the standard deviation of model g(m, r). The g(m, r) function is commonly known as the ground motion prediction equation (GMPE). The cumulative distribution function for magnitude is given by According to (9), the probability of the occurrence of discrete magnitudes is defined as The truncated Gutenberg-Richter recurrence law was used to quantify the frequency-magnitude distribution: where m max is the maximum magnitude and b is the slope of the Gutenberg-Richter distribution. Equation (7) integrates our knowledge about the rates of occurrence of earthquakes, the possible magnitudes and distances and the models that predict the ground motion parameters of those earthquakes.
Under the additional assumption that all seismic sources follow the Poissonian occurrence process, the probability of exceeding the intensity a during the next T years is given by For small probabilities, the value of λ is small compared to unity, and therefore, the probability in Equation (12) is approximately λT. In other words, the annual probability is approximately equal to the mean annual frequency.
The transition from the probability of exceeding the intensity to the probability of occurrence of the given intensity is as follows: Thus, the total probability of seismically-induced landslide occurrence in the next T years is computed by substituting Equations (6) and (13) into (1).

Occurrence Hazard Deaggregation
The deaggregation of seismically-induced landslide hazard, as well as seismic hazard deaggregation, serves several purposes: (1) The estimation of the contribution to the total hazard level from different magnitude-distance bins, and (2) the selection of the earthquake scenario (design earthquake) which is used for the ground motion selection.
In the case of PSHA, the deaggregation procedures require that the probability of exceeding the intensity should be expressed as a function of magnitude and/or source-to-site distance. In other words, the probability of exceeding the intensity is computed on a 2D magnitude-distance grid with increments (∆m, ∆r). This means that for each fixed (m, r) bin, Equation (7) is applied [29].
The seismic hazard deaggregation is produced in terms of exceeding the intensity or occurrence of a given intensity. The exceeding deaggregation approach is most commonly used in seismology.
At the same time, the use of the exceeding deaggregation approach as a target in ground motion selection is not consistent with the objective of the structural analysis, which is aimed at the determination Sustainability 2020, 12, 4977 7 of 14 of the seismic response of a structure to a seismic loading [30]. Fox et al. [30] have found that, for this form of assessment, one should use the occurrence deaggregation approach. The seismic slope stability assessment is closely related to the structural analysis as these approaches deal with the seismic response.
In the case of the occurrence deaggregation approach, the band of peak ground accelerations a − ∆a 2 , a + ∆a 2 should be considered. The results of hazard deaggregation, in that case, are expressed as It has been shown [30] that the width of the seismic intensity band ∆a is not very sensitive to the results of the deaggregation analysis. Thus, Equation (14) was used in the current study for the computation of the earthquake scenario.

Studied Area
For the case study, a local slope susceptible to landsliding was selected. The target area is located in the south of Sakhalin Island near the highway from Yuzhno-Sakhalinsk to Kholmsk (Figure 2). Shallow landslides were triggered here by human activity around the solid waste landfill (Figure 2). The block slide of the slope was covered by waste (Figure 3a), and snow melting in April 2018 caused the inundation of the slope bottom and refuse bank, forming a sliding surface and finally triggering a landslide. Furthermore, the deposition of a heavy load of concrete blocks on the top of the slope contributed to landsliding on April 15-16 2018. The landslide activity caused tension cracks along the highway (Figure 3b). For the case study, a local slope susceptible to landsliding was selected. The target area is located in the south of Sakhalin Island near the highway from Yuzhno-Sakhalinsk to Kholmsk (Figure 2). Shallow landslides were triggered here by human activity around the solid waste landfill (Figure 2). The block slide of the slope was covered by waste (Figure 3a), and snow melting in April 2018 caused the inundation of the slope bottom and refuse bank, forming a sliding surface and finally triggering a landslide. Furthermore, the deposition of a heavy load of concrete blocks on the top of the slope contributed to landsliding on April 15-16 2018. The landslide activity caused tension cracks along the highway (Figure 3b).

Geomechanical Slope Model
In terms of geology, the core part of the studied area is formed by Paleogene deposits (Figure 4). Soils of the Quaternary period from the core part of the studied area and include both human-made and Deluvial/Colluvial deposits [31]. The soils in a natural state have poor filtration characteristics, a high soaking and swelling capacity, a low internal friction angle and variable degrees of heaving.

Geomechanical Slope Model
In terms of geology, the core part of the studied area is formed by Paleogene deposits (Figure 4). Soils of the Quaternary period from the core part of the studied area and include both human-made and Deluvial/Colluvial deposits [31]. The soils in a natural state have poor filtration characteristics, a high soaking and swelling capacity, a low internal friction angle and variable degrees of heaving. For the hazard assessment, the geomechanical slope model was built for the local area where the landslide was triggered in April 2018. Their characteristics and parameters are presented in Table 1. Material parameters are selected from the geological and geophysical reports of the studied area [32], and should be considered as a simplified model. For the given model (Table 1), the critical acceleration in accordance with Equations (3,4) is 0.023g.   For the hazard assessment, the geomechanical slope model was built for the local area where the landslide was triggered in April 2018. Their characteristics and parameters are presented in Table 1. Material parameters are selected from the geological and geophysical reports of the studied area [32], and should be considered as a simplified model. For the given model (Table 1), the critical acceleration in accordance with Equations (3) and (4) is 0.023g.

Seismic Source Model and GMPEs
Area and line sources were aggregated into one seismic source model (Figure 1). The area sources around the studied area were selected from the regional area source model for Sakhalin Island [33]. Basically, the regional area sources were defined by analyzing the distribution of the shallow seismicity and known faults. Area sources describe the background seismicity pattern of small to moderate events up to magnitudes of Mw = 6. The area sources are considered in this study as volume sources with a uniform depth distribution from 5 to 15 km.
Line sources were selected from the source model of the general seismic zonation map of the Russian Federation [34]. The dynamic potentials of the line sources in our model are Mw = 6.5, 7 and 7.5 ( Table 2). Seismicity parameters for each of the considered sources are summarized in Table 2. The choice of the GMPE is very important for the hazard assessment because of its great influence on the final results. Regional PGA attenuation curves were used in this study to predict the ground motion parameters for the area sources [35]: log PGA = 0.87 * Mw − log R rup + 0.006 * 10 0.5 * Mw − 0.0038 * R rup − 1.726 ± 0.336, (15) where R rup is the rupture distance (km), Mw is the moment magnitude, and the last term is the standard deviation of the model. The peak ground acceleration in (15) is given in cm/s 2 . The average V S30 in the area was measured~300 m/s [35].
For the linear sources, the Next Generation Attenuation (NGA) West-2 GMPE model [36] was applied per the general results in Reference [35]. The NGA-West-2 GMPE model was corrected for the local soil conditions.

Results
The seismic hazard curve was computed using CRISIS 2015 software [37]. All calculations are performed for accelerations exceeding the critical value of 0.023g ( Figure 5).
The design ground shaking intensity for the 475-year return period (10% exceeding probability) corresponds to a PGA of 0.3g (Figure 5a). In terms of macroseismic intensity, it is of the order of VIII-IX MSK64 and in good agreement with the more superficial estimates based on the maps of general seismic zonation [34].
The distribution of the occurrence probability has a clear peak corresponding to the PGA value in the order of 0.07g (VI-VII MSK64 intensity) (Figure 5b). This is the most probable strong ground motion scenario for the studied area in the next 50 years, according to the given seismic source model and ground motion prediction models.
The probability of slope failure in relation to the ground shaking level is shown in Figure 5c. The curve is plotted according to Equation (6). The same figure shows the cumulative curve, according to Equation (1). The probability of slope failure in the next 50 years is obtained by summing the landslide occurrence probability p ij for the wide range of ground motion scenarios. The total probability of seismically-induced landslide occurrence in the next 50 years appeared to be about 16%. (Figure 5c). Thus, a high probability shows a considerable hazard of seismically-induced landslides in terms of the application of civil engineering. In order to estimate the most probable ground motion scenario of seismically-induced landslide, the discrete probabilities from Equation (1) were plotted versus the corresponding ground motion level, contributing to the total probability of the seismically-induced landslide hazard from each ground motion level (Figure 5d). The distribution has a modal form, with the peak corresponding to the PGA value of 0.15g (VII-VIII MSK64 intensity). When estimating the ground motion level that most probably causes landslides within a given time window, it is possible to perform an occurrence deaggregation of seismic hazard considering the bandwidth of ground shaking intensity. This helps to identify the parameters of the design In order to estimate the most probable ground motion scenario of seismically-induced landslide, the discrete probabilities p ij from Equation (1) were plotted versus the corresponding ground motion level, contributing to the total probability of the seismically-induced landslide hazard from each ground motion level (Figure 5d). The distribution has a modal form, with the peak corresponding to the PGA value of 0.15g (VII-VIII MSK64 intensity). When estimating the ground motion level that most probably causes landslides within a given time window, it is possible to perform an occurrence deaggregation of seismic hazard considering the bandwidth of ground shaking intensity. This helps to identify the parameters of the design earthquake (magnitude and source-to-site distance). Knowing the ground shaking level, magnitude, and source-to-site distance makes the earthquake scenario selection possible. The occurrence hazard deaggregation, according to (14) is found to be within the acceleration range of 0.1-0.2g, which corresponds to the median PGA level of 0.15g with a variance of ±0.05g. The probability distribution in relation to the magnitude-distance bins is shown in Figure 6. The peak of the distribution corresponds to earthquake design parameters of Mw = 5.5 and R rup = 5 km. Therefore, for the given seismic and slope models, the most probable scenario of an earthquake that would trigger a landslide in the next 50 years is defined by the parameters PGA = 0.15g, Mw = 5.5 and Rrup = 5 km.

Discussion
In the framework of the given models, a seismic event that triggers a slope failure in the next 50 years is most likely to have the parameters PGA = 0.15g, Mw = 5.5 and Rrup = 5 km. Regarding the earthquake design parameters, the cumulative Newmark displacement, according to (5) is about 16 cm. The 475-year shaking map intensity predicts a PGA of 0.3g (Figure 5a) and the corresponding Newmark displacement is about 53 cm, which is almost three times higher than the design sliding displacement predicted by the fully probabilistic approach.
One more modal distribution in the hazard deaggregation plot could be recognized in Figure 6. The peak has a smaller amplitude and corresponds to a seismic event with Mw = 6.5 and Rrup = 40 km. According to the seismic source model (Figure 1), this event is most likely to correspond to the thrust fault type event generated by a line source. The authors suggest that high-magnitude events are likely to have a long duration. The duration of ground shaking has an influence on the slope stability; thus, the earthquake design parameters of PGA = 0.15g, Mw = 6.5 and Rrup = 40 km should be considered Therefore, for the given seismic and slope models, the most probable scenario of an earthquake that would trigger a landslide in the next 50 years is defined by the parameters PGA = 0.15g, Mw = 5.5 and R rup = 5 km.

Discussion
In the framework of the given models, a seismic event that triggers a slope failure in the next 50 years is most likely to have the parameters PGA = 0.15g, Mw = 5.5 and R rup = 5 km. Regarding the earthquake design parameters, the cumulative Newmark displacement, according to (5) is about 16 cm. The 475-year shaking map intensity predicts a PGA of 0.3g (Figure 5a) and the corresponding Newmark displacement is about 53 cm, which is almost three times higher than the design sliding displacement predicted by the fully probabilistic approach.
One more modal distribution in the hazard deaggregation plot could be recognized in Figure 6. The peak has a smaller amplitude and corresponds to a seismic event with Mw = 6.5 and R rup = 40 km. According to the seismic source model (Figure 1), this event is most likely to correspond to the thrust fault type event generated by a line source. The authors suggest that high-magnitude events are likely to have a long duration. The duration of ground shaking has an influence on the slope stability; thus, the earthquake design parameters of PGA = 0.15g, Mw = 6.5 and R rup = 40 km should be considered as an additional earthquake scenario that may lead to slope instability in the next 50 years.
Constants in Jibson's landslide probability model [9] could differ in other regions if the strengths of the geologic medium, topography or soil moisture conditions are significantly different from those in southern California. In this case, the shape of the curve (Figure 5c) would be different from the predicted shape (6). For the probability hazard analysis, variable constants from Equation (6) should be proposed.
The current practice of seismic slope stability analysis does not take into account the spatial and seasonal variability of soils. The weakened sites of the soils may be missed during the geotechnical investigation, which may result in the risk of slope failure for those sites. In other words, the more accurate assessment of seismic slope stability requires the uncertainty analysis of geomechanical slope models. Such analysis may be performed with the simulation of the variability of the material parameters [20,38], which helps to obtain the rational hazard value. In future studies, we will analyze the uncertainties associated with the spatial and seasonal variability of soils and their influence on the risk of slope failure in the framework of the fully probabilistic approach.

Conclusions
An improved fully probabilistic technique for seismic slope stability assessment is suggested in this study. The method requires four crucial stages for data selection and processing: (1) Building the geomechanical slope models and its uncertainties; (2) building the seismic source models and its uncertainties; (3) building the ground motion prediction models and its uncertainties; and (4) calibrating the landslide probability model. Finally, the occurrence deaggregation technique is proposed for the selection of the ground shaking level, magnitude and source-to-site distance of the design earthquake that would most probably trigger a slope failure within the time interval of interest.
The method was applied considering a local slope in the south of Sakhalin Island, which was known as a seismically active and landsliding region. The total probability of slope failure under seismic loading in the next 50 years was computed to be in the order of 16%. The shaking intensity level estimated from the disaggregated earthquake-induced landslide hazard appeared to be 0.15g, while the 475-year seismic hazard curve predicts a value of 0.3g. The significant difference between ground shaking levels in terms of the 475-year seismic hazard curve and considered fully probabilistic approach suggests that the seismic landslide hazard could be underestimated or overestimated when using 475-year seismic hazard maps as a target for the selection of the scenario triggering condition. The approach given in this article follows the rational risk management idea that handles all possible scenarios.