Structure of the northwestern North Anatolian Fault Zone imaged via teleseismic scattering tomography

S. Rost ,1 G.A. Houseman,1 A.W. Frederiksen,2 D.G. Cornwell ,3 M. Kahraman,4,5 S. Altuncu Poyraz,5 U.M. Teoman,5 D.A. Thompson,1,* N. Türkelli,5,† L. Gülen,6 M. Utkucu6 and T.J. Wright7 1School of Earth and Environment, The University of Leeds, Leeds, LS29JT, United Kingdom. E-mail: s.rost@leeds.ac.uk 2Department of Geological Sciences, University of Manitoba, Winnipeg, Manitoba, Canada 3Department of Geology and Geophysics, School of Geosciences, University of Aberdeen, King’s College, Aberdeen, AB243UE, United Kingdom 4Eurasia Institute of Earth Sciences, Istanbul Technical University, Istanbul, Turkey 5Kandilli Observatory and Earthquake Research Institute, Department of Geophysics, Boğaçizi University, 34684 Cengelköy, Istanbul, Turkey 6Department of Geophysical Engineering, Sakarya University, Esentepe Campus, 54187 Sakarya, Turkey 7COMET, School of Earth and Environment, The University of Leeds, Leeds, LS29JT, United Kingdom


I N T RO D U C T I O N
The North Anatolian Fault Zone (NAFZ) is one of the longest continuous continental strike slip fault systems on Earth, posing considerable hazard to Northern Anatolia and especially the mega-city of Istanbul towards its western end. Here, we use a novel imaging Scattering tomography NW NAFZ 923 Barka 1992;Şengör et al. 2005;Reilinger et al. 2006) driven by the gradient of gravitational potential energy from the Anatolian plateau to the Hellenic Trench (England et al. 2016).
The NAFZ ruptured in a series of M > = 6.7 earthquakes during the 20th century from east to west (Stein et al. 1997), interpreted as stress transfer along the strike of the NAFZ from one earthquake bringing the next segment closer to failure. The two most recent events in the current series occurred in 1999 with epicentres iṅ Izmit (M = 7.6) and Düzce (M = 7.2, Barka et al. 2002;Gülen 2002), with the fault rupture extending into the Sea of Marmara and the next anticipated event in the series posing a pronounced risk to the city of Istanbul.
While the deformation at the surface is localized on faults (Bürgmann & Dresen 2008;Hussain et al. 2016), the distribution of deformation throughout the crust and into the mantle remains unclear (Bürgmann & Dresen 2008;Vauchez et al. 2012;Moore & Parsons 2015). Understanding the structure and dynamics of fault zones, especially at depth, is essential for our understanding of continental deformation and seismic hazard.
Here, we aim to better understand the structure of the NAFZ, especially in the middle and lower crust and into the upper mantle, using data from temporary seismic stations and exploiting the scattered seismic wavefield following the P-wave arrivals of teleseismic events (Frederiksen & Revenaugh 2004). We use data from the 18-month DANA deployment (DANA 2012) across the NAFZ in the region of the 1999 ruptures (Fig. 1a). The P-wave coda contains energy from P-to-P and P-to-S scattering at small-scale heterogeneities along the ray-paths. Structure can be recovered from the scattered seismic energy through migration approaches ranging from common-conversion-point or common-scattering-point stacking (e.g. Dueker & Sheehan 1997) to full depth migration (e.g. Ryberg & Weber 2000). Here we are using a tomographic waveform approach based on linear inverse theory of the scattered wavefield (Ji & Nataf 1998;Frederiksen & Revenaugh 2004) to resolve the structure of the lithosphere and potential shear zones. The scattering tomography builds on different principles than the more established P-receiver function method previously applied to this data set (Kahraman et al. 2015). With different resolution, limitations and trade-offs this study will contribute to our knowledge of the fine-scale structure of the lithosphere beneath this major fault zone.
We find that the two strands of the NAFZ evident in the shallow structure coincide with main interfaces and interface terminations throughout the crust and into the upper mantle indicating that the fault zone structure may extend to depths of at least ∼75 km in this region. We find evidence for small-scale variation of structure in the vicinity of the strands that might indicate the detection of heterogeneity related to past deformation along the present day fault.

