Crustal imaging across the North Anatolian Fault Zone from the autocorrelation of ambient seismic noise

Seismic images of active fault zones can be used to examine the structure of faults throughout the crust and upper mantle and give clues as to whether the associated deformation occurs within a narrow shear zone or is broadly distributed through the lower crust. Limitations on seismic resolution within the crust and difficulties imaging shallow structures such as the crust‐mantle boundary (Moho) place constraints on the interpretation of seismic images. In this study we retrieve body wave reflections from autocorrelations of ambient seismic noise. The instantaneous phase coherence autocorrelations allow unprecedented ambient noise images of the North Anatolian Fault Zone (NAFZ). Our reflection profiles show a Moho reflected P wave and additional structure within the crust and upper mantle. We image a distinct vertical offset of the Moho associated with the northern branch of the NAFZ indicating that deformation related to the fault remains narrow in the upper mantle.


Introduction
The behavior of fault zones within continental settings is more complex than simple plate tectonics would predict. While deformation in oceanic environments is typically localized on plate boundaries, deformation within continents occurs across regions that can be thousands of kilometers wide (e.g., the Alpine-Himalayan Belt). It has been proposed that ductile flow within the lower crust causes distributed deformation that otherwise occurs within narrow shear zones in the continental upper crust and upper mantle [Bürgmann and Dresen, 2008]. The formation and evolution of major continental transform faults, like the North Anatolian Fault, appears to be a delicate balance between accommodation of the dominant tectonic strain field and the exploitation of weak preexisting geological features such as suture zones.
In particular, understanding how structural changes throughout a fault zone can be used to infer the mechanism of strain localization at lower crustal depths has important implications for our knowledge of lithospheric rheology and our understanding of the earthquake cycle. Seismic imaging of active fault zones has provided two contrasting arguments: some authors, such as Wilson et al. [2004], argue that strain is more diffuse in the lower crust and upper mantle, forming a broad shear zone and citing pervasive seismic anisotropy in the lower crust and, in the case of the Marlborough Fault Zone, the lack of visible offsets in the Moho discontinuity beneath the surface signature of the fault.
In contrast, other studies on continental strike slip faults [Zhu, 2000;Weber et al., 2004] argue that variations in crustal thickness, correlated with the location of active faults, show that the Moho has been offset during faulting, therefore showing that strain may still be localized in a narrow band of ductile deformation to at least upper mantle depths. Receiver function analysis of the Tibetan Plateau [Zhang et al., 2014] has also revealed apparent Moho depth offsets that correspond to surface traces of major faults and sutures, in the region of the Earth's thickest continental crust. Their images suggest that narrow zones of localized strain structures can reach to at least upper mantle depths, even though crustal thickening and anisotropy imply that diffuse strain is also important.
Here we aim to image the crustal structure of the North Anatolian Fault Zone (NAFZ) in Turkey in the region of the 1999 Izmit and Düzce earthquakes. The NAFZ is a ∼1500 km long right-lateral strike-slip system along which the Anatolian Block adjoins Eurasia. The NAFZ formed in the middle Miocene and propagated westward to reach the Sea of Marmara ∼200 ka ago [Sengör et al., 2005] and extends now into the North Aegean trough. Since 1939 the NAFZ has displayed a westward migrating pattern of large earthquakes (M w ≥ 7), resulting most recently in the magnitude 7.6 and 7.2 Izmit and Düzce earthquakes of 1999 [Barka et al., 2002] and continues to pose significant seismic hazard to the region. In this region, the NAFZ divides into a northern and southern strand [Karimi et al., 2014], with the northern strand showing more seismic activity [Altuncu Poyraz et al., 2015]. The NAFZ appears to follow the trend of accretionary complexes (Figure 1) formed during the closing of the Tethys Ocean [Sengör et al., 2005]. The NAFZ roughly follows the Intra-Pontide Suture between the continental Istanbul Zone and the accretionary Sakarya Terrane in the south [Okay, 2008]. In between the northern and southern strands of the NAFZ lies the Armutlu-Almacik block (Figure 1), which formed in a similar fashion to the Sakarya Terrane, as a set of accretionary complexes during the closure of the Tethys Ocean [Sengör et al., 2005;Bekler and Gurbuz, 2008].
Between May 2012 and October 2013, the University of Leeds, Kandilli Observatory and Earthquake Research Institute (KOERI) and Sakarya University installed a 66 station network of 3-component seismometers in a roughly rectangular grid ( Figure 1) with an approximate station spacing of ∼7 km known as the Dense Array for Northern Anatolia [DANA, 2012]. DANA was supplemented by three permanent stations of the KOERI network (GULT, SAUV, and SPNC) and seven stations to the east of the array forming a semicircle, with larger station spacing. Each northward trending line of DANA is referred to by a code DA-F, with DA being the westernmost, and F being the easternmost line. The network was deployed across the rupture zone of the 1999 Izmit earthquake [Barka et al., 2002] and crossed both the northern and southern strands of the NAFZ to sample all three crustal terranes: the Istanbul Zone to the north, the Armutlu-Almacik block bounded by the two fault strands, and the Sakarya Terrane to the south ( Figure 1).

