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

Abstract

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.

1. 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² > 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.

2. Study Area

2.1. 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² (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.

2.2. 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 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)

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 seismic parameters related to maximum magnitude in study area.

Magnitude, M	Earthquake Hypocentral Depth, D (km)	MSK-64 Intensity, Im	Peak Ground Horizontal Acceleratio, PGA (g)	Horizontal Seismic Coefficient, Kh	Vertical Seismic Coefficient, Kv
7.3	15	9.4	0.44	0.22	0.07
7.3	20	8.8	0.30	0.15	0.05
7.3	25	8.4	0.20	0.1	0.033
7.3	30	8.0	0.16	0.08	0.026

.

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.

3. Data and Methods

3.1. 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.

3.2. Slope Stability Numerical Modelling

3.2.1. 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 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.

3.2.2. 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³, cohesion is 10 kPa and friction angle is 27°. As for intact mudstone, the unit weight is taken to be 22 kN/m³, 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.

3.2.3. 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:

$${\mathrm{F}}_{\mathrm{s}-\mathrm{f}}=\frac{{\displaystyle \sum (N_i\tan {\varphi }_i\cos \alpha +c_i{l}_i\cos \alpha })}{{\displaystyle \sum N_i\sin \alpha }},$$

(3)

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:

$${\mathrm{F}}_{\mathrm{s}-\mathrm{m}}=\frac{({\displaystyle \sum N_i\tan {\varphi }_i+{\displaystyle \sum c_i{l}_i}})R}{{\displaystyle \sum G_{i}x-{\displaystyle \sum N_{i}f}}},$$

(4)

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).

3.2.4. 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:

$$F_h=\frac{a_h}{g}G=K_{h}G,$$

(5)

$$F_v=\frac{a_v}{g}G=K_{v}G,$$

(6)

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.

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.

3.2.5. 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:

a_c = (Fs − 1)·g sinα,

(7)

where a_c is the critical acceleration, g is gravity acceleration, Fs is the static safety factor, and α is slope angle.

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)

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 I_a = −1.079 + 0.31 I_m,

(9)

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:

lg D_n = 0.162 I_m − 1.993 log a_c − 3.187,

(10)

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:

lg D_n = −1.993 log a_c − 1.67,

(11)

As for earthquake hypocentral depth 20 km:

lg D_n = −1.993 log a_c − 1.76,

(12)

As for earthquake hypocentral depth 25 km:

lg D_n = −1.993 log a_c − 1.83,

(13)

As for earthquake hypocentral depth 30 km:

lg D_n = −1.993 log a_c − 1.89,

(14)

4. Results

4.1. 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°, 25°, 30°, 35°, 40°, 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² > 0.96). The resulting regression equations between the safety factor and slope angle are expressed as shown in

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

Condition of Numerical Modelling	Magnitude, M	Earthquake Hypocentral Depth, D (km)	Regression Equations between Safety Factor (Fs) and Slope Angle (α)
Static	None	-	Fs = 0.0011α2 − 0.110α + 3.6686; R2 = 0.984 *
Seismic	7.3	15	Fs = 0.0005α2 − 0.051α + 1.9516; R2 = 0.993
Seismic	7.3	20	Fs = 0.0005α2 − 0.054α + 2.1955; R2 = 0.975
Seismic	7.3	25	Fs = 0.0008α2 − 0.075α + 2.6613; R2 = 0.994
Seismic	7.3	30	Fs = 0.0007α2 − 0.073α + 2.7166; R2 = 0.980

* R² means a coefficient of determination of the regression equation.

.

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.

4.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. 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.

4.3. 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.

5. 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.