Seismic hazard assessment of the central North Anatolian Fault (Turkey) from GPS-derived strain rates and b-values

ABSTRACT The North Anatolian Fault (NAF) represents one of the most seismically active transform zones on Earth. It is characterized by high rates of crustal deformation that generate destructive earthquakes. These rates are induced by convergence of the northward-migrating Arabian and African plates with respect to the stable Eurasian plate. Therefore, the NAF represents a natural earthquake laboratory with a wide range of earthquake sizes (M ≤ 7.9) to investigate by using interdisciplinary approaches (seismological, magnetism, geological, gravitational, and geodetic studies). In this study, we compare the results of an analysis of b-values from seismicity and GPS (Global Positioning System) measurements of the strain rate to understand their coupling in terms of faulting and earthquake hazard implications. In particular, this comparison allows investigation of the spatial correlation between b-value and strain rate maps and is thus able to locate fault segments that have a high potential of generating large earthquake(s). b-Values range from 0.5 to 1.5 along the central NAF. The maximum principal strain rates are positive (tensile), and the minimum principal strain rates are negative (compressive). The surface strain is positive, showing that tensile strain is predominant in areas with high strain rates, consistent with the trend of the corresponding stresses.


Introduction
Turkey is one of the most seismically active regions in the world. Tectonic evidence suggests that Anatolia is moving westward relative to Eurasia (McKenzie 1972). The tectonic structure of Turkey is dominated by the Hellenic arc, North Anatolian Fault (NAF) and East Anatolian Fault (EAF). The NAF is a 1200-km-long strike-slip fault running from Karliova to the west (Şeng€ or et al. 2005). The NAF has a long history of seismicity and is characterized by dextral strike-slip motion. Although many studies have used seismic data and geodetic methods along the NAF, the central NAF is less studied using the combination of these datasets (b-values and strain).
Seismological and geodetic data provide key information for elucidating the kinematics of crustal motion. Seismological data can be a stress indicator, while GPS (Global Positioning System) observations can provide a high temporal sampling of deformation, depending on the time span of the data. Furthermore, a comparison of two datasets can provide insights into the seismic potential of the region of interest. The aim of this study is to investigate the central NAF seismic potential, which is inferred from GPS strain rates and b-values.
A b-value represents the relation between large and small earthquakes and is assumed to be related to the tectonic heterogeneity of the asperity region (Mogi 1962). The asperities create locked segments that resist faulting. In the case of the 1999 Izmit Mw = 7.4 mainshock, Westerhaus et al. (2002) had identified a low b-value volume only 10 km from the epicentre as the most likely location for a mainshock a year before it occurred. Gorgun et al. (2009) studied the area east of the Izmit mainshock rupture and found low b-values 25 days before the 1999 Duzce Mw = 7.2 mainshock. According to Gorgun et al. (2009), low b-values are related to a high Coulomb stress increase in the Duzce mainshock region. The reason for a low b-value is likely the state of stress near an asperity. Low b-values in some areas prior to a major seismic event can be followed by high b-values after the rupture occurs. The change in b-value indicates a coseismic stress drop if the b-value is used as a stress-meter, as suggested by Zang et al. (1998), Lei (2003), and . Ahmad et al. (2017) used an earthquake catalogue of Syria and GIS data to develop a seismic hazard map based on seismicity parameters (a-and b-values), concepts of local probability and recurrence time. Mapping asperities promises to be useful in identifying locations that may become crucial in future large earthquakes. Thus, high b-values indicate low stress, while low b-values indicate high stress within asperities.

