# Depth-Dependent Anisotropy Along Northwest Segment of the North Anatolian Fault Zone: Evidence for Paleo-Tectonic Structures Contributing to Overall Complexity

Derya Keleş, Tuna Eken, Andrea Licciardi, Andrew Frederiksen, Tuncay Taymaz

► To cite this version:

Derya Keleş, Tuna Eken, Andrea Licciardi, Andrew Frederiksen, Tuncay Taymaz. Depth-Dependent Anisotropy Along Northwest Segment of the North Anatolian Fault Zone: Evidence for Paleo-Tectonic Structures Contributing to Overall Complexity. Journal of Geophysical Research: Solid Earth, 2024, 129 (6), pp.e2023JB028014. 10.1029/2023JB028014. hal-04921809

HAL Id: hal-04921809

https://hal.science/hal-04921809v1

Submitted on 31 Jan 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# JGR Solid Earth

## RESEARCH ARTICLE

10.1029/2023JB028014

**Key Points:**

* Vertically stratified seismic anisotropy exists in the crust and upper mantle along the north-western part of the North Anatolian Fault Zone (NAFZ)

* The NAFZ does not appear to control the deformation beyond ~20–25 km depth

* Unusually high anisotropy below the IZ reflects imprints of the frozen fabric associated with the paleo-tectonic Odessa Shelf

**Supporting Information:**

Supporting Information may be found in the online version of this article.

**Correspondence to:**

D. Keleş,
kelesde@itu.edu.tr

**Citation:**

Keleş, D., Eken, T., Licciardi, A., Frederiksen, A. W., & Taymaz, T. (2024). Depth-dependent anisotropy along northwest segment of the North Anatolian Fault Zone: Evidence for paleo-tectonic structures contributing to overall complexity. *Journal of Geophysical Research: Solid Earth*, 129, e2023JB028014. https://doi.org/10.1029/2023JB028014

Received 6 OCT 2023
Accepted 11 JUN 2024

**Author Contributions:**

**Conceptualization:** Derya Keleş, Tuna Eken, Andrea Licciardi, Andrew W. Frederiksen, Tuncay Taymaz
**Data curation:** Derya Keleş
**Formal analysis:** Derya Keleş, Andrea Licciardi, Andrew W. Frederiksen

**Funding acquisition:** Tuncay Taymaz
**Investigation:** Derya Keleş, Tuna Eken, Andrea Licciardi, Andrew W. Frederiksen

**Methodology:** Derya Keleş, Tuna Eken, Andrea Licciardi, Andrew W. Frederiksen

**Project administration:** Tuncay Taymaz

**Resources:** Tuncay Taymaz

**Software:** Derya Keleş, Tuna Eken, Andrea Licciardi, Andrew W. Frederiksen

© 2024. The Author(s).
This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.

<figure><img src="image_2.png" /></figure>

# Depth-Dependent Anisotropy Along Northwest Segment of the North Anatolian Fault Zone: Evidence for Paleo-Tectonic Structures Contributing to Overall Complexity

Derya Keleş¹ ID, Tuna Eken¹ ID, Andrea Licciardi² ID, Andrew W. Frederiksen³ ID, and Tuncay Taymaz¹ ID

¹Department of Geophysical Engineering, Faculty of Mines, Istanbul Technical University, Istanbul, Türkiye, ²Université Côte d'Azur, CNRS, OCA, IRD, Géoazur, Sophia Antipolis, France, ³University of Manitoba, Earth Sciences, Winnipeg, MB, Canada

**Abstract** The North Anatolian Fault Zone (NAFZ) is a prominent tectonic structure with a significant impact on the observed active deformation in Türkiye. Detailed knowledge of the seismic anisotropy in the crust and mantle along this nascent shear deformation zone provides insights into the kinematics associated with past and present tectonic events. We employed teleseismic earthquakes observed by the Dense Array North Anatolia seismic network to map 3-D variations in crustal and mantle anisotropy in/around the NW segment of the NAFZ. To achieve this, we first performed a harmonic decomposition analysis of P-receiver functions. The results were then used as a priori information to conduct an anisotropic receiver function inversion with the Neighborhood Algorithm that enabled imaging of the actual orientation and geometry of anisotropic structures. SKS splitting measurements are further used to make a comparison between the anisotropic behavior of crustal and mantle structures. Crustal anisotropy parameters estimated in our analyses/models well identify the signature of deformation caused by accumulated strain in the earthquake cycle through the strike of shallow cutting faults in the brittle crust beneath the NAFZ. Diffuse intense anisotropic energy at lower crustal depths was attributed to lattice preferred orientation of crystals or partially molten lenses elongated along the shear direction. Strong harmonic energy variations beneath the northern part of the Istanbul Zone likely reflect imprints of LPO-originated frozen fabric at shallow depths (0–20 km) associated with the palaeotectonic Odessa Shelf, Intra-Pontide Suture Zones or remnants of the Tethys Ocean.

**Plain Language Summary** We analyze teleseismic observations of distant earthquakes to constrain past/present geodynamic events that have shaped the crust and upper mantle structure beneath the NAFZ and its vicinity. By means of established seismological methods, we discovered detailed patterns of crustal and mantle kinematics that give clues about how the Earth's layers deform. In particular, we found that faults in the region respond to earthquakes and are likely to be concentrated in the brittle part of the crust. The unique patterns beneath the north of the NAFZ indicates the signature of a different lithospheric domain, particularly associated with the Odessa Shelf in the lower crust and upper mantle. The crustal deformation caused by the NAFZ shows how strain builds up during earthquake cycles, leading to the creation of faults in the crust. The details of lateral and vertical variations of seismic anisotropy resolved in the study region helps us to better understand the geological processes that shape the Earth's structure in this particular area.

## 1. Introduction

The westward movement of Anatolia is driven by a combination of multiple geodynamic forces (Figure 1). Multiple subduction processes and active fault zones beneath the Anatolian plate produce complicated geological features, such as the NAFZ, a major nascent plate boundary accommodating ongoing continental shear deformation (e.g., Taymaz et al., 2007; Şengör & Yılmaz, 1981). The history of crustal and mantle rocks that have undergone deformation at different depth ranges along the NAFZ is still puzzling within this complex tectonic setting. The surface deformation pattern resulting from the NAFZ is represented by a nearly 20 mm/yr with westward movement which is well correlated with active seismicity (McClusky et al., 2000; Reilinger et al., 2006). Among the best-known earthquakes in the western part of NAFZ are the 17 August 1999 Mw7.4 İzmit and the 12 November 1999 Mw7.2 Düzce 1999 earthquakes (e.g., Çakır et al., 2003).
**Supervision:** Tuna Eken

**Validation:** Derya Keleş, Tuna Eken, Andrea Licciardi, Andrew W. Frederiksen, Tuncay Taymaz

**Visualization:** Derya Keleş

**Writing – original draft:** Derya Keleş, Tuna Eken

**Writing – review & editing:** Derya Keleş, Tuna Eken, Andrea Licciardi, Andrew W. Frederiksen, Tuncay Taymaz

The western portion of the NAFZ splits into northern and southern branches separating the Istanbul Zone (IZ), Armutlu-Almacık Block (AAB) and Sakarya Zone (SZ) from north to south, respectively (Okay & Tüysüz, 1999) (Figure 1). Owing to recent deployments of seismic networks in the west of NAFZ, our knowledge of isotropic and anisotropic characteristics in the crust and mantle between these tectonic units has considerably increased. Early studies focusing on isotropic features in the crustal structure that are mainly based on receiver function analyses and models (e.g., Frederiksen et al., 2015; Jenkins et al., 2020; Kahraman et al., 2015; Saunders et al., 1998; Vanacore et al., 2013), in particular, have indicated a N-S increasing pattern in crustal thickness estimates from the SZ to IZ. Although a strong correlation between the style of active tectonics and crustal complexities throughout shear deformation zones has been documented in and around the NAFZ (Çubuk-Sab-uncu et al., 2017; Eken, Bohnhoff, et al., 2013; Izgi et al., 2020; Koulakov et al., 2010; Licciardi et al., 2018; Turunctur et al., 2023), how significant past tectonic events contribute to the observed seismic anisotropy is not well understood. Hence, detailed investigations/modeling of the structure and distribution of crust and mantle kinematics can provide key constraints on the past/present geodynamic forces and associated tectonic features shaping the region.

A detailed quantification of seismic anisotropy, that is, directional dependence of seismic wave properties, is greatly helpful in distinguishing the effects of past and present tectonic events on the deformation. The alignment of joints or microcracks, presence of a lithospheric fabric developed due to present-day or early (frozen) lattice preferred orientation (LPO) of anisotropic minerals, or highly foliated metamorphic rocks can lead to seismic anisotropy in the crust (Babuska & Cara, 1991; Karato, 1998; Montagner, 1998; Savage, 1999; Sherrington et al., 2004). In addition to the LPO, shape-preferred orientation (SPO) due to the presence of geometrical patterns of impedance contrasts can lead to a preferential fast and slow direction of seismic wave propagation in the crust and mantle (Fouch & Rondenay, 2006). For instance, in the crust, SPO fabrics may develop from fluid-filled cracks that form in response to stress caused by tectonic forces (e.g., Crampin et al., 1984). At depths within the mantle, the presence of melt-filled cracks, lenses, or compositional lamellae (e.g., Kendall, 2000; Walker et al., 2004) may serve as a potential source for the SPO. Two main sources of observed seismic anisotropy in the upper crust are structure-induced anisotropy and stress-induced anisotropy (Boness & Zoback, 2006). Structure-induced anisotropy can occur due to shear-induced deformation in the vicinity of active faults (Eken, Bohnhoff, et al., 2013; Zhang & Schwartz, 1994; Zinke & Zoback, 2000), sedimentary bedding planes (Kern & Wenk, 1990), and preferred mineral orientations (Sayers, 1994). Stress-induced anisotropy can result either from the extensive dilatancy of fluid-filled microcracks or the preferential closure of fractures by the in-situ stress field (Crampin, 1978). The latter results in a parallel alignment of the orientation of fast waves to the maximum horizontal stress (SHmax).

Our knowledge of seismic anisotropy in/around the NAFZ at crustal depths initially relied on local shear wave splitting measurements (e.g., Eken, Bohnhoff, et al., 2013; Hurd & Bohnhoff, 2012; Peng & Ben-Zion, 2004). Lateral variations of the estimated fast polarization azimuths in these studies exhibit a strong spatial correlation between the fast axis of shear waves and the strike of faults. Local S-wave polarization properties along north-western segments of the NAFZ suggests that both stress- and structure induced mechanisms cause seismic anisotropy in the upper 8–10 km of the crust (Eken, Bohnhoff, et al., 2013; Hurd & Bohnhoff, 2012). The joint interpretation of receiver function analyses (e.g., harmonic decomposition, and NA-based receiver functions modeling), for instance, can reveal the actual orientation and geometry of anisotropic structures at various depths of the crust and uppermost mantle (Bianchi & Bokelmann, 2019; Frederiksen & Bostock, 2000; Licciardi & Piana-Agostinetti, 2016). Licciardi et al. (2018) reported an asymmetry in the variation of crustal anisotropy parameters across the central NAFZ using a harmonic decomposition approach that exploits directional dependence of receiver functions. Furthermore, high-resolution 3-D local P-wave arrival time tomography of the crust and uppermost mantle beneath the Anatolia has revealed extension-parallel fast velocity azimuths in the lower crust and uppermost mantle in western Anatolia (Wang et al., 2020). Early core-mantle refracted SKS phase splitting measurements (e.g., Biryol et al., 2010; Confal et al., 2016; Evangelidis et al., 2011; Olive et al., 2014; Paul et al., 2014; Sandvol et al., 2003; Yolsal-Çevikbilen, 2014) and two-layer mantle anisotropy models (e.g., Erman et al., 2022; Merry et al., 2021) showed regionally consistent NE-SW fast polarization orientation (FPO) in Anatolia. These observations were primarily explained by LPO of olivine in the asthenospheric mantle due to the mantle flow induced by slab roll-back and/or extension beneath the Aegean Sea. Anisotropic variations of Pn wave velocities (Al-Lazki et al., 2003, 2004; Confal et al., 2020; Mutlu & Karabulut, 2011) and tomographic models of anisotropic P-waves at crustal and uppermost mantle depths based on the inversion of local earthquake
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Detailed tectonic map of north-western Anatolia. The major tectonic blocks are represented by different colors (see the legend for further explanation). NAFZ: North Anatolian Fault Zone. NNAF: Northern Branch of North Anatolian Fault Zone (NAFZ). SNAF: Southern Branch of NAFZ. IPSZ: Intra-Pontide Suture Zone. IAESZ: Izmir-Ankara-Erzincan Suture Zone. The active fault data presented here is taken from Emre et al. (2018). The tectonic units are modified after Okay and Tüysüz (1999).</figcaption></figure>

travel time residuals (Wang et al., 2020) indicate that P-wave anisotropy in the uppermost mantle is coherent with both SKS splitting directions and GPS velocities in the Aegean Sea, favoring vertically coherent deformation of the Anatolian lithosphere.

Our main aim is to provide significant insights into the geodynamic processes that shape the orientation of anisotropy in the crust and mantle in the study area. To achieve this aim, the present work employs three advanced seismic anisotropy imaging approaches. Firstly, we conducted harmonic decomposition of receiver functions to reveal the localization of anisotropy at crustal and uppermost mantle depths. Secondly, we utilized these findings as a priori information in a Neighborhood Algorithm (NA) based RF inversion in order to resolve the exact geometry of the seismic anisotropy and its relation with active seismicity and tectonic features, especially around the NAFZ. Finally, we analyzed station-averaged single-layer anisotropy via teleseismic core-mantle refracted SKS phase splitting measurements to further understand the mechanical coupling between crust and mantle.

## 2. Geological Setting

The NAFZ is located close to the Intra-Pontide Suture Zone (Intra-Pontide Suture Zones (IPSZ)) and Izmir-Ankara-Erzincan Suture Zone (IAESZ). The evolution of this ~1,500 km-long dominantly dextral strike slip fault (Okay & Tüysüz, 1999) started 13–11 Ma ago. Formation of the fault continued from the Karlıova triple junction in the east toward west to the Sea of Marmara. The ΙΖ located along the Odessa Shelf between the Moesian Platform and Crimea consists of continental-rift fragments derived from the south of the Odessa Shelf in late Cretaceous caused by the tearing of a Hercynian continental sliver (Okay et al., 2015; Figure 2). The rifted ΙΖ and its small fragments moved to the southern side of Black Sea, and were bounded by two large transform faults known as the dextral Black Sea fault in the west (Figure 2) and the sinistral West Crimean Fault (WCF) in the east. The basement of the ΙΖ is characterized by pre-Early Palaeozoic-Late
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Geological map of the study area that is modified after Taylor et al. (2019) and Akbayram, Sorlien, and Okay (2016). Tectonic zones across the branches of North Anatolian Fault Zone are annotated with gray transparent text. The dashed black line represents the passive thrust fault associated with the Intra-Pontide Suture within the Armutlu Peninsula and Almacık Mountains. AB and PB indicate the location of Adapazarı and Pamukova Basin in Istanbul Zone and Armutlu Peninsula, respectively.</figcaption></figure>

