Seismicity, focal mechanism, and stress tensor analysis of the Simav region, western Turkey

1 Introduction

The Hellenic and Cyprian trenches are major structures that have dominated the active tectonics of western Turkey. The western Anatolian plate shows EW trending horst and graben structures [1]. These structures and crustal thickening in western Anatolia occurred as a result of rapid extension and magmatism, which followed the closure and the destruction of the Neotethyan oceans [2]. Western Turkey is dominated by the northward subduction of this oceanic lithosphere between early Cenozoic and late Eocene [3]. Accordingly, the crust has been stretched [4] and thinned from 45–50 to 28–32 km [5,6]. In the western Anatolia, extension started at late Oligocene and still exists [4,7]. The observed crustal thickness of the region is between 30 and 34 km [8].

The Simav region is located in the Akşehir–Simav fault system, which is one of the most active fault systems in western Anatolia with a dominant normal oblique structure. The Simav graben has a WNW-ESE trend located between the Alaçam and the Eğrigöz mountains at the north and the Simav Mountain at the south. It is approximately 50 km long and has a ∼10 km wide EW trending basin and ∼1,100 m elevation above the sea level [9]. There are active faults located in and around the Simav region. The faults were mapped as a part of a project by Emre et al. [10] (Figure 1). The Simav Fault Zone, with its right-lateral strike-slip mechanism, marks the western edge of the Simav Graben showing NE-SW normal trending fault mechanism with average dip of 45°–50° [9,10]. The Naşa Fault Zone (NFZ) is another active fault system, trending NW-SE and 20 km in length and 5 km wide [10] (Figure 1).

Figure 1

Active major fault zones and topography map of the study area. SFZ, Simav Fault Zone; NFZ, Naşa Fault Zone; EGFZ, Emet-Gediz Fault Zone. Blue star symbols show epicenter location of moderate earthquakes [10].

Before the May 19, 2011, mainshock (M_w 5.9), there were three moderate (magnitude greater than 6.0) earthquakes that occurred in the study area in the twentieth century: May 2, 1928 (M 6.2) earthquake [11], June 25, 1944 (M 6.0) earthquake [11], and March 28, 1970 (M 7.1) Gediz earthquake [12] (Figure 1). The distribution of seismicity cannot be associated with the system of faults in the region, which has been of interest to many researchers in geology, geochemistry, and paleoseismicity [9,13,14,15,16,17,18,19]. A detailed study in western Anatolia, published by Akyol et al. [20], in 2006 observed lower than average crustal seismic velocities that may be associated with high crustal temperature. They also observed seismicity to a depth of about 10 km. The research of İnan et al. [17] concentrated on precursory anomalies before earthquake occurrence, by monitoring hot spring water chemistry changes and radon gas release in western Turkey. No precursory anomaly was observed in this region before the May 19, 2011, 20:15 UTC Simav earthquake. They indicated that this may be a result of pre-earthquake strain in the Simav block not transferred to the adjacent blocks. Yolsal-Çelikbilen et al. [21] indicated a centroid depth of about 9 km. Demirci et al. [18] did a study on rupture analysis, moment tensor (MT) inversion, damage distribution, and failure stress analysis. They observed a 70- to 80-degree dipping WNW-ESE normal fault, after the May 19, 2011, Simav mainshock, as the result of the active extension regime in western Turkey. Karasözen et al. [22] reassessed the previous moderate earthquakes occurred in the past half-century and imaged the May 19, 2011, Simav earthquake by using the interferometric synthetic aperture radar. They derived a north-dipping fault segment at about 7–9 km depth with little or no slip in the upper few kilometers by using seismic waveform and InSAR modeling of May 19, 2011, Simav earthquake.

This article reports on the results obtained from the analysis of seismicity and stress regime of the Simav region. There is no clear alignment on the rupture plane observed after the relocation of the earthquake sequence after the May 19, 2011, Simav mainshock. Local stress field and focal mechanism characteristics are consistent with the tectonic regime, crustal deformation, and seismicity pattern of the study area.

