Seismotectonics of Malatya Fault, Eastern Turkey

1 Introduction

Precise assessment of seismic risk requires not only the understanding of kinematics and structure of major plate boundary zones but intraplate deformation zones as well. In intraplate settings, rate of strain accumulation is relatively low and the associated fault rupture occurs over longer recurrence intervals (e.g. several hundreds to thousands of years; [1, 2]) than moderate or large earthquakes occurring on plate boundaries. However, the abrupt release of accumulated stress over years may lead to destructive earthquakes (e.g.Mexico, Bulnay [3, 4]).Moreover, major plate boundary zones are monitored by quite dense sensor networks (seismic, GPS-Global Positioning Systems-and strong motion), whereas intraplate settings, like the region of interest in this study, attract less attention even though they possess a high seismic risk.

Eastern Mediterranean is an unprecedented region with intriguing tectonic processes taking place: a continental collision between stable Eurasia and northward migrating Arabian Plate along Bitlis-Zagros-Suture Zone, another collision between stable Eurasia and Nubian Plate along Hellenic-Cyprus Arc and the lateral escape of Anatolian Plate as a consequence. To compensate these motions, Anatolian plate is accommodated by two main plate boundaries; dextral strike-slip North Anatolian Fault Zone (NAFZ), sinistral strike-slip East Anatolian Fault Zone (EAFZ) (Figure 1a). Both fault zones are mapped during the 20th century and studied extensively since then [5, 6, 7, 8, 9, 10, 11, 12, 13]. Many devastating earthquakes occurred on these faults resulting in several numbers of causalities and damage. In Figure 1b, seismicity between years 1900-2019 in a greater area, can be seen. Main target region of this study, the Malatya Fault, is located to the west of EAFZ covering east of the Anatolian Plate. According to studies based on geodetic sensors, Anatolian Plate moves westward with an insignificant amount of internal deformation -less than 2mm/year- [14, 15, 16, 17] whereas geological studies and seismic activity suggest the opposite, that Anatolian Plate accommodates a considerable amount of deformation [18, 19, 20, 21, 22, 23, 24, 25, 26]. The eastern half of the Anatolian Plate is dominated by several active, dextral to sinistral strike-slip faults and related pull-apart basins [15]. Several studies exist in order to identify seismic hazard related to the intraplate deformation in eastern Anatolian Plate, however, studies related to Malatya Ovacık Fault Zone is limited to mainly geological, morphochronological and paleoseismological studies [17, 27, 28, 29].

Figure 1

a) Tectonic map of the study area. Main fault zones in the region are; NAFZ: North Anatolian Fault Zone, EAFZ: East Anatolian Fault Zone, CAFZ: Central Anatolian Fault Zone, MOFZ: Malatya Ovacık Fault Zone. Blue inverted triangles represent seismic sensors and blue squares represent strong motion sensors operated by AFAD. Seismic sensors operated by KOERI are shown as green triangles [35]. Black lines represent the active faults [30]. b) Seismicity map of the study area. Earthquakes with magnitude above 3.0, between the years 1900 and 2019 are plotted. (Data provided from KOERI).

The 275 km long NE striking sinistral Malatya Ovacık Fault Zone (MOFZ) is one of the most important structures in Anatolian Plate since it is responsible for the initiation and deformation within the intraplate [27, 28, 29]. The northern branch, Ovacık Fault (OF) starts from Erzincan Basin and extends to Arapgir (Malatya) of about 110 km. It branches into several segments. The southern branch, NE striking, approx. 165 km-long Malatya Fault starts from Adatepe forming a junction with OF and terminates at Doğanşehir (Figure 1a). Malatya Fault is described in three segments in [30], namely Akçadağ, Arguvan and Doğanşehir segments, whereas [29] defined Malatya Fault in five segments, naming FS1 to FS5. Several authors [28, 29], suggest that Malatya Fault is an active fault partly accommodating the intraplate deformation of Anatolia on the contrary to the hypotheses of Malatya Fault being inactive since mid-Pliocene [31, 32]. According to recent paleoseismological studies, the average earthquake interval on Malatya Fault is 2275 ± 605 years and it may produce a M7.4 earthquake following the empirical equations of [28, 33]. Moreover, another study on the seismic potential of Malatya Fault reports that the fault passes through Quaternary units in some of parts of Akçadağ Segment; and thus, is capable to produce large earthquakes [34].

In this study, we analyze datasets obtained from seismic and strong motion sensor networks in order to identify the seismotectonics of Malatya Fault (MF). We estimate a 1-D crustal velocity model for the region. We have improved accuracy of earthquake hypocenter locations by merging data from two local networks and using all possible sensors in the region. Hypocenter locations of low-magnitude-earthquakes cannot be determined due to limitations in station distribution, therefore our results are constrained. Nevertheless, obtained focal mechanisms and stress directions are in good agreement with the tectonics of region. Moreover, we identify an active morphological component, which is not present in the active fault map of Turkey. Thus, this study has a significant contribution to the seismic hazard assessment of a region prone to large earthquakes.