Data and Methods
In order to image the structure of the NAFZ throughout the entire crust and upper mantle, we constructed stacked autocorrelograms of continuous seismic records from DANA. It has been shown that the autocorrelation of surface displacement at a seismometer yields information on the reflection response of structure near the receiver, scaled by the spectrum of the sources exciting the ground motion and by amplitude of the waves arriving at the structure beneath the instrument [Gorbatov et al., 2013]. Thus, by autocorrelating continuous records of ground motion, it is possible to obtain Earth's reflection response at the seismometer location and to image crustal structure, including free surface reverberations, for a colocated source and receiver [Claerbout, 1968]. This technique has been applied to image crustal structure beneath the continental United States [Tibuleac and von Seggern, 2012] and Australia [Gorbatov et al., 2013;Kennett et al., 2015;Kennett, 2015] and offers complementary information to the widely used receiver function technique. We constructed autocorrelograms for all stations of the DANA network by splitting the data set into 1 h long time windows of the vertical component of ground motion. In contrast to previous studies using autocorrelation interferometry [Tibuleac and von Seggern, 2012;Gorbatov et al., 2013], here we employ a phase autocorrelation method to enhance the signal-to-noise ratio, and ensure that coherent signals such as body wave reflections are retrieved independent of amplitude [Schimmel, 1999].
Phase correlations have been rigorously compared with amplitude correlation methods by Schimmel [1999] and Schimmel et al. [2011]. Schimmel et al. [2011] utilized phase correlations to obtain empirical Green's functions from seismic array data at a local array in southern Spain, as well as the GEOSCOPE global seismic network. Schimmel et al. [2011] noted close correspondence between the retrieved Green's functions obtained using phase and amplitude correlation, with those computed through phase correlation showing higher signal-to-noise ratio. The correspondence between phase and amplitude correlation implies that we can use phase autocorrelation to retrieve the reflection response of the Earth, even though the mathematical basis for this correspondence may be unclear. To calculate the phase correlograms, the 1 h data segments are Hilbert transformed and the analytical signal (a complex valued time series in which the amplitude represents the envelope of the seismogram, and the argument the instantaneous phase) is calculated [Bracewell, 1965]. From the analytical signal, the time series of instantaneous phase is extracted, and these records are autocorrelated using [Schimmel, 1999]: where C pac (t) is the phase autocorrelation, is instantaneous phase, and is the correlation lag time. The autocorrelation windows are then linearly stacked in order to enhance signal-to-noise ratio. Phase correlation has the advantage that it is amplitude independent and that coherent signals such as crustal reflections generated from the ambient seismic wavefield can be better retrieved than through an amplitude correlation, which is often dominated by high-amplitude surface waves [Schimmel et al., 2011]. The phase correlation also removes the need to remove coherent signals such as earthquake arrivals through, e.g., temporal normalization [Bensen et al., 2007], as the instantaneous phase of arrivals from such sources is amplitude independent and of comparable magnitude to arrivals from transient low-energy sources such as the ambient noise field. Thus, by stacking autocorrelations over a sufficiently long period of time, these transient energetic signals are rendered negligible in the final reflection image.
From the stacked phase autocorrelations it is possible to produce images consisting primarily of subsurface reflections, in contrast to many studies based on seismic interferometry that analyze the surface waves that dominate the ambient wavefield [e.g., Shapiro et al., 2005;Ren et al., 2013, andEkström, 2014].

