**Earthquake Environmental E**ff**ects of the 1992 MS7.3 Suusamyr Earthquake, Kyrgyzstan, and Their Implications for Paleo-Earthquake Studies**

### **Christoph Grützner 1,2,\*, Richard Walker 3, Eleanor Ainscoe 3, Austin Elliott <sup>3</sup> and Kanatbek Abdrakhmatov <sup>4</sup>**


Received: 20 May 2019; Accepted: 19 June 2019; Published: 21 June 2019

**Abstract:** Large pre-historical earthquakes leave traces in the geological and geomorphological record, such as primary and secondary surface ruptures and mass movements, which are the only means to estimate their magnitudes. These environmental earthquake effects (EEEs) can be calibrated using recent seismic events and the Environmental Seismic Intensity Scale (ESI2007). We apply the ESI2007 scale to the 1992 MS7.3 Suusamyr Earthquake in the Kyrgyz Tien Shan, because similar studies are sparse in that area and geological setting, and because this earthquake was very peculiar in its primary surface rupture pattern. We analyze literature data on primary and secondary earthquake effects and add our own observations from fieldwork. We show that the ESI2007 distribution differs somewhat from traditional intensity assessments (MSK (Medvedev-Sponheuer-Karnik) and MM (Modified Mercalli)), because of the sparse population in the epicentral area and the spatial distribution of primary and secondary EEEs. However, the ESI2007 scale captures a similar overall pattern of the intensity distribution. We then explore how uncertainties in the identification of primary surface ruptures influence the results of the ESI2007 assignment. Our results highlight the applicability of the ESI2007 scale, even in earthquakes with complex and unusual primary surface rupture patterns.

**Keywords:** earthquake environmental effects; Suusamyr earthquake; Kyrgyzstan; Tien Shan; surface rupture; landslide; digital elevation model (DEM); Structure-from-Motion

#### **1. Introduction**

Paleo-earthquake magnitudes are usually estimated from the lengths of mapped primary surface ruptures and the single-event offsets of geomorphological markers. Empirical relationships between magnitude, primary surface rupture length, and average/maximum offset then allow calculating magnitudes (e.g., References [1,2]). Those relationships are derived from recent surface-rupturing earthquakes with well-constrained magnitudes. Using this approach in paleoseismological studies comes with several limitations: (i) It can only be applied to earthquakes with primary surface ruptures. Large events without significant surface offset like the 2015 Gorkha MW7.8 Earthquake [3] will go unnoticed in the paleoseismological record; (ii) if earthquake recurrence intervals on a fault are long, erosion or sedimentation may have intensely altered the primary surface ruptures or even completely eradicated them (e.g., Reference [4]); and (iii) complex primary surface ruptures on several faults

may not be interpreted as a single earthquake. For example, the 2016 Kaikoura Earthquake in New Zealand with a magnitude of MW7.8 ruptured at least twelve major crustal faults and produced highly variable primary surface ruptures [5]. It is unlikely that such complexity can be understood with paleoseismological methods, and the magnitude of the paleo-earthquake will probably be under-estimated. (iv) Incoherent and anomalous primary surface ruptures can not only be hard to detect, but may also cause problems when it comes to applying the empirical relationships. For example, the 1992 MS7.3 Suusamyr Earthquake, Kyrgyzstan, broke the surface in two short sets of primary ruptures with a 25 km gap in between [6,7].

One approach to overcome the problem of estimating paleo-earthquake magnitudes based solely on primary surface ruptures is the application of the Environmental Seismic Intensity Scale, ESI2007 [8–10]. In contrast to classical intensity scales, such as MSK (Medvedev-Sponheuer-Karnik) or MM (Modified Mercalli), the ESI2007 only uses effects on the environment to assign earthquake intensities. It, therefore, avoids the influence of building styles and the problem of saturation at high intensities, and allows applying the scale to paleo-earthquakes. The conversion of ESI2007 intensities to magnitudes, however, can only be based on a large set of modern case studies, which are currently being collected in the Earthquake Environmental Effects (EEE) catalogue hosted by ISPRA [11]. Preferably, this set of case studies should include earthquakes with different mechanisms, magnitudes, depths, tectonic and geological settings, distributed all across the globe. A glance at the events included in the catalogue so far makes clear that this is not yet the case.

In this paper, we apply the ESI2007 scale to the 1992 Suusamyr Earthquake. There are several reasons for doing this: (i) The majority of entries in the EEE catalogue are currently from Europe and NW South America. Only a few Central Asian earthquakes are included. Thus, this study contributes to extending the entries for central Asia and compressional tectonic settings; (ii) the Suusamyr Earthquake was special in its primary surface rupture pattern, and the exercise of applying the ESI2007 scale to it thereby points out a general problem; (iii) the earthquake occurred in a sparsely populated area, which hampers the application of traditional intensity scales that mainly focus on damage to human-made infrastructure; and (iv) we document secondary cracks of the Suusamyr Earthquake with a high-resolution digital elevation model (DEM) computed from drone aerial imagery and the Structure-from-Motion (SfM) technique [7,12,13].

Here, we first provide an overview of the 1992 Suusamyr Earthquake and review the existing intensity assessments. Then we evaluate the published EEEs from the literature, we include our field observations, and we present an ESI2007 intensity map. Finally, we discuss ESI2007 intensities vs. MM and MSK intensities and the implications of our study for the application of the ESI2007 scale on paleo-earthquake in general.

#### **2. The 1992 Suusamyr Earthquake**

An earthquake with a magnitude of MS7.3 hit the Suusamyr Basin in Kyrgyzstan on 19 August 1992, 02:04 GMT (Figure 1). The Suusamyr Basin is an east-west elongated intramontane basin in the Kyrgyz Tien Shan, surrounded by Paleozoic bedrock. The thickness of the Cenozoic basin fill reaches several hundreds of meters in its widest part near the town of Suusamyr ([14,15] and references therein). A thick Quaternary cover blankets most of the Neogene rocks in the basin center. The basin is bound by E-W striking thrust faults to the north and the south (Figure 1c). North−South shortening is accommodated by the oppositely-vergent thrust. Reference [7] showed that several surface-rupturing earthquakes occurred in the Late Quaternary on these faults, testifying to the tectonic activity of the area. Several ~E-W trending anticlinal ridges in the basin also record ongoing shortening. The Chet Korumdy ridge, which will be discussed in this paper, is a 250 m high, 7 km long ridge paralleling the Suusamyr River (Figure 1). It is made up of steeply north-dipping sediments of Pliocene-Quaternary age [6]. At its western tip, prominent wind gaps record Late Quaternary uplift [7,15,16].

**Figure 1.** (**a**) The Suusamyr Basin in the Kyrgyz Tien Shan. Kyrg.: Kyrgyzstan. (**b**) Overview of the Suusamyr Basin and the 1992 earthquakes. Large white circles: MS7.3 main shock with the Harvard CMT solution [17] and the strongest aftershock (M6.6) located with data from the Kyrgyz regional network [18]. All other aftershocks are from the compilation of Reference [19]. Elevation data: SRTM1. Note the >20 km gap in between the two sets of primary surface ruptures. (**c**) The Suusamyr fault system (from Reference [7]) and location of the Bishkek-Osh Highway.

The Suusamyr Earthquake ruptured an east-west-trending, south-dipping thrust fault in the Suusamyr Basin with the epicenter beneath the Aramsu Range. Reference [18] reports a centroid depth between 5–21 km and a fault dip of 50◦ ± 13◦ based on the teleseismic body-wave inversion of the main shock and an aftershock study using broadband stations and temporal and regional networks. Reference [20] analyzed body waves and surface waves, and reported a hypocenter at a depth of 14 ± 2 km, a fault dip of 49◦ ± 6◦ to the south, and a rake of 105◦ ± 3◦. Within two hours of the main shock, three large aftershocks with magnitudes of MS6.6, MS6.6, and mb6.0 occurred west of the Aramsu Range, but their mechanisms and depths are unknown [18]. The more than 900 aftershocks that were registered in the months after the Suusamyr Earthquake occurred at depths shallower than 15 km. They occurred within a rather narrow, south-dipping zone at depths between 5 and 15 km and showed a more diffuse pattern near the surface [18].

The most surprising effect of the earthquake was its primary surface rupture pattern. Two sets of relatively short primary ruptures were identified. At the eastern tip of Chet Korumdy ridge, a short primary rupture occurred in the Suusamyr River bed [6]. This site is referred to as the 'eastern primary surface ruptures'. The second set of primary surface ruptures formed more than 20 km further west (the 'western primary surface ruptures'), with no hints for surface breaks in between these two (Figure 1). This peculiarity and its implications for ESI2007 intensity assignment will be discussed later in this paper. Apart from these primary surface ruptures, the Suusamyr Earthquake also caused secondary cracks, often referred to as secondary ruptures in the literature. Here we will use the term secondary ruptures to summarize all types of cracks and ruptures that cannot be attributed to primary faulting. The 1992 Suusamyr Earthquake also caused widespread mass movements, mud eruptions, and jumping stones.