T E C T O N I C S A N D P R E V I O U S G E O P H Y S I C A L S U RV E Y S
The study region (Fig. 1b) is an amalgam of continental and subduction-related oceanic fragments that remain after the closing of the Tethyan Ocean in the late Tertiary (e.g. Robertson & Ustaömer 2004;Okay et al. 2008). The DANA network samples three tectonic blocks situated in the study region (i) The Istanbul-Zonguldak Zone (IZ) to the north of the northern strand of the NAFZ, (ii) the Sakarya zone (SZ) to the south of the southern strand and (iii) the Armutlu and Almacık blocks (AA). The NAFZ splays into a northern (NNAFZ) and southern (SNAFZ) strand west of about 30.65 • , with the NNAFZ and the SNAFZ following the northern and southern boundary of the AA. The AA is interpreted as the cause for the splay (e.g. Akbayram et al. 2016). The NAFZ is colocated with the intra-Pontide suture for at least half its onshore length and the boundaries between the SZ, AA and IZ likely represent the suture in this locale (e.g. Okay et al. 2008). The NAFZ is believed to have developed ∼11 Ma in eastern Anatolia with strain localization propagating westward and reaching the Sea of Marmara, and our study region, before 3.9 Ma (Akbayram et al. 2016). Slip on the northern and southern NAFZ strands has been estimated to be approximately 16-25 and 5-19 mm yr -1 , respectively (Stein et al. 1997;Meade et al. 2002;Flerit et al. 2003). The northern branch of the NAFZ in our study area last ruptured in the 1999İzmit earthquake (Tibi et al. 2001; Barka et al. 2002) and still shows active slip at the surface (Hussain et al. 2016) although current seismicity is not focussing beneath either fault strand (Altuncu Poyraz et al. 2015).
A study using local earthquake waveforms (Horasan et al. 2002) finds a Moho depth of 32 km in the Marmara region. They find crustal discontinuities at 4 and 17 km depth with upper mantle velocities of 8.0 and 4.6 km s -1 for V P and V S , respectively and mantle densities of 3.4 g cm -3 . Upper crustal velocities are found to be 5.8 and 3.4 km s -1 , lower crustal velocities of 6.2 and 3.6 km s -1 and near-surface velocities of 3.5 and 2.2 km s -1 for V P and V S , respectively (Horasan et al. 2002).
Using fault zone head waves and fault zone reflected waves from seismicity and stations close to the fault zone it is possible to resolve the interface properties of the fault zone. Ben-Zion & Sammis (2003) and Ben-Zion et al. (2003) imaged a shallow fault zone extending about 3-4 km into the crust with velocity reductions of up to 50 per cent but a width of only 100 m, well below the resolution of most other seismic methods. In the area of the 1999 Düzce andİzmit earthquakes Bulut et al. (2012) and Najdahmadi et al. (2016) image the material properties across the fault zone using trapped waves  (Farr et al. 2007). Stations are indicated by yellow circles (permanent stations in red). Mapped faults (red lines) provided by Emre et al. (2018) and rupture of the 1999İzmit and Dücze earthquakes (yellow) provided by Gülen (2002). Dashed north-south and east-west lines indicate location of depth profiles shown in Figs 6 and 7 and are approximate locations of depth profiles provided by Kahraman et al. (2015). (b) Simplified geological map of the region outlining the three main tectonic blocks and geological areas. After Taylor et al. (2019a). and find a bimaterial interface down to the base of the seismogenic crust with an average velocity contrast of 3.4 per cent (Najdahmadi et al. 2016) and 6 per cent with the southern block being fast (Bulut et al. 2012).
Studies of local earthquakes detect upper crustal anisotropy around the NAFZ, limiting it to the upper 8 km (Hurd & Bohnhoff 2012) or 4 km (Peng & Ben-Zion 2004), likely due to aligned cracks in the uppermost crust. Detected splitting times are on the order of 10 ms not likely influencing the results of this study. SKS analysis detects asthenospheric anisotropy (Biryol et al. 2010;Paul et al. 2014;Legendre et al. 2021) with lag times between fast and slow direction of typically 1.5 ± 0.4 s with fast polarization directions smoothly varying from NNE-SSW in northern Turkey to NE-SW in eastern Turkey. A combined study of shear-wave splitting and anisotropic receiver functions show complex anisotropy especially in our study region (Lemnifi et al. 2017). Shear wave splitting tomography using aftershocks of the 1999İzmit and Düzce earthquakes show a 3 km wide anisotropic zone extending to 5 km depth with distinct asymmetry being related to damage from the unilateral eastward propagation of the 1999İzmit rupture (Li et al. 2014).
P-wave receiver functions (PRFs) east of the Sea of Marmara indicate a deepening of the Moho from west (29-32 km) to east (34-35 km) (Zor et al. 2003;Vanacore et al. 2013). The average crustal V P /V S in our study region is ∼1.75 (Vanacore et al. 2013). PRFs of the DANA data set (Kahraman et al. 2015) find crustal thickness and V P /V S variation in both EW and NS directions with the crust deepening from 36.5 km (V P /V S = 1.73) to 40 km (V P /V S = 1.73) in the IZ, a constant crustal thickness of ∼ 37 km (V P /V S = 1.69 to 1.70) in the AA, and a slight thinning from ∼ 35 km (V P /V S = 1.73) in the west to ∼ 34 km (V P /V S = 1.85) in the east of the SZ (Fig. 1). Combining data from several permanent stations and temporary station deployments, including DANA data, Jenkins et al. (2020) determined Moho depths across the Sea of Marmara region finding thick crust of up to 41 km in the IZ, with a shallower Moho (32-34 km) in the AA and SZ with evidence of discontinuous structure across the NAFZ. The transition also shows complex Moho structure around the NNAFZ. Additionally, Jenkins et al. (2020) find eastwest variation with a general deepening of the Moho towards the east.
Previous studies evidence strong crustal heterogeneity on scales of less than 10 km with sharp truncations of subhorizontal interfaces coinciding with the surface locations of the northern and southern NAFZ strands. The northern strand seems to penetrate deeper into the crust and may extend into the upper mantle based on analysis of receiver functions and ambient noise cross-correlations (Kahraman et al. 2015;Taylor et al. 2016;Jenkins et al. 2020). Both structural changes in north-south and east-west direction have been reported (e.g. Kahraman et al. 2015;Ç ubuk-Sabuncu et al. 2017;Jenkins et al. 2020).
Using P-wave transfer functions and a grid-search inversion approach (Frederiksen et al. 2015) detected a sharp change of crustal thickness across the northern NAFZ which is believed to follow the trace of the Intra-Pontide suture in this location and a change of the V P /V S ratio across the southern branch indicating a change in basement composition. The IZ shows thick crust (40-45 km) but low topography indicating that it is in isostatic disequilibrium or underlain by thicker lithosphere, a result supported by Jenkins et al. (2020). The transfer functions also provide evidence for thick sediments in the Sakarya and Pamukova basins in agreement with ambient noise analysis (Taylor et al. 2019a).
Tomographic studies using traveltimes from local, regional and teleseismic events (e.g. Bariş et al. 2005;Salah et al. 2007;Koulakov et al. 2010;Bakırcı et al. 2012;Beyhan & Alkan 2015;Polat et al. 2016;Papaleo et al. 2017Papaleo et al. , 2018 and full waveform information (e.g. Fichtner et al. 2013;Govers & Fichtner 2016; Ç ubuk-Sabuncu Downloaded from https://academic.oup.com/gji/article/227/2/922/6318864 by University of Aberdeen user on 30 July 202130 July et al. 2017Blom et al. 2020) find strong velocity contrasts in the crust along the NNAFZ and the SNAFZ. These are interpreted as the fault zone exploiting the sutures between the Istanbul zone, Armutlu block and Sakarya Zone as result of the closure of the Neo-Tethys ocean. Nonetheless, Fichtner et al. (2013) note the lack of a low velocity fault zone signature at depth west of ∼32 • E due to the absence of a well-localized suture and insufficient strain localization due to the young age of the fault zone. In contrast, Koulakov et al. (2010) using local tomography note the juxtaposition of high-velocity, low attenuation blocks (e.g. Armutlu block) to lower velocity areas. A smaller scale, full waveform tomographic study has been performed by Ç ubuk-Sabuncu et al. (2017) noticing strong lateral and vertical velocity variations and strong radial anisotropy of the crust in agreement with the active tectonics of western Turkey. Salah et al. (2007) focus on the rupture area of the 1999İzmit and Düzce earthquakes using local tomography to resolve P-and S-wave velocity as well as the Poisson ratio and detect prominent low velocity zones down to depths of 25 km as well as well as a high-velocity anomaly at a depth of 8 km between 30.0 and 30.4 • E along the strike of the fault zone. Seismicity is more prevalent in the high-velocity region although it also occurs in low velocity regions. Similarly, traveltime tomography resolves narrow subvertical low velocity zones coinciding with the surface expression of the SNAFZ and the NNAFZ (Papaleo et al. 2017(Papaleo et al. , 2018. Magnetotelluric (MT) data show differences in the crustal conductivity from south to north across the NAFZ (Tank et al. 2005) with a high resistivity (≥1000 m) crustal basement in the IZ to the north and a less resistive crustal basement (500 m) in the SZ. The MT data resolve a localized conductive zone (30-50 m) within the AA that extends into the upper mantle that has been attributed to partial melts or pore fluid flow from the upper mantle beneath the NAFZ.

M E T H O D
We apply the teleseismic scattering tomography approach by Frederiksen & Revenaugh (2004) to the DANA data set to resolve the small-scale structure beneath the array. The scattered seismic wavefield is more sensitive to short-wavelength variations in material properties than is the path-integrated sensitivity of transmitted phases such as used in, for example seismic traveltime tomography. The P-to-p and P-to-s scattered energy in the coda of teleseismic P waves travelling along different paths to the main arrival can uniquely determine Earth structure if the sampling of the seismic wavefield is dense enough to avoid spatial aliasing. In the tomographic approach some aliasing can be accepted without introducing issues with non-uniqueness of the solution due to the regularization of the problem. Several approaches to use the scattered coda energy to image the subsurface have been developed, forming a continuous spectrum of method complexity. The common approach of receiver function analysis uses stacked records of P-to-s (or S-to-p) scattered (converted) energy (Vinnik 1977;Langston 1979) which may be binned according to their common conversion point to improve signal-to-noise ratio (Dueker & Sheehan 1997) and mapped to depth. The method assumes a 1-D stratified seismic structure which is often violated in practice (Rondenay 2009). Lateral variation of structure leads to diffraction of the seismic wavefield and diffraction stacking, a backprojection of the diffracted energy along its traveltime hyperbola, can be used to image small-scale perturbations on the order of the seismic wavelength of the structure at depth. These methods are widely used in controlled-source applications (Yilmaz 2001), and are commonly described as migration techniques (Rondenay 2009) but implementation requires dense spatial sampling of the seismic wavefield. General improvements and densification of recent passive seismic deployments make the application of more complex methods, such as traveltime stacking of the scattered wavefield (Revenaugh 1995) or the application of inversion or backprojection operators in a 2-D or 3-D model space (Bostock & Rondenay 1999) possible and allow higher resolution of detail. For a full review of these methods see Rondenay (2009).
For a more complete treatment of the scattering problem, the scattering image problem can be formulated as a tomographic inversion (Ji & Nataf 1998). Using a waveform inversion, Frederiksen & Revenaugh (2004) have developed a linear tomographic inversion of the scattered seismic wavefield which we apply here. A full description is given in Frederiksen & Revenaugh (2004) and we outline only the main points of this approach here.
The Born approximation, a common approximation of the full scattering process, assumes a weakly scattering medium and single scattering (Sato et al. 2012), which is a good approximation for this application where the amplitudes of the scattered wavefield are much smaller than those of the direct wave. In this approximation the scattering properties are represented as perturbations in elastic parameters (δλ, δμ, δρ) to a background model (λ 0 , μ 0 , ρ 0 ). The seismic equation of motion for the displacement u in an isotropic medium is given by: with λ and μ being the Lamé parameters and ρ the density. Eq (1) can be expanded, using small perturbations to the elastic properties (δλ, δμ, δρ) around a background medium with elastic properties λ 0 , μ 0 , ρ 0 (Frederiksen & Revenaugh 2004), to: The wavefield can then be divided into a primary (background) and scattered component (u = u 0 + δu) with the unperturbed wavefield satisfying the unperturbed wave equation Assuming that the scattered wavefield is much weaker than the unperturbed wavefield this gives the first-order Born approximation by discarding higher-order terms: with Q i being a term of the unperturbed wavefield and the perturbed model parameters which is given by eq. (13.22) in Aki & Richards (2002): with u 0 being a solution for the unperturbed medium. Assuming Rayleigh scattering, where the wavelength of the incident wavefield is much larger than the scale of the heterogeneity, the scattering problem reduces to a point scatterer and the full scattered wavefield is approximated by that of an array of point scatterers. Following Wu & Aki (1985), it is possible to derive expressions for the equivalent point source in Rayleigh scattering. These expressions also contain the directivity of the radiation of the scattered wavefield, and are provided as eqs (7) S. Rost et al. (2004). This gives us the ability to compute both the amplitude and radiation pattern of scattering from small-scale heterogeneities in solving the forward problem of the waveform inversion.
We assume the incident P wave to be planar ( Fig. 2) with a known slowness vector, a condition well met for teleseismic records. The scattered wavefield is derived from the seismic observations by computing the 3-component receiver functions relative to the first arriving P wave. The considered input seismic wavefield includes the direct incident P and the free-surface reflections (Pp and Ps), producing forward-and back-scattering in the volume, respectively. The forward and backscattering of the input wavefield produces six possible scattered phases (where • indicates the scattering event along the ray path) at small-scale elastic heterogeneities: P•p, P•s, Pp•p, Pp•s, Ps•p and Ps•s. In the forward modelling, we consider every possible combination of perturbed parameter [P-and S-wave velocity perturbation (δα, δβ) and density perturbation (δρ)], incident wave (forward scattering P and backscattering free surface reflection Pp and Ps) and station location. The inclusion of the free surface backscattered energy as well as the forward scattered direct wave increases the resolution of the study volume and allows us to resolve a 3-D perturbation model, here represented as a regular grid of perturbed cells. We use ray tracing in a 1-D velocity model to determine traveltimes to and from the scattering heterogeneity and to calculate incidence and refraction angles. We use eqs (7)-(10) of Frederiksen & Revenaugh (2004) including a geometrical spreading factor for a layered medium to determine the amplitudes of the scattered energy in an elastic velocity model.
The Born approximation prescribes that single scattered waves propagate in the unperturbed medium and do not interact with heterogeneities again. Therefore, the scattered wavefields from individual heterogeneities are independent. The complete scattered wavefield T can therefore be obtained through the summation of the contributions of individual heterogeneities: with t ij representing the time series representing the scattering contribution of the jth perturbed parameter of the ith scatterer (Frederiksen & Revenaugh 2004). The medium beneath the array is parametrized into a 3-D grid of cells with each cell potentially containing a perturbation of elastic parameters. The perturbation for all cells can be collapsed in an δβ, δρ]. Summing over all contributing elements we obtain the N-element vector d i with the number of displacement samples depending on N = samples × stations × components × events. The dependence of the full scattered wavefield on arbitrary model m is then described as A is an N × M matrix describing the sensitivity of each data point to each model parameter, that is each column of A represents a differential seismogram for a perturbation of a single parameter in a single cell of the perturbed model. Eq. (7) is linear and can therefore be solved using linear inverse theory. To pose this problem as a damped inversion the inverse problem is formulated as a minimization: with I being an M × M identity matrix and λ a weighting factor, representing uniform damping. We use the LSQR method (Paige & Saunders 1982) to solve for the material properties in m. For the inversion of real data it has been found that regularization by smoothing is preferable to damping (Frederiksen & Revenaugh 2004) as it provides results with higher coherence. Using LSQR, the model is smoothed by posing m = Sx with S being a matrix containing a Gaussian smoother. We use B = AS and the minimization in which we solve for x rather than m (Van der Lee & Nolet 1997; Frederiksen & Revenaugh 2004). For all recovery tests and real data inversions, we apply a moving Gaussian smoother with a standard deviation of one model element in the horizontal directions, but we do not smooth in the vertical direction. No smoothing is applied beyond three standard deviations. This choice biases the recovered model towards lateral coherence, making recovered lateral changes more coherent in our study region where we expect strong lateral changes across the NAFZ. The smoothing limits lateral resolution of structures to about 10 km. Vertical resolution of 2 km is defined by the model discretization.
The model space is parametrized as a regular grid with 5 km horizontal grid spacing and 2 km vertical grid spacing with 30 (0-29) cells in horizontal directions and 60 (0-59) in vertical direction. Each cell is treated as a point scatterer with vertical and horizontal locations at depths 2 · j km (j = 0, . . . , 59) and longitude/latitude location of 5 · k (k = 0, . . . , 29), respectively. The maximum grid size is controlled by the maximum memory required to invert the data set (see below). We tested the method with doubled lateral and vertical grid spacing and do not find noticeable differences in the general structure of the solutions except for obvious impacts on the maximum possible resolution of the solutions.

DATA
We use passive seismic data from stations of the Dense Array for Northern Anatolia (DANA) that were installed across the NAFZ in the region of the 1999İzmit and Düzce ruptures (DANA 2012). DANA was deployed between May 2012 and October 2013 and stations were arranged in a quasi-rectangular region of 35 km by 70 km with a nominal station spacing of 7 km ( Fig. 1). Stations were aligned along seven north-south oriented lines (labelled A-F) and 11 east-west lines (labelled 01-11). Seven additional stations were installed in an eastern semi-circle with a radius of about 60 km. Three permanent stations (SPNC, SAUV, GULT) of Bogaziçi University and Kandilli Observatory and Earthquake Research Institute/National Earthquake Monitoring Centre (BU-KOERI/NEMC) located within the DANA network grid were included in the analysis. Stations were equipped mainly with Güralp CMG-6TD and CMG-3T medium broadband and broadband three-component instruments [full information on the network can be found in DANA (2012)]. Data were sampled at 50 Hz.
We use earthquakes within the deployment period with m b >5.5 from the catalogue of the National Earthquake Information Centre (NEIC) and angular distances of 30 • to 90 • . For the permanent stations we add events from 2009 onwards (in total 47 additional events contributing typically a single 3-component seismogram (ZRT) to the data set). Low frequency noise was suppressed by applying a 2-way, 2-pole high-pass filter with cut-off frequency of 0.1 Hz. We calculate 3-component receiver functions (RFs) with a maximum Sketch of the scattered phases included in the full waveform inversion. An incoming planar P-wave wave front interacts with a cell with a parameter disturbance (δα, δβ, δρ) either from the direct wavefront (forward scattered) or the back-scattered wave from the free-surface reflection. The wave type can convert upon scattering from P to S.
frequency of 1.2 Hz using the time domain iterative deconvolution approach by Ligorría & Ammon (1999) deconvolving the Z component from the vertical, radial and transverse components. The calculated receiver functions were visually inspected to select events following these criteria: (1) transverse RFs show lower or comparable amplitudes than radial RFs, (2) the direct P-wave arrival is close to the predicted travel time for a 1-D earth model and (3) no evidence for large amplitude ringing. The pre-processing used to obtain the receiver functions is similar to the method used by Kahraman et al. (2015), but applied to all three components (vertical, radial and transverse) of the traces in our analysis. To remove the first arrival, which does not contain any additional structural information, we mute the first 2.5 s of each trace following the theoretical P-wave arrival. In total, we use 1396 distinct source-receiver pairs (with 3 components) from 176 events in our analysis. The distribution of sources is shown in Fig. 3. Traces were cut and tapered to 100 s and downsampled from the original 50 Hz sampling to 5 Hz . Despite the downsampling, the matrix to invert is very large which limits the achievable resolution and model depth. Typical storage requirements for the matrix inversion using the sparse storage method are ≈338 Gb for a model space dimension (x × y × z) 145 × 145 × 118 km 3 with an element size of 5 × 5 × 2 km 3 and 1396 100 s long traces, sampled at 5 Hz. We are able to invert the full data set without recourse to inverting subsets of data and stacking the resulting images (Frederiksen & Revenaugh 2004;Zhang & Frederiksen 2013) leading to improved image quality of our results.

R E C OV E RY T E S T S
We tested several 1-D background velocity models for data inversion and synthetic data recovery including models by Karahan et al. (2001), Bekler & Gürbüz (2008) and Horasan et al. (2002) and models including constant velocity and linear vertical gradients.
The background models are used for ray tracing to determine traveltimes of the incident and scattered wavefield. While timing of arrivals changes slightly for all realistic velocity models, the overall recovered structure in our tests does not depend significantly on the choice of background model, although depths of interfaces change due to changes in the traveltimes. We chose to use the model by Karahan et al. (2001) for all inversions presented here (Table 1, Fig. 4). This velocity model is derived from seismic experiments in the study area and has been used in previous studies using this data set (Kahraman et al. 2015;Altuncu Poyraz et al. 2015). Fig. 5 shows the result of an inversion of the full (1396 individual source-receiver combinations) noisy synthetic data generated through the perturbation model (Fig. 4b). A subset of the synthetic traces used in this inversion, that is the stations recording event 20123211812 and used in the data inversion, are shown in Figs 4(c) and (d). Synthetic data were generated using ray tracing through the background velocity model with the addition of the scattered wavefield (i.e. the summation of all contributions of the single scatterers in the model). We use a 0.25 s wide Gaussian wavelet as the source time function. Synthetic tests use the source-receiver combinations for each event in the data set, therefore recreating the same resolution as the recorded data set. For comparison we show the recorded and deconvolved data in Fig. 4a) with the first arrival muted. Scattered phases can be seen coherently across the traces. The synthetic traces (Figs 4c and d) show similar structure although clearly are not able to capture the full complexity of the data due to the simplicity of the model (Fig. 4b). Noise is added to the synthetic data through a random number generator (Marsaglia & Bray 1964) using 10 per cent RMS amplitude variation Gaussian noise compared to the direct wave amplitude to produce this noisy synthetic data set (Fig. 4d).
The synthetic model is parametrized with 5 km cell spacing horizontally and 2 km vertically. The model contains a V P = +0.5 km s -1 and V S = +0.3 km s -1 anomaly for a single depth element (2 km) starting at 34 km depth and a V P = -0.5 km s -1 and V S = -0.3 km s -1 anomaly with thickness of 6 km starting at 78 km (Figs 5a and c).   (Karahan et al. 2001).  (d) show the recovery of the input model. Because the inversion uses the Born approximation, which generates signals from localized perturbation, the recovered model will be a bandpass filtered version of the input. We apply a Wiener optimum filter to minimize the effects of the inversion process, mainly to reduce sidelobes to aid interpretation. The optimization filter, as described for example by Gubbins (2004), is obtained by minimising the residual between the desired output g t (Figs 5a and c) and the signal obtained by convolution of the filter f 0 t with the actual output x t (Figs 5b and d) The effect of the inversion and the filter terms acting on a single trace of the synthetic model are shown in Fig. S1. Although the input model in this test does not contain any density (δρ) heterogeneity, Fig. 5(e) shows that the inverted model for the density structure is affected by cross-talk between the different   parameters (more examples given in Fig. S2). However, relative amplitudes ρ in this model are small and the effect is most prominent in areas with velocity anomalies. Tests with models including ρ show that density structure can be resolved. Complete input and output models for this recovery test and further tests are shown in the Figs S2-S8. These tests show that the recovery of velocity and density anomalies is variable within the model volume due to the relative sampling of the model volume by the ray configuration of the data set. Peripheral regions are generally less well resolved than the centre of the volume and areas of reduced resolution show poor recovery of amplitudes in the inverted model (Fig. 5). Within the central zone we do not observe strong depth or amplitude variations of the recovered model, adding confidence to our interpretation. Areas of the model space that are not well resolved are masked in all following figures (and Supporting Information) and the approximate limits of the well-resolved volume are shown in the N-S profiles (dot-dashed vertical lines in Fig. 5f), to which we limit our interpretation. These areas are estimated from our recovery tests as shown in Fig. 5   also away from them, which are not fully removed by the Wiener filtering. These are likely due to the sampling of the volume by the data set as well as the noise added to the synthetic data. Care has been taken when interpreting recorded data inversions to not interpret such artefacts as lithospheric structure.
Changing the depth extent of the inverted model space between 48 and 118 km (in 20 km steps) does not lead to strong changes in the inverted model. A comparison between a 48 and 118 km deep model containing the same structure for V P and V S is shown in Figs S3(a) and (b). This holds even when synthetic traces were generated including structure below the inverted volume (Fig. S3c) showing that heterogeneities underneath the volume are not erroneously mapped into the model volume. In the following section, we show models down to depths of 118 km (60 nodes with 2 km spacing) in a trade-off between achievable resolution, model size and required computer memory. The horizontal smoothing leads to some smearing of energy in horizontal directions. Nonetheless, Fig. 5 shows that terminating discontinuities can be accurately located within 1-2 horizontal elements (i.e. 5-10 km) in the central region of the model space. We also performed recovery tests using other structural models including velocity and density heterogeneities to better understand the performance of the method (for these further recovery tests please see the Supporting Information).

