Microearthquakes preceding a M4.2 Earthquake Offshore Istanbul

A primary hurdle in observing small foreshocks is the detection-limit of most seismic networks, which is typically about magnitude M1-1.5. We show that a start-up test of a borehole-based seismic network with a much lower detection limit overcame this problem for an Mw4.2 earthquake. This earthquake occurred offshore of Istanbul, Turkey, on a fault system that is likely to rupture in an M > 7 event in the coming decades. In the three days before and two after, a total of 62 or more earthquakes, including at least 18 foreshocks, came from the mainshock source area. The signal similarity of the foreshocks shows a clear increase during the hours before the Mw4.2 mainshock. Similar foreshock sequences have recently been reported for a few well monitored M > 7 plate-boundary earthquakes. The sequence surrounding the Mw4.2 gives the impression of stochastic failures that ended up interactively unloading stress concentrations. The Mw4.2 mainshock then resulted from the accumulated release of significantly smaller events, as suggested by other field and laboratory studies.

The TESV station is part of the 7-station Geophysical Observatory at the North Anatolian Fault zone (GONAF). These borehole seismic stations surround the eastern Sea of Marmara 13 . Each of the GONAF sites include vertical and 3-C seismic sensors distributed at ~75 m intervals along 300-m deep boreholes. Consequently, they function as an array with locally very low magnitude-detection threshold -down to M~0. As  14,21 . Stars (M > 6) and years (M > 6.8) mark large earthquakes along the NAFZ since 1912, including the 1999 Izmit and Düzce events 28,40 . During the last century, the entire NAFZ was activated except for the Marmara section offshore Istanbul (indicated by the red line), where the last large (M > 7) earthquake occurred in 1766 16,24,25 . The focus area of the study is the eastern portion of the Marmara section of the NAFZ indicated by the black square (enlarged in b). This figure was created using GMT (Generic Mapping Tools) version 4.5 available at http://gmt.soest.hawaii. edu/projects/gmt/wiki/Download. (b) Eastern Marmara region, GONAF stations and location of the M w 4.2 event. Location of the June 25, 2016, M w 4.2 earthquake (yellow star) below the eastern part of the pull-apart Cinarcik Basin where the NAFZ branches into the Armutlu fault and the Princes Islands segment. Locations of GONAF borehole-geophone arrays are indicated in red 13 . The schematic sketch on the left shows a cross section of the TESV sensor distribution with four levels of 1, 2, and 15 Hz vertical and 3-component geophones. The black rectangle is enlarged in the lower right showing epicenters of seismicity during the last decade (black dots) and the eight strongest events (red dots) of the seismic sequence framing the M w 4.2 mainshock. Epicenters are local seismicity from the preceding decade and determined from the permanent regional seismic network operated by the Turkish Disaster and Emergency Management Presidency of Turkey (AFAD) 31 . This figure was created using GMT (Generic Mapping Tools) version 4.5 available at http://gmt.soest.hawaii.edu/projects/gmt/ wiki/Download. a network for regional low-noise monitoring of the eastern Marmara target area it allows for accurate hypocenter locations down to M~1. Cross-correlating TESV waveforms of all 62 events suggests spatial clustering of the entire sequence within an area of 1 km² -about the source size of the M w 4.2 mainshock. Calculating running averages of the cross-correlation coefficients shows a well-defined increase during the hours before the mainshock -reminiscent of the lab results referred to above.