#### **3. Macroseismic E**ff**ects of the Suusamyr Earthquake**

The Suusamyr Earthquake resulted in more than 50 deaths, most of which were due to collapsed buildings throughout the valley and mass movements that occurred at Belaldy River [6,21,22]. The epicentral area was surveyed by field crews and from a helicopter during the months following the earthquake. The results of these surveys were reported in great detail in References [6,21–23]. Reference [7] revisited the Suusamyr Basin and reported additional information on the earthquake effects. In the following, we briefly report on the published intensity assessments and review the earthquake environmental effects of the Suusamyr Earthquake in light of the ESI2007 scale.

Reference [6] published a Modified Mercalli (MM) intensity map of the earthquake (Figure 2). They used data from 41 individual sites, assessing the damage to human-made structures and also environmental effects. They report intensities of MM ≥ X for two sites. The first one is an E-W elongated zone around the eastern primary surface ruptures and the Chet Korumdy ridge, the second one is an NNE-SSW elongated zone in the Belaldy River valley region. The first site is characterized by primary and secondary surface ruptures, large landslide volumes, and strong peak ground accelerations (PGA). The second site, however, does not encompass primary or secondary ruptures, but a very large mass movement in the Belaldy Valley. The combined area of MM <sup>≥</sup> X is 62 km2. Intensities of MM = IX were assigned to an ENE-WSW elongated area of 835 km2, surrounding the primary surface ruptures and encompassing the northern part of the area where numerous mass movements occurred. An E-W elongated area of 3000 km2 was assigned intensity MM = VIII. This isoseismal includes nearly all mass movements. The area of intensity MM = VII is NNW-SSE elongated (roughly parallel to the orogenic trend) and encloses 13,000 km2 (Table 1).

**Figure 2.** Earthquake environmental effects (EEEs) of the 1992 Suusamyr Earthquake. Primary surface ruptures are from References [6,7,21]. Secondary ruptures are from References [21,23]. Landslides, the location of the Belaldy rock avalanche, mud eruptions, and jumping rocks are from References [19,23]. Modified Mercalli (MM) intensities are from Reference [6]; Medvedev-Sponheuer-Karnik (MSK) intensities are from Reference [19]. Hillshade is based on SRTM1 elevation data.


**Table 1.** Comparison of the intensity distribution of the 1992 Suusamyr Earthquake.

