*Article* **Comparison of Earthquake-Triggered Landslide Inventories: A Case Study of the 2015 Gorkha Earthquake, Nepal**

#### **Sansar Raj Meena \* and Sepideh Tavakkoli Piralilou**

Department of Geoinformatics—Z\_GIS, University of Salzburg, 5020 Salzburg, Austria; sepideh.tavakkoli-piralilou@stud.sbg.ac.at

**\*** Correspondence: sansarraj.meena@sbg.ac.at

Received: 13 September 2019; Accepted: 9 October 2019; Published: 10 October 2019

**Abstract:** Despite landslide inventories being compiled throughout the world every year at different scales, limited efforts have been made to critically compare them using various techniques or by different investigators. Event-based landslide inventories indicate the location, distribution, and detected boundaries of landslides caused by a single event, such as an earthquake or a rainstorm. Event-based landslide inventories are essential for landslide susceptibility mapping, hazard modeling, and further management of risk mitigation. In Nepal, there were several attempts to map landslides in detail after the Gorkha earthquake. Particularly after the main event on 25 April 2015, researchers around the world mapped the landslides induced by this earthquake. In this research, we compared four of these published inventories qualitatively and quantitatively using different techniques. Two principal methodologies, namely the cartographical degree of matching and frequency area distribution (FAD), were optimized and applied to evaluate inventory maps. We also showed the impact of using satellite imagery with different spatial resolutions on the landslide inventory generation by analyzing matches and mismatches between the inventories. The results of our work give an overview of the impact of methodology selection and outline the limitations and advantages of different remote sensing and mapping techniques for landslide inventorying.

**Keywords:** mass movements; inventory map; amalgamation; earth observation (EO); spatial resolution

#### **1. Introduction**

Landslides are the most frequent hazards of mountain regions throughout the world [1]. Given landslides' variable characteristics, they cause enormous damage to human life and infrastructure [2]. Landslides are usually caused by a trigger, like an earthquake or rainfall, and these two phenomena are considered to be the common physical triggers for event-based landslides [3]. In the Himalayan region, rainfall in the monsoon period triggers several massive and small landslides every year [4]. However, landslides triggered by earthquakes are severely destructive. For instance, the recent global landslide triggering events of the Wenchuan earthquake (2008), China [5], the Gorkha earthquake (2015), Nepal [6], and the Bihar earthquake (2002), India [7], resulted in a large number of casualties and severe damage to private and public infrastructures. However, there are several reasons that make it difficult to extract information about the exact location of landslides in an area, such as difficulty in accessing the hazard-affected remote areas [1]. A landslide inventory map, including the exact location and the exact boundaries along with the distribution, is the prerequisite for landslide analysis, susceptibility assessment, and mapping [8,9]. Furthermore, for the case of event-based landslides, detailed and state-of-the-art information about the landslides is critical. To better understand the triggering factors in an event-based landslide, different aspects of tracking, recording, and analysing data must be considered [10]. Several definitions are available regarding landslide inventories in the

literature. According to [11], landslide inventory maps are the basis for obtaining records about the date of occurrence, location, and type of slope movements. A dataset of a landslide inventory can provide information on the time of the event, the location of occurrence, the type of landslide, and the extent of the landslide [12].

Landslide inventory maps are the basis of determining hazard and risk assessment [13,14]. Golovko et al. [15] analyzed multiple sources of slope failures for the establishment of a comprehensive multitemporal landslide inventory. They also described landslide inventories as the prerequisite to enable landslide hazard assessment. Landslide inventories are also crucial in carrying out a risk analysis by analyzing the impacts of past landslide events and relating them to the present criteria to predict future landslide-prone areas [16]. Thus, it is crucial to record required information on landslide occurrence to link this with triggering factors.

Landslide inventory maps can give information about probable threatened areas to disaster management authorities whcih can be used for reconstruction planning after an earthquake event [6]. Inventories are also a basis for training and validating various knowledge-based, machine learning, and deep learning methodologies related to automatic landslide detection [1,14].

Event-based landslide inventory maps can be prepared from various sources. Recently, the availability and use of high-resolution remote sensing optical images has been very useful in the identification of landslides [14]. For instance, triggering events like earthquakes trigger thousands of landslides in remote areas. Therefore, remote sensing and earth observation (EO) data play a significant role in mapping and analyzing inventories. The EO data of pre- and post-landslide events are required for conducting classification and interpretation of the hazard-affected area. There have been several attempts to map landslides using expert-based approaches, such as manual rule-based, automatic, and semi-automatic classification techniques. Two main approaches for the classification and extraction of landslides from the EO data, namely object-based and pixel-based, were distinguished [1]. However, the spatial resolution of the available EO data plays a critical role in the quality of the resulting landslide inventory maps [17,18]. Although most studies so far relied on EO data with a single scale for landslide extraction, some works considered multi-resolution satellite imageries and EO data [18]. Even using a single satellite image frame, some studies performed multi-scale methodologies and observed higher landslide extraction accuracies compared to single-scale performances [17]. Scientific progress toward pixel-based automatic identification of landslides was made using remote sensing imageries with different resolutions [19]. Considerable progress toward pixel-based automatic identification of landslides was achieved using deep-learning convolution neural networks (CNN) [1,14]. Ghorbanzadeh et al. [1] applied CNN and other machine learning models to identify landslides in RapidEye data from the Rasuwa district in Nepal. The extracted landslides from different models were then tested using Global Position System (GPS) data along with a manually detected landslide inventory of image spectral features from RapidEye data; topographic input was also used, including digital elevation models (DEM) from Advanced Land Observing Satellite (ALOS) data. In another pixel-based study [20], unsupervised classification resulted in the detection of about 60% of manually extracted and mapped landslides. There were also several object-based studies for landslide extraction, such as [21] and [22], which compared their results with manually extracted results. Therefore, manual landslide extraction and mapping is considered a standard technique for receiving the most detailed and accurate inventory. As manual detection and extraction is the most reliable technique for inventory generation, it is preferable for testing state-of-the-art models, such as deep-learning CNN models. Therefore, the quality of the manually detected landslide inventory is critical, as it usually considered to be the ground truth [14]. There are a limited number of landslide and mass movement detection studies evaluating the quality of their applied inventory dataset. In a specific study [14], the frequency area distribution (FAD) method was applied to evaluate the three available landslide inventory maps, with one of them being selected for validation of the results.

On one hand, the quality of a landslide inventory map easily affects the overall accuracy of any landslide detection, susceptibility assessment, or hazard- and risk- mapping study. On the other

hand, there are a limited number of landslide studies that tackled the problem of comparison of two or more inventories [3,23]. Pellicani and Spilotro [23] compared archives and surveyed inventories for the Daunia region in Italy. They compared two landslide inventory maps to determine the corresponding quality through direct comparison. Another work by [3] compared photo-interpreted and semi-automatic landslide inventory maps in the Pogliaschina catchment, Italy. They compared the quality of rainfall event-based inventories cartographically and statistically.

Quality and completeness levels of a landslide inventory depend upon the accuracy, certainty, and type of information included in the map [24]. Criteria for the assessment of inventories are lacking from previous landslide research [3,25,26]. For the comparison of inventories, two or more inventories are needed to compare the quality of landslide maps. Despite the great importance of establishing the quality of a landslide inventory for scientific investigations, the number of such studies is limited. Comparing two or more inventories does not occurs often due to the limited availability of two or more inventories. However, there are some studies that were carried out regarding landslide mapping, and also several databases of landslides were compiled earlier by [27,28]. After the 1989 Loma Prieta, California, earthquake (Mw 6.9) a total of 1046 landslides were mapped using field investigations and aerial photographs in an area of about 15,000 km2 in central California. The spatial distribution of the landslides was investigated statistically using one-way analysis of variance (ANOVA) and regression techniques. Correlations of landslide occurrence with distance from the earthquake source, slope steepness, and rock type were determined [29]. A comprehensive database of devastating landslides caused by catastrophic earthquakes that took place all over the world was compiled by Rodríguez et al. [28], covering the period 1980–1997. Another work by Esposito et al. [30] discussed and described the ground effects and landslides triggered by the 1997 Umbria–Marche seismic sequence in an area of 700 km2. The environmental, seismic intensity (ESI) scale, and earthquake hazard were studied by [31,32]. The ESI scale is a measurement which defines the earthquake intensity by considering the size and spatial distribution of earthquake environmental effects. Lekkas et al. [33] used the ESI scale and its correlation with geological structures for seismic hazard estimation of the 2008 Mw 7.9 Wenchuan, China, earthquake. In another study by [34], correlations between ESI-07 intensity, slope, and lithology were discussed regarding landslides triggered by the 2016 Mw 7.8 Pedernales, Ecuador, earthquake.

Furthermore, Ferrario et al. [35] investigated the role of earthquake environmental effects within seismic sequences from the 2018 Lombok (Indonesia). Statistical analysis was carried out for three nearly complete landslide inventories triggered by the 12 May 2008, Wenchuan Mw 7.9 earthquake of China [36]. Correlations of landslide occurrence with topographical factors and seismic parameters were studied for three inventories. This literature review shows that, over the past decades, there was an improving trend regarding both the documentation and statistical evaluation of the earthquake-induced landslides in some important works.

Several of the nowadays commonly applied methodologies were mostly developed by the Italian groups for the census of the effects induced by earthquakes of moderate magnitude. During the Emilia Romagna (northern Italy) 2012 earthquake sequence, for processing and real-time data sourcing, new approaches and technologies were developed by the INGV EMERGEO working group [37]. Just after the earthquake event, the EMERGEO working group surveyed the epicentral area for co-seismic geological hazards. Later, they organized the records and processed them with the EMERGEO Information System (siE). The EMERGEO working team of Civico et al. [38] presented a 1:25,000 scale map of the surface ruptures after the 30 October 2016 Mw 6.5 Norcia earthquake, central Italy. They used 11,000 oblique photographs taken from helicopter flights, which were verified with field datasets. They also provided the datasets through a database of the co-seismic effects following the Norcia earthquake [39].