Study Region and Data Base
The NAFZ separates the Anatolian and Eurasian plates, extending for 1200 km between the Karliova triple junction in eastern Anatolia and the Gulf of Saros, Northern Aegean Sea [14][15][16] (Fig. 1a). The westward movement of Anatolia has developed in the tectonic framework of the northward moving Arabian plate 17,18 . It is connected to southward rollback of the Hellenic subduction zone, where the African lithosphere is subducted below the Aegean plate 19,20 . It has an average GPS slip rate of 20-25 mm/yr, increasing towards its western end 17 . A dominantly strike-slip fault zone along the bulk of the NAFZ turns into a more complex transtensional system of fault branches in NW Turkey 15,21 .
Starting in 1939, seven M~7 earthquakes occurred between 1939 to 1999 whose epicenters progressed sequentially westward along the NAFZ, arriving at the eastern end of the Sea of Marmara in 1999 15,[22][23][24][25] (Fig. 1a). The 1999 M w 7.4 Izmit and M w 7.1 Düzce mainshocks resulted in the death of >20.000 persons 26,27 . This left the Marmara section as the only segment that has not produced a large earthquake since 1766. The average recurrence rate on this section is around 200-250 years 16 . The cumulative moment release in the Marmara region has doubled since the 1999 events 25 . Nonetheless activity along the main fault branch below the Sea of Marmara is  Table for  sparse. Few M > 4 earthquakes have occurred, and several aseismic fault patches were identified that could serve as nucleation points for the pending M > 7 event 25,28,29 .
The 5-day TESV-site equipment test captured the seismic activity surrounding the magnitude M w 4.2 earthquake at 05:40:15.18 UTC on 25 June 2016. This includes at least 18 foreshocks that immediately preceded the M w 4.2, the largest regional earthquake in several years.
To test for foreshocks of much lower magnitude than an M w 4.2 earthquake, we need a corresponding earthquake-catalogue completeness magnitude M c . In this case, M c should be ideally on the order of −1, preferably even lower. Burying seismometers 300 m underground can readily lower a networks detection threshold by as much as 2 or more magnitude levels, depending on local conditions. This type of installation results in an order of 5 to 50-fold increase in detection of small earthquakes 12,13,30 . Except for the SHTH site, which for our study period had only a surface station in operation, the data discussed here come from the multi-level arrays reaching this depth (see Fig. 1b). Especially important were the bottom 1 Hz and 2 Hz 3-C sensors, as these recorded clear S-waves, thereby determining hypocentral distances.
The sampling rate of our TESV data was set at 500 Hz. Except for a half hour gap starting at 12:30 UTC on 23 rd June 2016, the sensors of the TESV array recorded continuously from 08:00 UTC on 22 nd June to 08:00 UTC on 27 th June 2016 -a full 5-day span. The M w 4.2 earthquake occurred in the middle of the TESV test data set, at 05:40:15.181 UTC on 25 th June 2016.
The TESV borehole array allowed us to record approximately 7 times the number of events in the M w 4.2 sequence as the combined, surface-based, Turkish national networks. During the 5-day start-up test at least 61 more events with magnitudes ranging from M w = 0 to M w = 3.5 were detected at TESV as coming from the M w 4.2 rupture area. Some of the larger events were also recorded at other GONAF borehole arrays in operation at that time (ESNK and BOZB) and a few of these were also recorded at the island-based SHTH surface sensor. Including the mainshock, a total of 9 events -3 before the M w 4.2 -were strong enough to be located by national seismic networks 21 (red dots in the inset of Fig. 1b).

Methods and Results
The 62 earthquakes were located in an area with diffuse background seismicity in the preceding decade 31 (black dots in the inset of Fig. 1b). To compare the space-time relationships and waveforms of these events to the background activity, we applied both statistical and signal processing methods. Our statistical analysis compared inter event times and distance for the 62 events with all other events surrounding the epicenter as listed in the Turkish national catalogue of AFAD 31 . According to the ANOVA analysis of variance test 32 , their clustering is significantly different from the background seismicity in the same area during the rest of 2016.  We filtered out electrical and low-frequency seismic noise from our seismograms by applying a fourth order Butterworth band pass filter between 3 and 45 Hz. We also inspected the spectrogram of these recordings as a function of sensor depth 13 . In this way, we visually identified 110 earthquakes in the five days of data. Of these, we found a total of 61 events with sufficiently high signal-to-noise ratio to accurately measure S-P differential arrival times. With these data it was possible to pick S-P times relative to that of the M w 4.2 to within 0.1 s or less. The average S-P time of all 62 events is 2.00 s +/− 0.07 s. Their waveforms as recorded at the 290-m deep 1 Hz vertical seismometer of the TESV array are shown in Fig. 2.
To assign magnitudes for the 62 events, we compared TESV P-wave amplitudes to the nine events for which magnitudes were determined by the national network. These magnitudes ranged between 4.2 and 1.2. Aligning the nine reported magnitudes against their amplitudes recorded at the TESV, we determined the relation M = log(A TESV ) − 1.74 and used it to estimate the magnitudes of the remaining events (Fig. 3a). This resulted in magnitudes for the entire sequence ranging from 0.1 to 4.2 (see Table 1). Plotting the magnitudes in chronological order shows a very rough trend toward increasing values leading up to the M w 4.2 mainshock. This is followed by a typical aftershock sequence, with the largest aftershock being about one magnitude step (M ~ 3) smaller than the mainshock (Fig. 3b).
We do not have enough recordings with good azimuthal coverage for accurately determining the hypocenters of each of the foreshocks. We were, however, able to estimate epicentral distances using S-P differential travel times recorded at the GONAF borehole stations BOZB, ESNK, and TESV. S-P data was also available from the   Shown from left to right are the vertical seismogram, the EW-depth, and the NS-depth particle motion, respectively. (d) Relative back azimuths for events with magnitudes M > 1 as determined from the polarization analysis based on particle motions. A consistent value of 145° for the back azimuth is obtained for the larger events while the scatter starts to increase for the smallest events due to reduced signal-to-noise ratios. Squares, crosses and circles indicate foreshock, mainshock and aftershocks, respectively. The color is encoded with the rectilinearity value obtained for each event. The grey rectangle frames one standard deviation of the azimuths obtained for the events with M > 1.3. The general indication is that due to the consistent back azimuth together with the uniform S-P differential travel time the entire 62-event sequence presented in this study could originate from the same fault patch that was activated during the M w 4.2 mainshock.  SHTH surface station on the Princes Island Sivriada (Fig. 4). The epicentral circles are estimated to have a radial precision of about 350 m. Their intersections are concentrated around the epicenter associated with the mainshock. Based on the scaling relations of Bohnhoff et al. 33 they cover a patch on the order of 1 km². The projection of these circles to the estimated hypocentral depth (11 km) of the M w 4.2 further reduces the space containing the 61 additional events.
To help constrain this source volume, a polarization analysis to the sequence's P wave particle motions was done. For this study, the 2 Hz 3-C borehole sensors of the TESV array was used. Stable back-azimuths were found for M > 1 events (Fig. 5). The consistent back azimuth together with the uniform S-P differential travel times for the 62-events suggests that they all originated from the same 1-2 km long zone that failed during the M w 4.2 mainshock.
To quantify their similarity, we cross-correlated event waveform pairs as recorded on TESV's 290 m deep 1 Hz vertical sensor. We de-trended and tapered the waveforms in 4.1 s long windows, starting 0.1 sec prior to the P-wave onset. These windows thus included both the P-and S-wave arrivals and their codas. Their cross-correlation maxima were arranged in a time sequential, square, 62 × 62-element, matrix with their autocorrelations lying along its diagonal axis 34 . The values of the resulting 1,891 coefficients range between 0.08 and 0.92 (Fig. 6a). The resulting matrix contains a high-correlation sub-matrix of events just before the M w 4.2, and a less extensive one during the aftershocks. Running averages of the coefficients with varying window lengths shows their time-dependent trends. The averages were calculated in a retrospective manner: the results include only events prior to the time point shown (Fig. 6b). These averages increased about 20 hours prior to the mainshock, reaching a maximum 10 minutes before the 4.2 mainshock (Fig. 6c). The results also illustrate how our correlation method might be implemented in an earthquake forecasting system involving real-time signal processing and waveform cross-correlation.
The events at the beginning of the sequence and those towards the end show larger differences in S-P times and smaller cross-correlation coefficients (Fig. 7). This indicates larger spatial differences in their hypocentral location. In contrast, the events surrounding the mainshock show the smallest S-P differences and highest cross-correlations, suggesting their close spatial relationship with the M w 4.2 earthquake.