Intensity data were also compiled within the framework of the SNSF Project No IB7320-110694 (http://www.kyrgyzstan.ethz.ch) by Reference [19], who published Medvedev-Sponheuer-Karnik (MSK) intensities (Figure 2). The shape of the isoseismals differs from the MM intensities reported by Reference [6], but the affected locations and extent are similar. Highest intensities of MSK = X were also assigned to two separate locations, covering 51 km2. The one near the eastern primary surface ruptures is located south of the Suusamyr River and has a shape similar to the MM ≥ X intensities, although shifted to the south by ~1 km. The other one near the Belaldy River valley is elongated in NW-SE direction and also located further south compared to MM ≥ X. Intensities of MSK = IX were assigned to two individual patches with a total of 600 km2. The first one runs parallel to the Suusamyr River and is stretched in E-W direction, the second one is elongated NW-SE (parallel to the orogenic trend), covering the high mountains that border the Suusamyr Valley to the southwest. The zone of intensity MSK = VIII overlaps with the MM = VIII intensities, but is more E-W elongated and also followed the orogenic trend. The total area of MSK = VIII is 3200 km2. No MSK = VII intensities were reported by Reference [19] (Table 1).

In the following, we discuss all earthquake environmental effects described in the literature and present new data from our fieldwork. We assign ESI2007 intensities for all EEEs, which are the basis for our intensity map.

#### *3.1. Primary Surface Ruptures*

The eastern primary surface ruptures occurred at the eastern tip of Chet Korumdy ridge next to the Bishkek-Osh highway (162 km marker; Figures 1–3). As a result of the main primary surface rupture, the Suusamyr River temporarily changed its course. The ruptures were trenched by Reference [24] who confirmed that the fault reached the surface and that the surface deformation is due to thrust motion on an S-dipping fault plane. Reference [21] report a 400 m long, N-facing scarp of up to 2.7 m height and up to 0.3 m of lateral slip. Here, Reference [6] report a fold-scarp geometry and 500 m of primary ruptures with a height of 2.7 m. They estimate a total slip of 4.2 m assuming no horizontal component of motion. Both authors mention numerous rupture segments that also affected the nearby highway, but state that the scarp in the river bed is by far the dominating feature. Furthermore, the authors state that the other segments are likely secondary ruptures. Reference [7] surveyed the eastern primary surface ruptures again in 2015/16 with a drone and produced a high-resolution DEM with the SfM technique. They document a 600 m-long primary surface rupture with up to 3.1 m vertical displacement and a total slip of ~3.6 m that occurred in the Suusamyr River bed. No evidence for the other, smaller scarps on the highway was preserved in 2015/16, due to road repairs, but small secondary scarps were still visible south of the road. Reference [21] reported additional ~4 km of N-facing scarps on the hills further to the east, which could be a continuation of the large main scarp in the river valley. A 300 m wide deformation zone encompasses numerous fault strands with vertical offsets of 0.1–1.05 m height. These scarps face north and are at least in one location formed by a steeply S-dipping fault. References [6,23] do not mention these additional ruptures.

**Figure 3.** Earthquake environmental effects of the 1992 Suusamyr Earthquake. (**a**) Map with the photo locations. (**b**) The central section of the western primary surface ruptures in summer 2015. View to the south. (**c**) The highest section of the eastern primary surface ruptures in summer 2016. View to the west. (**d**) Secondary ruptures on top of the Chet Korumdy ridge in summer 2015. View to the west. (**e**) The 1992 Chet Korumdy landslide that damaged the Bishkek-Osh highway (visible in the background) in summer 2015. View to the west. (**f**) Oblique satellite image of the Belaldy rock avalanche. Source: GoogleEarth/CNES/Airbus (2019), imagery from16 August 2013. View to the northeast. (**g**) Mud eruption site W of the Aramsu range. Note that the eruption site sits at the landslide toe. View to the southwest. (**h**) Mud eruption site in the W of the valley. This site also sits on top of a landslide. View to the west. (**i**) Close-up of the cracks from which the mud erupted.

The western primary surface ruptures occurred on six segments with a total length of 7 km, spanning a 12 km long rupture zone [6,7,21] (Figures 2 and 3). N-facing scarps are typically 0.7–1.2 m high, but locally reach up to 2 m. Reference [7] also report that some scarps show up on pre-1992 satellite imagery, but their pre-1992 height in unknown. This is why the total vertical offset measured from high-resolution DEMs of marker surfaces has to be taken as a maximum value. However, References [6,21] report 0.9–1.8 m high fresh ruptures scarps that formed in 1992, and we take this information as the best available data on their height.

The two sets of primary surface ruptures with their >20 km separation and the unusual scarp length-to-height ratio of the eastern primary surface ruptures can be treated in three different ways in an ESI2007 intensity assignment [8,9]. (i) If we consider the ruptures separately, the western primary surface ruptures correspond to intensity ESI = IX because of their length of 7 km. Using the 12 km long rupture zone reported by Reference [7] would correspond to ESI = X. The surface displacement also results in ESI = X. The eastern primary surface ruptures are classified as ESI = VIII if we use the length of the main scarp only (600 m). The additional ~4 km of scarps reported by Reference [21] would lead to ESI = IX. The maximum surface displacement of 3.6–4.2 m results in ESI = XI, since 3 m offset is recommended as the threshold between ESI = X and ESI = XI. (ii) The second option is to add up the primary surface rupture lengths of the western and the eastern sites and to ignore the spatial separation. Using the 7 km primary rupture length in the west and the 600 m length in the east gives ESI = IX. Using the 12 km long primary rupture zone in the west leads to ESI = X, regardless whether the additional 4 km in Reference [21] from the east are taken into account or not. Similarly, using the additional 4 km from Reference [21] leads to ESI = X, using either the 7 km primary rupture length or the 12 km primary rupture zone in the west. (iii) The third option is to measure a total primary rupture zone length from the western tip of the western primary surface ruptures to the eastern tip of the eastern primary surface ruptures. This leads to a length of 35 km without and 40 km with the additional 4 km of primary surface ruptures reported by Reference [21], corresponding to ESI = X.

#### *3.2. Secondary Ruptures*

Extensive ground cracks were caused by the Suusamyr Earthquake throughout the Suusamyr Basin and along with the surrounding mountain ranges (Figures 2 and 3). These effects were mapped in detail after the earthquake and published in References [6,21,22]. Reference [21] report that secondary ruptures covered an area of more than 4000 km2. Reference [6], however, estimates that the total area affected by secondary earthquake effects (including secondary ruptures, mass movements, and mud eruptions) is only ~2500 km2, based on the report in Reference [25]. We used the map from Reference [22] and calculated a total secondary rupture length of 114 km, covering an area of 2520 km2.

Most secondary cracks occurred in the E-W running Suusamyr River valley between the western and the eastern primary ruptures, but secondary cracks were also reported from SE of the Aramsu Range and the area between Toluk, Sarysogat, and the Belaldy rock avalanche [21,22] (Figure 2). Along the flanks of the mountain ranges, seismic shaking led to gravitational cracks that mostly paralleled the valley. These reached individual lengths of several tens of meters, spanning fissure zones of several kilometers length. The Chet Korumdy ridge was particularly affected by secondary cracks, which mainly formed as E-W elongated grabens on top of the ridge as a result of extension (Figures 3 and 4) and S-facing scarps on its southern side. We surveyed parts of the secondary ruptures on the Chet Korumdy ridge with a drone and created a high-resolution DEM, shown in Figure 4.

**Figure 4.** Secondary ruptures at Chet Korumdy ridge included small grabens that formed at the top of the ridge and mainly south-facing scarps on its southern side. (**a**) Overview of the western part of Chet Korumdy ridge. (**b**)Detail of the small grabens and secondary cracks that formed on top of the ridge. Also see Figure 3d. (**c**) Secondary cracks on top of the ridge and south-facing crown of a landslide that formed in the 1992 earthquake. The hillshade is from drone survey imagery and has a resolution of 2.91 cm/pixel.

Secondary ruptures were—to a lower extent—also reported from some flat areas in the Suusamyr basin and west of Chet Korumdy ridge. Most of them occurred close to rivers and are likely related to the presence of unconsolidated, saturated sediments subject to lateral spread.

We assign intensities of ESI = X to the Chet Korumdy area and the eastern primary surface ruptures, based on the clustering and the dimensions of the secondary ruptures. Near the western primary surface ruptures, we assign ESI = IX also because of the number, the size, and the extent of the phenomena. Intensities of ESI = VIII were assigned to the remaining, more isolated occurrences of secondary cracks.

#### *3.3. Mass Movements*

The 1992 Suusamyr Earthquake caused widespread mass movements, including landslides, rock falls, and rock avalanches. These mass movements have been intensively studied after the earthquake [6,19,21,22,26–31], and we compiled the intensity assignment from these studies.

The largest mass movement was the Belaldy rock avalanche (Figures 2 and 3) with an estimated volume of 107–108 m3 [6,22,30]. Additional large landslides were reported from the downstream section of the Belaldy river, with volumes of a few thousand to a few tens of thousand cubic kilometers [21]. Those deposits were re-activated in the year after the earthquake and caused a catastrophic mud flow that affected several villages [22]. At the southern flank of Chet Korumdy ridge, landslides with up to 0.5–1 <sup>×</sup> 106 m3 volume were mobilized and damaged the Bishkek-Osh highway [26]. A large cluster of

mass movements occurred in the Toluk-Sarysogat area, where also intense secondary ruptures were observed (Figure 2). Based on the maps in References [6,19,21,22], we compute a total affected area of 2336 km2.

The volume of the largest mass movement corresponds to intensities of ESI = X-IX, smaller slides are associated with ESI = VIII-IX. Little information is available on the volume of individual mass movements apart from the examples discussed above. For ESI2007 intensity assessments, the total affected area is also an input parameter. The southern cluster of landslides in the Toluk-Sarysogat area covers ~500 km2, which is much higher than the threshold of ~100 km2 recommended for ESI = VIII zoning. We, therefore, assigned intensities of ESI = IX to areas where mass movements cluster and intensities of ESI = VIII to more isolated mass movements. The isolated mass movements in the far west and the far east of the study area were assigned intensity ESI = VII. Using the total affected area (2636 km3) as an input value leads to a maximum intensity of ESI = X.

#### *3.4. Mud Eruptions*

No liquefaction was reported from the Suusamyr Earthquake, bud mud eruptions in the days following the earthquake were described from four locations [19,21]. Reference [21] report that the explosive ejections occurred at 'boggy sites on young active faults' and that they were even accompanied by gas combustion. Reference [6] cites eye witnesses who heard 'loud noise' and saw 'smoke that bore a pungent aroma'. They also report radially scattered rocks and mention the connection underlying active faults. Reference [21] mention that one of the eruptions resulted in a 400 m long and 100 m wide mud flow. We visited two of the sites in 2015 (Figure 3) and found muddy bogs related to landslides. At the central site near the western primary surface ruptures a muddy layer was found at a landslide toe (Figure 3g). We excavated a 25 cm deep pit into the deposits and found them to be dark, uniform, saturated soils without any hints of burned layers. We did not find any radially scattered rocks or similar signs of an explosive event. Instead, the site appeared like an over-pressured landslide mass from which water was seeping. At the westernmost site near the highway (Figure 3h,i), we encountered swampy areas on top of a landslide deposit. Puddles of several meters length contained fissures with seeping water (Figure 3i). No hints for explosive eruptions could be identified in 2015/16. If those mud eruptions happened as reported by the eyewitnesses, their traces were entirely gone 23 years after the earthquake and the sites looked like ordinary, yet very wet landslides.

These effects are hard to evaluate with the ESI2007 scale, since they are rather unusual. We decided to assign them intensity ESI = IX to account for the fact that they are unusual phenomena related to the seismic shaking, and that they had a considerable impact on the natural environment. However, the observation that they occurred after the earthquake and not co-seismically complicates the assessment. With the ESI = IX assignment, they fall well within the range of values of the neighboring EEEs.

#### *3.5. Jumping Rocks*

At several sites, jumping rocks were found during the post-earthquake surveys [21]. These are boulders of considerable size sitting on flat or nearly flat surfaces that were displaced as a result of the strong ground motion, not to be confused with rock falls or other types of gravitational effects. These effects require >1 g peak ground acceleration. Reference [21] describe rocks with a weight of 40–70 kg that moved 2–5 m on a flat surface west of the Aramsu range and an ~8 ton boulder that moved 30 cm uphill at the same site (Figure 2). Near the eastern primary surface ruptures, a ~1.3 ton block reportedly moved 2 m and two boulders of ~10 tons sitting on flat ground underwent displacement of 25–30 cm and 10◦ rotation. In the same area, an undisclosed number of rocks were moved laterally 85 cm [21].

We visited the Chet Korumdy ridge in 2015/16 and looked for the jumping rocks, but we failed to identify the boulders that moved in 1992. Most likely the footprints of the pre-earthquake position were already completely altered in the 23 years after the earthquake, so that there is no way to distinguish between a boulder that had moved and ones that did not jump.

The ESI2007 guidelines assign intensities of ESI = IX-XII to the phenomenon of jumping rocks, depending on their size, roundness, travel distance, and the local slope. Based on the values reported by Reference [21], we assign ESI = X to all the sites with large jumping boulders.

#### *3.6. Summary of the ESI2007 Assignment*

We used 132 data points to assign ESI2007 intensities (Figure 5; Table 1; see Supplementary Materials for a list with the details of each point).

**Figure 5.** Map of ESI2007 intensities compiled in this study. The location of the main shock and the strongest well-located aftershock are shown. See the Supplementary Data for details on each data point.

ESI2007 intensity ESI = XI was assigned to two sites only: To the eastern primary surface ruptures because of the amount of offset, and to the Belaldy rock avalanche because of the large volume of this mass movement.

Intensities of ESI = X cover an area of 317 km2 and are mainly based on the evaluation of primary surface ruptures, clusters of secondary ruptures, jumping rocks, and large mass movements. This isoseismal has an E-W elongated shape paralleling the fault strike and includes the Suusamyr river valley and the slopes south of it.

ESI = IX intensities cover an area of 1735 km<sup>2</sup> and include secondary ruptures, clusters of mass movements, and mud eruption sites. This isoseismal has an E-W elongated, rectangular shape.

Intensities of ESI = VIII were found in an E-W elongated ellipse with an area of 2700 km2. This intensity is based on rather isolated secondary ruptures and individual mass movements.

For intensity ESI = VII, only a few data points were available in the far east and the far west of the study area. These intensity values are based on four data points only, two mass movements and two secondary cracks. Thus, they come with high uncertainty, and the ESI = VII isoseismal could only be roughly estimated to cover ~6500 km2. This value is highly uncertain. In contrast, drawing a concave hull around all secondary effects combined (cracks, mud eruptions, landslides, rock falls, jumping stones) results in an area of 3254 km2 only, which is little larger than the ESI = VIII isoseismal.

#### **4. Discussion**

We combined our intensity assessment with the macroseismic effects on the environment reported by Reference [6], who list 41 individual sites with macroseismic observations. We followed the procedures described in detail in References [8–10].

#### *4.1. The ESI2007 Map of the Suusamyr Earthquake and Comparison with Other Intensity Scales*

The assignment of ESI2007 values for the primary surface ruptures is not straight forward as pointed out in Section 3.1. Depending on how the spatial separation between the two sets of primary surface ruptures is treated, one option is to plot two areas with different intensities around either location (ESI = X in the west and up to ESI = XI in the east). The other option is to average the observations and to plot one or two zones with ESI = X. Taking into account the unusual scarp length-to-height ratio of the eastern primary surface ruptures, for which the underlying mechanism is not yet well understood, we decide to treat them as an exception that does not justify assigning ESI = XI. Based on the scarp length, there we assign ESI = X.

Secondary effects allow drawing a more coherent image of the intensity distribution. Although the Belaldy rock avalanche stands out, due to its huge volume that would justify intensity ESI = XI, we treat it as an isolated event that fits more into the general pattern of ESI = X mass movements in the epicentral area. In general, we found that the highest ESI intensities can be assigned, due to the size of individual effects and their spatial distribution. Intensity values of ESI = IX and ESI = VIII, on the contrary, can be distinguished based on the spatial clustering of secondary effects only. This is mainly due to the lack of detailed descriptions for each locality, such as landslide volume or length of individual secondary ruptures and cracks.

The highest ESI2007 intensity X covers an E-W elongated area that roughly includes MM = X and MSK = X intensities from previous studies (Figures 5 and 6). In contrast to the ESI2007 intensities, MM and MSK intensities were assigned to two isolated patches surrounding the eastern and western primary surface ruptures, respectively, and to the Belaldy rock avalanche area (Figure 6). Intensity MSK = IX was also attributed to two isolated patches by Reference [19], whereas MM = IX was assigned to an E-W elongated patch [6]. Our intensity assessment results in a larger patch of ESI = IX, which reaches further south than MM and MSK intensities. This is due to a large cluster of mass movements in the Sarysogat-Toluk area and probably reflects the fact that this area is sparsely populated, but sits on the hanging wall above the actual fault plane. The mud eruption sites were also assigned ESI = IX, although they occurred after the earthquake [21]. However, they are related to the Suusamyr event. In any case, they do not have a strong influence on the overall ESI2007 map because they do fall in areas where other effects were mapped with ESI = IX, too. Further examples of similar phenomena are needed to understand these effects and their potential for ESI assignments better. Intensity ESI = VIII was found in an area smaller than that assigned MM = VIII and MSK = VIII, and the same is true for intensity ESI = VII. This can be explained by the fact that minor shaking can still cause damage to houses and infrastructure, while it is not sufficient anymore to cause significant earthquake environmental effects, such as landslides or secondary ruptures.

The observations on the ESI2007 intensity distribution and the comparison with MM and MSK intensity assessments highlight three peculiarities of the ESI2007 scale. First, it allows distinguishing high intensities in places were traditional scales saturate because of the total collapse of human-made structures that are not reinforced. Second, it allows assigning intensities in sparsely populated places, where traditional scales cannot be applied, due to the lack of observations. Third, it rather captures high intensities (>VII), which cause significant environmental effects. Thus, for modern earthquakes, it is best applied in combination with the traditional scales to get the full picture of earthquake effects. For paleo-earthquakes, however, the insensitivity of small to moderate events is useful because it makes it less likely to confuse moderate earthquakes with significant events that release a large seismic moment.

**Figure 6.** Comparison of ESI2007 intensities with MM isoseismals [6] and MSK intensities from Reference [19]. Note the shift to the south of the ESI = IX, which we attribute to the sparse population in the high mountains of the Aramsu range, which hampers traditional MM and MSK intensity assessments.

#### *4.2. Treating Uncertainties in the Application of the ESI2007 Scale*

Uncertainties in intensity assignment mainly arise from the unknown volumes of most landslides. However, the maximum ESI2007 intensity is well-constrained by the primary surface ruptures and the Belaldy rock avalanche. No mass movements of similar size have been reported, so maximum intensities are unlikely to be higher. The lowest intensities are defined by the extent of documented environmental effects, which is also well-known. Therefore, the shape of the isoseismals between maximum and minimum may somewhat vary, but the overall picture and the maximum intensity are unlikely to be affected by these uncertainties.

Another source of uncertainty is the unknown height of the pre-1992 scarps at the western primary surface ruptures that [7] identified on satellite imagery from 1968. The post-earthquake surveys report a height of 0.9–1.8 m [6,21], while the survey in 2015/16 found heights of 0.7–1.2 m with some outliers reaching 2 m height [7]. The imagery of the pre-1992 scarps does not reveal a free-face (Figure 5 in Reference [7]), which indicates that the height of the fresh ruptures reported by References [6,21] should be reliable because the height of the free-face is easy to measure. These heights would correspond to ESI = X. If the pre-1992 vertical offset was larger than 0.5–1.4 m, a lower intensity of ESI = IX would have to be assigned [8,10]. Given the state of preservation of the free-face in 2015/16, we consider this very unlikely (Figure 3b, Reference [7]). If the length of the western primary surface ruptures is taken into account, ESI2007 intensities do not change of course. Therefore, the height of the 1992 western primary surface ruptures does not have a large effect on the overall ESI2007 assignment. While the problem of pre-existing scarps is usually a source of uncertainty hard to deal with in paleoseismology studies, this effect does not pose a significant problem in the case of the 1992 Suusamyr earthquake. Extending the use of scarp-diffusion modelling to multi-event scarps [32–34] may help to improve the interpretation of multi-event scarps.

From a paleoseismology perspective, the peculiar offset-to-length ratio of the eastern primary surface ruptures is interesting. Although the nature of the eastern primary surface ruptures of the 1992 Suusamyr Earthquake is well understood [6,7,21,23–25], such unusual primary effects are likely hard to

interpret when attributed to a paleo-earthquake. One obvious interpretation of such an unusal feature would be that a longer primary scarp once existed, and that erosion has eradicated the prolongation of the scarp along strike. In that case, using the vertical offset only would result in ESI = X-XI. Using the mapped length to assign ESI2007 intensities would result in ESI = VIII only. In that case, the earthquake intensity, and thus, also the magnitude, would be dramatically underestimated.

If in case of a paleo-earthquake, the eastern and western primary ruptures would be assigned to the same event depends on the precision of the dating methods and the interpretation of the individual researcher. Reference [7] have shown in the Suusamyr case, this may indeed be challenging. Another possibility would be that, due to the unusual offset-to-length ratio of the eastern primary surface ruptures, they are not even identified as such, but attributed to other effects, such as being a landslide toe or being a terrace riser. Paleoseismological trenching could help to solve this problem. However, in the case of false attribution, there would be no intensity ESI = XI in this area. The clustering of landslides and secondary cracks could probably be ascribed to topography and preservation effects, and the epicenter of the earthquake would be located further to the west near the western surface ruptures.

#### *4.3. The Problem of Strong Aftershocks*

A significant source of uncertainty in the ESI2007 assignment is related to the three strong aftershocks that occurred within two hours of the main event. These earthquakes had magnitudes of MS6.6, MS6.6, and mb6.0, and they occurred west of the Aramsu Range [18], although their mechanism and location are not well constrained. Given their magnitudes, they may have even caused the western primary surface ruptures [7] or at least contributed to their total lengths. The same holds true for the secondary effects. There is virtually no possibility to distinguish between the effects of the main shock and those of the three strong aftershocks. The earthquakes occurred so close in time, and in such a sparsely populated region, that there are no eyewitness accounts that could help to solve this problem. Perhaps the shape of the isoseismals with the highest intensities could support the interpretation of a combined effect of the main shock and the aftershocks, but we will most likely never know.

This problem usually arises in paleoseismology, when due to the limits of Quaternary dating methods, it is impossible to distinguish between one strong earthquake and two smaller ones, or similar combinations. However, the Suusamyr case illustrates that even for modern events this may be an issue, especially in sparsely populated regions and when significant aftershocks occur right after the main shock before any documentation of the earthquake effects has been undertaken.

#### *4.4. Lessons for ESI2007 Assignments*

The peculiar primary surface rupture pattern of the Suusamyr Earthquake poses a problem for ESI2007 assignments as discussed above. Thanks to the widespread secondary effects, we were able to produce an intensity map that relies on a large number of observations. If only information on the primary surface ruptures were available, for example from a similar paleo-earthquake, the picture would look different. This is important from a paleoseismological point of view. For many, if not most pre-historic earthquakes, only data on the primary surface ruptures are derived from field mapping and/or trenching studies. Reference [7] have shown that this may be misleading in the Suusamyr Basin. Another example is the MW6.6 1998 Fandoqa Earthquake that occurred on the Gowk Fault in Iran [35]. This earthquake re-ruptured a section of the fault that had already ruptured in an MW7.1 earthquake in 1981. Primary surface ruptures were observed in both cases, but the offsets caused by the smaller earthquake were several times larger than those of the stronger event. Reference [35] report that the two earthquakes ruptured the fault at different depths and rightly emphasized the problem this causes for paleoseismological studies in similar settings.

Similarly, complex primary surface ruptures from paleo-earthquakes may not be entirely understood and, therefore, lead to a wrong ESI2007 assignment [36]. Modern examples of such complex primary surface ruptures are among others the 1889 Chilik Earthquake [37], the 1911 Chon Kemin Earthquake [38], the 2010 El Mayor Cucapah Earthquake [39], and the 2016 Kaikoura Earthquake [5].

Modern earthquake ruptures can be studied extensively in the field [40], using seismological methods, and with (space) geodesy [41]. Their secondary effects can also be documented in great detail. Most studies of paleo-earthquakes, however, mainly focus on the primary surface ruptures, because this is often the only definitively seismogenic evidence left in the landscape. Notwithstanding, with the Suusamyr example in mind, we emphasize the importance of collecting as much information on the environmental effects of such earthquakes as possible.

Unfortunately, secondary effects may also be problematic for intensity assessment. The Belaldy rock avalanche formed a landslide lake which burst ten months after the Suusamyr Earthquake in June 1993. Triggered by snowmelt, a debris flow mainly incorporating material from the 1992 landslides reached and damaged the village of Torkent 30 km downstream [22,42]. Such phenomena may be easily confused with coseismic effects. In case of the Suusamyr Earthquake, the mud eruption sites that we visited a few years after the earthquake were virtually indistinguishable from ordinary landslides. Not only will they be interpreted as such in the case of a paleo-earthquake, but mud eruptions may also happen in any future landslide without a seismic trigger. From a paleoseismological perspective, temporal clustering of earthquake environmental effects can help to identify seismic events, but as can be seen from the 1998 Fandoqa Earthquake [35] or the Central Italy earthquake series [43], earthquakes can occur very close in time on the same fault, which is beyond the resolution limit of any Quaternary dating method (perhaps excluding varves).

#### **5. Conclusions**

We have compiled an ESI2007 intensity map of the 1992 Suusamyr Earthquake based on earthquake environmental effects only. This assignment was based on published data on primary and secondary earthquake effects, and in our field mapping. The ESI2007 intensity distribution is roughly comparable to the MM and MSK intensities reported, although we note important differences mainly in the shape of the isoseismals (Figure 6), but also in the area they encompass (Table 1). One major advantage of the ESI2007 scale is that it allows incorporating data from uninhabited areas, which led to a shift to the south of the overall intensity pattern.

Due to the peculiar pattern of 1992 primary surface ruptures, the assignment of ESI2007 intensities to these effects is not straight forward and allows different approaches. Depending on the input data used, the resulting ESI2007 intensities vary significantly. When combined with the secondary effects of the earthquake, however, a coherent picture emerges.

We show that this may not be the case when dealing with a paleo-earthquake of similar characteristics. The unusual offset-to-length ratio of the eastern primary surface ruptures and the gap between the western and eastern primary surface ruptures are prone to false or simplified interpretations which may hamper the correct assignment of ESI2007 intensities. This is especially important as the secondary effects, such as gravitational cracks and mass movements are likely to be short-living features that could also happen without a seismic trigger. The eastern primary surface ruptures occurred in a river bed and are also likely to be eroded within a relatively short time span. Thus, a paleo-earthquake that resembles the 1992 event will likely be described incompletely with the ESI2007 scale. This is an important lesson for the application of this scale to paleo-earthquakes in similar environments, and to events with a complex or anomalous primary surface rupture pattern.

**Supplementary Materials:** The following are available online at http://www.mdpi.com/2076-3263/9/6/271/s1.

**Author Contributions:** Conceptualization, C.G.; methodology, C.G.; validation, C.G., R.W., E.A., A.E. and K.A.; formal analysis, C.G.; investigation, C.G., R.W., E.A., A.E. and K.A.; resources, C.G., R.W., E.A., A.E. and K.A.; data curation, C.G., E.A. and A.E.; writing—original draft preparation, C.G.; writing—review and editing, E.A.; project administration, R.W.; funding acquisition, R.W.

**Funding:** This study was funded by the Earthquakes without Frontiers project, NERC and ESRC grant code: EwF\_NE/J02001X/1\_1), by the Center for Observation and Modelling of Earthquakes, Volcanoes and Tectonics (COMET, GA/13/M/031), and by Looking inside the Continents from Space (LiCS) (NE/K011006/1).

