The Influence of the North Anatolian Fault and a Fragmenting Slab Architecture on Upper Mantle Seismic Anisotropy in the Eastern Mediterranean

The eastern Mediterranean hosts, within the span of a few hundred kilometers, extensional, strike‐slip, and collision tectonics above a set of fragmenting subducting slabs. Slab roll‐back, toroidal flow, and lithospheric dripping/delamination processes are also believed to be operating. Associated asthenospheric flow and lithospheric deformation are expected to manifest as seismic anisotropy, measurable via study of SKS shear wave splitting. Surprisingly, previous SKS splitting investigations have resolved only long wavelength patterns of anisotropy in the region, interpreting them as large‐scale asthenospheric flow; moreover, no anisotropic signature has been associated with the North Anatolian Fault (NAF), unlike other major strike‐slip plate boundaries worldwide. We present a 29‐year record of SKS splitting observations, revealing hitherto‐unrecognized short‐length‐scale variations in anisotropy, and backazimuthal variations of splitting parameters that attest to multi‐layered anisotropy. Lithospheric anisotropy beneath the NAF exhibits fast directions either fault‐parallel or intermediate between the principle extensional strain rate axis and fault strike, diagnostic of a relatively low‐strained transcurrent mantle shear zone. Elsewhere, anisotropy is consistent with asthenospheric flow through tomographically imaged slab gaps, and driven by Hellenic trench retreat. Evidence for westward flow of asthenosphere driving Anatolian plate motion is lacking. Shorter splitting delay times and nulls in central Anatolia suggest weaker azimuthal anisotropy in the asthenosphere, supporting models that invoke vertical mantle flow patterns (lithospheric dripping/asthenospheric upwelling). Thus, we conclude that the signal of mantle anisotropy more closely reflects the lithospheric deformation, complex slab architecture and geodynamic diversity of the region than previously recognized.

Such tectonic and geodynamic processes each have implications for the strain field in the mantle, and thus may be expected to result in seismic anisotropy via the development of lattice preferred orientation (LPO) in mantle minerals such as olivine (e.g., Zhang & Karato, 1995). Thus, observations of seismic anisotropy in this area should help elucidate how the lithospheric mantle deforms beneath a young strike-slip plate boundary, how the asthenosphere may record a signal of processes such as slab tearing and lithospheric dripping/delamination, and whether Anatolia's tectonic motions are mirrored at depth by a large-scale asthenospheric flow.
Azimuthal seismic anisotropy in the upper mantle is commonly observed via shear-wave splitting of core-refracted seismic phases such as SKS, SKKS, and PKS (hereafter XKS), which record no source-side anisotropy. When an XKS wave traverses an anisotropic medium, it splits into two orthogonally polarized waves, one traveling faster than the other. The azimuth ( E ) of the fast polarization direction, and the delay time ( E t) between the fast and slow arrivals at a seismic station, provide a simple means to characterize the anisotropy (e.g., Silver & Chan, 1991). In Anatolia, most previous XKS shear-wave splitting investigations (Biryol et al., 2010;Paul et al., 2014;Sandvol et al., 2003;Yolsal-Çevikbilen, 2014) have reported remarkably consistent NNE/SSW to NE/SW  E observations over large length scales (Figure 2a), and thus have attributed their results to sub-plate anisotropy caused by large-scale mantle flow. However, localized NW/SE  E observations in southwest Anatolia have been attributed to localized mantle flow through a slab gap (Paul et al., 2014), indicating shorter length-scale asthenospheric processes may be recorded too. Despite this, some have even suggested the anisotropy mainly arises in the mantle transition zone, with the upper mantle being only very weakly anisotropic (Legendre et al., 2021).
Major strike-slip faults are generally associated with highly anisotropic mantle shear zones exhibiting fault-parallel  E (e.g., Duclos et al., 2005;Özalaybey & Savage, 1995;Silver & Chan, 1991;Vauchez et al., 2012, see Table 1), yet previous XKS studies have found no such signal along the NAF (Biryol et al., 2010;Mutlu & Karabulut, 2011;Paul et al., 2014;Sandvol et al., 2003). Possible explanations for this are not limited to the hypothesis that the NAF lithospheric mantle is decoupled from the crustal motions. For example, single XKS splitting observations are the path-integrated result of all anisotropic layers from the core-mantle boundary to the surface. Single, or station-averaged, XKS measurements (e.g., Biryol et al., 2010;Paul et al., 2014;Sandvol et al., 2003) will therefore not necessarily parallel lithospheric or asthenospheric anisotropic fabrics, if both exist. Multiple anisotropic layers yield characteristic backazimuthal variations in XKS shear-wave splitting parameters from which the anisotropic characteristics of different layers (e.g., the lithosphere and the asthenosphere) can be constrained in the presence of good earthquake coverage (Silver & Savage, 1994). Additionally, the fact that the NAF has accrued relatively little slip ( E 85 km; Hubert-Ferrari et al., 2002) compared to other strike-slip faults worldwide (Table 1) may mean it has a comparatively low-strained mantle shear zone. In this case, LPO-derived fast directions may be intermediate between fault strike and the axis of maximum extensional strain ( 1 E ) (Zhang & Karato, 1995).
Here, we use the increased number and deployment time of broadband seismic instruments in the eastern Mediterranean ( Figure 2b), with over 150 new stations since the most recent plate-wide XKS splitting studies of Paul et al. (2014) and Lemnifi et al. (2017), to provide the most comprehensive study of teleseismic shear-wave splitting in the region to date. We assess 6,066 measurements of split and null XKS arrivals in the period 1991-2020, with the specific goal of identifying possible short length-scale asthenospheric flow patterns, and the enigmatic signature of lithospheric anisotropy associated with the North Anatolian Fault.
MERRY ET AL.

