Source Mechanisms and stress field of the 2017 Ayvacık/Çanakkale earthquake sequence in NW Turkey

Ayvacık district of Çanakkale province was hit by four moderate size earthquakes on 6 February 2017 03:51 M w 5.2, 6 February 2017 10:58 M w 5.1, 7 February 2017 02:24 M w 5.2, 12 February 2017 13:48 M w 5.3. Centroid moment tensor solutions for 60 events with moment magnitudes (M w ) between 3.4 and 5.3 are computed by applying a waveform inversion method on data from the Turkish and Greek broadband seismic networks. The time span of data covers the period between 1 January and 28 February 2017. The mainshock is a shallow focus normal event with strike-slip component at a depth of 10 km. Focal depths of aftershocks range from 3 to 24 km. The seismic moment (M o ) of the mainshock is estimated as 7.62 × 10 16 Nm. The calculated rupture duration of the Ayvacık mainshock is 3 s. The focal mechanisms of the aftershocks are mainly normal faulting with a minor strike-slip component. The geometry of focal mechanisms reveals a normal faulting regime with NE-SW trending direction of T -axis in the entire activated region. We perform stress tensor inversion to acquire more accurate picture for the regional stress field. The stress tensor inversion results indicate a predominant normal stress regime with a WNW-ESE oriented maximum horizontal compressive stress (S H ). With respect the newly determined focal mechanisms, the effect of the propagation of the North Anatolian Fault into Aegean Sea is very clearly pronounced. According to the double-difference relocation algorithm for the 2017 Ayvacık seismic sequence, ten clusters are revealed. In order to find these clusters, we use HypoDD software. According to the relative relocations from HypoDD, we define 10 clusters along the five profiles. The aftershock activity in the observation period between 1 January and 28 February 2017 extends from W to E direction. Seismic cross-sections indicate that a complex pattern of the hypocenter distribution with the ruptured segments. The western cluster is associated with a fault plane trending mainly N-S and dipping vertical, while the eastern is related to a fault plane trending NE-SW and dipping towards NW. The best constrained focal depths operate in the approximate depth range from 3 to 24 km. Coulomb stress analysis is performed to calculate the cumulative effect of all the earthquakes of the swarm and correlate it with the activated region. Positive lobes with stress more than 3 bars are obtained, these values are large enough to trigger failure towards NW-SE of the swarm activity.