2 Data

Turkey, a region exhibiting an intense seismic activity (Figure 1b), is continuously monitored with seismic, geodetic and strong motions sensors. Several M+7 earthquakes occurred in the region; therefore, a dense sensor network is necessary to investigate such a complex and seismically active landscape in order to better understand active tectonics.Moreover, precisely determined earthquake parameters are important not only for understanding the seismicity and seismotectonics of a region but also for providing an accurate estimation of seismic hazard and risk parameters. Local agencies; Kandilli Observatory and Earthquake Research Institute (KOERI) [35] and Disaster and Emergency Management Authority Presidential of Earthquake Department (AFAD) operate seismic and strong motion sensor networks (Figure 1a). Both agencies are utilizing similar type of seismic sensors mostly recording at 100 sps. In this study, we employ seismic and strong motion sensors mainly located in eastern Turkey. A total number of 123 stations are included in earthquake hypocenter determination. Even though the sensor networks seem to be robust, due to low signal to ratio waveforms or GPS clock errors, some stations are not employed in the analysis.

3 Estimation of crustal velocity of the region and surroundings

In order to determine a minimum 1-D crustal velocity model for the study area, an earthquake catalog is constructed. By merging the earthquake catalogs of KOERI and AFAD between 2008-2018, approximately 1000 earthquakes have been compiled and revised for this study in the region (37.10^∘–39.10^∘ E, 37.30^∘–39.10^∘ N). Merged catalogs are revised and relocated by using the earthquake location code HYPOCENTER [36].

A minimum 1-D velocity model can be defined as the velocity model that produces the smallest possible uniform error for the set of events with well-constrained locations [37]. For the identification of an optimum 1-D velocity model, a widely acknowledged software package VELEST code [38, 39] is used. The method simultaneously inverts P and S wave travel times for a 1-D velocity model, station delays and hypocenters of earthquakes. Simply, it solves the forward problem by tracing direct reflected and refracted rays from source to receiver and then uses the standard damped least squares in order to solve the inverse problem [37].

Generally, events are selected according to the following criteria;

• The number of P and S phases should be more than 8 with at least one S phase,

• The azimuthal gap should be less than 180^∘,

• • The RMS residual should be less than 0.5,

• The distance to nearest station is less than one focal depth.

After eliminating the earthquakes according to the given criteria, 495 earthquakes are selected for inversion. In the computation of final velocity model, nine iterations performed for each test by trying 100 different initial models. We also test a number of different starting 1-D velocity models [40, 41]. Out of several iterations, using different initial velocity models, a final velocity model chosen with the minimum rms value reduced to 0.35 (Table 1).

4 Earthquake Hypocenter Determination

Precise hypocenter locations of earthquakes are essential in understanding the tectonic regime in a region. In this study, waveforms recorded on both seismic and strong motion sensors operated by two local networks; KOERI and AFAD are taken into account. P- and S-wave arrival times of a dataset of 146 events are re-picked manually. At this stage, re-picking of data is performed by one experienced person in order to provide consistency. We included a re-picked data set of events with minimum 40 P-wave (1/3 ratio S-wave) observations per event. Integrating different kinds of sensors from two different local networks significantly reduced the azimuthal gap to 20-40 degrees, which is an important parameter to provide accuracy. Even though the azimuthal gap is ensured to be as low as possible to obtain accurate hypocenter locations, still another important parameter, focal depth, needs to be constrained. Several rules of thumb exist to obtain reliable depth estimates. At least one P-wave phase is required at a station within a close proximity of the source, approximately 1.5 times the focal depth. In this study, the minimum distance between the source and station-a P-wave phase can be picked- is approximately 20 km, which is one broadband seismic sensor. For the rest of the sensors, minimum distance to source is 50 km. To determine hypocenter locations, we applied NonLinLoc (Non Linear Localization), which is an absolute probabilistic localization algorithm [42]. This algorithm follows the probabilistic inversion approach defined by [43], and the earthquake location methods of [44, 45]. Detailed information on NonLin-Loc can be found in [42].

As it can be seen in Figure 2, merging data of seismic and strong motion sensors from two networks, significantly improved the hypocenter locations and the depth resolution. Improved seismic activity focuses along Malatya Fault (MF), mainly the southern portion of MF, around Kızık Basin and near Hekimhan.

Figure 2

Earthquake hypocenter locations and depth cross-sections. Green dots represent hypocenter locations from AFAD catalog, red dots represent hypocenter locations from KOERI catalog. Black dots represent the hypocenter locations re-located by combining the entire data from KOERI and AFAD seismic networks.