Results
It is only possible to retrieve reflections from interfaces representing relatively strong contrasts in material properties, due to the low amplitude of body wave reflections in the ambient noise field [Gorbatov et al., 2013]. Our autocorrelogram profiles show strong evidence for Moho reflections and evidence for internal crustal reflectors ( Figure 2). We have visually interpreted the reflections seen in our autocorrelation profiles, and a sketch of Earth structure thus inferred beneath the DA and DF lines of DANA is shown in Figure 3. All autocorrelograms are band-pass filtered between 0.2 and 0.4 Hz. Analysis of noise spectra show that this band contains high power from sources such as secondary ocean microseisms [Stehly et al., 2006], and teleseismic earthquakes [Seriff et al., 1965;Gupta, 1965]. Higher-frequency noise sources would be able to increase vertical resolution of this approach, but the ambient noise field recorded at DANA contains very little power at frequencies ≥ 1 Hz.
A coherent reflection is visible in all autocorrelograms at ∼12 s, in close correspondence with the theoretical arrival time of a P wave reflection from the Moho (PmP) at 39 km depth (12.53 s) as calculated for the crustal model of Karahan et al. [2001], which we interpret as the Moho reflection. The arrival time of this reflection has a maximum variation between stations of ± 1 s, which is likely the combined result of variations in P wave velocity structure in the crust, and topography on the Moho. Kahraman et al. [2015] and Frederiksen et al. [2015] utilized P receiver functions and a transfer function method, respectively, to estimate crustal thickness across DANA as varying from 32 km in the Sakarya Terrane up to 45 km in the Istanbul Zone.
Evidence of shallower structure throughout the crust is also visible. Prominent arrivals are visible at ∼2.5 s and ∼7 s in all autocorrelations. Both of these arrivals display variation in arrival times on the order of 2 s. The arrival at ∼2.5 s may correspond to a reflection from a crystalline basement, which has been found beneath the Istanbul Zone [Okay, 2008], overlain by Tethyside accretionary complexes. An arrival is also apparent at ∼7 s, indicating a potential middle crustal reflector. This reflector may correspond to an interface at 20 km depth also evident in receiver function images of structure beneath the DANA array [Kahraman et al., 2015]. The arrival at ∼7 s also appears to be systematically offset to a deeper location within the Istanbul Zone, particularly in profiles sampling the western half of the study area (DA, DB, and DC).
A further arrival is visible at ∼5 s. This arrival is most prominent at stations located on the Istanbul Zone; however, it is also visible at some stations within the Sakarya Terrane but not within the Armutlu Block. This arrival generally shallows toward the north of the array, arriving on average ∼1 s earlier in the Istanbul Zone than in the Sakarya Terrane. This northward shallowing reflector could be indicative of a southward dipping reflector or of faster seismic velocities in the upper crust of the Istanbul Zone.
Finally, a coherent arrival is also seen below the Moho at ∼17.5 s. The arrival times of this phase vary by ±1 s. A reflector at this depth has also been identified in wide-angle refraction studies in the continental United States [Green and Hales, 1968], as well as being identified in a previous reflection study using ambient noise in northern Africa [Ruigrok et al., 2011]. The autocorrelation method may have the potential to retrieve reflections from deeper upper mantle interfaces, but signal-to-noise ratio appears to decrease with depth.