In this study, we evaluated four manually extracted landslide inventory maps of the Gorkha earthquake 2015, Nepal. Four different research teams from different parts of the world carried out landslide inventories related to this earthquake and published their resulting maps. These studies used various EO data sources, such as very high-resolution WorldView imagery and coarser-resolution imageries. We compared these inventory maps quantitatively and qualitatively using the cartographical degree of matching and frequency area distribution (FAD) methods [40]. The reasons for the differences are discussed here, outlining the limitations and advantages of different mapping techniques.

#### **2. Study Area**

The epicentre of the earthquake was in the Gorkha district and its aftershock was about 140 km in the Dholaka district [41]. However, landslides triggered by the earthquake and several aftershocks were scattered around large areas. Seven districts that were severely affected by the earthquake were selected for inventory mapping by different researchers. The study area covered a 14,502 km2 region, which spread over most of central Nepal. The area had two primary drainage systems, namely Narayani and Saptakoshi [42]. Our case study area lay in central Nepal country within the fold and thrust zone of Himalaya. This zone was caused by the collision of the Eurasian plate with the Indian Plate [43] (see Figure 1). The collision resulted in extensive crustal shortening and upheaval, leading to the formation of the quintessential collided orogen, Himalaya. All evaluations of the present study were done based on an overlapping region of four considered inventory maps.

**Figure 1.** Area covered by the different investigators during mapping landslide inventories of the 2015 Gorkha earthquake.

#### **3. Landslide Inventory Maps of Gorkha Earthquake (National Scale)**

#### *3.1. Inventory A*

Authors carried out a field survey and manual interpretation of Landsat-8 Enhanced Thematic Mapper (ETM) images with a resolution of 15 m, Gaofen-1 (GF-1) images with resolution of 2 m, GF-2 images with resolution of 0.8 m, and Google Earth imagery to prepare the polygon-based inventory [42]. A total of 15,456 km<sup>2</sup> area was selected for investigation based on previous studies and updates during the digitation process. A total of 3716 co-seismic landslides were mapped. The largest landslide mapped was 9983 m<sup>2</sup> and the smallest landslide was less than 50 m2. The co-seismic landslides covered a total area of 15.93 km2.

#### *3.2. Inventory B*

Earthquake-induced landslides were mapped by comparing pre- and post-event very highresolution satellite imageries. DigitalGlobe WorldView-2 and -3 imagery, along with Pleiades satellite data, were used [44]. The spatial resolution of imagery used in mapping landslides varied from 30–50 cm in most of the areas. Images were acquired from 26 April to 15 June 2015, with most images collected between 2 May and 8 May 2015. This team was the only group working on Gorkha earthquake who differentiated between landslide source and deposits areas. They were able to map more than 25,000 landslides in the area affected by the earthquake.

#### *3.3. Inventory C*

Authors used Google Earth imagery, which was updated after the Gorkha earthquake, to prepare a polygon-based inventory. A total of 4000 km<sup>2</sup> area was selected for investigation based on previous studies and updates during the digitizing process [45]. In most of the region, the satellite imagery of pre-earthquake was from December 2014 and post-earthquake imagery from 2–4 May, which was around one week after the main earthquake. A total of 17,000 co-seismic landslides were mapped by [45]. Results also showed the spatial correlation of topographical parameters with landslide occurrence.

#### *3.4. Inventory D*

