Seismic Tomographic Imaging of the Eastern Mediterranean Mantle: Implications for Terminal‐Stage Subduction, the Uplift of Anatolia, and the Development of the North Anatolian Fault

The Eastern Mediterranean captures the east‐west transition from active subduction of Earth's oldest oceanic lithosphere to continental collision, making it an ideal location to study terminal‐stage subduction. Asthenospheric‐ or subduction‐related processes are the main candidates for the region's ∼2 km uplift and Miocene volcanism; however, their relative importance is debated. To address these issues, we present new P and S wave relative arrival‐time tomographic models that reveal fast anomalies associated with an intact Aegean slab in the west, progressing to a fragmented, partially continental, Cyprean slab below central Anatolia. We resolve a gap between the Aegean and Cyprean slabs, and a horizontal tear in the Cyprean slab below the Central Anatolian Volcanic Province. Below eastern Anatolia, the completely detached “Bitlis” slab is characterized by fast wave speeds at ∼500 km depth. Assuming slab sinking rates mirror Arabia‐Anatolia convergence rates, the Bitlis slab's location indicates an Oligocene (∼26 Ma) break‐off. Results further reveal a strong velocity contrast across the North Anatolian Fault likely representing a 40–60 km decrease in lithospheric thickness from the Precambrian lithosphere north of the fault to a thinned Anatolian lithosphere in the south. Slow uppermost‐mantle wave speeds below active volcanoes in eastern Anatolia, and ratios of P to S wave relative traveltimes, indicate a thin lithosphere and melt contributions. Positive central and eastern Anatolian residual topography requires additional support from hot/buoyant asthenosphere to maintain the 1–2 km elevation in addition to an almost absent lithospheric mantle. Small‐scale fast velocity structures in the shallow mantle above the Bitlis slab may therefore be drips of Anatolian lithospheric mantle.