5 Determination of Focal Mechanisms

Focal mechanism solutions are determined using a merged dataset of waveforms from AFAD and KOERI networks. The number of earthquakes employed at this stage reduced significantly due to low signal to noise ratio of waveforms and limitations of station distribution. Waveform data for 57 earthquakes (M>3.2) are selected. Focal mechanisms are calculated by using an exhaustively tested and widely used ISOLA software package [46, 47, 48]. This technique utilizes an iterative deconvolution inversion [49] to solve for the best single- or multiple-point source representation of each earthquake. Moment tensors are estimated by a least-squares minimization of mis-fits between observed and synthetic waveforms; subevent positions and relative times are calculated through grid search [47]. In order to represent a full wavefield at regional and local distances, Green’s functions are calculated using the discrete wavenumber method of [50] and [51]. Employing an accurate velocity model in computation of the Green’s function plays a crucial role in inversion process. Significance of velocity model and the corresponding variations are analyzed in several studies [52, 53, 54, 55]. Hence, our first step was to compute a site-specific model for the target region (Section 2, Table 1). Green’s functions are estimated utilizing this local velocity model.

In Figure 3a,Moment Tensor (MT) solution of the M4.9 Hekimhan-Malatya earthquake and in Figure 3b observed versus synthetic waveforms can be seen.Waveform data of 17 stations from AFAD and KOERI network are merged in order to improve station distribution. Presented results exhibit a variance reduction of 0.91 and a consistent shallow depth at 6.7 km (hypocentral depth) and 10.0 km (centroid depth).

Figure 3

a) Moment Tensor Solutions of the M4.9 Hekimhan-Malatya earthquake and b) the observed versus synthetic waveform fits.

MT solutions obtained for the target region are presented in Figure 4 (Appendix A). Moreover, the results are plotted as a ternary diagram [56]. In Figure 5, the distribution of computed focal mechanisms can be seen in a triangular representation. At the vertices of the triangle, mechanisms with vertical T, P and B axes are plotted for thrust, normal and strike-slip faulting, respectively.

Figure 4

Focal mechanisms of 57 earthquakes.

Figure 5

Ternary diagram, a graphical representation of focal mechanism orientations. The vertices of the triangle represent mechanisms with vertical T, P and B axes. Focal mechanism solutions obtained in this study cluster around strike-slip and normal faulting.

6 Stress Inversion Analysis

To provide a more reliable determination of the direction of released stress during the earthquake, we applied a stress inversion method by using focal mechanism solutions obtained in the previous section. The most frequently used technique is defined by [57]. This method determines the orientation of three principle stress axes and aims to find the best fit between the stress tensor and observed focal mechanisms. In this linear inversion method, nodal plane knowledge is not required. This approach performs a bootstrap re-sampling to choose fault plane randomly from nodal planes. The method assumes a uniform and invariant stress field in space and time, and slip occurs in the direction of maximum shear stress [58]. In this study, Michael’s method is preferred due to its appropriate estimate of the uncertainty and accuracy of results even with low signal to noise ratio data.

The MT solutions are inverted to obtain the principal stress directions and shape ratio. A STRESSINVERSE code by [59] is used to analyze the statistical properties of focal mechanisms with respect to tectonic stress and failure conditions. The iterative stress inversion based on Michael’s method [57] can be easily implemented with this code [59]. Stress inversion results are represented by three principal stress axes; a "o" for σ1, a “x” symbol for σ2, and a “+” symbol for σ3 (Figure 7a). The pressure (P) and tension (T) axes are well clustered (Figure 7b). The azimuth angle is measured from north and the plunge angle from the horizontal plane. The shape ratio (or stress magnitude ratio) (R) is 0.24 and describes the relative magnitudes of the principal stresses, (R = (σ1−σ2)/(σ1−σ3)) [60] (Figure 5). Inversion is performed in an iterative mode, by using grid search through the principal stress and with friction ranging from 0.4 to 1, in 0.05 increments. Accuracy of the stress inversion depends on the number of focal mechanisms and signal to noise ratio of data.

Figure 6

Stress inversion solutions. (a) Distribution of P-T axes orientations of individual nodal planes, (b) P and T axes shown in stereographic projection.

Figure 7

Focal mechanisms of earthquakes with magnitude above 3.8 are plotted within the shear zone identified with respect to the seismic activity focused near Hekimhan. Solid red line indicates the corresponding source fault of the observed seismicity.