The landslide inventory map was manually mapped by the authors in the immediate aftermath of the earthquake using a range of EO data sources, including web-hosted high-resolution optical data in Google™ Crisis Response (e.g., United Kingdom - Disaster Monitoring Constellation-2 (UK-DMC2), Disaster Monitoring Constellation for the International Charter (DMCii), Worldview, Digital Globe Inc., SPOT National Centre for Space Studies (CNES), imagery accessed via the Disaster Charter, imagery available from United States Geological Survey (USGS) Hazards Data Distribution System (HDDS Explorer), and imagery specifically tasked over regions of interest (e.g., Pleiades CNES). A total of 2117 co-seismic landslides were mapped by [46] (see Figure 2).

**Figure 2.** Landslide mapped after the 2015 Gorkha earthquake for the affected region. (**A**) Inventory A, (**B**) inventory B, (**C**) inventory C, and (**D**) inventory D.

#### **4. Landslide Inventory Comparison Methodologies**

#### *4.1. Spatial Distribution*

Landslide inventories for the Gorkha event were collected from various sources. Detailed information on inventories is given in Table 1. Out of the inventories, we used four that were available and prepared by research teams from various parts of the world. To analyze the landslide density per square kilometre, point-based inventory data were used. Figure 3 shows the landslide density distribution of selected inventories for the Gorkha event; all of the inventories were co-seismic and prepared using post-earthquake satellite imageries.

Information regarding the data used to prepare the inventories is given in Table 1. There was a significant variation in the number of landslides mapped for the Gorkha event, ranging from *NLT* = 2117 to *NLT* = 24, 915. None of these inventories classified the types of landslides as proposed by [47]. The data sources used for mapping varied from high-resolution imagery of about 30 cm resolution to Google Earth images. The landslide density analysis results, illustrated in Figure 3, showed that most of the landslides were concentrated in the Gorkha and Sindhupalchock districts for most of the inventories; the southern part of the Rasuwa district was also severely hit, as can be interpreted from the Figure 3. Differences in the density of landslides in an area depended on factors such as the effect of amalgamation while mapping, the purpose of the mapping, and the data sources used for mapping. Areas near the Dholaka and Sindhupalchok districts showed high landslide density ias a result of the aftershock on 12 May 2015 in the Dholaka district. The higher density of landslides showed that inventories were prepared after the aftershock of 12 May 2015 and hence had a larger number of landslides compared to inventories prepared for the first earthquake event. Another factor that affected the density of landslides was the different coverage of the mapping area. Comparing the different inventories at larger scale is not an easy task as their coverage areas are different. To compare landslide inventories prepared by different interpreter, we chose a commonly mapped area that was

mapped by most of the inventories. For the analysis, we chose four inventories from various sources which were polygon-based, and compared them statistically for smaller regions.

**Table 1.** Landslide inventories prepared after the 2015 Gorkha earthquake, Nepal.


**Figure 3.** Landslide inventory densities per km<sup>2</sup> for the Gorkha event. (**A**) Inventory A, (**B**) inventory B, (**C**) inventory C, and (**D**) inventory D.

To compare the inventories for the Gorkha earthquake, we chose four polygon-based inventories that were available. The commonly mapped area was covered in all four inventories, which allowed for the statistical comparison of the inventories. We compared the four inventories and examined the differences and similarities of the total area covered, the number of landslides mapped, and the size of the landslides.

#### *4.2. Cartographical Degree of Matching*

For the comparison of four inventories available for the common mapped area in the Rasuwa district, an attempt was made to determine the cartographic matching and mismatch (see Figure 4).

**Figure 4.** (**A**) Example of a union with features within a feature class that overlaps. (**B**) Example of an intersection with a feature within a feature class that overlaps.

Carrara et al. (1993) [52] provided a method to evaluate the degree of matching and mismatch between two inventory maps. The mismatch index, *E*, was given by

$$E = \frac{(A\_1 \cup A\_2) - (A\_1 \cap A\_2)}{(A\_1 \cup A\_2)}, \ 0 \le E \le 1\tag{1}$$

$$M = 1 - E,\\ 0 \le M \le 1 \tag{2}$$

where E is the mismatch and M is the matching between the inventories. Union and intersection are used in Equation (1) to observe the matching and mismatch between the inventories.

This equation was valid for up to two inventories. However, in our case, we wanted to compare four inventories. Thus, we formulated an equation for the comparison of four inventories (see Equation (3)).

$$E = 1 - \frac{(A \cap B) \cup (A \cap C) \cup (A \cap D) \cup (B \cap C) \cup (B \cap D) \cup (C \cap D)}{(A \cup B \cup C \cup D)} \tag{3}$$

The statistics of the four landslide inventories can be seen in Table 2. The total number of landslides varied from 33 to 49 for the same area. Moreover, there was a significant difference in the total area of landslides. Cartographical mapping differences in a single landslide for the same event can be seen in Figure 5. Interpreters mapped the landslide boundary differently, and the causes of such differences are discussed in Section 5. Based on Equation (3), we observed cartographical match and mapping errors, as presented in Table 3. The pairwise comparison of inventories and comparison of inventories mapped by different interpreters are represented in Figures 6 and 7, respectively.

**Table 2.** Statistics for the four event-based landslide inventories.


**Figure 5.** Mapping by different interpreters for the Gorkha earthquake of the same landslides. (**A**,**B**) represent examples of the mapped landslides.

**Figure 6.** Pairwise comparison of inventories mapped by different interpreters.

**Figure 7.** Comparison of inventories mapped by different interpreters.


**Table 3.** Cartographical matching and mismatch between inventories.

#### *4.3. Frequency Area Distribution (FAD)*

Landslide inventories were statistically analyzed using frequency area distribution (FAD) curves, in which the landslide areas were plotted versus the cumulative and non-cumulative landslide frequencies. In the study by [53], observations showed that the power law was valid for medium and large landslides. The probability of occurrence of landslide size can be given by the power-law equation.

$$p(\mathbf{x}) = \varepsilon \mathbf{X}^{-\beta},\tag{4}$$

where *X* is the observed values, *c* is a normalization constant, and β is the power-law exponent.

Figure 8 shows the power-law distribution for medium to massive landslides and divergence from the power-law toward lower frequencies with a rollover point, where frequency decreased for smaller landslides. The trend of the FAD of most landslide inventories diverged from the power-law for small landslides [53–56]. The point where this divergence began was defined as the cut-off point [56,57]. For non-cumulative probability density distributions of landslide areas, the peak point of the probability distribution curve, after which the probability value began to decrease for smaller landslides following a positive power-law decay, was referred to as the rollover point [58]. According to [58], in a power-law distribution, the slope of the distribution was defined with a power-law exponent. The part that was represented by large events was referred to as the power-law tail, as shown in Figure 9 (with a scaling parameter, β). Malamud et al. [53] investigated four well-documented landslide events and concluded that rollover was a real phenomenon for landslide-event inventories, depending upon the bias and under-sampling of the smaller landslides. They modeled the FAD for these four inventories and established theoretical curves to estimate the total landslide area triggered by an earthquake or rainfall event.

**Figure 8.** Schematic representation of the main components of a non-cumulative frequency area distribution (FAD) for a landslide inventory.

Further distributions were used to fit the frequency area distribution of landslides. The double Pareto model described the majority of the data well, but [59] indicated that the same model was less good at the tails of the distribution. Another method was proposed by [60], who showed that the entire FAD of landslides could be explained by a three-parameter inverse-gamma distribution (equation). This approach also described a way to estimate the landslide event magnitude (*mLs*). The *mLS* is the indication of the size of the landslide triggering event and gives an indication of the severity of the event in terms of landslide occurrence in a particular area for an event.

$$p(A\_L; \rho, a, s) = \frac{1}{a\Gamma(\rho)} \left[\frac{a}{A\_L - s}\right]^{\rho + 1} \exp\left[-\frac{a}{A\_L - s}\right] \tag{5}$$

where ρ is the parameter primarily controlling power-law decay for medium and large values, Γ(ρ) is the gamma function of ρ, *AL* is landslide area, *a* is the location of rollover point, *s* is the exponential decay for small landslide areas, and −(ρ + 1) is the power-law exponent. Malamud et al. [60] provided a best fit for the power-law exponent and showed that −(ρ + 1) = 2.4.

Table 4 shows that the power-law exponent of four analysed inventories ranged from 2.27 to 2.48, which was consistent with the literature describing an interval having a central tendency of around 2.3–2.5 [55,58]. The minimum landslide area mapped ranges from 35.50 m<sup>2</sup> to 500.53 m<sup>2</sup> for the inventories. Also, the largest landslide mapped ranges from 764,038 m<sup>2</sup> to 157,265.25 m2. The rollover points ranged from 256.74 m<sup>2</sup> to 1258 m2 for the inventories.


**Table 4.** Comparison of the frequency area statistics of the landslide area.

**Figure 9.** Landslide frequency size distribution, representing the dependence of landslide probability density p on the landslide area.

There was a zigzag pattern in the plotted figure of landslide probability density against the inverse gamma fit. The differences in the probability distribution and inverse gamma fit may have been the result of gaps regarding mapped landslides for given inventories, indicating that some landslides were missing or not mapped by interpreters due to various reasons, such as rapid mapping or the amalgamation of smaller landslides into single landslide features. The rollover points for all inventories differed to each other. For the inventory by interpreter B, the rollover point toward the smaller landslides was 85.13, which was smaller in comparison to other inventories, as there was a wider distribution of small landslides in this inventory.

#### **5. Discussion and Conclusions**

The availability of four independent landslide inventories for the same triggering event in the same geographical area allowed us to quantitatively compare and outline strengths and weaknesses of the methods used to prepare the inventories. The inventory by interpreter B (Figure 10b) portrayed more landslides (498) than the other inventories. This difference was significant, as the other interpreters could not detect most of the smaller landslides visually using Google Earth imagery, which was used to prepare the three inventories by interpreters A, C, and D. The more significant number of landslides in the inventory by interpreter B with the very high-resolution Worldview imagery used for the visual identification of the landslides, compared to the spatial resolution of the satellite imagery used for the other inventories, is explained here. The amalgamation of the landslide information merged multiple very small landslides into a single larger landslide, reducing the total number of the mapped event landslides. Considering only the number of landslides, quantitative identification could not be performed because of both the amalgamation of adjacent landslides and the subjectivity of the landslide extraction process. Many factors caused the amalgamation of landslides in a landslide inventory map. The manual extraction of landslide borders and representation with polygons is a subjective process that is affected by the applied method, the preferences of the experts and interpreters, and how much time and effort are invested into the inventory generating process [61]. Amalgamation is the mapping of nearby smaller landslides as a single polygon, which may lead to possibly severe distortion of the statistical analysis of these inventories. An adjacent landslide polygon is usually described as a single polygon if the runouts or scars overlap in areas. Thus, differentiating between them is difficult. Low image resolution, the working scale, and the contrast between affected and unaffected areas are all other reasons for amalgamation [62]. In some regions, landslides can be very condensed, and several contiguous landslides may join runout areas. Amalgamation is often due to errors resulting from a lack of expertise of the interpreter. It may also happen when landslide inventory mapping is carried out using (semi)-automated classification and change detection based on optical satellite images, e.g., [63].

There are different numbers of mapped landslides for the same area interpreted by different interpreters, which related to the scale of working, applied EO spatial resolution, personal perspective, and the method of mapping that was adopted. Inventory C (Figure 10C) and inventory D (Figure 10D) landslides were mapped using Google Earth imagery; 197 and 336 landslides were mapped, respectively. In Figure 10A, inventory A had 144 landslides for the same common area using Google Earth imagery. Figure 10B shows 498 landslides in inventory B, which was the result of using high-resolution Digital Globe WorldView 2–3 satellite imagery by the author. However, it should be mentioned that very high resolution (VHR) optical satellite imagery can have significant distortions and georeferencing errors in high-relief areas if they are orthorectified using only provided rational polynomial coefficient (RPC) models and not by using manual ground control points (GCPs). Also, this can be more severe, as many of the VHR rapid emergency acquisitions over a disaster area-of-interest can have small incidence angles, thus increasing automatic orthorectification errors. Differences in the mapped outlines of landslides can be due to this. The offset from automatic RPC processing is no less than 5 m for the newer sensors (WorldView 2- 3, Pleiades) and even less for older satellites (GeoEye-1, etc). High mountain relief and low-resolution digital elvetaion model (DEM) used for orthorectification (e.g., Shuttle Radar Topography Mission (SRTM), Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) or Advanced Land Observing Satellite (ALOS) may add more possible errors to these offsets [64].

**Figure 10.** Effect of amalgamation for the same area mapped by four different people in the area near the Trishuli river, (**A**) Inventory A, (**B**) inventory B, (**C**) inventory C, and (**D**) inventory D.

The results of our mismatch index in this work confirm that the difference between the four-event inventories was significant. However, the mismatch index was in the range of the differences measured by other investigators in previous studies that compared landslide inventories in similar physiographical settings (e.g., [25,52]). In addition to the causes for the mismatch discussed by these investigators, in our test case, the difference was also the result of the period of landslides mapped and of the different spatial resolutions of the satellite imagery. Despite the mismatch index, E, visual inspection of the four inventories revealed a similar spatial distribution of the event landslides in the four landslide maps. This matter was confirmed (i) by the spatial correlation between the four inventories and (ii) by the similarities of the double Pareto density functions for the four inventories (see Figure 10). The results of the power-law exponent of the four analysed inventory ranged from 2.27 to 2.48, which was consistent with the literature that showed that the power-law exponent interval had a central tendency around 2.3–2.5 [55,58].

According to the results of the present study, we conclude that both comparison methods, i.e., the cartographical degree of matching and FAD, are essential for evaluating the quality of landslide inventory productions in terms of their similarities in the total number, area size, and spatial density of landslides. However, the cartographical degree of matching method is suitable more for the evaluation and validation of the location and boundaries of the landslide affected areas. Therefore, the cartographical degree of matching method is essential for landslide extraction/annotation studies, especially those using machine learning models for this aim. The accuracy of the resulting landslide extraction from remote sensing data by using the machine-learning models is directly related to the quality of the training data of landslide inventory datasets [1]. Thus, the cartographical degree of the cartographical matching method plays a vital role in enhancing the accuracy of any landslide extracting studies using machine learning models.

On the other hand, although the FAD is also advantageous for landslide detection studies, it is a crucial method for studies which analyse the landslide areas and their distribution in a case study area. Landslide susceptibility mapping, landslide-prone region analysis, and landslide modeling and risk assessment are considered to be studies, where the FAD method should be considered to evaluate the applied inventory datasets.

For the Gorkha earthquake-affected region in central Nepal, we compared four independent event landslide inventories showing landslides triggered by a high-intensity earthquake that hit the area on 25 April 2015. The first inventory was obtained through the visual interpretation of Google Earth imagery. The second inventory was obtained by exploiting a semi-automatic procedure applied to a VHR resolution worldview satellite imagery. We compared the four inventories by exploiting methods already present in the literature and by proposing new qualitative and quantitative criteria. Comparison of the four independent event inventory maps led to the conclusion that the mismatch between the four inventories was significant but consistent with differences measured by other investigators in similar physiographical areas. The mismatch was attributed to (i) different spatial resolutions of the satellite images, (ii) the amalgamation of smaller landslides into larger landslides while delineating the boundaries of landslides; the minimum landslide area was 35.50 m2, which was mapped using VHR resolution worldview satellite imagery, whereas, the landslide mapped using Landsat 8 ETM imagery was 500.53 m2, and (iii) the inability of the operator to recognize landslides in shadowed areas. For our future work, we will apply these inventories along with some data-driven models to generate landslide susceptibility maps and compare the resulting susceptibility accuracies.

**Author Contributions:** Conceptualization, S.R.M.; methodology, S.R.M., and S.T.P.; validation, S.R.M.; data curation, S.R.M.; writing—original draft preparation, S.R.M., and S.T.P.; writing—review and editing, S.R.M.; visualization, S.R.M. All authors read and approved the final manuscript.

**Funding:** This research is partly funded by the Open access fund by University of Salzburg

**Acknowledgments:** We would like to thank three anonymous reviewers for their valuable comments.

**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* **Dynamics of the Zones of Strong Earthquake Epicenters in the Arctic–Asian Seismic Belt**

**Lyudmila P. Imaeva 1,2,\*, Valery S. Imaev 1,3 and Boris M. Koz'min <sup>3</sup>**


Received: 28 February 2019; Accepted: 9 April 2019; Published: 12 April 2019

**Abstract:** Our comprehensive study of the Russian Arctic region aims to clarify the features and types of seismotectonic deformation of the crust in the Arctic–Asian Seismic Belt, specifically in the zones of strong earthquakes in the Laptev Sea Segment, the Kharaulakh Segment, and the Chersky Seismotectonic Zone. We have analyzed modern tectonic structures and active fault systems, as well as tectonic stress fields reconstructed by tectonophysical analysis of the Late Cenozoic faults and folds. The investigated neotectonic structures are ranked with respect to the regional classification principles. Changes in the crustal stress–strain state in the lithospheric plate boundaries between the Eurasian, North American, and Okhotsk Sea Plates are analyzed, and regularities of such changes are discovered. A set of models has been constructed for the studied segments of plate boundaries with account of the dynamics of the regional geological structures. The models can give a framework for the assessment of potential seismic risks of seismogenerating structures in the Russian Arctic region.

**Keywords:** Arctic–Asian seismic belt; regional segment; active fault; paleoseismogenic structure; Late Cenozoic deformation; earthquake mechanism; seismotectonic deformation; potential seismicity

#### **1. Introduction**

The Russian Arctic region is covered by comprehensive studies combining geological, geophysical, and seismological methods that have identified the Arctic–Asian seismic belt, which includes the spreading Gakkel Ridge, the system of rift basins in the Laptev Sea shelf, and the seismogenerating structures of the continental crust in the Chersky Seismotectonic Zone [1]. The geodynamic processes in this belt are studied as indicators of types of seismotectonic deformations of the crust in the contact zone where the Eurasian, North American, and Okhotsk Sea lithospheric plates interact and move relative to each other (Figure 1).

This paper reviews published results regarding the Laptev Sea Segment, Kharaulakh Segment, and the Chersky Seismotectonic Zone [2,3], and presents new research results. It discusses the dynamics of the formation of neotectonic structures and the types of the crustal stress–strain state in the zones of strong earthquake epicenters in the Arctic–Asian seismic belt. In each of the above-mentioned segments, the structural–tectonic positions of the modern structures and the systems of active faults are identified and analyzed on the basis of the field geological and geostructural data collected by the authors and the literature data.

Our research provides data for discovering the tectonic positions and the structural dynamic pattern of the main fields of earthquake epicenters in the study area and allows identifying the blocks that acts as tectonic stress concentrators. Changes in the crustal stress–strain state in different tectonic settings at the margins of the Eurasian, North American, and Okhotsk Sea lithospheric plates are analyzed, and the regularities of such changes are discovered.

**Figure 1.** Geodynamic activity of neotectonic structures of the Arctic–Asian seismic belt (modified after [1]). 1—classes of geodynamic activity of the domains: 1–2—low activity, 3–5—medium activity, 6–9—high activity; 2—earthquake epicenters sized by magnitude: ≤ 4.0, 4.1–5.0, 5.1–6.0, 6.1–7.0, ≥ 7.0; 3—kinematics of active faults: (a) thrusts, (b) normal faults, (c) strike-slip faults; 4—Balagan–Tas volcano; 5—focal mechanisms of earthquakes (lower hemisphere; the principal stress axes of compression and extension are marked by black and white dots, respectively). Lithospheric plates: EU—Eurasian, NA—North American, OK—Okhotsk. Inset—location of the Arctic–Asian seismic belt. Segments: I—spreading Gakkel Ridge, II—Laptev Sea rift, III—Kharaulakh, IV—Chersky Seismotectonic Zone.

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

Seismotectonic studies are based on the concept of the structural and dynamic uniformity of the geophysical medium and regularities in the development of seismogeodynamic processes, which is elaborated by the Institute of Physics of the Earth of the Russian Academy of Sciences (IPE RAS) [4,5] and developed in [1–3,6,7]. The research method used in our studies of the Russian Arctic includes several stages. The first stage aims at establishing the general trends in the neotectonic development of the study area. Modern (Late Cenozoic) structures in the study area are analyzed considering the recent tectonics as a structural framework comprising active faults and other features of modern tectonic activity related to regional seismicity. Neotectonic structures of the Russian Arctic are ranked by their degree of activity, according to the regional principles of the classification described in [6,7].

In our study, a domain is a neotectonic geodynamic taxon of the territorial rank, which is considered to be a spatially localized integral object with a multifactorial interaction of its main components in the profile of the earth's crust. The classification of domains is a multi-level system, including nine levels (i.e., classes) of activity of modern geodynamic processes that lead to the formation of neotectonic structures. Each activity class is characterized by its specific set of primary and additional features pertaining to tectonics (geodynamic settings), geophysics (seismicity, heat flow, field gravity anomalies, and crustal thickness), morphostructure (terrain elevation, difference between the highest

and lowest elevations, and rates of vertical and horizontal movement of the ground surface), material composition, deformation indicators, and GPS measurements. Additionally, the inherited dynamics of neotectonic structures are considered with respect to the conditions in the previous stages of the domain development. To assess the modern geodynamic activity and specify the class of each domain, we interpreted both its primary and additional features. The nine classes of domains are grouped by the degree of the modern geodynamic activity: low (classes 1–2), moderate (classes 4–6), and high (classes 7–9) activity. The characteristics of the classes and the data on active faults are shown in Figure 1.