Introduction
The Ayvacık earthquake (EQ) occurred at 03:51:40.0 GMT on 6 February 2017 in Çanakkale province. The mainshock was a moderate size (M w =5.2) event at a depth of 10 km. The Ayvacık EQ is directly located in the westernmost part of Biga Peninsula (BP) along the western North Anatolian Fault (NAF) strike-slip fault system ( Figure 1a). Biga Peninsula is one of the most seismically active parts of the northwestern Anatolia. It is dominated by a series of strike-slip structures bounded by normal or oblique faults [McKenzie, 1978;Şengör et al., 1985;Taymaz et al., 1991;Le Pichon et al., 1995;Barka et al., 1997;Kreemer et al., 2004;Gürer et al., 2006;Kürçer et al., 2008;Ateş et al., 2009;Le Pichon and Kreemer, 2010;Fichtner et al., 2013;Selim and Tüysüz, 2013;Vanacore et al., 2013;Kind et al., 2015]. The main seismotectonic structures of the BP region are the Edremit Fault Zone (EFZ), Çan-Biga Fault (ÇBF), Evciler Fault (EF), Gülpınar Fault (GF) and Kestanbol Fault (KF) that compose transform plate boundaries where the Eurasian plate (EU) to the north is sliding relative to the Anatolian (AT) and Aegean Sea plates (AS) to the south (inset in Figure 1a and b). The NAFZ is one of the most active transform faults worldwide extending along 1600 km from eastern Anatolia to the Aegean Sea ( Figure 1a). The direction of slip corresponds well with the Global Positioning System (GPS) derived 25-30 mm yr −1 westward motion of the Anatolian Block with respect to Eurasia [McClusky et al., 2000;Reilinger et al., 2006]. The deformation zone associated with the transform is very wide from the Thrace to Aegean Graben system [Taymaz et al., 1991;Kreemer et al., 2004;Nyst and Thatcher, 2004;Gürer et al., 2006;Selim and Tüysüz, 2013]. Approximately 45° NE trend is observed as a strike-slip deformation along the BP. Strike-slip motion is transformed into normal on southwest part of the BP (see Figure 1b) [Kiratzi, 2002;Kreemer et al., 2004;Gürer et al., 2006;Aydoğan, 2007]. BP is located SW of Marmara region describing the southern branch of NAFZ ( Figure 1). Figure 1b indicates an overview on the seismotectonic settings of BP region with Aegean Sea.
BP forms the southwestern part of the Marmara Sea. In this region, it is correlated a wide strike-slip fault zone with an oblique motion [Erdik et al., 2004;Bohnhoff et al., 2016]. The focal mechanisms of recent moderate-size Ethem Görgün et al.  Anatolia [McClusky et al., 2000;Reilinger et al., 2006]. Seismically active faults are shown by red lines [Şaroğlu et al., 1992]. Blue, green and yellow triangles with station codes depict locations of the KOERI, AFAD and AUTH broadband seismic stations, respectively. The epicenter of the 2017 Ayvacık M w 5.2 EQ at 03:51 is indicated by the red star. The inset in the above left shows whole Turkey. Significant active faults are indicated by the red lines in Biga Peninsula [Emre et al., 2013]. Two closely spaced red single arrows display shear sense of major faults (NAT: North Aegean Trough, NAFZ: North Anatolian Fault Zone, EAFZ: East Anatolian Fault Zone). Boundaries (heavy colored red lines) of the Aegean Sea (AS) and Anatolian (AT) plates, which are surrounded by the African (AF), Arabian (AR) and Eurasian (EU) plates (Bird, 2003). The solid black rectangle shows the study area, which is enlarged in Fig. 1a. b) For reference, focal mechanisms of the previous significant earthquakes are plotted [Kalafat et al., 2009]. Black and white circles on beach-balls exhibit P and T-axes, respectively. Red star represents the 6 February 2017 M w 5.2 03:51 mainshock. earthquakes that have occurred since 1985 around the BP is shown in Figure 1b [Kalafat et al., 2009]. The distribution of these epicenters indicates a low activity compared with Marmara Sea [Erdik et al., 2004;Bohnhoff et al., 2016].
These earthquakes are mostly upper plate events (depth ≤ 25 km) [Kalafat et al., 2009]. The seismic pattern of BP and surrounding area is generally characterized by moderate and large size events [Kalafat et al., 2009]. The focal mechanisms of some moderate recent events are related to the NAFZ strike-slip and oblique system combination.
We interpret their orientation and spatial distribution in terms of the major strike-slip fault with normal motions. Ganas et al. [2018] stated that Ayvacık EQ has normal slip kinematics, dimensions of 6 by 6 km (length, width) and dips at 45° according to geodetic inversion results. Their INSAR analysis also showed that there is a uniform slip of 28 cm. Özden et al. [2018] also pointed out that analysis of the recent Ayvacık earthquake sequence indicates the domination of the normal faulting stress regime in the southern part of Biga Peninsula with a consistent NNEtrending 3 axis. Mesimeri et al. [2018] indicated that spatial distribution of Ayvacık EQ sequences revealed a south-dipping causative fault along with secondary and smaller antithetic segments. In this respect, the seismicity started at the western part and migrated to the eastern part to the activated area. Svigkas et al. [2019] also stated that geodetic observations constrained the rupture distribution of the main events and they found a single fault with striking N110°E and dipping ~40° to SW.
In this study, a stress tensor inversion of earthquake focal mechanism data is performed to obtain a more accurate picture of the regional stress field. For this purpose, seismic waveforms at local and regional distances are used to calculate source parameters of the 60 events (5.3 ≥ M w ≥ 3.4) of the 2017 Ayvacık EQ sequence using the waveform inversion method [Nakano et al., 2008]. This provides additional information on the stress field that may improve kinematic models for the Ayvacık EQ area and thus develop the understanding of the Biga Peninsula local and regional tectonics. Furthermore, the double-difference relocation algorithm [Waldhauser and Ellsworth, 2000] and the Coulomb stress analysis [Toda et al., 2005] are applied to determine the expanded spatial distribution of the Ayvacık EQ sequence. Determination of accurate stress field orientations and event locations using data from local and regional seismic networks are crucial for investigations of the seismotectonics in and around southwestern Marmara Sea region. Thessaloniki (AUTH) broadband seismic stations used in this study. We use seismic records obtained from 24 KOERI stations (blue triangles), 2 AFAD stations (green triangles) and 4 AUTH (yellow triangles) stations in this study (see Figure 1a).