R E S U LT S
The results of the tomographic scattering inversion of the DANA data set are shown in Figs 6 and 7. Fig. 8 presents an interpreted section of the results. Slices in Figs 6 and 7 were extracted from the three-dimensional inversion volume along north-south (Fig. 6) and east-west (Fig. 7) profiles at locations shown in Fig. 1. The locations of the profiles were chosen to be in similar locations to those shown in Fig. 6 of Kahraman et al. (2015) [for an equivalent display to Kahraman et al. (2015) see Figs S9, S11 and S13]. Full solutions are presented in the form of animated GIFs in Figs S15 and S16. A kml file is also provided to display the profiles in their correct geographical location. The model is filtered with the Wiener optimization filter as discussed above. In the following, we report depths at the top of a heterogeneity in the filtered sections which gives the  best agreement of input and recovered model on the recovery tests (e.g. Fig. 5).
Generally, the S-wave images show greater amplitude and are better constrained. The S-wave tomographic images also seem to show more fine scale structure likely related to the shorter wavelength. The density ( ρ) profiles show some of the major structure and are shown in Figs S13 and S15 but suffer from cross-talk as shown in Fig. 5. As the interpretation of the ρ profiles is more difficult and there is no independent constraint on the density structure we do not discuss this parameter further in the text.