2 Seismicity and earthquake relocation

The moderate-size (M_w 5.9) Simav earthquake occurred on May 19, 2011, at 20:15:23.00 UTC time (mentioned as “mainshock” hereafter) at 8 km depth in southern edge of the NFZ. A Mw 3.6 earthquake occurred almost at the same location only 15 min before the mainshock. Significant aftershock activity lasted over 6 months with more than 5,000 earthquakes with M > 0.5. The largest aftershock (M_w 4.9) occurred 9 days after the mainshock. On May 3, 2012, another moderate earthquake with M_w 5.0 occurred to the northeast of the mainshock. A total number of 21,508 earthquakes were recorded by both KOERI (Boğaziçi University Kandilli Observatory and Earthquake Research Institute – Regional Earthquake-Tsunami Monitoring Center) [23] and TU (National Seismic Network of Turkey (DDA)) [24] networks in the region (latitude: 39.7°–38.5°N and longitude: 28.30°–30.0°E) between January 2005 and August 2018. These catalogs were merged by author and earthquake distribution and statistics of the region analyzed by using the ZMAP 7.1 software [25]. The seismicity rate is calculated for background seismicity between January 2005 and mainshock and for seismic distribution after mainshock till the end of August 2018 (Figure 2).

Figure 2

Seismicity rate for background seismicity and for the period after mainshock.

Evolution of the seismic network after 2008 in Turkey [26] and temporary networks located after Simav mainshock also played a significant role in increasing the detection rate of small-magnitude earthquakes. The majority of the background seismicity change between 2.5 and 2.9 magnitudes (Figure 2). As a result of increasing number of local seismic stations and temporary networks (such as TU, KOERI, Çanakkale Onsekiz Mart University (COMU), and Tübitak Marmara Research Center (MRC)), magnitude detection threshold decreased (Figure 2). Magnitude distribution of the 21,508 earthquakes is shown in Figure 3.

Figure 3

Magnitude distribution of the 21,508 earthquakes.

One of the critical parameters for seismicity studies is magnitude of completeness of the earthquake catalog (M_c), which is also critical for the determination of the parameter b. In the earthquakes detected between January 2005 and May 19, 2011, a number of 5,197 events were used to calculate M_c and the b value. Increase in the distribution of background seismicity is represented in Figure 4b with respect to cumulative moment release (Figure 4a). From the frequency–magnitude distribution of earthquakes (Figure 4b), the magnitude of completeness was determined as 2.7 and the resultant b-value is 1.3.

Figure 4

(a) Cumulative moment release determined from the background (1.1.2005–19.5.2011 until Simav mainshock) seismicity, plotted as a function of time. (b) Frequency–magnitude distribution of background seismicity (1.1.2005–19.5.2011 until Simav mainshock) with magnitude completeness of M_c = 2.7.

Earthquakes with available waveform data (both P and S phase arrivals) are repicked manually by using the PQLX code [27] and are relocated by using the HYPOINVERSE Earthquake Location code [28] in order to minimize combined catalog data uncertainties. Body wave phase readings are used to apply double-difference algorithm, compute a 1-D crustal model, and for the focal mechanism analysis in the study area. The results of HYPOINVERSE were used as starting locations for a new and more accurate earthquake relocation using the double-difference algorithm HypoDD [29]. The station distribution is quite dense with the total number of 86 seismic stations (Figure 5). All waveform data from three-component broadband seismometers belonging to permanent and temporary networks are processed. Waveform data from seismic stations at epicentral distances closer than 50 km are used to reduce location residuals for earthquakes with magnitudes smaller than 2.0. A total number of 10,727 earthquakes from KOERI and TU catalogs was revised as mentioned above. (This catalog is mentioned as “revised catalog” hereafter.) Earthquake locations of the revised catalog with respect to earthquake magnitudes are shown in Figure 6.