Data and waveform inversion method
Centroid moment tensor (CMT) solutions of earthquakes in the Ayvacık region are computed using the waveform inversion method developed by Nakano et al. [2008]. In this approach, if the location yielding a minimum residual lies at the edge of a search area, the grid is extended to surround the location of the minimum residual. When the location of minimum residual lies within a search area, a new search is performed around the location using a reduced grid spacing to find a detailed source location. The seismic moment and rupture duration are estimated from the deconvolved form of the moment function [Nakano et al., 2008;.
Three-component seismograms are used for the inversion of the Ayvacık EQ sequence. Seismograms with good data quality (high signal-to-noise ratio) are selected. The average number of waveforms used for the CMT analysis is 10. The observed velocity seismograms are corrected for instrument response and then integrated in time to obtain the displacement seismograms. Waveforms are bandpass filtered between 20 and 50 s and decimated to a sampling frequency of 0.5 Hz. A total data length of 512 s (256 data points in each channel) is used for the inversion.
Green's functions are generated using the discrete wavenumber method [Bouchon, 1979], assuming the crustal structure model of Akyol et al. [2006] for the calculation. The synthetics are calculated for a horizontally layered structure given in Table 1. Green's functions are computed for every 10 km of epicentral distance up to 1500 km. At each of these radius steps, three-component displacements are calculated at every 1° in azimuth for each basis of moment tensor. Green's functions are also calculated at every 5 km for source depths (hypocentral depth) shallower than 100 km. For spatial grid search, hypocenter locations estimated by KOERI are used as an initial location.
Adaptive grid spacings, in which the grid spacings are gradually decreased in each step of the search, are also applied.

4
Spatial grid search is started with a horizontal grid spacing of 0.5° and a vertical grid spacing of 10 km. In the next step, the grid spacing is reduced to 0.2° horizontally and 5 km vertically. Finally, the horizontal grid spacing is reduced to 0.1° and vertical is 5 km. At each grid point, the fault and slip orientation parameters (strike, dip and rake angles) are searched in 5° steps. For each combination of source location, fault and slip orientation parameters, the waveform inversion is carried out to estimate the best-fitting source parameters (see for details Nakano et al., [2008] and Nakano et al., [2010]).
The CMT solutions are obtained using data not only from stations in Turkey but also from Greece (yellow triangles in Figure 1a), and therefore azimuthal coverage is sufficient. The inversion method of Nakano et al. [2008] uses the double-couple constraint, which stabilizes the inversion solution and reduces the trade-off between source location and non-double-couple components (see for detail Nakano et al. [2008]). Furthermore, all stations used are located in similar crustal structures as Thrace and Greece thus; possible effects of the structural contrast may be minimized. Determination of accurate source parameters, especially source locations, using data from the local and regional seismic network are crucial for investigations of the seismotectonics in and around Turkey and Greece.

Stress tensor inversion
The stress tensor has six unknowns, either three principal stresses and orientations, or three normal and three shear stress components (e.g., Zang and Stephansson [2010]). Four of the unknowns are resolved by the inversion of the stress tensor, the fifth unknown is calculated by the assumption that slip occurs in the direction of maximum shear stress [Wallace, 1951;Bott, 1959] and the sixth unknown is usually resolved using the assumption that the stress tensor is homogeneous and constant in the binning region throughout the time interval of interest.
In this study, the technique of Michael [1984;1987] is applied to the selected 60 events (residuals are smaller than 0.4 according to Nakano et al., [2006] criteria; Table 2). The Michael-approach provides a more appropriate estimate of uncertainty, compare to Gephart and Forsyth [1984] approach [Hardebeck and Hauksson, 2001] The algorithm uses the statistical method of bootstrap re-sampling and allows determining the orientation of the three principal stresses ( 1 = maximum principal compressive stress, 2 = intermediate and 3 = minimum) as well as the stress ratio R = ( 2 -3 )/ ( 1 -3 ), also called relative stress magnitude [Bott, 1959]. The R is defined using the standard geologic/geophysical notation with compressive stress positive and 1 > 2 > 3 [Zoback, 1992]. The stress ratio (R) ranges from 0 to 1. Values of R < 0.5 and R > 0.5 indicate a transpressional and transtensional regime, respectively. All parameters are determined by finding the best fitting stress tensor to the observed focal mechanisms.
Assumptions that must be fulfilled by the input data are: (1) stress is uniform in the area of interest during the observed time interval, however this assumption cannot be entirely valid given the diversity and complexity of the structures, (2) earthquakes are shear-dislocations on pre-existing faults, (3) similar shear stress magnitude are present on each fault and (4) slip occurs in the direction of the resolved shear stress on the fault plane. When analyzing real data sets, the uncertainties of the retrieved stress directions, the shape ratio, the fault selection and overall friction on faults are verified by the bootstrap method as proposed by Michael [1987] and Vavryčuk [2014].
Here the uncertainties are calculated as the maximum differences between the results of the inversion for noisefree and noisy data with 1000 noise realizations. This method is convenient for several reasons. First, it can take into  account the case when some nodal planes are more uncertain than the others by specifying differently noise levels for the fault orientations and slip directions. And secondly, the obtained uncertainties are more realistic. The efficiency of the iterative inversion was tested on the Ayvacık data. The method worked well and produced a significantly more accurate shape ratio than the inversion of the focal mechanisms with randomly selected faults.
To quantify the misfit between the best stress tensor and the data, the angle between the calculated slip vector from stress tensor inversion and observed slip vector from fault plane solutions is used. This angle is referred to as . The angle refers to the mean value of for the data in a single inversion [Michael, 1987]. A synthetic control study showed that the amount of heterogeneity in the stress field could be characterized by the average misfit between the observed and predicted slip directions ( ) [Michael et al., 1990]. If ≤ 33°, stress tensors are spatially uniform. If > 33°, the inversion result is interpreted in terms of a spatially heterogeneous state of stress (Michael, 1991). This 33° marks the transition from an uniform stress field to a heterogeneous stress field. In our study, stress tensor inversion results give angles smaller than 33°. Heterogeneity of the stress field was documented in the average misfit level of the inversion. For each stress inversion, 2000 bootstrap iterations are performed.
Michael's algorithm is quite fast and accurate when compensating directions of the principal stresses. It gives a reasonable accuracy even for randomly selected fault planes in focal mechanisms. For example, Vavryčuk [2014] indicated that the principal stress directions are determined with a similar accuracy regardless of whether the faults are randomly selected from nodal planes or correctly identified by the inversion. The behavior of the shape ratio is, however, remarkably different. Data with randomly selected faults yield a broad distribution of the shape ratio with its maximum significantly biased as compared to the inversion with correctly selected faults. The maximum error is 44 and 12 percent for randomly and correctly selected faults, respectively. He calculated 1000 realizations of noise which produced an accuracy similar to that of the retrieved focal mechanisms (mean error of 6° in the fault normal and slip direction). However, the accuracy of the stress ratio is significantly lowered when the fault planes are not correctly chosen. The stress ratio is more sensitive to the correct choice of the fault plane than the principal stress directions and substituting the faults by the auxiliary nodal planes introduces high errors [Lund and Slunga, 1999;Vavryčuk, 2014]. This difficulty is removed by modifying Michael's algorithm and inverting jointly for stress and for fault orientations [Vavryčuk, 2014]. The fault orientations are determined by applying the fault instability constraint and the stress is calculated in iterations. As a by-product, overall friction on faults is determined. Vavryčuk [2014] showed that the new iterative stress inversion is fast and accurate, and performs much better than the standard linear inversion. Since the iterative stress inversion is based on Michael's method, it can easily be implemented in these codes enhancing their accuracy. Vavryčuk [2014] proposed the MATLAB code of this inversion called STRESSINVERSE that is provided on the web page (http://www.ig.cas.cz/stress-inverse).

Relocation of hypocenters
We locate all earthquakes using the 1-D (P-and S-wave) velocity model of Akyol et al. [2006]. We choose only events with at least 5 stations with 5 P-and 2 S-phases, hypocenter parameters with RMS errors < 0.5 s, azimuthal gap < 180° and horizontal and vertical errors < 1 km, using the HYPOCENTER earthquake location program [Lienert and Havskov, 1995]. 3140 events which are recorded by at least five stations for two months period following the date of 1 January 2017. We calculate the absolute hypocentral parameters of aftershocks with the hypocenter location algorithm described by Lienert and Havskov [1995]. The catalogue covers the period from 1 January to 28 February 2017 for shallow earthquakes with hypocentral depth < 25 km. The catalogue consists of 3140 events. The average horizontal and vertical uncertainties of the 3140 events (M l ≥ 2.0) are found to be 3 and 4 km, respectively. In minimizing the errors in the location parameters, network geometry, phase reading quality and uncertainties in the crustal structure are restricting factors. Relative earthquake location methods can improve absolute hypocenter locations. For this purpose, the double-difference relocation algorithm HypoDD, Waldhauser and Ellsworth, [2000] is used. The algorithm assumes that the difference in travel times for two close events observed at one station can be connected to the spatial offset between the events with high accuracy. Absolute location parameters of aftershocks and P/S travel time differences between the event pairs are used for the inputs. Relative location errors of the 2326 relocated events are smaller than 1 km in horizontal and vertical directions.
We divided the catalogue into two parts. The first part of the catalogue ranges from 1 January to 6 February. This part exhibits the seismicity prior to Ayvacık EQ sequence. The second part of the catalogue includes aftershocks of the 2017 Ayvacık EQ sequence. Figure 2 and Figure 3 display the seismicity before and after the 6 February 2017 Ayvacık EQ sequence, respectively. From the spatial distribution of the epicenters (Figure 2), it is observed that the epicenters indicate a prominent NW-SE alignment. The along strike dimension of the activated zone is ~30 km and its width is greater than 10 km, prior to the February 2017 Ayvacık mainshock (M w =5.2). This observation is also represented in the depth cross-section along N-S direction (Figure 2, profile A-A', indicated by blue). The event depths range from 3 to 24 km.
A number of aftershocks are located off-shore (A-A' and B-B' profiles) with respect to the main aftershock cluster.
The swarm activity in the observation period between 6 and 28 February extends from the W to the E direction Ethem Görgün et al.
6 indicate that the aftershock sequence is mainly confined in the crust (depth < 25 km) and is running in the approximate depth range from 3 to 24 km. In order to investigate evaluation of aftershocks along the Ayvacık EQ area, we examine spatial distribution of aftershocks along 5 profiles. Aftershocks are clearly activated in the ~50 km length and ~10 km width of fault rupture zone (Figure 3). Along these 5 profiles, we detect ten secondary fault segments in the study area based on double-difference algorithm. Ten secondary fault segments are related to visual inspection after application of double difference algorithm. We indicate main active faults (red lines in Figure 3) along the Ayvacık earthquake region and Figure 3 also shows a possible Ayvacık earthquake rupture zone (white dashed line) according to aftershock distribution. We display 3D view of aftershock distribution for Ayvacık EQ sequence in Figure 3 in accordance with to the 5 profiles.
In the westernmost end of the main cluster, the seismicity occurred in a zone of ~10 km wide and the depths of aftershock are found to be between 5 and 20 km (Figure 3,  Furthermore, Figure 4a and b indicates cumulative number of EQ as a function of time before and after the February 6 mainshock, respectively. Figure 4 also present clearly temporal variations of the events occurred before and after the February 6 mainshock for the entire EQ catalogue.

7
The 2017 Ayvacik Earthquake Sequence First, the source centroid location and focal mechanism of these events are estimated. Figures 5 and 6 show the best fitting source location of 60 events obtained from the waveform inversion. Using the hypocenter location determined by KOERI as an initial source location, the waveform inversion searched for the best-fitting source parameters. Revised source locations vary in average 4 km horizontally and 5 km vertically with respect to the initial locations of KOERI, but this variation does not change the general pattern in seismicity. Location differences depend on velocity model, station coverage and number of stations use to calculate the hypocenter. Data from some of near stations (distance < 50 km) are not used because the waveforms are clipped or saturated. The focal mechanism of the 6 February 2017 M w 5.2 03:51 EQ (event #4, Table 2) obtained at the best fitting source location shows a normal motion with significant strike-slip component (Figure 7), characterized by a nodal plane with strike, dip and rake 333°, 76° and -63°, respectively. Waveforms recorded by four stations (BZC, green triangle; EZN and GADA, blue triangle and SIGR, yellow triangle; Figure 1a), which are located at epicentral distances ranging from approximately 30 to 80 km are used. The seismic moment of the mainshock is estimated as M o = 7.62 × 10 16 Nm, and the corresponding moment magnitude is M w = 5.2. The seismic moment function obtained from the waveform inversion is shown in Figure 7 (upper right). The rupture duration estimated from the moment function is 3 s (Figure 7, upper right side). Waveform fits between observed and synthetic seismograms calculated for the best fitting source parameters are also displayed in Figure 7 (Nakano et al., [2008] for details on the approach).
Waveform fits are in good agreement with a normalized residual (R) of 0.05 (Figure 7, lower left).
The 6 February 2017 10:58 M w 5.1 event is investigated (event #8, Table 2, Figure 8). Waveforms recorded by four stations which are located at epicentral distances ranging from 35 to 75 km used for the previous event. The seismic moment of the 10:58 event is estimated as M o = 5.70 × 10 16 Nm, and the corresponding moment magnitude is M w = 5.1. This event is located at the SW of the 03:51 event. The estimated focal mechanism is indicated in Figure   8; a nodal plane corresponds to 300° strike, 60° dip and -60° rake angles. The focal mechanism of this event is normal faulting. This indicates that this event does not have the similar focal mechanism as the 03:51 event.

Ethem Görgün et al.
8  The compressional axis is NW-SE. The fits between observed and synthesized seismograms are also displayed in The 7 February 2017 02:24 M w 5.2 event is also examined using waveform inversion technique (event #8, Table 2, Figure 9). The seismic moment of the 02:24 event is estimated as M o = 6.70 × 10 16 Nm. This event is located at the SE of the 03:51 and 10:58 events. The focal mechanism of this event is computed as a pure normal fault with a compressional axis is E-W direction. The small residual (R = 0.05) exhibits that the observed seismograms are in good agreement with the synthetic ones (Figure 9, red and black seismograms).
The 12 February 2017 13:48 M w 5.3 event is the largest EQ in the catalogue. We examined using waveform inversion technique same as former events (event #43, Table 2, Figure 10). The depth of this event is shallow and comply with the former M w > 5.0 events. The seismic moment of the 13:48 event is estimated as M o = 1. 16 × 10 17 Nm. This EQ is located at the same coordinates with 02:24 event. We estimate its focal mechanism as a normal fault with a compressional axis is NW-SE direction. The waveform fits between observed and synthesized seismograms are also shown in Figure 10 (red and black seismograms). We calculate a small residual (R = 0.05) for this event that observed seismograms correspond with the synthetic ones ( Figure 10).

Analysis of stress tensors
The results of the stress tensor inversion are presented in Figure 11. The measure of the misfit values is evaluated by the areas of 95% confidence limit (see Figure 11). Stress tensor inversion of the Ayvacık EQ aftershocks reveals a substantial high level of uniform stress for the entire catalog ( = 23°) according to Michael (1991). The maximum principal stress, 1 , trends N30°W with a plunge of 4° and the minimum principal stress, 3 , trends N182°E with a plunge of 84° (Figure 11). The resulting stress tensor corresponds to a transtensional stress regime (R = 0.75). The orientation of the maximum horizontal compressional stress (S H ) is WNW-ESE and normal faulting stress regime is found for the study region.
Low variance (0.10) indicates a good fit of homogeneous stress tensor to the observed focal mechanism, and shows low heterogeneity of the stress field [Michael et al., 1990;Lu et al., 1997;Wiemer et al., 2002;Görgün et al., 2010]. The alignment of 1 is NW-SE. Deformation regimes obtained from the stress tensor inversion indicate the predominant normal faulting regime with strike-slip component in the mainshock area ( Figure 10). Wiemer et al. [2002] concluded that a low degree of heterogeneity in stress field (variance < 0.2) indicates uniformity. There is no significant variance in stress tensor orientations (variance < 0.2). Thus, to first order, the Ayvacık EQ region is characterized by a homogeneous intraplate stress field.

Calculation of the Coulomb stress change
The static Coulomb stress change caused by a mainshock has been widely applied to assess areas of subsequent off-fault aftershocks [Reasenberg and Simpson, 1992;Toda et al., 2008;2011a]. The Coulomb stress change is defined as CFF = + μ , where is the shear stress on the fault (positive in the inferred direction of slip), is the normal stress (positive for fault-normal tension) and μ is the apparent friction coefficient. Failure is promoted if CFF is positive and inhibited if negative; both increased shear and unclamping of faults are taken to promote failure, with the influence of unclamping controlled by fault friction [Toda et al., 2011a]. After the occurrence of an earthquake, the areas with positive values of CFF are loaded with stress. To resolve the Coulomb stress change on a 'receiver fault' requires a source model of the earthquake fault slip, as well as the geometry and slip direction on the receiver [Toda et al., 2011a]. We can find the receiver faults at every point in the Coulomb stress increase that maximize the earthquake stress change and the tectonic stress [King et al., 1994], termed the 'optimally-oriented' Coulomb stress change [Toda et al., 2011a].
The source parameters deduced from waveform modelling are used to calculate the Coulomb stress change for the 4 largest events ( Table 2). The stress tensor principal axes ( 1 , 2 and 3 ) are also added to compute the spatial Ethem Görgün et al.
16 Figure 11. a) P/T-axes for the events in Fig. 4 (Table 2) with retrieved principal stress directions (green signs on the figure).
The P and T axes are marked by red circles and blue plus signs, respectively b) Confidence limits of the principal stress directions retrieved by the iterative method [Vavryčuk, 2014]. Red, green and blue colours correspond to the 1 , 2 and 3 stress directions, respectively. The contours around the best model define the 95% confidence region of each stress axis. c) Histograms of the stress ratio (R) which is found to be 0.75. This value indicates that the study region is described by transtensional stress regime. distribution of the Coulomb stress change (Figure 12). The stress changes are resolved at 10 km depth based on the estimated geometry of major active faults. Increased shear stress in the rake direction and unclamping on surrounding "receiver" faults are interpreted to promote failure [Toda et al., 2008]. We resolve the shear and normal (clamping/unclamping) components of the stress change on a 2D grid of points or on specified "receiver" fault planes. We make a regular grid with receiver faults every 5 km. Receiver faults are planes with a specified strike, dip, and rake, on which the stresses imparted by the source faults are resolved; we call these "specified faults" in the Coulomb stress calculation. The shear stress increase or decrease is dependent on the position, geometry, and slip of the source fault and on the position and geometry of the receiver fault, including its rake. The normal stress change (clamping or unclamping) is independent of the receiver fault rake [Toda et al., 2011b]. The Coulomb 3.3 software (www.coulombstress.org) is used to project normal faulting with significant strike-slip component mechanism at a depth of 10 km with a strike/dip/rake of 333°/ 76°/ -63° and a friction coefficient μ = 0.4 in an elastic half space with uniform isotropic elastic properties [Lin and Stein, 2004;Toda et al., 2005;2011a].
The spatial distributions of the Coulomb stress change due to the occurrence of the 2017 Ayvacık EQ sequence M w > 5.0 events are presented in Figure 12. We analyze four mainshocks (the same as we described in the preceding section) with M w bigger than 5.0. According to our analysis, four mainshocks indicate similar Coulomb stress failure patterns about Ayvacık EQ sequence. Stress increase regions brought closer to failure are represented by the red lobes. the Ayvacık EQ sequence ruptures are brought to 3-5 bars closer to failure. Stress changes of ≥ 1 bars are generally observed to influence seismicity rates [King et al., 1994]. Prior to February 6 M w 5.2 mainshock, there is another moderate event occurred on 14 January 2017 with a M w 4.5. This event has identical focal mechanism with four Mw > 5.0 earthquakes. 14 th January event is the foreshock just before occurring February 6 moderate size earthquakes. This event is important for Coulomb stress distribution since M w > 5.0 earthquakes are quite close to each other.
To further compare the Coulomb stress change and Ayvacık EQ aftershocks, a cross-section is presented in Figure   13 along the route A-B (perpendicular to the fault geometry) shown in Figure 12. Most of Ayvacık EQ aftershocks (M l ≥ 3.0) lie in the region of < -1 bar Coulomb stress decrease (Figure 12 and 13). Nevertheless, some of aftershocks (especially 6 February 03:51 and 7 February 02:24 events, Figure 13, A-B profile) M l ≥ 3.0 lie above the fault plane, where the Coulomb stress is calculated to have increased.

Discussion and conclusions
The spatio-temporal distribution and source characteristics of the February 2017 Ayvacık EQ sequence are investigated. The distribution of relocated hypocenters and focal mechanisms clearly indicate the activation of a N182°E trending normal faulting system. The stress tensor inversion results reveal a predominant normal stress regime with an almost horizontal (E-W, for entire 60 events) and NW-SE (according to four M w >5.0 events) orientation the maximum horizontal compressional stress (S H ). The entire study region is rather homogeneous according to the stress tensor inversion results (variance < 0.2, Figure 11). The relative stress magnitude (R) is estimated to be 0.75 (transtensional stress regime) for the Ayvacık EQ region. Therefore, to first order, the Ayvacık EQ region is characterized by a homogeneous intraplate stress field as far as the orientation of stress is concerned (Figure 11).
High-resolution aftershock locations reveal a complex pattern of the hypocenter distribution with the activation of ten different fault segments (Figure 3, defined by red and black dashed-lines). The westernmost segment (blue, see Figure 3) is related to a fault plane trending NE-SW and dipping towards SSW, while the eastern segments are Ethem Görgün et al.
18 Figure 13. Cross-sections of Coulomb stress distribution for four mainshocks. Profiles are shown in Fig. 10. Red line donates coseismic ruptures of the mainshocks according to fault geometry based on this study (Table 2). White circles represent earthquake hypocenters. The size of the solid white circles is proportional to magnitude. related to a fault plane trending E-W and dipping towards almost horizontal and SE. There are several smaller-scale transtensional structures along faults in off-shore Biga Peninsula where close to the 2017 Ayvacık EQ sequence occurred. According to the seismic reflection data from Boztepe-Güney et al. [2001] and Yaltırak et al. [2012], offshore Biga Peninsula is active, because it exhibits sea-bottom scarps and cuts across the western prolongation of the NAFZ. The dimensions of the activated region do not remain constant during the evolution of the seismic sequence, even though the seismicity gradually decreased. To explain the observed extended spatial distribution, the double-difference relocation algorithm is performed, using the absolute locations as an input calculated in the framework of this study. The analysis reveals prominent eastern and western segments of the activated seismogenic region. The westernmost cluster (blue, Figure 3) has relatively small magnitude events with M l ≤ 3.5. The eastern segments (orange and purple, see Figure 3), where the mainshock occurred, indicates the seismicity along the east direction and to cause the observed elongation of the activated region.  (Figure 1b). In a broad context, the Aegean Sea region (including west of Anatolia) is a domain characterized by diffuse deformation (widespread plate boundary deformation between Eurasia and Anatolian plates), which is affected by both the Aegean Sea back-arc extension, approximately N-S, and by the shear transferred from east or north, due to the movement of the North Anatolian Fault Zone [Nyst and Thatcher, 2004;Le Pichon and Kreemer, 2010;Reilinger et al., 2010]. BP faults located in well-defined extensional basins are controlled by NW-SE trending normal and NE-SW trending strike-slip movements. These five events are example of a NE-SW trending strike-slip faulting with normal component (Figure 1b). This type of strike-slip faulting is predominantly observed in the northern Aegean Sea region and further to the south (Figure 1). The sequence investigated here is a typical example of pure NW-SE trending normal motions.

19
The 2017 Ayvacik Earthquake Sequence Figure 14. Topographic map of Ayvacık EQ sequences after the 6 February 2017 M w 5.2 03:51 mainshock. Blue circles indicate aftershock epicenters following the mainshock. Red lines represent the active fault segments in vicinity of the EQ swarm region [Emre et al., 2013;Özden et al., 2018]. White dashed line shows the 2017 Ayvacık EQ rupture. Red arrows depict maximum (S H ) and minimum (S h ) horizontal compressive stresses.
During the 2017 seismic sequence in the Ayvacık area, an epicenter migration is also observed for the 4 larger magnitude earthquakes. The spatio-temporal distribution of Ayvacık seismic activity exhibits that the master event of the 6 February 03:51 is followed by three other equivalent earthquakes within the rupture zone ( Figure 14). This migration appears to follow a direction of progressively increasing stress, from the north-western part of the EQ cluster, towards the south-eastern part (Figure 3). In agreement with these results is also the Coulomb stress redistribution due to the occurrence of the first two events that indicate a maximum stress transfer from the northwest, towards the epicenters of the following two large shocks that occurred subsequently in the southeastern part of the rupture zone.
In this study, all aftershocks are shallower than 25 km. Furthermore, similar findings were published by Akyol et al. [2006]. Accordingly, it can be speculated that the thickness of brittle seismogenic crust in Biga Peninsula area is about 25 km. This aftershock study also indicates that the majority of aftershocks occur in the crust shallower than 20 km with the exception of few events about 23 km (in the main cluster, profiles B-B' and C-C'). Fault and then the aftershock activity spreads out NW-SE direction. This rupture zone is indicated by the white dashed line. This can be interpreted as the rupture zone is parallel to the Tuzla Fault. This zone is also located between the Tuzla and Kocaköy Faults. We state that this rupture zone that belongs to February 2017 EQ swarm activity, is an independent and unknown segment from the existing fault lines. Kandilli Observatory and Earthquake Research Institute (KOERI) intensity data distribution is correlated with this new fault segment in the Ayvacık EQ region. According to KOERI intensity map, event locations along this fault zone are related to high intensity values (www.koeri.boun.edu.tr). This new segment has not only on-shore but also off-shore parts. Almost half of the EQ epicenters are located along the off-shore sections of this rupture zone ( Figure 14). Livaoğlu et al. [2018] also indicated a damage distribution map related to the Ayvacık EQ sequences. According to their findings, damage ratio of structures is higher in villages close to epicenters along the new segment that found in this study.
To conclude, the Ayvacık EQ sequence that occurred in the north of Aegean Sea off-shore Biga Peninsula is clearly indicating that active NE-SW trending ( 3 ) normal faulting systems are wide-spread in the Aegean region.
These tectonic structures impose a threat to the nearby districts and provinces along the northern Aegean Sea region and western Anatolia. Additionally, the Ayvacık EQ sequence helps to improve the knowledge of the offshore Biga Peninsula and Edremit Gulf seismotectonics and Coulomb stress change. It also may carry important information on the nature of future earthquakes along the Biga Peninsula and near tectonic structures in Mytilini and NW Anatolia.