Western profile
Profiles for V P (Fig. 6 a) and V S (Fig. 6 b) have been extracted along a longitude of 30.2 • E. Areas with limited resolution as determined from the recovery tests (Fig. 5) have been masked in this profile in transparent grey. The V P profile (Fig. 6a) is dominated by a velocity increase at ∼40 km depth for most of the profile, which we associate with the Moho. The Moho velocity increase bifurcates south of ∼40.4 • with a shallower velocity increase located at ∼32 km depth deepening to 40 km at ∼40.4 • N. The anomaly also seems to fade, that is showing less of a velocity anomaly, south of about 40.3 • N. The point of bifurcation coincides with the surface expression of the southern strand of the NAFZ. A similar Moho signal is observed in the S-wave anomaly at ∼40 km, shallowing to about 38 km within the Armutlu block, which shows a thickening of this interface. The S-wave anomaly does not show the same shallow branch observed in the P waves but shows lower amplitudes south of ∼40.4 • N, that is south of the southern NAFZ strand.
Observed crustal structure includes a weak high V P anomaly at ∼18 km in the Armutlu block with weak, complex V S structure in the Sakarya zone. Complex structure starting at ∼32 km depth (positive and negative anomalies) in V P and V S can be seen in the vicinity of the northern strand (40.7 • N) just overlying the Moho. The V P model shows less structure in the upper crust except a fast anomaly to depths of ∼5 km around the southern branch and a slow (also seen in V S ) overlying fast anomaly between ∼10 km and ∼20 km depth in the Armutlu block.
The high velocity anomaly at 40 km depth is underlain by a strong low V P and V S anomaly at depths of ∼50 km. This anomaly shows lower amplitudes in the Sakarya Zone with the change coinciding with the surface expression of the southern NAFZ strand. We also identify a velocity increase in V P and V S at ∼64 and ∼66 km depth, respectively, around 40.4 • N (southern strand) and ∼74 km in V P beneath the Istanbul zone (with a termination at the northern strand). The V S anomaly shows a low velocity anomaly at ∼68 km depth just north of the northern strand changing to a high velocity anomaly at ∼74 km depth north of 40.9 • N. At greater depths we observe a fast anomaly in V S at ∼78 km depth and a fast anomaly in V P and V S at ∼92 km but showing depth variation in V S . The 78 km anomaly seems to merge with the deeper anomaly in the Istanbul zone.