The next stage of the study aims at identifying the most probable areas of the recent activity for a more detailed investigation and search for reference objects. Large-scale morphostructural and structural–dynamic mapping is carried out to provide the two components of morphotectonic analysis. The main concept of this analysis is consistency between terrain features and corresponding rates and types of endogenous processes. The relative movements of crustal blocks during neotectonic activation cycles create the main features of the terrain, specifically the morphostructures bordered by active faults. The types of endogenous geodynamics are reflected in the features of modern geodynamic activity in the blocks of different ranks and linear fault zones between the blocks.

The database, including geological, geophysical, and seismological data on the study area, provides the basis for investigating the structural–dynamic features of the main fields of earthquake epicenters characterized by the maximum seismic potential. To this end, this study stage includes the collection and interpretation of large-scale remote data and laser-scanning images of the isoseismal areas of strong earthquakes, mapping active faults, fault kinematic analysis, and selection of areas to be covered by detailed field surveys. Attention is given to additional indicators of recent fault activity, such as displacements identified by repeated geodetic surveys, earthquake epicenters confined to the zones of dynamic influence of faults, focal mechanisms of earthquakes as indicators of the dynamics and directions of crustal movements, geothermal and gas-hydrochemical anomalies that give evidence of an increased permeability of the crust, seismic profiling, seismic survey data, gravimetric data, and electrical survey data. Field studies collect and clarify information on the deformation and displacement of young relief elements and sediments and discover evidence of strong paleo- and modern earthquakes. Trenching is performed on seismogenic deformation sites. The field database helps to clarify the kinematic types of Late Cenozoic folds and faults and the structural parageneses of active faults.

Reference objects (i.e., zones of earthquake epicenters) are selected in the seismogeodynamic zones and studied in detail. In such zones, we identify linear fault zones and blocks that act as tectonic stress accumulators, which may have high seismic potential; determine the kinematic types of contact zones of the main seismogenerating structures; and create models showing regional structures and dynamics. The experience of some Russian and international joint research projects shows that the above-described sequence of data collection and processing ensures that the resultant datasets provide a complete and reliable investigation of modern seismogeodynamic processes.

The research results reported in this paper are based on an updated and more comprehensive regional database consolidated by the authors, as well as the data on geology, tectonics, geophysics, seismogeology, and hydrogeology of the study area from publications and sources provided by industrial companies and research organizations.

#### **3. Results**

#### *3.1. Laptev Sea Segment*

#### 3.1.1. Structure and Tectonics of Laptev Sea Segment

The Laptev Sea shelf is located at the northern margin of the Eurasian Plate (Figure 1). In this area, the SE flank of the mid-oceanic Gakkel Ridge is traced orthogonally to the continental slope (Figure 2). The oceanic basin is in contact with the shelf areas of the East Arctic seas, which represent a system of

structures between the continent and the ocean. This area experienced several phases of collisional deformation in the Late Paleozoic and Mesozoic [8,9].

**Figure 2.** Seismotectonics of the Laptev Sea segment (modified after [7]). 1—Gakkel spreading ridge; 2—continental slope; 3—grabens at the bottom of the Laptev Sea: L—Lyakhov, B—Belkovsky–SvyatoyNos, S—Shirostonsky, Ch—Chondonsky, Y—Ust–Yana, U—Ust–Lena; 4—conventional boundaries of the Laptev Sea plate: S—southern, W—western, E—eastern; 5—boundaries of large troughs and uplifts; 6—earthquake epicenters sized by magnitude: ≤3.0, 3.5, 4.0, 4.5, 5.0, 5.5, 6.0, 6.5, 7.0; 7—focal mechanisms, dates and magnitudes of earthquakes (lower hemisphere; the principal stress axes of compression and extension are marked by black and white dots, respectively). Lithospheric plates: EU—Eurasian, NA—North American.