Figure 5

Seismic stations distribution operated by different networks. RETMC: Kandilli Observatory and Earthquake Research Institute Regional Earthquake-Tsunami Monitoring Center (21 stations), ERD: Prime Ministry Disaster and Emergency Management Presidency, Department of Earthquake (24 stations), MRC: Tübitak Marmara Research Center (38 stations), COMU: Çanakkale Onsekiz Mart University Geophysics Department (2 stations). Study area is framed as orange.

Figure 6

Map shows the raw locations of events before revision (gray) and double-difference (blue) relocated events of the revised catalog. The P-wave velocity model (this study) depicted at the right-bottom corner (depths in km, velocity in km/s). The NS and EW cross-sections show the depth distribution of the seismicity. Red star represents the mainshock. Black lines refer to the active faults [10]. Earthquakes with magnitudes greater than 0.7 are shown. Most of the earthquakes are relocated at depths shallower than ∼18 km.

Determining earthquake locations and origin times with minimum errors and maximum accuracy is highly dependent on the velocity structure of the study area. A regional velocity model for the study area is constructed using the VELEST inversion code [30]. Earthquakes located at less than 150 km epicentral distances with minimum of 15 P-phase and 3 S-phase readings are selected from the revised catalog in order to construct the velocity model. The earthquakes recorded by less than 15 seismic stations and having an azimuthal gap greater than 60 degrees are eliminated. Earthquakes with location errors greater than 2 km both at latitude and longitude and earthquakes with depth errors greater than 5 km were ignored. These section criteria were chosen in order to construct the 1-D velocity model with minimum location errors. The velocity model, computed by Akyol et al. [20], was chosen as an initial model for the relocation process. A total number of 515 earthquakes, with 9,981 P-phases and 1,605 S-phases high-quality picks were relocated using the initial velocity model, as an input to the iterative 1-D inversion algorithm. Station corrections were performed with respect to a reference seismic station. The Wadati diagram showed quality data, and the resultant V_p/V_s ratio was 1.77 (Figure 7). As a result of iterations, a 1-D velocity model that performed minimum RMS (root mean square) errors and best trade-off with station corrections was selected. The Moho discontinuity was observed at 33 km depth. The final model is similar to the velocity model of Karasözen et al. [22], obtained by a similar methodology but using different data sets (Figure 8). The RMS error is decreased from 0.493 to 0.197 at the end of iterations for the revised catalog, after relocation with the new 1-D velocity model. Earthquake locations and depths differ from the revised catalog with initial velocity model in minor changes at latitude (∼0.096 km), longitude (∼0.12 km), and depth (∼5 km). The seismicity in the study area is mainly observed within ∼18 km depths (Figure 6).

Figure 7

Wadati diagram for revised catalog. The V_p/V_s ratio computed as 1.773.

Figure 8

Velocity models from different studies (initial model is depicted as green line) and a model resulted from the final inversion (red line).

The improved earthquake locations show that the seismicity mainly occurred around the mainshock location, with a small clustering in the Eğriöz Mountain (Figure 6). A dense activity was located in northwestern Tavşanlı, and a sparse activity took place along the EGFZ associated with the NW-SE trending structure of this fault zone. On the other hand, the aftershock activity in the study area is not clearly associated with the main Naşa fault system or the rest of the faults in this region. As observed, they do not correspond to any geomorphological structure. In other words, the aftershock activity cannot be unambiguously related to a tectonic structure or to any of active fault systems in the study area (Figure 6). Cumulative number of earthquakes dramatically increases after the Simav mainshock. No time-dependent earthquake clustering is observed, which might be evidence of fluid flow and dynamic triggering in the study area. Neither swarm nor depth-dependent activity migration is observed (Figures 9a–c).

Figure 9

The cross-section of revised catalog earthquakes with latitude (°) (a), longitude (°) (b), and depth (km) (c) as a function of time.