Eastern profile
The eastern north-south profile at 30.51 • E (Figs 6 c and 6d) shows more structure, especially in the crust, than the western profile despite the close proximity of the two profiles.
We observe a strong, fast V P anomaly at a depth of ∼34 km terminating halfway through the Armutlu block and re-emerging at a depth of ∼42 km just north of the northern strand in the Istanbul zone. In V S we observe a more continuous structure with a high velocity anomaly at ∼36 km depth in the south, stepping to ∼42 km at ∼40.7 • N coinciding with the northern strand. The V P anomaly is 932 S. Rost et al. Figure 7. As Fig. 6 showing west-east oriented slices through the inverted V P (left-hand panels) and V S (right-hand panels) structure. Slices are located in the Istanbul zone (IZ -a,b) at latitude 40.81 • , the Armutlu-Almacık block (AA -c,d) at 40.58 • and Sakarya zone (SZ -e,f) at 40.36 • . Areas with limited resolution as determined from the recovery tests (Fig. 5) masked in grey. Black circles indicate local seismicity as determined by Altuncu Poyraz et al. (2015) within a ±5 km corridor projected onto the profile. The top panel in each subpanel shows SRTM topography along the profile (Farr et al. 2007). weak in the Armutlu block on this profile and seems to terminate at 40.6 • E, while the V S anomaly is more continuous, but also weakens in this region. The amplitude variation of these anomalies cannot be explained by the limitations of the sampling (see Fig. 5).
Especially striking in this profile is the complex V S structure in the Sakarya Zone down to depths of about 30 km manifesting as series of fast and slow anomalies between ∼10 and 32 km (see Fig. S7). The V P structure is similar but weaker than V S . The structure terminates abruptly at the southern strand with little crustal structure in the Armutlu block. The Adapazarı basin (centred at about 40.7 • N) is representing as a low velocity anomaly between 40.6 • N and 40.7 • N to depths of about 6 km (V S ).
Similar to the western profile we identify a slow anomaly in both V P and V S at depths of ∼56 and ∼52 km, respectively. The V S anomaly seems to show more complexity. We identify a weak slow anomaly at ∼76 km depth in the Sakarya zone in V P which appears detectable but much weaker in V S . This anomaly seems to terminate at the southern branch. Fast anomalies are detected at ∼92 km in V P and V S across the profile with shallower fast anomalies for V P and V S at ∼76 km depth beneath the Istanbul zone and the Armutlu block. In V P there is evidence of this interface splitting into a deeper interface deepening to ∼102 km across the southern strand.

