Shallow Landslide Susceptibility Mapping in Sochi Ski-Jump Area Using GIS and Numerical Modelling

The mountainous region of Greater Sochi, including the Olympic ski-jump complex area, located in the northern Caucasus, is always subjected to landslides. The weathered mudstone of low strength and potential high-intensity earthquakes are considered as the crucial factors causing slope instability in the ski-jump complex area. This study aims to conduct a seismic slope instability map of the area. A slope map was derived from a digital elevation model (DEM) and calculated using ArcGIS. The numerical modelling of slope stability with various slope angles was conducted using Geostudio. The Spencer method was applied to calculate the slope safety factors (Fs). The pseudostatic analysis was used to compute Fs considering seismic effect. A good correlation between Fs and slope angle was found. Combining these data, sets slope instability maps were achieved. Newmark displacement maps were also drawn according to empirical regression equations. The result shows that the static safety factor map corresponds to the existing slope instability locations in a shallow landslide inventory map. The seismic safety factor maps and Newmark displacement maps may be applied to predict potential landslides of the study area in the case of earthquake occurrence.


Introduction
The Caucasus region is among the landslide research hotspots with medium to high landslide susceptibility classification [1].The mountainous area of Greater Sochi, where the Olympic ski-jump facilities have been built, is located in the northern Caucasus.Although the Olympic Winter Games were over, the construction of the ski-jump area related to tourism is still being carried out.During the initial construction of the area, the main loss Quaternary deposits like sand and clay have been cleared out.However, the weathered surface of the Jurassic mudstone, almost underlying the whole area, is widely distributed [2,3].Compared with other sorts of bedrock, the mudstone always features low strength [4,5].The low strength weathered mudstone can easily cause shallow landslides.Another important natural factor contributing to the formation of landslides is the high seismic activity of the region [6].Therefore, the analysis of the slope stability of the study area in either static or seismic states is of significance.
Quantitative methods based on GIS for evaluating landslide susceptibility include statistically-based models and physically-based models [7,8].In previous work, there have been several statistically-based models for landslide susceptibility evaluation, such as the certainty factor (CF) model, artificial neural network (ANN) model, and logistic regression (LR) model, which may be used to assess landslide susceptibility over large areas [9][10][11][12].Physically-based models based on the limit equilibrium method are often used for small area landslide susceptibility, such as the 3D model and infinite slope stability model [7,13].Compared with statistically-based models, physically-based models can obtain the slope safety factors, which can better reflect the mechanism of slope instability.As for the Sochi region, Zakharov et al. [14] and Kuzin et al. [15] have used a geographic information system (GIS) for landslide susceptibility mapping.However, their mapping processes did not concern the local slope stability setting.In numerical calculation of the slope stability, the safety factor is a very crucial index for showing whether a slope is stable or unstable.In a previous study by other scholars, concerning the slope stability analysis in GIS, the safety factor (Fs) values were usually achieved based on empirical equations [16][17][18].In this work, the safety factor (Fs) values for several representative slope angles were calculated by numerical modelling, and a good correlation (R 2 > 0.9) between Fs and slope angle was found.The safety factor (Fs) values for all slope angles were estimated by regression equations for the study area.Regarding slope stability analysis, numerical methods include the limit equilibrium method [19], finite element method [20,21], finite difference method [22], discrete element method [23,24] and so on.Some scholars have used numerical modelling to analyze the stability of some slopes in this region using different methods [2,[25][26][27].In this study, the Spencer method [28], which is a limit equilibrium method, was applied to calculate the slope safety factors (Fs) with the help of commercial program GeoStudio.The region includes several active faults, with the predicted maximum magnitude of earthquakes M max = 7.3 [29,30].Due to high intensity of the territory seismicity, the seismic effect on slope stability should be taken into account.Pseudostatic analysis was used to compute safety factors (Fs) considering seismic effect.The regression curves of Fs and slope angles were built for this specific area.Intersecting these data sets in ArcGIS landslide safety-factor maps were presented.The Newmark method is also widely used to evaluate the seismic landslide susceptibility [31][32][33][34].At the end of the study landslide Newmark-displacement maps were drawn according to the empirical equations proposed by Jibson et al. [16].The Jibson's empirical equations include Arias intensity, which is different from MSK-64 intensity, adopted in Russia.The transformed equations for Newmark-displacement related to MSK-64 intensity were first introduced in this work.Due to lack of literature considering seismic effect on slope stability of the study area, it is necessary to conduct landslide susceptibility maps to guide the future road construction related to tourism.

General Geographical and Geological Background
The Olympic ski-jump complex in Sochi, Russia was selected as the study area.The area is located in the northern Caucasus on the east Black Sea coast (Figure 1a).The area is on the left bank of the Mzymta River and it covers an area of 0.53 km 2 (Figure 1b).
In recent years, some construction related to tourism has been planned.Some work has been carried out, so that the topography relief and the geological deposits have been changed.Some loss Quaternary deposits like sand and clay have been cleared out.Regarding geological conditions (Figure 2), most of the area is covered by Jurassic mudstone as outcropping bedrock.The upper layer of mudstone bedrock is widespread weathered mudstone.The intact mudstone of bedrock is of low permeability (0.001 m/d), while the weathered layer is more permeable (1 m/d).The thickness of the weathered mudstone all over the area is about 2-3 m.Due to the low strength of weathered layer, the area is supposed to be prone to landslide.Through the field investigation, some shallow landslide warning signs, such as new cracks and unusual bulges, could be found as shown in Figure 3.In this study, an earthquake is believed to be the primary trigger factor of landslide in Sochi ski-jump area, rather than rainfall.Therefore, rainfall-induced landslide has not been considered yet.In recent years, some construction related to tourism has been planned.Some work has been carried out, so that the topography relief and the geological deposits have been changed.Some loss Quaternary deposits like sand and clay have been cleared out.Regarding geological conditions (Figure 2), most of the area is covered by Jurassic mudstone as outcropping bedrock.The upper layer of mudstone bedrock is widespread weathered mudstone.The intact mudstone of bedrock is of low permeability (0.001 m/d), while the weathered layer is more permeable (1 m/d).The thickness of the weathered mudstone all over the area is about 2-3 m.Due to the low strength of weathered layer, the area is supposed to be prone to landslide.Through the field investigation, some shallow landslide warning signs, such as new cracks and unusual bulges, could be found as shown in Figure 3.In this study, an earthquake is believed to be the primary trigger factor of landslide in Sochi ski-jump area, rather than rainfall.Therefore, rainfall-induced landslide has not been considered yet.In recent years, some construction related to tourism has been planned.Some work has been carried out, so that the topography relief and the geological deposits have been changed.Some loss Quaternary deposits like sand and clay have been cleared out.Regarding geological conditions (Figure 2), most of the area is covered by Jurassic mudstone as outcropping bedrock.The upper layer of mudstone bedrock is widespread weathered mudstone.The intact mudstone of bedrock is of low permeability (0.001 m/d), while the weathered layer is more permeable (1 m/d).The thickness of the weathered mudstone all over the area is about 2-3 m.Due to the low strength of weathered layer, the area is supposed to be prone to landslide.Through the field investigation, some shallow landslide warning signs, such as new cracks and unusual bulges, could be found as shown in Figure 3.In this study, an earthquake is believed to be the primary trigger factor of landslide in Sochi ski-jump area, rather than rainfall.Therefore, rainfall-induced landslide has not been considered yet.

Geotectonic and Seismic Conditions of the Region
In terms of seismic hazard, the study area is of high seismic activity.According to the existing map of general seismic zoning (GSZ) of the Russian Federation (OSR-97), the territory of Sochi falls in the zone prone to seismic intensity of 8-9 [29].The last strong earthquake, which triggered widespread distributed landslides, was thought to occur about 600 years ago [35].The active faults of the Sochi region include Main Caucasian, Bekisheiskii, Krasnopolyanskii, Monastyrskii and Sochinskii faults [30].All these active faults are reverse faults, except that the Main Caucasian is a thrust fault.According to the empirical relationship between the earthquake magnitude and the displacement on the fault, the fault displacement with an amplitude of 1.2 m corresponds to the magnitude Mw = 7.3, which is the maximum magnitude of the Sochi region [35].A slope instability map of the study area.The main existing landslide warning signs were found in the western part of the area; some also could be found in the southern and northeastern parts.

Geotectonic and Seismic Conditions of the Region
In terms of seismic hazard, the study area is of high seismic activity.According to the existing map of general seismic zoning (GSZ) of the Russian Federation (OSR-97), the territory of Sochi falls in the zone prone to seismic intensity of 8-9 [29].The last strong earthquake, which triggered widespread distributed landslides, was thought to occur about 600 years ago [35].The active faults of the Sochi region include Main Caucasian, Bekisheiskii, Krasnopolyanskii, Monastyrskii and Sochinskii faults [30].All these active faults are reverse faults, except that the Main Caucasian is a thrust fault.According to the empirical relationship between the earthquake magnitude and the displacement on the fault, the fault displacement with an amplitude of 1.2 m corresponds to the magnitude Mw = 7.3, which is the maximum magnitude of the Sochi region [35].
As for the ski-jump area, the potential earthquake may be sourced from the Main Caucasian fault zone, Bekisheiskii fault zone, and Krasnopolyanskii fault zone (Figure 4a).These active faults in the study zone generally run from east-southwest to west-northeast and they run parallel to each other.As mentioned above, the Main Caucasian fault is a thrust fault; the Bekisheiskii and Krasnopolyanskii faults are reverse faults.In this paper, we consider the Krasnopolyanskii fault as the source of the potential earthquakes, which may represent the most negative seismic condition (Figure 4b).As for the ski-jump area, the potential earthquake may be sourced from the Main Caucasian fault zone, Bekisheiskii fault zone, and Krasnopolyanskii fault zone (Figure 4a).These active faults in the study zone generally run from east-southwest to west-northeast and they run parallel to each other.As mentioned above, the Main Caucasian fault is a thrust fault; the Bekisheiskii and Krasnopolyanskii faults are reverse faults.In this paper, we consider the Krasnopolyanskii fault as the source of the potential earthquakes, which may represent the most negative seismic condition (Figure 4b).A slope instability map of the study area.The main existing landslide warning signs were found in the western part of the area; some also could be found in the southern and northeastern parts.

Geotectonic and Seismic Conditions of the Region
In terms of seismic hazard, the study area is of high seismic activity.According to the existing map of general seismic zoning (GSZ) of the Russian Federation (OSR-97), the territory of Sochi falls in the zone prone to seismic intensity of 8-9 [29].The last strong earthquake, which triggered widespread distributed landslides, was thought to occur about 600 years ago [35].The active faults of the Sochi region include Main Caucasian, Bekisheiskii, Krasnopolyanskii, Monastyrskii and Sochinskii faults [30].All these active faults are reverse faults, except that the Main Caucasian is a thrust fault.According to the empirical relationship between the earthquake magnitude and the displacement on the fault, the fault displacement with an amplitude of 1.2 m corresponds to the magnitude Mw = 7.3, which is the maximum magnitude of the Sochi region [35].
As for the ski-jump area, the potential earthquake may be sourced from the Main Caucasian fault zone, Bekisheiskii fault zone, and Krasnopolyanskii fault zone (Figure 4a).These active faults in the study zone generally run from east-southwest to west-northeast and they run parallel to each other.As mentioned above, the Main Caucasian fault is a thrust fault; the Bekisheiskii and Krasnopolyanskii faults are reverse faults.In this paper, we consider the Krasnopolyanskii fault as the source of the potential earthquakes, which may represent the most negative seismic condition (Figure 4b).
where I m is the MSK-64 intensity, M is magnitude, D is distance between the hypocenter and the object location in kms, and PGA is the maximum horizontal acceleration in gs.
Note that D may respond to the hypocentral depth when the object is in the fault zone.In this study, the area was assumed located upon the earthquake source, as shown in Figure 4b. the distance D between the epicenter and the object location refers to the depth of the hypocenter.The depth of the hypocenter is assumed to range from 15 km to 30 km.According to the Equations ( 1) and ( 2), the assumed I m and PGA of the study area were summarized, as shown in Table 1.The study aims to analyze shallow landslide susceptibility of the area under different seismic conditions.As the previous text shows, the conditions include earthquakes with hypocenters of different depth.

Topographic Data and Slope Map
A 5-m digital elevation model (DEM) with the 0.5-m error obtained from the aerial scanning data was used for topographic analysis.Topographic features were presented in ArcGIS (Figure 5).The slope map was created by using a simple algorithm from the DEM that compares the elevations of adjacent cells and calculates the maximum slope.Figure 5 showed the slopes with large slope angles (>40 • ), including two northern artificial scarps.The following analysis aims to assess natural slope susceptibility; therefore, the stability results of the northern artificial slopes will be excluded.According to the slope distribution map, the statistics of the study area with different slope angles were analyzed, as shown in Figure 6.The slope-angle analysis showed that the steepest slope is about 45 • and the most distributed slope is between 20-30 • .Therefore, the variation of the slope angles was achieved, which can be applied to classify the slope angle of slope models.

Slope Models
According to slope-angle distribution (Figure 6), six simplified mudstone slope models were set up in the GeoStudio software, as shown in Figure 7.The slope angles of models were set to be 20  , 40 • , and 45 • .The slope instability of the slope angle below 20 • was not considered due to the results of initial calculation, in which the Fs of the slope of slope angle below 20 • is hardly below 1.0 even in the most negative seismic condition.The slope model height was 20 m; the thickness of weathered mudstone was 2-3 m.Through the field investigation, the steepest slopes have the minimum thickness of weathered mudstone, while the thickness of weathered mudstone in the slope with a slope angle 20 • is about 3 m.Therefore, the thickness of weathered mudstone in 45 • slope was set to be 2 m, and the weathered layer thickness in 20 • slope −3 m.The weathered mudstone thickness value in other slope models (with slope angle 25 • , 30 • , 35 • , 40 • ) was assigned by linear interpolation.In addition, the tension cracks were set in the weathered mudstone.These models aim to analyze the shallow slope instability.

Slope Models
According to slope-angle distribution (Figure 6), six simplified mudstone slope models were set up in the GeoStudio software, as shown in Figure 7.The slope angles of models were set to be 20°, 25°, 30°, 35°, 40°, and 45°.The slope instability of the slope angle below 20° was not considered due to the results of initial calculation, in which the Fs of the slope of slope angle below 20° is hardly below 1.0 even in the most negative seismic condition.The slope model height was 20 m; the thickness of weathered mudstone was 2-3 m.Through the field investigation, the steepest slopes have the minimum thickness of weathered mudstone, while the thickness of weathered mudstone in the slope with a slope angle 20° is about 3 m.Therefore, the thickness of weathered mudstone in 45° slope was

Slope Models
According to slope-angle distribution (Figure 6), six simplified mudstone slope models were set up in the GeoStudio software, as shown in Figure 7.The slope angles of models were set to be 20°, 25°, 30°, 35°, 40°, and 45°.The slope instability of the slope angle below 20° was not considered due to the results of initial calculation, in which the Fs of the slope of slope angle below 20° is hardly below 1.0 even in the most negative seismic condition.The slope model height was 20 m; the thickness of weathered mudstone was 2-3 m.Through the field investigation, the steepest slopes have the minimum thickness of weathered mudstone, while the thickness of weathered mudstone in the slope with a slope angle 20° is about 3 m.Therefore, the thickness of weathered mudstone in 45° slope was

Physical-Mechanical Parameters of Mudstone
Based on the experimental result, the typical physical-mechanical parameters of weathered mudstone are taken as follows: the unit weight is 19 kN/m 3 , cohesion is 10 kPa and friction angle is 27 • .As for intact mudstone, the unit weight is taken to be 22 kN/m 3 , cohesion-25 kPa and friction angle is 28 • .Note that the intact mudstone refers to the lower layer close to the weathered mudstone.The above parameters were, then, used as input in numerical calculation of slope stability analysis.
set to be 2 m, and the weathered layer thickness in 20° slope −3 m.The weathered mudstone thickness value in other slope models (with slope angle 25°, 30°, 35°, 40°) was assigned by linear interpolation.In addition, the tension cracks were set in the weathered mudstone.These models aim to analyze the shallow slope instability.

Physical-Mechanical Parameters of Mudstone
Based on the experimental result, the typical physical-mechanical parameters of weathered mudstone are taken as follows: the unit weight is 19 kN/m 3 , cohesion is 10 kPa and friction angle is 27°.As for intact mudstone, the unit weight is taken to be 22 kN/m 3 , cohesion-25 kPa and friction angle is 28°.Note that the intact mudstone refers to the lower layer close to the weathered mudstone.The above parameters were, then, used as input in numerical calculation of slope stability analysis.

Spencer Method of Slope Stability Analysis
The limit equilibrium method for analysis of slope stability has been applied in geotechnical engineering for many decades.According to the method, the landslide body is divided into many slices in order to analyze mechanical state.The main concept of limit equilibrium method is to get safety factor with respect to horizontal force equilibrium equation or moment equilibrium equation.The safety factor is calculated by comparing the resistance with the sliding strength of the slope.When the safety factor is above 1.0, the slope is supposed to be stable.When the safety factor is below 1.0, the slope is believed to be unstable.The slope is more stable as the safety factor is bigger.
The safety factor equation with respect to force equilibrium is: where ci is cohesion of i-th slice, φi is friction angle of i-th slice, α is inclination of slice base, Ni is slice base normal force and li is geometric parameter.The safety factor equation with respect to moment equilibrium is:

Spencer Method of Slope Stability Analysis
The limit equilibrium method for analysis of slope stability has been applied in geotechnical engineering for many decades.According to the method, the landslide body is divided into many slices in order to analyze mechanical state.The main concept of limit equilibrium method is to get safety factor with respect to horizontal force equilibrium equation or moment equilibrium equation.The safety factor is calculated by comparing the resistance with the sliding strength of the slope.When the safety factor is above 1.0, the slope is supposed to be stable.When the safety factor is below 1.0, the slope is believed to be unstable.The slope is more stable as the safety factor is bigger.
The safety factor equation with respect to force equilibrium is: where c i is cohesion of i-th slice, ϕ i is friction angle of i-th slice, α is inclination of slice base, N i is slice base normal force and l i is geometric parameter.The safety factor equation with respect to moment equilibrium is: where G i is i-th slice weight; R, x and f are geometric parameters; others are the same as Equation (3).Spencer developed a method accessing safety factor satisfying both force and moment equilibrium [28].In the Spencer method, interslice shear to normal ratio is a constant, which makes the two safety factors from Equations ( 3) and (4) equal.As a type of limit equilibrium method, the Spencer method has mathematically more rigorous formulations rather than traditional methods like the Janbu method [36] and Bishop [37] method.The detailed content can be referred to Spencer's study [28].In a word, the safety factor can be achieved through the Spencer method based on the material unit weight, cohesion, friction angle, and slope angle (the slope model geometry).

Pseudostatic Analysis
The pseudostatic analysis shows the seismic effects by an earthquake's ground accelerations that leads to inertial forces.The direction of seismic forces can be horizontal and vertical at the centroid of each slice while using the Spencer method (Figure 8).The horizontal and vertical seismic forces are defined as: where a h and a v are horizontal and vertical pseudostatic accelerations, g is the gravitational acceleration constant, G is slice weight, K h and K v are horizontal and vertical seismic coefficients.
material unit weight, cohesion, friction angle, and slope angle (the slope model geometry).

Pseudostatic Analysis
The pseudostatic analysis shows the seismic effects by an earthquake's ground accelerations that leads to inertial forces.The direction of seismic forces can be horizontal and vertical at the centroid of each slice while using the Spencer method (Figure 8).The horizontal and vertical seismic forces are defined as: where ah and av are horizontal and vertical pseudostatic accelerations, g is the gravitational acceleration constant, G is slice weight, Kh and Kv are horizontal and vertical seismic coefficients.In this study, the selection of seismic coefficients was based on Eurocode 8 [38]: ah is considered as half of PGA, while, av is equal to 0.33 ah.Therefore, Kh and Kv can be computed on the basis of PGA values responding to different earthquake conditions, as shown in Table 1.

Newmark Deformation Method
In 1965, Newmark introduced a method of analysis for estimating deformations of an artificial embankment dam owing to strong shaking caused by an earthquake [31].This method was called the Newmark method and then used in natural slope analysis.At present, the Newmark method is usually applied to evaluate the performance of landslides induced by an earthquake.The main concept is to estimate accumulated displacement of the slope through doubly integrating the acceleration over the critical acceleration in the sliding block, as shown in Figure 9.In this study, the selection of seismic coefficients was based on Eurocode 8 [38]: a h is considered as half of PGA, while, a v is equal to 0.33 a h .Therefore, K h and K v can be computed on the basis of PGA values responding to different earthquake conditions, as shown in Table 1.

Newmark Deformation Method
In 1965, Newmark introduced a method of analysis for estimating deformations of an artificial embankment dam owing to strong shaking caused by an earthquake [31].This method was called the Newmark method and then used in natural slope analysis.At present, the Newmark method is usually applied to evaluate the performance of landslides induced by an earthquake.The main concept is to estimate accumulated displacement of the slope through doubly integrating the acceleration over the critical acceleration in the sliding block, as shown in Figure 9.
The sliding block has a critical acceleration, a c , which represents the threshold to overcome resistance to starting to slide.Newmark pointed out that the critical acceleration is a function of the static safety factor and slope angle, defined as: where a c is the critical acceleration, g is gravity acceleration, Fs is the static safety factor, and α is slope angle.
resistance to starting to slide.Newmark pointed out that the critical acceleration is a function of the static safety factor and slope angle, defined as: where ac is the critical acceleration, g is gravity acceleration, Fs is the static safety factor, and α is slope angle.
where Dn is Newmark displacement in centimeters, Ia is the Arias intensity in meters per second, and ac is the critical acceleration in gs.Since Russia adopts MSK-64 intensity, the regression equation of Newmark displacement requires to be a function of critical acceleration and MSK-64 intensity.Based on the study of Margottini et al. [39], Arias intensity can be transformed into MSK-64 intensity through an equation as follows: where Ia is the Arias intensity, Im is the MSK-64 intensity Combining Equations ( 8) and ( 9 where D n is Newmark displacement in centimeters, I a is the Arias intensity in meters per second, and a c is the critical acceleration in gs.Since Russia adopts MSK-64 intensity, the regression equation of Newmark displacement requires to be a function of critical acceleration and MSK-64 intensity.Based on the study of Margottini et al. [39], Arias intensity can be transformed into MSK-64 intensity through an equation as follows: lg where I a is the Arias intensity, I m is the MSK-64 intensity. Combining Equations ( 8) and ( 9), we can get the regression equation of Newmark displacement as a function of critical acceleration and MSK-64 intensity: As mentioned in the previous context, four different seismic conditions are considered with respect to magnitude 7.3 and earthquake hypocentral depth 15-30 km.The MSK-64 intensity can be calculated by Equation (1), as shown in Table 1.Once intensity is known, Newmark displacements D n in these conditions can be expressed as follows, As for earthquake hypocentral depth 15 km: As for earthquake hypocentral depth 20 km: As for earthquake hypocentral depth 25 km: As for earthquake hypocentral depth 30 km:

Relationships of Safety Factor to Slope Angle
The variation of physical-mechanical parameters may influence the results of Fs.However, the variation of parameters is limited for the study area.Therefore, the mean values were selected to use in numerical modelling.As the physical-mechanical parameters were fixed, the relationship of safety factor to slope angle was specifically analyzed.
Numerical modelling was conducted for slope models from 6 different slope angles (20  , and 45 • ), as shown in Figure 7. Static and seismic safety factors were computed using Spencer method.Applying pseudostatic analysis, seismic slope stability modelling included 4 conditions, in which seismic coefficients F h and F v were used, referring to Table 1.
The results of safety factor in either static or seismic states, plotted in Figure 10, turned out to fit polynomial curves.It could be found that the regression equations were well constrained (R 2 > 0.96).The resulting regression equations between the safety factor and slope angle are expressed as shown in Table 2.

Static and Seismic Safety Factor Mapping
Since the slope map (Figure 5) derived from DEM has been achieved, the safety factor maps can be drawn based on the regression equations between the safety factor and slope angle.The mapping process was carried out using the raster calculation tool in ArcGIS program.
Figure 11 shows the static safety factor map. From the figure, we could see the western part of the area suffered the lowest safety factors, the value of which is near 1.0.As mentioned above, slope stability or instability is based on the interplay between driving forces and resisting forces: when driving forces overcome resisting forces (Fs < 1.0), the slope is unstable and it results in slope failure.Since the slope map was made, the relationship of the safety factor to slope angle helped to conduct safety factor maps.Then, the equations were implemented based on the ArcGIS program in order to produce safety factor maps.

Static and Seismic Safety Factor Mapping
Since the slope map (Figure 5) derived from DEM has been achieved, the safety factor maps can be drawn based on the regression equations between the safety factor and slope angle.The mapping process was carried out using the raster calculation tool in ArcGIS program.
Figure 11 shows the static safety factor map. From the figure, we could see the western part of the area suffered the lowest safety factors, the value of which is near 1.0.As mentioned above, slope stability or instability is based on the interplay between driving forces and resisting forces: when driving forces overcome resisting forces (Fs < 1.0), the slope is unstable and it results in slope failure.According to the field investigation, the slopes possessing the lowest safety factors turned out to be the most vulnerable to slope failure, with landslide warning signs investigated.Because the safety factor is above 1.0, the entire slope failure has not occurred yet.Compared with the landslide warning sign locations in Figure 3, the static safety factor map of Figure 11 corresponds to the actual potential landslide map, which could verify the model validation.

Static and Seismic Safety Factor Mapping
Since the slope map (Figure 5) derived from DEM has been achieved, the safety factor maps can be drawn based on the regression equations between the safety factor and slope angle.The mapping process was carried out using the raster calculation tool in ArcGIS program.
Figure 11 shows the static safety factor map. From the figure, we could see the western part of the area suffered the lowest safety factors, the value of which is near 1.0.As mentioned above, slope stability or instability is based on the interplay between driving forces and resisting forces: when driving forces overcome resisting forces (Fs < 1.0), the slope is unstable and it results in slope failure.According to the field investigation, the slopes possessing the lowest safety factors turned out to be the most vulnerable to slope failure, with landslide warning signs investigated.Because the safety factor is above 1.0, the entire slope failure has not occurred yet.Compared with the landslide warning sign locations in Figure 3, the static safety factor map of Figure 11 corresponds to the actual potential landslide map, which could verify the model validation.
The seismic safety factor maps were drawn with respect to hypocentral depth 30 km (Figure 12a), hypocentral depth 25 km (Figure 12b), hypocentral depth 20 km (Figure 12c), and hypocentral depth 15 km (Figure 12d).Results of seismic safety factor maps indicate that the area of Fs below 1.0 expands, as hypocentral depth decreases: when hypocentral depth is 30 km, the slope instability occurs only in small parts of western and eastern slopes; when hypocentral depth is 15 km, the slope failure appears all over the area except for northern and central area.The seismic safety factor maps were drawn with respect to hypocentral depth 30 km (Figure 12a), hypocentral depth 25 km (Figure 12b), hypocentral depth 20 km (Figure 12c), and hypocentral depth 15 km (Figure 12d).Results of seismic safety factor maps indicate that the area of Fs below 1.0 expands, as hypocentral depth decreases: when hypocentral depth is 30 km, the slope instability occurs only in small parts of western and eastern slopes; when hypocentral depth is 15 km, the slope failure appears all over the area except for northern and central area.

Newmark Deformation Mapping
As indicated above, conducting a Newmark displacement map requires the critical acceleration map according to Equations ( 1)-( 4).First, the critical acceleration map was made, based on a function of the static safety factor and slope angle, the values of which had been known.The critical acceleration map is shown in Figure 13.Finally, the predicted Newmark displacements under potential earthquakes were estimated in each raster cell, which made up the Newmark displacement map (Figure 14).From the Newmark displacement maps, we can see that the largest displacement (>10 cm) occurs in the western slopes, correlated with the part of the lowest safety factors.The result also indicates that Newmark displacement increases, as the hypocentral depth decreases.Under condition of earthquake with hypocentral depth of 15 km, the entire area suffers the ground movement, however, most parts of the area will move less than 0.35 cm.The most dangerous slopes, which may undergo the largest displacement, are believed to be the western part of the area.

Newmark Deformation Mapping
As indicated above, conducting a Newmark displacement map requires the critical acceleration map according to Equations ( 1)-( 4).First, the critical acceleration map was made, based on a function of the static safety factor and slope angle, the values of which had been known.The critical acceleration map is shown in Figure 13.Finally, the predicted Newmark displacements under potential earthquakes were estimated in each raster cell, which made up the Newmark displacement map (Figure 14).From the Newmark displacement maps, we can see that the largest displacement (>10 cm) occurs in the western slopes, correlated with the part of the lowest safety factors.The result also indicates that Newmark displacement increases, as the hypocentral depth decreases.Under condition of earthquake with hypocentral depth of 15 km, the entire area suffers the ground movement, however, most parts of the area will move less than 0.35 cm.The most dangerous slopes, which may undergo the largest displacement, are believed to be the western part of the area.

Newmark Deformation Mapping
As indicated above, conducting a Newmark displacement map requires the critical acceleration map according to Equations ( 1)-( 4).First, the critical acceleration map was made, based on a function of the static safety factor and slope angle, the values of which had been known.The critical acceleration map is shown in Figure 13.Finally, the predicted Newmark displacements under potential earthquakes were estimated in each raster cell, which made up the Newmark displacement map (Figure 14).From the Newmark displacement maps, we can see that the largest displacement (>10 cm) occurs in the western slopes, correlated with the part of the lowest safety factors.The result also indicates that Newmark displacement increases, as the hypocentral depth decreases.Under condition of earthquake with hypocentral depth of 15 km, the entire area suffers the ground movement, however, most parts of the area will move less than 0.35 cm.The most dangerous slopes, which may undergo the largest displacement, are believed to be the western part of the area.

Discussion and Conclusions
The method combining GIS and numerical modelling were proposed.After numerical calculation of 6 slope models with different slope angles, the relationships between safety factor and slope angle were acquired.From the relationships of safety factor to slope angle, it can be concluded that safety factor decreases as slope angle increases.Also, safety factor is proved to decrease as seismic effect (Kc) increases.The results are easy to understand, however, we found that under either static or seismic conditions the relationships of safety factor to slope angle can be fitted by polynomial curves.Using the regression equations, it was easier to conduct safety factor maps in GIS.
According to the static safety factor map, the western slopes of the area were believed to be the most vulnerable to slope failure.The result of the lowest Fs in the static safety factor map is correlated with field performance, where the landslide warning signs occurred.As the hypocentral depth decreases, the seismic effect increases.The seismic safety factor maps give an approach to investigate the safety factor variation due to the change of seismic effect.
The transformed Newmark-displacement equations based on MSK-64 intensity were presented.The result shows that the Newmark displacement is supposed to increase, as the seismic effect increases.In addition, the slopes with the largest Newmark displacements correspond to the ones with the lowest safety factors.Newmark displacement maps provide an index to predict the landslide susceptibility.
The shortcomings of this study are lack of comparison of the predicted Newmark displacement with the actual landslide performance.This is due to the fact that the study area has been rebuilt and the last earthquake inducing the widespread landslides in this region occurred 600 years ago.Future work will expand the scope of the study area, considering the inventory of landslides triggered by historical earthquakes.

Discussion and Conclusions
The method combining GIS and numerical modelling were proposed.After numerical calculation of 6 slope models with different slope angles, the relationships between safety factor and slope angle were acquired.From the relationships of safety factor to slope angle, it can be concluded that safety factor decreases as slope angle increases.Also, safety factor is proved to decrease as seismic effect (K c ) increases.The results are easy to understand, however, we found that under either static or seismic conditions the relationships of safety factor to slope angle can be fitted by polynomial curves.Using the regression equations, it was easier to conduct safety factor maps in GIS.
According to the static safety factor map, the western slopes of the area were believed to be the most vulnerable to slope failure.The result of the lowest Fs in the static safety factor map is correlated with field performance, where the landslide warning signs occurred.As the hypocentral depth decreases, the seismic effect increases.The seismic safety factor maps give an approach to investigate the safety factor variation due to the change of seismic effect.
The transformed Newmark-displacement equations based on MSK-64 intensity were presented.The result shows that the Newmark displacement is supposed to increase, as the seismic effect increases.In addition, the slopes with the largest Newmark displacements correspond to the ones with the lowest safety factors.Newmark displacement maps provide an index to predict the landslide susceptibility.
The shortcomings of this study are lack of comparison of the predicted Newmark displacement with the actual landslide performance.This is due to the fact that the study area has been rebuilt and the last earthquake inducing the widespread landslides in this region occurred 600 years ago.Future work will expand the scope of the study area, considering the inventory of landslides triggered by historical earthquakes.

15 Figure 1 .
Figure 1.Location of the study area: (a) the area is located in northern Caucasus on the east Black Sea coast (based on Google Earth); (b) aerial scanning photo of the study area with the Sochi Olympic skijump facilities.

Figure 2 .
Figure 2. A typical slope profile of the study area.The main geological bedrock is Jurassic mudstone, including two layers: (a) the upper layer is weathered mudstone with thickness of around 2 m; (b) the lower layer is intact mudstone.

Figure 1 .
Figure 1.Location of the study area: (a) the area is located in northern Caucasus on the east Black Sea coast (based on Google Earth); (b) aerial scanning photo of the study area with the Sochi Olympic ski-jump facilities.

Figure 1 .
Figure 1.Location of the study area: (a) the area is located in northern Caucasus on the east Black Sea coast (based on Google Earth); (b) aerial scanning photo of the study area with the Sochi Olympic skijump facilities.

Figure 2 .
Figure 2. A typical slope profile of the study area.The main geological bedrock is Jurassic mudstone, including two layers: (a) the upper layer is weathered mudstone with thickness of around 2 m; (b) the lower layer is intact mudstone.

Figure 2 .
Figure 2. A typical slope profile of the study area.The main geological bedrock is Jurassic mudstone, including two layers: (a) the upper layer is weathered mudstone with thickness of around 2 m; (b) the lower layer is intact mudstone.

Figure 3 .
Figure 3.A slope instability map of the study area.The main existing landslide warning signs were found in the western part of the area; some also could be found in the southern and northeastern parts.

Figure 4 .
Figure 4.The tectonic schemes of the region: (a) the main active faults around the study area; (b) the schematic plot of the active faults with assumed earthquake focuses of different depth.

Figure 3 .
Figure3.A slope instability map of the study area.The main existing landslide warning signs were found in the western part of the area; some also could be found in the southern and northeastern parts.

Figure 3 .
Figure 3.A slope instability map of the study area.The main existing landslide warning signs were found in the western part of the area; some also could be found in the southern and northeastern parts.

Figure 4 .
Figure 4.The tectonic schemes of the region: (a) the main active faults around the study area; (b) the schematic plot of the active faults with assumed earthquake focuses of different depth.

Figure 4 .
Figure 4.The tectonic schemes of the region: (a) the main active faults around the study area; (b) the schematic plot of the active faults with assumed earthquake focuses of different depth.As for the Caucasus region, the empirical relation of intensity (I m ) to magnitude (M) and hypocentral distance (D) and the empirical relation of peak ground acceleration (PGA) to intensity (I m ) are shown as follows [14]: I m = 1.5M − 4.7 lgD + 4.0, (1) lgPGA = 0.301 I m − 3.1,(2)

Figure 5 .
Figure 5. Slope map derived from the digital elevation model (DEM) of the study area (location shown in Figure1).The steepest slope of 45° appears in the western part of the area.The two artificial scarps, also showing steep slopes, are located in the northern part of the area (see the legend of Figure3).

Figure 6 .
Figure 6.The statistics concerning area distribution of different slope intervals: the steepest slope is about 45° and the most distributed slope is between 20-30°.

Figure 5 .
Figure 5. Slope map derived from the digital elevation model (DEM) of the study area (location shown in Figure 1).The steepest slope of 45 • appears in the western part of the area.The two artificial scarps, also showing steep slopes, are located in the northern part of the area (see the legend of Figure 3).

15 Figure 5 .
Figure 5. Slope map derived from the digital elevation model (DEM) of the study area (location shown in Figure1).The steepest slope of 45° appears in the western part of the area.The two artificial scarps, also showing steep slopes, are located in the northern part of the area (see the legend of Figure3).

Figure 6 .
Figure 6.The statistics concerning area distribution of different slope intervals: the steepest slope is about 45° and the most distributed slope is between 20-30°.

Figure 6 .
Figure 6.The statistics concerning area distribution of different slope intervals: the steepest slope is about 45 • and the most distributed slope is between 20-30 • .

Figure 7 .
Figure 7. Slope models for stability analysis from 6 different slope angles (20°, 25°, 30°, 35°, 40°, and 45°) with static safety factors.The orange part responds to intact mudstone; the yellow part responds to weathered mudstone; the green part means the calculated landslide body.

Figure 7 .
Figure 7. Slope models for stability analysis from 6 different slope angles (20 • , 25 • , 30 • , 35 • , 40 • , and 45 • ) with static safety factors.The orange part responds to intact mudstone; the yellow part responds to weathered mudstone; the green part means the calculated landslide body.

Figure 8 .
Figure 8. Schematic map showing the directions of seismic forces in the slope slice in pseudostatic analysis.Note that for briefly presenting the seismic forces, the interslice forces and forces under the slice base are not shown in the stretch.

Figure 8 .
Figure 8. Schematic map showing the directions of seismic forces in the slope slice in pseudostatic analysis.Note that for briefly presenting the seismic forces, the interslice forces and forces under the slice base are not shown in the stretch.

Figure 9 .
Figure 9. Schematic diagram indicating integration of acceleration to determine landslide block displacement (adapted from [16]).(A) Ground motion acceleration versus time with critical acceleration (horizontal dashed line).(B) Velocity of the sliding-block versus time.(C) Displacement of the sliding-block versus time.

10 )Figure 9 .
Figure 9. Schematic diagram indicating integration of acceleration to determine landslide block displacement (adapted from [16]).(A) Ground motion acceleration versus time with critical acceleration (horizontal dashed line).(B) Velocity of the sliding-block versus time.(C) Displacement of the sliding-block versus time.To get Newmark displacement requires knowing the static safety factor, the slope angle and selecting an earthquake ground strong-motion record.After considering a large number of strong-motion records, Jibson and et al. developed an empirical regression equation to estimate Newmark displacement [16].The regression equation of Newmark displacement is a function of critical acceleration and Arias intensity, expressed as: lg D n = 1.521 lg I a − 1.993 lg a c − 1.546,(8) ISPRS Int.J. Geo-Inf.2019, 8, 148 11 of 15

Figure 10 .
Figure 10.Relationships between safety factor and slope angle with respect to static condition and four seismic conditions.

Figure 10 .
Figure 10.Relationships between safety factor and slope angle with respect to static condition and four seismic conditions.

Figure 10 .
Figure 10.Relationships between safety factor and slope angle with respect to static condition and four seismic conditions.

Figure 11 .
Figure 11.Static safety factor map of the study area (location shown in Figure 1).Figure 11.Static safety factor map of the study area (location shown in Figure 1).

Figure 11 .
Figure 11.Static safety factor map of the study area (location shown in Figure 1).Figure 11.Static safety factor map of the study area (location shown in Figure 1).

Figure 12 .
Figure 12.Seismic safety factor maps of the study area (location shown in Figure 1), with respect to (a) hypocentral depth 30 km, (b) hypocentral depth 25 km, (c) hypocentral depth 20 km, and (d) hypocentral depth 15 km.

Figure 13 .
Figure 13.The critical acceleration map related to Newmark deformation calculation.The values represent the threshold of ground acceleration to overcome resistance for a landslide block to start to slide.

Figure 12 .
Figure 12.Seismic safety factor maps of the study area (location shown in Figure 1), with respect to (a) hypocentral depth 30 km, (b) hypocentral depth 25 km, (c) hypocentral depth 20 km, and (d) hypocentral depth 15 km.

Figure 13 .
Figure 13.The critical acceleration map related to Newmark deformation calculation.The values represent the threshold of ground acceleration to overcome resistance for a landslide block to start to slide.

Figure 13 .
Figure 13.The critical acceleration map related to Newmark deformation calculation.The values represent the threshold of ground acceleration to overcome resistance for a landslide block to start to slide.

Figure 14 .
Figure 14.The predicted Newmark displacement maps: (a) in case of M 7.3 earthquake with hypocentral depth 30 km; (b) in the case of M 7.3 earthquake with hypocentral depth 25 km; (c) in case of M 7.3 earthquake with hypocentral depth 20 km; (d) in case of M 7.3 earthquake with hypocentral depth of 15 km.

Figure 14 .
Figure 14.The predicted Newmark displacement maps: (a) in case of M 7.3 earthquake with hypocentral depth 30 km; (b) in the case of M 7.3 earthquake with hypocentral depth 25 km; (c) in case of M 7.3 earthquake with hypocentral depth 20 km; (d) in case of M 7.3 earthquake with hypocentral of 15 km.

Table 1 .
The seismic parameters related to maximum magnitude in study area.

Table 2 .
The relationships of safety factor to slope angle in deferent conditions.