The orientations of P/ T axes are nearly vertical and nearly horizontal (Figure 5), respectively, suggesting an extensional regime. In addition to this, nearly subhorizontal σ3 and nearly subvertical σ1 indicate the existence of an ongoing extensional regime within the study area. This may imply the existence of an ancient fragmented fault in the region. There are few earthquakes with pure strike slip mechanisms however our results indicate strike-slip mechanism with a normal component. Our results coincide with [18], which suggests northeast-southwest shortening and corresponding northwest-southeast extension for the Central “ova” neotectonic province of Anatolia. Nevertheless, more comprehensive analysis needs to be performed with enhanced instrumental and paleoseismological studies in the vicinity of the fault.

7 Discussion and Conclusion

In this study, seismotectonics of the Malatya Fault, eastern Turkey, is investigated by using a unified dataset of local agencies, KOERI and AFAD. A local crustal velocity model is estimated to be further employed in the analysis. Earthquakes in the study area are relocated in order to obtain more precise hypocenter parameters. Moment tensor solutions are computed utilizing these improved hypocenter locations. Focal mechanism solutions and stress inversion results of this study depict a well recovered image of the active tectonic structure in the target region.

Malatya Fault, the southern branch of MOFZ, is argued to be the plate boundary between Anatolian and Arabian Plate until the development of NAFZ and EAFZ in mid-Pliocene [31] and labeled as inactive since then. In a recent study, [28] found evidence of four earthquakes from pre-instrumental period according to their paleoseismological observations in the region with a recurrence interval of 2275+-605 years. They stated that another earthquake with a surface rupture is due. These findings are in good agreement with geological and morphological evidence [20, 26, 29] highlighting the Malatya Fault as an active fault, partly accumulating intraplate deformation of Anatolian Plate [17]. According to geodetic studies within the region, most of the strain accumulation due to westward movement of Anatolia is accommodated by EAFZ rather than Malatya Fault [17].

In the light of these arguments, we analyze seismic activity in the region between the years 2008-2019. Seismicity mainly focuses in two zones; one along the northern branch of Malatya Fault, and the second one near Hekimhan, almost oblique to Malatya Fault. The highest magnitude events are M5.1 occurred on 29.11.2015 (near Hekimhan) and another event, magnitude 4.8, occurred on 25.03.2019 (near Arguvan).We determine hypocenter locations of 146 events. In earthquake hypocenter determination, we merge datasets of seismic and strong motion sensor networks in order to reduce azimuthal gap, an important parameter to ensure accuracy. However, instrumental problems in sensors and low signal to noise ratio data limit our calculations. Focal mechanisms of only 57 events are calculated, reflecting normal and strike slip faulting. Relatively shallow depths obtained in hypocenter determinations, might be an evidence of insignificant amount of accumulated strain. In a recent study, it is reported that Tohoku-Oki earthquake rupture is terminated by the seismic activity in the adjacent segments [61]. Reliable geodetic surveys need to be conducted in this region in order to better understand whether Malatya Fault has accumulated any amount of strain to precisely determine the seismic risk.

The most striking result of this study is the identification of a NW-SE trending deformation zone, which becomes evident following the seismic activity near Hekimhan (Figure 2). In the interior parts of the Anatolian Block, all tectonic regimes including such a NW-SE trending one are suggested to be formed as a result of the activation of pre-weakness belts due to the westward motion of the Anatolian Block [18, 62]. Major shear zones such as the Malatya Ovacık Fault, are often associated with a structural pattern of faults arranged as en échelon [63, 64]. The observed seismicity may be related to an anti Riedel structure within the shear zone of the Malatya Fault (Figure 7). The solid red line in Figure 7 may be good candidate for the source fault of these seismic events, which also makes a lineament on the morphology. As a property of anti Riedel structures, a dextral strike-slip movement is expected [18], since the kinematics of the principal displacement zone is sinistral. According to the results of shear-box experiments of [63], the inclination of anti Riedel shears is approximately 78^∘. We determine an inclination of 55-60^∘ of the observed anti Riedel structure (Figure 7). There may be several reasons of this low inclination angle. Tchalenko [63] performed these experiments using homogeneous kaolin clay whereas the study region consists of highly heterogeneous distribution of geological units, as expected. Another reason may be the amount of data employed in this study. Limited number of earthquakes and related focal mechanisms may not be enough to address the state of stress in such a complex tectonic region.

Even though the target region is monitored with a high number of stations from two different networks, the results of this study are mainly constrained by the limitations of sensor networks along the Malatya Fault. Moreover, there is not any consensus on the active tectonic regime of the region. Therefore, it is necessary to increase the density of the station distribution along this fault and analyze the region in more detail. The results of this study underline the importance of deploying dense networks of different types of sensors to monitor earthquake prone regions as well as providing maintenance of these sensors regularly. Despite all limitations, defining an active morphological element, which is not included in recent active fault map of Turkey, prove to be noteworthy for seismic risk assessment of the target region. More detailed studies on the geological and morphological evidence of this deformation zone need to be conducted.