Sakarya zone
The west-east profile for V P and V S (Figs 7e and f, respectively) has been extracted along 40.36 • N and is fully located within the Sakarya zone. The Sakarya zone is the southernmost tectonic block in the study region. The inverted scattering tomography model shows a positive anomaly at depths of ∼38 km. In V P this interface shallows to ∼32 km around 30.6 • E. This anomaly seems rather complex and might be discontinuous. We also identify a laterally limited fast anomaly at ∼30 km between 30.4 • E and 30.6 • E. A deeper slow anomaly at about 54 km depth can be seen that shows a slight step down to about 60 km (V P ) at about 30.4 • E and seems complex in V S . The western part of the profile shows a fast anomaly at ∼68 km, with a slow anomaly at ∼78 km in the east. A fast anomaly at ∼98 km depth (94 km in V P ) is identified which terminates at 30.2 • E in V P .

Armutlu block
In contrast to the Sakarya zone, the Armutlu block (Figs 7c and d for V P and V S , respectively) shows more structure down to depths of 40 km. A fast anomaly at ∼40 km terminates around 30.4 • E and appears as shallow as 30-32 km further east in V P . V S also shows the termination but a less pronounced step. The step around 30.6 • E seems to coincide with the profile moving from the Armutlu block to the Almacık mountains. West of ∼30.4 • E, this interface is underlain by a slow anomaly at ∼50 km showing a step to ∼58 km at 30.4 • E in V P . Overall V S seems more complex. We identify several small scale fast and slow anomalies in the crust, the strongest at ∼14 km around 30.4 • E in V S . Slow anomalies shallower than 40 km are indicated between 30.5 • E and 30.9 • E.
A fast anomaly at ∼94 km stretches across most of the profile in V P , with comparable but more complex structure in V S . the V S section also shows more localized structures at depths greater than 80 km.

Istanbul Zone
The Istanbul zone (Figs 7a and b for V P and V S , respectively) shows very little structure down to depths of about 40-42 km where a strong fast anomaly can be detected in V P and V S . This fast anomaly seems to terminate around 30.6 • E for V P but continues across the IZ in V S . A slow anomaly is visible in V S at depths less than 10 km between 30.6 • E and 30.9 • E and a fast anomaly between 30.2 • E and 30.5 • E.
The strong Moho signal is underlain by a slow anomaly around 52 km depth again terminating around 30.6 • E for V P . V P shows a fast anomaly at ∼74 km depth, which like the Moho signal in this block, terminates at about 30.6. • E; the corresponding structure in V S is weaker and discontinuous. A strong fast anomaly at ∼92 km depth can be seen in V P and V S , and again the V S structure is complex.

D I S C U S S I O N
The scattering tomography results show changes in the structure over distances of 10 km in the lithosphere. The smoothing process implemented in the inversion means that more abrupt changes present in the actual structure would also appear smoothed over that distance. These changes can be related to the different structure of the tectonic blocks Bulut et al. 2012) and manifest in, for example the north-south profiles but can be detected in the full data volume (see Supporting Information). Nonetheless, we observe structural changes also in east-west direction similar to those detected earlier (e.g. Bariş et al. 2005;Beyhan & Alkan 2015;Ç ubuk-Sabuncu et al. 2017) where more continuous structure within the tectonic blocks might be expected.
The profile across the Armutlu block follows the strike of the NAFZ east of about 30.7 • E where it leaves the Armutlu block (Fig. 1). The depth slices through the model shown in Figs 6 and 7 show strong changes between the two north-south trending profiles despite their close proximity. Interpreted NS sections are shown in Fig. 8. We have performed recovery tests for the dominant interpreted structure in the tomographic model (Fig. 9). We include complex crustal structure in the Sakarya zone, a Moho step and lithospheric structure in this complex model. We are able to recover the input structure very well in V S and V P with stronger recovered anomaly amplitudes in V S (Figs 9a and b). Some low amplitude spurious signals due to noise and the inversion volume sampling are visible in these models as discussed earlier, but are generally of lower amplitude than the recovered model. These are mainly visible close to the input anomalies (Fig. 9a). The difference in V P and V S recovery indicates that a joined interpretation of V P and V S might be necessary for robust interpretation of the scattering tomographic images. To highlight the most coherent part of the model we stack the depth profiles in longitude and divide these at 30.4 • E to show western and eastern stacks in Fig. 10 for both V P and V S . Fig. 11 shows a schematic of the dominant structure in the tomographic images. Elev.
[m] Comparing the individual slices and the stacked velocity-depth profiles shows that many features are coherent along stretches of the profile, but can change on short scale-lengths in both N-S and E-W directions.

