Seismicity, focal mechanisms and active stress ﬁeld around the central segment of the North Anatolian Fault in Turkey

We analysed locations and focal mechanisms of events with magnitude ≥ 3, which are recorded by 39 broad-band seismic stations deployed during the North Anatolian Passive Seismic Experiment (2005–2008) around central segment of the North Anatolian Fault (NAF). Using P - and S -arrival times, earthquakes are relocated and a new 1-D seismic velocity model of the region is derived. Relocated events in the area are mainly limited to a depth of 15 km and present seismicity in the southern block indicates widespread continental deformation. In the next step, focal mechanisms are derived from ﬁrst motions ( P, SH ) and amplitude ratios ( SH/P ) using a grid-search algorithm in an iterative scheme. Analysis of our well-constrained focal mechanisms indicate mainly strike-slip motions apart from some normal and few thrust events that are related to complex local fault geometry. Calculated pressure/tension axes are mainly subhorizontal and maximum horizontal stress directions ( SH max) are oriented predominantly in NW–SE direction which corresponds well with the slip character of NAF and its splays. In the east, E–W trending splays show right-lateral strike-slip mechanisms similar to the main strand whereas in the west, antithetic N–S trending faults show left lateral strike-slip motions. The seismic cluster that converged near C¸orum after relocation indicates a dominant right-lateral strike-slip mechanism along the E–W trending fault. These focal mechanisms are used to perform stress tensor inversion across the region to map out the stress ﬁeld in detail. Overall, maximum ( σ 1 ) and minimum ( σ 3 ) principal stresses are found to be subhorizontal and the intermediate principle stress ( σ 2 ) is vertically orientated, consistent with a dominant strike-slip regime. These directions point to the clockwise rotation of stress trajectories from N to S where NW–SE directed σ 1 in the north turns towards N–S in the south away from the NAF. Moreover, the 200-km-long Ezinepazar–Sungurlu Fault which is previously mapped as an active strike-slip fault is characterized by minor seismic activity and trends perpendicular to the computed maximum stress direction in the southwest away from the main strand of NAF suggesting that the Sungurlu segment is either compressional in nature or inactive.


I N T RO D U C T I O N
The tectonic framework of the eastern Mediterranean is characterized by convergence of the northward moving Arabian and African plates with respect to the relatively stable Eurasian Plate (Dewey et al. 1986). In the east, continental collision between African and Eurasian plates led to a compressional regime, whereas in the west, rollback of subducting African lithosphere beneath Aegean Trench resulted in a backarc extensional regime (Bozkurt 2001). This gave palaeotectonic structures created by closure of Tethys ocean during the early Cenozoic. All indications suggest that the NAF is a newly coalescing continental transform plate boundary which forms the northern margin of the Anatolian Plate (i.e. Biryol et al. 2010).
After 1999 Izmit (M w = 7.4) and Düzce (M w = 7.1) earthquakes, many seismic studies were performed to obtain accurate earthquake locations and focal mechanisms along the western segments of NAF. Bohnhoff et al. (2006), Görgün et al. (2010) and Pınar et al. (2010) investigated aftershock focal mechanisms of Izmit earthquake to reveal fault segmentation and local stress conditions along the fault zone.Örgülü (2011) analysed microseismicity and focal mechanisms across the major basins in Marmara Sea to identify active stress field and style of deformation.
Although, its geological and geomorphologic features are well defined and extensively studied, crustal deformation and associated seismicity around central segment of the NAF is relatively less known. Since, this segment has an important role in detecting the seismic risk in Ankara (Fig. 2) and its surroundings, it is essential to characterize the seismicity and segmentation along the central section of NAF. In this study, earthquakes recorded by the North Anatolian Passive Experiment (2005)(2006)(2007)(2008) have been relocated using an average 1-D velocity model computed from the analysed seismic data. Their focal mechanisms are then determined and used to resolve the current local stress field. Finally, these results are correlated with the local geology to characterize the active tectonics of the region.