Discussion
The clear P wave Moho reflection in Figure 2 shows the imaging potential of the autocorrelation technique. The Moho reflection is strong and coherent in the Sakarya Terrane and the Armutlu Block but is less clear within the Istanbul Zone (e.g., lines DB and DF and north of latitude 40.7).
Our results show apparent topography on the Moho, though Moho depth variation may trade off against average crustal velocity. By varying the seismic velocities of the crustal model of northwestern Anatolia derived from seismic refraction [Karahan et al., 2001] by ± 5%, we calculate an upper limit of ± 0.6 s of travel time variation can be attributed to lateral variations in seismic velocities. This leaves ∼50% of our observed travel time variations that could be attributed to topography on reflectors. Notably, Kahraman et al. [2015] detect a deepening of the Moho from west (DA, DB, and DC) to east (DE and DF) within the Istanbul Zone. Figure 2 shows that we also observe a delay in the Moho reflection in the northern part of profiles DE and DF, possibly indicating a deeper Moho toward the East.
We note a significant offset and disruption of the Moho reflection associated with the northern strand of the NAFZ, between the Istanbul Zone and Armutlu Block (station latitude north of ∼40.8). The Moho reflection is visible as a strong arrival in every autocorrelogram but disappears for stations close (within 7 km) to the northern NAFZ ( Figure 2). The weak reflection in this region may indicate that deformation associated with the active fault has reduced the reflectivity of the crust-mantle boundary to the point that no reflections are generated. A potential mechanism for the reduced reflectivity of the upper mantle beneath the northern branch of the NAFZ is the serpentinization of the upper mantle, due to the percolation of fluids down through the shear zone. Serpentinization would reduce the seismic velocity of upper mantle rocks, causing the velocity contrast at the Moho to become more gradational, reducing the reflectivity. This process has been previously invoked to explain the absence of a Moho discontinuity beneath the Mid-Hungarian Shear Zone in a receiver function study of the Pannonian Basin [Hetényi et al., 2015].
The autocorrelation images show a variety of crustal reflections of varying amplitude. The phase correlation technique is sensitive to coherence, so a lower amplitude arrival (e.g., the Moho arrival within the Istanbul Zone) is indicative of reduced reflection coherence. The loss of coherence could be due to a loss of reflectivity but could also be due to the existence of small scale structure at the reflector. A further consequence is that the Moho, which would be expected to be the most prominent discontinuity within the crust and upper mantle, is not higher in amplitude than any other reflector in our profiles. This phenomenon can be attributed to the amplitude independence of phase correlation.
We note that the region of absent Moho reflections is narrow, suggesting fault-related strain in the lower crust may be localized to within 10 km of the NAFZ. The diminished strength of the reflected signal is consistent with receiver function analysis of the DANA data set [Kahraman et al., 2015] that shows a clear reduction in amplitude of P to S conversions at the Moho in the vicinity of the northern strand of the NAFZ. The receiver 10.1002/2016GL067715 functions also show a somewhat less clear Moho amplitude reduction near the southern strand that is not observed in our study. We therefore infer that the damage zone of the northern branch of the NAFZ at least reaches the Moho and P wave reflections from the boundary, as well as P to S conversions, are suppressed by a more gradual increase of velocity with depth across the Moho. The restricted spatial extent of the effect supports the contention that localized strain is occurring at lower crustal depths [Weber et al., 2004]. Below the crust, full-waveform inversion of the velocity field by Fichtner et al. [2013] has also shown that the mantle lithosphere beneath the NAFZ is generally marked by a relatively narrow band of low shear velocity that extends down to the asthenosphere and presumably defines a locus for more focussed deformation at crustal levels.
Our observation of offsets in Moho depth associated with the northern branch of the NAFZ indicates that it is likely that localized shear associated with the NAFZ has reached the upper mantle. Narrow shear zones that extend into the continental mantle have also been inferred at other large continental fault systems, such as the Dead Sea Transform [Weber et al., 2004], and within the Tibetan Plateau [Zhang et al., 2014]. The apparent continuity of the Moho reflection across the southern strand of the NAFZ (Figure 2) might then be interpreted as implying less strain accumulated on this structure.
Below the Moho, the autocorrelograms show a further coherent arrival at ∼17 s for most stations, with travel time variations on the order of 1 s. A similar arrival has been detected in autocorrelograms by Tibuleac and von Seggern [2012] in Nevada and has been interpreted as an S wave arrival from the Moho (SmS). Vertical component autocorrelograms, however, should not be sensitive to vertically traveling S waves, so we interpret this arrival as a P wave reflection off an upper mantle discontinuity, likely the Hales discontinuity, which has previously been attributed as a phase transition from spinel to garnet in aluminous peridotite [Hales, 1969].
An alternative explanation for the Hales discontinuity is that it is the upper limit of a zone of enhanced seismic anisotropy due to mineral fabric beneath the lithospheric mantle lid that may appear as a seismic discontinuity to low frequency seismic waves [Levin and Park, 2000]. A potential mechanism for the formation of this fabric is the westward movement of Anatolia [Reilinger et al., 2006;Nocquet, 2012] as it "escapes" toward the Aegean Sea from the collision of Arabia with Eurasia [Jackson and McKenzie, 1984]. The motion of the strong lithospheric lid over a weaker asthenosphere may have the potential to induce a coherent fabric in olivine crystals, forming the Hales discontinuity. Such continental escape processes have been invoked to explain the appearance of the Hales discontinuity in Arabia [Levin and Park, 2000] and the Basin and Range Province of the Sierra Nevada Mountains [Jones and Phinney, 1998].
Stations located in the Istanbul Zone show a strong reflector at ∼5 s. It appears to terminate at the location of the surface trace of the northern strand of the NAFZ. The Intra-Pontide Suture between the Istanbul Zone and Armutlu Block provides a simple explanation of discontinuous structure. This arrival is also visible at some stations within the Sakarya Terrane in the south but is not detected within the Armutlu Block. Receiver function analysis also reveals a strong discontinuity at a depth of 15 km within the Istanbul Zone, also terminating at the northern branch of the NAFZ [Kahraman et al., 2015]. It is likely that the reflection visible at ∼5 s corresponds to this discontinuity. While this structure is observed in receiver functions to be stronger in the eastern part of the study region, our autocorrelograms show a reflector that is equally visible throughout the entire Istanbul Zone.