3 Regional focal mechanisms

More than 80 earthquakes with magnitude greater than 3.5 occurred in the region, between January 2005 and September 2018; however, due to poor azimuthal coverage and the low signal to noise ratio (SNR) values, many of these events were not used for focal mechanism computations. MT solutions of 59 earthquakes were calculated to determine the kinematics of the earthquake sequence and stress field orientation. For these calculations, the ISOLA software (ISOLated Asperities) [31] was used, which is a regional centroid MT inversion algorithm proposed by Kikuchi and Kanamori [32]. MT solutions were calculated as point sources. The time and position of point sources were obtained by the grid search. The discrete wavenumber method [33,34] was used to calculate the Green’s functions. A 1-D velocity model (mentioned as this study in Figure 8) was used for the calculations of the Green’s functions. The velocity model plays a significant role in estimating the Green’s functions. It is important to select stations with good azimuthal coverage rather than the number of stations. For each MT solution, minimum of five stations were used. Seismic data are filtered in different frequency bands, between 0.06 Hz (±0.02) and 0.18 Hz (±0.02) for different sizes of earthquakes to obtain an optimum solution [31]. Stations with clipped waveform were eliminated from the database.

The ISOLA code provides a number of parameters in order to assess the solution quality. The variance reduction (VR) is a measure of a residual between the observed and synthetic waveforms. Other indicators are the condition number (CN: ratio of maximum to minimum singular values of Green’s functions matrix G), the SNR, the double-couple (DC) percentage, the focal mechanism variability (FMVAR) index, and the space-time variability (STVAR) index [35]. The MT solutions with VR > 0.5, CN < 10, FMVAR < 30, and STVAR < 0.3 are acceptable according to the ISOLA user’s guide [31]. Five earthquake focal mechanism solutions, which did not fulfill the acceptability criteria, are eliminated (listed in Table S1 as supplementary). The MT solutions are obtained for 54 earthquakes as plotted in Figure 10. The generalized linear inversion is used to calculate full MTs and decomposed into DC and non-DC components. VR values range between 0.5 and 0.7.

Figure 10

The focal mechanisms determined in this study. Beachball symbols numbered for each MT solution up to 54. The number 8 (red beachball) represents the May 19, 2011, Simav mainshock solution. The dimension of circles represents magnitude of events. Black lines represent active faults in the study area [10] (see supplementary Table S1).

In this study, the MT of the May 19, 2011, Simav mainshock is determined as normal faulting for auxiliary planes with strike, dip, and rake angles of 265°, 68°, −109° (selected as fault plane) and 127°, 29°, −52°, respectively. The mean of dip direction for a total of 54 MT solutions is 59.4°. Figure 10 shows that fault-plane solutions dominantly exhibit normal faulting (−140° < rake < −90°).

A ternary diagram shows a distribution of types of focal mechanisms (Figure 11). The majority of the solutions show a normal faulting motion, which is consistent with the splay’s movement of the study area and an extensional regime of the western Aegean. Also a few oblique fault solutions with strike-slip components are observed but no thrust faulting was detected.

Figure 11

Rake-based ternary diagram of 54 MT solutions. No thrust fault mechanism is observed. Most of the solutions belong to the normal faulting mechanism (MATLAB code by M. Hallo [36] is used to plot the diagram).

4 Stress inversion analysis

Many techniques have been proposed to determine tectonic stress orientation and stress ratio (R) from focal mechanisms of earthquakes [37]. One of the most common techniques is developed by Michael [38]. This method determines the orientation of the three-principle stress axes and aims to find the best-fitting stress tensor to the observed focal mechanisms. This is a linear inversion method and defining of the fault plane is not required. On the other hand, this method performs a bootstrap resampling to choose a fault plane randomly from nodal planes. This method assumes a uniform and invariant stress field in space and time, and slip occurs in the direction of maximum shear stress [39].