Overview
The fate and geometry of subducting slabs during the final stages of subduction are debated. Several studies have shown that slabs can segment and tear as they descend (e.g., Portner, Delph, et al., 2018;Wortel & Spakman, 2000), particularly in continental collision settings where subduction termination is marked by detachment (e.g., Duretz & Gerya, 2013). The interplay between terminal-stage slab evolution and overriding plate deformation can be complex; for example, several postcollisional processes, such as continental lithospheric delamination, may cause magmatism and uplift. Slab segmentation and/or break-off produce gaps thought to facilitate upwelling and decompression melting of hot asthenosphere, leading to magmatism (e.g., Davies & von Blanckenburg, 1995). The Eastern Mediterranean (Figure 1) is an ideal study locale for late-stage subduction processes and terminal-stage slab morphology: It captures the transition from subduction of Earth's oldest oceanic lithosphere in the west to collision-dominated convergence in the east.
Debate surrounds the causes of magmatism and uplift in central and eastern Anatolia, and if they have a common geodynamic evolution relating to consequences of terminal-stage subduction. Some petrological and geodynamic modeling studies favor slab steepening/peel-back of initially flat subduction and eventual break-off as the cause, suggesting it is required to explain N-S discrepancies in volcanic alkalinity (Bartol & Govers, 2014;Keskin, 2003;Şengör et al., 2008) and satisfy the 1-2 km plateau uplift (Göğüş & Pysklywec, 2008;Göğüş & Ueda, 2018;Memiş et al., 2020). Other geochronologic and geodynamic studies favor lithospheric removal/thinning in the form of gravitational instabilities developed during collision (Göğüş, Pysklywec, & Faccenna, 2017;Schleiffarth et al., 2018), advocating it is required to explain the rapid surface uplift and volcanism along plateau edges. Better constraints for lithospheric thickness and the presence of thermal or melt derived anomalies below central and eastern Anatolia may help discriminate between the two scenarios.
The Neo-Tethyan suture zone, formed from repeated collisions of multiple subduction zones (e.g., Okay & Tüysüz, 1999), has reactivated as a continent-scale transform to accommodate the westward motion of Anatolia: the North Anatolian Fault Zone (NAFZ; Figure 1). The Dead Sea Transform (DST; Figure 1; e.g., Weber et al., 2009), which connects to the East Anatolian Fault Zone, has formed along other preexisting crustal weak zones to accommodate differential motion between Africa and Arabia after continental collision. Both are major active plate boundaries, yet consensus is lacking on their expression in the lithospheric mantle (e.g., Çubuk-Sabuncu et al., 2017;Fichtner et al., 2013;Koulakov et al., 2006;Papaleo et al., 2017Papaleo et al., , 2018, specifically, whether there are differences in properties across the fault that leave an imprint on seismic structure to mantle lithospheric depths. To address these outstanding tectonic and geodynamic questions, teleseismic P and S wave relative arrival-time tomography is used to probe the upper-mantle seismic structure of the Eastern Mediterranean. Seismic tomography can help constrain upper-mantle structure, temperature, and the presence of melt. At uppermost-mantle depths, we will test whether slow wave speeds are diagnostic of thinned lithosphere and/or the presence of decompression melt. In midmantle depths, the integrity of slab structure when traversing the west-to-east subduction-collision transition will be determined, revealing any variations in physical properties between the subducting slabs. The deepest parts of our models will reveal how the slabs interact with the mantle transition zone (MTZ). Our study benefits from larger data volumes than previous studies, spanning the time period 1997-2019, including data from a seismograph network that operated in Cyprus from March 2016 to March 2019 (Bastow et al., 2017), and constrains the first S wave body-wave model for the region. Our results provide insight into slab character at terminal-stage subduction zones and have important implications for regional magmatism and uplift.  Hayes et al. (2018). Black arrows show plate motion relative to Eurasia in mm year −1 (Reilinger et al., 2006).

Subduction
Northward subduction of the oldest, Neo-Tethyan, oceanic lithosphere (e.g., Granot, 2016), is ongoing beneath the Aegean and Cyprean trenches (e.g., Lefebvre, 2011). Slab rollback and southward Aegean trench migration since 30-35 Ma  exerts some of the highest extensional strain rates (>100 × 10 −9 year −1 ) in the world within western Anatolia (Kreemer et al., 2014;Reilinger et al., 2006). Both trenches are associated with Wadati-Benioff zones terminating at 200 km depth, with the Cyprean trench much less seismically active (e.g., Duman et al., 2018;Wdowinski et al., 2006). The arrival of a Gondwanan continental fragment-the Eratosthenes seamount ( Figure 1; Cowie & Kusznir, 2012;Granot, 2016)-at the Cyprean trench signifies the transition to collisional underthrusting of African continental lithosphere (e.g., Feld et al., 2017;Glover & Robertson, 1998). Until recently, it was widely accepted that the Cyprean trench connects to the Aegean trench via the Pliny-Strabo Transform Edge Propagator (STEP) fault system (Figure 1; e.g., Govers & Wortel, 2005). However, Bocchini et al. (2018) reject this model from a seismicity perspective, arguing there is a clear subducting slab dipping NW to the north of Rhodes, and instead propose a very shallow tear between the two slabs. Kaymakçı et al. (2018) observe a lack of differential rotations in NE directions, through a paleomagnetic approach, to come to a similar conclusion. On the surface, the meeting point between the two arcs is the Isparta Angle ( Figure 1a), a tectonically complex zone experiencing a mixture of extensional and transtensional regimes (e.g., Över et al., 2016). The Isparta Angle is bounded by a NE-trending left-lateral strike-slip fault that splays into multiple subparallel faults in the west and NW-SE trending faults in the east (e.g., Över et al., 2016;Yagmurlu et al., 1997). These faults bound a series of SW-and SE-trending rock units emplaced from the Late Cretaceous to Late Miocene (e.g., Bozkurt, 2001;Robertsonet al., 2003). The submarine Anaximander seamounts, bound the Isparta Angle to the south (Figure 1; e.g., Zitter et al., 2003). The Cyprean trench then joins the Bitlis-Zagros suture through the East Anatolian Fault via a diffuse boundary beneath the Adana Basin (Figure 1; e.g., Abgarmi et al., 2017).
The structure of the DST is debated, with Koulakov et al. (2006) finding no seismic signature extending into the mantle in contrast to Laske et al. (2008) who suggest, based on Rayleigh wave dispersion tomography, that the upper mantle east of the DST has been altered enough by an elevated geotherm to produce a seismic signal. In the mantle lithosphere below the NAFZ, Berk Biryol et al. (2011) constrained a 4% faster P wave velocity in the north of the fault compared to the south, whereas the full-waveform tomography models of Govers and Fichtner (2016) indicate no significant velocity contrasts across the NAFZ west of 33°E, where the fault bifurcates.

Seismic Data and Networks
Relative arrival-time residuals were determined using teleseismic recordings from 553 temporary and permanent broadband stations (see Table S8 in the supporting information for station details) in operation between 1997 and 2019 ( Figure 2a). Each is sensitive to a broad range of frequencies (0.008-25 Hz), recording data at ≥50 samples per second. 12,285 earthquakes of m b ≥ 5.5 in the teleseismic epicentral distance range 30°<Δ<100°were visually inspected for a good signal-to-noise ratio. Station-earthquake pairs with Δ ≤ 30°were excluded to avoid nonunique raypaths caused by MTZ-associated triplications. A total of 2,448 earthquakes were ultimately analyzed for direct P phases. Similarly, 1,440 earthquakes yielded direct Sand SKS phases (Figures 2b and 2c). Earthquake coverage is best from the NE (Himalayas and the western Pacific subduction zones: 30-120°). PKP, PKIKP, and PKiKP phases were omitted because multiple phase arrivals generally overlapped temporally and thus could not be attributed to a specific phase or the correct raypath.

Method of Relative Arrival-Time Determination
The first identifiable coherent peak or trough of body-wave energy across the network was manually picked on waveforms filtered using zero-phase two-pole Butterworth bandpass filters with corner frequencies of 0.4 and 2 Hz and 0.04 and 0.2 Hz for P and S waves, respectively. This allows for the inclusion of higher frequency signals, satisfying the infinite frequency approximation (ray theory) assumed in the inversion process by forming a trade-off between optimizing high-frequency signals and limiting high-frequency noise in the data set (a method also adopted by previous studies; e.g., Bastow et al., 2005;Boyce et al., 2016). Direct P and S phases were picked on the vertical and tangential component seismograms, respectively; SKS phases were picked on the radial. The manual alignment of traces was refined through multichannel cross correlation (MCCC; VanDecar & Crosson, 1990) allowing for accurate computation of relative arrival-time residuals.
The MCCC procedure cross-correlates waveforms of a given phase for a single earthquake in a time window of 3 s for P waves and 12 s for S waves, avoiding interference from secondary phase arrivals, to search for the cross-correlation maximum position, τ max ij , between trace pairs. These are subsequently optimized using the least squares approach to determine the most appropriate cross-correlation derived relative delay time, Δt ij , for each station pair (i, j) for each earthquake: where t P i and t P j are preliminary-picked-arrival estimates for the ith and jth traces, respectively. Rays with a mean cross-correlation coefficient <0.85 were rejected. Cycle skips and GPS timing errors were systematically identified and eliminated. The best fitting station delay times are then calculated by finding the least squares solution of a system of overdetermined linear equations formulated from the set of relative time shifts. The least squares minimization procedure provides more accurate relative arrival-time measurements for each station by forcing the mean to be zero.
Error associated with each arrival in the MCCC method is expressed as the standard deviation (σ) of distributed arrivals associated with that trace: where i and j represent pairs of the n stations. Using this method, relative arrival times for the P and S wave data set have a σ of order 0.01 and 0.08 s, respectively. In line with comparable studies (e.g., Bastow et al., 2005;Boyce et al., 2016), these errors were considered optimistic, especially given the trace sampling rate of ≥0.01 s.
Relative arrival-time residuals, t rel , were then computed using The back azimuthal and distance distribution of earthquakes used in this study between 1997 and 2019. (b) m b ≥ 5.5 P wave event epicenters (circles). (c) m b ≥ 5.5 S wave (circles) and SKS (crosses and red circles) event epicenters. The concentric circles mark 30°intervals in epicentral distance from the center of the network at 32°N, 32°E (black star). Note that although some earthquakes are located within 30°from the study area on the figure, no arrivals within 30°from individual stations have been included in the analysis to avoid triplications associated with the mantle transition zone.
where t i is the phase pick time for each station i, t ei is the expected IASP91 (Kennett & Engdahl, 1991) traveltime for station i, and t e is the mean IASP91 expected traveltime for the specific earthquake across the network. The final data sets comprise 69,108 P wave relative arrival times, significantly more than any previous tomographic studies of the region (e.g., Berk Biryol et al., 2011;Lei & Zhao, 2007;Portner, Delph, et al., 2018;Zor, 2008), and 47,419 S wave relative arrival times.

Analysis of Relative Arrival-Time Residuals
MCCC is a relative arrival-time method that calculates anomalies with respect to the background mean of the data set by completely removing it from the resulting models. Any fast and slow anomalies are thus relative to the regional average arrival time (e.g., Bastow, 2012). The network used covers regions of both globally fast-along the Aegean subduction zone-and slow-in eastern Anatolia-anomalies when compared to various global tomographic models (e.g., Ritsema et al., 2011;Simmons et al., 2012).
Both P and S wave mean relative arrival times (Figures 3a and 3b) reveal a trend of earlier-than-expected arrivals along the Aegean subduction zone (δt p ≈ −1 s; δt s ≈ −3 s) and later arrivals in the EAP (δt p ≈ 1 s;

10.1029/2020GC009009
Geochemistry, Geophysics, Geosystems and NW likely reflect fast, cold subducting lithosphere, with later arrivals from the SE (Figure 4). In eastern Anatolia (e.g., station PTK), most arrivals are slow. Similar patterns are obvious along the DST (e.g., station KZIT), despite lacking ray coverage between 135°and 270°azimuths. Generally, back azimuthal and epicentral distance coverage for the whole data set is good (Figures 2b, 2c, 3d, and 3f). Station ANTB, in the center of the Isparta Angle, shows the largest peak-to-peak variation of δt p =−0.51± 0.98 s (Figure 4), with a trend more pronounced in the S wave data set, highlighting S wave sensitivity to partial melt, perhaps due to hot asthenospheric upwelling through a slab window (e.g., Berk Biryol et al., 2011). Station DB08, along the NAFZ, has the least variation (δt p =−0.16± 0.14 s).

Model Parameterization, Inversion Procedure, and A Priori Crustal Model 2.4.1. A Priori Crustal Model
Steep incidence angles and the lack of crossing raypaths at the shallowest depths mean teleseismic tomography is largely incapable of resolving crustal structure (e.g., Waldhauser et al., 2002). Given the thick (<9 km) sedimentary basins present throughout Anatolia (e.g., 5 km-thick Tuz Gölü basin in Central Anatolia; Cemen et al., 1999) and the DST (e.g., Garfunkel & Ben-Avraham, 1996), accounting for delays through low-speed sediment is also important before inverting for upper-mantle velocity perturbations. As well as 3-D crustal velocity variations, laterally varying Moho topography can enhance the impact on delay times.
Corrections assumed straight raypaths, taking into account the slowness, and constant P and S wave velocities through sediment, crust, and uppermost mantle above the station with the deepest Moho (see Text S1 for calculation and Table S1 for velocities used). Sediment thickness beneath each station was extracted from global model CRUST 1.0 (Laske et al., 2013). Crustal thicknesses produced using H-κ stacking data from Vanacore et al. (2013), Abgarmi et al. (2017), and Mohsen et al. (2005) enabled Moho corrections for all stations. For any stations lacking H-κ stacking data, crustal thicknesses were extrapolated based on the aforementioned models.
The contribution of the crustal model to the relative traveltime residuals ranges from −0.50 to 0.75 s for P waves and −1.0 to 1.5 s for S waves. The corrected residuals will more accurately reflect the proportion of delay times mapped into parts of the model where interpretations will be made. After inversion, the crustal model has had a noticeable effect at <200 km depth, beyond which differences in the model become indistinguishable. To explore the effect of errors in our a priori crustal model, we have created both checkerboard and synthetic tests introducing ± 1 and ± 3 km Moho depth uncertainties into the models (see Figures S8  and S9). We conclude that deterioration of the model resolution caused by crustal uncertainties at 100 km depth is minimal and ≥300 km is negligible.

Tomographic Inversion Procedure
The regularized, least squares inversion method of VanDecar et al. (1995) was used to invert the relative arrival-time residuals for mantle velocity perturbations. An equal number of conjugate gradient model iterations were used for both the P and S wave models. Model slowness is parameterized, identically for both P and S wave models, using B-splines under tension (Cline, 1981) over a dense grid of 101,400 knots in total, covering 26 knots between 0 and 1,400 km in depth, 75 knots between 14°E and 57°E in longitude, and 52 knots between 20°N and 51°N in latitude. The best resolved, central portion of the study area is sampled at 30 km depth intervals, restricted to above 800 km depth and 0.3°in latitude/longitude (Table S2). Beyond this region coarser spacing is used, creating zones where unwarranted or spurious structure can map without being incorporated into the interpreted parts of the model. Results outside the central part of the study area are thus not considered robust enough for interpretation.
To address the underdetermined nature of the tomographic inversion, with fewer rays than model unknowns, and to prevent overinterpretation of finer-scale structures, regularization is applied to suppress primary and secondary spatial gradients. These correspond to flattening and smoothing parameters, respectively. A preferred model was chosen by investigating the trade-off between RMS (root mean square) residual reduction (data fit) and model roughness (model complexity) controlled by the regularization parameters. Concurrently, we avoid models achieving greater residual reduction than justifiable from MCCC-defined data noise levels, which we consider unrealistically low, in line with previous studies (e.g., Tilmann et al., 2001). An Occam's razor approach was adopted to determine the simplest model with the least structure required to fit the data. The preferred models, chosen empirically, are located near the knee of the trade-off curves ( Figures S4 and S5) and fit 88.55% and 82.86% of the RMS estimated relativedelay times for the P and S wave models, respectively.

Resolution Testing
To assess linear resolving power of the P and S wave inversions, a standard checkerboard analysis was conducted ( Figure 5), using alternating δV p = δV s = ± 5%, 50 km diameter spherical anomalies, defined by Gaussian functions across their diameter (see Figure S6 for the P wave checkerboard test). They were placed at 200 km depth intervals from 100 to 700 km depth, and every 2°in latitude/longitude. The structures were inverted using identical model parameterization and regularization to the original data set inversion. Gaussian residual time error components with standard deviations of 0.03 and 0.1 s were added to theoretical traveltimes for P and S waves, respectively.
Amplitude recovery is >70% at depths >300 km although anomalies experience increasing lateral smearing with depth ( Figure 5). Shallower layers recover ∼50% of amplitudes with very high lateral resolution since ray density is highest at <300 km depth. Although vertical smearing is minimal and largely observed at boundaries, structure in the lower 600-800 km of the model does not smear into shallower regions.
Checkerboard tests are poor for identifying possible model artifacts, typically representing extreme limits of model resolution. More realistic synthetic tests are thus used to assess the ability of the data set to retrieve tectonically/geodynamically plausible structures. Several structural models were created ( Figure 6) to address outstanding debates, including variable slab dip and segmentation (Portner, Delph, et al., 2018) Geochemistry, Geophysics, Geosystems terrane boundaries (Fichtner et al., 2013). See Table S5 for synthetic test parameters and Figure S7 for further synthetic tests.
Gaussian residual time errors were added and the data inverted as with the checkerboard tests. The Pwave model recovers >70% of anomaly amplitudes (Figure 6), increasing with depth in central Anatolian regions and decreasing around the Aegean. Long-wavelength structures are easier to retrieve with smoothing parameterization than steep gradients, hence the greater amplitude recovery compared to checkerboard tests (Figures 5). There is minimal lateral smearing, and vertical boundaries are well constrained with some vertical smearing occurring at model edges; the Aegean slab extends for an extra ∼100 km depth, with the Cyprean slab smeared to a lesser extent due to a higher density of crossing ray paths (Figure 6; A-A′). The S wave model also achieves ∼70% amplitude recovery (Figure 6), experiencing marginally greater smearing than the P wave model. Figures S8 and S9 show the effect of crustal thickness errors in our a priori crustal model on several of the resolution tests in Figures 5 and 6. Ray hit count maps also indicate the best resolved parts of the model (Figures S10 and S11).

Tomographic Results
The velocity perturbations in our relativearrival-time derived P and S wave models (Figures 7 and 8) are relative to the regional mean (Bastow, 2012). Based on global tomographic model LLNL-G3Dv3 (Simmons et al., 2012), a mean fast wave speed anomaly of δV p =0.42% exists below 38.59°N, 34.99°E in the upper 700 km of the model, where Station AT22 (Figure 1) has near-zero relative arrival-time residual in the P wave data set. Similarly for S waves, in the SL2013sv model (Schaeffer & Lebedev, 2013) a global mean δV s =0.07% exists below station AT06 (39.67°N, 35.66°E) in the upper 700 km of the model. Thus, the P wave model could be considered to be on a pedestal that is 0.42% slower than the global average; the S wave model may be close to the global mean. The wide-aperture network used therefore covers regions of both genuinely fast and slow wave speeds when benchmarked against global tomographic studies. Subducting slabs in the Eastern Mediterranean are faster than the global mean; areas such as the EAP are slower.

Geochemistry, Geophysics, Geosystems
As well as deviations from the global zero mean, further uncertainties in amplitude estimation stem from the underdetermined nature of the smoothed inversions. Resolution testing implies amplitude recovery varies throughout the model from ∼50-70%, with long-wavelength anomalies recovered best (Figures 5 and 6). Uncertainties in δV p,s are also revealed during trade-off curve analysis (Figures S4 and S5) when regularization parameters are chosen.

Causes of Seismic Heterogeneity
When interpreting tomographic models, it is important to note that the δV p,s = 0% contour in relative traveltime tomography is not necessarily the same as the global zero mean (e.g., Bastow, 2012). The δV p,s = 0% contour in Figures 7 and 8 differs from the global mean by δV p = 0.42% (Simmons et al., 2012) and δV s = 0.07% (Schaeffer & Lebedev, 2013), respectively. However, the relative high and low velocity anomalies we present, therefore, genuinely span fast and slow mantle structure. Various factors can affect upper-mantle wave speeds, but temperature exerts greater influence than composition: for compositions ranging from refractory dunite to harzburgite and fertile peridotite, estimated velocity variations are usually <1% (e.g., Cammarano et al., 2003;Goes et al., 2000). We thus begin by exploring the mantle thermal structure that would be expected if we assume our seismic velocity anomalies are purely due to temperature variations.

Thermal Structure From Seismic Anomalies
To obtain absolute wave speeds for a thermodynamic conversion to temperature, we express our S wave velocity anomalies with respect to IASP91 (Kennett & Engdahl, 1991). We assume our relative S wave anomalies are true anomalies since our δV s = 0% contour is close to the global zero mean. The conversion was made following methods described in Cobden et al. (2008) and was based on phase compositions and elastic parameters computed with PerPleX (Connolly, 2005) using the database from Stixrude and Lithgow-Bertelloni (2011) and the anelasticity model Q7g (Styles et al., 2011). Western Anatolia (31°E, 38.5°N) was set as a reference point assumed to fall along the 1,300°C adiabat for a pyrolitic composition (Styles et al., 2011) below 100 km depth. Although basaltic oceanic crust has a compositionally distinct seismic signature, it would be too thin to resolve in our models, so for an estimate of thermal anomalies, we consider subducting slabs as peridotitic in composition. Uncertainties in temperature estimates inherent in the conversion method are a few tens of degrees (Cammarano et al., 2003;Goes et al., 2000). Additional uncertainty arises from the dependence of anomaly amplitudes on resolution and inversion regularization parameter choices (Figures S4 and S5). An RMS uncertainty of 0.5% in anomaly amplitude, when moving down the trade-off curve by 2% in RMS residual reduction, produces an RMS uncertainty of ± 73°C in 10.1029/2020GC009009 Geochemistry, Geophysics, Geosystems inferred temperature across the study area at 100 km depth. Near the base of the lithosphere, a further uncertainty arises from the limitations on depth resolution, so that variations in lithospheric thickness will map into variations in velocity and hence by conversion into variations in temperatures.
The highest temperatures, 1,600-1700°C at 100 km depth, are located in eastern Anatolia (Figure 9b), strongly indicative of hot asthenosphere/thin lithosphere. Temperatures of 1,600-1700°C exceed both the dry solidus (e.g., Rey, 2015) and petrological estimates of mantle potential temperature below the region (1,400°C; e.g., McNab et al., 2018), implying melt is likely present in the eastern Anatolian upper mantle (Figures 9b and 9c). Even considering uncertainties in these temperature estimates, temperatures in eastern Anatolia would still exceed both adiabatic mantle (by 100-400°C; Figure 9c) and dry mantle solidus temperatures (by 100-200°C; Figure 9c). With melt present, our temperatures would be overestimates. Assuming the actual mantle is as hot as inferred from petrology, and the rest of the seismic anomaly is due to melt, small fractions of melt (depending on geometry) could suffice to explain the additionally slow seismic velocities. In central Anatolia, temperature estimates of (1,400-1,500°C) only slightly exceed petrological estimates, and little to no melt would be required. Lithosphere thinner than 50-70 km would facilitate decompressional melting, a point further discussed in section 4.3.

Geochemistry, Geophysics, Geosystems
Assessing traveltimes rather than velocity perturbations circumnavigates amplitude recovery issues associated with smoothed, underdetermined tomographic velocity models. Several studies have cited the ratio of S and P wave traveltimes ( δt s δt p ) for common earthquake-station pairs as a diagnostic tool for causes of seismic anomalies (e.g., Bastow et al., 2005;Bolton & Masters, 2001;Civiero et al., 2016). A δt s δt p of 3.6-3.8 defines the upper bounds of the thermal range (e.g., Cammarano et al., 2003). A direct comparison of residual subsets in central and eastern Anatolia indicates a generally higher ratio in the east ( δt s δt p ∼ 4:01; Figure 9a), suggestive of the presence of mantle melt as a cause for the seismic anomalies. In central Anatolia, the ratio lies within the thermal regime ( δt s δt p ∼ 3:54; Figure 9a). Bastow et al. (2005) attribute a δt s δt p of between ∼5 and ∼10 to significant fractions of mantle melt beneath the Main Ethiopian Rift; larger ratios than eastern Anatolia, implying that if mantle melt is present, it is likely less voluminous. Since traveltimes are path-integrated values (i.e., sample 3-D heterogeneity along the entire ray path), and Swaves are generally longer wavelength than P waves, thus sampling a greater portion of the mantle, these ratios are likely to be underestimates of local ratios (e.g., Schmandt & Humphreys, 2010).
Differences are evident in the seismic and inferred thermal structure of the Aegean (δV p =3%; Figures 7 and 8; ∼800°C; Figure 9) and Cyprean (δV p =1.5%; Figures 7 and 8; ∼1,000°C; Figure 9) slabs. A simple explanation for this is the slower subduction rate below Cyprus (present-day <10 mm year −1 ; Reilinger et al., 2006) compared to the Aegean (30 mm year −1 ; Reilinger et al., 2006), which will render the Cyprus slab warmer and seismically slower in the uppermost mantle, as is observed (Figure 9). This, as well as differences in slabcharacter, where the Cyprean slab is perhaps partially continental further east (Figure 1), may add to the along-strike velocity differences that characterize slab morphology; we revisit this hypothesis in section 4.2.

Slab Versus Drip Anomalies
Diffuse fast wave speed anomalies below Tibet have been proposed to correspond to large blocks of foundered lithosphere (Chen et al., 2017), while small-scale (tens to hundreds of kilometers wide) fast anomalies below the Western United States have been attributed to lithospheric drips (e.g., van Wijk et al., 2008). Short length-scale (20-40 km wide) fast wave speed anomalies beneath central and eastern Anatolia (Figures 7 and  8; F-F′), previously interpreted as Bitlis slab fragments may, in fact, be lithospheric drips given the scale and relatively low amplitudes of the anomalies (e.g., van Wijk et al., 2008). Continental lithospheric drips are likely to be imaged as lower amplitude anomalies than oceanic slab fragments since their temperatures are at most a few hundred degrees lower than the surrounding mantle (van Wijk et al., 2008); oceanic slab likely has larger deviations, and thus faster wave speeds. These observations, in conjunction with a thinned or absent lithospheric mantle below Anatolia, motivate our interpretation that these anomalies are lithospheric drips. Synthetic tests ( Figure 6) show that high-amplitude (>3%), >30 km wide structures can be recovered from the tomographic inversions; however, we cannot completely preclude the possibility that they are model artifacts.

Contrasts Across Transform Faults
Observed peak-to-peak variations of δV p ≈ 5% and δV s ≈ 7%, corresponding to temperature contrasts of ∼700°C (Figure 9c), exist at 100 km depth across the western segment of the NAFZ (Figures 7 and  8; B-B′, C-C′), too large to be explained by lateral contrasts in lithospheric composition. The sharp, near-vertical, boundary is thus likely the result of juxtaposed asthenospheric material in the south and possibly refractory (Artemieva & Mooney, 2001), relatively cold, Precambrian lithosphere in the north (the Istanbul Zone; e.g., Okay, 2008). Such a big step in lithospheric thicknesses would require recent thinning of Anatolian lithosphere south of the East European Craton (∼200 km-thick Istanbul Zone; Figure 1b; Schaeffer & Lebedev, 2013). Associated with slow wave speeds, and in agreement with thinner lithosphere in the south, is the Miocene-Quaternary Galatia Volcanic Province (GVP; Figure 1), where geochemical studies suggest volcanism is primarily the result of asthenospheric decompression melting (e.g., Varol et al., 2014;Wilson et al., 1997). Observations of crustal structure and isostatic calculations (Frederiksen et al., 2015;Jenkins et al., 2020) corroborate the suggestion of thicker lithosphere beneath the Istanbul Zone compared to Sakarya south of the NAFZ (Figure 1). Subduction-related hydration 10.1029/2020GC009009

Geochemistry, Geophysics, Geosystems
and magmatic processes may have weakened the overriding Anatolian lithosphere (Arcay et al., 2006) inducing lithospheric thinning via delamination south of the more refractory East European craton, perhaps akin to the Basin and Range (e.g., Humphreys et al., 2003). Alternatively, the NAFZ may mark the northern limit of metasomatic modification, regardless of initial contrasts between Precambrian and Phanerozoic lithospheric structure.
Unlike the NAFZ, the DST lacks significant wave speed and temperature contrasts between the African and Arabian plates (Figures 7-9b, and 9c). Instead, low wave speeds (δV p ≈ δV s ≈ −1.5%) and high estimated temperatures (1,300-1,500°C) characterize shallow depths (∼100 km) along strike of the DST, which could be consistent with asthenospheric material flowing from the Afar hot spot along DST base-of-the-lithosphere topography (e.g., Ebinger & Sleep, 1998;Faccenna et al., 2013). Alternatively, Chang and van der Lee (2011) have suggested the existence of a mantle plume below Jordan based on their own S wave tomographic model of Arabia and East Africa. A shallow (<68 km deep) LAB has been constrained by receiver function studies (Mohsen et al., 2005), consistent with our observed low wave speeds at shallow depths (Figures 7 and 8).

Slab Morphology Beneath the Eastern Mediterranean 4.2.1. Shallow Slab Regions
A striking observation in Figures 7 and 8 at 100 km depth is the fast wave speed material at the Cyprean and Aegean trenches, readily interpreted as subducting African lithosphere. A slow δV p,s slab window exists separating the two zones below the Isparta Angle, implying the presence of a vertical tear; consistent with previous studies (Berk Biryol et al., 2011;Paul et al., 2014;Piromallo & Morelli, 2003;Portner, Delph, et al., 2018;Spakman et al., 1993). The discontinuity parallels the Strabo and Pliny Transforms, termed a Subduction Transform Edge Propagator (STEP; Govers & Wortel, 2005;Özbakır et al., 2013). Berk Biryol et al. (2011) suggest tearing nucleated at 300 km depth where the Aegean and Cyprean slab anomalies are furthest apart in their P wave models. In Figures 7 and 8, however, the window is approximately equidistant 200-410 km depth. The existence of the Aegean-Cyprean slab window (Figures 7 and 8; 100-500 km depth slices) is also consistent with shear-wave splitting studies (e.g., Evangelidis, 2017;Paul et al., 2014;Wei et al., 2019), which interpret strong azimuthal mantle anisotropy patterns as flow through the slab gap. If the flow includes a component of upwelling, it may explain Neogene magmatism within the Isparta Angle (Figure 1), which is dominated by anorogenic potassic or ultrapotassic alkaline volcanic rocks (e.g., Dilek & Altunkaynak, 2010;Karaoğlu & Helvacı, 2014;Yagmurlu et al., 1997): uncharacteristic subduction zone lithologies.
Pliny-Strabo STEP tear development hypotheses generally involve continental lithosphere: oceanic versus continental subduction between the Aegean and Cyprean slabs (Wortel & Spakman, 2000), the arrival of continental fragments at the trench in the Late Eocene (the collision of the continental Central Anatolian Crystalline Complex with the Pontides; Figure 1b; Govers & Fichtner, 2016), or the present-day arrival of the attenuated-continental Anaximander seamount (Figure 1a; Delph et al., 2015;Zitter et al., 2003). Although Granot (2016) place the current continental-oceanic transition further east below Cyprus (Figure 1), the lower temperatures (Figures 9b and 9c) and wave speeds (Figures 7 and 8) that characterize the Cyprean slab may reflect buoyancy contrasts that drive variable sinking rates and thus the observed slab tearing.
At ∼38°N, slow wave speeds (δV p <−1.5%) in the eastern portion of the Cyprean slab are suggestive of a horizontal tear at 200 km depth (Figure 7; C-C′). Although the velocity contrast is subtle, synthetic tests ( Figure 6) indicate the separation is resolvable. Buoyancy contrasts within the Cyprean slab, and the plate geometry that makes subduction increasingly restricted east of the Aegean trench, likely contribute to the difference in convergence rate along the Cyprean and Aegean trenches (<10 and ∼30 mm year −1 , respectively; Reilinger et al., 2006). By tearing, both vertically below the Isparta Angle, and horizontally in the eastern portion of the Cyprean slab, these differential rates can be accommodated within the subducting plate. Further east, transitioning from active subduction to continental collision, fast wave speed anomalies are confined to the transition zone, suggesting the Bitlis slab has completely broken off (e.g., Berk Biryol et al., 2011;Keskin, 2003;Lei & Zhao, 2007;Portner, Delph, et al., 2018).

Intermediate Depth Slab Geometry
Eastern Mediterranean seismicity terminates at 200 km depth (Papazachos et al., 2000), so we follow Portner and Hayes (2018) to create an upper mantle slab geometry model ( Figure 10) from our P wave tomography model by assuming the center of the fast anomaly reflects the dipping slab. Figure 10 needs to be treated cautiously below Cyprus. The non-slab-like nature of the Cyprean fast wave speed anomaly (Figures 7 and 8; C-C′) could be a consequence of increasing slab buckling resulting from slower African plate sinking rates, and perhaps added drips of removed, Anatolian, continental lithosphere sitting atop the slab surface (we revisit this concept in section 4.3). The larger thickness of the Cyprean slab could also reflect its complex geometry around a tightly curved trench (Figure 1), where the slab is resolved tomographically as a broad feature. Figure 10 reveals an intact Aegean slab to 600 km depth. In contrast, the Cyprean slab is increasingly structurally complex eastward, traversing the transition from subduction to continental collision. A kink in the slab isodepth contours NW of Cyprus aligns with the dextral Paphos Transform Fault (Figure 10a). Where oceanic lithosphere is currently impinging on this part of the trench, west of the fault, the Cyprean slab exhibits simpler geometry, as previously noted by Portner, Delph, et al. (2018). Line C-C′ in Figure 10a corresponds to cross-section C-C′ in Figure 7, highlighting the correlation between slab tearing and the CAVP at ∼38°N.

Slabs in the MTZ: Implications for Slab Sinking Rates and Bitlis Break-off Age
The Aegean and Cyprean slabs traverse the 410 km discontinuity uninhibited. However, both flatten (dips of 20°and 40°, respectively) and/or thicken as they penetrate the MTZ (Figures 7 and 8; B-B′, C-C′). Geodynamic modeling (e.g., Karato et al., 2001) coupled with seismic tomography constraints indicate viscosity contrasts at the 660 km discontinuity can alter slab geometry (e.g., Fukao & Obayashi, 2013;Goes et al., 2017;van der Hilst et al., 1991). This is difficult to test on the Cyprean slab given the non-slablike nature of the fast wave speed anomaly in the tomographic images (Figures 7 and 8; C-C′). The Aegean and Cyprean slabs sink independently below 660 km depth (Figure 7; A-A′). This challenges previous tomographic studies (Berk Biryol et al., 2011;Bijwaard et al., 1998) that suggest the slabs merge and pond on the 660 km discontinuity but agrees with Portner, Delph, et al. (2018) who extended the work of Berk Biryol et al. (2011) with a deeper model parameterization. Also, in agreement of slab segmentation in the MTZ are receiver function studies (Kaviani et al., 2018), imaging varying MTZ thicknesses from west to east Anatolia.
Below the eastern Anatolian collision zone, fast wave speed anomalies appear within the MTZ (Figures 7  and 8; F-F′). These are located directly below the Bitlis-Zagros suture zone. The simplest explanation is that they belong to the Bitlis slab rather than the Pontides slab, which has been proposed by some geochronological studies to have broken off in Paleocene times (72-52 Ma; e.g., Parlak et al., 2013;Schleiffarth et al., 2018) and is thus likely to be located deeper in the mantle compared to the Bitlis slab. Alternatively, it has been suggested that the Bitlis slab broke off first and Pontides slab subduction persisted until the Mid-Miocene (e.g., Keskin, 2003;Memiş et al., 2020) before peeling off beneath eastern Anatolia. We cannot preclude the possibility that some of the shallow anomalies are Pontides slab, but it would require it to have been fragmented during peel-back. Below eastern Anatolia, the fast wave speed anomalies, interpreted as the Neo-Tethyan Bitlis slab, appear within the MTZ at progressively deeper depths moving eastward, at ∼400 km (Figures 7 and 8; D-D′) and ∼500 km depth (Figures 7 and 8; F-F′) at 40°E and 44°E, respectively. This supports the hypothesis of gradual detachment as a tear beginning below Iran and propagating west (e.g., Hafkenscheid et al., 2006;Omrani et al., 2008;Rabayrol et al., 2019). Southeastern Anatolian basalts, younging toward the west (Reid et al., 2019), are consistent with this interpretation.
Observations of detached slabs in the mantle allow us to explore slab sinking rates. In settings where mantle viscosity governs both density-driven sinking velocity and lateral plate motion, vertical sinking rates are found to be comparable to plate convergence rates (e.g., Butterworth et al., 2014;Schellart & Spakman, 2012). Globally, upper mantle sinking rates of ∼10-50 mm year −1 are estimated (e.g., Billen, 2010;Goes et al., 2017;Hafkenscheid et al., 2006), with the upper values reflective of present-day Pacific sinking rates (e.g., Goes et al., 2011). Sinking decelerates at the MTZ and rates are estimated to be ∼10 mm year −1 in the lower mantle (Butterworth et al., 2014;Hafkenscheid et al., 2006). In the case of Anatolia, Arabia-Eurasia convergence rates have roughly remained constant at 20 ± 5 mm year −1 in the last 56 Ma Jolivet & Faccenna, 2000). The ∼500 km-deep Bitlis slab should thus have broken off in Oligocene times (26± 6 Ma).
This break-off age may lie on the earlier end of estimates, particularly since fast anomalies in tomographic models tend to smear downward ( Figure 6). If Arabia-Eurasia convergence rates, and therefore sinking rates, were slightly faster at 20-30 mm year −1 , the timing of slab break-off (25-16 Ma) would overlap with the exhumation of the Bitlis Massif (18-13 Ma) as determined from Apatite fission-track dating (Okay et al., 2010). A 10-25 Ma delay between the start of collision and break-off is often predicted (van de Zedde & Wortel, 2001;van Hunen & Allen, 2011). Similarly, a delay of a few Ma may exist between break-off and resulting surface uplift since the slab originally exerts a downward pull on the overlying plate (e.g., Buiter et al., 2002). Given that the first signs of collision between Arabia and Eurasia occurred 35 Ma (e.g., Allen & Armstrong, 2008), break-off would have occurred ≤35 Ma, consistent with our predictions for break-off. These values differ from the 10 Ma estimate based on similar calculations for upper-mantle fast wave speed anomalies imaged in the Anatolian tomographic model of Portner, Delph, et al. (2018), which they interpret as fragments of the Bitlis slab. Our estimates differ because we interpret these shallower small fast anomalies not as fragments of the Bitlis slab, but as delaminated Anatolian lithosphere (see section 4.1.2). Thus, our interpretation of a deeper Bitlis slab leads to an earlier estimate of slab break-off.
Studies of Mid-Miocene to present-day lavas have proposed a 10 Ma slab break-off age for the Bitlis slab (e.g., Keskin, 2003;Şengör et al., 2003). However, these geochronologically defined ages, which assume 10.1029/2020GC009009 Geochemistry, Geophysics, Geosystems zero time lag between break-off and volcanism, imply a Pacific-style ∼50 mm year −1 sinking rate, which far exceeds Arabia-Eurasia convergence. These rates, therefore, are likely too fast, especially given that the slab is likely to have been stalled by the MTZ. The 10 Ma magmatism is thus probably too late for detachment of a slab currently residing at 500 km depth (Figures 7 and 8) but may instead be a consequence of break-off on the westernmost end of Bitlis-Zagros Suture, where the slab resides at shallower depths (Figures 7 and 8), or of alternative magma-generating mechanisms. Agreeing with a later break-off in the west are geochemical signatures of Quaternary volcanic rocks in the Karlıova triple junction (KTJ; Figure 1), located along the westernmost Bitlis-Zagros Suture, which are primarily subduction-related, unlike further east, where volcanism is increasingly asthenospherically derived (Karaoğlu et al., 2020). We next explore lithospheric removal as an alternative mechanism for magmatism across Anatolia, investigating whether some of the fast wave speed mantle anomalies are removed Anatolian lithosphere, rather than subducted African oceanic lithosphere.

Implications for Causes of Magmatism in Central and Eastern Anatolia
The eastern termination of the Cyprean slab is marked by extensive and continuous low wave speed anomalies (δV p ≈−3%; δV s ≈−5%; Figures 7 and 8) extending to 300 km depth beneath the EAP and 200 km beneath the CAP. To first order, our tomographic images do not support a common evolution between central and eastern Anatolia. Eastern Anatolia likely has thinner lithosphere, evidenced by low velocities (Figures 7 and 8) and high temperatures (Figures 9) to shallow depths, with a completely broken off slab; in central Anatolia, slab material still exists above the MTZ. Surface wave studies agree with the presence of thinned lithosphere, particularly beneath the EAP where the mantle lithospheric component has been characterized as absent (e.g., Gök et al., 2007;Şengör et al., 2003). Similar geophysical studies in Central Anatolia (e.g., Vinnik et al., 2014) find a 20-30 km-thick lithospheric mantle underlying a 35 km-thick crust, consistent with the 60 km-deep lithosphere-asthenosphere boundary determined from studies of recent lavas (Reid et al., 2019). In addition to differences in lithospheric structure, central and eastern Anatolia also varies in terms of geology and volcanism, with eastern Anatolia often described as an accretionary complex with much more extensive basaltic volcanism (e.g., Keskin, 2003). Several mechanisms have been cited as responsible for volcanism in central and eastern Anatolia: slab steepening of initially flat Pontides slab subduction (Keskin, 2003), Bitlis slab break-off, and lithospheric removal in the form of dripping of continental, Anatolian, lower lithosphere (e.g., Göğüş, Pysklywec, & Faccenna, 2017). Each process is expected to generate a different surface expression, petrological signature, and mantle wave speed structure, each of which we explore here (McNab et al., 2018;Reid et al., 2017;Schleiffarth et al., 2018).

Slab Break-off and Tearing:
Postcollisional magmatism, often interpreted as the surface manifestation of break-off, has been attributed to asthenospheric upwelling and associated decompression melting in response to slab steepening and detaching from the overriding lithosphere (Davies & von Blanckenburg, 1995;van de Zedde & Wortel, 2001). Models show that strong temperature increases leading to melting occur with a delay of a few Myr, and magmatism is expected if break-off occurs at depths within the lithosphere (van de Zedde & Wortel, 2001). In this light, EAVP magmatism of varying ages has been used to infer the timing of Bitlis slab break-off (e.g., Kocaarslan & Ersoy, 2018;Rabayrol et al., 2019). However, the 10 Ma age cited by some studies (e.g., Şengör et al., 2008) is probably too late, as argued in section 4.2.3. Although some studies have proposed break-off may have progressed amagmatically (Niu, 2017;Schleiffarth et al., 2018), the petrological and geochemical characteristics of Early-Miocene calc-alkaline lavas in central-to-east Anatolia (Dargahi et al., 2010;Ghasemi & Talbot, 2006;Kocaarslan & Ersoy, 2018;Rabayrol et al., 2019) would be consistent with a 26± 6 Ma break-off age, as estimated using slab sinking rates that mirror convergence rates. Late-Miocene (10 Ma) EAVP magmatism may instead be related to a later lithospheric removal event (e.g., Rabayrol et al., 2019).
The upper mantle below the CAVP is generally characterized by low wave speeds (δV p =−2%; δV s =−4%), similar to the EAVP. The striking NE-SW linearity of the volcanic centers (Figures 1 and 10) and the calc-alkaline affinity of older magmas (12-3 Ma) suggest it was formed as a subduction-related volcanic arc (Govers & Fichtner, 2016;Schleiffarth et al., 2018). We image a tear in the Cyprean slab below the CAVP, and the slab is likely now too deep (∼200 km; Figures 7 and 8; C-C′) to feed dehydration-related volcanism. Asthenospheric melt has been proposed to explain increasing alkalinity in the most recent 10.1029/2020GC009009 Geochemistry, Geophysics, Geosystems (Quaternary; 3 Ma-12 ka) CAVP volcanism (Aydin, 2008;Reid et al., 2017). Our slab tear may provide the necessary route to the surface for this melt (e.g., Uslular & Gençalioğlu-Kuşcu, 2019). The onset of more alkaline volcanism can therefore be interpreted as a response to the initiation of slab tearing. Unless the tear started shallow enough (<75 km; Rey, 2015), adiabatic upwelling of asthenosphere alone would be insufficient to cause melting: The mantle would need to have elevated mantle potential temperatures to melt. Supporting this hypothesis, Nikogosian et al. (2018) suggest Neogene Anatolian intraplate lavas (Figure 1) have a prominent mantle plume component that results from larger-scale mantle flow affecting the wider Arabian-North African region, impacting Anatolia through gaps and tears in the subducting African lithosphere.
Interestingly, there is little present-day volcanism between the CAVP and EAVP (volcanoes in Figure 1; e.g., Schleiffarth et al., 2018;Yılmaz et al., 1998), despite the area being underlain by slow wave speeds and solidus-exceeding estimated temperatures (∼1,700°C) to shallow depths (100 km). Between 36°N and 40°N , significant volcanism last occurred during the Pliocene (10-20 Ma) and has been linked to the western propagation of Bitlis slab tearing (Rabayrol et al., 2019). This would be a further argument for two distinct causes for magmatism: break-off related magmatism (possibly aided by influx of hotter mantle) starting below the EAVP around 20 Ma, and is currently active below the CAVP, and subsequent EAVP volcanism associated with lithospheric removal. This is consistent with Reid et al. (2019), who suggest Miocene-aged volcanism in central and eastern Anatolia can be attributed to two distinct causes: the westward propagation of slab tearing in central Anatolia and lithospheric removal by dripping in eastern Anatolia.

Delamination of Anatolian Lithosphere
Lithosphere delamination is often cited as a mechanism to cause tectonic uplift and magmatism in postcollisional settings. It can take the form of lithospheric drips that remove the lower lithosphere (Ducea et al., 2013;Göğüş, Pysklywec, & Faccenna, 2017;Houseman et al., 1981), potentially aided by the production of melt (Hernlund et al., 2008), or removal of the lower lithosphere as a large coherent block along a detachment surface (Göğüş, Pysklywec, Şengör, et al., 2017;Memiş et al., 2020). The Anatolian plate is a good candidate for the first process (e.g., Göğüş, Pysklywec, & Faccenna, 2017;Schleiffarth et al., 2018) because its lithosphere has been weakened by extensive hydration and magmatism during protracted African plate subduction. But it has also been suggested that progressive removal/peel-back of the Pontides slab led to EAP lithospheric thinning, magmatism, and uplift (e.g., Memiş et al., 2020). Low wave speeds and elevated temperatures ( Figure 9) at shallow depths (<100 km) below the CAVP and EAVP are consistent with thinned lithosphere. Contrary to past consensus, Topuz et al. (2017) suggest eastern Anatolia is a continental block, not an accretionary complex, implying that lithospheric thinning there may thus not be due to oceanic slab steepening and delamination beneath an accretionary complex. Instead, lithospheric removal in the form of lithospheric drips would explain spatially sporadic present-day volcanism (e.g., Schleiffarth et al., 2018), in a region that has been influenced by prolonged subduction-related hydration and magmatism. Complex fast wave speed anomalies below both central and eastern Anatolia may thus be a combination of fragmenting subducted African plate, and dense, drips of delaminated Anatolian lower lithosphere. Lithospheric removal through dripping of continental Anatolian lithospheric material rather than peel-back of an oceanic Pontide slab could, like shallow break-off of a continental Bitlis slab, also satisfy the continental crust contamination present in CAVP and EAVP basalts (e.g., Kocaarslan & Ersoy, 2018).

Implications for Uplift in Anatolia-Isostatic Considerations
Warm asthenosphere impinging on the crust in response to lithospheric delamination may not only play a role in magmatism but also in regional uplift. According to various seismic studies, central and eastern Anatolia have normal thickness crust (e.g., Abgarmi et al., 2017;Vanacore et al., 2013;Zor, 2008) but elevated topography (1-2 km; e.g., Komut, 2015). The plateaus must thus be isostatically compensated by low-density or thin lithospheric roots or dynamically uplifted by mantle flow.
Residual topography, defined as the difference between the observed elevation and that expected from crustal buoyancy alone (e.g., Gvirtzman & Nur, 2001;Gvirtzman et al., 2016), can reveal how much of a lithospheric mantle root is required to balance crustal buoyancy. We calculated residual topography at each seismograph station using crustal thicknesses from the same, receiver function-derived, a priori crustal model used in our tomographic inversions (see Text S2 for details). Densities of 3,200 and 2,840 kgm −3 (average crustal density for continents; Karabulut et al., 2019) were adopted for asthenosphere and 10.1029/2020GC009009 Geochemistry, Geophysics, Geosystems crust, respectively. Residuals were reduced to a mid-ocean ridge reference level where mantle lithospheric buoyancy is zero (Lachenbruch & Morgan, 1990). Following Yu et al. (2016), it is instructive to define the residual topography expected from a 35 km-thick continental crust under isostatic equilibrium at zero elevation, a common reference for Airy isostasy (Figure 11; Molnar et al., 1993). For crustal buoyancy alone, this reference crust would have a negative residual topography of −1.5 km, implying a 96 km-thick, negatively buoyant (assuming a density of 3,250 kgm −3 ; see Text S3 for calculation) lithospheric mantle is needed to maintain the top of the crust at sea level. Lower values of residual topography imply a thicker mantle lithosphere and/or a downward pull by flow, while higher values imply a thinner lithosphere and/or a component of dynamic uplift.
Results reveal significant lateral variations across central and eastern Anatolia ( Figure 11), with the largest contrast across the NAFZ. Residuals just north of the NAFZ suggest it is close to Airy isostasy (−1.5 km). The lowest residual topographies (<−2 km), located further north, are likely the result of thinned and/or denser-than-reference crust below the Black Sea, which is part oceanic and part heavily intruded continental crust (e.g., Shillington et al., 2009). Outside the Black Sea, mantle lithospheric thicknesses need to be ∼120 km, increasing to 180 km in the northernmost part of the study area, to counteract the negative buoyancy and satisfy observed elevation ( Figure S12). The change in lithospheric thickness across the NAFZ is estimated to be 40-60 km. The low negative residuals (−2 to −3 km) in southernmost Anatolia, below the CAP, could be a consequence of subduction.
The CAP and EAP interiors have maximum values of residual topography (−0.3 to +0.5 km), suggesting the area deviates from Airy isostatic equilibrium and may require other geodynamic processes to maintain elevated topography. Residuals near zero indicate crustal buoyancy alone satisfies the observed elevation; thus, very thin to completely absent mantle lithosphere is required under normal asthenospheric temperatures (akin to mid-ocean ridges). Positive residuals in the innermost portion of the plateaus (stars in Figure 11) require buoyant material at depth or additional dynamic support to satisfy the elevation. Changing the density contrasts between asthenosphere, lithospheric mantle and crust impacts our calculated residual topographies and lithospheric thicknesses: ± 30 kgm −3 in crustal density causes ± 0.37 km uncertainty in residual topography and ± 21 km in lithospheric mantle thicknesses. Our interpretations, however, remain robust with some stations in the CAP and EAP still characterized by positive residuals (Figure 11).
Asthenospheric support, perhaps from hotter-than-average mantle replacing a convectively removed lower lithosphere plus melt buoyancy, could thus be required to fully account for the uplift in parts of the CAP and EAP. Lithospheric thicknesses derived from both surface wave inversions (Şengör et al., 2003) and joint inversion of receiver functions and surface waves (Gök et al., 2007) also indicate a significant amount of mantle lithospheric material to be missing, and almost absent beneath the highest elevation areas in the EAP (e.g., Gök et al., 2007;Şengör et al., 2003). Karabulut et al. (2019) find a low-density contrast between the crust and mantle (∼315 kgm −3 ) from isostatic calculations, which they attribute to the thermal state of the mantle arising from an anomalously thin lithospheric lid, that could also be supporting the high topographies of the CAP and EAP. Local subduction-collision processes explain adequately most Anatolian tectonics, and although studies modeling dynamic topography (Komut, 2015;Şengül Uluocak et al., 2016) attribute temperature variations in the near-field upper mantle as a primary controlling factor of elevation, they do not preclude the possibility of far-field contributions, such as flow from the Afar mantle plume(e.g., Faccenna et al., 2013).
Our low wave speed anomalies indicating elevated mantle temperatures and the presence of melt below central and eastern Anatolia (Figures 9b and 9c) imply that a combination of thinned lithospheric mantle and elevated asthenospheric temperatures is likely responsible for large deviations in isostatic compensation.

Geochemistry, Geophysics, Geosystems
When reviewed in light of positive long-wavelength free-air gravity anomalies and negative Bouguer anomalies concurrent with mass deficit (Bonvalot et al., 2012), our observed slow P and S wave anomalies (Figures 7 and 8) are consistent with the presence of buoyant material within the uppermost mantle beneath Anatolia (Bakırcı et al., 2012;Berk Biryol et al., 2011;Govers & Fichtner, 2016).

Conclusions
We present the highest-resolution P wave, to date, and first S wave relative arrival-time tomography models of the Eastern Mediterranean. Results reveal the complex structure of the subducting African lithosphere, suggesting slabs do not necessarily remain cohesive: They can instead tear as they sink toward the lower mantle, likely in response to lateral variations in slab properties. Figure 12 illustrates our main conclusions.
Abrupt lateral wave speed contrasts to 100-150 km depth are observed across the NW segment of the NAFZ, associated with fast wave speed (δV p = 3%;δV s = 5%), cold (∼800°C), likely refractory Precambrian lithosphere to the north and slower (δV p = −1%; δV s =−3) asthenosphere below Anatolia to the south. Slow structure (δV p = −1%; δV s = −2%) along strike of the DST may be indicative of thinned lithosphere, which may provide a pathway for hot spot material flowing from Afar as far north as eastern Anatolia.
The Aegean and Cyprean slabs, dipping at 60°and 40°, respectively, are separated by slow anomalies (δV p = −2%; δV s = −3%) defining gaps and windows that permit asthenospheric upwelling beneath the tectonically complex Isparta Angle. Where continental lithosphere is subducting east of Cyprus, horizontal slab tearing at ∼38°N at ∼200 km depth correlates with the CAVP. The Cyprean slab is therefore in early detachment and may develop akin to the Bitlis slab, which sits at 500 km depth. Observations either support a 50± 8 mm year −1 slab sinking rate for the Bitlis slab based on a Miocene (10 Ma) break-off age or are consistent with an earlier Oligocene (∼26 Ma) break-off age with sinking rates more in line with convergence rates (∼20 mm year −1 ). All three slabs interact with the MTZ, descending into the lower mantle, flattening or thickening slightly as they do so.
The slowest wave speeds (δV p = −3%; δV s = −5%) at <200 km depth are found below active volcanic centers in central and eastern Anatolia and, when considered in light of residual topography calculations, are interpreted as asthenospheric material below thinned continental lithosphere (0-20 km-thick mantle lithosphere below the crust). Relative traveltime residuals and velocity-derived temperatures at 100 km depth in excess of those from petrological estimates require the presence of mantle melt. Based on these observations,the

10.1029/2020GC009009
Geochemistry, Geophysics, Geosystems distribution of deformational domains throughout Anatolia correlates directly with upper-mantle slab features such as tears, edges, and windows.