Conclusions
We have shown that it is possible to extract images of body wave reflections by constructing stacked autocorrelograms of ambient seismic noise. The dense station spacing of DANA enables high lateral resolution of the crustal structure in the vicinity of the NAFZ. While these autocorrelograms are only sensitive to strong contrasts in physical properties in the subsurface, it is still possible to image crustal discontinuities such as the Moho using this technique. The dense station spacing of DANA greatly aids in identifying and interpreting coherent reflections from subsurface interfaces beneath the array. We have shown that it is possible to use microseismic noise sources in the frequency band 0.2-0.5 Hz in order to construct these autocorrelograms, negating the need for energetic sources such as earthquakes. Our reflection profiles show that similar to receiver functions [Kahraman et al., 2015], P wave reflections from the Moho fade out beneath the northern NAFZ. In addition, offsets in Moho depth are visible within ∼10 km of the northern strand of the NAFZ, while the southern strand has a less clear effect on the Moho. We interpret that deformation associated with the northern branch of the NAFZ has at least reached the Moho, in the form of a narrow shear zone. These indications of localized Moho variation beneath the NAFZ are consistent with localized low-velocity regions in the lithospheric mantle beneath [Fichtner et al., 2013], as well as observations from other large continental fault zones, such as the Dead Sea Transform [Weber et al., 2004], San Andreas Fault [Zhu, 2000], and Tibetan Plateau [Zhang et al., 2014], which show narrow bands of deformation associated with active faults which extend downward at least as far as the Moho. In contrast, the southern branch of the NAFZ either does not cut through the entirety of the crust or the deformation in the upper mantle likely occurs by a more diffuse shear zone, as interpreted for the Marlborough Fault System, New Zealand [Wilson et al., 2004]. Our observations demonstrate that the autocorrelation method may provide reliable estimates of important parameters such as crustal thickness and the ability to map heterogeneity in the crust and upper mantle independently of other established seismological techniques, such as receiver functions.