S E I S M O T E C T O N I C S E T T I N G
The main strand of NAF had a record of large devastating earthquakes in the last century (Stein et al. 1997). The central northwards convex bend of the NAF includes segments that ruptured during large earthquakes (M > 7) occurred in 1939, 1942, 1943and 1944). The measured long-term geological slip rate along the main strand of NAF is 20 ± 5 mm yr −1 (Kozacı et al. 2007) which is comparable to GPS studies indicating 24 mm yr −1 uniform slip rate along the NAF (McClusky et al. 2000;Reilinger et al. 2006). In particular, the overlap zone between the two largest events (1939 and 1943) is characterized by southward bifurcating right lateral splay faults extending further westwards into the interiors of the Anatolian Plate, which can be compared with synthetic Riedel shears. Main splays from SE to NW include the Central Anatolian Fault (Koçyigit & Beyhan 1998), the Almus Fault (Bozkurt & Koçyigit 1996), Ezinepazarı-Sungurlu Fault (Kaymakçı et al. 2003) and Laçin Fault (Kaymakçı 2000;Bozkurt 2001). These splays create wedge-like blocks rotating in a counterclockwise sense (Piper et al. 1986;Tatar et al. 1995;Taymaz et al. 2007). Unlike other E-W oriented splays, the Ezinepazarı-Sungurlu Fault changes its direction (36 • E) from E-W to NE-SW towards west and strikes trough the Anatolian Plate for more than 200 km (Şaroglu et al. 1992). During the 1939 Erzincan earthquake (M w = 7.9), the easternmost segment of the Ezine-Sungurlu Fault was ruptured but measured offsets along the fault decreases westwards as the fault turns southwards and becomes distant to the NAF (Erturaç & Tüysüz 2012). Recent GPS observations based on a local network in the region, also did not show any significant movement along the southwestern section of the Ezine-Sungurlu Fault (Yavaşoglu et al. 2011). Apart from subparallel splays, there are also two major approximately N-S trending faults, which are analogues to antithetic Riedel shears. The Dodurga Fault is a left-lateral strike-slip fault with considerable normal slip component (Koçyigit et al. 2000) whereas the Eldivan Fault possibly resulted from reactivation of a Neotethyan suture displays left-lateral strike-slip motion with minor reverse component (Kaymakçı 2000;Kaymakçı et al. 2000Kaymakçı et al. , 2003. Moreover, the northern tip of the NW-SE trending Tuz Gölü Fault, which is a right-lateral strike-slip fault with a mapped length of about 200 km (Dirik & Göncüoglu 1996;Koçyigit & Beyhan 1998) and its northward continuation, the Afşar Fault (Tan et al. 2010) lies southwest of the study area. Fig. 3 shows the distribution of instrumental seismicity (M ≥ 3) recorded between 1995 and 2009 by Kandilli Observatory and Earthquake Research Institute (KOERI). The observed seismicity is shallow and is not only restricted to the main strand of NAF. The widespread earthquake distribution in the southern block correlates well with the faults mapped in the region and provides evidence for ongoing internal deformation across the Anatolian Plate. During this time period, several moderate earthquakes (M > 5) were detected in the study area. Mescitözü earthquakes occurred sequentially during 1996 and 1997 along the right-lateral splay Laçin Fault with a possible triggering effect. In the west, the 2000 Orta earthquake (M w = 6) originated from the activation of the Dodurga Fault (Koçyigit et al. 2000) and formed a dense seismic cluster across Downloaded from https://academic.oup.com/gji/article-abstract/196/1/405/587190 by MIDDLE EASTTECHNICAL UNIVERSITY user on 24 July 2020  1939, 1942, 1943 and 1944. Faults are compiled and modified from Şaroglu et al. (1992), Koçyigit & Beyhan (1998), Kaymakçı (2000), Kaymakçı et al. (2000), Bozkurt & Koçyigit (1996) and Koçyigit (2009 Table 1. the fault. Most recently, the Bala earthquake sequence in 2005 resulted due to a conjugate strike-slip faulting (Koçyigit 2009) and was followed by the 2007 sequence, which occurred on the NW-SE trending Afşar Fault (Tan et al. 2010). In the region, there are a total of 41 available fault plane solutions of earthquakes recorded between 1938 and 2010 ( Fig. 4 and Table 1). The rake-based ternary diagram given as an inset in Fig. 4, clearly displays that most of these solutions have strike-slip mechanism in accordance with the NAF and its splays. Extensional mechanisms are mainly occurring along minor normal faults that exist in transtensional zones between strike-slip faults segments. The 2000 Orta earthquake (M w = 6) is the only major event displaying a significant normal slip component. The body-wave inversion results (Taymaz et al. 2007) however suggest larger component of left-lateral slip than the Harvard-CMT solution. Except the 1968 Bartın earthquake, which indicates almost pure thrusting, the interiors of the northern block are characterized by lack of major events. In the south, the focal mechanisms indicate active faulting mainly along the NW-SE trending right-lateral Afşar Fault.

E A RT H Q UA K E R E L O C AT I O N
In this study, earthquakes recorded by 39 broad-band seismic stations during the North Anatolian Passive Experiment (2005-2008Gans et al. 2009) are analysed. Initially processed data set included close to 200 events with magnitude ≥3. During the phase picking process, seismic records are rotated to radial and transverse components and together with the vertical component are picked using SAC software (Goldstein et al. 2003). P waves, which are a combination of P g and P n arrivals, are picked from the vertical components, and SH waves are picked from transverse. Earthquake relocation was carried out by using HYPOCENTER location algorithm that applies adaptively damped least squares (Lienert et al. 1986;Lienert & Havskov 1994). We used the Wadati Diagram in Fig. 5(a) in order to select the reference V p /V s ratio of 1.73. A total of 151 events with higher quality picks (3745 P-and 973 S-phase) are relocated using velocity models of Toksöz et al. (2003) and IASP91 (Kennett & Engdahl 1991) as initial trial models in the relocation process (Fig. 5b).
Earthquake relocation is considered to be a classical inverse problem in seismology. Locating an earthquake and finding its origin time depends highly on the velocity structure of the region since it controls the ray paths and hence, the arrival times of recorded seismic waves at various stations (Stein and Wysession 2003). In order to constrain a better velocity model for the region, relocated events which have a good azimuthal coverage with a maximum gap of 120 • , relatively low rms and at least 13 P-wave readings are selected. As a result, a total of 98 events with 1666 P-phase picks are input to the iterative 1-D inversion algorithm, VELEST. In the single event mode, this algorithm inverts for hypocentre locations and associated station corrections along with the optimum 1-D crustal velocity model (Kissling 1995). The station corrections are calculated with respect to an assigned reference station. It is worth to note that the inversion is limited to the first arriving P phase in order to minimize embedded picking errors that are likely more significant in S-phase picks.
The inversion processes are strongly dependent on the starting model and damping factors. Thus, a careful sensitivity analysis including numerous inversion tests is carried out by using different starting velocity models and various combinations of damping factors (Kissling 1995). After sensitivity analysis, two-stage inversion scheme, which does not allow presence of low-velocity zones as suggested by Kissling (1995), is chosen since it provided the most stable results for the velocity model and station corrections. In this scheme, earthquake location parameters of the relocated events are overdamped at both stages whereas damping of station corrections are kept higher than the velocity model in the first inversion to resolve the preferred velocity model better. During the second inversion, both velocity model and station corrections are damped equally to determine the required station corrections. In order to minimize the effect of the starting model, the first inversion is applied separately for two different starting models (Toksöz et al. 2003) and IASP91 (Kennett & Engdahl 1991) and their output velocity models along with an average model are then used during the second inversion. The inversion process was stopped when the station corrections and the velocity values did not vary significantly in subsequent iterations. At the end, all three resultant velocity models showed converging trends. The velocity model that gives the minimum rms misfit and best trade-off with station corrections is selected as the 1-D P-wave velocity model of the region together with associated station corrections. According to the obtained model, the crust has an average P-wave velocity of 5.7 km s −1 down to 10 km depth and 6.3-6.4 km s −1 from 10 to 35 km. Moho discontinuity is at 35 km depth and P-wave velocity of uppermost mantle is 7.7 km s −1 , which is slower than the global average ( Fig. 5b). Overall, our model is also consistent with the 1-D velocity model developed using similar methodology but slightly different data set for the local tomography study of the region (Yolsal-Ç evikbilen et al. 2012). Station corrections that are calculated with respect to a reference station reflect highly variable near-surface site conditions and partly account for spatial validity of the 1-D velocity model across the study area (Kissling 1995). The positive station corrections indicate true velocities slower than the model or vice versa. In this study, station KGAC located near the centre of the network (lat: 40.941, long: 34.323) with the highest number of records is selected as the reference station. The obtained final station corrections mostly vary between −0.3 and 0.3 s. Negative corrections indicating faster Finally, high-quality events are relocated by the HYPOCENTER program using the resultant 1-D P-wave velocity model and final station corrections which significantly reduced the overall rms error and provided a total of 152 well-relocated events (Fig. 6a). The improved earthquake locations show only minor variations in terms of latitude and longitude but rather significant changes in hypocentre depths, which is critical for focal mechanism analysis (Billings et al. 1994). According to our results, the studied events are mainly limited to the top 15 km and widely distributed suggesting a broad zone of deformation especially in the southern block. Moreover, relocation process revealed the existence of a cluster located near Ç orum. Close up view of this cluster clearly shows earthquakes that converged together after relocation (Fig. 6b).

F O C A L M E C H A N I S M
A focal mechanism solution describes the geometry and mechanism of the faulting during an earthquake. In this study, focal mechanisms of the relocated earthquakes are determined using the grid-search program FOCMEC (Snoke et al. 1984), which searches the focal sphere for potential solutions consistent with polarity and/or amplitude ratio data. Data picks includes polarity, impulsiveness and quality of P and SH waves. For accurate determination, the azimuth and take-off angles are computed by ray tracing program TauP (Crotwell et al. 1999) using the improved event locations and 1-D velocity model obtained during the relocation process.
There are multiple potential solutions that fit the polarity data of an event. To obtain more unique solutions, polarities of P and SH phases are used along with SH/P amplitude ratios since high ratios define a point located along the strike direction of a nodal plane or vice versa. The maximum allowed log10 ratio which is the maximum limit of the difference between theoretical and observed (Snoke 2003), is selected as 0.6 and values greater than this limit are assigned as amplitude error. In this step, an iterative grid-search scheme is applied by increasing the number of data errors to find solutions with minimum polarity and/or amplitude errors (Anderson et al. 2007). An example of our results obtained with and without using SH/P amplitude data are given in Fig. 7 along with P-and S-wave polarities plotted over the expected radiation pattern of the selected solution. Note that the input of amplitude ratios in addition to the polarities strongly limits the number of possible solutions and results in better constrained solutions.
In the presence of multiple solutions, the similarity of the solutions are also quantified with precision parameter κ of P-and T-axes, calculated using Fisher statistics (Fisher 1953). The high values of κ represent solutions with similar mechanisms and after visual inspection poorly resolved earthquakes with low κ values are excluded from further analysis. At the end, solutions with minimum polarity and amplitude ratio errors and lowest rms misfits are chosen as the best-fitting focal mechanisms, which led to a total of 109 well-constrained focal mechanism solutions (Table 2 and Fig. 8a).
Most earthquakes correlate well with the trend and slip character of NAF and its splays. Rake-based ternary diagram of resultant focal mechanisms (Fig. 8b) shows that apart from some normal and few thrust faulting, strike-slip motion is dominant in the region. In the west, NAF is characterized by right-lateral strike-slip mechanism with some thrust component which is in accordance with the arcshaped fault geometry and towards south observed focal mechanism are highly variable and relates to the N-S trending antithetic faults. In the east, overlapping right-lateral splays of NAF forms a zone growing wide towards west where some active normal faulting also occurs at the transtensional fault steps. Near Ç orum, a higher magnitude (M = 4.5) event is followed by many aftershocks forming a cluster that converged together after relocation. At this location, our seismic network has an excellent coverage and resulted in numerous reliable focal mechanisms indicating a dominant right-lateral strike-slip mechanism along an E-W trending fault (Fig. 8c).
Calculated pressure (P) and tension (T) axes orientations of our focal mechanisms and solutions from previous studies are shown together in Fig. 9. Most of the P-and T-axes are subhorizontal and indicate a strike-slip regime with NW-SE directed compression and NE-SW directed dilatation. In our earthquake database (Table 2) and the one compiled from previous studies (Table 1), the maximum horizontal compressive stress directions (SH max) of the individual focal mechanisms are calculated according to Lund & Townend (2007) and stress regime following the World Stress Map Project standard (available online at www.world-stress-map.org) which is first presented by Zoback (1992).
In order to compare and correlate with mapped faults, SH max orientations of individual events are plotted along with their type in Downloaded from https://academic.oup.com/gji/article-abstract/196/1/405/587190 by MIDDLE EASTTECHNICAL UNIVERSITY user on 24 July 2020  To investigate local tectonic variations in detail, the study area is subdivided based on their geographical location and tectonic setting. Fig. 11 shows the rake-based ternary diagrams of earthquakes occurred in each subregion along with calculated P-and T-axes orientations. Region 1, a narrow elongated area containing Gerede and Tosya segments of the NAF, displays a transpressional character. On the other hand, region 2 which includes Niksar-Erbaa pull-apart basin in the east indicates localized extension due to transtension. In the centre, region 3 (Amasya-Aydıncık) and region 4 (Ç orum-Sungurlu) are both dominated by strike-slip events with subhorizontal P-and T-axes. Region 5 located on the N-S trending antithetic faults (Dodurga-Eldivan) includes both strike-slip and normal earthquakes resulting in more vertical and diffuse P-axes cluster. In the south, events that lie within region 6 over the Afşar Fault (Bala-Afşar) are all strike-slip in character and trend of their P-axes indicate some rotation towards N-S.

S T R E S S T E N S O R I N V E R S I O N
In the presence of pre-existing faults, P-and T-axes could deviate significantly from the principal stress directions (Zoback et al. 1987). Therefore, numerous stress inversion techniques have been proposed for the determination of stress field from focal mechanism solutions (Angelier 1984;Gephart & Forsyth 1984;Michael 1984Michael , 1987. These techniques are developed to obtain the best-fitting stress tensor to the observed focal mechanisms, by minimizing the discrepancy between the shear stress and slip direction for all earthquakes in the data set. They have similar assumptions but vary on how they define the best model and the fault plane ambiguity (variance and misfit ;Angelier 1984;Michael 1987;Gephart 1990;Michael 1991;Bohnhoff et al. 2004;Görgün et al. 2010). Angelier (1984) developed a non-linear inversion technique which solves focal mechanisms by treating both nodal planes as fault planes (Michael 1987). This is an unphysical assumption because it can only work if the stress field is uniaxial. Another limitation of this technique is that, it would fail to determine the shape of the stress ellipsoid if the actual fault and auxiliary plane is not known. Gephart and Forsyth's (1984) and Michael's (1984Michael's ( , 1987 methods are widely used because they developed an algorithm which can correctly distinguish fault planes from the auxiliary plane during Downloaded from https://academic.oup.com/gji/article-abstract/196/1/405/587190 by MIDDLE EASTTECHNICAL UNIVERSITY user on 24 July 2020  the inversion (Michael 1987). Gephart and Forsyth's method (1984) attempts to find the stress tensor using a grid search over the stress field parameter space (Gephart & Forsyth 1984;Gephart 1990). Linear inversion method of Michael (1984Michael ( , 1987 solves the stress tensor using a linear least-squares method. This method performs a bootstrap approach where the fault plane is chosen randomly from the two nodal plane. All these methods assume that (1) the stress in the study is uniform and invariant in space and time, and (2) earthquake slip occurs in the direction of maximum shear stress (Delvaux & Barth 2010). If both conditions exist one may determine the directions of the principal stresses (σ 1 − σ 3 ) and a relative stress magnitude (R). The first assumption indicates a homogeneous stress field. For highly heterogeneous stress fields with a high average misfits, the inversion method may not give meaningful results. The second assumption implies isotropic fault planes. Because the real faults may not be isotropic, for the inversion process, additional random noise is used in the the fault set. Hardebeck & Hauksson (2001) compared these two methods in detail, applying them on synthetic focal mechanism data sets. Based on this comparison, both models can determine stress parameters accurately and increasing the data set can improve the accuracy of results. The method of Gephart & Forsyth (1984) results in more reliable stress orientations with high-quality data sets, whereas the method of Michael (1984Michael ( , 1987 provides a more appropriate estimate of the uncertainty and is more accurate with very noisy data sets (Hardebeck & Hauksson 2001).
In this study, we used the linearized stress inversion technique of Michael (1987) in the ZMAP software package (Wiemer 2001), which performs bootstrap resampling where the fault plane is chosen randomly from the two nodal planes. For each stress inversion, Downloaded from https://academic.oup.com/gji/article-abstract/196/1/405/587190 by MIDDLE EASTTECHNICAL UNIVERSITY user on 24 July 2020   (Table 2) and ones compiled from literature (Table 1). principal stress orientations (σ 1 , σ 2 , σ 3 ) where σ 1 > σ 2 > σ 3 and the stress shape ratio (Angelier 1979) are determined using 2000 bootstrap iterations, The shape parameter or stress magnitude ratio R describes the relative magnitude of the principal stresses and hence constrains the shape of deviatoric stress ellipsoid. Inversion of the entire data set indicates predominantly strike-slip deformation (Fig. 12). Both σ 1 and σ 3 are subhorizontal and trends N36 • W and N55 • E, respectively, whereas σ 2 is almost vertical. The resultant stress shape ratio of R = 0.62 suggests that the magnitude of σ 2 is close to the average of σ 1 and σ 3 defining a triaxial stress ellipsoid. The stress field heterogeneity is reflected by variance of the inversion and average misfit angle β between observed and predicted slip directions. Variance less than 0.2 (Wiemer et al. 2002) and β less than 33 • (Michael 1991;Görgün et al. 2010) is accepted as a relatively uniform stress field. High variance and β values are inconsistent with a single uniform stress tensor and hence imply spatially heterogeneous stress field within the analysed volume. In general, the resultant variance (0.18) and β (38 • ) are reasonable considering the size of the study area and errors present in focal mechanism solutions (Michael 1991) but some systematically larger values are observed locally reflecting the inherit spatial variations in the stress field. For proper detection of the stress variations present in the study area, stress tensors of the previously defined geographic regions are also inverted independently. The results are listed in Table 3 and shown in Fig. 13 along with the expected Riedel shear pattern of the region given as an inset figure.
Along Gerede and Tosya segments (region 1), the subhorizontal σ 1 rotates 25 • counterclockwise with respect to the stress field calculated for the entire area and trends N61 • W. The σ 2 becomes more horizontal while σ 3 moves towards vertical suggesting rightlateral strike-slip and some thrust faulting. Moreover, the low R value (0.2) calculated for the region indicates closer magnitudes for σ 2 and σ 3 and support local transpression with uniaxial compression in strike-slip stress regime (Michael et al. 1990). This resultant stress tensor has relatively higher variance and σ values indicating complex faulting and correlates well with the main strand of NAF and NE-SW trending thrust faults mapped in the south (Fig. 13). Across Niksar-Erbaa pull-apart basin, σ 1 trends N24 • W in accordance with the arc-shaped fault geometry of NAF and is closer to vertical indicating active normal faulting. The stress tensor of this region also displays the highest R value (0.76) and reveals localized uniaxial extension in a transtensional strike-slip stress regime. This is consistent with the orientation of normal faults located at the strike-slip fault steps (Fig. 13).
The central part away from the main strand of NAF, displays relatively uniform strike-slip stress regime with subhorizontal σ 2 and σ 3 . Region 3 (Amasya-Aydıncık) which includes splays of NAF, is characterized by N34 • W directed σ 1 . Region 4 (Ç orum-Sungurlu) contains 65 focal mechanisms solutions including the 'Ç orum cluster' (39 events). In order to eliminate effects of temporal stress associated possibly with the cluster, stress inversions are performed separately for all the data, only for the cluster and without the cluster (Table 3). The obtained results are comparable to each other with only minor differences. The trend of σ 1 is between N30 • W and N44 • W. In the absence of 'Ç orum cluster', the plunge of σ 1 is slightly higher implying minor normal faulting. For both regions (3 and 4), computed R values are similar (0.53-0.65) and define a triaxial stress ellipsoid. The resultant stress tensors indicate predominantly right-lateral strike-slip faulting along ∼E-W trending faults and compatible with most of the mapped faults. However, the western section of Ezinepazarı-Sungurlu Fault strikes in NE-SW direction, which is perpendicular to the computed trend of σ 1 (Fig. 13).
Downloaded from https://academic.oup.com/gji/article-abstract/196/1/405/587190 by MIDDLE EASTTECHNICAL UNIVERSITY user on 24 July 2020 Figure 11. Tectonic map of the study area showing rake-based ternary diagrams and stereographic P-and T-axes plots of earthquakes located in each subregion (dashed).

Figure 12.
Results of stress tensor inversion for the entire database. Stereographic projection of the resultant principal stresses with contours defining the 95 per cent confidence region is on the left, histogram of the stress shape ratio R is on the right. Across the antithetic faults (region 5: Dodurga-Eldivan), the inversion results are characterized by high variance and β values pointing to the heterogeneous stress field. Even though, mean σ 1 and σ 3 are subhorizontal, the confidence areas of both σ 1 and σ 3 are elongated along N30 • -40 • W with highly variable plunge amounts suggesting presence of both strike-slip and normal faulting in the region (Fig. 13). Relatively high R value (0.73) also reflects the uniaxial extension occurring in a strike-slip regime. This complex faulting can be also due to temporal stress variations associated with the 2000 Orta earthquake (M w = 6) which has an oblique strike-slip mechanism with significant normal slip component. In this respect, our results for this region are in agreement with mapped left-lateral strike-slip faults. It is also worth to note that a further subdivision can be applied but is avoided here considering the spatial clustering and limited number of events available in the region.
In the south, the events located mainly on the Afşar Fault (region 6: Bala-Afşar) and reveal a uniform stress field with low variance and β values. The moderate R value (0.43) is also an indicative of a triaxial stress ellipsoid. According to the inversion results, both σ 1 and σ 3 are horizontal and σ 2 is vertical, which Table 3. Results of stress tensor inversion applied for different regions. Trend (tr) and plunge (pl) of each stress axes (σ 1 , σ 2 , σ 3 ) is given along with the number of focal mechanisms, stress regime, stress shape ratio (R), average misfit angle (β) and variance.  suggest pure strike-slip faulting. The resolved trend of σ 1 is N12 • W and indicates a ∼20 • clockwise rotation in stress trajectories with respect to the north (Fig. 13).

D I S C U S S I O N A N D C O N C L U D I N G R E M A R K S
The central segment of NAF makes a northward convex bend and characterized at the surface by complex splays and Riedel shear structures. Throughout the study area, the seismogenic zone extend to a depth of about 15 km. The observed seismicity is not limited to a narrow band near the main strand of NAF rather widespread in the southern block suggesting a broad zone of internal deformation within the Anatolian Plate. Our 1-D velocity model of the region which indicates slow uppermost mantle velocities (7.7 km s −1 ) and on average 35-km-thick crust is within the range of previously conducted seismic refraction studies (Kuleli et al. 2001;Gürbüz et al. 2004). Apart from some normal and few thrust events related to local fault geometry, resultant focal mechanisms indicate mainly strikeslip motions in the region. Calculated pressure axes are mainly subhorizontal and maximum horizontal stress directions (SH max) are oriented predominantly in NW-SE direction. Similarly, the resultant stress tensor computed using all available earthquake data indicates a dominant strike-slip regime with a triaxial stress ellipsoid developed under NW-SE compression and NE-SW extension. According to our results, the near vicinity of NAF (between latitudes 40 • and 41 • ) is characterized by E-W trending right-lateral, N-S trending left-lateral and NW-SE trending normal active faults. The seismic cluster that converged near Ç orum after relocation also indicates a dominant right-lateral strike-slip mechanism along an E-W trending fault. On the other hand, south of latitude 40 • , active faulting occurs mainly along NW-SE trending right-lateral and NE-SW trending left-lateral strike-slip faults.
Based on stress inversions from slip vectors measured from fault planes, Bellier et al. (1997) proposed temporal stress change from transpression before middle Pleistocene to present-day transtension for this region. Separate stress analysis of regions defined according to their geographic location and tectonic setting in our study area revealed notable spatial variations in present stress state where arc shaped geometry of NAF and overlapping splays with extensional steps result in local transpressional (Gerede-Tosya) and transtensional (Niksar-Erbaa, Dodurga-Eldivan) stress regimes.
Since orientation of faults with respect to the stress trajectories defines the fault activity and type, combining the fault pattern with the present-day stress field is essential to gain further insight on the spatial distribution of earthquakes. In this respect, the stress trajectories for the maximum principle stress σ 1 have been manually reconstructed using stress orientations inferred from earthquake focal mechanisms at various locations that are connected, extrapolated and smoothened to represent the orientation of the stress field across the study area (Fig. 14). The reconstructed stress trajectories display clockwise rotation from N to S where NW-SE directed σ 1 in the north turns towards N-S in the south away from the NAF. The effect of this rotation is also evident on the fault orientations. The length-based rose diagrams constructed separately for the faults located north and south of 40 o N (inset in Fig. 14) reveal a prominent change in fault trends conformable with the detected rotation in the stress field.
The study area lies northeast of the capital city, Ankara. In the last century, the main strand of NAF along the central segment has ruptured in an eastward progressive fashion by a sequence of large earthquakes, including events in 1939 (M w = 7.9), 1942 (M w = 7.1), 1943 (M w = 7.4) and 1944 (M w = 7.3). The main strand which accommodates 90 ± 5 per cent of the slip rate (Yavaşoglu et al. 2011), is more than 100 km away from Ankara, but there are also other faults located closer to the city. Among them, Dodurga Fault recently activated during the 2000 Orta earthquake (M w = 6) and Afşar Fault during the 2005 and 2007 Bala earthquakes (M w ≤ 5.7). These events are felt strongly in Ankara located 50-80 km away from the events. Although inherited uncertainty related to the source parameters of the 1938 Kırşehir earthquake (M w = 6.7), which occurred before the instrumental era, is rather large, the event still stands as the largest magnitude event recorded close to Ankara (80 km away) and should also be taken into account for seismic hazard related studies. The 200-km-long Ezinepazar-Sungurlu Fault which is previously mapped as an active strike-slip fault is another major tectonic element that deserves close attention. The eastern Ezinepazarı segment which trends in E-W direction, is seismically active and ruptured during the 1939 (M w 7.9) event whereas in the west, the Sungurlu segment which strikes in NE-SW direction away from the main strand of NAF, displays minor seismic activity and trends almost perpendicular to the resolved maximum stress direction. Mature pre-existing strike-slip faults with low frictional strength and/or elevated pore pressure cannot sustain high shear stresses and thus can explain the maximum stress Downloaded from https://academic.oup.com/gji/article-abstract/196/1/405/587190 by MIDDLE EASTTECHNICAL UNIVERSITY user on 24 July 2020 observed at a high angle to the fault strike (e.g. Zoback et al. 1987). However, both GPS measurements and measured geological offsets along the Sungurlu segment show no significant displacement (Yavaşoglu et al. 2011;Erturaç & Tüysüz 2012) and the segment display minor seismic activity supporting that the fault is not an active weak strike-slip fault but either compressional in nature or inactive.

A C K N O W L E D G E M E N T S
This research was partly supported by NSF grant EAR0309838. The seismic instrumentation for the North Anatolian Fault Passive Seismic Experiment was provided by the Incorporated Research Institutions for Seismology (IRIS) through the PASSCAL Instrument Center. Most of the figures were created using the Generic Mapping Tools (GMT) software (Wessel & Smith 1995). The authors thank Nuretdin Kaymakçı, Bora Rojay and Erdin Bozkurt for their active support and constructive comments.