Tectonic and Geodynamic Formation and Evolution
The eastern Mediterranean lies on the Alpine-Himalayan orogenic belt, in a zone of north-south convergence between Eurasia, Africa, and Arabia (Figure 1). Regional tectonics have been dominated by this convergence throughout the Cenozoic, involving the near-complete subduction of the Paleo-and Neo-Tethys oceans and the accretion of continental fragments and ophiolites onto the southern margin of Eurasia (e.g., Okay & Tüysüz, 1999;Şengör & Yilmaz, 1981). Anatolia is formed almost entirely of this accreted material, Figure 2. (a) All previous published splitting data, retrieved from the shear-wave splitting database (Wüstefeld et al., 2009). Bars are oriented along the fast direction  E and have length proportional to  E t . In all cases apart from those of Lemnifi et al. (2015Lemnifi et al. ( , 2017 these bars represent station-averaged splitting parameters; those of Lemnifi et al. (2015Lemnifi et al. ( , 2017  manifest as approximately east-west trending terranes, separated by suture zones (Figure 1a: e.g., Okay & Tüysüz, 1999;Şengör & Yilmaz, 1981).
The northernmost of these terranes is the Rhodope-Pontide fragment (Şengör & Yilmaz, 1981), which includes Paleotethys suture zones from collisional events prior to the mid-Jurassic (e.g., Şengör & Yilmaz, 1981) and is bound to the south by the Intra-Pontide suture, a Neotethyan suture zone of Eocene age (Akbayram et al., 2016). South of the Intra-Pontide suture is the Sarkarya continent, which also contains a Paleotethys suture zone. The Neotethyan İzmir-Ankara-Erzincan suture zone separates the Sarkarya continent and Rhodope-Pontide fragment from the Kırşehir Block and Menderes-Taurus Block to the south (Figure 1a). These two blocks underlie most of western and central Anatolia, and it is debated whether they originate from a continuous continental fragment or from two fragments previously separated by now-subducted oceanic lithosphere (e.g., Görür et al., 1984;Şengör et al., 2019;van Hinsbergen et al., 2016). Eastern Anatolia lacks these units, and is dominated by subduction-accretion complexes (Şengör et al., 2008;Şengör & Yilmaz, 1981).
Miocene-to-Recent volcanism and uplift in eastern Anatolia has been attributed to slab steepening and break-off processes, resulting in the removal of almost all lithospheric mantle (Bartol & Govers, 2014;Şengör et al., 2003). Uplift of the Central Anatolian Plateau, with volcanism occurring at its northern and southern margins, has been explained by invoking a lithospheric drip centered beneath the Kırşehir Block ( Figure 1a; Göğüş et al., 2017), or a continuation of slab-steepening and break-off as in eastern Anatolia (Bartol & Govers, 2014). Additionally, very rapid uplift rates  E 3 mm/yr at the southern margin of central Anatolia since 0.5 Ma may be a response to incipient break-off of the Cyprus slab (Öğretmen et al., 2018), which has been imaged in seismic tomography (Kounoudis et al., 2020;Portner et al., 2018). Finally, a transition in western Anatolia from subduction zone-style volcanism to intraplate-style volcanism from ca. 10 Ma indicates an increasing asthenospheric source of melt in this area (Dilek & Altunkaynak, 2010;McNab et al., 2018). This may be related to the acceleration of Hellenic trench retreat at this time, with the opening of the gap between the Cyprus and Aegean slabs beneath the Isparta Angle (Figure 1b), which is evident today in seismic tomographic models (Biryol et al., 2011;Kounoudis et al., 2020).

Previous Investigations of Seismic Anisotropy and Mantle Deformation
In the Aegean and far-western Anatolia, XKS splitting investigations ( Figure 2a: Evangelidis et al., 2011;Hatzfeld et al., 2001;Paul et al., 2014;Schmid et al., 2004), anisotropic Pn tomography (Mutlu & Karabulut, 2011), anisotropic P-wave tomography Wei et al., 2019) and anisotropic surface wave tomography (Endrun et al., 2011) have shown fast directions that parallel ongoing crustal stretching directions and mineral lineations in exhumed shear zones (L. Jolivet et al., 2009), consistent with interpretations that deformation is vertically coherent throughout the crust and mantle lithosphere (L. Jolivet et al., 2009Jolivet et al., , 2013Kreemer et al., 2004). This is in contrast with the rest of Anatolia, where early shear-wave splitting investigations indicated consistently NE/SW-oriented  E , despite the complex surface strain field (Biryol et al., 2010;Sandvol et al., 2003). In particular, these authors concluded that the major faults (NAF, EAF) had no apparent effect on upper mantle anisotropy, and interpreted their results as representing an asthenospheric flow field. 7 of 26 A subsequent regional XKS splitting study by Paul et al. (2014) augmented these results with observations throughout the rest of Anatolia, with  in south-west Anatolia was interpreted as toroidal flow through the gap between the Cyprus and Aegean slabs (Paul et al., 2014), as subsequently modeled by Confal et al. (2018). In the rest of Anatolia, Biryol et al. (2010) had attributed their XKS splitting results to south-westward asthenospheric flow induced primarily by rollback of the Aegean and Cyprus slabs. Paul et al. (2014) point out that evidence for significant rollback of the Cyprus slab is lacking, and that Aegean rollback alone is unlikely to induce south-westward asthenospheric flow in eastern Anatolia. Paul et al. (2014) prefer an interpretation following the mantle convection models of Faccenna et al. (2013), who suggest that Hellenic trench retreat and flow from the distant Afar plume upwelling combine to create a toroidal asthenospheric flow regime matching the tectonic motions. A northward advance of Afar plume material beneath Arabia and the Levant has been suggested to explain a variety of geological and geophysical evidence, including northward progression of volcanism and uplift throughout western Arabia (Camp & Roobol, 1992;Wilson et al., 2014), N/S oriented  E in shear-wave splitting studies on the Arabian plate (Hansen et al., 2006;Kaviani et al., 2013), and a continuous region of low upper-mantle seismic velocities (e.g., Hansen et al., 2006), while westward flux of this material beneath Anatolia is supported by westward-propagating uplift since the Miocene (McNab et al., 2018).
At the NAF, in addition to the lack of fault-parallel  E observations in the aforementioned XKS splitting studies (Figure 2a), fault-oblique uppermost mantle anisotropic fast directions have been observed in anisotropic Pn tomography (Mutlu & Karabulut, 2011), as has a sharp change in fast directions across the Moho in anisotropic P-wave tomography . Given the evidence for a thinned lithosphere, the apparent lack of influence of either the NAF or EAF on upper mantle anisotropy, and the discrepancy between plate motion and anisotropic fast directions, it has been suggested that the Anatolian crust and mantle are mechanically decoupled (L. Jolivet et al., 2013;Mutlu & Karabulut, 2011;Wang et al., 2020). On the other hand, low-velocity upper mantle structures beneath the NAF have been imaged in tomographic studies and interpreted as signatures of a plate-scale shear zone (Fichtner et al., 2013;Papaleo et al., 2017Papaleo et al., , 2018, favoring a hypothesis of coherent crust and lithospheric mantle deformation. More recent shear-wave splitting work by Lemnifi et al. (2017), reanalyzing many of the stations in Anatolia used by Paul et al. (2014), suggests more complexity is present, with evidence of backazimuthal variations in splitting parameters in northern Turkey. The authors propose two anisotropic layers, with the upper contained entirely within the crust. Depth-dependent anisotropy is also suggested beneath the central NAF by receiver function analysis (Vinnik et al., 2016) and there is evidence for backazimuth-dependent splitting parameters at station ISP in south-west Anatolia, interpreted as two anisotropic layers by Şapaş and Boztepe-Güney (2009). Contrastingly, recent anisotropic Rayleigh wave tomography is used by Legendre et al. (2021) to suggest that anisotropy is weak throughout the Anatolian upper mantle, and XKS splitting therefore arises in the mantle transition zone.

Seismological Data Used in the Study
We analyze XKS shear-wave splitting at publicly available seismic stations in the eastern Mediterranean from 1991 to 2020, and the temporary TROODOS network which operated on Cyprus in 2017-2019 (Bastow et al., 2016). While many of these stations have been used in previous shear-wave splitting studies (Figure 2b; Biryol et al., 2010;Confal et al., 2016;Evangelidis et al., 2011;Hatzfeld et al., 2001;Kaviani et al., 2011Kaviani et al., , 2013Lemnifi et al., 2017;Olive et al., 2014;Paul et al., 2014;Sandvol et al., 2003;Yolsal-Çevikbilen, 2014), we re-analyze all data to ensure a consistent approach and, for long-running, permanent stations, to extend the temporal coverage of constraints. Additionally, a recent analysis of KOERI seismic stations (network KO in Figure 2b), measuring polarizations of P waves and Rayleigh waves (Büyükakpınar et al., 2021), shows that many of the permanent stations in Turkey included in our analysis are, or have been at some stage, misaligned. For all stations with a misalignment estimate of   5 , we have corrected horizontal component azimuths prior to analysis.

of 26
We present, for the first time, XKS splitting results using temporary networks MARSite (Ergintav et al., 2013) and DANA (Altuncu Poyraz et al., 2015) in north-west Anatolia, the CD-CAT network (Sandvol, 2013) in Central Anatolia, as well as many stations from the SIMBAAD network (Paul et al., 2013) in western Anatolia ( Figure 2b). There have also been several additions to the permanent networks we use (notably the KOERI network) since previous studies were published. We also process data from the Georgia, Azerbaijan and Iraq national networks, which are not included in previous shear-wave splitting studies.

Data Processing
We identified 36,792 XKS arrivals by visual inspection from earthquakes of magnitude  E 5.9 at epicentral distances   88 E from 684 stations (Figure 4a) throughout the region in the period 1991-2020. All waveforms were filtered using a zero phase Butterworth bandpass filter with corner frequencies of 0.04 and 0.3 Hz.
We performed splitting analysis on each identified XKS arrival, using the semi-automated workflow of Teanby et al. (2004), based on the Silver and Chan (1991) eigenvalue minimization method. Horizontal component seismograms are rotated and time-shifted to minimize the second eigenvalue of the covariance matrix for particle motion within a time window around the shear wave arrival, effectively linearizing the particle motion and minimizing tangential-component shear wave energy (Figure 3a-3c). If the particle motion is already linear and there is no tangential-component energy, this is called a 'null' measurement and indicates either that the fast direction is parallel or perpendicular to the radial shear wave direction, or that there is no azimuthal anisotropy.
Traditional XKS splitting analysis takes a single, manually picked, shear-wave analysis window. Following the Teanby et al. (2004) method, however, splitting analysis is performed on 100 different time windows and cluster analysis is used to select splitting parameters that are stable over many windows (see Teanby et al., 2004, for details) (Figures 3e and 3f). We accept only those results with stable clusters, visibly linearized particle motion, and errors of   15 in  t. Errors in each final individual station measurement are calculated via an F-test to obtain the 95% confidence interval (Figure 3d).

Summary of Results and Comparisons to Previous Studies
Accepted splits were recorded at 516 of the 684 stations analyzed, including 4296 SKS, 661 PKS and 340 SKKS arrivals (Figure 4a). We recorded 769 null results across 251 stations; there were 38 stations where only null results were observed. The average number of accepted splits at each station was nine, with a maximum of 127 at station EIL. Across the network, most splits have 0.5 s    E t 2.5 s. Each accepted split is recorded in Table S2, and each null in Table S3.
For stations with multiple measurements, we stack the associated error surfaces (Figure 3d), weighted according to the signal-to-noise ratio, to obtain a single set of splitting parameters for each station (see Wolfe & Silver, 1998). At stations where two or more splits are recorded, we include any null results in the stacking process, since they help to constrain 

E
. The results of the stacking process are shown in Figure 4b and  Figures 2a and 4). However, with the greater number of stations and lengthened earthquake catalog, we do find additional evidence for lateral variations and backazimuth dependence in splitting parameters, especially in northwest and south-central Anatolia. Thus, where previous studies (Biryol et al., 2010;Paul et al., 2014;Sandvol et al., 2003) have proceeded on the assumption of a single, horizontal, homogeneous anisotropic layer, we can consider more complex scenarios. Next, we summarize the observations that can be used to assess the relative contributions of lithospheric and asthenospheric anisotropy, informing subsequent two-layer modeling work described in Section 4.3, and finally describe where our results agree with other predictors of seismic anisotropy.

Lateral Variations in Splitting Parameters
The observation that  E t varies throughout the region, with   E t 1 s through much of central Anatolia and   E t 1.5 s in parts of northwest Anatolia (Figure 5a), indicates that there are changes in strength of azimuthal anisotropy and/or anisotropic layer thickness between these regions. The magnitude of  E t, in light of lithospheric thickness estimates, provides a first-order constraint on whether the splits can be produced within the lithosphere. Assuming an upper mantle shear-wave velocity of 4.48 km/s (Kennett et al., 1995) and 4% anisotropy (Savage, 1999), a 100 km thick anisotropic layer would produce   E t 0.9 s. Thus, lithospheric thickness estimates of 50-100 km (Figure 5b: e.g., Hoggard et al., 2020;Kind et al., 2015;McNab et al., 2018;Reid et al., 2019) imply that the lithosphere cannot generally be the sole source of splitting, with   E t 1 s throughout most of the region, aside from central Anatolia (Figure 5a). On the other hand, the increased spatial density of observations, especially in northwest Anatolia with the addition of the DANA and MARSite networks (Figure 2), allows us to present evidence for lateral variations in splitting parameters unidentified by previous investigations. The best example of this is around the North Anatolian Fault (NAF) at   30 E, where we observe changes in  E of   60 over  E 20 km ( Figure 6). Following Rümpker and Ryberg (2000), the minimum XKS Fresnel zone width at 60 km depth (the likely minimum LAB depth; e.g., McNab et al., 2018;Reid et al., 2019) is  E 50 km. Thus, changes in splitting parameters observed on a horizontal scale of  E 20 km must arise in the lithosphere, since Fresnel zones overlap in the asthenosphere.
Some of the lithospheric delay times may arise in the crust, but prior investigations have revealed crustal anisotropy to be relatively weak in Anatolia. Local earthquake shear-wave splitting has resulted in delay times commonly    (Licciardi et al., 2018). Other seismological methods have revealed shear-wave crustal anisotropy to be E 2% on average, or E 3% in the upper crust (Legendre et al., 2020;Taylor et al., 2019;Vinnik et al., 2016). Thus, crustal anisotropy is unlikely to substantially affect XKS splitting parameters, with a maximum contribution of   E t 0.2 s.

Backazimuthal Variations and Two-Layer Modeling
Stacking results at each station ( Figure 4b, Section 4.1) implicitly imposes the assumption of a single, horizontal, homogeneous anisotropic layer. Where this assumption is invalid, we expect to see splitting parameters that vary according to the backazimuth of the incoming shear wave (Silver & Savage, 1994). This is observed in some locations, notably around the Dead Sea Transform, many parts of the NAF, and in some areas of southern Anatolia (Figure 7). In the case of two anisotropic layers, we expect a characteristic  90 E periodicity in splitting parameters as a function of backazimuth (Silver & Savage, 1994), and where we observe this we model such variations using the equations of Silver and Savage (1994), following Walker and Wookey (2012). We use a grid-search misfit minimization approach similar to Liddell et al. (2017) Kaviani et al. (2013), where  E is the standard deviation in observed splitting parameters (see Text S1 for further details). We assume an anisotropic strength of 4% (e.g., Savage, 1999) in order to test thickness; however, we acknowledge that anisotropic strength trades off against layer thickness.
An excellent example of a characteristic two-layer variation is in the region of the Dead Sea Transform, where the station EIL has recorded 127 splits across a wide range of backazimuths (Figures 7a and S1).   Figure S3). We also model station KZIT, around 90 km west of the DST, and while the lower layer  l E is less well constrained, our  u =  10 E remains similar to regional values ( Figure S4).
Best-fitting two-layer models at all stations where modeling was possible are shown in Figures S1-S38, summarized in Table S1; equivalent splitting parameters for upper layers are displayed in Figure 8a and lower layers in Figure 8e. At most stations across Anatolia, there is either little evidence for backazimuthal variations in splitting parameters, or insufficient backazimuthal coverage to model the variations. In particular, instruments from temporary deployments (open polygons in Figure 2b) lack sufficient backazimuthal coverage to resolve multi-layer anisotropy, notably meaning we are unable to perform this analysis for most of the central and eastern NAF, with the exception of stations ILGA (which we combine with nearby temporary stations BEDI, SYUN, and ALIC into the hybrid station CNAF; Figures 7b and S5), ERZN ( Figure S6) and YEDI ( Figure S7). Nevertheless, backazimuthal variations are consistent with two layers of anisotropy at all stations with sufficient coverage for modeling around the NAF ( Figures S5-S38), and we are able to constrain best-fitting two layer models at many stations (e.g., Figures 7b and 7c), supporting the presence of both lithospheric and asthenospheric anisotropy here (Section 5.2). Figure 9 shows both station-stacked results and best-fitting two-layer models for stations associated with the NAF, defined here as those within 30 km (i.e., approximately half a minimum Fresnel zone width at lithospheric mantle depths: Rümpker & Ryberg, 2000) of the North Anatolian Shear Zone (NASZ), defined by Şengör et al. (2005) as the crustal extent of NAF-associated dextral motion.
More complex multi-layered anisotropy beneath the central NAF has been suggested from joint inversion of shear-wave splitting and receiver functions across a wide area of northern Anatolia (Vinnik et al., 2016), with an E/W-oriented lower asthenosphere beneath N/S upper asthenosphere and NE/SW lithosphere, but the backazimuthal variations we observe are inconsistent with this (Figure 7b, Text S1.1). We favor a simpler two-layer model with  u , similar to other stations at the NAF.
By contrast, backazimuthal variations evident at some stations in southern and southwestern Anatolia do not show the classic  90 E periodicity associated with two-layer anisotropy (e.g., station ISP, Figure 7d, and stations ANTB, ELL, DALY, and IKL, Figures S39-S42). This may be caused by lateral variations at depth (e.g., Alsina & Snieder, 1995), suggesting complex asthenospheric anisotropy may be present here (see Section 5.5). At station ISP, the two-layer model previously proposed by Şapaş and Boztepe-Güney (2009) yields

Comparison With Surface Strain and Other Anisotropic Constraints
Given the evidence outlined in Sections 4.2 and 4.3 for lithospheric anisotropy, we investigate possible relationships between ongoing tectonic strain and anisotropy by comparing our data with published models of

10.1029/2021GC009896
14 of 26 current interseismic strain rate at the Earth's surface, calculated using geodetic data (Figures 6, 8a, and 9). We use two models: the Global Strain Rate Model (GSRM), a geodetic model that relies mainly on GPS data (Kreemer et al., 2014), and the more recent regional model of Weiss et al. (2020), henceforth referred to as Weiss20, which additionally uses satellite radar interferometry data to obtain more complete coverage in Anatolia. We generally use Weiss20 where possible due to its increased data coverage. The primary ) and lower ( l E ) anisotropic layers from two layer models. The approximate strike of the North Anatolian Fault is marked, as are the average axes of maximum extensional strain rate within the same area (Weiss et al., 2020).
MERRY ET AL.

10.1029/2021GC009896
16 of 26 Regional asthenospheric anisotropy is constrained by the anisotropic P-wave tomography of Wei et al. (2019), who use their model to predict XKS splitting parameters (Figure 8e). Regions where observed  E agrees with their predicted  E include eastern Anatolia, where  E is remarkably consistent (Figure 4) despite the highly variable active strain and multiple terranes, and the northern Aegean and northwest Anatolia (Figure 8f). The best-fitting lower layer  l E in the DST region (Section 4.3; Kaviani et al., 2013) also matches well with those of Wei et al. (2019) (Figure 8e). Those regions where our results disagree with those of Wei et al. (2019)-that is, much of north, central and south-west Anatolia-are likely to have either significant lithospheric anisotropy, or asthenospheric anisotropy that varies on shorter length scales than are resolved in their model. These regions represent the locations of, respectively, the rapidly straining NAF, the complex Cyprus slab, and the gap between the Cyprus and Aegean slabs.

Discussion
The results presented here require a re-examining of the broad consensus (Biryol et al., 2010;Legendre et al., 2021;Paul et al., 2014;Sandvol et al., 2003;Yolsal-Çevikbilen, 2014) that regional-scale, sub-plate processes dominate the anisotropy beneath Anatolia. On the one hand, when considering station-stacked results (Figure 4b), the observations that led previous authors to discount lithospheric anisotropy as a significant source of splitting remain: splitting parameters for the most part are relatively undisturbed by major lithospheric terrane boundaries and faults. On the other hand, the increased spatial coverage of observations in northwest Anatolia reveals the importance of shorter-scale processes around the NAF (Section 4.2, Figures 4  and 6), where backazimuthal variations support the presence of both lithospheric and asthenospheric anisotropy (Section 4.3). Furthermore, our improved coverage of central Anatolia reveals variable  E ( Figure 4) and a  E t that is markedly lower than surrounding regions ( Figure 5). Such observations are incompatible with the homogeneous, sub-plate source of anisotropy generally invoked by previous studies, and particularly with the suggestion by Legendre et al. (2021) that the entire upper mantle is only weakly anisotropic.
In the following sections, we consider first how well lithospheric anisotropy can explain the observed splitting in actively straining regions. We then consider the influence of the complex environment of subducting slabs on splitting parameters, and the implications of our results for ideas surrounding small-versus broadscale asthenospheric flow.

Western Anatolia and the Aegean
In the extensional region of western Anatolia and the central Aegean, splitting parameters are generally independent of backazimuth (e.g., stations DKL, KULA, and SMG, Figures S43-S45), consistent with a co-orientation of  E in the lithosphere and asthenosphere. In extensional regimes, lithospheric anisotropy as a result of vertically coherent deformation is expected to result in  E paralleling the axis of maximum extensional strain rate,  1 E , at the surface (Savage, 1999;Silver & Chan, 1991). Thus, the observation that  E is parallel with  1 throughout this region (Figures 8a and 8b) is explained well by vertically coherent deformation related to back-arc extension, as has already been suggested (L. Jolivet et al., 2009Jolivet et al., , 2013. However, anisotropy is unlikely to be confined to the lithosphere, given   E t 1.2 s throughout the region (Figure 5a). A similar  E in the asthenosphere can be attributed readily to rollback-induced flow toward the retreating Hellenic trench (e.g., Confal et al., 2018;Evangelidis et al., 2011). This also agrees well with the lower layer of two-layer models in the westernmost part of the NAF, just to the north of this region (Figure 9). The southward decrease in  E t (Figure 5a) may indicate either a decreasing anisotropic layer thickness, or a rotation of fast axes toward the vertical, as asthenospheric flow becomes more vertical closer to the sinking slab.

The North Anatolian Fault
One of the key observations motivating previous authors to only consider sub-plate anisotropy has been the lack of parallelism of  E with the NAF (Biryol et al., 2010;Paul et al., 2014;Sandvol et al., 2003); this is generally reproduced in our station-stacked results (Figure 4b). By contrast, fault-parallel  E is seen at almost

Influence of Complex Slab Geometry
In southern Anatolia, our increased data coverage indicates much more variable splitting parameters than previously recognized, with a clear influence of the complex geometry of subducting slabs. The edges of slab-related fast wavespeed anomalies coincide with abrupt variations in splitting parameters: superimposing splits at their 150 km depth piercing points on the 150 km depth-slice of the S-wave tomographic model of Kounoudis et al. (2020) demonstrates that  E diverts around slab edges and through slab gaps ( Figure 10). This includes the NW/SE fast directions in south-west Anatolia, already proposed by Paul et al. (2014) to indicate toroidal flow through the Aegean-Cyprus slab gap, evident in the Kounoudis et al. (2020) model as slow wavespeeds between the faster anomalies of the Cyprus and Aegean slabs ( Figure 10).
Lateral variations around slab edges are evident even at individual stations such as ISP, above the western edge of the Cyprus slab, where splitting parameters vary depending on whether the backazimuth is to the east (slab) or west (slab-gap) (Section 4.3, Figures 7d and 10). The observed sharp lateral changes in  E around slab edges could be an indication of toroidal flow, or could reflect more complex anisotropic characteristics of either the slab itself or the hydrated mantle wedge.
A region of predominantly null results from several backazimuths around the Cappadocia Volcanic Province (CVP), central Anatolia (Figure 10), suggests azimuthal anisotropy in the asthenosphere is weak or absent. This coincides with a horizontal tear imaged in the Cyprus slab (Kounoudis et al., 2020), suggesting azimuthally isotropic, vertical flow through the opening gap. River profile inversion indicates that, along with the Isparta angle, this area has hosted some of the highest uplift rates ( E 0.25 mm/yr) in Anatolia over the last 4 Myr (Figure 10 inset: McNab et al., 2018). Additionally, very rapid uplift (>3 mm/yr) of central Anatolia's southern margin (pink star in Figure 10) is evident over the last 0.5 Myr from uplifted marine terraces (Öğretmen et al., 2018), while a switch in basalt geochemistry from calc-alkaline in the Neogene to Na-alkaline in the Quaternary in the CVP has been attributed by Aydin et al. (2014) to a change in mantle melt source from predominantly lithospheric to predominantly asthenospheric. Collectively, the coincidence of null splitting with focused rapid uplift, volcanism, recent changes in basalt geochemistry, and a tomographically imaged slab tear (Kounoudis et al., 2020) constitute compelling evidence for upward flow of sub-slab asthenosphere through the tear beneath the CVP. Additional local vertical flux of asthenosphere through the slab may be indicated beneath the south coast of Anatolia, where null results from a range of western (but not eastern) backazimuths recorded at station IKL (Figures 10 and S42) coincide with both a slower region in the tomography (Kounoudis et al., 2020) and the rapid uplift rates recorded by Öğretmen et al. (2018).  (Figures 4 and 5a). This region coincides with the southwards-younging Kır-ka-Afyon-Isparta volcanic fields that are proposed to be related to the development of the Cyprus-Aegean slab gap (Dilek & Altunkaynak, 2009;Karaoğlu & Helvaci, 2014), which is clearly present at the southern edge of this region ( Figure 10). Thus, the strong anisotropic signal is likely to be an expression of the development of this tear and subsequent lateral mantle flow.
Splitting results in the forearcs of the Cyprus and Aegean subduction zones are highly variable and backazimuth-dependent ( Figure 10; see Figure S48 for an example at station CSS in Cyprus). This may reflect differing subslab, intraslab and supraslab anisotropy, combined with localized features such as the layer of serpentinized mantle suggested in the Aegean forearc by Olive et al. (2014). Resolution of such details, however, is best achieved in conjunction with high-resolution seismic images as in Olive et al. (2014), or other data with greater vertical resolution than XKS studies, including local earthquake shear-wave splitting (e.g., Evangelidis, 2017), so we do not explore these stations further here.

Implications for Asthenospheric Flow and Anatolian Plate Driving Forces
The rapid westward motion of Anatolia has been attributed by some to a toroidal asthenospheric flow field: northward through western Arabia, westward through Anatolia and then southwestward in the Aegean (e.g., Faccenna et al., 2013;Le Pichon & Kreemer, 2010). A hypothesis of SW-directed flow in the Aegean matches our results (Section 5.1). To the east,  E in eastern Anatolia and the lower layer of the Dead Sea Transform consistently parallels the trend of an elongate low-wavespeed anomaly at 150 km depth ( Figure 11; Schaeffer & Lebedev, 2013) that has been associated with 30 Ma-Recent base-of-lithosphere asthenospheric flow from the Afar hotspot (e.g., Camp & Roobol, 1992) and/or the putative Jordan hotspot (Chang & Van der Lee, 2011). However, rather than rotating counterclockwise from Arabia to Anatolia,  E rotates in the opposite sense to become NE/SW-oriented (Figure 8c), implying that, if there is a large scale northward flow, it is likely oriented northeast toward the Caucasus (Figure 11). While Faccenna et al. (2013) indicate that counterclockwise toroidal flow can induce NE-SW fast directions in parts of central and eastern Anatolia, the widespread consistency of  E in eastern Anatolia is not fully reproduced in their models.
Evidence for westward-directed mantle flow in Anatolia is lacking. E/W oriented  E is only observed locally beneath the Kırşehir Block in Central Anatolia (Figure 4). Seismological and petrological evidence for a relatively thick lithosphere here (Fichtner et al., 2013;Reid et al., 2019) means that the associated short delay times (  1 E t s) could be accrued entirely within the lithosphere. While the backazimuthal coverage is poor here, there is no evidence for backazimuth-dependence of splitting parameters and thus no reason to invoke multiple anisotropic layers (e.g., stations YOZ and AT06, Figures S49 and S50). A local sub-parallelism with suture zones also lends credence to a lithospheric anisotropy hypothesis (Figure 8a). To the north, in the central NASZ, an eastward decrease in  E t coincides with station-stacked  E becoming somewhat closer to the fault strike (Figures 9b and 9c), perhaps indicating weaker asthenospheric splitting and thus 10.1029/2021GC009896 21 of 26 a greater relative contribution from shear-zone, lithospheric anisotropy (Section 5.2). Additionally, there are stations with predominantly null results from a range of backazimuths, both in the south (CVP; see Section 5.4) and northeast of central Anatolia (Figure 4a). Small  E t and null results imply weaker azimuthal anisotropy in the asthenosphere, consistent with geodynamic models proposing asthenospheric upwellings or sinking lithospheric drips (e.g., Göğüş et al., 2017), and vertical flow through the Cyprus slab tear in the south (Section 5.4).
Finally, we argue that a lack of westward asthenospheric flow is compatible with the uplift and volcanism history of Anatolia, despite arguments on the basis of the westward progression of uplift over the last 15 Ma and westward-decreasing mantle potential temperature (McNab et al., 2018). In fact, the areas of greatest recent uplift and volcanism in Anatolia are those where we infer an introduction of subslab material into the Anatolian asthenosphere, either vertically through the Cyprus slab tear or laterally through the slab gap ( Figure 10, Section 5.4). This suggests that any hot material being introduced originates below the African slabs, perhaps part of a regional North African/Mediterranean upwelling (Nikogosian et al., 2018) or simply a background 'spoke-pattern' mantle convection pattern that has been suggested to occur in the region (Lees et al., 2020;McKenzie, 2020).
Given the lack of evidence for westward flow of Anatolian asthenosphere, and the compelling explanations for azimuthal anisotropy (or lack thereof) from localized geodynamic processes, we suggest that the westward motion of Anatolian lithosphere is not driven by drag from the asthenosphere. Thus other forces such as gradients in gravitational potential energy and plate boundary forces are likely the dominant forces driving the tectonics here (Capitanio, 2016;England et al., 2016;Özeren & Holt, 2010).

Conclusions
We have used a 29-year record of split and null XKS phases to constrain azimuthal seismic anisotropy in the eastern Mediterranean that is more variable, both laterally and vertically, than previously recognized. While asthenospheric anisotropy is likely the most important source of anisotropy in much of the region (as previously surmised), with the benefit of denser station coverage and greater backazimuthal coverage we are able to identify a shallower contribution to splitting that provides new constraints on lithospheric dynamics in the region. Below the North Anatolian Fault, we find compelling evidence for lithospheric anisotropy, and in many cases a characteristic variation in splitting parameters with backazimuth that indicates multiple anisotropic layers. Lithospheric fast polarization directions ( E ) lie either close to fault strike, or between fault strike and the axis of maximum extensional strain rate ( 1 E ), consistent with a relatively lowstrain strike-slip mantle shear zone. There is thus no need to invoke a decoupled crust-mantle system in the lithosphere. Along-strike variations in the relationship between splitting parameters and  1 E /fault strike may reflect a variation in shear zone width and strain localization.
These results, in combination with the additional variability in splitting parameters revealed by the increased coverage of data in central and southern Anatolia, diminish the importance of large-scale asthenospheric fabrics previously invoked. Our results support SSW-directed asthenospheric flow in western Anatolia and the Aegean toward the retreating Hellenic trench, with corresponding lithospheric deformation, while in the east of the region, they are consistent with a northward flow of asthenosphere in Arabia, both as predicted in mantle flow models (Confal et al., 2018;Faccenna et al., 2013). However, in between, we find no evidence for westward flow in central Anatolia, contradicting models that invoke toroidal asthenospheric flow driving the plate motions. Instead, we suggest asthenospheric flow in Anatolia is dominated by smaller-scale processes. High  E t in west-central Anatolia is readily associated with mantle flow between the Aegean and Cyprus slabs, and  E is diverted around the edges of the slabs. Lower values of  E t below Central Anatolia, and null results in the Cappadocia Volcanic Province and the north-central coast of Anatolia, corroborate studies that invoke vertical mantle flow fields, due to dripping of Anatolian continental lithosphere and/or flow through a developing tear in the Cyprus slab. Anisotropy in eastern Anatolia suggests either NE-or SW-directed flow. Given the lack of evidence for a toroidal flow of asthenosphere matching the Anatolian APM, we conclude plate motions in this region may be driven by other forces, such as gradients in gravitational potential energy (e.g., England et al., 2016) and/or boundary forces (e.g., Capitanio, 2016).

Data Availability Statement
Seismic data from network codes marked in Figure 2b as AB, GO, IU, MP, XG, XH, YB, YH, and YL are available from the IRIS DMC (https://ds.iris.edu/ds/nodes/dmc); CQ, HC, HL, and HT from NOA (http://eida. gein.noa.gr); GE, IS, JS, ZR, Z3, and Z4 from GEOFON (http://geofon.gfz-potsdam.de); KO and 6G from KOERI (http://eida-service.koeri.boun.edu.tr); MN from INGV (http://webservices.ingv.it); TU and YF from ORFEUS (http://www.orfeus-eu.org); XW, XY, and YI from RESIF (http://ws.resif.fr). Data from the TROODOS network (Bastow et al., 2016) in Cyprus (code XG in Figure 2b) will become available in 2022 via IRIS. All seismic data were processed using the Seismic Analysis Code (Helffrich et al., 2013). Shear wave splitting source code is available at https://github.com/jwookey/sheba (Wüstefeld et al., 2010). Two-layer modeling was performed using the MSAT toolkit (Walker & Wookey, 2012), available at https://github.com/ andreww/MSAT. Prior shear-wave splitting results displayed in Figure 2a were retrieved from http://splitting.gm.univ-montp2.fr/DB/ (Wüstefeld et al., 2009), and our results will be added to this database. gave advice and example codes for grid-search inversions using the MSAT toolkit. We thank F. McNab for providing the uplift model used in Figure 10 inset and advice on its representation. We are also grateful to C. Kreemer, W. Wei and J. Weiss for their helpful responses to data requests. We are indebted to the numerous people involved in deploying and maintaining the seismic networks used in this study, and making the data publicly available. We created figures in this paper using the Generic Mapping Tools software (Wessel et al., 2019). We thank R. Govers and an anonymous reviewer for helpful comments and advice.