**Acknowledgments:** We thank three anonymous reviewers for their constructive comments that helped to improve the manuscript. We thank the editor Alicia Wang for handling our paper.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Earthquake-Induced Landslide Risk Assessment: An Example from Sakhalin Island, Russia**

#### **Alexey Konovalov \*, Yuriy Gensiorovskiy, Valentina Lobkina, Alexandra Muzychenko, Yuliya Stepnova, Leonid Muzychenko, Andrey Stepnov and Mikhail Mikhalyov**

Sakhalin Department of Far East Geological Institute, Far Eastern Branch, Russian Academy of Sciences, Yuzhno-Sakhalinsk 693023, Russia

**\*** Correspondence: a.konovalov@geophystech.ru

Received: 30 April 2019; Accepted: 8 July 2019; Published: 11 July 2019

**Abstract:** Damages caused by earthquake-induced ground effects can be of the order or significantly exceed the expected damages from ground shaking. A new probabilistic technique is considered in this study for earthquake-induced landslide risk assessment. A fully probabilistic technique suggests a multi-stage hazard assessment. These stages include the determination of seismic hazard curves and landslide probabilistic models, a vulnerability assessment, and geotechnical investigations. At each of the stages, the uncertainties should be carefully analyzed. A logic tree technique, which handles all available models and parameters, was used in the study. The method was applied considering child education facilities located at the foot of a natural slope in the south of Sakhalin Island which is known as an active seismic and land sliding area. The significant differences in the ground motion scenario in terms of the 475-year seismic hazard map and the fully probabilistic approach considered suggests that seismic landslide risk could be underestimated or overestimated when using the 475-year seismic hazard map for risk assessment. The given approach follows the rational risk management idea that handles well all possible ground motion scenarios, slope models, and parameters. The authors suggest that the given approach can improve geotechnical studies of slope stability.