Proterozoic metamorphic and metavolcanic rocks (Akbayram et al., 2013; Chen et al., 2002; Yiğitbaş et al., 2004). This basement is unconformably overlain by a sequence of Ordovician-Carboniferous sedimentary deposits (Görür et al., 1997). The extensively distributed outcrops of the Armutlu Peninsula include thrust contacts which are formed by high-grade metamorphic rocks and accretionary complex (Figure 2). These metamorphic basement rocks are covered by younger sedimentary sequences formed in the Upper Cretaceous to Lower Eocene (Akbayram et al., 2013; Okay, 2008; Şengör & Yılmaz, 1981; Figure 2). The Almacık Mountains are comprised of metamorphic ophiolites and volcanic rocks overlain by an Early-Cenozoic sedimentary succession (Gülmez et al., 2013; Robertson & Ustaömer, 2004; Yiğitbaş et al., 2004; Yılmaz et al., 1981). The basement of the SZ is formed of strongly deformed and locally metamorphosed active margin units of the Paleo-Tethys including the Variscan orogen and low-grade metamorphic sequences (e.g., Robertson & Ustaömer, 2004). The sedimentary sequence overlying the whole basement of SZ is relatively young and more diffuse compared to the IZ (Akbayram, Sorlien, & Okay, 2016). Parts of central and western Sakarya Basin are formed by nearly 5 km thick sediments (Figure 2) (Akbayram, Şengör, & Özcan, 2016).

# 3. Data and Method

## 3.1. Data

We used a teleseismic data set recorded by 72 broadband, three-component stations of the Dense Array North Anatolia (DANA) temporary passive seismic imaging array initiated within the framework of the FaultLab project (Figure 3). This array has 72 stations in total; 63 out of them were operated for 18 months (2012–2013) (Figure 3c). We extracted the remaining teleseismic data from six temporary and three permanent seismic stations that are operated by the Kandilli Observatory and Earthquake Research Institute (KOERI). The inter-station spacing is about 7 km except for 7 temporary stations from the KOERI (Figure 3c). In our work, we used two different data sets for receiver functions and core-mantle refracted SKS splitting analyses. For the calculation of receiver functions, we extracted digital waveforms of teleseismic events with Mw ≥ 5.2 with epicentral distances between 30° and 90°. Following a visual inspection procedure, we selected 12,758 high-quality P-wave signals, out of 641 events and 35,869 waveforms, for further P-receiver function computations (Figure 3a). We achieved a
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Station and teleseismic earthquake distribution used in this study. (a, b) Azimuthal equidistant projection centered on the Dense Array North Anatolia (DANA) array displaying the location of teleseismic earthquakes (color coded by event focal depth) employed to generate the RF and SKS phase data sets, respectively. Only events from 30° to 140° of epicentral distance have been used. Size of circles is proportional to event magnitude (M) as per legend. (c) Map showing station distribution (triangles) in the study area. Blue triangles indicate temporary stations within the DANA Network. Green and purple triangles are temporary and permanent stations obtained from the KOERI, respectively.</figcaption></figure>

good azimuthal coverage at most of the stations even though most of the events in our data set approach the array from the first quadrant at all stations.

The selected raw waveform data prior to generating receiver functions in this study were windowed from 20 s before to 180 s after the theoretical P-wave onset calculated based on the 1-D IASP91 radial earth model (Kennett & Engdahl, 1991). Later, we employed a frequency domain deconvolution method using a water-level parameter of 0.1 (Di Bona, 1998; Langston, 1979) to calculate receiver functions. We used band-pass filtering between 0.02 and 1.25 Hz with a Gaussian width of 2.2 to avoid high frequency noise from the raw data.

When a fully isotropic medium exists below the station, the tangential component of receiver functions ideally is not expected to display any signal since the P wave converts only into the SV phase (Sav- age, 1999). The real Earth is rather more complicated and a detailed look into overall azimuthal variation of both radial and tangential receiver functions can enable us to detect anisotropic signatures associated with crustal structure. Anisotropic structures with different orientations can lead to either a two-lobed (π peri- odicity) or a four-lobed (2π periodicity) pattern of directional dependency for the energy variation of P-SH converted phase signals on tangential receiver functions (e.g., Audet, 2010; Bianchi et al., 2010; Lamarque & Julià, 2019; Langston, 1977; Levin & Park, 1997; Licciardi & Piana Agostinetti, 2016; Savage, 1999). In
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Examples from receiver function data set. (a) Locations of all seismic stations (gray triangles) including the selected stations (purple triangles) whose receiver function examples are shown in this figure. The faults and suture zones are represented by black lines and red lines with triangles, respectively. NAFZ: North Anatolian Fault Zone, SNAF: Southern branch of North Anatolian Fault Zone (NAFZ), NNAF: Northern branch of NAFZ, AB: Adapazarı Basin, PB: Pamukova Basin, IPS: Intra-Pontide Suture, IAESZ: Izmir-Ankara-Erzincan Suture Zone. The back-azimuthal variations of all receiver functions generated at stations (b) DD10, (c) DB06 and (d) DE04. Positive and negative polarities are marked by blue and red colors, respectively. The dashed line on radial receiver functions for each station indicates 0 s, corresponding to the direct P phase, while the solid black line shows the onset of the Moho Ps converted phase.</figcaption></figure>

Figure 4, we present azimuthal variations of receiver functions at three selected stations (DD10, DB06 and DE04) that are located in the IZ, the Armutlu-Almacık region near the northern branch of NAFZ and the SZ (near the southern branch of NAFZ), respectively. Figure 4 shows azimuthally varying time shifts and polarities of converted phases on tangential receiver functions obtained at each selected station. The positive
and negative polarities on the traces are represented by blue and red colors, which indicates an increase and/or decrease in absolute velocity contrast between two different layers with increasing depth, respectively. Radial receiver functions provide insight into the Moho signal, which is fairly notable around 4.0–4.75 s (Figures 4b and 4d for the stations DD10, DB06, DE04, respectively). To better understand the origin of periodic signals observed on the tangential receiver functions, a harmonic decomposition analysis on the receiver function data set was performed. Prior to the harmonic decomposition, the radial and tangential receiver functions (RRF and TRF) computed from the three-component seismograms of a single event were migrated from the time domain into the depth domain using the IASP91 radial earth model (Kennett & Engdahl, 1991). The depth conversion process was performed using the mean slowness value at each station.

For the SKS splitting measurements, we extracted a total of 17,732 digital waveform data from 642 teleseismic events with Mw > 5.5 which occurred within the epicentral distance range of 65°–140° (Figure 3b). We calculated theoretical arrival times of SKS phases and other possible contaminating phases (e.g., S, ScS, SS or sS) based on the 1-D IASP91 radial earth model (Kennett & Engdahl, 1991). Then we selected only the events with an analysis window of 40 s (15 s before and 25 s after the SKS phase onset) that is not contaminated by other phases. Finally, we performed a visual inspection on the remaining data set to ensure good quality data (i.e., high signal-to-noise ratio). In this way we obtained 9,755 waveforms from 250 events suitable for further SKS splitting measurements.

## 3.2. Method

The present work focuses on lateral and vertical variation of anisotropic structures in the crust and mantle by implementing three different methodologies. To detect anisotropic signatures in the crust and uppermost mantle and to model the actual orientation of possible anisotropic structures, we perform the harmonic decomposition technique and NA-based inversion on the high-quality receiver functions. In addition, we perform SKS splitting measurements to elucidate the deeper continuation of anisotropic behavior in the mantle lithosphere or asthenosphere.

### 3.2.1. Harmonic Decomposition Analysis

The variations caused by complexities in the real earth (e.g., dipping velocity heterogeneity or anisotropy) affect both the RRF and the TRF and contribute to increase the level of complexity of the observed waveforms, making first-order qualitative interpretations from raw RF data sets difficult. To isolate depth-dependent 3-D effects from the isotropic earth response, we employ the harmonic decomposition method developed by Farra and Vinnik (2000) and established through applications in various tectonic settings (e.g., Audet, 2015; Bianchi et al., 2010; Licciardi et al., 2018; Shiomi & Park, 2008; Tarayoun et al., 2017). This technique allows a quantitative analysis and interpretation of the RRF and TRF simultaneously (Audet, 2015; Bianchi et al., 2010) and offers an advantage over classical RF analyses by particularly separating effects of isotropy and anisotropy with horizontal axis of symmetry on the observed RF signals.

The method approximates the N-S and E-W components of RF in depth (z) domain as a series of cos[k Ø(z)] and sin[k Ø(z)] terms with k = 0, 1, 2. Introducing event back azimuth (Ø) enables the isolation of azimuthally varying features from isotropic ones. The k = 0 term is the azimuthally constant term associated with isotropic features. By stacking the RRF and TRF with a positive shift in back-azimuth (+π/2k) Shiomi and Park (2008) showed that enhancement in signal quality of two-lobed or four-lobed periodicity on k = 1 and k = 2 harmonics, respectively, could be achieved. The two-lobed periodicity (k = 1) is related to dipping velocity contrasts or hexagonal anisotropy with plunging axis of symmetry, whereas the four-lobed periodicity (k = 2) is solely an indication of anisotropic structure with horizontal axis of symmetry (Maupin & Park, 2007; Shiomi & Park, 2008). The following linear system of equations is solved for five coefficients, A (z), B (z), C (z), D (z), E (z) that represent the harmonic terms of the RF data set (Bianchi et al., 2010; Licciardi et al., 2018):
$$ \begin{bmatrix} RF_1^R(z) \\ \vdots \\ RF_n^R(z) \\ RF_1^T(z) \\ \vdots \\ RF_n^T(z) \end{bmatrix} \qquad (1) $$

$$ = \begin{bmatrix} 1 & \cos(\emptyset_1) & \sin(\emptyset_1) & \cos(2\emptyset_1) & \sin(2\emptyset_1) \\ \vdots & \vdots & \vdots & \vdots & \vdots \\ 1 & \cos(\emptyset_n) & \sin(\emptyset_n) & \cos(2\emptyset_n) & \sin(2\emptyset_n) \\ 0 & \cos(\emptyset_1 + \pi_2) & \sin(\emptyset_1 + \pi_2) & \cos(2\emptyset_1 + \pi_4) & \sin(2\emptyset_1 + \pi_4) \\ \vdots & \vdots & \vdots & \vdots & \vdots \\ 0 & \cos(2\emptyset_n + \pi_2) & \sin(2\emptyset_n + \pi_2) & \cos(2\emptyset_n + \pi_4) & \sin(2\emptyset_n + \pi_4) \end{bmatrix} \times \begin{bmatrix} A(z) \\ B(z) \\ C(z) \\ D(z) \\ E(z) \end{bmatrix} $$

where *n* is the total number of events at each station, *RF*<sub>*i*</sub><sup>*R*</sup>(*z*) and *RF*<sub>*i*</sub><sup>*T*</sup>(*z*) are the *i*th radial (R) and tangential (T) receiver functions generated for the event approaching from the back azimuth specified by Ø<sub>*i*</sub>. A(z) represents the harmonic coefficient corresponding to *k* = 0 harmonics. cos[kØ(*z*)] and sin[kØ(*z*)] represent orthogonal N-S and E-W components for *k* = 1 harmonics, respectively. B(z) and D(z) represent the cos[kØ(*z*)] terms while C(z) and E(z) indicate the sin[kØ(*z*)] terms for the *k* = 1 and *k* = 2 harmonics. cos[kØ(*z*)] and sin[kØ(*z*)] terms correspond to the N-S and E-W oriented orthogonal components for *k* = 1. In order to suppress the influence of dipping heterogeneity and anisotropy, Equation 1 can be solved with the implementation of a negative phase shift (−π/2k) (Liu & Park, 2017). Resultant coefficients reflect unmodeled components associated with complex 3-D structures, scattering, or incomplete geometrical coverage and are considered to be an estimate of the error of the harmonic decomposition (e.g., Park & Levin, 2016). Applying the technique at each station involves stacking RF data sets based on Ø and Δ (epicentral distance) with 10° and 20° bin-widths. To minimize bias, events in each bin are selected with zero overlap, preventing the reuse of the same RF in the analysis.

Figures 5–7 display the harmonic decomposition analysis obtained for stations DD10, DB06 and DE04, respectively. The *k* = 0 harmonics (Figures 5a–7a) illustrate the isotropic features at each station. Inferred from the location of a strong Ps converted phase on the signal of *k* = 0 harmonic coefficient, the Moho boundary can be detected considering the IASP91 velocity model mentioned in the depth migration step (e.g., for stations DD10, DB06, and DE04, the Moho depth estimates are approximately 42 km, 37 and 35 km, respectively). In addition, the *k* = 1 and *k* = 2 harmonic coefficients of the N-S and E-W components provide insight into the orientation and strength of deformation (i.e., anisotropy) in the crust and upper mantle. Visual inspection of the amplitudes of the signal, and the positive-negative phase changes on the *k* = 1 and *k* = 2 harmonic terms (see the subfigures in Figures 5–7 for each selected station) suggests dominance of the signal of the *k* = 1 harmonic coefficient over the *k* = 2 coefficient, implying a low significance of the horizontal axis of symmetry, which can be disregarded.

### 3.2.2. Energy Maps and Orientation of Anisotropy

In order to reveal the horizontal and vertical extent of anisotropic structures and/or dipping velocity heterogeneities at crustal depths and to estimate possible anisotropic orientations (fast/slow axes) we map the distribution of energy values and anisotropic axes of symmetry over the study area. We compute the energy values using the following equation (Equation 2) which represents the sum of squares of the depth converted modeled part of the *k* = 1 harmonic coefficient (cos (Ø) and sin (Ø) terms) (after Licciardi & Piana Agostinetti, 2016).

$$ E_{k=1} = \sum_{z=1}^{n} (B(z)^2 + C(z)^2) \qquad (2) $$

where $E_{k=1}$ indicates the total energy of the *k* = 1 harmonic calculated from its N-S and E-W components.
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Harmonic decomposition results for station DD10. (a) *k* = 0 harmonic coefficient. (b) the distribution of the analyzed events in harmonic decomposition analysis with respect to back-azimuthal angle and epicentral distance. Modeled and unmodeled part of *k* = 1 (with φ) and *k* = 2 (with 2φ) harmonic coefficients represented in panels (c, d), respectively. Blue and red colors in the waveforms represent positive and negative polarity, respectively.</figcaption></figure>

Figures 8 and 9 presents the cross-section profiles including depth-migrated *k* = 0 and *k* = 1 harmonic coefficients that are used to gain insights into Moho depth and possible anisotropic layering beneath each profile.

A particle motion analysis is further applied to the depth converted *k* = 1 harmonics to estimate the orientation (i.e., trend) of the symmetry axis of anisotropy (or the dip direction of a dipping velocity contrast). The resulting directions can be mapped on different horizontal slices at selected depths (Bianchi et al., 2010). In brief, orientation information can be extracted from the slope of the particle motion following a simple linear regression procedure. Interpretation of the estimated orientations are conditional on the assumptions about the origin of the observed signals: seismic anisotropy or dipping velocity contrast. In the first case, considering hexagonal anisotropy, the retrieved orientation corresponds to the trend of the plunging symmetry axis, whose plunging direction is determined by the nature of the symmetry axis itself (fast or slow), which is not known a priori (e.g., Sherrington et al., 2004). In the second case, the retrieved angle is the dip direction of a dipping velocity contrast. Even if we are not able to discriminate between anisotropy and dipping velocity contrasts in the upper crust, we note that the presence of strongly dipping (more than 30°) velocity contrasts is unlikely at larger depth and we