Tectonic background
Although Central Anatolia has been regarded as a deformation-free region in many studies, Anatolia and the surrounding regions encompass a wide variety of tectonic phenomena, including transform strike-slip faulting (North and East Anatolia Faults), continental collision and major thrust faulting (Bitlis-Zagros, Caucasus), subduction (Nubia, Arabia), contraction (Caucasus, Marmara Sea), extension (Western Anatolia) and numerous small-scale processes (McKenzie 1972;Jackson and McKenzie 1984;Taymaz et al. 1991;Barka and Reilinger 1997;Kocyigit and Beyhan 1998). In addition, recent geodetic studies reveal that a gradual westward increase in velocities has been occurring (Reilinger et al. 2006) along the NAFZ. On the other hand, the results of the dense GPS velocity field presented by Aktug et al. (2013) indicate an ongoing internal deformation within Anatolia.
The fault motion of these large earthquakes suggests that they belong to the broad North Anatolian Fault Zone (NAFZ), where a transform plate boundary exists between the Anatolian and Eurasian plates in northern Anatolia. GPS vectors by Reilinger et al. (1997), McClusky et al. (2000, and Reilinger et al. (2006) indicate that northern Anatolia moves westward along the NAF. This motion of the Eurasian plate results into westward escape of the Anatolian plate along two prominently large transform faults: the NAF and EAF (Figure 1). These studies focused on the general tectonic frame of Turkey and lack sufficient spatial resolution to identify individual segments of the NAF without taking the elastic strain accumulation into account. However, recent regional studies provide the best estimates of the crustal thickness, deformation rates and locking depths, which are vital for estimating the strain accumulation rate and for further investigations of earthquake hazard models. Receiver function modelling suggests that the crustal thicknesses are thinner (25-30 km), relative to the high elevation of the Anatolian plateau (Zor et al. 2003;Gok et al. 2011). The recent seismic activity, the thickness of the crust and the deformation rates of the central NAF have been studied by Yavasoglu et al. (2011), Tatar et al. (2012, Kaneko et al. (2013), Ozener et al. (2013), Aktug et al. (2013), andPoyraz (2016). Recently, Yavasoglu et al. (2011) found that the locking depth in the central segment of the NAF near Ladik is 16 km, which is greater than the results of Reilinger et al. (2006) and Aktug et al. (2015), and provide an average locking depth along the NAF of 14.3 § 1.7 km, which is less than the estimates of previous studies (Peyret et al. 2012). Poyraz (2016) also suggests an average slip rate of 21 mm/yr and a locking depth of 12.72 km for the easternmost part of the NAF.
The region of interest in this study is bounded by the longitudes 31 -42 and latitudes 39 -41 . Figure 2 displays the seismic activity in the region. The activity is concentrated along the  easternmost part of the NAF and to the west. The largest earthquakes at the central NAF were the 1942 Erbaa (M 7.1) and the 1943 Tosya (M 7.3) earthquakes (Ambraseys 1970;Barka 1996).
Moreover, many geophysical and geological studies already have been conducted in the study area, other than the geodetic ones mentioned above. These studies determined that the crustal structure and active stress field are related to the variations in the strain rate (e.g. Yolsal-Çevikbilen et al. 2012;Warren et al. 2013;Karas€ ozen et al. 2013). In fact, interseismic deformation models cannot be differentiated by using solely geodetic data since the crustal structure and rheology that control the deformation at depth cannot be uniquely determined. Furthermore, combining the fault pattern with the present-day stress field is essential to gain further insight into the spatial distribution of earthquakes.

Data and method
GPS velocity field GPS measurements provide subcentimetre resolution of the displacement of the GPS stations relative to a stable reference frame. The continuously increasing GPS datasets allow quantification of more accurate slip rates, which can generate more precise strain maps. In this study, GPS data come from an analysis of a new GPS dataset and similar combination of available datasets covering the region of interest (Aktug et al. 2015). Using a dense GPS network that is widespread in the study area, Aktug et al. (2015) have analyzed the variation in kinematics and compared the results with geological studies. A GPS velocity field is obtained by combining different survey-mode and permanent GPS network solutions using the Kalman filter method in the GLOBK software (King et al. 2009). The reference frame was defined by estimating a 12-parameter transformation in ITRF2008. The combination strategy is similar to Nocquet (2012). The combined velocity field comprises 84 sites ( Figure 3).
Although these datasets have been analyzed by Aktug et al. (2015), here, we combine them into a single velocity field with GPS data and apply a strain analysis for a further comparison with b-values. Clearly, GPS measurements allow the dense monitoring of coseismic and postseismic strain transients, which can contribute to understanding the physics of the earthquake process; therefore, the capability of an interseismic GPS velocity field from medium-aperture geodetic networks can be evaluated.
It has been shown that slip rates along a long fault system can vary depending on the fault orientations and variations in elasticity, even if the driving velocity is constant (Ch ery 2008). This study shows that the westward increase in the GPS velocities within Central Anatolia due to a counterclockwise rotation confirms the results of Aktug et al. (2013). The southern part of the NAF has a pattern of an increasing rate of velocity from the east to west, which is indicated by an increase in the velocity of 15.3 mm/yr, to 20.4 mm/yr. The divergence in the directions of the velocity vectors west of the Sungurlu branch (SF) suggests a north-south extension caused by the southwestern motion of the Aegean plate. Horizontal slip along the NAF is not pure dextral strike-slip but rather transtensive (Dhont et al. 1998).

Frequency-magnitude relations (Gutenberg-Richter law)
The frequency-magnitude distribution (FMD) describes the relationship between the number of earthquakes events and their magnitude (Gutenberg and Richter 1944): where N(M) refers to the frequency of earthquakes with magnitudes greater or equal to M. In a semi-logarithmic plot, the constants a (zero-offset) and b (slope) can be determined. A high b-value indicates an increased material heterogeneity or crack density, whereas a low b-value indicates a high applied shear stress. The b-value has been interpreted to represent the presence of an asperity, the stress level and the material heterogeneities along the fault plane (Mogi 1962;Scholz 1968;Amitrano 2003;. The b-value is calculated using the maximum likelihood method (MLM) published by Utsu (1966Utsu ( , 1999.
b ¼ log 10 ðeÞ ðMÞ À ðM c À DM bin 2 Þ Â Ã ; Here, (M) is the mean of the binned magnitudes, M c is the magnitude of the completeness, and DM bin is the binning width of the catalogue.

Seismological data and magnitude of completeness (M c )
We use the catalogue from the Kandilli Observatory and Earthquake Research Institute (KOERI). We take into account all events from 1 January 2000 to the present. In general, we locate all the earthquakes using the 1-D (P-and S-waves) velocity model of Pinar et al. (2007). We choose events with at least four P-phases and one S-phase, hypocentre parameters with RMS errors <0.5 s, an azimuthal gap <180 and horizontal and vertical errors <1 km by using the HYPOCENTER earthquake location program (Lienert and Havskov 1995). The catalogue covers the period from 1 January 2000 to the present for shallow earthquakes with hypocentral depths <32 km. Figure 2 indicates the epicentral distribution along the NAF. The KOERI catalogue consists of 11,767 events (Figure 4(a)). The FMD method is fully capable of resolving the temporal changes in magnitude of completeness (M c ), as we demonstrate here for the case of the NAF seismicity ( Figure 2). We determine M c as a function of time for the entire period using a moving window approach. M c is calculated with a sliding window of 250 events, in steps of 50 events starting from 1 January 2000. Completeness immediately after the mainshock increases over that calculated during the background time by up to 2 magnitude units because small events are hidden in the coda of larger ones and because analysts are often overwhelmed by the considerable number of events. This well-known behaviour was documented, for example, by Wiemer and Wyss (2002) and Wiemer and Katsumata (1999). Detailed knowledge of M c is necessary to correctly assess seismic hazard (Reasenberg and Jones 1989;Wiemer ans Wyss 2002). Using the FMD as a guideline, we can derive with some confidence the temporal behaviour of M c as a function of time (Figure 4(b)). M c ranges from 1.8 to 3.3 for the entire catalogue.
To determine a complete catalogue, we compute the b-value during the seismicity as a function of depth ( Figure 5(a)) and time ( Figure 5(b)), by using an overlapping moving window approach and applying the entire magnitude range (EMR) method (Woessner and Wiemer 2005). On the basis of this analysis, we define b-values for the entire catalogue, ensuring that the catalogue is largely complete in both space and time (Woessner and Wiemer 2005;Woessner et al. 2006).
A bootstrap approach is used to estimate the errors in b and M c (Chernick 1999). Amor ese et al.

Spatial mapping of b-values
The variations in b-value are calculated using the public software package ZMAP by Wiemer (2001). Projecting all events onto the fault plane, we calculate the maximum likelihood b-values (Utsu 1966(Utsu , 1999 for each node of a 0.01 £ 0.01 grid from all events with M M c within a 5-km search radius, requiring a minimum of 50 events. Figure 5 presents the mean b-values calculated from 1000 mappings at each node. Repeating the mapping and using the mean of the different mappings substantially reduce the deviation from the original magnitudes. Although in many cases the deviations seem to concentrate on nodes along the margins of the section and in high b-value areas, we also find substantial changes at the centre of low b-value patches. We classify the amplitude of the observed variations by comparing it with the Amor ese et al. (2010) error. For 50 events (the minimum required number for our b-value calculations), the intrinsic error of the b-values is approximately 0.12. All deviations that are larger than this error are thus not expected from bootstrapping the original data at that point, and hence we regard these deviations as significant. According to this definition, the maximum observed errors indicate large areas with insignificant deviations. We represent high b-value areas along the NAF as depth sections in Figure 6. According to these three profiles, which show low b-values in the Bolu,Çorum, Tokat and Erzincan provinces, the BB 0 profile exhibits remarkably low b-values for the entire area. However, the AA 0 and CC 0 profiles also depict high b-values, especially at the western rims at distances between 0 and 20 km, as shown in Figure 7.

Results and discussion
There are many related examples of comparisons between geodetic data and seismicity from different parts of Earth (e.g. Ward 1998;Shen-Tu et al. 1998;Kreemer et al. 2000;Masson et al. 2005;Rontogianni 2010). These comparisons show that the geodetic strain rates are higher than the seismicity rates. In the literature, there are also many studies that have discussed the correlation between seismotectonic variables and GPS strain along active transform fault zones in Turkey (e.g. Oncel and Wilson 2004;Oncel and Wilson 2006;Aktug et al. 2013;Yavasoglu 2015). The results by Wilson (2004, 2006) reveal a positive correlation between average dilatation and the b-value and a negative correlation between the spatial fractal dimension D 2 and the GPS-derived geodetic strain in their region of interest. Regarding their results, the correlation between the b-value and GPS-derived dilatation suggests that the regions in compression have an increased probability of a rupture of greater magnitude. These studies involve different stages of the seismic cycle. However, this study investigates only the interseismic stage of the seismic cycle. We used a dense GPS velocity field along the central NAF with a large data span . We found a consistency in the results obtained with geodetic and geophysical methods, similar to Aktug et al. (2013), Aktug et al. (2015) and Yavasoglu (2015).
Strain rates are calculated from the geodetic data, which are based on position estimates using two or more GPS observations. The density and the distribution of the GPS sites might have an effect on the results of the strain rates. However, an overall comparison between the seismic parameters and geodetic findings provides insight into the seismic hazard in the target region. It is also important to explain the observed surface deformation. Biryol et al. (2011) investigated the strain rates and the dynamic interaction between lithosphere and asthenosphere (plate margin). Their results indicate a variation in the principal infinitesimal shortening and extensional strain axes from east to the west, following the counter-clockwise rotation of Anatolian plate. Biryol et al. (2011) suggest that the associated lithospheric and asthenospheric strain fields are the result of similar plate scale forces. Our findings agree with the results of Biryol et al. (2011), regarding the overall tectonic deformation of the Anatolian plate at this plate boundary. Figure 5 displays the spatial distribution of the b-values in map view. The areas with low b-values are located in the Bolu,Çorum, Amasya, Tokat and Erzincan provinces along the NAF. These low b-values are associated with high strain values. We plot low b-value regions as depth sections in Figure 6. These low b-value areas range from depths of 2-20 km. At these depths, we see a spatial correlation between several low b-value and high-strain-rate regions. However, this result does not reflect an overall correlation pattern for the entire target region.
We use a GPS velocity field to determine the crustal strain rate field as the gradient of displacement. SSPX software (Cardozo and Allmendinger 2009) is used to determine the strain accumulation in the area. Since crustal strain is assumed to occur in the horizontal plane, we use horizontal  velocities to infer the strain rates. The GPS sites are located at both sides of the NAF; therefore, the velocities used in this study are relative to a fixed Eurasian plate. The GPS sites in the study area are not distributed uniformly. To extract a regional tectonic pattern from the data, we used the grid-distance weighted method with a 20-km grid distance. The results are consistent with the long-term deformation of the area. The strain shows counter-clockwise rotation, which is consistent with the orientation of the regional geological structure. The results indicate that a significant strain accumulation occurs around Niksar (north of the intersection of the main branch of the NAF and the Sungurlu fault), where low b-values are prominent (Figure 8). The high strain values along borders of the study area can be artificial due to insufficient GPS data.
It can be seen from Figure 8 that the extension and shortening axes are approximately 45 to the dextral fault zone, as expected. This gridded homogenizing strain method is very useful when working with heterogeneous regional GPS networks, as in this study. On the other hand, the important features of the heterogeneous geology in such a relatively large region might be missed with this method. Many fault zones in this region splay from the NAF, running south-west across the Anatolian Plate for hundreds of kilometers (Bozkurt 2001). However, the results of the strain accumulation agree with the seismicity and geology in the region of interest.
According to the GPS velocity field, there is no deformation at Sungurlu fault (SF). However, strain has been accumulating where the SF and main NAF merge (Yavasoglu 2015). From the GPS results, such tectonic movements relatively oblique to the NAF imply that this fault does not strictly obey a model of an intracontinental transform fault, along which the upper crustal block motions must be parallel to its trace (Dhont et al. 1998). A number of basins occur within the NAF. One of them is the Niksar basin, which lies between the 1939 (M 7.9) and 1942 (M 7.1) earthquake fault segments of the NAF and does not present pure strike-slip movement (Tatar et al. 1995;Tatar 1996). Barka et al. (2000) suggest that the high-strain/low b-value area identified in this study has a distinct morphological expression of these basins. The results of Karasozen et al. (2013) also reveal that the Niksar-Erbaa pull-apart basin indicates localized extension due to transtension with notable spatial variations in the stress state. According to Yolsal-Çevikbilen et al. (2012), the lack of correlation and distributed seismicity along several segments of the NAF in this region might be an indicator of an interseismic stage, during which the strain energy is accumulating along locked portions of the fault.

Conclusions
Here, we use an improved velocity field for the region and two decades of geodetic data in an attempt to produce a strain rate map and an improved seismicity catalogue by using a new 1-D Figure 8. Strain rate field: principal strain axes and shear strain. velocity model (Karasozen et al. 2013). Although the seismic activity and the geological structure of an area are the essential inputs for seismic hazard maps, the strain rates produced by geodetic data are important for understanding the seismic hazard of that region. GPS data provide the deformation pattern as a strain rate map that is consistent with the longÀterm history of the geology in the region. The density and the distribution of the GPS stations provide an improved strain rate map. Clearly, a longer observation time span and additional GPS sites are needed for further research. Therefore, a new homogenous velocity field across all of Turkey for three decades is the work in progress, and we are currently studying the integrated strain map from this homogenous velocity field for the next step of this study.
The GPS and seismological data are used to generate various maps of regions that will potentially have large earthquakes in the future. These maps allow us to evaluate seismic hazards through spatial correlation between the GPS measurements and seismicity catalogues. The present study demonstrates the benefit of combining GPS and seismological data to improve seismic hazard assessment along the NAF and to thus differentiate between low-and high-potential regions for future earthquakes. Improved earthquake hazard maps will allow fine-tuning microzonation of these areas for city infrastructure planning.