**Keywords:** earthquake-induced landslide; fully probabilistic technique; Newmark's method; Sakhalin Island; risk

#### **1. Introduction**

Large earthquakes affecting urban areas are one of the most destructive natural hazards and can lead to significant impacts on the built and human environment. Generally, earthquake loss models consider ground shaking and ground failure (such as landslides, liquefaction, and faulting) hazards. Damages caused by earthquake-induced ground effects, in some cases, significantly exceed the damages from direct ground shaking [1,2]. Damages related to seismically-induced landslides can be considerable due to the full collapse or loss in functionality of facilities, roads, pipelines, and other lifelines [3–7].

There are numerous causative factors for seismically-induced gravity-slope processes on Sakhalin Island, which is recognized as an area with a high level of geohazards. A total of 70% of the South Sakhalin territory is susceptible to landslide activity [8]. According to general seismic hazard maps, Sakhalin Island is a seismically active area, with an 8–9 MSK-64 macroseismic intensity for the 475-year return period [9].

As a recent example of earthquake-induced landslides, the Mw = 6.2 2 August 2007 Nevelsk earthquake should be noted. The Nevelsk earthquake was followed by aftershock sequences with a relatively high productivity level [10]. Focal mechanisms of the mainshock and aftershocks indicates the west dipping (38–40◦) fault planes.