Mohorovičic`discontinuity
In the west, the Mohorovičic`discontinuity (the Moho) is visible in both V P and V S as a dominant fast velocity at depths of ∼40 km with variations in V P in the south and in V S in the Armutlu block. In the east the Moho is shallower in the south [34 km (V P ), 38 km (V S )] but shows a step to greater depths (42 km) between 40.6 • N to 40.7 • N at 30.51 • E.
The deepening of the Moho might indicate the existence of a shear zone at the location of the northern branch which has also been indicated in teleseismic tomography in this region (Papaleo et al. 2017(Papaleo et al. , 2018. The Moho at 30.2 • E is overlain by a slow anomaly in V P between 40.7 • N and 40.9 • N (Fig. 8). In the V S model there is weak evidence for a similar, but weak and intermittent, structure between 40.8 and 41.0 • N .
In the east, the Moho seems much weaker and discontinuous across all three tectonic blocks. The strongest change in Moho depth can be identified around 40.8 • N in the eastern profile where we observe a step from 32 to 40 km coinciding with the surface expression of the northern branch. For the northern branch the discontinuous structure seems to extend into the mantle as discussed later.
While there are strong north-south changes in the profiles in Figs 6-10 we also observe strong east-west changes, for example in the Sakarya zone with a complex Moho structure around 30.5 • E and the weakening of the Moho east of 30.4 • E or the change in the Armutlu block at ∼30.5 • E. The latter might be related to the stepover structure of the NAFZ related to the differential movement of the Armutlu and Almacık blocks and the trend of the suture zones between the tectonic blocks. We also observe a pronounced change of the Moho depth between the Armutlu block and the Almacık mountains at around 30.6 • E (Figs S14 and S15), indicating strong contrasts in crustal structure from the Armutlu block to the SNAFZ shear zone.

Crustal structure
We find evidence for strong crustal structure variation along some of the profiles. The most striking structure is the apparent strong crustal layering south of the southern branch in the Sakarya zone (Figs 6c and d) for both V P and V S (best visible in V S ). The crustal heterogeneity is clearly truncated by the surface location of the southern branch and forward models indicate that it consists of a series of high and low velocity anomalies (e.g. Fig. 9) perhaps related to emplacement of magmatic sills during the Tethys closure (Karabulut et al. 2003). The crust in the Armutlu block on the other hand is relatively homogeneous, adding to the stark difference across the southern NAFZ branch.
Overlying the Moho in the area of the northern strand along the eastern profiles we detect small-scale, complex Moho structure. Modelling indicates that it could be related to a heterogeneity with limited extent approximating a point scatterer perhaps related to the material property changes in the fault zone. The lateral smoothing inherent to our inversions leads to a lack of resolution in this case.
We detect evidence for the Adapazarı basin as low velocity anomalies between 40.6 • N and 40.8 • N in the eastern profiles. Our method does not allow the necessary depth resolution at these depths for conclusion on the depth of the basin. The high velocity Iznik metamorphics (Taylor et al. 2019a) can be detected between 40.4 • N and and 40.6 • N in the western profile.
Areas in the proximity of the surface expressions of the northern and southern strands show more heterogeneous structures than areas further away, perhaps related to increased damage around the fault zone (Ben-Zion & Sammis 2003). We detect a few localized crustal heterogeneities in the Sakarya zone and Armutlu block. There is evidence for a more continuous low velocity anomaly at ∼10 and ∼25 km depth in the Armutlu block and the Sakarya zone that is best visible in the V P models (Fig. 7). The scattering tomography shows less heterogeneity in the Istanbul Zone than in the neighboring tectonic units of the Armutlu block and the Sakarya Zone, which could be related to the reported absence of metamorphism and the lack of major deformation (Okay 1989).

Subcrustal structure
Below the Moho we identify a dominant low velocity layer at depths between ∼50 and ∼60 km in the north-south profiles for V P and V S (Fig. 10). The low velocity layer weakens but remains observable around 40.6 • N and is possibly linked to the surface expression of the northern strand. The weakening is more pronounced in V P than V S . The interface to the anomaly is slightly deeper (52 km) in the stacked eastern profile but also shows changes in the extent of the reduced seismic velocities from traveltime tomography (Papaleo et al. 2017(Papaleo et al. , 2018. The continuity of this structure beneath all tectonic blocks, although with possible depth and structural variations, indicates that it is related to lithospheric structure post-dating the amalgamation of northern Anatolia and the development of the suture zones. It is similar to a signal detected by Kahraman et al. (2015).
The fast anomalies detected at depths greater than 60 km show changes in depth and structure in the vicinity of the surface locations of the NAFZ branches although slightly offset to the north, possibly consistent with shear zones that dip to the North (Kahraman et al. 2015;Papaleo et al. 2017Papaleo et al. , 2018. The tectonic implications of of such a northerly dip remain unclear. There is little evidence for a coherent deep low velocity anomaly in our model that can be interpreted as a lithosphere-asthenosphere boundary (LAB). The lower part of the models seems dominated by high velocity anomalies, although there is weak evidence for a low velocity anomaly in V S between 110 and 120 km depth. Results from previous studies suggesting shallow LAB depths between 80 and 100 km are confirmed in the entire region outside the subduction zones (Kind et al. 2015). Therefore we cannot confirm a detection of the LAB in our models. The LAB might be too gradational to show up as signal in the P-wave coda and to be imaged using our method.

Shear zones
In our scattering tomographic model we see the strongest evidence for the NAFZ shear zone in the abrupt changes of crustal and subcrustal structures. We see crustal structures that terminate on or near both fault strands, most clearly in the changes of the crustal structure transitioning from the Sakarya Zone to the Armutlu block (i.e. across the SNAFZ) which we can trace to Moho depths (e.g .  Figs 6c and d). In general, the AA shows almost no heterogeneity in the crust. At the northern boundary of the AA, coinciding with the NNAFZ we detect energy that seems consistent with the existence of a laterally very limited heterogeneity (Fig. 8) that might be related to the damage zone of the fault zone (Ben-Zion & Sammis 2003).
The Moho step detected in the eastern profiles (e.g. Figs 6c and d) seems to coincide with the surface expression of the NNAFZ and might indicate a localized subvertical shear zone extending deeper than the Moho and into the mantle. Some interfaces in the lithospheric mantle (e.g. Figs 6a and c) also show terminations coinciding with the NNAFZ indicating sub-Moho structure related to the shear zone. Willis et al. (2019) showed that localization of a shear zone in the lower crust can be produced by thermal activation in this tectonic setting if the crust has a rheology comparable to that of dry plagioclase. Furthermore, thermally activated shear zone localization of the upper 10 km or so of the mantle is also possible for a dry peridotite rheology.
Other continental transform faults such as the San Andreas Fault system (SAF), the Alpine Fault (AF) and the Dead Sea Transform (DST) show similar structures (e.g Stern & McBride 1998;Weber et al. 2004;Mohsen et al. 2005;Ford et al. 2014) indicating localized shear throughout the crust. The SAF in southern California (Yan & Clayton 2007) and the DST along the Aravia fault (Mohsen et al. 2005) seem to offset the Moho in close proximity to the surface expression of the fault similar to the eastern profiles across the NAFZ. There is evidence that the SAF also offsets the LAB even in the upper mantle (Ford et al. 2014).
Due to the intra-Pontide suture zone that juxtaposes tectonic blocks of different provenance in the study area, it is difficult to separate the residual signature of a suture zone from the shear zone. The NAFZ seems to exploit a crust weakened by the presence of sutures. Nonetheless, our results provide first evidence that the southern branch might extend throughout the crust (Fig. 6). We also see evidence in the crust indicating small-scale heterogeneity coinciding with the location of the southern and northern strands .

Comparison to receiver function structure
Our results allow a comparison with crust and upper mantle structures resolved in the region using other approaches. Direct comparison with the P-wave receiver function study of (Kahraman et al. 2015) using the same receiver array combines the comparative strength of scattering tomography in imaging lateral and vertical changes with the receiver function sensitivity to vertical discontinuities. Kahraman et al. (2015) noted pronounced variations in crust and upper mantle structure and properties both in north-south and east-west directions in agreement with this and other studies (e.g. Beyhan & Alkan 2015;Ç ubuk-Sabuncu et al. 2017). Kahraman et al. (2015) imaged lateral terminations in key subhorizontal discontinuities beneath the southern and northern fault zones. Here we Downloaded from https://academic.oup.com/gji/article/227/2/922/6318864 by University of Aberdeen user on 30 July 2021 constrain: (1) complex crustal layering in the Sakarya zone towards the east and (2) a mid-crustal feature in the Armutlu block in the west of the study region, which appears to be confined by NAFZ fault branches extending deeper into the crust. Kahraman et al. (2015) find a deepening of the Moho from north to south in the east which is not as clearly seen in this study and they do not show evidence for a step in Moho depth roughly along the surface expression of the northern NAFZ, although there is a high-amplitude receiver function signal beneath the Moho in the Istanbul Zone, which may indicate a subcrustal anomaly. Kahraman et al. (2015) show evidence for structure in the lithospheric mantle in some parts of their profiles that is broadly in agreement with structures resolved in this study. In particular, their anomaly at ∼60 km depth beneath the Istanbul Zone is colocated with a sub-Moho low-velocity zone seen here (Fig. 8). A northwards dipping shear zone, interpreted through the termination of interfaces in the crust and lithospheric mantle in Kahraman et al. (2015), is not evident in the scattering tomography results. Instead, we observe a Moho step beneath the northern NAFZ surface expression and terminations of lithospheric features beneath the NNAFZ could indicate a subvertical extension of the fault zone into the lithospheric mantle. In contrast, the SNAFZ appears to terminate at the Moho in the scattering tomography results, in agreement with Kahraman et al. (2015).

C O N C L U S I O N
We have used data from a dense deployment of seismometers over the actively deforming NAFZ in the region of the 1999İzmit and Düzce ruptures to analyse the scattered seismic wavefield following teleseismic P-wave arrivals. Extending the analysis of the scattered seismic wavefield to a tomographic inversion (Frederiksen & Revenaugh 2004) we detect crustal and mantle heterogeneities that can be linked to the structure and tectonics of the lithosphere beneath the North Anatolian Fault Zone in northwest Turkey. Our highresolution images from the scattering tomography down to depths of 120 km allow unprecedented insight into the lithospheric-scale structure of a major continental strike-slip fault. We show complex structure in crust and lithospheric mantle that can be linked to modern active tectonic processes as well as the structure of the crustal terranes that form the region (Fig. 11).
Our tomographic models show complex crustal structure in the Sakarya zone terminating at the southern branch of the NAFZ and terminations of crustal discontinuities at the northern branch. The terminations of crustal structure are sharp within the resolution of our approach. We observe a step in Moho depth coinciding with the surface location of the northern branch of the NAFZ across most of the study region. Terminations of subhorizontal structures beneath the Moho might indicate that the shear zone extends into the upper mantle to depths of at least 75 km. We detect changes in lithospheric structure perpendicular and parallel to the NAFZ indicating the imprint of complex tectonic history of the region onto the lithospheric structure.
We show that scattering tomography in conjunction with dense recordings of the seismic wavefield is able to provide deeper insight into crustal and mantle structure and the fine-scale structure around fault zones adding to our knowledge of lithospheric structure around fault zones. Strain associated with the NAFZ seems to be localized through the crust and into the mantle. The NAFZ likely exploits weaknesses due to old sutures in this region following the northwards subduction of the Tethys during the amalgamation of Anatolia.    Figure S2. Spurious density structure obtained from inversion of synthetic model described in Fig. 5. The synthetic model includes velocity anomalies, with a homogeneous density structure, so the signal shown here represents cross-talk in the inversion process between velocity and density components of the model. Model is displayed without interpolation between model cells. Figure S3. Comparison of synthetic model with 48 and 118 km depth for V P (a) and V S (b). Slices are at 30.2 • and 30.51 • longitude. Input model contains a 4-km-thick low velocity (−0.5 and −0.3 km s -1 for V P and V S , respectively) layer starting at 16 km depth and a 2-km-thick high velocity anomaly (+0.5 km s -1 ; +0 km s -1 for V P and V S ) starting at 40 km depth, followed by 2 km low velocity anomaly (−0.5 km s -1 ; −0.3 km s -1 for V P and V S ) at 42 km depth. (c) Inversion of model containing high velocity interface at 80 km depth in a 48 km deep model space. Parametrization of all models is 5 km laterally and 2 km vertically. Model is displayed without interpolation between model cells. Please note the difference in colour scale for models in panel (a) and (b). Figure S4. Recovery test for synthetic model containing a single cell (5×5×2 km 3 positive and negative (V P = ±0.3 km s -1 , V S = ±0.18 km s -1 , ρ = ±0.1 kg m -3 ) anomaly. Input model shown on the left with the recovered model on the right. Recovered models are unfiltered and show the effect of the inversion with the Born approximation leading to side-lobes. Model is displayed without interpolation between model cells. Location and velocity change of heterogeneity can be recovered although the full magnitude of the anomaly might not be recovered. Input and retrieved model are shown on the same colour scale. Figure S5. Recovery test for synthetic model containing a block of high velocity (V P = 0.5 km s -1 , V S = 0.3 km s -1 , ρ = 0.2 kg m -3 ) starting at 30.4 • longitude. Top row shows V P input model and bottom ow recovered model. Slices taken from 3-D volume at 40.36 • latitude and 30.51 • longitude. The volumetric velocity anomaly and the boundaries can be resolved although there is some bleeding of the vertical boundary and not the full magnitude of the velocity anomaly can be recovered. Recovered odel contains some sidelobes and inversion artefacts. Input and retrieved model are shown on the same colour scale. Input and retrieved model are shown on the same colour scale. Figure S6. Recovery test for synthetic model containing a continuous 2 km thick negative and positive anomaly at at 36 and 42 km depth (V P = ±0.3 km s -1 , V S = ±0.18 km s -1 , ρ = ±0.1 kg m -3 ). Input model shown on the left with the recovered model on the left. Recovered models are unfiltered and show the effect of the inversion with the Born approximation leading to side-lobes. Model is displayed without interpolation between model cells. Input and retrieved model are shown on the same colour scale. Figure S7. Recovery test for synthetic model containing a continuous V P anomaly that is 2 km thick at at 38 km depth. V P reduction is −0.25 km s -1 with no V S anomaly (V S = 0.0 km s -1 , ρ = 0.0 kg m -3 ). Input model shown on the left with the recovered model on the left. Recovered models are unfiltered and show the effect of the inversion with the Born approximation leading to side-lobes. The inverted V S model shows some anomalies of low magnitude (colour scale provided in recovered V S profile at longitude 30.52 • ) due to cross talk between the model parameters. Model is displayed without interpolation between model cells. Please note the different colour scale between then input and retrieved model. Input and retrieved model are shown on the same colour scale. Figure S8. Recovery test for complex synthetic model containing V P and V S anomalies in the crust, a stepped Moho like structure and a deeper low velocity anomaly. Similar to Fig. 8 but opposite crustal interface velocities in the Sakarya zone complex crust and lithospheric interface. Figure S9. Recovered V S structure of the data along the profile locations from Kahraman et al. (2015). Figure S10. Recovered V S structure along selective north-south gridlines within the well resolved model space. Figure S11. Recovered V P structure of the data along the profile locations from Kahraman et al. (2015). Figure S12. Recovered V P structure along selected north-south gridlines within the well resolved model space. Figure S13. Recovered density (ρ) structure of the data along the profile locations as in Kahraman et al. (2015). Figure S14. Recovered density (ρ) structure along all north-south gridlines. Figure S15. Animated longitudinal slices through the P-wave velocity model (V P ). The full model space is shown. Insert shows profile location, DANA station locations and active faults in the area. Figure S16. Animated latitudinal slices through the P-wave velocity model (V P ). The full model space is shown. Insert shows profile location, DANA station locations and active faults in the area. Figure S17. Animated longitudinal slices through the S-wave velocity model (V S ). The full model space is shown. Insert shows profile location, DANA station locations and active faults in the area. Figure S18. Animated latitudinal slices through the S-wave velocity model (V S ). The full model space is shown. Insert shows profile location, DANA station locations and active faults in the area. Figure S19. KML file of depth profiles for Google Earth.

A C K N O W L E D G E M E N T S
Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the paper.