Conclusions
We report on the borehole-based detection of a 62-event sequence framing an M w 4.2 mainshock along the Marmara seismic gap offshore of Istanbul in NW Turkey. The similarity of the earthquakes reported here indicates that they occurred within several hundred meters of the mainshock -in other words within its estimated source area 33 . We found three lines of evidence supporting this conclusion. First, the epicentral circles determined by S-P times at four GONAF sites intersect within a few hundred meters of the M w 4.2 location. Second, the Pwave particle motions of best resolved M > 1 events consistently point to the same back azimuth as the M w 4.2. Third, the waveforms of these events are more strongly correlated in time and space than other events in the 1 st and 2 nd halves of 2016.
In part we observe a set of foreshocks with increasing waveform similarity during the hours before the M w 4.2 earthquake that show a similar behavior as foreshocks repeatedly observed during laboratory rock deformation test and more recently before large plate-bounding earthquakes.
The exact locations of the events discussed here could not be fully constrained since their magnitudes were too small for them to be registered at the regional surface-networks or the other operational GONAF stations. Cross-correlation coefficients for event pairs plotted with differential S-P time. Cross-correlation coefficients for all event pairs plotted with their difference in S-P time. Orange and red dots show event pairs for early (events 1-11) and immediate (events 12-18) foreshocks. Blue and gray dots show event pairs for immediate (events 20-23) and later (from event 24) aftershocks. The event pairs with the highest crosscorrelation coefficients also have the smallest S-P time differences indicating that they could origin from the same fault patch activated by the mainshock. This is interpreted to reflect the emergent failure process leading to the M w 4.2 mainshock. The lack of even more sensitive or nearer stations precluded testing, for example, the inverse Omori law for foreshocks 35 , where the rate of earthquakes before a mainshock increases according to a power law. The same applies for both (1) a decrease in b-values, as posited for near-offset events preceding a mainshock and (2) for the migration of foreshocks toward the mainshock. Well-documented field evidence for foreshock behavior in nature is still sparse. Our example is one of only a few field-based observations of much smaller, near-hypocenter seismicity preceding a mainshock. Consequently, there are a number of explanations for the relationship of foreshock to mainshocks. These include, for example, the controversial cascade model where earthquakes trigger aftershocks larger than themselves 5,36 . Foreshock laws are partly still seen as statistical in nature, observable when averaging over a large number of sequences, but not systematically for every event 35 .
In this light, the progress in monitoring instrumentation and its operation in a seismically active region that is discussed here may hold promise for hazard and risk reduction. Going underground with seismic monitoring, as it were, can work even in highly urbanized areas. Adding to this technology the sort of analysis methods we present here can contribute towards refining operational earthquake forecasting, even perhaps helping plan activities such as critical evacuations [37][38][39] .