The largest aftershock with a magnitude of Mw = 5.9 followed about two hours after the mainshock origin time. They generated tsunami waves up to 3.2 m high. The ground shaking effects in Nevelsk caused by the mainshock and largest aftershock corresponded to a 7–8 MSK-64 macroseismic intensity (Figure 1). Massive damages of buildings (>200), bridges, railway and roads were found during macroseismic inspection [10].

The 2007 Nevelsk earthquakes caused massive release of methane from the coal beds in the coastal zone of 40 km in length and an uplift of benches in the area of the Nevelsk sea port [10].

**Figure 1.** Contour map of peak ground acceleration (% *g*) for the 2 August 2007 Nevelsk earthquake (Mw = 6.2). Filled boxes indicate the settlements with MSK-64 felt reports. The mainshock (Mw = 6.2) is shown by the red filled circle, the largest aftershock (Mw = 5.9) by the blue filled circle.

As a result of significant ground shaking, subsidence cracks and shallow landslides (up to 200–300 m3) were widely recorded (Figure 2) within the Nevelsk urban area (16–21 km from the epicenter). The 2007 Nevelsk earthquake occurred in a relatively dry period. There was a recorded 69 mm of rain precipitation, representing 45% of the mean annual value, for the two months before the mainshock [11]. Therefore, seismically-induced landslides remain a major natural hazard on Sakhalin Island that should be considered in the risk assessment strategy.

**Figure 2.** An example of shallow landslides caused by the 2007 Nevelsk earthquake.

Most hazard assessment techniques are designed for producing susceptibility or likelihood maps on large scales (regional or global) [2,12,13]. Because material parameters are difficult to identify in detail for large areas, slope parameters are estimated from topographic, geologic, and other geospatial information. These maps are most commonly used for an estimation of the general hazard level and specifying sites where detailed geotechnical investigations are needed.

The aim of this study is to apply the earthquake-induced landslide risk assessment technique at a local scale. This paper proposes a fully probabilistic approach that handles all available ground motion scenarios and geomechanical slope models well. The important issue is the estimation of the uncertainties from the spatial and temporal variability of soil parameters.

The child education facilities (Figure 3) located under the natural slope in Nevelsk (Sakhalin Island, Russia) were considered in this study as an example of the given approach. This site was proposed for application of the methodology due to available slope material parameters previously obtained from geotechnical investigations.

**Figure 3.** Map of the slope fragment with an expected landslide boundary.

#### **2. Materials and Methods**

#### *2.1. Fully Probabilistic Risk Assessment Technique*

The seismic hazard maps in terms of the 475-year ground shaking intensity are most commonly used as a triggering condition for analyzing the slope stability under seismic loading [2,13,14]. Uncertainties about the 'best' ground motion scenario for land sliding prediction are still being discussed [15].

A fully probabilistic risk assessment technique that handles all available ground motion scenarios well is considered in this study [16]. The given approach follows a rational risk management idea that deals with all possible models, scenarios, and uncertainties.

The total risk of an earthquake-induced landslide hazard to the facility located at the foot of the slope in the next *T* years can be expressed in terms of probability and vulnerability. Formally, the probability of the occurrence of *N* seismically-induced landslides within a certain time interval should be considered. In practice, the probability of the triggering of two or more slope failures within 50 years is extremally small.

A fully probabilistic approach suggests determination of the occurrence probabilities of various ground shaking levels and probabilities of landslide triggering for these ground motion scenarios. That gives the total probability of earthquake-induced slope failure in a certain time interval considering all possible ground motion scenarios:

$$P\_{sf}(T) = \sum\_{i} \sum\_{i} w\_{j} P(\text{PGA} = a\_{i}|T) \cdot Pr\_{j}(a\_{i}) = \sum\_{j} \sum\_{i} w\_{j} p\_{i j \prime} \tag{1}$$

where *P*(PGA = *ai*|*T*) is the probability of the occurrence of ground shaking level *ai* in the next *T* years (ground motion scenario), *Prj*(*ai*) is the probability that seismic loading *ai* will trigger a landslide in the frame of the geomechanical slope model *j*, *pij* is the probability of slope failure in the next *T* years under ground motion scenario *ai* in the frame of the geomechanical slope model *j*, and PGA is the peak ground acceleration at the given site.

The summation in (1) is carried out for all available geomechanical slope models ranked by weights *wj*, where:

$$\sum\_{j} w\_{j} = 1.\tag{2}$$

The resistance of the facility to prevent the land sliding can be expressed in terms of the structural damage risk (physical vulnerability). The vulnerability is defined in this study as the degree of loss of a given element at risk and ranging from 0 (no loss) to 1 (total loss).

Vulnerability depends on a number of factors [17,18], such as the sliding mass volume and its velocity, position of landslide initiation, etc. For an accurate risk assessment, all of these factors should be considered. In many case studies, the vulnerability assessment is still characterized by large uncertainties.

In a simple case, vulnerability depends on the intensity of the landslide that can hit the facility located at some distance from the landslide source [18,19]. It means that a certain geomechanical slope model can be associated with the expected intensity level of the landslide.

Then, the total risk *R*(*T*) of facility damage caused by seismic land sliding within *T* years can be written as:

$$R(T) = \sum\_{j} \sum\_{i} w\_{j} v\_{j} p\_{i j \prime} \tag{3}$$

where *vj* is the vulnerability of the facility in the frame of the geomechanical slope model *j*.

Basically, Equations (1) and (3) suggest a multi-stage and multi-hazard approach. These stages include site-specific probabilistic seismic hazard analysis, vulnerability assessment, geotechnical investigations, and landslide probability calibration.

#### *2.2. Causative Factors of Studied Area*

Seismically-induced landslide risk assessment is complicated by the numerous factors and conditions contributing to slope failure. This section discusses the major causative factors of the studied area, helping to constrain the slope models and parameters for the probabilistic hazard analysis (Figure 4).

**Figure 4.** Side view of the location of the natural slope and school considered in this study. The red line indicates the expected landslide breakout wall. The red arrow shows the path direction of the expected landslide.

#### 2.2.1. Geology

From the basis of the tectonic terrains scheme, the studied area is related to the southern part of the West Sakhalin Terrane located to the north of the Lopatino Cape. The target area consists of Late Cretaceous and Cenozoic outcropping rocks (Figure 5).

Quaternary and Neogene age deposits are widely presented in the studied area. The Quaternary soils of low-thickness are presented by deluvial-eluvial deposits. The Neogene soils of the Nevelskaya formation are presented by volcanic sandstones, siltstones, tuffs, and tuffits. A geologic map with detailed formation lithology is given in Figure 5 [20].

From the view of general engineering geology, the following layers are highlighted in the lithological slope section: 1—topsoil (0.2 m); 2—fine-gravelly with sandy filler up to 50%, dense (deluvial-eluvial) (0.2 m); 3—tuff-sandstone on clayey cement with low-strength siltstone interlayers (3.6 m).

**Figure 5.** Geological structural map of studied area compiled according to [20]: 1—Upper Miocene Kurasiiskay formation: lower subformation (kr1)—siliceous claystones and upper subformation (kr2)—alternation of sandstones and siltstones; 2—Middle Miocene Verhneduiskaya formation (vd)—alternation of sandstones and siltstones and claystones, hards coal; 3—Lower Miocene Nevenlskaya formation: lower subformation (nv1)—alternation of sandstones and siltstones and upper subformation (nv2)—alternation of sandstones and siltstones; 4—synclineaxis (a) anticlineaxis (b); 5—fault zone; 6—studied area.

#### 2.2.2. Geomorphology

Nevelsk is located on the western flanks and sea terrace's foot, which are western branches of the South Kamyshoviy range. Terraces have an absolute elevation till 200 m and are dissected by narrow, deeply incised V-shaped valleys of rivers and rills [21].

Landscape dissection depth reaches 200 m, and the slope steepness is 45◦ and more. Terrace surfaces are inclined to the sea side.

Abrasion sea terrace scarps, rivers, and rill valley slopes are well-sodded and covered by grassy-shrub vegetation, and in some places are forested. Vegetation does not prevent landslides.

River valley and sea terrace slopes are complicated by old block slides, which make slopes stepped. The child education facilities are located close to the foot of the high sea terrace (Figures 3 and 4).

#### 2.2.3. Climatic Settings

The Nevelsk district is located along the Tatar Strait of the Japanese Sea, and is strongly affected by a warm Tsushima current, so it witnesses the warmest summer and winter of the entire island. The territory is characterized by a large precipitation quantity during the second part of the summer.

The mean annual precipitation quantity in Nevelsk is 911 mm [22]. Precipitation falls inhomogenously. The major precipitation period is related to the second half of the summer and early autumn, and the minimum period is the second half of winter.