At the end of the Cretaceous and Cenozoic, intensive rifting took place due to opening of the Eurasian spreading basin [8,9]. The Laptev Sea Rift has been described in detail using the data from marine multichannel seismic profiling by reflected waves (CDP–SRM) [10–12]. According to these surveys, the entire Laptev Sea area is framed by the Mesozoic fold systems that are traced for considerable distances into the sea. In seismic profiles, the fold systems are recorded as an acoustic basement and can be viewed as a tectonic base of the Laptev Sea Rift. The basement top is disturbed only by normal faulting, which gives evidence of the post-folding age of this surface and suggests that the sedimentary cover may be dated to the Cenozoic or Late Cretaceous.

The large thickness of the sedimentary cover (to 10 km) is indicative of the presence of both Cenozoic and Cretaceous deposits in the Laptev Sea Rift. The basement includes the Paleozoic and Early Mesozoic formations of the Verkhoyansk–Kolyma fold system, wherein folding was completed by the Middle Cretaceous. During the Cenozoic, subsidence occurred on both the continental shelf and coastal lowlands. The profiles of the exhumed Cenozoic sediments on several coastal lowlands and islands have been studied in detail. On the Anzhu Islands, the Upper Cretaceous and Cenozoic deposits, and even the Miocene sediments, are disturbed by linear folds in thrusts that are non-conformably overlain with horizontally deposited Upper Pliocene sediments [13]. It can thus be suggested that Cenozoic extension was interrupted by a compression episode at the end of the Miocene. According to the seismic profiling, a system of narrow NW-striking grabens, troughs (Ust–Lena, Ust–Yana, Chondon, Bel'kovsky–Svyatoy Nos, etc.), and associated submarine horsts have been identified in this area (Figure 2) [11,14–16]. The troughs are up to 200–250 km long and 40–60 km wide. The most remarkable element of this rift system is the Ust–Lena graben, traced for a distance of 400–420 km north of the southern termination of Buor–Khaya Bay; its northern part is 150–170 km wide, and gradually narrows to 30–40 km towards the south. Faults bordering the grabens control the details of their interior structure. Two types of faults are typical of this area: sublatitudinal and NW-striking normal faults. Strike-slips faults oriented sub-orthogonally to the normal faults are also observed. Such faults form the E–NE, NE, and NW-trending fault systems that have horizontally displaced the deposits in the depressions by almost two kilometers.

#### 3.1.2. Seismotectonics of Laptev Sea Segment

The Laptev Sea Segment is a zone of diffuse seismicity, which spatially corresponds to the system of rift basins in the Laptev Sea shelf and occupies the area between the Tajmyr Peninsula, Lena River delta, and Novosibirsk Islands (Figure 2). This zone includes several NW-striking subzones of more densely spaced earthquake hypocenters. Strong earthquakes mainly tend to occur in the zone that extends from the Gakkel Ridge to the Yana Bay of the Laptev Sea and follows the boundary between the Eurasian and North American lithospheric plates. In this zone, the earthquake sources are either concentrated in the basins of the Laptev Sea rift system or occur at the sides of the depressions. Furthermore, two zones without much seismic activity are identified along the boundaries of the Laptev Sea shelf. The Lena–Tajmyr zone stretches across the Lena River delta, along the coast of Olenek Bay of the Laptev Sea, and towards the Tajmyr Peninsula; and another zone extends sublongitudinally from the East Siberian Sea between the Faddeev and New Siberian islands and can be traced further northwards (Figure 2).

In our previous studies [1,7], seismotectonic deformation parameters were calculated from the focal mechanisms of local earthquakes recorded at the southeastern termination of the Gakkel Ridge, the continental slope and the rift depressions of the Laptev Sea shelf (Figure 2). Seismic moment tensors were taken from the Global Centroid Moment Tensor (CМТ) Catalog and the International Seismological Centre (ISC) Bulletin [17]. According to our calculations, the seismotectonic setting of extension is evidently dominant in the above-mentioned areas of the Russian Arctic. The principal stress axis of extension is subhorizontal and oriented in a NE–SW direction across the strike of the main structures [7].

Considering the regular coincidence of earthquake epicenters with the tectonic structures of the Laptev Sea shelf, we distinguish the three most seismically active portions of the Laptev Sea Segment, which are termed 'eastern', 'southern', and 'western' (Figure 2). The eastern portion is traced along the Belkovsky–Svyatoy Nos and Lyakhovsky grabens and coincident with the northern flank of the Mesozoic fold zone belonging to the Verkhoyansk–Kolyma fold system. The southern portion is traced along the latitudinal Mesozoic fold branch of the Olenek sector in the portion of the Lena–Anabar regional suture (Figure 3). The western portion is traced from the Tajmyr fold system along and towards the abyssal basin. In plan, these seismogenerating structures contour the margins of the Laptev Sea Microplate [1,3,7]. According to the seismological data, the western and eastern boundaries of the Laptev Sea Microplate are subjected to compression (Figures 1 and 2). In this area, thrusting caused local earthquakes in the contact zones between the continental and rift structures.

From the south, the Laptev Sea Rift is separated from the Siberian platform by the Mesozoic faults of the Olenek section of the Lena–Anabar regional suture (Figure 3). In this area, a series of W–NW-trending folds formed in the Mesozoic due to sublatitudinal left-lateral shearing along the northern margin of the Siberian platform. Such folds may continue on the Laptev Sea shelf, wherein their possible extension is limited by an extended zone of high-gradient gravity anomalies. The latter is detected as a sublatitudinal system of alternating small-size linear positive and negative anomalies of varying intensity [18]. In this zone, the seismic records show several areas of maximum seismic activity values, which correlate with the above-mentioned system of gravity anomalies (Figure 3).

**Figure 3.** Seismotectonics of the Olenek sector of the Lena–Anabar regional suture (modified after [7]). 1—Lower Cretaceous continental deposits; 2—anticline axis; 3—syncline axis; 4—isoheights of the Lower Cretaceous bottom; 5—axis of consedimental bar; 6—axis of consedimental basin; 7–9—kinematics of active faults: 7—thrust, 8—normal fault, 9—strike-slip fault; 10—focal mechanisms, dates and magnitudes of earthquakes (lower hemisphere; the principal stress axes of compression and extension are marked by black and white dots, respectively); 11—earthquake epicenter; 12—density of earthquake epicenters (number of seismic events per 1◦ × 1◦): (a) 26–30, (b) 21–25, (c) 16–20, (d) 11–15, (e) 6–10, (f) 1–5.

In the zone of the Lena–Anabar regional suture, the seismic process develops in conditions of both extension (Lena River delta, shores of Olenek, and Anabars bays) and compression (Taimyr Peninsula). Focal mechanism solutions for this zone are varying, as well as the fault kinematics, including normal, reverse, and strike-slip faults, and their modifications. The analysis of seismotectonic deformation structures shows the dominance of crustal stretching with a small shear component [6,7]. The extension stress azimuth in this zone differs in orientation from that of the Laptev Sea Rift. The principal stress axes are directed across the strike of the main tectonic elements; those with low dip angles are NW–SW-trending (Figures 1 and 3). Most focal mechanisms of the local earthquakes show normal faulting, with the exception of the 1990 М<sup>s</sup> 4.9 Taimyr earthquake in the western margin of the Laptev Sea shelf and the 2015 М*<sup>s</sup>* 4.3 earthquake in the NW part of the Taimyr Peninsula (Figure 2) [3,7]. Their focal mechanism solutions significantly differ from those of all other seismic events recorded in the area along the Laptev Sea coast (Figures 1 and 3).

The 1990 Taimyr earthquake, which occurred at depth of 15 km, caused thrust-and reverse-type displacements, along the gently sloping and subvertical planes of the NW-striking and submeridional

faults, respectively. It should be noted that the conclusion on compression in the Taimyr Peninsula, which is based on the focal mechanisms, is consistent with the results of the geostructural studies based on the geological survey data [19]. It is also supported by the rates of visible uplifting of the Laptev Sea coast, which were estimated from long-term measurements of the global ocean level. Over a 10-year measurement period, the rates of modern tectonic uplifting in the Taimyr Peninsula were +(1–2) mm/yr on average [20].

Thus, the modern geodynamic setting of the Laptev Sea Segment is determined by the seismogenerating structures located in the area where the riftogenic zone of the Eurasian basin of the Arctic Ocean extends up to the shelf area. This predetermines the style of seismotectonic deformation, the locations of zones with specific tectonic settings, and the crustal stress-strain state, as well as the structural and dynamic features of the main fields of earthquake epicenters in the study area. In general, the seismotectonic deformation parameters reflect the trends in the distribution of crustal stress fields reconstructed from the structural data.

#### *3.2. Kharaulakh Segment*

#### 3.2.1. Structure and Tectonics of Kharaulakh Segment

The Kharaulakh Segment is the northern flank of the Verkhoyansk fold–thrust belt (Figure 4). In the Neo-Proterozoic, this segment originated at the reworked margin of the Siberian platform and developed as a passive continental margin. Its evolution is reflected in the structures and kinematic types of the observed dislocations. During the Cenozoic, the rotation pole of the North American and Eurasian plates changed its position several times and thus caused alternating extension and compression periods, as confirmed by the structural tectonics of this area [7,9,21]. In the Kharaulakh Segment, the Cenozoic megacomplex is represented mainly by Paleocene-Eocene continental deposits that occur with a sharp angular unconformity on different horizons of the deformed Precambrian–Mesozoic megacomplex [6,7]. These deposits fill a series of sublongitudinally oriented depressions, e.g., Kengdei, Kunga, Sogin, Bykov, etc. (Figure 4), which formed in the Paleogene during the earliest rifting stage along the continental continuation of the spreading Gakkel Ridge. At some locations, the Paleogene deposits are folded and cut by thrusts and reverse faults.

All the above-mentioned observations provide evidence of compression in the Cenozoic. The structural studies of the Kharaulakh Segment (Figure 4) suggest sublatitudinal compression in the Middle Miocene, according to the Cenozoic profiles of this segment and adjacent areas [21].

The next episode in the Cenozoic evolution of the Kharaulakh Segment was extension in the Pliocene–Quaternary. Young normal faults that displaced the Neogene weathering crust are observed along the coast of Buor–Khaya Bay (Figure 4). The extension axis was either sublatitudinal or NE-trending. The Quaternary evolution of the near-coastal zone of the Kharaulakh Segment was sharply different from the development of its continental part. Normal faulting dominated in this area and defined the block structure of this territory, as shown by differences in the hypsometric positions of the Late Quaternary and Holocene deposits. The data on normal faults give evidence of the extension phase and suggest that the extension axes were sublatitudinal and NE-oriented [7]. Thus, the Kharaulakh Segment is a transition zone where the mid-oceanic and continental crustal structures are conjugated [9,16,21].

**Figure 4.** Seismotectonic map of the Kharaulakh segment (modified after [21]). 1—Cenozoic depression: (a) Kengdey, (b) Khorogor, (c) Kunga, (d) Kharaulakh, (e) Nyaybinskaya; 2—seismodislocations; 3—locations of seismogravity effects; 4—rock fracturing diagrams, positions of vectors of the principal stress axes and fault planes; 5—focal mechanisms, dates and magnitudes of earthquakes (lower hemisphere; the principal stress axes of compression and extension are marked by black and white dots, respectively); 6–8—kinematics of active faults: 6—strike-slip fault, 7—thrust, 8—normal fault; 9—earthquake epicenters sized by magnitude: ≥6.8, 6.7–5.0, 4.9–4.0, ≤3.9. Systems of active faults: I—Primorskaya, II—West Verkhoyansk, III—Kharaulakh, IV—Buor–Khaya.

#### 3.2.2. Seismotectonics of Kharaulakh Segment

In our study, we use the data collected from seismic and geostructural surveys of the Kharaulakh Segment in combination with the information on tectonic stress fields reconstructed by tectonophysical analysis of the Late Cenozoic faults and folds in the study area. The available geological and geophysical datasets are also used to identify the systems of regional and local faults that were active in the Cenozoic. Their kinematics are confirmed by corresponding fracturing diagrams and earthquake focal mechanisms (Figure 4). In the zones of their dynamic influence, we identify seismogenerating structures of various sizes, which correlate with seismic events of M ≥ 7.5. Four main groups of fault systems are distinguished with respect to spatial locations, lengths, and kinematics: (I) Primorskaya (normal and strike-slip faults), (II) West Verkhoyansk (thrusts), (III) Kharaulakh (normal and strike-slip faults), and (IV) Buor–Khaya (normal faults).

The most active strike-slip faults are located in the central zone of the Kharaulakh Segment, which is traced as a system of closely spaced subparallel longitudinal faults arranged 'en echelon' to each other (Figure 4). In the map, this zone coincides with a field of minimum gravity values, which extends far to the south beyond the study area [7,21]. In the aerial photographs, morphological features show that the Kharaulakh faults are active in the Quaternary; numerous troughs, grabens, landslides, and avalanches are related to these faults. Furthermore, the available seismic records give evidence of active faulting in the Kharaulakh Segment. In 1927–1928, five strong earthquakes (М<sup>s</sup> 5.8–7.0) took place near the village of Bulun (Yakutia) at the southern termination of the fault zone in the same crustal block.

In the southern (most active) flank of this area, the sublongitudinal faults cross the western slopes of the Kharaulakh Ridge and run parallel to its axial line for a distance of 15 km. In the aerial photograph (Figure 5), a straight fault line displaces numerous river channels and cuts the watershed ridges.

**Figure 5.** Seismodislocation Beris (after [21]). (**a**) Photograph of surface ruptures resulting from the Bulun earthquakes of 1927–1928 (M ≥7.0); (**b**) aerial photograph of the Beris fault scarp, arrows point to the fault; (**c**) fragment of the deciphered image (1—strike-slip fault; 2—watershed axis); (**d**) photograph of seismogenic extension ruptures in the Kharaulakh fault zone (original location of the Beris fault).

According to the field structural and geological observations in the zone of dynamic influence and the rock fracturing analysis, this is a right-lateral strike-slip fault with a normal component. These fault kinematics are confirmed by the focal mechanism solution of the Bulun earthquake (1927; М<sup>s</sup> 7.0). The dip and strike of the fault plane are consistent in the reconstructions based on both the structural and seismic data. Recent activation features are clearly observed in the 'diagonal link' zone, which includes more than 20 paleo and modern seismodislocations caused by gravity and tectonic events. Some of the seismodislocations are marked by outcropped fault segments with horizontal displacement of 5–7 m (Figure 5).

Another seismically active area is the Buor–Khaya normal fault zone along the western coast of the Buor–Khaya Bay (Figure 4), which includes normal faults observed in the coastal outcrops from the Bykov Channel to the Kharaulakh depression for a distance of more than 160 km. The faults cut the basement of the rift structure, and many of them penetrate into the upper horizons of the sedimentary cover, which reveals their young age (Pliocene-Quaternary) [7,21]. The same age is determined from the field geological and geomorphological observation data [22] and multichannel seismic profiling [23]. Some of the faults are traced on land and to the Buor–Khaya Bay bottom and are clearly reflected in the sea bottom relief [11,12]. Their activity is evidenced by morphological features detected in the satellite images, as well as by their relationship to earthquakes (М<sup>s</sup> 3.5–5.4) and seismodislocations.

Considering the structural dynamics of the Kharaulakh Segment, we distinguish several areas that differ in types of the crustal stress state (Figure 4). The zones of strong earthquake epicenters differ in seismic parameters from one area to another, as shown by the data from the field geological and structural observations and the seismotectonic crustal deformation regimes reconstructed from the seismic data. According to our analysis of the state of crustal stresses, the Kharaulakh Segment is a unique transition region where extension is replaced with compression [7,21]. Similar regions on the globe are the Afar Rift in East Africa and the Northern California region. In the Kharaulakh Segment, stresses are mainly concentrated in the systems of sublongitudinal strike-slip faults with a normal component, which compose a block-concentrator that has a maximum seismic potential. Faulting in this area was influenced by the zone of the left-lateral strike-slip displacements with a normal component in the Olenek sector of the Lena–Anabar regional suture, as well as by the dynamically conjugated northwestern system of faults traced from the Chersky mountains (Figures 1 and 3). This conclusion is confirmed by the general sublongitudinal strike of the tectonic structures in the Kharaulakh Segment. The flanks of the dome structures show the northwestward deviation. In this area, the Tuora–Sis Range (Figure 4) extends to the Lena River left bank (Chekanov Range). Other facts in support of the above findings are the Cenozoic basins bordered by the strike-slip faults with a normal component from the west and the thrust faults from the east. The structural and seismological data show that the Kharaulakh duplex, represented in the west by ramp anticlines of the Tuora–Sis Range, is the main seismogenerating structure in the Kharaulakh Segment [21].

#### *3.3. Segment of the Chersky Seismotectonic Zone*

#### 3.3.1. Structure and Tectonics of the Segment

The Chersky Seismotectonic Zone belongs to a collage of terranes, which differ in structural and geological evolution features, comprising the Verkhoyansk fold–thrust belt [24]. In the Middle Jurassic, most of these terranes amalgamated through various episodes of subduction, collision and extension into a single larger tectonic block termed the Kolyma–Omolon Superterrane. The southeastern segment of the Chersky Seismotectonic Zone constitutes a continental fragment of the Okhotsk crustal plate. Its basement is composed of Archaean and Early Proterozoic crystalline rocks overlain with continental volcanic rocks of the Okhotsk–Chukotka belt [24,25]. The structure of the Chersky Seismotectonic Zone is complex and includes numerous linear fold zones and crustal blocks characterized by a high degree of seismic activity [26]. In this zone, the Yana–Indigirka and Indigirka–Kolyma regional sectors

are recognized by various features, including specific parageneses of active structures, anomalies of potential fields, and seismicity patterns (Figure 6).

#### 3.3.2. Seismotectonics of the Segment

Structural, morphotectonic, and seismic data were analyzed to investigate the kinematic features of active faults in the Chersky Zone (Figure 6). The detailed analysis shows that most of the active faults were formed due to horizontal compression and crustal shortening. This conclusion is supported by the structural observations of movements in the fault zones, field studies of tectonic fracturing, intense dislocations of Cenozoic sediments, as well as geological maps [26]. Compression is also confirmed by the focal mechanism solutions of strong earthquakes recorded in the zones of influence of active faults.

**Figure 6.** Seismicity and active faults of the Chersky Seismotectonic Zone (modified after [26]). 1—Balagan–Tas volcano; 2—focal mechanisms, dates and magnitudes of earthquakes (lower hemisphere; the principal stress axes of compression and extension are marked by black and white dots, respectively); 3–5—kinematics of active faults: 3—strike-slip fault, 4—thrust and reverse faults, 5—normal fault; 6—inferred fault; 7—earthquake epicenters sized by magnitude: ≤3.0, 3.5, 4.0, 4.5, 5.0, 5.5, 6.0, ≥6.5. Active faults: 1—Lena–Anabar regional suture, 2—Verkhoyansk marginal suture, 3—Omoloi, 4—Ege–Khaya, 5—Nelkan–Kyllakh regional suture, 6—Burkhala, 7—East Sette–Daban, 8—Yudoma, 9—Bilyakchan, 10—Ketanda, 11—Nyut–Ulbei, 12—Yana, 13—Adycha–Taryn, 14—Ulakhan, 15—Darpir, 16—Inyali–Debin, 17—Chai–Yureya, 18—Bryungade, 19—Ilin–Tas, 20—Arga–Tas, 21—Myatis, 22—Yarhodon, 23—Korkodon, 24—Chelomdzha–Yamsk, 25—Polousnensky, 26—Yarkan (South Anyui), 27—Khetachan. Lithospheric plates: EU—Eurasian, NA—North American, OK—Okhotsk. Sectors: A—Yana–Indigirka, B—Indigirka–Kolyma.

The Ulakhan strike-slip fault system is the largest in the NE regions of Russia (Figures 6 and 7). It has a complete set of characteristic of similar fault systems located in other regions of the world: a major shear suture represented by a straight-line fault, a system of echelon ruptures, an asymmetrical geometric pattern of folds and fractures feathering the major suture, an intrusion belt at a distance from the major suture, and a wide stratigraphic range of rocks [24]. Specific structural features of this fault system are an echelon series of left-lateral strike-slip faults and a chain of young pull-apart mini-troughs that formed in the crustal stretching segments of this fault. The Ulakhan fault crosses the Rassokha and Omulevka rivers. During the Middle Pleistocene–Holocene, the Ulakhan fault shifted the river channels to the left for about 24 km [26,27], which shows that an average rate of horizontal tectonic movements along the fault amounted to 5–7 mm/yr.

**Figure 7.** Structural dynamics of the Omulevka block. 1—active strike-slip fault; 2—earthquake epicenter; 3—focal mechanisms, dates and magnitudes of earthquakes (lower hemisphere; the principal stress axes of compression and extension are marked by black and white dots, respectively); 4—seismodislocation Eryun–Tas–Takh. Arrows show directions of plate movements. Lithospheric plates: NA—North American, OK—Okhotsk.

In the Ulakhan fault zone, many traces of ancient and modern seismodislocations are observed, including the most spectacular one—an almost 50 m high landslide dam in the upper reaches of the Tirekhtyakh River. This natural dam occurred due to a catastrophic earthquake hundreds or thousand years ago [28] (Figure 8).

**Figure 8.** Dam on the Eryun–Tas–Takh River (after [25,26]). (**A**) front view (the headwater segment the Eryun–Tas–Takh River is blocked by the dam); (**B**) top view (the river valley and the detachment wall are blocked by collapsed rocks); (**C**) transverse profile across the seismogenic dam (after [28]): 1—coarse-grained material: (a) granite gneiss, (b) carbonates, 2—fine-grained sandstone material, 3—erosion line of the rock collapse, 4—line before the rock collapse, 5—reconstructed line of the rock collapse wall, 6—lower boundary of the dam lake, 7—hypsometric level before the rock collapse, 8—observation points, 9—absolute elevation; (**D**) rock fracturing diagram for the active fault zone: 1—active fault plane, 2—fracture density isolines.

The Ulakhan fault conjugates with the Darpir fault (Figure 7). This pair of active conjugated faults reflects the style of tectonic deformation at the SE boundary between the North American and Okhotsk plates. Both faults are clearly detectable on satellite images. In topographic maps, these faults are shown as lengthy lineaments of the NW strike, which converge at an angle of 20–25◦. The faults border the uplifted Omulevsky block of Paleozoic rocks, which is elevated relative to the Mesozoic structures by 450–550 m. The Omulevsky block is viewed as a terrane that formed in the Mesozoic structural frame during collisional and post-collisional transformations of the Verkhoyansk–Kolyma belt. The background seismicity in the Omulevsky block suggests that this rootless terrane [24] was detached from the neighboring rocks as a result of horizontal displacements.

The Chai–Yureya fault is a significant structure in the Chersky system of active left-lateral strike-slip faults (hereafter the Chersky fault system), which is traced from the Okhotsk Sea coast to the Indigirka River as a series of separate fault segments typical of left-lateral strike-slip faults (Figures 6 and 9).

**Figure 9.** Geological model of the 1971 Artyk earthquake epicenter zone [29]. 1—modern alluvial deposits; 2—Middle Upper Quaternary water–glacial deposits; 3—Neogene–Quaternary deposits of the Verkhnener basins; 4—Jurassic sediments; 5—Triassic sediments; 6—lines of confirmed and inferred tectonic faults; 7—Artyk earthquake epicenter; 8—area of maximum values of ground surface deformation; 9—focal mechanism, dates and magnitudes of earthquakes (lower hemisphere; the principal stress axes of compression and extension are marked by black and white dots, respectively).

The entire length of this fault is detected as a zone of intense rock dislocations, which is marked by a linear magnetic anomaly and a sharp gravity gradient jump [29]. The Artyk earthquake (1971, М<sup>s</sup> 6.6), which took place in the zone of the Chai–Yureya fault, was one of the strongest earthquakes in NE Asia. During 12 months after the main shock, more than 1500 aftershocks were recorded in a 60 km long area extending in the NW direction from the epicenter of the main seismic event along the Chai–Yureya fault. The distribution of the aftershocks with depth and focal solution of the main shock show strike-slip movements on the nearly vertical fault plane (Figure 9) [2].

In the north–west, the Chersky fault system is conjugated with the Polousnensky thrust zone of sublatitudinal strike. The latter is represented by a series of subparallel thrust and reverse faults with planes dipping to the south and south–west (Figure 6). The focal mechanism of the Irgichan earthquake (1962, М<sup>s</sup> 6.2) also shows thrusting in the zone of the Polousnensky faults.

The Adycha–Taryn zone of active thrust faults borders the entire Chersky fault system from the south–west. This zone includes a series of 'en echelon' faults shifted by straight transverse strike-slip faults. Thrusting is confirmed by the focal mechanism solutions of the earthquakes (1951, М<sup>s</sup> 5.2–6.4) that occurred near the Upper Adycha depression [29]. Several seismodislocations are located in the area where the Triassic sandstones are thrust onto the Holocene terrace. In the middle course of the Adycha River, the Triassic rocks are thrust onto the Middle Pleistocene and Middle Quaternary rocks. The planes of all the above-mentioned thrust faults dip to the northeast underneath the Chersky Ridge. The northeastern boundary of the Chersky Seismotectonic Zone is the Myatiss thrust fault that is traced over 700 km along the NE foot of the Moma Ridge in the zone of its junction with the Indigirka–Zyryansky trough (Figure 6). At the NW flank of this fault, there is the Andrey–Tas area of maximum seismicity, as evidenced by the Uyanda (1984, М<sup>s</sup> 5.6), Andrey–Tas (2008, М<sup>s</sup> 6.1), and Ilin–Tas (Abyi) (2013, М<sup>s</sup> 6.9) earthquakes, as well as numerous minor seismic events (Figure 10). The main shock of the Ilin–Tas earthquake and its aftershocks occurred in the duplex compressional structure defined by the tectonic deformation data. Compression is confirmed by the focal mechanism solutions [30,31].

**Figure 10.** Strong earthquakes in the Andrei–Tas area of maximum seismic activity values (after [30,31]). 1—earthquake epicenters sized by magnitude: 3.0–3.9, 4.0–4.9, 5.0–5.9, 6.0–6.9; 2—focal mechanisms, dates and magnitudes of earthquakes (lower hemisphere; the principal stress axes of compression and extension are marked by black and white dots, respectively); 3—Andrei-Tas earthquake epicenter of 22 June 2008 (inset **A**) and Ilin–Tas (Abyi) earthquake epicenter of 14 February 2013 (inset **B**); 4—zones of seismogravity effects The insets show seismotectonic deformation structures resulting from the earthquakes.

The Andrey–Tas block bordered with thrust faults is a compact rectangular structure experiencing intense uplift (up to 2500 m high) and subjected to faulting of various ranks and kinematic types. Shearing is observed mainly in its central part, and a fan-shaped series of thrust faults develops in the NW flank. In the images from Roskosmos Geoportal, seismotectonic dislocations are detectable in the epicenter areas of the Andrey–Tas and Ylin–Tas earthquakes (Figure 10). The strike and morphological features suggest reverse faulting and shearing along this fault. The features detected from the seismotectonic and seismogravitational effects are possibly related to recent seismic events recorded in the Andrei–Tas block and the paleoearthquakes. It is noted that the areas with seismotectonic deformations and the epicenter zones of the seismic events are spatially coincident, and the shapes of the ground surface dislocations are consistent with the movements reconstructed from the corresponding focal mechanism solutions.

The Chersky fault system is sublongitudinally adjacent to the Ketandino–Ulbey system of right-lateral strike-slip faults (Figure 6). In the space images and aerial photographs, the latter are detected from the segments of straight fault lines that disturb the watershed parts of the ridges. Displaced granite intrusions and geological boundaries along these faults suggest that the right-lateral shear displacements amounted to 20 km [26,29]. This type of shearing is also confirmed by the focal mechanism solutions of strong earthquakes that occurred in the zones of dynamic influence of the faults. The Chelomdzha–Yamsky fault traced near the Okhotsk Sea coast is a boundary between the outer and inner zones of the Okhotsk–Chukotka volcanogenic belt (Figure 6). The magnetic and gravity fields are detected from magnetic anomalies and a sharp gravity gradient jump, respectively [29]. The fault is accompanied by a chain of intermountain depressions filled with Neogene and Quaternary sediments. The fault kinematics are defined by straightness and the presence of folded and thrust dislocations in the Cenozoic sediments of the depressions influenced by the fault, which is a thrust with a left-lateral strike-slip component.

Seismic activity increases at the sutures bordering the Chersky Seismotectonic Zone. The Sette–Daban earthquake (1951, М<sup>s</sup> 6.5) occurred in the Kyllakh block, which is the area of dynamic influence of the Nelkano–Kyllakh marginal suture (Figures 6 and 11). In this area, the deformation field is complex due to overlapping of the seismogenerating structures that belong to the Arctic–Asian, Okhotsk–Chukotka, and Baikal–Stanovoy seismic belts. The structure is characterized by lystric thrust faults that steeply dip near the surface and continue eastwards as low-angle detachments at depth. In the western part of the Nelkano–Kyllakh area, the sediments are thrust onto the subhorizontal beds of the Jurassic and Cretaceous deposits of the Siberian platform [24,25]. Due to thrusting, the eastern wing of the anticline is overlain with the Vend–Cambrian deformed beds of the East Sette–Daban tectonic zone. The E–NE orientated compression is evidenced by two fault planes reconstructed from the seismological data for the Sette–Daban earthquake, the NE plane (right-lateral strike-slip), and the NW plane (left-lateral strike-slip) [27,32]. In general, the identified seismotectonic deformation parameters reflect the trends in the distribution of crustal stress fields reconstructed from the structural data.

In order to investigate the crustal stress state of the Chersky Seismotectonic Zone, we analyzed the focal mechanism solutions of earthquakes recorded in the frontal zone, where the Kolyma–Omolon superterrane interacts with the Eurasian plate (Figure 12). These seismic events took place in conditions of steady NE-oriented compression. The principal compressional axis is subhorizontal (dip from 3◦ to 44◦) and acts across the strike of the structural elements in this zone (Figure 6). The extension stress axis is often coincident with the fault strike and oriented either horizontally or subvertically (dip from 2◦ to 85◦). The spatial orientations of the intermediate stress axis are chaotic, and the dip angles vary in a wide range (from 0◦ to 82◦). The above-described stress orientations reconstructed from the focal solutions are dominant across the entire Chersky zone.

**Figure 11.** Structural dynamics of the Kyllakh block, South Verkhoyansk sector (modified after [27,32]). 1–7—ages of deposits: 1—Cretaceous, 2—Jurassic, 3—Carboniferous–Permian, 4—Ordovician–Silurian–Devonian, 5—Vendian–Cambrian, 6—Middle–Upper Riphean, 7—Lower Riphean; 8–10—kinematics of active faults: 8—thrust; 9—strike-slip; 10 – normal fault; 11—syncline axis; 12—anticline axis; 13—Kyllakh earthquake epicenter; 14—focal mechanism, dates and magnitudes of earthquakes (lower hemisphere; the principal stress axes of compression and extension are marked by black and white dots, respectively). Faults: B—Burkhala, K—Kyllakh.

According to the kinematics of the principal tectonic stresses inferred from the seismic data, most of the earthquake focal mechanisms in the Chersky Seismotectonic Zone yield reverse (40%), strike-slip (30%), and thrust (20%) solutions, and about 10% show a combination of strike-slip and normal faulting (the percentage is calculated for 24 solutions; see Figure 6). These data clearly demonstrate that in the zone of interaction between the Kolyma–Omolon superterrane and the Eurasian plate, the seismic process develops under compression along the system of major faults that are conjugated with the marginal thrust and reverse faults.

Based on the structural pattern of the major seismogenerating zones and their dynamics, we proposed a model for the entire Chersky Seismotectonic Zone (the inset in Figure 12). This model shows transpression settings initiated by the interacting frontal structures in the contact zone between the Eurasian and North American lithospheric plates moving relative to each other at different rates [7,29]. Such a setting is possible if the Kolyma–Omolon block (in the frontal part of the North American plate) operated as an active indenter during convergence of the lithospheric plates in the NE direction. Due to its impact, a fan-shaped set of NW-trending left-lateral and SE-trending right-lateral strike-slip faults was formed. Seismogenerating zones of reverse faulting and thrusting, which developed at the terminations of the above-mentioned faults, have high seismic potential. Specific features of the

modern geodynamics of the Chersky Seismotectonic Zone are reflected in the regular pattern of the fields of local earthquake epicenters.

**Figure 12.** Scheme showing contemporary dynamics and earthquake epicenters in Yana–Indigirka segment of Chersky seismotectonic zone (after [26]). 1–2—kinematics of active faults: 1—strike-slip, 2—thrusts and reverse; 3—directions of movements of blocks; 4—focal mechanisms of earthquakes (lower hemisphere; the principal stress axes of compression and extension are marked by black and white dots, respectively); 5—Balagan–Tas volcano; 6—earthquake epicenters sized by magnitude: ≤3.0, 3.5, 4.0, 4.5, 5.0, 5.5, 6.0, ≥6.5.

Due to compression, the continental part of the Okhotsk plate is shifted in the E–SE direction, which contributes to the left-lateral shearing along the Ulakhan faults and right-lateral displacements in the Ketandino–Ulbey zone (Figure 6). This conclusion is supported by the development of the Pliocene–Quaternary extension zones along the northern coast of the Okhotsk Sea, which resulted from shifting of the Okhotsk Sea plate to the southeast. The E–SE thrusting of the continental part of the Okhotsk Seaplate is confirmed by the topographic and geomorphological observations that show uplifting in the mountainous area along the Okhotsk Sea coast. The watershed line between the Okhotsk Sea and the Arctic Ocean is disproportionately close to the Okhotsk Sea. In this region, shifting to the southeast relative to the Eurasian plate at a rate of 2–4 cm/yr is detected by satellite

geodetic monitoring of the reference points installed by the US Geological Survey in Magadan and Petropavlovsk–Kamchatsky [26,29]. The model showing the structure and geodynamics of the Chersky Seismotectonic Zone allows us to discover the regularities in the occurrence of active structures of various types and identify the blocks that act as tectonic stress concentrators that have high seismic potential (М<sup>s</sup> ≥ 6.5).

#### **4. Discussion and Conclusions**

The territory of the Arctic–Asian seismic belt is influenced by a variety of geodynamic processes: spreading (Gakkel Ridge), rifting (Laptev Sea shelf), transtension (Kharaulakh Zone), and transpression (Chersky Seismotectonic Zone). The fault–block structures differ by kinematic types in the belt segments but reflect the regularities in the development of the major seismogenerating structures in the entire Arctic–Asian seismic belt.

Spreading takes place in the Gakkel Ridge that represents the Arctic segment of the boundary between the Eurasian and North American lithospheric plates. In this area, the dominating seismotectonic setting is crustal stretching, and the principal stress axes are sublatitudinal across the strike of the main structural elements. South of the continental slope is the part on the Laptev sea microplate (Figure 2). The seismotectonic deformation parameters in this area show transverse subhorizontal extension. According to the seismological data, compression takes place at the western and eastern margins of the Laptev Sea microplate, as well as in the Lena River delta, in response to rifting in the Laptev Sea shelf. The Olenek sector (a structural border between the Laptev Sea rift system and the Siberian platform) is also subjected to compression, with a small left-lateral shear component. Zones differing in the crustal stress state (stretching, compression, and a variety of stretching–compression combinations) are observed in the Kharaulakh Segment, where the mid-oceanic and continental structures are conjugated. Southeast of the transtension area, in the Chersky Seismotectonic Zone, the North American and Eurasian plates approach each other at an oblique angle. In the Chersky zone, the steady NE–trending compression / transpression is confirmed by the structural–tectonic and seismological data and evidenced by the seismogenerating structures of thrust, reverse and shear types.

A change in the geodynamic setting from tension to compression can be satisfactorily explained in the global plate tectonics theory, considering that an assumed pole of rotation of the North American and Eurasian plates is located near the Buor–Khaya Bay [21,29]. It can, therefore, be concluded that stretching (extension) is currently typical of the neotectonic structures located to the north of the rotation pole, while compression takes place in the structures located to the south and southeast of the rotation pole.

The geodynamic processes dominating in the regional segments of the Arctic–Asian seismic belt predetermine the types of focal zones of strong earthquakes. The transtension zones, where the mid-oceanic and continental crustal structures are conjugated, demonstrate the highest seismic potential. In the transpression zones, increased seismic potential is typical for the blocks that acts as tectonic stress concentrators. Such blocks are formed at the terminations of the wings of active shear systems and associated with 'pull-apart' and/or duplex compression structures. In remote images, these blocks are clearly detectable from the geomorphological and morphodynamic features of the modern terrain. The flanks of the Chersky zone, which are confined to the reactivated structures of the marginal sutures, are also characterized by an increased level of seismic activity. Strong earthquake focal zones are formed in accordance with the dynamics of the major seismogenerating structures located in the central zone of the Arctic–Asian seismic belt.

Our results of the seismogeodynamic analysis maybe used for assessments of potential seismic risks of seismogenerating, and potentially active, structures in the Russian Arctic region, which are located at the lithospheric plate boundaries between the Eurasian, North American and Okhotsk Sea Plates.

**Author Contributions:** L.P.I. carried out the preparation and analysis of the structural-geological and morphodynamic situation in the seismic zone of Chersky, and she also had the idea of determining the geodynamic activity of neotectonic structures. V.S.I. conducted field work on the study of the pleistoseist regions of strong earthquakes of the Chersky seismic belt. He was also responsible for determining the kinematics of active faults and the paleoseismological situation in the seismic zone. B.M.K. was engaged in catalogues of earthquake epicenters and focal mechanisms of earthquakes in the Chersky belt, as well as materials of seismological data.

**Funding:** This study was supported by the grant of the Russian Foundation for Basic Research (RFBR) No. 19-05-00062. The study was carried out at the financial support of RFBR in the framework of the scientific project No. 16-05-000224, a project of the Institute of the earth's crust Siberian branch of RAS (No. 346-2018-0001) and the Institute of diamond Geology and noble me-for metal Siberian branch of RAS (No. 0381-2616-0001) and programs of the Government of republic Yakutia for study of region in 2016-2020-glare of Sakha (Yakutia) on complex study territory of the Republic for 2016-2020.

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

MDPI St. Alban-Anlage 66 4052 Basel Switzerland Tel. +41 61 683 77 34 Fax +41 61 302 89 18 www.mdpi.com

*Geosciences* Editorial Office E-mail: geosciences@mdpi.com www.mdpi.com/journal/geosciences

MDPI St. Alban-Anlage 66 4052 Basel Switzerland

Tel: +41 61 683 77 34 Fax: +41 61 302 89 18

www.mdpi.com

ISBN 978-3-0365-1878-7