The MT solutions of 54 earthquakes were inverted for principal stress directions and for the shape ratio. The MATLAB public-open STRESSINVERSE code by Vavryčuk [40] was used to analyze the statistical properties of focal mechanisms with respect to tectonic stress and failure conditions. The iterative stress inversion can be easily implemented with this code, which is based on the Micheal’s method [38] modified by Vavryčuk [40].

Stress inversion results are represented by the three principal stress axes. The pressure (P) and tension (T) axes are shown (Figure 12a). The principal stress directions (Figure 12b) are as follows (azimuth/plunge): σ1 = 205°/79°, σ2 = 100°/3°, σ3 = 10°/11°. 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.56. It describes the relative magnitudes of the principal stresses (R = (σ1 − σ2)/(σ1 − σ3)) [41]. The inversion was run repeatedly, using a grid in searching through the principal stress orientations, and with the friction coefficient ranging from 0.4 to 1, in 0.05 increment. The optimum friction coefficient computed is 0.7. The stress inversion accuracy varies depending on the number of focal mechanisms and noise level of data. Six iterations, with 1,000 realizations of noise, were performed. Average fault normal and slip direction errors calculated as ±6°.

Figure 12

(a) Distribution of P–T (red-blue, respectively) axes and principal stress axes (green), (b) confidence regions of principal stress axes, shown in the stereographic projection.

5 Conclusion and discussion

According to the dense seismic station coverage, seismicity of the study area was observed in detail between 2004 and 2018. Statistical analysis performed in the area with a combined catalog of earthquakes outlines characteristic seismicity features. The seismicity after the May 19, 2011, Simav mainshock was not clearly compatible with the main Simav fault, but it’s compatible with NFZ.

The Simav region has a morphology shaped with the EW trending high-angle normal faults. In this study, strike, dip, and rake angles for Simav mainshock are determined as 265°, 68°, and −109°, respectively. In the Global Centroid Moment Tensor catalog [42], the Simav mainshock MT solution values are listed as 286°, 46°, 85° and 98°, 44°, −96° as strike, dip and rake, respectively. Normal faulting mechanisms with a small strike-slip component as a result of teleseismic inversion after May 19, 2011, Simav mainshock in the region were observed by Yolsal-Çelikbilen et al. [21]. They determined strike, dip, and rake angles of the fault plane as 287°, 58°, and −94°, respectively, from P-wave first-motion polarities. Karasözen et al. [22] published a study on the Simav region. By using teleseismic body waves and regional MT analysis, they observed MT inversion solution for the mainshock fault plane as 289°, 54°, and −93° as strike, dip, and rake angles, respectively. They observed normal faulting at 7–9 km beneath surface without any visible rupture at the surface consistent with the extensional regime of the region. As listed in Karasözen et al. [22], different observations with various techniques consistently indicate WNW-ESE normal faulting in the study area with dip greater than 50°. Although the dip direction of the mainshock in this study is larger than in previous studies, the average dip angle is consistent with the NNE-SSW regional trending extensional regime in the region. Demirci et al. [18] determined the stress tensor in the Simav region by using the inversion code developed by Carey [43] and Carey-Gailhardis and Mercier [44], with 23 focal mechanisms of earthquakes with M > 3.5. They observed the N12°E extensional stress regime, which is coherent with the NNE-SSW regional trending extensional regime and the average shape ratio as 0.8.

The Simav region is one of the major metamorphic regions in the western Turkey, with high-angle north-dipping normal mechanism in NNE-SSW extension direction. The weakness zone in the western Anatolia is substantially in the NS direction [45,46], because the dominant characteristic of the region is extension in the NS direction. The stress inversion results in this study reveal north-dipping mechanism with a general trend of the NNE-SSW extension for the Simav region and dominant normal fault characterization in the WNW-ESE direction is compatible with the orientation of σ₁ observed in this study. As a further topic, it would be desirable to study the strain rate in order to know the velocity gradient and rates of displacements in the Simav region.