Two landslide activation periods related to soil moisture conditions are typical for the studied area.

The first period (end of May to early June) is related to melting snow which leads to rapid soil humidification and often, to shallow (1 m deep) landslides. Rainfall can trigger landslides during this period.

The second period (August to October) is related to high cyclone activity, leading to a large rain precipitation quantity.

The long-term annual average precipitation quantity during the warm period is 579 mm. The rainiest month is September, with a precipitation quantity of ~111 mm. The precipitation quantity during a cyclone event (in August–September) can exceed the monthly average level. The maximum precipitation quantity (211 mm) was registered during the Phillis cyclone event (2–7 August 1981).

#### 2.2.4. Soil Moisture Conditions

The soil moisture conditions of a potentially sliding mass represent one of the most important factors for landslide prediction.

Three soil moisture models are assumed. The models are constrained by generalizing the climatic settings in the studied area. These simplified models are characterized by the following properties:


#### *2.3. Materials*

#### 2.3.1. Ground Motion Scenarios

The main goal of probabilistic seismic hazard analysis (PSHA) is to determine the probability of exceeding a certain ground shaking level at a given site within the time interval of interest [23]. The result of such analysis is a seismic hazard curve, which demonstrates the relationship between the probability of exceeding *P*(PGA > *a*|*T*) and ground shaking level *a*. Under ground shaking conditions, we consider PGA.

The occurrence probability *P*(PGA = *ai*|*T*) is recognized in this study as a ground motion scenario. The transition from the probability of exceeding to the discrete probability is given as:

$$P(\text{PGA} = a\_i | T) = P(\text{PGA} > a\_i | T) - P(\text{PGA} > a\_{i+1} | T). \tag{4}$$

Therefore, Equation (4) defines the probability of the occurrence of a certain ground shaking level in the next *T* years. It was substituted into Equations (1) and (3).

The seismic hazard curve for a given site was computed using CRISIS 2015 software [24]. The authors used the regional seismic source models and ground motion prediction equations tested in the previous PSHA studies [25].

The average 30-m ground layer shear-wave velocity VS30 at the site is of the order of 300 m/s. It was used for site correction within the PSHA stage considering the site-correction term in the ground motion prediction equations.

The seismic hazard curve used in this study is shown in Figure 6. The 10% exceeding probability (or 475-year return period probability) corresponds to the PGA value of 0.38 *g* (Figure 6). Peak ground acceleration of 0.38 *g* corresponds to a ground shaking intensity of 9 MSK64. It is in good agreement with the estimates based on the maps of general seismic zonation [9].

**Figure 6.** The probability of exceeding the given peak ground acceleration. Red dotted line indicates the 10% probability level.

#### 2.3.2. Probability of Landslide Triggering

Jibson et al. [12] calibrated the slope failure probability by analyzing landslide field data after the 1994 Northridge earthquake and predicted Newmark displacement (*DN*). The authors used the Newmark sliding block model [26] to quantify the probability of landslide triggering given the PGA value (*a*) and critical slope acceleration (*ac*). The Jibson probability model follows the Weibull distribution:

$$Pr(D\_N) = 0.335 \left[ 1 - \exp\left(-0.048 D\_N^{1.565}\right) \right] \tag{5}$$

where

$$
\log D\_N(a) = -0.215 + \log \left[ \left( 1 - \frac{a\_c}{a} \right)^{2.341} \left( \frac{a\_c}{a} \right)^{-1.438} \right] \pm 0.51. \tag{6}
$$

The critical slope acceleration *ac* of a sliding mass is considered as a simple function of the static factor of safety *FS* and the slope geometry:

$$a\_{\mathbb{C}} = (F\_{\mathbb{S}} - 1)g \sin \alpha\_{\prime} \tag{7}$$

where *g* is the gravity factor and α is the dip angle of the potential sliding surface.

According to the limit equilibrium theory, the static factor of safety *FS* is defined as the relationship between the force keeping the sliding mass on the slope and the force moving the sliding mass down the slope [12]:

$$F\_S = \frac{c'}{\gamma z \sin \alpha} + \frac{\tan \phi'}{\tan \alpha} - \frac{m \gamma\_w \tan \phi'}{\gamma \tan \alpha},\tag{8}$$

where *c* is the effective cohesion, *z* is the slope-normal thickness of the potential sliding mass, γ is the material unit weight, γ*<sup>w</sup>* is the unit weight of ground water, ϕ is the effective friction angle, and *m* represents the fractional depth of the water table with respect to the total slide depth. The sliding mass is stable if *FS* > 1, and unstable when *FS* < 1.

Soil moisture conditions in the studied area and southern California are significantly different. Virtually no rain had fallen prior to the 1994 Northridge earthquake, which was used for calibrating the seismically-induced landslide probability [12]. The pore-water pressure term was dropped from Equation (8) in the Jibson probability model. The soils in the studied area are expected to be saturated most of the year, so the authors paid great attention to the third term of Equation (8). Several models with a non-zero pore-water pressure term are hypothesized in this study.

Therefore, the Jibson probabilistic model for California (5) was imported into the Equations (1) and (3).

#### 2.3.3. Geomechanical Slope Models and Logic Tree

The physico-mechanical parameters of the slope soils were obtained during geotechnical investigations at the site (see Acknowledgments Section). Geotechnical studies included sampling undisturbed cores and laboratory tests. The material parameters of the soil mass are given in Table 1.


**Table 1.** Physico-mechanical parameters of the soil mass.

Identifying the slab thickness *z* of potential sliding mass requires high-quality seismic profile, geological, and geophysical data. These can give additional information about the layer structure that helps to fix the slope-normal thickness or to identify the preexisting landslides. Such data is not available for the studied area. Therefore, the slab thickness variability was hypothesized through the logic tree approach, which is commonly used in PSHA. Four models with varied slab thicknesses (1, 2, 4, and 8 m) were considered. The corresponding weights are equal to 0.3, 0.3, 0.3, and 0.1 (Figure 7). The weights were determined by an expert's view.

The given difference in weighting reflects our view that the landslide initiation depth of an order of 8 m seems to be unrealistic. However, the authors are aware of the uncertainties connected with this choice in weighting scheme.

The accurate weighting of models with slab thickness uncertainties can be realized through the corresponding probability density distribution. The probability that the landslide has a certain slip surface depth value is defined by its statistical distribution. A significant amount of observations in the studied area are required for developing such a probability model.

The second variable parameter is the proportion of the sliding mass thickness that is saturated (*m* in Equation (8)). The water saturation parameter was defined as a simplified relation given by soil moisture conditions (see Section 2.2.4):

$$m = \frac{\text{depth of water inversion}}{\text{slab} - \text{normal thickness}}.\tag{9}$$

The proportion between the period duration of water invasion into the corresponding soil depth and the annual period was used for the weighting. For instance, for dry soil conditions, the corresponding weight was defined as 5/12~0.4 (Figure 7).


**Figure 7.** Logic tree for handling the epistemic uncertainties in the considered models and parameters. The weights of the branches are given in circles.

A full logic tree that handles uncertainties in the considered models and parameters is shown in Figure 6. The resulting critical slope acceleration for each branch was calculated according to (7). For the unstable case (*FS* < 1), the critical slope acceleration was defined as the minimum acceleration in the given seismic hazard curve (0.005 *g*).

#### 2.3.4. Vulnerability Assessment

The physical vulnerability term is widely used in many scientific and engineering fields [27]. It defines the probability that a given element at risk will be damaged under a certain impact. Recent seismic risk analyses have usually dealt with seismic fragility curves [28]. The fragility curve defines the conditional probability of exceeding the certain damage state of the building for a given level of ground shaking intensity. Increasing the research volume of earthquake damage data allows researchers to calibrate seismic fragility curves in many seismically-active regions.

Unfortunately, accurate calibrating of the fragility curves is still complicated by limited damage data for most types of landslides. Moreover, different physical mechanisms are associated with different types of landslides (rockfalls, block slides, debris flow etc.). Vulnerability depends on a number of factors [17,18], such as the sliding mass volume, landslide velocity, landslide source, etc.

A physical vulnerability assessment is characterized by uncertainties that can be either epistemic or aleatory [29]. Epistemic uncertainties can be associated with the simplification of the landslide intensity estimation, the characterization of elements at risk, the vulnerability models, expert view, etc. Aleatory uncertainties can be associated with the spatial variability of parameters [29], such as landslide

source location. For this type of uncertainty, the expected structural damage differs, depending on whether the facility is located on the foot or landslide body [19].

Russian regulation norms have a specific definition of landslide hazard, depending on its type and volume, but quantitative estimations of the structural response are not available.

Due to the reasons mentioned above, vulnerability assessment still remains somewhat subjective and is usually performed in small study areas (local scale). The main issue of the vulnerability assessment in our study is to give a suitable estimation of expected damage using different sliding mass models.

Available data that defines the physical loss of buildings for the wide range of sliding processes were used in this study [18,19]. Most available data for physical damage estimates come from debris flow studies. However, the authors are aware of the uncertainties connected with the choice of vulnerability values.

Different vulnerability values in connection to the expected intensity level of a landslide and facility type are given in [18]. Jaiswal et al. [18] used the landslide volume for intensity estimation (M-I, M-II and M-III). The M-I intensity class corresponds to shallow landslides with a volume of less than 1000 m3. If a landslide becomes bigger, with a volume ranging from 1000 m3 to 10,000 m3, the intensity follows the M-II class. For landslides with a volume greater than 10,000 m3, the intensity reaches the M-III class.

For similarity, the landslide volume was estimated as *V* = *AL*×*z*, where *AL* is the sliding area. The equation *AL* = min(706 × *z*, *ALmax*) links the normal slope thickness and sliding area [19]. According to field studies, the maximum expected sliding area *ALmax* is of the order of 2335 m2 (Figure 3). Despite the fact that the given relationships have significant uncertainties, they help to constrain the landslide intensity class in the studied area.

According to [19], the landslide intensity is associated with the slip surface depth, which is directly related to the slab thickness of the sliding mass. There are [12] five different landslide intensity classes that correspond to the given slip surface depths of 1 m, 2 m, 6 m, 10 m, and 20 m. The estimated vulnerabilities for slip surface depths of 1 m, 2 m, 6 m, and 10 m were imported to our model with corresponding slab thicknesses of 1 m, 2 m, 4 m, and 8 m.

The reinforced concrete building type (Type-4 or SBT4) was used in this study for the vulnerability assessment. This building type is closely related to the child education facilities considered by its material strength properties.

The vulnerability exhibits average, minimum, and maximum values. For dry soil conditions, we use average vulnerability values, and for water-saturated soil conditions, we use maximum vulnerability values, as we expect that a water-saturated sliding mass has a relatively higher landslide velocity.

We used the estimated vulnerability of elements at risk located within run-out paths of a landslide. Table 2 contains the vulnerability estimations for the considered facilitates and slope models in full accordance with [18,19]. The final weights for each of the considered slope models in Figure 6 are given with respect to the vulnerability values from Table 2.


**Table 2.** Simplified vulnerability model for the facilities located within landslide run-out paths.

#### **3. Results and Discussion**

Since the critical slope acceleration becomes available for each of the considered slope models, it is substituted into the Jibson probabilistic model (4). Next, the probability of occurrence of a landslide (*Pr*(*DN*(*ai*))) under the given seismic loading is multiplied by the probability of the occurrence a certain ground shaking level in the next 50 years (*P*(PGA = *ai*|*T*)). This gives the probability (*pij*) of slope failure in the next 50 years under ground motion scenario *ai* in the frame of the geomechanical slope model *j*. Then, the summation of discrete probabilities *pij* is carried out for all available ground motion scenarios *ai* and geomechanical slope models ranked by weights *wj*. Since the final weights for each of the considered slope model are calculated with respect to the vulnerability (Figure 7), the summation of *pij* multiplied by its final weights *wj* gives the total risk value.

The total risk of an earthquake-induced landslide hazard for the child education facility located at the foot of the slope within the next 50 years was computed according to Equation (3) and the weighting scheme illustrated in Figure 7. The corresponding risk value appeared to be 7.4%. It reflects the high level of seismic and landslide hazards and high vulnerability value for the facilities located at the foot of the natural slope. Therefore, the high risk value shows a considerable hazard of seismically-induced landslides in terms of civil engineering.

In order to estimate the most probable ground motion scenario that will cause damage to the facility by an earthquake-induced landslide, the discrete probabilities *pij* from Equation (3) were plotted against the corresponding ground motion level. This gives the contribution to the total risk from each slope model and each ground shaking level (Figure 8). As can be seen from Figure 8, the distributions have a modal form.

**Figure 8.** Average weighted (1) and unweighted (2–8) curves of probability of slope failure in the next 50 years as a function of peak ground acceleration: 2—*ac* = 0.695 *g*; 3—*ac* = 0.533 *g*; 4—*ac* = 0.247 *g*; 5—*ac* = 0.166 *g*; 6—*ac* = 0.086 *g*; 7—*ac* = 0.023 *g*; 8—*ac* = 0.005 *g*.

The average weighted probability density function computed according to the logic tree scheme (Figure 7) has a peak at PGA = 0.15 *g* (7–8 MSK64 intensity level). The PGA value of the order of 0.15 *g* can be recognized as the 'best' ground motion scenario for the considered slope models and parameters. In terms of PSHA, the given ground shaking level is closely related to the 100-year return period (according to Figure 6).

At the same time, Figure 7 shows the contribution to the total risk from each slope model. The saturated soil's model (*ac* = 0.005 plot in Figure 8) should be considered as the most probable scenario for the considered area. Generally, it makes the slope unstable, regardless of the triggering conditions.

The slope model corresponding to the slip surface depth of *z* = 4 m and dry soil conditions (*m* = 0) has a clear peak around PGA = 0.2–0.3 *g* (*ac* = 0.023 *g* plot in Figure 8). This ground shaking value is closely related to the 475-year probability.

It should be noted that a significant ground shaking level corresponding to a 7–8 MSK-64 intensity caused shallow landslides within the Nevelsk urban area in 2007. The PGA value corresponding to the 'best' ground motion scenario is the same order as the ground shaking intensity in the studied area generated by the 2007 Nevelsk earthquake. No facilities were damaged by the earthquake-induced landslide. This may be explained in terms of soil moisture conditions. The point is that several months prior to the Mw = 6.2 2 August 2007 Nevelsk earthquake, no rain had fallen. As a result, the 2007 Nevelsk earthquakes occurred in a relatively dry period. For such untypical conditions, the corresponding ground motion scenario is the order of PGA = 0.2–0.3 *g* or an 8–9 MSK64 intensity.

When generalizing results of the study, several assumptions and models should be mentioned.

Constants of the Jibson landslide probability model in the target area are generally not the same as those in southern California. The shape of the probability curve depends on the geological, geomorphological, and soil moisture conditions. Jibson et al. [12] argue that the variability of constants in Equation (4) should be proposed if a significant difference between the target area and southern California is found. This means that the logic tree will be complicated by uncertainties from the landslide probability model.

One more source of epistemic uncertainty comes from the spatial and seasonal variability of slope material parameters. Geotechnical studies of rocks and soils within the large area of the natural slope within different seasons should significantly reduce the uncertainties associated with the variability of geomechanical slope models (Figure 8).

#### **4. Conclusions**

A fully probabilistic technique is considered in this study for an earthquake-induced landslide risk assessment in a relatively small area. The given method suggests a multi-stage and multi-hazard approach. These stages include a site-specific probabilistic seismic hazard analysis, vulnerability assessment, geotechnical investigations, and landslide probability calibration.

As a case study, the child education facility located under the natural slope was considered in this study. The total risk of earthquake-induced landslide hazard to the child education facility located at the foot of the slope within the next 50 years was of the order of 7.4%.

A significant difference between the ground motion scenario in terms of the 475-year seismic hazard map and considered fully probabilistic approach suggests that seismic landslide risk could be underestimated or overestimated when using the 475-year seismic hazard map for landslide risk assessment. The given approach follows the rational risk management idea that handles well all possible ground motion scenarios and geomechanical slope models.

An important factor that leads to an increase of the total risk is the saturated soil mass. Geotechnical studies of rocks and soils within large areas of a natural slope within different seasons should significantly reduce the uncertainties associated with the variability of geomechanical slope models.

The aim of future research is to produce regional seismically-induced landslide hazard maps using the fully probabilistic approach.

**Author Contributions:** Conceptualization, A.K.; methodology, A.K.; software, A.S.; validation, A.S.; formal analysis, Y.G.; investigation, Y.G., V.L., A.M., Y.S., L.M., A.S. and M.M.; resources, A.K. and Y.S.; writing—review and editing, A.K. and L.M.; visualization, Y.S. and A.M.; project administration, A.K. and Y.G.

**Acknowledgments:** This study was supported through computational resources provided by the Shared Facility Center "Data Center of FEB RAS" (Khabarovsk) [30]. Material parameters of the studied area were supported by Biolit LLC (Yuzhno-Sakhalinsk).

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