<figure><img src="image_8.png" /><figcaption>**Figure 6.** Same as Figure 5 for harmonic decomposition results obtained at station DB06.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Same as Figure 5 for harmonic decomposition results obtained at station DE04.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 8.** Locations of six N-S (A-A', B-B', C-C', D-D', E-E', F-F') and one NW-SE (G-G') oriented profiles (see thin red lines) used for vertical cross-section images presented in Figures 9a–9g. The map shows station distribution (triangles) in the study area. Blue triangles indicate the temporary stations within Dense Array North Anatolia Network. Green and purple stations are temporary and permanent stations obtained from the KOERI, respectively. Black lines are active faults.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>**Figure 9.** (a) Vertical cross-section of harmonic decomposition results taken along profile A-A'. Topography data, the location of tectonic zones and branches of North Anatolian Fault Zone along this profile are shown on the top panel. Depth variation of *k* = 0 and *k* = 1 harmonic coefficients and their interpolated distribution along the profile are given by three cross-sections at the center of the figure. The bottom displays depth variation of interpolated power of components for *k* = 1. (b) Same as panel (a) for harmonic decomposition results obtained for profile B-B'. (c) Same as panel (a) for harmonic decomposition results obtained for profile C-C'. (d) Same as panel (a) for harmonic decomposition results obtained for profile D-D'. (e) Same as panel (a) for harmonic decomposition results obtained for profile E-E'. (f) Same as panel (a) for harmonic decomposition results obtained for profile F-F'. (g) Same as panel (a) for harmonic decomposition results obtained for profile G-G'.</figcaption></figure>

therefore attribute the observed signal to anisotropy only in the lower crust and upper mantle. Following Schulte-Pelkum and Mahan (2014), we plot both the retrieved orientations (gray) together with the strike of the plane perpendicular to it (white) (Figure 10). The latter identifies the strike of foliation or of the plane of preferred cracks or mineral alignment, assuming hexagonal anisotropy with a slow symmetry axis or the strike of the equivalent dipping velocity contrast. Our maps, which highlight the lateral variation of harmonic energy with the orientation of 3D features, provide spatial information about azimuthally varying structure. This information serves as a priori constraint for the implementation of NA-based receiver function inversion.
<figure><img src="image_12.png" /><figcaption>Figure 9. (Continued)</figcaption></figure>

3.2.3. Neighborhood Algorithm Inversion

The results obtained from harmonic decomposition analysis are a good indicator to understand how the anisotropic pattern beneath the receiver side at different depth ranges varies. However, they cannot be used alone to interpret the actual orientation of an anisotropic structure identified by its plunge and trend (fast or slow axis orientation) and its strength. Thus, we performed a receiver function inversion using the NA to examine the 3-D crustal structure (Frederiksen et al., 2003; Sambridge, 1999a, 1999b). Before the inversion we built an initial model that incorporates both isotropic (Delph et al., 2015) and anisotropic features of the crust in the modeled depth extent.

The NA is a kind of global minimization method. This algorithm begins by choosing some number of models randomly from a multidimensional bounded model parameter space described within user-specified ranges. We initially created 100 samples from the parameter space, we further performed 350 iterations with 50 new
<figure><img src="image_13.png" /><figcaption>Figure 9. (Continued)</figcaption></figure>

samples in each iteration. Overall, we obtained 17,600 models to specify the depth dependent anisotropy. Forward calculations of synthetic waveforms were achieved using the RAYSUM software that was originally developed by Frederiksen and Bostock (2000). The complexity of the Earth makes modeling RF data with anisotropic features challenging. For simplicity, the NA inversion approach assumes hexagonal symmetry based on the formulation of Farra et al. (1991). The approach that allows the symmetry axis orientation to vary arbitrarily was originally designed for the mantle where LPO of olivine is widely accepted to be the source of anisotropy (Frederiksen et al., 2003). Therefore, the fast-symmetry-axis case is assigned a positive sign. Nevertheless, an SPO fabric with slow symmetry axis (negative sign of anisotropic strength in the NA parametrization) can develop due to fracturing in the crust as a consequence of the presence of active faults. The fitting process was considered successful when the cross-correlation coefficients between the observed and synthetic receiver functions were maximized. We determined our optimal model parameters below each selected station incorporating our harmonic decomposition results and an isotropic velocity model optimized
<figure><img src="image_14.png" /><figcaption>**Figure 9.** (Continued)</figcaption></figure>

in Delph et al. (2015). The misfit for each inversion iteration was calculated by subtracting the cross-correlation coefficient from one (1-c). Our aim was to select the best-fitting model by minimizing the misfit of all estimated parameters. To assess reliability, we calculated the standard deviation of model parameters from an ensemble of models with misfits less than 1.1 times of the misfit of the best-fit model (Bianchi & Bokelmann, 2019).

### 3.2.4. SKS Splitting Measurements

In addition to quantifying crustal anisotropy, we measured path integrated SKS splitting parameters to investigate average anisotropic pattern and its lateral variation beneath the upper mantle part of the same region. To achieve this, we performed the tangential energy minimization approach in a grid search manner (Long & Silver, 2009;
<figure><img src="image_15.png" /><figcaption>**Figure 10.** Lateral variation of the total energy calculated from the *k* = 1 harmonic coefficient at depths between (a) 0–20 km, (b) 20–35 km, and (c) 35–60 km. The gray and white bars on each station represents the parallel and perpendicular direction to horizontal axis of symmetry calculated by particle motion. Its length is correlated with the energy value. The black lines indicate active faults in each figure. The red line with triangles indicates the Intra-Pontide Suture Zones. The energy values are normalized with respect to highest energy value at each depth. The blue and red colors represent the lowest and highest energy values, respectively. The black dots in panel (a) indicate local seismicity taken from Altuncu Poyraz et al. (2015).</figcaption></figure>

Savage, 1999; Silver & Chan, 1991). We tested several pairs of splitting parameters varying by 1° of increment for the fast polarization orientations (FPOs) and time delays (TD) up to a maximum delay of 5 s at 0.05 s intervals. We visually checked each splitting result to ensure the linearity of the particle motion after correction for optimum splitting parameters. In total, we obtained 990 high-quality splitting examples from 142 events. Uncertainty estimates for splitting measurements were derived using a simplified inverse F-test approach (Silver & Chan, 1991). The absolute range of the error in splitting measurements is rather less reliable since it depends on the assumed number of degrees of freedom of the data, which is difficult to estimate. This test therefore provides a reasonable estimate of relative errors. In the present work one degree of freedom per second of waveform was
<figure><img src="image_16.png" /><figcaption>**Figure 11.** Lateral variation of station-averaged SKS splitting measurements in the study area. Red bars at each station (blue triangle) indicate the fast polarization orientation and are scaled by station averaged splitting time delay. The purple line with triangles represents the Intra-Pontide Suture Zones. Solid black lines indicate the location of active faults.</figcaption></figure>

considered, as previously used in Silver and Chan (1991). After selecting high-quality data and performing transverse energy minimization for each event at each station, we took the average of all individual apparent splitting parameters measured at each station using the Von Mises approach (Cochran et al., 2003) to calculate station average FPOs. This averaging method, similar to the arithmetic mean, considers the 180° symmetry of FPOs.

## 4. Results

Below we will describe our findings from the harmonic analyses, and anisotropic inversion of the RF data set. For this purpose, we plot harmonic decomposition analysis results along seven profiles (A-A', B-B', C-C', D-D', E-E', F-F' and G-G') (see Figure 8 for location of profiles and Figures 9a-9g for details) to reveal the Moho architecture and to identify the signature of possible anisotropic features. Energy maps generated including orientation of seismic anisotropy are used to examine the location of anisotropic contrasts at different depths (Figure 10). Furthermore, we present the SKS splitting measurements to constrain the deformation in the upper mantle (Figures 11 and 12.). The vertical extent of anisotropic structures and their strengths resolved through measurements and models inferred from both RFs and SKS data in this study will allow further interpretations on the origin of seismic anisotropy and deformation processes responsible for it in the crust and mantle part of the study area.

### 4.1. Moho Topography and Anisotropic Crustal Structure

We exploited seismological observations from the DANA array's dense station distribution with an interstation spacing of about 7 km. This allowed us to take depth profiles along lines of seismic stations to discuss 2D variations of isotropic/anisotropic properties along each profile (see Figure 8 for the location of profiles). We built six N-S, and one NW-SE oriented profiles crossing the stations (Figure 9). Nearly 4 km vertical resolution in the results of the *k* = 0 harmonic coefficient gives information about the Moho structure (Piana Agostinetti & Malinverno, 2018). From the *k* = 0 coefficient (isotropic Moho features) estimated for N-S profiles, we detected relatively deep Moho for the stations located in the IZ compared to the stations of SZ (Figures 9a-9f). IZ's crust thickens eastward (38-45 km). Within the Armutlu-Almacık Block, the shallower discontinuity is relatively thin (~34-35 km) compared to the crust beneath the Almacık Mountains, although the Armutlu Peninsula exhibits relatively complicated crustal structure with a notable doubled Moho boundary feature. To the south, the SZ is identified by the lowest crustal thickness (e.g., between 31 and 35 km), which shows a slight eastward increase with respect to other tectonic units (Figures 9a-9d). Overall, the depth sections of RF profiles suggest crustal thickening from SZ to IZ in the entire study area. At the same time, a localized thinning of the crust is notable in and around the Armutlu Peninsula (Figures 9a-9c).

We examine the characteristics of depth-dependent deformation by analyzing *k* = 1 harmonics and energy distributions along specific profiles. In the western part, our findings suggest that the anisotropy is mostly restricted to the crust and a few km into the uppermost mantle (down to depths of ~35–40 km) in the SZ (Figures 9a and 9b). No strong energy on either horizontal component of *k* = 1 harmonics is observed below the Moho in this region. Within the Armutlu Peninsula, anisotropic signals are strong at depths of 15–20 km in the seismogenic zone (Figures 9a and 9b). Below this depth, down to 80 km, less significant anisotropic energy is observed. Profile C-C' shows significant anisotropic energy at uppermost mantle depths (55–60 km) in the SZ (Figure 9c). Toward the west the anisotropy is confined to the crust. In the Armutlu Block and Almacık Mountains region, anisotropic signals can be traced down to the uppermost mantle depths. Depth variation of anisotropy from the SZ to the Armutlu Block indicates a shifted position of anisotropic contrast at shallow depths (~15–20 km) beneath the Armutlu Block (Figure 9). To the north, almost all N-S profiles (except for F-F'; Figures 9a-9e, and 9g) clearly show the anisotropy extending into the uppermost mantle of the IZ with strongest signals along A-A', B-B' and E-E' (Figures 9a, 9b, and 9e).
<figure><img src="image_17.png" /><figcaption>**Figure 12.** Lateral variation interpolated station averaged SKS splitting measurements (a) Interpolated time delays ranging from 1.0 to 2.0 s. (b) Interpolated station averaged fast polarization orientations (FPOs) ranging from 0° to 70°. (c) Interpolated discrepancies between absolute plate motion (APM) and FPOs. APM values are calculated in the No-Net Rotation frame using the NOVEL-1 plate motion model (Argus & Gordon, 1991).</figcaption></figure>
## 4.2. Variations in Depth-Dependent Energy and Anisotropic Orientation

Mapping the distribution of energy values and orientation of possible anisotropic structures at each station for the selected layers reveals the horizontal and vertical extent of anisotropic structures and/or dipping velocity heterogeneities at crustal depths. Figure 10 presents lateral variation of energy values and orientation of 3-D features for the upper crust (0–20 km), lower crust (20–35 km) and uppermost mantle (35–60 km). High energy values (0.8–1.0) near the SNAF suggest significant seismic anisotropy in the upper crust, with particle motion perpendicular to the fault (Figure 10a). Along the NNAF segment and Armutlu Peninsula, relatively high energy values (0.5–0.7) likely indicate deformed zones. At the center of the fault, high energy values tend to be focused into much narrower zones within the upper crust compared to the SNAF (Figure 10a). The symmetry axis of the anisotropy along this part of fault are oriented NNE-SSW (Figure 6). Beneath the Almacık Mountains, a strong and local energy variation is notable. However, insufficient station deployment and large interstation spacing (~40 km) makes the spatial resolution of our estimates rather poor beneath this region compared to tectonic elements in other regions of the study area (Figure 10a). In the northernmost part of the IZ, strong anisotropic variations exist in energy values (~1.0) with nearly N-S orientation represented by gray color in Figure 10a. To the south in the SZ and surrounding area, the symmetry axis of the anisotropy in the upper crust is NNE-SSW oriented with energy values varying between 0.4 and 0.65 along the SNAF. However, the region is generally represented by small values (~0.1–0.3) with varying axis of symmetry (Figure 10a). In the lower crust (25–30 km) we observe more diffuse variation with high energy of anisotropy spread over a larger region than the upper crust (Figure 10b). In particular, the anisotropy within the lower crust of the IZ is characterized with strongly varying energy values (e.g., between 0.65 and 1.0) and with complicated spatial variation of anisotropic orientations. Similar diffuse patterns can be seen in the Armutlu Peninsula and the SZ, with lower energy values (0.35–0.65) (Figure 10b).

In the uppermost mantle (between 35 and 60 km), the largest energy values (~0.9–1.0) are observed beneath the IZ. Near the NNAF, the symmetry axis of the anisotropy exhibits NNE-SSW orientation, while, to its further north at the IZ, anisotropic orientations are rotated to nearly N-S (Figure 10c) (Figure 10c). We excluded the data from stations (DA08, DA10, DD09, DE08, DE09, DF07, DF08, DF09, DF10) at the IZ since they show strong undesired effects of a thick sedimentary basin. In the SNAF and to the south of the Armutlu Peninsula, the observed energy values (e.g., varying between ~0.5 and 0.9) are restricted to a small localized area and do not show any particular coherency with the geometry of active tectonic features (Figure 10c).

## 4.3. Lateral Distribution of Mantle Anisotropy

We present directional variation of SKS splitting measurements at each station (Figures S31–S98 in Supporting Information S1) and a table of station-averaged SKS splitting parameters (see Table S1 in Supporting Information S1) in the supplement. These results indicate coherent FPOs aligned NE-SW with split time delays (dt) generally varying between ~1.1 and 1.8–2 s (Figure 11). In the SZ, the FPOs vary between N30°E and N45°E and time delays are estimated to be between 1.25 and 1.4 s (Figures 11 and 12a). Beneath the Armutlu Peninsula including the Pamukova Basin and near-fault regions in the north and south, splitting measurements generally show sub-parallel anisotropic orientations to the NNAF and SNAF (e.g., FPOs measured from stations DC05, DC06, DB07, DA07, DA08, DB08 vary between 55° and 60° ± 4.0°) (Figure 11). The lateral variation of the FPOs for the entire study area appears to deviate beneath station DC04 located on the boundary between the SZ and Armutlu Peninsula, where the average FPO of ~61° ± 4.07° is aligned almost parallel to the main strike of the SNAF. In the IZ to the north, the anisotropic parameters are characterized by average FPOs oriented along N30°E ± 5° and 1.5 ± 0.22 s delay time (Figure 11). In general, split time delays tend to increase from south to north indicating a notable contrast between stations located at the north and south of the NAFZ (Figure 11a).

Figure 11c compares FPOs with Absolute Plate Motion (APM) directions with respect to a No-Net Rotation reference frame using the NOVEL-1 plate motion model (Argus & Gordon, 1991). Discrepancies between FPO and APM directions increase from the SZ to the IZ. In the SZ, the angular differences are up to 15°, while they reach 30–35° for the Armutlu Almacık Block and a maximum of 50° for the IZ (Figure 11c).
## 4.4. Anisotropic Modeling

### 4.4.1. Istanbul Zone

We applied the NA inversion technique to stations DB08, DC09 and DD10 located in the IZ (Figures S1–S9 in Supporting Information S1). Stations DB08 and DC09 located in the IZ indicate a negative anisotropy for the upper crust with N0°E ± 23.0° and N179°E ± 15.1° oriented symmetry axis (see Figure S1 in Supporting Information S1 for station DB08, Figure S4 in Supporting Information S1 for station DC09, Table 1). Also, the plunge angle of the anisotropic orientation is ~81.5° ± 5.45° beneath these stations. However, the anisotropic model resolved beneath station DD10, located just to the northeast of the IZ, indicates an anisotropic fabric with a positive P-wave anisotropy represented by an axis of symmetry trending N118.0°E ± 16.8° and plunging 34.0° ± 9.2°. The resultant station-averaged fast wave azimuth beneath the station DD10 tends to rotate from N-S to nearly WNW-ESE direction (Figures 13a–13f).

In the lower crust, the trends of symmetry axes for stations DD10, DC09 and DB08 are N1.0°E ± 12.6°, N178.0°E ± 38.2°, and N31.0°E ± 21°, respectively. The plunge estimates for these stations are 21.0° ± 10.0°, 79.0° ± 13.8°, and 47.0° ± 12.8°, respectively. The P- and S- wave anisotropic strengths in the lower crust beneath stations DB08, DC09 and DD10 are given in Table 1. Beneath the station DB08 we observe the largest anisotropic strength, calculated as 18.7 ± 4.7 and 17.7 ± 3.8% for the P- and S- wave, respectively. At greater depths (~35–60 km), stations DB08 and DC09 indicate an anisotropic strength with positive sign (Table 1). NNW-SSE fast axes of symmetry oriented along N175.0°E ± 25.9° and N4.0°E ± 13.9°, are observed for these two stations, respectively. The axis plunges beneath these stations were estimated as 75.0° ± 13.7°, and 38.0° ± 11.8°, respectively. For station DD10, the anisotropic orientation and its plunge angle are N115.0°E ± 22.17° and 45.0° ± 15.8°, respectively (see Figures 13a–13f). Of the three modeled stations, DD10 yields the most robustly recovered model parameters with a misfit value of ~0.60, compared to DC09 (Figure S6 in Supporting Information S1) and DB08 (Figure S3 in Supporting Information S1) in the IZ.

### 4.4.2. Armutlu-Almacık Block

In this block, we performed the NA inversion using four stations (DA04, DA05, DB05, DB06). Our modeling results indicate a N-S anisotropic orientation in about 17 km of the crust corresponding to the seismogenic zone, except for station DB06 located in close proximity to the NAFZ. This station exhibits a N116°E ± 16.5° trending fast wave axis with 71.0° ± 8.9° of plunge angle (Table 1) that is fairly consistent with the main strike of the NAFZ segment (Figures 13h and 13i). Anisotropic strength with negative sign was found for station DB06 (Figure S10 in Supporting Information S1) at this depth (Table 1). For other selected stations in this block, we estimate a negative anisotropy with nearly N-S (varying between N0°E and N4°E) anisotropic trend. In the lower crust, we observed three stations (DA04, DA05, and DB06) with similar orientation (N135°E ± 30.4°, N111°-E ± 19.7°, and N50°E ± 34.5°, respectively) to the NAFZ (Table 1) (see Supplementary section). Inversion results at uppermost mantle depth ranges (down to a depth of 55–60 km) for the stations DA04, DB05, and DB06 show very similar N-S oriented symmetry axes compatible with those resolved for the IZ (Table 1). By contrast, station DA05 gives N114°E ± 22.8° oriented anisotropy with 2.0° ± 8.1° of plunge direction and with positive anisotropy strength, which deviates from the general pattern in this zone.

### 4.4.3. Sakarya Zone

Stations DE04, DC01 and DF04 (Figures S22, S25, and S28 in Supporting Information S1, respectively) located nearby the SNAF show weak S-wave anisotropy with strengths varying between 0.1% and 0.4% (Table 1). The symmetry axes modeled for the upper crustal depths (15–16 km) beneath these stations range from as N160°E to N175°E ± 75.0° with negative anisotropy strengths (see for example station DE04, Figures 13m–13r). In the lower crust, the fast axes of anisotropy vary between ~N68°E and ~N98°E ± 23.0° with ~42–45° ± 11.5° of plunge direction and with negative sign of anisotropy (Table 1). The anisotropic strengths for both P- and S-waves at this depth for all inverted stations vary from 3.9% to 7.8% ± 2.0 with negative sign and 7.6%–11.0% ± 1.8 with negative sign, respectively. At upper mantle depths (~54–55 km), inversion results indicate anisotropic structure identified with positive strength and with the average orientation along N25°E ± ~25° direction. This orientation is consistent across the study area, with a slight rotation from N-S to NE-SW toward the south. For stations DE04, DC01 and DF04 our inverse modeling resulted in the plunge directions of anisotropic orientation with 44.0° ± 9.98°, 49.0° ± 12.6°, and 81.0° ± 10.65°, respectively. Station DE04 yielded the
<table><caption>**Table 1**
Inverted Parameters Obtained From NA Inversion at Selected Stations</caption><thead><tr><th rowspan="2">Tectonic zone</th><th rowspan="2">Stations</th><th colspan="4">Upper crust</th><th colspan="4">Lower crust</th><th colspan="4">Uppermost mantle</th></tr><tr><th>H (km)</th><th>T (°)</th><th>P (°)</th><th>KP (%)</th><th>KS (%)</th><th>H (km)</th><th>T (°)</th><th>P (°)</th><th>KP (%)</th><th>KS (%)</th><th>H (km)</th><th>T (°)</th><th>P (°)</th><th>KP (%)</th><th>KS (%)</th></tr></thead><tbody><tr><td rowspan="2">Istanbul Zone</td><td>DB08</td><td>16.1 ± 0.4</td><td>0.0 ± 23</td><td>79.0 ± 5.5</td><td>-14.9 ± 1.23</td><td>-0.3 ± 1.51</td><td>15.6 ± 0.53</td><td>31.0 ± 21</td><td>21.0 ± 10</td><td>18.7 ± 4.7</td><td>17.7 ± 3.8</td><td>24.1 ± 0.4</td><td>175.0 ± 25.9</td><td>75.0 ± 13.7</td><td>9.9 ± 1.36</td><td>3.7 ± 1.89</td></tr><tr><td>DC09</td><td>18.4 ± 0.22</td><td>176.0 ± 15.1</td><td>83.0 ± 5.4</td><td>-13.6 ± 1.5</td><td>0.0 ± 2.04</td><td>14.7 ± 0.4</td><td>178.0 ± 38.2</td><td>79.0 ± 13.8</td><td>14.4 ± 3.35</td><td>16.3 ± 3.37</td><td>25.1 ± 0.11</td><td>4.0 ± 13.9</td><td>38.0 ± 11.8</td><td>8.9 ± 1.7</td><td>10.9 ± 2.8</td></tr><tr><td rowspan="2">Armutlu Almak Block</td><td>DD10</td><td>17.9 ± 0.25</td><td>118.0 ± 16.8</td><td>34.0 ± 9.2</td><td>20.0 ± 1.9</td><td>0.1 ± 3.8</td><td>14.0 ± 0.7</td><td>1.0 ± 12.6</td><td>47.0 ± 12.8</td><td>13.4 ± 2.65</td><td>10.7 ± 3.33</td><td>25.4 ± 0.2</td><td>115.0 ± 22.17</td><td>45.0 ± 15.8</td><td>7.3 ± 2.57</td><td>10.0 ± 2.88</td></tr><tr><td>DA04</td><td>18.5 ± 0.3</td><td>0.0 ± 33</td><td>89.0 ± 11</td><td>-12.8 ± 3.7</td><td>-14.9 ± 2.15</td><td>16.5 ± 0.53</td><td>135.0 ± 30.4</td><td>86.0 ± 23.6</td><td>-9.3 ± 2.97</td><td>-7.6 ± 2.24</td><td>25.5 ± 0.1</td><td>6.0 ± 25.2</td><td>21.0 ± 11.7</td><td>0.3 ± 2.85</td><td>17.1 ± 3.14</td></tr><tr><td rowspan="2">Sakarya Zone</td><td>DA05</td><td>18.5 ± 0.15</td><td>4.0 ± 10.7</td><td>80.0 ± 6.22</td><td>-10.3 ± 2.7</td><td>-15.0 ± 1.43</td><td>16.8 ± 0.37</td><td>111.0 ± 19.7</td><td>29.0 ± 10.4</td><td>10.4 ± 3.09</td><td>15.2 ± 3.9</td><td>25.4 ± 0.1</td><td>114.0 ± 22.8</td><td>2.0 ± 8.1</td><td>0.1 ± 4.5</td><td>18.0 ± 6.04</td></tr><tr><td>DB05</td><td>15.8 ± 0.3</td><td>0.0 ± 18</td><td>84.0 ± 5.5</td><td>-15.8 ± 2.45</td><td>-20.0 ± 2.58</td><td>16.6 ± 0.4</td><td>160.0 ± 42</td><td>64.0 ± 15.4</td><td>-0.2 ± 2.52</td><td>-13.7 ± 3.8</td><td>25.8 ± 0.3</td><td>174.0 ± 44.7</td><td>17.0 ± 16.4</td><td>0.6 ± 1.7</td><td>1.5 ± 2.3</td></tr><tr><td rowspan="2">DF04</td><td>DB06</td><td>16.5 ± 0.5</td><td>116.0 ± 16.5</td><td>71.0 ± 8.9</td><td>-14.1 ± 2.17</td><td>-19.7 ± 4</td><td>16.0 ± 0.3</td><td>50.0 ± 34.5</td><td>81.0 ± 6.37</td><td>-5.0 ± 2.05</td><td>-12.5 ± 2.93</td><td>25.4 ± 0.3</td><td>1.0 ± 24</td><td>41.0 ± 14.8</td><td>11.9 ± 2.9</td><td>4.2 ± 2.65</td></tr><tr><td>DC01</td><td>17.5 ± 0.35</td><td>179 ± 18.35</td><td>74.0 ± 9.3</td><td>-5.0 ± 0.6</td><td>-0.4 ± 1.15</td><td>15.1 ± 0.3</td><td>68.0 ± 20.3</td><td>42.0 ± 10.7</td><td>-7.5 ± 1.2</td><td>-7.6 ± 1.5</td><td>23.9 ± 0.55</td><td>27.0 ± 18.4</td><td>49.0 ± 12.6</td><td>7.6 ± 1.74</td><td>7.3 ± 2.13</td></tr><tr><td rowspan="2">DF04</td><td>DE04</td><td>17.0 ± 0.26</td><td>161.0 ± 15.02</td><td>81.0 ± 5.53</td><td>-19.8 ± 1.71</td><td>-0.1 ± 1.54</td><td>15.1 ± 0.45</td><td>81.0 ± 22.45</td><td>54.0 ± 12.1</td><td>-3.9 ± 2.39</td><td>-11.1 ± 2.88</td><td>25.1 ± 0.27</td><td>22.0 ± 17.2</td><td>44.0 ± 9.98</td><td>11.9 ± 1.4</td><td>3.2 ± 2.5</td></tr><tr><td>DF04</td><td>14.0 ± 0.44</td><td>178 ± 47.6</td><td>78.0 ± 9.1</td><td>-4.9 ± 0.7</td><td>-0.1 ± 0.75</td><td>15.0 ± 0.2</td><td>97.0 ± 23.1</td><td>30.0 ± 11.2</td><td>-7.8 ± 1.35</td><td>-8.0 ± 1.56</td><td>23.7 ± 0.6</td><td>19.0 ± 38</td><td>81.0 ± 10.65</td><td>11.6 ± 2.73</td><td>8.6 ± 2.84</td></tr></tbody></table>

Note. All results from the stations are clustered according to tectonic zones. Each column represented by blue color shows the upper crust, lower crust and uppermost mantle from left to right, respectively. These columns include inverted parameters which are thickness of layers, horizontal axis of symmetry representing the angle from the north, dip angle of the horizontal axis of symmetry, strength of P wave anisotropy and S wave anisotropy.
<figure><img src="image_18.png" /><figcaption>**Figure 13.** Results of the NA inversion of receiver functions for stations DD10 (top row), DB06 (middle row), DE04 (bottom row). Each graph except (f), (l), and (r) represents variations of model parameters within the upper crust (light blue-filled bar), lower crust (purple-filled bar), and upper mantle (turquoise-filled bar). The best model parameter in each subfigure is marked by red line on the box chart. (f), (l) and (r) exhibit overall variation of misfit values with varying model numbers for stations DD10, DB06, and DE04, respectively. The star in (f), (l) and (r) represents the best model misfit achieved; the misfit variation is color-coded by gray scale ranging between 0.6 and 1.0. See the supplementary materials for the details.</figcaption></figure>

anisotropic strength estimate for the uppermost mantle with the largest positive sign for both P- and S- waves (11.9% ± 1.4 and 3.2% ± 2.05, respectively). Station DE04 has the highest positive anisotropy strength for both P- and S-waves (11.9% ± 1.4 and 3.2% ± 2.05, respectively). This station is considered as the one with the most reliable model parameters given the lowest misfit value of 0.69 (see Figures 13m–13r).

# 5. Discussion

## 5.1. Seismic Anisotropy

In this section, we present a joint interpretation of the depth-varied anisotropy, its origin, and how it relates to the past and present tectonic and deformation history of the region.

### 5.1.1. Upper Crust

The major question that arises when explaining the observed anisotropic signal at upper crustal depth ranges is whether the source of anisotropy is structure- or stress-controlled (Boness & Zoback, 2006). In case of structure-controlled anisotropy, the alignment of macroscopic features associated with existing faults (e.g., Tadokoro et al., 2002), sedimentary bedding planes (Kern & Wenk, 1990), and preferred mineral alignments (Sayers, 1994), can result in fast shear waves parallel to the strike of these structural fabrics. The extensive dilatancy of fluid-filled microcracks (Crankin, 1987) or preferential closure of fractures by the in-situ stress field (Boness & Zoback, 2006) are often considered to be the main cause for stress-induced anisotropy. In this case, vertically propagating shear waves move fastest along the maximum horizontal compressive stress (SHmax).

The seismogenic areas and the geometry of relevant fault structures have been earlier well characterized by precise relocations of micro-seismic activity (e.g., Altuncu Poyraz et al., 2015) within the brittle upper crust beneath our study area (Figure 10a; Hurd & Bohnhoff, 2012; Peng & Ben-Zion, 2004, 2005). Radial anisotropy variation within the crust following a full-waveform inversion of local/regional earthquakes by Çubuk-Sabuncu et al. (2017) revealed strong negative anisotropy close to the weakness zone related to the NAFZ and around the Armutlu Peninsula in the upper crust (~20 km). Previous shear wave splitting analyses conducted along various segments of the NAFZ show that the deformation is usually localized around the fault zone within the upper
10 km (Eken, Bohnhoff, et al., 2013; Hurd & Bohnhoff, 2012; Peng & Ben-Zion, 2004, 2005). However local shear wave-derived anisotropic parameters suffer from poor vertical resolution. Anisotropic inversion of ambient noise-derived dispersion curves (Taylor et al., 2019) suggested a consistent upper crustal anisotropy parameter in this area.

Our work that exploits azimuthal variation of converted waves has for the first time revealed 3-D variations of strong deformation zones identified by high energy values. These zones appear to be restricted to a narrow portion of the fault zone. Further validation of our harmonic energy estimates comes from Rost et al. (2021) that suggested a more homogeneous crustal structure in the Armutlu Peninsula, which is confirmed by our low harmonic energy distribution resolved for upper crustal depths for this region. The NA inversion in this study, yielded well-constrained models with anisotropic trend and strike of fast plane estimates consistent with geomorphological observations and existing knowledge of neo-tectonic features. We found the highest strength of anisotropy with negative sign (slow-axis perpendicular to the strike of the fast plane) following our NA-based anisotropic inversion at stations near the SNAF and NNAF branches. This implies a clear role for the SPO fabric, which can develop due to fracturing in the crust as a consequence of the presence of active faults. The SPO in this region is thought to have developed as a consequence of the presence of active NNAF and SNAF branches, which may have led to the formation of impedance contrasts that facilitate anisotropy in the upper crust. This is achieved by the development of oriented fluid-filled cracks and microcracks in response to tectonic forces (Fouch & Rondenay, 2006; Savage, 1999). At locations away from the SNAF and NNAF within the Sakarya and IZs, our models tend to correlate well with SHmax in the region, indicating a dominant role of the stress regime for the observed anisotropic signal on converted phases. Beneath the northern part of the IZ, another zone of strong harmonic energy variation becomes significant farther away from the NNAF at upper crustal depths (e.g., 0–20 km) (Figure 10b). However, structure-induced seismic anisotropy does not seem to be the primary source. This observation mainly depends on the NA-based model outcomes suggesting a strong anisotropy (similar to harmonic decomposition results) with positive sign below station DD10 (Figures S7–S8 in Supporting Information S1 and Table 1). This is most likely the signature of a frozen fabric of LPO origin at these shallow depths, presumably associated with palaeotectonic features, that is, the Odessa Shelf, IPSZ or remnants of the Tethys Ocean. In fact, geological records (e.g., Akdoğan et al., 2021; Okay et al., 1994) indicate the presence of a different lithospheric domain in this part of the study area, based on strong similarities between stratigraphic data from the IZ and the Moesian Platform. Further evidence for the morphological fit of the IZ to Odessa Shelf between Moesia and Crimea supports the idea that the IZ represents the southern margin of the Odessa Shelf north of the Black Sea. Akdoğan et al. (2021) claimed that the IZ is contiguous with the Moesian Platform prior to the opening of the Black Sea, and later it moved southward to its current position during the Late Cretaceous opening of the Western Black Sea Basin. Recently reported thick crust estimates based on Common Conversion Point (CCP) stacking of receiver functions in Jenkins et al. (2020; Figure 14) or from the inversion of teleseismic P-coda autocorrelation functions (Eken et al., 2021) further implied changes in lithospheric properties in this region. The actual orientation of anisotropic structure resolved for the station DD10 after the NA inversion is compatible with the NNW-SSE trend of Odessa Shelf further supports the imprint of relic petro-fabric related to this paleo-tectonic feature (Figures 10 and 14).

### 5.1.2. Lower Crust

Our harmonic decomposition analysis reveals a diffuse anisotropic energy pattern at lower crustal depths (20–35 km). In a recent teleseismic tomography experiment conducted for the same region, Papaleo et al. (2017) attributed variation of low velocity zones imaged within the crust and upper mantle beneath the NAFZ to a narrow shear zone widening in the upper mantle, where its width reached ~50 km. Based on 3-D finite element modeling of the earthquake cycle in an elastic crust above Maxwell viscoelastic layers, Yamasaki et al. (2014) observed a spatial variation of crustal viscosity consistent with pre-seismic strain localization within about 25 km of the fault and a rapid post-seismic transient evidenced from temporal behavior of geodetic data. These authors showed that relatively sharp weak-zone rheology bounded by abrupt material property variations could occur due to lithological contrast, grain size reduction, fabric development, or water content, in addition to shear heating. Houseman (2018) suggested that this sub-vertical zone of shear localization remained aseismic, likely moving steadily at slow rates below the seismogenic zone in the upper crust and down to Moho depths.

Additional proof of this weak zone at lower crustal depths is provided by 3-D variations of geoelectrical properties found from magnetotelluric imaging along the western part of the NAFZ (e.g., Kaya et al., 2009; Tank
<figure><img src="image_19.png" /></figure>

**Figure 14.** Cartoon illustration of the sketch model inferred from our receiver function modeling and SKS splitting measurements for the study area. The topography is extracted from GeoMapApp. The color bar represents variation in the topography. The black solid lines show active major faults in the region. The boundary between the lower crust and uppermost mantle is marked with three significant tectonic areas known as IZ, AAB, and SZ. These zones are also referred to as the Istanbul Zone, Armutlu-Almacık Block, and Sakarya Zone, respectively. The double-reversed directed arrows indicate the right lateral motion of the North Anatolian Fault Zone. On the right side the angular representation corresponds to the horizontal axis of symmetry and its plunge angle is given as text.

et al., 2005) (Figure 10b). These studies emphasized the possible role of partial melting in producing an observed conductive body at 30–50 km depth between two strands of NAFZ. The diffuse anisotropic energy pattern evident from our harmonic decomposition analysis at lower crustal depths (20–35 km) can be associated with this weak zone where the LPO of minerals during ductile flow and/or partially molten lenses are elongated along the shearing direction (e.g., Brownlee et al., 2017; Hirn et al., 1995; Weiss et al., 1999). The NA-based inverse modeling beneath station DA05 results in a positive anisotropic strength for these fabrics that may be explained by a possible LPO origin (Table 1 and Figures S16–S17 in Supporting Information S1) (Frederiksen et al., 2003).

### 5.1.3. Upper Mantle

The energy map for the uppermost mantle shows that relatively strong anisotropy is mostly localized in the IZ to the north of the NNAF (Figure 10c). In this region, the harmonic decomposition and NA modeling indicate a clear N-S horizontal axis of symmetry (Figure 14). This slightly deviates from our NE-SW directed SKS-derived fast polarization azimuths that are coherent with the regional pattern (e.g., Biryol et al., 2010; Confal et al., 2016; Erman et al., 2022; Evangelidis, 2017; Evangelidis et al., 2011; Hatzfeld et al., 2001; Kaviris et al., 2018; Merry
et al., 2021; Olive et al., 2014; Paul et al., 2014; Sandvol et al., 2003; Schmid et al., 2004; Yolsal-Çevikbi-
len, 2014) (see the bottom layer in Figure 14).

According to Confal et al. (2020), the results of isotropic P-wave tomography indicate that a relatively high-
velocity and thick lithosphere exists in the IZ. Jenkins et al. (2020) suggest that the thick crust imaged from
CCP stacked RFs in the IZ may be related to the motion of the Odessa Shelf accommodated by the combined
activity of the West Black Sea Fault (WBSF) and WCF which are one of the most important tectonic features
during the opening of the Black Sea (Figure 1). The strong anisotropy observed in the IZ at all depths, far from
active fault zones, may be the result of this paleo-tectonic activity. The N-S axis of symmetry resolved in our NA-
based inverse modeling supports the presence of an anisotropic fabric beneath the IZ and is well-correlated with
the main strike of the inactive WBSF and WCF (Kalvoda et al., 2008; Yalçın & Yılmaz, 2010) (see Black Sea
Fault location and its orientation in Figure 1). Besides, a comparison of the lateral variation in our SKS-derived
FPOs with APM directions exhibits an apparent discrepancy within an area spanning from the IZ to the south over
the Armutlu Peninsula (Figure 13c). This implies a possible remnant frozen-in anisotropic fabric with paleo-
subducting slab origin. At the south, beneath the SZ, however, very small differences between APM and
FPOs provide clear evidence for sub-lithospheric LPO considered to be the main source for the observed
anisotropic signal on core-mantle refracted SKS phases (Confal et al., 2018). This indicates a termination of this
paleo-structure effect on the observed anisotropic signals in this part of the study region (Figure 12c). The LPO
type of anisotropy beneath the SZ develops under a regional-scale asthenospheric mantle flow regime that
occurred due to a combination of various dynamic effects (e.g., the Nubian-Eurasian plate convergence with slab
rollback in the Aegean Sea, a tear in the African slab, and detachment of the Arabian plate) (Confal et al., 2018;
Faccenna et al., 2013). Notable discrepancies between APM and FPOs imply that the contribution of this local
paleo-structural effect of the uppermost mantle alone to the observed seismic anisotropy is likely to be isolated
from the rest of the mantle structure beneath the IZ and terminated beneath the SZ.

A comparison of the crustal anisotropy pattern with that in the mantle can provide useful insights into possible
crust-mantle interactions that are controlled by the temperature and pressure conditions of the mantle or tectonic
movements affecting the flow regime mantle. Kinematics in the crust and mantle revealed from harmonic
decomposition, NA-based modeling, and SKS splitting measurements imply vertically incoherent distribution of
deformation due to differences in the origin of deformation in different layers (Figure 14). For instance, to the
north, the presence of frozen fabrics associated with the Odessa Shelf beneath the IZ in both the upper and lower
crust results in a different anisotropic orientation from the shearing direction in the upper mantle. In the Almacık
Zone, between the branches of NAFZ, the upper and lower crust suggest a complex anisotropic pattern that is not
coherent with horizontal basal drag characterized by regional NE-SW oriented fast shear wave directions
(Figure 14). Previous RFs analyses (e.g., Licciardi et al., 2018; Vinnik et al., 2016) in the central part of the NAFZ
attributed the anomalous character of the observed anisotropy on converted waves at around depths of 35–60 and
60–90 km to peridotite-type rocks containing melt with fast orientations normal to the direction of shear in the
mantle. Under low-temperature and high-stress conditions, Jung and Karato (2001) and Holtzman et al. (2003)
demonstrated that fast orientations can develop normal to the direction of flow in water-rich olivine and molten
peridotite, respectively. Karato et al. (2008) observed that the strain divided between the melt-enriched bands and
melt-depleted lenses may influence the rock fabric and upper mantle anisotropy. This may imply that the
deformation was not transferred coherently through the lower crust and sub-Moho structures with FPO estimates
oblique to the main fault branches in our study region. A similar decoupling observed in our study area may
suggest that crustal deformation does not penetrate deeply into the mantle lithosphere. Yet, an answer to the issue
of whether a localized channel flow driven by pressure gradients in response to an upper mantle vertical loading
(Castellanos et al., 2020) controls deformation in the mid-lower crust requires further geodynamic modeling.

**5.2. Crustal Kinematics Associated With the NAFZ**

Anisotropic constraints for the upper crustal part of the central NAFZ (Licciardi et al., 2018) do not show a
complete match with those observed in the western NAFZ (this study). Licciardi et al. (2018) reported strong
anisotropy in the central NAFZ that is significantly concentrated in the southern part. However, to the west, we
observe that seismic anisotropy related to the NAFZ exists across the branches of the NAFZ. In addition to this,
our results highlight a considerable amount of anisotropy in the northern part of the IZ (Figures 9 and 14).
Approaching from the central to western NAFZ, the NAFZ splits into two branches. Variations along the
westward continuation of neo-tectonic structures and available paleo-tectonic features might have led to a relatively complex deformation pattern in the upper crust.

Another important debate is related to the question of how the NAFZ continues in the crust and lithosphere. To shed light on this issue, we investigate the depth-variation in anisotropic parameters that may result from a change in the dominant mineralogy at different depth layers including upper, lower crust and mantle (Bernard & Behr, 2017; Frothingham et al., 2023). These microscale variations can strongly control the regional lithospheric strength and depth distribution of strain. In the present study, however, we do not specifically aim to model possible changes in the mineralogical properties of the rock, which would be in the scope of a future study. Rtaher, we primarily examine depth-varying changes in macro-scale seismic anisotropy properties to understand the importance of lithospheric strain along the continental transform faults. Our motivation was to infer the depth extent of the NAFZ and the associated accumulated strain that still remains enigmatic. Early examples of two-layer SKS splitting measurements (Özalaybey & Savage, 1995; Savage & Silver, 1993) showed that lithospheric strain due to the San Andreas Fault (SAF) controls the seismic anisotropy parameters of the upper layer, while E-W oriented flow in the asthenosphere is consistent with the anisotropy of the lower layer. Zhang and Karato (1995) suggested that the amount of lithospheric strain would decrease with depth along the SAF because the asthenosphere is less coupled to the upper layers. In another example from New Zealand, the observed coupling between lithospheric and mantle anisotropic parameters implied a large and consistent fault-parallel motion across the width of the South Island where a long-term strain is acting (Audoine et al., 1997; Klosko et al., 1997). Integrating receiver function-derived anisotropic parameters with those from SKS measurements can provide a simple guide to interpreting the change in deformation with depth.

Along the NAFZ, the apparent mismatch between anisotropic orientations at varying depths can be interpreted as the lack of vertical extent of the fault into the ductile lower crust and further into the lithosphere. Seismic images of the lithosphere beneath the NAFZ from previous tomographic inversions revealed varying results in different studies. Biryol et al. (2011) and Fichtner et al. (2013a, 2013b) based on teleseismic P-wave arrival time and full waveform inversions suggested a strong velocity contrast down to 100–150 km depth across the NAFZ. A Rayleigh wave tomography study by Salaün et al. (2012), however, did not find significant velocity changes within the mantle lithosphere. Similarly, fast P-wave orientations in Wang et al. (2020) resolved for the upper-most mantle beneath the western NAFZ revealed variation of NE-SW FPOs, uniformly oblique to the strike of the major fault segments (Figure 11). Vertical continuation of the shear deformation zones is hinted at by, for instance, lateral variations of simple SKS splitting parameters. In northeast Tibet or in Pamir, the anisotropic orientations estimated in close proximity to the North and South Kunlun Faults or Talas-Ferghana and Karakorum faults, respectively, were found to be parallel to the strike of these faults (Eken, Tilmann, et al., 2013; Kufner et al., 2018; León Soto et al., 2012). This suggested the hypothesis that these faults extend deeper into the lithosphere. A recent SKS splitting study (Erman et al., 2022) that involved two-layer splitting measurements in western Anatolia reported a lithospheric-scale anisotropic effect on SKS phases. They explained the fault-parallel orientation of anisotropy (at ~120 km) at the western edge of the NAFZ by the signature of a lithospheric petro-fabric formed during recent deformation. However, for the part of the NAFZ we examine, strike-parallel azimuthal anisotropies are only observed within the upper crustal volume. This indicates that the shear deformation related to the accumulated strain in the earthquake cycle concentrates on shallow cutting faults within the brittle part of the crust. The decrease in strain with depth may be related to the fact that the NAFZ is a relatively young strike-slip fault with a smaller offset (~13–11 Ma and ~80–90 km; Şengör et al., 1985, 2005). Fast wave azimuths oblique to the main strike of the NAFZ observed in almost all single-layer SKS splitting measurements (e.g., Biryol et al., 2010; Paul et al., 2014 and this study) provide further evidence that the NAFZ is unlikely to be a lithospheric-scale shear zone and is only deforming the crust.

Early isotropic model results obtained from the inversions of teleseismic P-wave travel times and Ps converted phases (Kahraman et al., 2015; Papaleo et al., 2017) indicated a widening deformation zone at further depths below the Moho. 2-D variation of our high energy estimates along the N-S profiles from east to west (from FF' to AA'; see Figures 9a–9f) also tends to extend deeper. The high energies mapped from the SNAF to the IZ below the Moho (Figure 10c) can result from a sub-vertical widening deformation zone indicative of an aseismic shear localization zone steadily deforming with slow rates below the seismogenic zone.
### 5.3. Moho Topography

Owing to the dense station distribution and small interstation spacing, our Moho depth estimates from $k = 0$ harmonics along the selected profiles (Figure 9) clearly mark variations in crustal thickness between ~30 and 42 km. These estimates are generally in good agreement with previous studies (e.g., Delph et al., 2015; Eken et al., 2021; Frederiksen et al., 2015; Jenkins et al., 2020; Kahraman et al., 2015; Vanacore et al., 2013).

Different approaches were used to determine the Moho architecture in the study area primarily relying on receiver function analyses. Converted wave analyses in Frederiksen et al. (2015) revealed a discrepancy between Moho depth, Bouguer gravity anomalies, and topography. Similarly, the $k = 0$ harmonic coefficient sections in our study exhibit an inconsistency between crustal thickness and topography (Figure 9). Vertical cross-sections along N-S profiles indicate an undulating Moho boundary, particularly between the branches of NAFZ in the vicinity of the Armutlu Peninsula. We attribute this uneven Moho topography to complex block rotation and non-isostatic equilibrium (Frederiksen et al., 2015; Hergert et al., 2011). The Moho deepens toward the IZ while the topography remains low. The IZ has a notably thick crust of ~42–45 km, which is limited by two N-S striking faults: the WBSF in the west, and the northern strand of the North Anatolian Fault in the south (Figure 14). Jenkins et al. (2020) used CCP stacks of receiver functions to image crustal discontinuities across the Sea of Marmara Region in NW Türkiye. Their work showed Moho depths ranging between 30 and 40 km with an increasing pattern from Sakarya to IZ, which is in line with our Moho estimates (see profile E-E' in Figure 9e). Along the profile E-E', Moho depths reach approximately 34–37 km in the SZ, while tending to decrease to ca. 33 km in the east of the Armutlu Peninsula. In the northern region, the crust of the IZ was found to be thick, reaching a maximum depth of 47 km. These estimates display a general trend of crustal thickening from the SZ to the IZ. Isostatic calculations propose an increased crustal thickness and low elevation below the IZ that should be supported by a lithospheric layer most likely thickened by remnants of an ancient structure. New constraints on the crustal depths incorporated with depth-varying anisotropic features resolved by the harmonic decomposition analysis and model estimates beneath the IZ clearly align with this idea, highlighting the role of the Odessa Shelf (e.g., Akdoğan et al., 2021; Okay et al., 1994) on the origin of this zone.

### 6. Conclusions

Our anisotropic measurements and model results obtained through three different seismological methods implemented in the present study imply that deformation is localized around the fault zones. The orientation of anisotropy in the upper crust near the fault zone is similar to the strike of the NAFZ and past deformation from the IPSZ. When the new constraints on upper crustal anisotropy and existing knowledge of structural geology, seismotectonic properties and geo-chronology are interpreted together, this study favors the SPO anisotropy for Armutlu-Almacık Block, SZ and the LPO beneath the IZ. Relatively diffuse and significant anisotropy in the lower crust can be explained by a widening deformation that can be represented as a weak zone composed of, for example, mineral alignment and/or partially molten lenses elongated along the shearing direction leading to anisotropic signals.

Our analysis and model results for the upper and lower crust near the branches of the NAFZ indicates a complex anisotropic pattern. We found that the NAFZ does not control the observed anisotropic signals starting from ~20 to 25 km, which usually corresponds to the brittle-ductile transition zone in this region. In particular, the observed decoupling between fast axes, which tend to align with the strike of the NNAFZ and SNAFZ in the upper crust, and regional NE-SW oriented fast shear wave directions in the mantle suggests that the faults in this part of the study area do not penetrate deeply into the mantle lithosphere.

### Data Availability Statement

The present work uses data from various temporary and permanent seismic stations obtained from the Dense Array for North Anatolia (DANA, 2012; https://doi.org/10.7914/SN/YH_2012) network initiated within the FaultLab project. The permanent regional network data are operated by the Kandilli Observatory and Earthquake Research Institute-Regional Earthquake-Tsunamic Monitoring Center (KOERI-RETCM) (https://doi.org/10.7914/SN/KO). The digital waveform recordings were retrieved from IRIS DMC (IRIS Data Management Centre) (www.iris.edu). Teleseismic earthquake catalogue data were obtained from the United States Geological Survey (USGS). We generated all figures using Generic Mapping Tool (GMT) (Wessel et al., 2019) and plotting utilities
Acknowledgments

This study is a part of an ongoing Ph.D.
thesis by Derya Keleş under the
supervision of Assoc. Prof. Dr. Tuna Eken.
We have revisited harmonic
decomposition analysis results earlier
presented within the M. Sc. thesis of Derya
Keleş supervised by Prof. Dr. Tuncay
Taymaz. TT, TE and DK would like to
thank Istanbul Technical University-
Research Fund (ITU-BAP, Project
ID:42236 granted to Prof. Dr. Tuncay
Taymaz), International Graduate
Fellowship program of the Scientific and
Technological Research Council of
Türkiye (TUBITAK, EEEAG Project ID:
117E505) for their financial support. We
express our sincere appreciation to Editor
Michael Bostock, Associate Editor, and
two anonymous reviewers for their
invaluable efforts, which substantially
contributed to the refinement of this article.

written in Python language. Figures 1 and 2 were modified using InkScape (http://www.inkscape.org). Figure 14 is created in InkScape and the 3D topographic map in Figure 14 is extracted from GeoMapApp (http://www.geomapapp.org).

References

Akbayram, K., Okay, A. I., & Satir, M. (2013). Early Cretaceous closure of the Intra-Pontide Ocean in western Pontides (Northwestern Turkey). *Journal of Geodynamics*, 65, 38–55. https://doi.org/10.1016/j.jog.2012.05.003

Akbayram, K., Şengör, A. M. C., & Özcan, E. (2016). The evolution of the Intra-Pontide suture: Implications of the discovery of late Cretaceous–Early Tertiary mélanges. In R. Sorkhabi (Ed.), *Tectonic evolution, collision, and seismicity of southwest Asia: In honor of Manuel Berberian's Forty-five years of Research contributions* (Vol. 525). Geological Society of America Special Paper. https://doi.org/10.1130/2016.252518

Akbayram, K., Sorlien, C. C., & Okay, A. I. (2016). Evidence for a minimum 52 ± 1 km of total offset along the northern branch of the North Anatolian Fault in northwest Turkey. *Tectonophysics*, 668–669, 35–41. https://doi.org/10.1016/j.tecto.2015.11.026

Akdoğan, R., Hu, X., Okay, A. I., Topuz, G., & Xue, W. (2021). Provenance of the Paleozoic to Mesozoic siliciclastic rocks of the ΙΖ constrains the timing of the Rheic ocean closure in the eastern Mediterranean region. *Tectonics*, 40(12), e2021TC006824. https://doi.org/10.1029/2021TC006824

Al Azlki, A., Sandvol, E., Seber, D., Barazangi, M., & Turkelli, N. (2004). Pn tomographic imaging of mantle lid velocity and anisotropy at the junction of the Arabian, Eurasian and African Plates. *Geophysical Research Letters*, 30(24). https://doi.org/10.1111/j.1365-246X.2004.02355.x

Al-Azlki, A. I., Seber, D., Sandvol, E., Turkelli, N., Mohamad, R., & Barazangi, M. (2003). Tomographic Pn velocity and anisotropy structure beneath the Anatolian plateau (eastern Turkey) and the surrounding regions. *Geophysical Research Letters*, 30(24), 8043. https://doi.org/10.1029/2003GL017391

Altuncu Poyraz, S., Teoman, U., Türkelli, N., Kahraman, M., Rost, S., Houseman, G., et al. (2015). New constraints on micro-seismicity and stress state in the western part of the North Anatolian Fault Zone: Observations from a dense seismic array. *Tectonophysics*, 656, 190–201. TEC-TO10106R2. https://doi.org/10.1016/j.tecto.2015.06.022

Argus, D. F., & Gordon, R. G. (1991). No-net-rotation model of current plate velocities incorporating plate motion model NUVEL-1. *Geophysical Research Letters*, 18(11), 2039–2042. https://doi.org/10.1029/91GL01532

Audet, P. (2010). Temporal variations in crustal scattering structure near Parkfield, California, using receiver functions. *Bulletin of the Seismological Society of America*, 100(3), 1356–1362. https://doi.org/10.1785/0120090299

Audet, P. (2015). Solid earth layered crustal anisotropy around the San Andreas Fault near Parkfield, California. *Journal of Geophysical Research: Solid Earth*, 120(May), 1–17. https://doi.org/10.1002/2014JB011821

Audoine, E., Savage, M., & Gledhill, K. (1997). Mantle deformation and seismic anisotropy of the South Island, New Zealand (abstract). *Geological Society of New Zealand Miscellaneous Publication*, 95A, 10.

Babuska, V., & Cara, M. (1991). Seismic anisotropy in the earth. In *Modern approaches in geophysics* (Vol. 10). Springer. https://doi.org/10.1007/978-94-011-3600-6

Bernard, R. E., & Behr, W. M. (2017). Fabric heterogeneity in the Mojave lower crust and lithospheric mantle in Southern California. *Journal of Geophysical Research: Solid Earth*, 122(7), 5000–5025. https://doi.org/10.1002/2017JB014280

Bianchi, I., & Bokelmann, G. (2019). Probing crustal anisotropy by receiver functions at the deep continental drilling site KTB in Southern Germany. *Geophysical Prospecting*, 67(9), 2450–2464. https://doi.org/10.1111/1365-2478.12883

Bianchi, I., Park, J., Piana Agostinetti, N., & Levin, V. (2010). Mapping seismic anisotropy using harmonic decomposition of receiver functions: An application to Northern Apennines, Italy. *Journal of Geophysical Research*, 115(12), 1–14. https://doi.org/10.1029/2009JB007061

Biryol, C. B., Beck, S. L., Zandt, G., & Özacar, A. A. (2011). Segmented African lithosphere beneath the Anatolian region inferred from teleseismic P- wave tomography. *Geophysical Journal International*, 184(3), 1037–1057. https://doi.org/10.1111/j.1365-246X.2010.04910.x

Biryol, C. B., Zandt, G., Beck, S. L., Ozacar, A. A., Adıyaman, H. E., & Gans, R. C. (2010). Shear wave splitting along a nascent plate boundary: The North Anatolian Fault Zone. *Geophysical Journal International*, 181, 1201–1213. https://doi.org/10.1111/j.1365-246X.2010.04576.x

Boness, N. L., & Zoback, M. D. (2006). Mapping stress and structurally controlled crustal shear velocity anisotropy in California. *Geology*, 34(10), 825–828. https://doi.org/10.1130/G22309.1

Brownlee, S. J., Schulte-Pelkum, V. A., Mahan, R. J., Condit, C., & Orlandini, O. F. (2017). Characteristics of deep crustal seismic anisotropy from acompilation of rock elasticity tensors and their expression in receiverfunctions. *Tectonics*, 36(9), 1835–1857. https://doi.org/10.1002/2017TC004625

Çakır, Z., Barka, A., De Chabalier, J. B., Armijo, R., & Meyer, B. (2003). Kinematics of the November 12, 1999 (Mw = 7.2) Düzce earthquake deduced from SAR interferometry. *Turkish Journal of Earth Sciences*, 12(1), 105–118.

Castellanos, J. C., Perry-Houts, J., Clayton, R. W., Kim, Y., Stanciu, A. C., Niday, B., & Humphreys, E. (2020). Seismic anisotropy reveals crustal flow driven by mantle vertical loading in the pacific nw. *Science Advances*, 6(28), eabb0476. https://doi.org/10.1126/sciadv.abb0476

Chen, F., Siebel, W., Satir, M., Terzioglu, M., & Saka, K. (2002). Geochronology of the Karadere basement (NW Turkey) and implications for the geological evolution of the ΙΖ. *International Journal of Earth Sciences*, 91(3), 469–481. https://doi.org/10.1007/s00531-001-0239-6

Cochran, E., Vidale, J., & Li, Y.-G. (2003). Near-fault anisotropy following the Hector Mine earthquake. *Journal of Geophysical Research*, 10(B9), 2436. https://doi.org/10.1029/2002JB002352

Confal, J. M., Bezada, M. J., Eken, T., Faccenda, M., Saygin, E., & Taymaz, T. (2020). Influence of upper mantle anisotropy on isotropic P-wave tomography images obtained in the eastern Mediterranean region. *Journal of Geophysical Research: Solid Earth*, 125(8). https://doi.org/10.1029/2019JB018559

Confal, J. M., Eken, T., Tilmann, F., Yolsal-Çevikbilen, S., Çubuk-Sabuncu, Y., Saygin, E., & Taymaz, T. (2016). Investigation of mantle kinematics beneath the Hellenic-subduction zone with teleseismic direct shear waves. *Physics of the Earth and Planetary Interiors*, 261, 141–151. https://doi.org/10.1016/j.pepi.2016.10.012

Confal, J. M., Faccenda, M., Eken, T., & Taymaz, T. (2018). Numerical simulation of 3-D mantle flow evolution in subduction zone environments in relation to seismic anisotropy beneath the eastern Mediterranean region. *Earth and Planetary Science Letters*, 497, 50–61. https://doi.org/10.1016/j.epsl.2018.06.005

Crampin, S. (1978). Seismic wave propagation through a cracked solid: Polarization as a possible dilatancy diagnostic. *Geophysical Journal of the Royal Astronomical Society*, 53, 467–496.
Crampin, S. (1987). The geological and industrial implication of extensive dilatancy anisotropy. *Nature*, **328**(6130), 491–496. https://doi.org/10.1038/328491a0

Crampin, S., Chesnokov, E. M., & Hipkin, R. G. (1984). Seismic anisotropy—The state of the art: II. *Geophysical Journal International*, **76**(1), 1–16. https://doi.org/10.1111/j.1365-246X.1984.tb05017.x

Çubuk-Sabuncu, Y., Taymaz, T., & Fichtner, A. (2017). 3-D crustal velocity structure of western Turkey: Constraints from full-waveform tomography. *Physics of the Earth and Planetary Interiors*, **270**, 90–112. https://doi.org/10.1016/j.pepi.2017.06.014

Delph, J. R., Biryol, C. B., Beck, S. L., Zandt, G., & Ward, K. M. (2015). Shear wave velocity structure of the Anatolian plate: Anomalously slow crust in southwestern Turkey. *Geophysical Journal International*, **202**(1), 261–276. https://doi.org/10.1093/gji/ggv141

Dense Array for North Anatolia (DANA). (2012). International federation of digital seismograph networks. Dataset/seismic network [Dataset]. https://doi.org/10.7914/SN/YH_2012

Di Bona, M. (1998). Variance estimate in frequency-domain deconvolution for teleseismic receiver function computation. *Geophysical Journal International*, **134**(2), 634–646. https://doi.org/10.1111/j.1365-246x.1998.tb07128.x

Eken, T., Bohnhoff, M., Bulut, F., Can, B., & Aktar, M. (2013). Crustal anisotropy in the eastern Sea of Marmara region in northwestern Turkey. *Bulletin of the Seismological Society of America*, **103**(2 A), 911–924. https://doi.org/10.1785/0120120156

Eken, T., Tilmann, F., Mechie, J., Zhao, W., Kind, R., Su, H., et al. (2013). Seismic anisotropy from SKS splitting beneath northeastern Tibet. *Bulletin of the Seismological Society of America*, **103**(6), 3362–3371. https://doi.org/10.1785/0120130054

Eken, T., Tork Qashqai, M., Schiffer, C., Taymaz, T., & Saygin, E. (2021). New insights into crustal properties of Anatolia and its surroundings inferred from P-coda autocorrelation inversions. *Journal of Geophysical Research: Solid Earth*, **126**(12), e2021JB023184. https://doi.org/10.1029/2021JB023184

Emre, Ö., Duman, T. Y., Özalp, S., Şaroğlu, F., Olgun, Ş., Elmacı, H., & Çan, T. (2018). Active fault database of Turkey. *Bulletin of Earthquake Engineering*, **16**(8), 3229–3275. https://doi.org/10.1007/s10518-016-0041-2

Erman, C., Yolsal-Çevikbilen, S., Eken, T., Tilmann, F., Keleş, D., & Taymaz, T. (2022). Constraints on the lithospheric kinematics in the Aegean and western Anatolia unveiled by SKS splitting observations. *Journal of Geophysical Research: Solid Earth*, **127**(12), e2022JB025265. https://doi.org/10.1029/2022JB025265

Evangelidis, C. P. (2017). Seismic anisotropy in the Hellenic subduction zone: Effects of slab segmentation and subslab mantle flow. *Earth and Planetary Science Letters*, **480**, 97–106. https://doi.org/10.1016/j.epsl.2017.10.003

Evangelidis, C. P., Liang, W.-T., Melis, N. S., & Konstantinou, K. I. (2011). Shear wave anisotropy beneath the Aegean inferred from SKS splitting observations. *Journal of Geophysical Research*, **116**(B4), B04314. https://doi.org/10.1029/2010JB007884

Faccenna, C., Becker, T. W., Jolivet, L., & Keskin, M. (2013). Mantle convection in the Middle East: Reconciling Afar upwelling, Arabia indentation and Aegean trench rollback. *Earth and Planetary Science Letters*, **375**, 254–269. https://doi.org/10.1016/j.epsl.2013.05.043

Farra, V., & Vinnik, L. (2000). Upper mantle stratification by P and S receiver functions. *Geophysical Journal International*, **141**(3), 699–712. https://doi.org/10.1046/j.1365-246x.2000.00118.x

Farra, V., Vinnik, L. P., Romanowicz, B., Kosarev, G. L., & Kind, R. (1991). Inversion of teleseismic S particle motion for azimuthal anisotropy in the upper mantle: A feasibility study. *Geophysical Journal International*, **106**(2), 421–431. https://doi.org/10.1111/j.1365-246x.1991.tb03905.x

Fichtner, A., Saygin, E., Taymaz, T., Cupillard, P., Capdeville, Y., & Trampert, J. (2013b). The deep structure of the North Anatolian Fault Zone. *Earth and Planetary Science Letters*, **373**, 109–117. https://doi.org/10.1016/j.epsl.2013.04.027

Fichtner, A., Trampert, J., Cupillard, P., Saygin, E., Taymaz, T., Capdeville, Y., & Villaseñor, A. (2013a). Multiscale full waveform inversion. *Geophysical Journal International*, **194**(1), 534–556. https://doi.org/10.1093/gji/ggt118

Fouch, M. J., & Rondenay, S. (2006). Seismic anisotropy beneath stable continental interiors. *Physics of the Earth and Planetary Interiors*, **158**(2-4), 292–320. https://doi.org/10.1016/j.pepi.2006.03.024

Frederiksen, A. W., & Bostock, M. G. (2000). Modelling teleseismic waves in dipping anisotropic structures. *Geophysical Journal International*, **141**(2), 401–412. https://doi.org/10.1046/j.1365-246X.2000.00090.x

Frederiksen, A. W., Folsom, H., & Zandt, G. (2003). Neighbourhood inversion of teleseismic Ps conversions for anisotropy and layer dip. *Geophysical Journal International*, **155**(1), 200–212. https://doi.org/10.1046/j.1365-246X.2003.02043.x

Frederiksen, A. W., Thompson, D. A., Rost, S., Cornwell, D. G., Gülen, L., Houseman, G. A., et al. (2015). Crustal thickness variations and isostatic disequilibrium across the North Anatolian Fault, western Turkey. *Geophysical Research Letters*, **42**(3), 751–757. https://doi.org/10.1002/2014GL062401

Frothingham, M. G., Mahan, K. H., Schulte-Pelkum, V., Goncalves, P., & Zucali, M. (2023). Confronting solid-state shear bias: Magmatic fabricContribution to crustal seismic anisotropy. *Geophysical Research Letters*, **50**(6), e2022GL102399. https://doi.org/10.1029/2022GL102399

Görür, N., Monod, O., Okay, A. I., Sengör, C. A. M., Tüysüz, O., Yigitbas, E., et al. (1997). Palaeogeographic and tectonic position of the carboniferous rocks of the western Pontides (Turkey) in the frame of the Variscan Belt. *Bulletin de la Societe Geologique de France*, **168**, 197–206.

Gülmez, F., Genç, Ş. C., Keskin, M., & Tüysüz, O. (2013). A post-collision slab-breakoff model for the origin of the Middle Eocene magmatic rocks of the Armutlu-Almacik belt, NW Turkey and its regional implications. *Geological Society - Special Publications*, **372**(1), 107–139. https://doi.org/10.1144/SP372.12

Hatzfeld, D., Karagianni, E., Kassaras, I., Kiratzi, A., Louvari, E., Lyon-Caen, H., et al. (2001). Shear wave anisotropy in the upper mantle beneath the Aegean related to internal deformation. *Journal of Geophysical Research*, **106**(B12), 30737–30753. https://doi.org/10.1029/2001JB000387

Hergert, T., Heidbach, O., Bécel, A., & Laigle, M. (2011). Geomechanical model of the Marmara Sea region - I. 3-D contemporary kinematics. *Geophysical Journal International*, **185**(3), 1073–1089. https://doi.org/10.1111/1365-246X.2011.04991.x

Hirn, A., Jiang, M., Sapin, M., Diaz, J., Neressian, A., Lu, Q. T., et al. (1995). Seismic anisotropy as an indicator of mantle flow beneath the Himalayas and Tibet. *Nature*, **375**(6532), 571–574. https://doi.org/10.1038/375571a0

Holtzman, B., Groebner, N., Zimmerman, M., Ginsberg, S., & Kohlstedt, D. (2003). Stress-driven melt segregation in partially molten rocks. *Geochemistry, Geophysics, Geosystems*, **4**(5). https://doi.org/10.1029/2001GC000258

Houseman, G. A. (2018). Why earthquakes threaten two major European cities: Istanbul and Bucharest. *European Review*, **26**(1), 30–34. https://doi.org/10.1017/s1062798717000448

Hurd, O., & Bohnhoff, M. (2012). Stress and structure-induced shear-wave anisotropy along the 1999 Izmit Rupture, northwest Turkey. *Bulletin of the Seismological Society of America*, **102**(5), 2177–2188. https://doi.org/10.1785/0120110270

Izgi, G., Eken, T., Gaebler, P., Eulenfeld, T., & Taymaz, T. (2020). Crustal seismic attenuation parameters in the western region of the North Anatolian Fault Zone. *Journal of Geodynamics*, **134**, 101694. https://doi.org/10.1016/j.jog.2020.101694
Jenkins, J., Stephenson, S. N., Martínez-Garzón, P., Bohnhoff, M., & Nurlu, M. (2020). Crustal thickness variation across the Sea of Marmara region, NW Turkey: A reflection of modern and ancient tectonic processes. *Tectonics*, 39(7). https://doi.org/10.1029/2019tc005986

Jung, H., & Karato, S. I. (2001). Water-induced fabric transitions in olivine. *Science*, 293(5534), 1460–1463. https://doi.org/10.1126/science.1062235

Kahraman, M., Cornwell, D. G., Thompson, D. A., Rost, S., Houseman, G. A., Türkelli, N., et al. (2015). Crustal-scale shear zones and heterogeneous structure beneath the North Anatolian Fault Zone, Turkey, revealed by a high-density seismometer array. *Earth and Planetary Science Letters*, 430, 129–139. https://doi.org/10.1016/j.epsl.2015.08.014

Kalvoda, J., Babek, O., Fatka, O., Leichmann, J., Melichar, R., Nehyba, S., & Spacek, P. (2008). Brunovistulian terrane (Bohemian Massif, central Europe) from late proterozoic to late Paleozoic: A review. *International Journal of Earth Sciences*, 97(3), 497–518. https://doi.org/10.1007/s00531-007-0183-1

Karato, S. (1998). Seismic anisotropy in the deep mantle, boundary layers and the geometry of mantle convection. *Pure and Applied Geophysics*, 151, 565–587. https://doi.org/10.1007/978-3-0348-8777-9_19

Karato, S.-i., Jung, H., Katayama, I., & Skemer, P. (2008). Geodynamic significance of seismic anisotropy of the upper mantle: New insights from laboratory studies. *Annual Review of Earth and Planetary Sciences*, 36(1), 59–95. https://doi.org/10.1146/annurev.earth.36.031207.124120

Kaviris, G., Fountoulakis, I., Spingos, I., Millas, C., Papadimitriou, P., & Drakatos, G. (2018). Mantle dynamics beneath Greece from SKS and PKS seismic anisotropy study. *Acta Geophysica*, 66(6), 1341–1357. https://doi.org/10.1007/s11600-018-0225-z

Kaya, T., Tank, S. B., Tunçer, M. K., Rokityansky, I. I., Tolak, E., & Savchenko, T. (2009). Asperity along the North Anatolian Fault imaged by magnetotellurics at Düzce, Turkey. *Earth Planets and Space*, 61(7), 871–884. https://doi.org/10.1186/bf03353198

Kendall, J. M. (2000). Seismic anisotropy in the boundary layers of the mantle. In S. Karato, A. Forte, R. Liebermann, G. Masters, & L. Stixrude (Eds.), *Earth's deep interior: Mineral physics and tomography from the atomic to the global scale*. Vol. 117 of geophysical monographs (pp. 133–159). American Geophysical Union.

Kennett, B. L. N., & Engdahl, E. R. (1991). Traveltimes for global earthquake location and phase identification. *Geophysical Journal International*, 105(2), 429–465. https://doi.org/10.1111/j.1365-246X.1991.tb06724.x

Kern, H., & Wenk, H.-R. (1990). Fabric-related velocity anisotropy and shear wave splitting in rocks from the Santa Rosa Mylonite Zone, California. *Journal of Geophysical Research*, 95(B7), 11213–11223. https://doi.org/10.1029/JB095iB07p11213

Klosko, E., Wu, F., Anderson, H., Eberhart-Phillips, D., & McEvilly, T. (1997). Shear-wave splitting and seismic anisotropy in the New Zealand region (abstract). *Eos Trans. AGU*, 78(46). Fall Meet. Suppl., F488.

Koulakov, I., Bindi, D., Parolai, S., Grosser, H., & Milkereit, C. (2010). Distribution of seismic velocities and attenuation in the crust beneath the North Anatolian fault (Turkey) from local earthquake tomography. *Bulletin Seismological Society of America*, 100, 207–224. https://doi.org/10.1785/0120090105

Kufner, S.-K., Eken, T., Tilmann, F., Schurr, B., Yuan, X., Mechje, J., et al. (2018). Seismic anisotropy beneath the Pamir and the Hindu Kush: Evidence for contributions from crustal, mantle lithosphere, and asthenosphere. *Journal of Geophysical Research: Solid Earth*, 123(12), 10727–10748. https://doi.org/10.1029/2018JB015926

Lamarque, G., & Julià, J. (2019). Lithospheric and sub-lithospheric deformation under the Borborema Province of NE Brazil from receiver function harmonic stripping. *Solid Earth Discussions*, 1–20. https://doi.org/10.5194/se-2019-41

Langston, C. A. (1977). Corvallis, Oregon, crustal and upper mantle receiver structure from teleseismic P and S waves. *Bulletin of the Seismological Society of America*, 67(3), 713–724. https://doi.org/10.1785/bssa0670030713

Langston, C. A. (1979). Structure under mount Rainier, Washington, inferred from teleseismic body waves. *Journal of Geophysical Research*, 84(B9), 4749–4762. https://doi.org/10.1029/jb084ib09p04749

León Soto, G., Sandvol, E., Ni, J. F., Flesch, L., Hearn, T. M., Tilmann, F., et al. (2012). Significant and vertically coherent seismic anisotropy beneath eastern Tibet. *Journal of Geophysical Research*, 117(B5), n/a–n/a. https://doi.org/10.1029/2011JB008919

Levin, V., & Park, J. (1997). P-SH conversions in a flat-layered medium with anisotropy of arbitrary orientation. *Geophysical Journal International*, 131(2), 253–266. https://doi.org/10.1111/j.1365-246X.1997.tb01220.x

Licciardi, A., Eken, T., Taymaz, T., Piana Agostinetti, N., & Yolsal-Çevikbilen, S. (2018). Seismic anisotropy in central North Anatolian Fault Zone and its implications on crustal deformation. *Physics of the Earth and Planetary Interiors*, 277(January), 99–112. https://doi.org/10.1016/j.pepi.2018.01.012

Licciardi, A., & Piana Agostinetti, N. (2016). A semi-automated method for the detection of seismic anisotropy at depth via receiver function analysis. *Geophysical Journal International*, 205(3), 1589–1612. https://doi.org/10.1093/gji/ggw091

Liu, Z., & Park, J. (2017). Seismic receiver function interpretation: Ps splitting or anisotropic underplating? *Geophysical Journal International*, 208(3), 1332–1341. https://doi.org/10.1093/gji/ggw455

Long, M. D., & Silver, P. G. (2009). Shear wave splitting and mantle anisotropy: Measurements, interpretations, and new directions. *Surveys in Geophysics*, 30(4–5), 407–461. https://doi.org/10.1007/s10712-009-9075-1

Maupin, V., & Park, J. (2007). Theory and observations - Wave propagation in anisotropic media. In *Treatise on geophysics* (Vol. 1, pp. 289–321). https://doi.org/10.1016/B978-044452748-6.00007-9

McClusky, S., Balassanian, S., Barka, A., Demir, C., Ergintav, S., Georgiev, I., et al. (2000). Global Positioning System constraints on plate kinematics and dynamics in the eastern Mediterranean and Caucasus. *Journal of Geophysical Research*, 105(B3), 5695–5719. https://doi.org/10.1029/1999jb900351

Merry, T. A. J., Bastow, I. D., Kounoudis, R., Ogden, C. S., Bell, R. E., & Jones, L. (2021). The influence of the North Anatolian Fault and a fragmenting slab architecture on upper mantle seismic anisotropy in the eastern Mediterranean. *Geochemistry, Geophysics, Geosystems*, 22(9), 1–36. https://doi.org/10.1029/2021GC009896

Montagner, J. P. (1998). Where can seismic anisotropy be detected in the earth's mantle? In boundary layers.... *Pure and Applied Geophysics*, 151(4), 223–256. https://doi.org/10.1007/s000240050113

Mutlu, A. K., & Karabulut, H. (2011). Anisotropic Pn tomography of Turkey and adjacent regions. *Geophysical Journal International*, 187(3), 1743–1758. https://doi.org/10.1111/j.1365-246x.2011.05235.x

Okay, A., Altiner, D., & Kılıç, A. (2015). Triassic limestone, turbidites and serpentinite—the Cimmeride orogeny in the Central Pontides. *Geological Magazine*, 152(3), 460–479. https://doi.org/10.1017/S0016756814000429

Okay, A. I. (2008). Geology of Turkey: A synopsis. *Der Anschnitt*, 21, 19–42.

Okay, A. I., Şengör, A. M. C., & Görür, N. (1994). Kinematic history of the opening of the Black Sea and its effects on the surrounding regions. *Geology*, 22(3), 267–270. https://doi.org/10.1130/0091-7613(1994)022<0267:khotoo>2.3.co;2

Okay, A. I., & Tüysüz, O. (1999). Tethyan sutures of northern Turkey. *Geological Society - Special Publications*, 156(1), 475–515. https://doi.org/10.1144/GSL.SP.1999.156.01.22
Olive, J. A., Pearce, F., Rondenay, S., & Behn, M. D. (2014). Pronounced zonation of seismic anisotropy in the Western Hellenic subduction zone and its geodynamic significance. *Earth and Planetary Science Letters*, **391**, 100–109. https://doi.org/10.1016/j.epsl.2014.01.029

Özalaybey, S., & Savage, M. K. (1995). Shear wave splitting beneath western United States in relation to plate tectonics. *Journal of Geophysical Research*, **100**(B9), 18135–18149. https://doi.org/10.1029/95jb00715

Papaleo, E., Cornwell, D. G., & Rawlinson, N. (2017). Seismic tomography of the North Anatolian Fault: New insights into structural heterogeneity along a continental strike-slip fault. *Geophysical Research Letters*, **44**(5), 2186–2193. https://doi.org/10.1002/2017GL072726

Park, J., & Levin, V. (2016). Anisotropic shear zones revealed by backazimuthal harmonics of teleseismic receiver functions. *Geophysical Journal International*, **207**(2), 1216–1243. https://doi.org/10.1093/gji/ggw323

Paul, A., Karabulut, H., Mutlu, A. K., & Salaün, G. (2014). A comprehensive and densely sampled map of shear-wave azimuthal anisotropy in the Aegean–Anatolia region. *Earth and Planetary Science Letters*, **389**, 14–22. https://doi.org/10.1016/j.epsl.2013.12.019

Peng, Z., & Ben-Zion, Y. (2004). Systematic analysis of crustal anisotropy along the Karadere-Düzce branch of the North Anatolian fault. *Geophysical Journal International*, **159**(1), 253–274. https://doi.org/10.1111/j.1365-246X.2004.02379.x

Peng, Z., & Ben-Zion, Y. (2005). Spatio-temporal variations of crustal anisotropy from similar events in aftershocks of the 1999 M7.4 İzmit and M7.1 Düzce, Turkey, earthquake sequences. *Geophysical Journal International*, **160**(3), 1027–1043. https://doi.org/10.1111/j.1365-246x.2005.02569.x

Piana Agostinetti, N., & Malinverno, A. (2018). Assessing uncertainties in high-resolution, multifrequency receiver-function inversion: A comparison with borehole data. *Geophysics*, **83**(3), KS11–KS22. https://doi.org/10.1190/geo2017-0350.1

Reilinger, R., McClusky, S., Vernant, P., Lawrence, S., Ergintav, S., Cakmak, R., et al. (2006). GPS constraints on continental deformation in the Africa-Arabia-Eurasia continental collision zone and implications for the dynamics of plate interactions. *Journal of Geophysical Research*, **111**(5), 1–26. https://doi.org/10.1029/2005JB004051

Robertson, A. H. F., & Ustaömer, T. (2004). Tectonic evolution of the intra-Pontide suture zone in the Armutlu Peninsula, NW Turkey. *Tectonophysics*, **381**(1–4), 175–209. https://doi.org/10.1016/j.tecto.2002.06.002

Rost, S., Houseman, G. A., Frederiksen, A. W., Cornwell, D. G., Kahraman, M., Altuncu Poyraz, S., et al. (2021). Structure of the northwestern North Anatolian Fault Zone imaged via teleseismic scattering tomography. *Geophysical Journal International*, **227**(2), 922–940. https://doi.org/10.1093/gji/ggab265

Salaün, G., Pedersen, H. A., Paul, A., Farra, V., Karabulut, H., Hatzfeld, D., et al. (2012). High-resolution surface wave tomography beneath the Aegean-Anatolia region: Constraints on upper-mantle structure. *Geophysical Journal International*, **190**(1), 406–420. https://doi.org/10.1111/j.1365-246x.2012.05483.x

Sambridge, M. (1999a). Geophysical inversion with a neighbourhood algorithm - I. Searching a parameter space. *Geophysical Journal International*, **138**(2), 479–494. https://doi.org/10.1046/j.1365-246X.1999.00876.x

Sambridge, M. (1999b). Geophysical inversion with a neighbourhood algorithm--II. Appraising the ensemble. *Geophysical Journal International*, **138**(3), 727–746. https://doi.org/10.1046/j.1365-246X.1999.00900.x

Sandvol, E., Turkelli, N., Zor, E., Gok, R., Bekler, T., Gurbuz, C., et al. (2003). Shear wave splitting in a young continent-continent collision: An example from eastern Turkey. *Geophysical Research Letters*, **30**(24). https://doi.org/10.1029/2003gl017390

Saunders, P., Priestley, K., & Taymaz, T. (1998). Variations in the crustal structure beneath western Turkey. *Geophysical Journal International*, **134**(2), 373–389. https://doi.org/10.1046/j.1365-246x.1998.00571.x

Savage, M. K., & Silver, P. G. (1993). Mantle deformation and tectonics: Constraints from seismic anisotropy in western United States. *Physics of the Earth and Planetary Interiors*, **78**, 207–228.

Savage, M. S. (1999). Seismic anisotropy and mantle deformation: What have we learned from shear wave splitting? *Reviews of Geophysics*, **37**(1), 65–106. https://doi.org/10.1029/98RG02075

Sayers, C. M. (1994). The elastic anisotropy of shales. *Journal of Geophysical Research*, **99**(B1), 767–774. https://doi.org/10.1029/93jb02579

Schmid, C., van der Lee, S., & Giardini, D. (2004). Delay times and shear wave splitting in the Mediterranean region. *Geophysical Journal International*, **159**(1), 275–290. https://doi.org/10.1111/j.1365-246x.2004.02381.x

Schulte-Pelkum, V., & Mahan, K. H. (2014). Imaging faults and shear zones using receiver functions. *Pure and Applied Geophysics*, **171**(11), 2967–2991. https://doi.org/10.1007/s00024-014-0853-4

Şengör, A. M. C., Görür, N., & Şaroğlu, F. (1985). Strike-slip faulting and related basin formation in zones of tectonic escape: Turkey as a case study. In K. T. Biddle (Ed.), *Strike-slip deformation, basin formation, and sedimentation, Soc. Econ. Paleontol. Miner Spec. Publ. 37 (in honor of J.C. Crowell)* (pp. 227–264). N Christie-Blick. https://doi.org/10.2110/pec.85.37.0227

Şengör, A. M. C., Tüysüz, O., İmren, C., Sakınç, M., Eyidoğan, H., Görür, N., et al. (2005). the North Anatolian Fault: A new look. *Annual Review of Earth and Planetary Sciences*, **33**(1), 37–112. https://doi.org/10.1146/annurev.earth.32.101802.120415

Şengör, A. M. C., & Yılmaz, Y. (1981). Tethyan evolution of Turkey: A plate tectonic approach. *Tectonophysics*, **75**(3–4), 181–241. https://doi.org/10.1016/0040-1951(81)90275-4

Sherrington, H. F., Zandt, G., & Frederiksen, A. (2004). Crustal fabric in the Tibetan Plateau based on waveform inversion for seismic anisotropy parameters. *Journal of Geophysical Research*, **109**(B2), B02312. https://doi.org/10.1029/2002JB002345

Shiomi, K., & Park, J. (2008). Structural features of the subducting slab beneath the Kii Peninsula, central Japan: Seismic evidence of slab segmentation, dehydration, and anisotropy. *Journal of Geophysical Research*, **113**(10), 1–13. https://doi.org/10.1029/2007JB005535

Silver, P. G., & Chan, W. W. (1991). Shear wave splitting and subcontinental mantle deformation conversion at the boundary. *Journal of Geophysical Research*, **96**(B10), 16429–16454. https://doi.org/10.1029/91jb00899

Tadokoro, K., Ando, M., Bariş, Ş., Nishigami, K., Nakamura, M., Balamur Ücer, S., et al. (2002). Monitoring of fault healing after the 1999 Kocaeli, Turkey, earthquake. *Journal of Seismology*, **6**(3), 411–417. https://doi.org/10.1023/a:1020091610510

Tank, S. B., Honkura, Y., Ogawa, Y., Matsushima, M., Oshiman, N., Tunçer, M. K., et al. (2005). Magnetotelluric imaging of the fault rapture area of the 1999 İzmit (Turkey) earthquake. *Physics of the Earth and Planetary Interiors*, **150**(1–3), 213–225. https://doi.org/10.1016/j.pepi.2004.08.033

Tarayoun, A., Audet, P., Mazzotti, S., & Ashoori, A. (2017). Architecture of the crust and uppermost mantle in the northern Canadian Cordillera from receiver functions. *Journal of Geophysical Research: Solid Earth*, **122**(7), 5268–5287. https://doi.org/10.1002/2017JB014284

Taylor, G., Rost, S., Houseman, G. A., & Hillers, G. (2019). Near-surface structure of the North Anatolian Fault Zone from Rayleigh and Love wave tomography using ambient seismic noise. *Solid Earth*, **10**(2), 363–378. https://doi.org/10.5194/se-10-363-2019

Taymaz, T., Yılmaz, Y., & Dilek, Y. (2007). The geodynamics of the Aegean and Anatolia: Introduction. *Geological Society - Special Publications*, **291**, 1–16. https://doi.org/10.1144/SP291.1

Turunçtur, B., Eken, T., Chen, Y., Taymaz, T., Houseman, G. A., & Saygin, E. (2023). Crustal velocity images of north-western Türkiye along the North Anatolian Fault Zone from transdimensional Bayesian ambient seismic noise tomography. *Geophysical Journal International*, **234**(1), ggad082–649. https://doi.org/10.1093/gji/ggad082
Vanacore, E. A., Taymaz, T., & Saygin, E. (2013). Moho structure of the Anatolian Plate from receiver function analysis. *Geophysical Journal International*, 193(1), 329-337. https://doi.org/10.1093/gji/ggs107

Vinnik, L., Oreshin, S., & Erduran, M. (2016). Melt in the mantle and seismic azimuthal anisotropy: Evidence from Anatolia. *Geophysical Journal International*, 205(1), 523-530. https://doi.org/10.1093/gji/ggw021

Walker, K. T., Bokelmann, G. H. R., & Klemperer, S. L. (2004). Shear-wave splitting beneath the Snake River Plain suggests a mantle upwelling beneath eastern Nevada, USA. *Earth and Planetary Science Letters*, 222(2), 529-542. https://doi.org/10.1016/j.epsl.2004.03.024

Wang, H., Huang, Z., Eken, T., Keleş, D., Kaya-Eken, T., Confal, J. M., et al. (2020). Isotropic and anisotropic P-wave velocity structures of the crust and uppermost mantle beneath Turkey. *Journal of Geophysical Research: Solid Earth*, 125(12), e2020JB019566. https://doi.org/10.1029/2020jb019566

Weiss, T., Siegesmund, S., Rabbel, W., Bohlen, T., & Pohl, M. (1999). Seismic velocities and anisotropy of the lower continental crust: A review. In D. Gajewski & W. Rabbel (Eds.), *Seismic exploration of the deep continental crust: Methods and concepts of DEKORP and accompanying projects* (pp. 97-122). Birkhäuser.

Wessel, P., Luis, J. F., Uieda, L., Scharroo, R., Wobbe, F., Smith, W. H. F., & Tian, D. (2019). The generic mapping tools version 6. *Geochemistry, Geophysics, Geosystems*, 20(11), 5556-5564. https://doi.org/10.1029/2019gc008515

Yalcin, N. M., & Yilmaz, I. (2010). Devonian in Turkey—A review. *Geologica Carpathica*, 61(3), 235-253. https://doi.org/10.2478/v10096-010-0014-3

Yamasaki, T., Wright, T. J., & Houseman, G. A. (2014). Weak ductile shear zone beneath a major strike-slip fault: Inferences from earthquake cycle model constrained by geodetic observations of the western North Anatolian Fault Zone. *Journal of Geophysical Research: Solid Earth*, 119(4), 3678-3699. https://doi.org/10.1002/2013JB010347

Yiğitbaş, E., Kerrich, R., Yilmaz, Y., Elmas, A., & Xie, Q. (2004). Characteristics and geochemistry of Precambrian ophiolites and related volcanics from the Istanbul-Zonguldak Unit, Northwestern Anatolia, Turkey: Following the missing chain of the Precambrian South European suture zone to the east. *Precambrian Research*, 132(1-2), 179-206. https://doi.org/10.1016/j.precamres.2004.03.003

Yilmaz, Y., Gözübol, A. M., & Tüysüz, O. (1981). Geology of an area in and around the northern anatolian transform Fault Zone between Bolu and Akyazı. In A. M. Isıkara & A. Vogel (Eds.), *Multidisciplinary approach to earthquake Prediction: Earthquake Prediction Research, Friedr* (pp. 45-66). Vieweg and Sohn.

Yolsal-Çevikbilen, S. (2014). Seismic anisotropy along the Cyprean arc and northeast Mediterranean Sea inferred from shear wave splitting analysis. *Physics of the Earth and Planetary Interiors*, 233, 112-134. https://doi.org/10.1016/j.pepi.2014.05.010

Zhang, S., & Karato, S. (1995). Lattice preferred orientation of olivine aggregates deformed in simple shear. *Nature*, 375(6534), 774-777. https://doi.org/10.1038/375774a0

Zhang, Z., & Schwartz, S. Y. (1994). Seismic anisotropy in the shallow crust of the Loma Prieta segment of the San Andreas Fault system. *Journal of Geophysical Research*, 99(B5), 9651-9661. https://doi.org/10.1029/94JB00241

Zinke, J. C., & Zoback, M. D. (2000). Structure-related and stress-induced shear-wave velocity anisotropy: Observations from microearthquakes near the Calaveras fault in central California. *Bulletin of the Seismological Society of America*, 90(5), 1305-1312. https://doi.org/10.1785/0119990099